#include "../moldyn.h"
#include "../math/math.h"
-//#include "harmonic_oscillator.h"
+#include "harmonic_oscillator.h"
int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
t_3dvec force,distance;
double d,f;
double sc,equi_dist;
+ double energy;
- params=moldyn->pot2b_params;
+ params=moldyn->pot_params;
sc=params->spring_constant;
equi_dist=params->equilibrium_distance;
if(bc) check_per_bound(moldyn,&distance);
d=v3_norm(&distance);
if(d<=moldyn->cutoff) {
- moldyn->energy+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
+ energy=(0.5*sc*(d-equi_dist)*(d-equi_dist));
+ moldyn->energy+=energy;
+ ai->e+=0.5*energy;
+ aj->e+=0.5*energy;
/* f = -grad E; grad r_ij = -1 1/r_ij distance */
f=sc*(1.0-equi_dist/d);
v3_scale(&force,&distance,f);