/* add forces */
v3_add(&(ai->f),&(ai->f),&force);
- v3_sub(&(aj->f),&(aj->f),&force);
+ v3_sub(&(aj->f),&(aj->f),&force); // reason: dri rij = - drj rij
#ifdef DEBUG
if((ai==&(moldyn->atom[0]))|(aj==&(moldyn->atom[0]))) {
ni=*(exchange->n_i);
tmp=exchange->betaini*pow(exchange->zeta_ij,ni-1.0);
b=(1.0+exchange->zeta_ij*tmp);
- db=chi*pow(b,-1.0/(2*ni)-1.0);
+ db=chi*pow(b,-1.0/(2.0*ni)-1.0);
b=db*b;
db*=-0.5*tmp;
}
v3_scale(&force,&(exchange->dist_ij),df_a*f_c+f_a*df_c);
v3_scale(&force,&force,-0.5*b);
v3_add(&(ai->f),&(ai->f),&force);
- v3_sub(&(aj->f),&(aj->f),&force);
+ v3_sub(&(aj->f),&(aj->f),&force); // dri rij = - drj rij
#ifdef DEBUG
if((ai==&(moldyn->atom[0]))|(aj==&(moldyn->atom[0]))) {
/* virial */
//virial_calc(ai,&force,&dist_ij);
- /* derivatice wrt j */
+ /* derivative wrt j */
v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
/* force contribution */
//virial_calc(aj,&force,&dist_ij);
/* derivative wrt k */
- v3_scale(&force,&dist_ik,dfcg);
+ v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rij = - drj rij
v3_scale(&tmp,&dcosdrk,fcdg);
v3_add(&force,&force,&tmp);
v3_scale(&force,&force,pre_dzeta);
else {
tmp=betaini*pow(zeta,ni-1.0); /* beta^n * zeta^n-1 */
b=(1+zeta*tmp); /* 1 + beta^n zeta^n */
- db=chi*pow(b,-1.0/(2*ni)-1); /* x(...)^(-1/2n - 1) */
+ db=chi*pow(b,-1.0/(2.0*ni)-1); /* x(...)^(-1/2n - 1) */
b=db*b; /* b_ij */
db*=-0.5*tmp; /* db_ij */
v3_scale(&force,&(exchange->dzeta_ij),f_a*db);
// ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
// ATOM_ATTR_2BP|ATOM_ATTR_HB,
// &r,&v);
+ //r.z=0.27*sqrt(3.0)*LC_SI/2.0; v.z=0;
+ //r.x=(TM_S_SI+TM_R_SI)/4.0; v.x=0;
+ //r.y=0; v.y=0;
+ //r.x=0; v.x=0;
+ //add_atom(&md,SI,M_SI,0,
+ // ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
+ // ATOM_ATTR_2BP|ATOM_ATTR_HB,
+ // &r,&v);
+ //r.z=-r.z; v.z=-v.z;
+ //r.y=0; v.y=0;
//r.x=0; v.x=0;
- //r.y=0; v.y=0;
- //r.z=sin(M_PI*60.0/180.0)*(TM_S_SI+TM_R_SI)/4.0; v.z=0;
//add_atom(&md,SI,M_SI,0,
// ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
// ATOM_ATTR_2BP|ATOM_ATTR_HB,