int moldyn_shutdown(t_moldyn *moldyn) {
+ printf("[moldyn] shutdown\n");
moldyn_log_shutdown(moldyn);
link_cell_shutdown(moldyn);
- moldyn_log_shutdown(moldyn);
rand_close(&(moldyn->random));
free(moldyn->atom);
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;
while(count) {
- atom[count-1].element=element;
- atom[count-1].mass=mass;
- atom[count-1].attr=attr;
- atom[count-1].bnum=bnum;
+ moldyn->atom[count-1].element=element;
+ moldyn->atom[count-1].mass=mass;
+ moldyn->atom[count-1].attr=attr;
+ moldyn->atom[count-1].bnum=bnum;
count-=1;
}
moldyn->atom=ptr;
atom=moldyn->atom;
- atom->r=*r;
- atom->v=*v;
- atom->element=element;
- atom->bnum=bnum;
- atom->attr=attr;
+ 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 (E = 3/2 N k T), E: kinetic energy
*/
+
+ if(moldyn->t==0.0) {
+ printf("[moldyn] no velocity scaling for T = 0 K\n");
+ return -1;
+ }
+
e=0.0;
for(i=0;i<moldyn->count;i++)
e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
unsigned int e,m,s,v;
t_3dvec p;
t_moldyn_schedule *schedule;
+ t_atom *atom;
int fd;
char fb[128];
schedule=&(moldyn->schedule);
+ atom=moldyn->atom;
/* initialize linked cell method */
link_cell_init(moldyn);
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);
/* 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,
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;