small mods to support site energies and kinetic energies per atom

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

diff --git a/potentials/tersoff.c b/potentials/tersoff.c
index 6bfcf30..953b839 100644 (file)
--- a/potentials/tersoff.c
+++ b/potentials/tersoff.c
@@ -93,6 +93,7 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        int brand;
        double s_r;
        double arg;
+       double energy;
 
        printf("WARNING! - tersoff_mult_2bp is obsolete.\n");
        printf("WARNING! - repulsive part handled in 3bp/j2 routine.\n");
@@ -173,7 +174,10 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        virial_calc(ai,&force,&dist_ij);
 
        /* energy 2bp contribution */
-       moldyn->energy+=f_r*f_c;
+       energy=f_r*f_c;
+       moldyn->energy+=energy;
+       ai->e+=0.5*energy;
+       aj->e+=0.5*energy;
 
        return 0;
 }
@@ -348,6 +352,7 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        unsigned char brand;
        double ni,tmp;
        double S,R,s_r,arg;
+       double energy;
 
        params=moldyn->pot_params;
        exchange=&(params->exchange);
@@ -436,7 +441,9 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        exchange->pre_dzeta=-0.5*f_a*f_c*db;
 
        /* energy contribution */
-       moldyn->energy+=0.5*f_c*(b*f_a+f_r);
+       energy=0.5*f_c*(b*f_a+f_r);
+       moldyn->energy+=energy;
+       ai->e+=energy;
 
        /* reset k counter for second k loop */
        exchange->kcount=0;