small mods to support site energies and kinetic energies per atom

[physik/posic.git] / potentials / harmonic_oscillator.c

diff --git a/potentials/harmonic_oscillator.c b/potentials/harmonic_oscillator.c
index 0d933b6..b203453 100644 (file)
--- a/potentials/harmonic_oscillator.c
+++ b/potentials/harmonic_oscillator.c
@@ -17,7 +17,7 @@
 
 #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) {
 
@@ -25,8 +25,9 @@ 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;
 
@@ -37,7 +38,10 @@ int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        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);