v3_add(&(atom[i].r),&(atom[i].r),&delta);
v3_scale(&delta,&(atom[i].f),0.5*tau_square/atom[i].mass);
v3_add(&(atom[i].r),&(atom[i].r),&delta);
+//if(i==0) printf("und dann %.20f\n",atom[i].r.x*1e10);
check_per_bound(moldyn,&(atom[i].r));
+//if(i==0) printf("und danach %.20f\n",atom[i].r.x);
/* velocities */
v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
v3_add(&(atom[i].v),&(atom[i].v),&delta);
}
+moldyn_bc_check(moldyn);
/* neighbour list update */
link_cell_update(moldyn);
/* f_r_ij = f_r_ji, df_r_ij = df_r_ji */
f_r=A*exp(-lambda*d_ij);
- df_r=-lambda*f_r/d_ij;
+ df_r=lambda*f_r/d_ij;
/* f_a, df_a calc (again, same for ij and ji) | save for later use! */
exchange->f_a=-B*exp(-mu*d_ij);
moldyn->debug++; /* just for debugging ... */
db=0.0;
b=chi;
+ v3_scale(&force,dist_ij,df_a*b*f_c);
}
else {
tmp=betaini*pow(zeta,ni-1.0); /* beta^n * zeta^n-1 */
v3_scale(&temp,dist_ij,df_a*b);
v3_add(&force,&force,&temp);
v3_scale(&force,&force,f_c);
- v3_scale(&temp,dist_ij,df_c*b*f_a);
- v3_add(&force,&force,&temp);
- v3_scale(&force,&force,-0.5);
-
- /* add force */
- v3_add(&(ai->f),&(ai->f),&force);
}
+ v3_scale(&temp,dist_ij,df_c*b*f_a);
+ v3_add(&force,&force,&temp);
+ v3_scale(&force,&force,-0.5);
+
+ /* add force */
+ v3_add(&(ai->f),&(ai->f),&force);
/* add energy of 3bp sum */
moldyn->energy+=(0.5*f_c*b*f_a);
/* dVji */
zeta=exchange->zeta_ji;
if(zeta==0.0) {
- // do nothing ... (db=0, b=chi)
moldyn->debug++;
+ b=chi;
+ v3_scale(&force,dist_ij,df_a*b*f_c);
}
else {
tmp=betajnj*pow(zeta,nj-1.0); /* beta^n * zeta^n-1 */
v3_scale(&temp,dist_ij,df_a*b);
v3_add(&force,&force,&temp);
v3_scale(&force,&force,f_c);
- v3_scale(&temp,dist_ij,0.5*df_c*b*f_a);
- v3_add(&force,&force,&temp);
- v3_scale(&force,&force,-0.5);
-
- /* add force */
- v3_sub(&(ai->f),&(ai->f),&force);
}
+ v3_scale(&temp,dist_ij,df_c*b*f_a);
+ v3_add(&force,&force,&temp);
+ v3_scale(&force,&force,-0.5);
+
+ /* add force */
+ v3_sub(&(ai->f),&(ai->f),&force);
return 0;
}
return 0;
}
+
+/*
+ * debugging / critical check functions
+ */
+
+int moldyn_bc_check(t_moldyn *moldyn) {
+
+ t_atom *atom;
+ t_3dvec *dim;
+ int i;
+
+ atom=moldyn->atom;
+ dim=&(moldyn->dim);
+
+ for(i=0;i<moldyn->count;i++) {
+ if(atom[i].r.x>=dim->x/2||-atom[i].r.x>dim->x/2)
+ printf("FATAL: atom %d: x: %.20f (%.20f)\n",
+ i,atom[i].r.x*1e10,dim->x/2*1e10);
+ if(atom[i].r.y>=dim->y/2||-atom[i].r.y>dim->y/2)
+ printf("FATAL: atom %d: y: %.20f (%.20f)\n",
+ i,atom[i].r.y*1e10,dim->y/2*1e10);
+ if(atom[i].r.z>=dim->z/2||-atom[i].r.z>dim->z/2)
+ printf("FATAL: atom %d: z: %.20f (%.20f)\n",
+ i,atom[i].r.z*1e10,dim->z/2*1e10);
+ }
+
+ return 0;
+}
+