/* albe 3 body potential function (first k loop) */
int albe_mult_3bp_k1(t_moldyn *moldyn,
/* 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) {
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
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
- /* 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;
/* albe 3 body potential function (second k loop) */
int albe_mult_3bp_k2(t_moldyn *moldyn,
/* 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) {
/* cos_theta derivatives wrt i,j,k */
dijdik_inv=1.0/(d_ij*d_ik);
/* cos_theta derivatives wrt i,j,k */
dijdik_inv=1.0/(d_ij*d_ik);
v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
v3_add(&dcosdrj,&dcosdrj,&tmp);
v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
v3_add(&dcosdrj,&dcosdrj,&tmp);
v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
v3_add(&dcosdrk,&dcosdrk,&tmp);
v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
v3_add(&dcosdrk,&dcosdrk,&tmp);
v3_scale(&dcosdri,&dcosdri,-1.0);
/* f_c_ik * dg, df_c_ik * g */
v3_scale(&dcosdri,&dcosdri,-1.0);
/* f_c_ik * dg, df_c_ik * g */
- //v3_scale(&force,&force,-1.0);
- //virial_calc(ai,&force,&dist_ij);
+ v3_scale(&force,&force,-1.0);
+ virial_calc(ai,&force,&dist_ij);
- //v3_scale(&force,&force,-1.0);
- //virial_calc(ai,&force,&dist_ik);
+ v3_scale(&force,&force,-1.0);
+ virial_calc(ai,&force,&dist_ik);