X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=potentials%2Falbe.c;h=688712dc0c2da74ae1473d96de2604b18c43d555;hb=684bf7c398cdfa98549b0c7a1fa37e6dc5b35bea;hp=a79f58d3268834f91aaead5cfc26ba56c6334929;hpb=92ef07d77a4c879527180224acea73a3f6564497;p=physik%2Fposic.git

diff --git a/potentials/albe.c b/potentials/albe.c
index a79f58d..688712d 100644
--- a/potentials/albe.c
+++ b/potentials/albe.c
@@ -98,7 +98,7 @@ int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
 /* albe 3 body potential function (first k loop) */
 int albe_mult_3bp_k1(t_moldyn *moldyn,
-                        t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
+                     t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 
 	t_albe_mult_params *params;
 	t_albe_exchange *exchange;
@@ -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);
@@ -283,7 +284,7 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	}
 
 	/* force contribution */
-	scale=-0.5*(f_c*(df_r+b*df_a)+df_c*(f_r+b*f_a));
+	scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a));
 	v3_scale(&force,&(exchange->dist_ij),scale);
 	v3_add(&(ai->f),&(ai->f),&force);
 	v3_sub(&(aj->f),&(aj->f),&force); // dri rij = - drj rij
@@ -305,11 +306,13 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	/* virial */
 	virial_calc(ai,&force,&(exchange->dist_ij));
 
-	/* dzeta prefactor = - 0.5 f_c f_a db */
-	exchange->pre_dzeta=-0.5*f_a*f_c*db;
+	/* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
+	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;
@@ -319,7 +322,7 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
 /* albe 3 body potential function (second k loop) */
 int albe_mult_3bp_k2(t_moldyn *moldyn,
-                        t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
+                     t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 
 	t_albe_mult_params *params;
 	t_albe_exchange *exchange;
@@ -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++;