/* zero absolute time */
moldyn->time=0.0;
- /* debugging, ignre */
+ /* debugging, ignore */
moldyn->debug=0;
/* executing the schedule */
if(!(i%v)) {
visual_atoms(&(moldyn->vis),moldyn->time,
moldyn->atom,moldyn->count);
- printf("\rsched: %d, steps: %d, theta: %d",
+ printf("\rsched: %d, steps: %d, debug: %d",
sched,i,moldyn->debug);
fflush(stdout);
}
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);
/* 2bp post function */
if(moldyn->func2b_post) {
-printf("DEBUG: vor 2bp post\n");
moldyn->func2b_post(moldyn,
&(itom[i]),
jtom,bc_ij);
-printf("DEBUG: nach 2bp post\n");
}
}
}
+//printf("DEBUG: %.15f \n",itom[i].f.x);
}
return 0;
/* clear 3bp and 2bp post run */
exchange->run3bp=0;
exchange->run2bp_post=0;
+
+ /* reset S > r > R mark */
+ exchange->d_ij_between_rs=0;
/*
* calc of 2bp contribution of V_ij and dV_ij/ji
/* 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);
df_c=-0.5*sin(arg)*(M_PI/(s_r*d_ij));
/* two body contribution (ij, ji) */
v3_scale(&force,&dist_ij,-df_c*f_r-df_r*f_c);
+ /* tell 3bp that S > r > R */
+ exchange->d_ij_between_rs=1;
}
/* add forces of 2bp (ij, ji) contribution
/* Vij and dVij */
zeta=exchange->zeta_ij;
- 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); /* chi * (...)^(-1/2n - 1) */
- b=db*b; /* b_ij */
- db*=-0.5*tmp; /* db_ij */
- v3_scale(&force,&(exchange->dzeta_ij),f_a*db);
- v3_scale(&temp,dist_ij,df_a*b);
- v3_add(&force,&force,&temp);
- v3_scale(&force,&force,f_c);
+ if(zeta==0.0) {
+ 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 */
+ b=(1+zeta*tmp); /* 1 + beta^n zeta^n */
+ db=chi*pow(b,-1.0/(2*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);
+ 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);
/* add energy of 3bp sum */
moldyn->energy+=(0.5*f_c*b*f_a);
- /* add force (sub, as F = - dVij) */
- v3_sub(&(ai->f),&(ai->f),&force);
-
/* dVji */
zeta=exchange->zeta_ji;
- tmp=betajnj*pow(zeta,nj-1.0); /* beta^n * zeta^n-1 */
- b=(1+zeta*tmp); /* 1 + beta^n * zeta^n */
- db=chi*pow(b,-1.0/(2*nj)-1); /* chi * (...)^(-1/2n - 1) */
- b=db*b; /* b_ij */
- db*=-0.5*tmp; /* db_ij */
- v3_scale(&force,&(exchange->dzeta_ji),f_a*db);
- v3_scale(&temp,dist_ij,df_a*b);
- v3_add(&force,&force,&temp);
- v3_scale(&force,&force,f_c);
+ if(zeta==0.0) {
+ 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 */
+ b=(1+zeta*tmp); /* 1 + beta^n zeta^n */
+ db=chi*pow(b,-1.0/(2*nj)-1); /* x(...)^(-1/2n - 1) */
+ b=db*b; /* b_ij */
+ db*=-0.5*tmp; /* db_ij */
+ v3_scale(&force,&(exchange->dzeta_ji),f_a*db);
+ 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 (sub, as F = - dVji) */
+ /* add force */
v3_sub(&(ai->f),&(ai->f),&force);
return 0;
/* betajnj * zeta_jk ^ nj-1 */
tmp=exchange->betajnj*pow(zeta,(n-1.0));
tmp=-chi/2.0*pow(1+tmp*zeta,-1.0/(2.0*n)-1)*tmp;
- v3_scale(&temp1,&temp1,tmp*B*exp(-mu*d_jk)*f_c_jk);
+ v3_scale(&temp1,&temp1,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5);
v3_add(&(ai->f),&(ai->f),&temp1); /* -1 skipped in f_a calc ^ */
+ /* scaled with 0.5 ^ */
+ }
+
+ 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;
projects / physik / posic.git / blobdiff