!(jtom->attr&ATOM_ATTR_3BP))
continue;
- /*
- * according to mr. nordlund, we dont need to take the
- * sum over all atoms now, as 'this is centered' around
- * atom i ...
- * i am not quite sure though! there is a not vanishing
- * part even if f_c_ik is zero ...
- * this analytical potentials suck!
- * switching from mc to md to dft soon!
- */
+ /* neighbourhood of atom j is not needed! */
// link_cell_neighbour_index(moldyn,
// (jtom->r.x+moldyn->dim.x/2)/lc->x,
exchange=&(params->exchange);
exchange->run3bp=0;
+ exchange->run2bp_post=0;
/*
* we need: f_c, df_c, f_r, df_r
f_r=A*exp(-lambda*d_ij);
df_r=-lambda*f_r/d_ij;
- /* f_a, df_a calc + save for 3bp use */
+ /* f_a, df_a calc + save for later use */
exchange->f_a=-B*exp(-mu*d_ij);
exchange->df_a=-mu*exchange->f_a/d_ij;
exchange->f_c=f_c;
exchange->df_c=df_c;
- /* enable the run of 3bp function */
+ /* enable the run of 3bp function and 2bp post processing */
exchange->run3bp=1;
+ exchange->run2bp_post=1;
/* reset 3bp sums */
exchange->sum1_3bp=0.0;
params=moldyn->pot2b_params;
exchange=&(params->exchange);
+ /* we do not run if f_c_ij was dtected to be 0! */
+ if(!(exchange->run2bp_post))
+ return 0;
+
db_ij=&(exchange->db_ij);
f_c=exchange->f_c;
df_c=exchange->df_c;
df_a=exchange->df_a;
betan=exchange->betan;
n=*(exchange->n);
+ chi=exchange->chi;
dist_ij=&(exchange->dist_ij);
db_ij_scale1=(1+betan*exchange->sum1_3bp);
b_ij=chi*db_ij_scale1*help;
db_ij_scale1=-chi/(2*n)*help;
+ /* db_ij part */
v3_scale(db_ij,db_ij,(db_ij_scale1*db_ij_scale2));
v3_scale(db_ij,db_ij,f_a);
+ /* df_a part */
v3_scale(&temp,dist_ij,b_ij*df_a);
+ /* db_ij + df_a part */
v3_add(&force,&temp,db_ij);
v3_scale(&force,&force,f_c);
+ /* df_c part */
v3_scale(&temp,dist_ij,f_a*b_ij*df_c);
/* add energy of 3bp sum */
moldyn->energy+=(0.5*f_c*b_ij*f_a);
- /* add force of 3bp calculation */
+
+ /* add force of 3bp calculation (all three parts) */
v3_add(&(ai->f),&temp,&force);
return 0;
bracket=0.0;
bracket_n_1=0.0;
bracket_n=0.0;
- printf("Foo -> 0: ");
+ //printf("Foo -> 0: ");
}
else {
bracket=f_c_ik*g;
bracket_n_1=pow(bracket,n-1.0);
bracket_n=bracket_n_1*bracket;
- printf("Foo -> 1: ");
+ //printf("Foo -> 1: ");
}
//printf("%.15f %.15f %.15f\n",bracket_n_1,bracket_n,bracket);