int moldyn_shutdown(t_moldyn *moldyn) {
+ moldyn_log_shutdown(moldyn);
link_cell_shutdown(moldyn);
moldyn_log_shutdown(moldyn);
rand_close(&(moldyn->random));
strncpy(moldyn->sfb,fb,63);
break;
case VISUAL_STEP:
- moldyn->mwrite=timer;
+ moldyn->vwrite=timer;
strncpy(moldyn->vfb,fb,63);
visual_init(&(moldyn->vis),fb);
break;
}
int add_atom(t_moldyn *moldyn,int element,double mass,u8 bnum,u8 attr,
- t_3dvec r,t_3dvec v) {
+ t_3dvec *r,t_3dvec *v) {
t_atom *atom;
void *ptr;
perror("[moldyn] realloc (add atom)");
return -1;
}
-
- atom=ptr;
- atom->r=r;
- atom->v=v;
+ moldyn->atom=ptr;
+
+ atom=moldyn->atom;
+ atom->r=*r;
+ atom->v=*v;
atom->element=element;
atom->bnum=bnum;
atom->attr=attr;
return p_total;
}
-double estimate_time_step(t_moldyn *moldyn,double nn_dist,double t) {
+double estimate_time_step(t_moldyn *moldyn,double nn_dist) {
double tau;
- tau=0.05*nn_dist/(sqrt(3.0*K_BOLTZMANN*t/moldyn->atom[0].mass));
- tau*=1.0E-9;
- if(tau<moldyn->tau)
- printf("[moldyn] warning: time step (%f > %.15f)\n",
- moldyn->tau,tau);
+ /* nn_dist is the nearest neighbour distance */
+
+ if(moldyn->t==5.0) {
+ printf("[moldyn] i do not estimate timesteps below %f K!\n",
+ MOLDYN_CRITICAL_EST_TEMP);
+ return 23.42;
+ }
+
+ tau=(0.05*nn_dist*moldyn->atom[0].mass)/sqrt(3.0*K_BOLTZMANN*moldyn->t);
return tau;
}
t_linkcell *lc;
int i;
+ int fd;
+
+ fd=open("/dev/null",O_WRONLY);
lc=&(moldyn->lc);
lc->cells=lc->nx*lc->ny*lc->nz;
lc->subcell=malloc(lc->cells*sizeof(t_list));
- printf("initializing linked cells (%d)\n",lc->cells);
+ printf("[moldyn] initializing linked cells (%d)\n",lc->cells);
for(i=0;i<lc->cells;i++)
- list_init(&(lc->subcell[i]),1);
+ //list_init(&(lc->subcell[i]),1);
+ list_init(&(lc->subcell[i]),fd);
link_cell_update(moldyn);
perror("[moldyn] realloc (runs)");
return -1;
}
+ moldyn->schedule.runs=ptr;
moldyn->schedule.runs[count-1]=runs;
ptr=realloc(schedule->tau,count*sizeof(double));
perror("[moldyn] realloc (tau)");
return -1;
}
+ moldyn->schedule.tau=ptr;
moldyn->schedule.tau[count-1]=tau;
return 0;
}
if(v) {
if(!(i%v)) {
- visual_atoms(moldyn->visual,i*moldyn->tau,
+ visual_atoms(&(moldyn->vis),moldyn->time,
moldyn->atom,moldyn->count);
- printf("\rsteps: %d",i);
+ printf("\rsched: %d, steps: %d",sched,i);
fflush(stdout);
}
}
+
}
/* check for hooks */
}
/* neighbour list update */
-printf("list update ...\n");
link_cell_update(moldyn);
-printf("done\n");
/* forces depending on chosen potential */
-printf("calc potential/force ...\n");
potential_force_calc(moldyn);
//moldyn->potential_force_function(moldyn);
-printf("done\n");
for(i=0;i<count;i++) {
/* again velocities */
/* reset energy */
moldyn->energy=0.0;
+printf("DEBUG: count = %d\n",count);
for(i=0;i<count;i++) {
/* reset force */
/* 2 body pair potential/force */
if(atom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)) {
-
+
+printf("DEBUG: processing atom %d\n",i);
link_cell_neighbour_index(moldyn,
(atom[i].r.x+moldyn->dim.x/2)/lc->x,
(atom[i].r.y+moldyn->dim.y/2)/lc->y,
countn=lc->countn;
dnlc=lc->dnlc;
+printf("DEBUG: countn = %d - dnslc = %d\n",countn,dnlc);
for(j=0;j<countn;j++) {
this=&(neighbour[j]);
if((btom->attr&ATOM_ATTR_2BP)&
(atom[i].attr&ATOM_ATTR_2BP))
+printf("DEBUG: calling func2b\n");
moldyn->func2b(moldyn,
&(atom[i]),
btom,