small mods to support site energies and kinetic energies per atom

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

diff --git a/potentials/albe.c b/potentials/albe.c
index a79f58d..28008b0 100644 (file)
--- a/potentials/albe.c
+++ b/potentials/albe.c
@@ -226,6 +226,7 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        double d_ij,r0;
        unsigned char brand;
        double S,R,s_r,arg;
+       double energy;
 
        params=moldyn->pot_params;
        exchange=&(params->exchange);
@@ -309,7 +310,9 @@ int albe_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*(f_r+b*f_a);
+       energy=0.5*f_c*(f_r+b*f_a);
+       moldyn->energy+=energy;
+       ai->e+=energy;
 
        /* reset k counter for second k loop */
        exchange->kcount=0;
@@ -380,13 +383,13 @@ int albe_mult_3bp_k2(t_moldyn *moldyn,
 
        /* cos_theta derivatives wrt i,j,k */
        dijdik_inv=1.0/(d_ij*d_ik);
-       v3_scale(&dcosdrj,&dist_ik,dijdik_inv);
+       v3_scale(&dcosdrj,&dist_ik,dijdik_inv);         // j
        v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
        v3_add(&dcosdrj,&dcosdrj,&tmp);
-       v3_scale(&dcosdrk,&dist_ij,dijdik_inv);
+       v3_scale(&dcosdrk,&dist_ij,dijdik_inv);         // k
        v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
        v3_add(&dcosdrk,&dcosdrk,&tmp);
-       v3_add(&dcosdri,&dcosdrj,&dcosdrk);
+       v3_add(&dcosdri,&dcosdrj,&dcosdrk);             // i
        v3_scale(&dcosdri,&dcosdri,-1.0);
 
        /* f_c_ik * dg, df_c_ik * g */
@@ -428,8 +431,8 @@ int albe_mult_3bp_k2(t_moldyn *moldyn,
 #endif
 
        /* virial */
-       //v3_scale(&force,&force,-1.0);
-       //virial_calc(ai,&force,&dist_ij);
+       v3_scale(&force,&force,-1.0);
+       virial_calc(ai,&force,&dist_ij);
 
        /* derivative wrt k */
        v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
@@ -449,8 +452,8 @@ int albe_mult_3bp_k2(t_moldyn *moldyn,
 #endif
 
        /* virial */
-       //v3_scale(&force,&force,-1.0);
-       //virial_calc(ai,&force,&dist_ik);
+       v3_scale(&force,&force,-1.0);
+       virial_calc(ai,&force,&dist_ik);
        
        /* increase k counter */
        exchange->kcount++;