t_moldyn_schedule *schedule;
schedule=&(moldyn->schedule);
- count=++(schedule->content_count);
+ count=++(schedule->total_sched);
- ptr=realloc(moldyn->schedule.runs,count*sizeof(int));
+ ptr=realloc(schedule->runs,count*sizeof(int));
if(!ptr) {
perror("[moldyn] realloc (runs)");
return -1;
}
- moldyn->schedule.runs=ptr;
- moldyn->schedule.runs[count-1]=runs;
+ schedule->runs=ptr;
+ schedule->runs[count-1]=runs;
ptr=realloc(schedule->tau,count*sizeof(double));
if(!ptr) {
perror("[moldyn] realloc (tau)");
return -1;
}
- moldyn->schedule.tau=ptr;
- moldyn->schedule.tau[count-1]=tau;
+ schedule->tau=ptr;
+ schedule->tau[count-1]=tau;
return 0;
}
int moldyn_integrate(t_moldyn *moldyn) {
- int i,sched;
+ int i;
unsigned int e,m,s,v;
t_3dvec p;
- t_moldyn_schedule *schedule;
+ t_moldyn_schedule *sched;
t_atom *atom;
int fd;
char dir[128];
double ds;
- schedule=&(moldyn->schedule);
+ sched=&(moldyn->schedule);
atom=moldyn->atom;
/* initialize linked cell method */
moldyn->debug=0;
/* executing the schedule */
- for(sched=0;sched<moldyn->schedule.content_count;sched++) {
+ for(sched->count=0;sched->count<sched->total_sched;sched->count++) {
/* setting amount of runs and finite time step size */
- moldyn->tau=schedule->tau[sched];
+ moldyn->tau=sched->tau[sched->count];
moldyn->tau_square=moldyn->tau*moldyn->tau;
- moldyn->time_steps=schedule->runs[sched];
+ moldyn->time_steps=sched->runs[sched->count];
/* integration according to schedule */
visual_atoms(&(moldyn->vis),moldyn->time,
moldyn->atom,moldyn->count);
printf("\rsched: %d, steps: %d, debug: %d",
- sched,i,moldyn->debug);
+ sched->count,i,moldyn->debug);
fflush(stdout);
}
}
}
/* check for hooks */
- if(schedule->hook)
- schedule->hook(moldyn,schedule->hook_params);
+ if(sched->hook)
+ sched->hook(moldyn,sched->hook_params);
/* get a new info line */
printf("\n");
}
}
+#ifdef DEBUG
+printf("\n\n");
+#endif
return 0;
}
/* add forces of 2bp (ij, ji) contribution
* dVij = dVji and we sum up both: no 1/2) */
v3_add(&(ai->f),&(ai->f),&force);
+#ifdef DEBUG
+if(ai==&(moldyn->atom[0])) {
+ printf("dVij, dVji (2bp) contrib:\n");
+ printf("%f | %f\n",force.x,ai->f.x);
+ printf("%f | %f\n",force.y,ai->f.y);
+ printf("%f | %f\n",force.z,ai->f.z);
+}
+#endif
/* energy 2bp contribution (ij, ji) is 0.5 f_r f_c ... */
moldyn->energy+=(0.5*f_r*f_c);
/* add force */
v3_add(&(ai->f),&(ai->f),&force);
+#ifdef DEBUG
+if(ai==&(moldyn->atom[0])) {
+ printf("dVij (3bp) contrib:\n");
+ printf("%f | %f\n",force.x,ai->f.x);
+ printf("%f | %f\n",force.y,ai->f.y);
+ printf("%f | %f\n",force.z,ai->f.z);
+}
+#endif
/* add energy of 3bp sum */
moldyn->energy+=(0.5*f_c*b*f_a);
/* add force */
v3_add(&(ai->f),&(ai->f),&force);
+#ifdef DEBUG
+if(ai==&(moldyn->atom[0])) {
+ printf("dVji (3bp) contrib:\n");
+ printf("%f | %f\n",force.x,ai->f.x);
+ printf("%f | %f\n",force.y,ai->f.y);
+ printf("%f | %f\n",force.z,ai->f.z);
+}
+#endif
return 0;
}
v3_scale(&temp2,&temp2,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5);
v3_add(&(ai->f),&(ai->f),&temp2); /* -1 skipped in f_a calc ^ */
/* scaled with 0.5 ^ */
+#ifdef DEBUG
+if(ai==&(moldyn->atom[0])) {
+ printf("dVjk (3bp) contrib:\n");
+ printf("%f | %f\n",temp2.x,ai->f.x);
+ printf("%f | %f\n",temp2.y,ai->f.y);
+ printf("%f | %f\n",temp2.z,ai->f.z);
+}
+#endif
}
int main(int argc,char **argv) {
/* check argv */
- if(argc!=2) {
- printf("[sic] error: arg1 (vis/log/save location) ");
- printf("must be given!\n");
+ if(argc!=3) {
+ printf("[sic] usage: %s <logdir> <temperatur>\n",argv[0]);
return -1;
}
double tau;
/* testing location & velocity vector */
- //t_3dvec r,v;
+ t_3dvec r,v;
/* values */
tau=1.0e-15; /* delta t = 1 fs */
moldyn_bc_check(&md);
/* testing configuration */
- //r.x=2.45/2; v.x=0;
+ //r.x=2.95/2; v.x=0;
//r.y=0; v.y=0;
//r.z=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_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP,//|ATOM_ATTR_HB,
// &r,&v);
- //r.x=-2.45/2; v.x=0;
+ //r.x=-2.95/2; v.x=0;
//r.y=0; v.y=0;
//r.z=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_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP,//|ATOM_ATTR_HB,
// &r,&v);
/* setting a nearest neighbour distance for the moldyn checks */
set_nn_dist(&md,0.25*sqrt(3.0)*LC_SI); /* diamond ! */
/* set temperature */
- printf("[sic] setting temperature\n");
+ printf("[sic] setting temperature -> %f\n",273+atof(argv[2]));
//set_temperature(&md,273.0+1410.0);
- set_temperature(&md,273.0+450.0);
+ //set_temperature(&md,273.0+450.0);
//set_temperature(&md,273.0);
//set_temperature(&md,1.0);
//set_temperature(&md,0.0);
+ set_temperature(&md,atof(argv[2])+273.0);
/* set pressure */
printf("[sic] setting pressure\n");
/* create the simulation schedule */
printf("[sic] adding schedule\n");
- moldyn_add_schedule(&md,20001,1.0);
+ moldyn_add_schedule(&md,20000,.1);
+ moldyn_add_schedule(&md,10000,.2);
+ moldyn_add_schedule(&md,6667,.3);
+ moldyn_add_schedule(&md,5000,.4);
+ moldyn_add_schedule(&md,4001,.5);
/* activate logging */
printf("[sic] activate logging\n");
moldyn_set_log_dir(&md,argv[1]);
- moldyn_set_log(&md,LOG_TOTAL_ENERGY,200);
- moldyn_set_log(&md,VISUAL_STEP,200);
+ moldyn_set_log(&md,LOG_TOTAL_ENERGY,100);
+ moldyn_set_log(&md,VISUAL_STEP,100);
/*
* let's do the actual md algorithm now