int moldyn_shutdown(t_moldyn *moldyn) {
- link_cell_shutdown(moldyn);
+ printf("[moldyn] shutdown\n");
moldyn_log_shutdown(moldyn);
+ link_cell_shutdown(moldyn);
rand_close(&(moldyn->random));
free(moldyn->atom);
}
int set_temperature(t_moldyn *moldyn,double t) {
-
+
moldyn->t=t;
return 0;
return 0;
}
+int set_nn_dist(t_moldyn *moldyn,double dist) {
+
+ moldyn->nnd=dist;
+
+ return 0;
+}
+
int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z) {
if(x)
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 moldyn_log_shutdown(t_moldyn *moldyn) {
+ printf("[moldyn] log shutdown\n");
if(moldyn->efd) close(moldyn->efd);
if(moldyn->mfd) close(moldyn->mfd);
- if(moldyn->visual) visual_tini(moldyn->visual);
+ if(&(moldyn->vis)) visual_tini(&(moldyn->vis));
return 0;
}
int count;
int ret;
t_3dvec origin;
- t_atom *atom;
count=a*b*c;
- atom=moldyn->atom;
if(type==FCC) count*=4;
if(type==DIAMOND) count*=8;
- atom=malloc(count*sizeof(t_atom));
- if(atom==NULL) {
+ moldyn->atom=malloc(count*sizeof(t_atom));
+ if(moldyn->atom==NULL) {
perror("malloc (atoms)");
return -1;
}
switch(type) {
case FCC:
- ret=fcc_init(a,b,c,lc,atom,&origin);
+ ret=fcc_init(a,b,c,lc,moldyn->atom,&origin);
break;
case DIAMOND:
- ret=diamond_init(a,b,c,lc,atom,&origin);
+ ret=diamond_init(a,b,c,lc,moldyn->atom,&origin);
break;
default:
printf("unknown lattice type (%02x)\n",type);
}
moldyn->count=count;
+ printf("[moldyn] created lattice with %d atoms\n",count);
while(count) {
- atom[count-1].element=element;
- atom[count-1].mass=mass;
- atom[count-1].attr=attr;
- atom[count-1].bnum=bnum;
count-=1;
+ moldyn->atom[count].element=element;
+ moldyn->atom[count].mass=mass;
+ moldyn->atom[count].attr=attr;
+ moldyn->atom[count].bnum=bnum;
+ check_per_bound(moldyn,&(moldyn->atom[count].r));
}
+
return ret;
}
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;
- atom->element=element;
- atom->bnum=bnum;
- atom->attr=attr;
+ moldyn->atom=ptr;
+
+ atom=moldyn->atom;
+ atom[count-1].r=*r;
+ atom[count-1].v=*v;
+ atom[count-1].element=element;
+ atom[count-1].mass=mass;
+ atom[count-1].bnum=bnum;
+ atom[count-1].attr=attr;
return 0;
}
}
/* velocity scaling */
- scale_velocity(moldyn);
+ scale_velocity(moldyn,VSCALE_INIT_EQUI);
return 0;
}
-int scale_velocity(t_moldyn *moldyn) {
+int scale_velocity(t_moldyn *moldyn,u8 type) {
int i;
- double e,c;
+ double e,scale;
t_atom *atom;
+ int count;
atom=moldyn->atom;
/*
* - velocity scaling (E = 3/2 N k T), E: kinetic energy
*/
+
e=0.0;
+ count=0;
+ for(i=0;i<moldyn->count;i++) {
+ if(atom[i].attr&ATOM_ATTR_HB) {
+ e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
+ count+=1;
+ }
+ }
+
+ /* temporary hack for e,t = 0 */
+ if(e==0.0) {
+ if(moldyn->t!=0.0)
+ thermal_init(moldyn);
+ else
+ return 0;
+ }
+
+ /* direct scaling */
+ scale=(1.5*count*K_BOLTZMANN*moldyn->t)/e;
+ if(type&VSCALE_INIT_EQUI) scale*=2.0; /* equipartition theorem */
+ scale=sqrt(scale);
for(i=0;i<moldyn->count;i++)
- e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
- c=sqrt((2.0*e)/(3.0*moldyn->count*K_BOLTZMANN*moldyn->t));
- for(i=0;i<moldyn->count;i++)
- v3_scale(&(atom[i].v),&(atom[i].v),(1.0/c));
+ if(atom[i].attr&ATOM_ATTR_HB)
+ v3_scale(&(atom[i].v),&(atom[i].v),scale);
return 0;
}
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);
count2=27;
a=nx*ny;
-
cell[0]=lc->subcell[i+j*nx+k*a];
for(ci=-1;ci<=1;ci++) {
bx=0;
}
}
- lc->dnlc=count2;
+ lc->dnlc=count1;
lc->countn=27;
return count2;
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;
unsigned int e,m,s,v;
t_3dvec p;
t_moldyn_schedule *schedule;
-
+ t_atom *atom;
int fd;
char fb[128];
+ double ds;
schedule=&(moldyn->schedule);
+ atom=moldyn->atom;
/* initialize linked cell method */
link_cell_init(moldyn);
/* sqaure of some variables */
moldyn->tau_square=moldyn->tau*moldyn->tau;
moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff;
-
/* calculate initial forces */
potential_force_calc(moldyn);
+ /* do some checks before we actually start calculating bullshit */
+ if(moldyn->cutoff>0.5*moldyn->dim.x)
+ printf("[moldyn] warning: cutoff > 0.5 x dim.x\n");
+ if(moldyn->cutoff>0.5*moldyn->dim.y)
+ printf("[moldyn] warning: cutoff > 0.5 x dim.y\n");
+ if(moldyn->cutoff>0.5*moldyn->dim.z)
+ printf("[moldyn] warning: cutoff > 0.5 x dim.z\n");
+ ds=0.5*atom[0].f.x*moldyn->tau_square/atom[0].mass;
+ if(ds>0.05*moldyn->nnd)
+ printf("[moldyn] warning: forces too high / tau too small!\n");
+
/* zero absolute time */
moldyn->time=0.0;
-
for(sched=0;sched<moldyn->schedule.content_count;sched++) {
/* setting amount of runs and finite time step size */
}
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 */
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);
- v3_per_bound(&(atom[i].r),&(moldyn->dim));
+ check_per_bound(moldyn,&(atom[i].r));
/* velocities */
v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
}
/* 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 */
int potential_force_calc(t_moldyn *moldyn) {
int i,j,k,count;
- t_atom *atom,*btom,*ktom;
+ t_atom *itom,*jtom,*ktom;
t_linkcell *lc;
- t_list neighbour[27];
- t_list *this,*thisk,*neighbourk;
- u8 bc,bck;
+ t_list neighbour_i[27],neighbour_j[27];
+ t_list *this,*that;
+ u8 bc_ij,bc_ijk;
int countn,dnlc;
count=moldyn->count;
- atom=moldyn->atom;
+ itom=moldyn->atom;
lc=&(moldyn->lc);
/* reset energy */
moldyn->energy=0.0;
for(i=0;i<count;i++) {
-
+
/* reset force */
- v3_zero(&(atom[i].f));
+ v3_zero(&(itom[i].f));
/* single particle potential/force */
- if(atom[i].attr&ATOM_ATTR_1BP)
- moldyn->func1b(moldyn,&(atom[i]));
+ if(itom[i].attr&ATOM_ATTR_1BP)
+ moldyn->func1b(moldyn,&(itom[i]));
/* 2 body pair potential/force */
- if(atom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)) {
-
+ if(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)) {
+
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,
- (atom[i].r.z+moldyn->dim.z/2)/lc->z,
- neighbour);
+ (itom[i].r.x+moldyn->dim.x/2)/lc->x,
+ (itom[i].r.y+moldyn->dim.y/2)/lc->y,
+ (itom[i].r.z+moldyn->dim.z/2)/lc->z,
+ neighbour_i);
countn=lc->countn;
dnlc=lc->dnlc;
for(j=0;j<countn;j++) {
- this=&(neighbour[j]);
+ this=&(neighbour_i[j]);
list_reset(this);
if(this->start==NULL)
continue;
- bc=(j<dnlc)?0:1;
+ bc_ij=(j<dnlc)?0:1;
do {
- btom=this->current->data;
+ jtom=this->current->data;
- if(btom==&(atom[i]))
+ if(jtom==&(itom[i]))
continue;
- if((btom->attr&ATOM_ATTR_2BP)&
- (atom[i].attr&ATOM_ATTR_2BP))
+ if((jtom->attr&ATOM_ATTR_2BP)&
+ (itom[i].attr&ATOM_ATTR_2BP))
moldyn->func2b(moldyn,
- &(atom[i]),
- btom,
- bc);
+ &(itom[i]),
+ jtom,
+ bc_ij);
/* 3 body potential/force */
- if(!(atom[i].attr&ATOM_ATTR_3BP)||
- !(btom->attr&ATOM_ATTR_3BP))
+ if(!(itom[i].attr&ATOM_ATTR_3BP)||
+ !(jtom->attr&ATOM_ATTR_3BP))
continue;
link_cell_neighbour_index(moldyn,
- (btom->r.x+moldyn->dim.x/2)/lc->x,
- (btom->r.y+moldyn->dim.y/2)/lc->y,
- (btom->r.z+moldyn->dim.z/2)/lc->z,
- neighbourk);
+ (jtom->r.x+moldyn->dim.x/2)/lc->x,
+ (jtom->r.y+moldyn->dim.y/2)/lc->y,
+ (jtom->r.z+moldyn->dim.z/2)/lc->z,
+ neighbour_j);
+ /* neighbours of j */
for(k=0;k<lc->countn;k++) {
- thisk=&(neighbourk[k]);
- list_reset(thisk);
+ that=&(neighbour_j[k]);
+ list_reset(that);
+
+ if(that->start==NULL)
+ continue;
+
+ bc_ijk=(k<lc->dnlc)?0:1;
+
+ do {
+
+ ktom=that->current->data;
+
+ if(!(ktom->attr&ATOM_ATTR_3BP))
+ continue;
+
+ if(ktom==jtom)
+ continue;
+
+ if(ktom==&(itom[i]))
+ continue;
+
+ moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ijk);
+
+ } while(list_next(that)!=\
+ L_NO_NEXT_ELEMENT);
+
+ }
- if(thisk->start==NULL)
+ /* neighbours of i */
+ for(k=0;k<countn;k++) {
+
+ that=&(neighbour_i[k]);
+ list_reset(that);
+
+ if(that->start==NULL)
continue;
- bck=(k<lc->dnlc)?0:1;
+ bc_ijk=(k<dnlc)?0:1;
do {
- ktom=thisk->current->data;
+ ktom=that->current->data;
if(!(ktom->attr&ATOM_ATTR_3BP))
continue;
- if(ktom==btom)
+ if(ktom==jtom)
continue;
- if(ktom==&(atom[i]))
+ if(ktom==&(itom[i]))
continue;
- moldyn->func3b(moldyn,&(atom[i]),btom,ktom,bck);
+ moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ijk);
- } while(list_next(thisk)!=\
+ } while(list_next(that)!=\
L_NO_NEXT_ELEMENT);
}
v3_sub(&distance,&(ai->r),&(aj->r));
- v3_per_bound(&distance,&(moldyn->dim));
if(bc) check_per_bound(moldyn,&distance);
d=v3_norm(&distance);
if(d<=moldyn->cutoff) {
d=+h1-h2;
d*=eps;
v3_scale(&force,&distance,d);
- v3_add(&(ai->f),&(aj->f),&force);
+ v3_add(&(ai->f),&(ai->f),&force);
}
return 0;
* tersoff potential & force for 2 sorts of atoms
*/
+/* create mixed terms from parameters and set them */
+int tersoff_mult_complete_params(t_tersoff_mult_params *p) {
+
+ printf("[moldyn] tersoff parameter completion\n");
+ p->Smixed=sqrt(p->S[0]*p->S[1]);
+ p->Rmixed=sqrt(p->R[0]*p->R[1]);
+ p->Amixed=sqrt(p->A[0]*p->A[1]);
+ p->Bmixed=sqrt(p->B[0]*p->B[1]);
+ p->lambda_m=0.5*(p->lambda[0]+p->lambda[1]);
+ p->mu_m=0.5*(p->mu[0]+p->mu[1]);
+
+ return 0;
+}
+
/* tersoff 1 body part */
int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) {