t_stage *stage;
switch(type) {
+ case STAGE_DISPLACE_ATOM:
+ psize=sizeof(t_displace_atom_params);
+ break;
case STAGE_INSERT_ATOMS:
psize=sizeof(t_insert_atoms_params);
break;
char error[128];
char line[128];
char *wptr;
- char word[16][32];
+ char word[16][64];
int wcnt;
int i,o;
+ t_displace_atom_params dap;
t_insert_atoms_params iap;
t_continue_params cp;
t_anneal_params ap;
wptr=strtok(line," \t");
if(wptr==NULL)
break;
- strncpy(word[wcnt],wptr,32);
+ strncpy(word[wcnt],wptr,64);
wcnt+=1;
}
if(!strncmp(word[1],"lj",2))
mdrun->potential=MOLDYN_POTENTIAL_LJ;
}
+ else if(!strncmp(word[0],"continue",8))
+ strncpy(mdrun->continue_file,word[1],128);
else if(!strncmp(word[0],"cutoff",6))
mdrun->cutoff=atof(word[1]);
else if(!strncmp(word[0],"nnd",3))
csp.type|=CHSATTR_TCTRL;
}
if(!strncmp(word[i],"prelax",6)) {
- csp.dp=atof(word[++i]);
+ csp.dp=atof(word[++i])*BAR;
csp.type|=CHSATTR_PRELAX;
}
if(!strncmp(word[i],"trelax",6)) {
mdrun->visualize=atoi(word[1]);
else if(!strncmp(word[0],"stage",5)) {
// for every stage line, add a stage
- if(!strncmp(word[1],"ins_atoms",9)) {
+ if(!strncmp(word[1],"displace",8)) {
+ dap.nr=atoi(word[2]);
+ dap.dx=atof(word[3]);
+ dap.dy=atof(word[4]);
+ dap.dz=atof(word[5]);
+ add_stage(mdrun,STAGE_DISPLACE_ATOM,&dap);
+ }
+ else if(!strncmp(word[1],"ins_atoms",9)) {
iap.ins_steps=atoi(word[2]);
iap.ins_atoms=atoi(word[3]);
iap.element=atoi(word[4]);
- iap.element=atoi(word[4]);
iap.brand=atoi(word[5]);
for(i=0;i<strlen(word[6]);i++) {
switch(word[6][i]) {
break;
}
}
- // only rand mode by now
- if(word[8][0]=='t')
- iap.type=INS_TOTAL;
- else {
- iap.type=INS_REGION;
- iap.x0=atof(word[8]);
- iap.y0=atof(word[9]);
- iap.z0=atof(word[10]);
- iap.x1=atof(word[11]);
- iap.y1=atof(word[12]);
- iap.z1=atof(word[13]);
+ switch(word[7][0]) {
+ // insertion types
+ case 'p':
+ iap.type=INS_POS;
+ iap.x0=atof(word[8]);
+ iap.y0=atof(word[9]);
+ iap.z0=atof(word[10]);
+ break;
+ case 'r':
+ if(word[8][0]=='t') {
+ iap.type=INS_TOTAL;
+ iap.cr=atof(word[9]);
+ }
+ else {
+ iap.type=INS_REGION;
+ iap.x0=atof(word[8]);
+ iap.y0=atof(word[9]);
+ iap.z0=atof(word[10]);
+ iap.x1=atof(word[11]);
+ iap.y1=atof(word[12]);
+ iap.z1=atof(word[13]);
+ iap.cr=atof(word[14]);
+ }
+ default:
+ break;
}
add_stage(mdrun,STAGE_INSERT_ATOMS,&iap);
}
double delta;
- if(!(mdrun->sattr&SATTR_PRELAX))
+ if(!(mdrun->sattr&SATTR_PRELAX)) {
return TRUE;
+ }
delta=moldyn->p_avg-moldyn->p_ref;
return TRUE;
}
+int displace_atom(t_moldyn *moldyn,t_mdrun *mdrun) {
+
+ t_displace_atom_params *dap;
+ t_stage *stage;
+ t_atom *atom;
+
+ stage=mdrun->stage.current->data;
+ dap=stage->params;
+
+ atom=&(moldyn->atom[dap->nr]);
+ atom->r.x+=dap->dx;
+ atom->r.y+=dap->dy;
+ atom->r.z+=dap->dz;
+
+ return 0;
+}
+
int insert_atoms(t_moldyn *moldyn,t_mdrun *mdrun) {
t_insert_atoms_params *iap;
t_stage *stage;
t_atom *atom;
t_3dvec r,v,dist;
- double d,dmin;
+ double d,dmin,o;
int cnt,i;
double x,y,z;
double x0,y0,z0;
cr_check=FALSE;
v.x=0.0; v.y=0.0; v.z=0.0;
+ x=0; y=0; z=0;
+ o=0; dmin=0;
+
+ switch(mdrun->lattice) {
+ case CUBIC:
+ o=0.5*mdrun->lc;
+ break;
+ case FCC:
+ o=0.25*mdrun->lc;
+ break;
+ case DIAMOND:
+ case ZINCBLENDE:
+ o=0.125*mdrun->lc;
+ break;
+ default:
+ break;
+ }
switch(iap->type) {
case INS_TOTAL:
x=moldyn->dim.x;
- x0=0.0;
+ x0=-x/2.0;
y=moldyn->dim.y;
- y0=0.0;
+ y0=-y/2.0;
z=moldyn->dim.z;
- z0=0.0;
+ z0=-z/2.0;
cr_check=TRUE;
break;
case INS_REGION:
z0=iap->z0;
cr_check=TRUE;
break;
+ case INS_POS:
+ x0=iap->x0;
+ y0=iap->y0;
+ z0=iap->z0;
+ cr_check=FALSE;
+ break;
default:
printf("%s unknown insertion mode: %02x\n",
ME,iap->type);
while(cnt<iap->ins_atoms) {
run=1;
while(run) {
- r.x=(rand_get_double(&(moldyn->random))-0.5)*x+x0;
- r.y=(rand_get_double(&(moldyn->random))-0.5)*y+y0;
- r.z=(rand_get_double(&(moldyn->random))-0.5)*z+z0;
+ if(iap->type!=INS_POS) {
+ r.x=rand_get_double(&(moldyn->random))*x;
+ r.y=rand_get_double(&(moldyn->random))*y;
+ r.z=rand_get_double(&(moldyn->random))*z;
+ }
+ else {
+ r.x=0.0;
+ r.y=0.0;
+ r.z=0.0;
+ }
+ r.x+=x0;
+ r.y+=y0;
+ r.z+=z0;
+ // offset
+ if(iap->type!=INS_TOTAL) {
+ r.x+=o;
+ r.y+=o;
+ r.z+=o;
+ }
run=0;
if(cr_check==TRUE) {
dmin=1000; // for sure too high!
}
add_atom(moldyn,iap->element,pse_mass[iap->element],
iap->brand,iap->attr,&r,&v);
- printf("%s atom inserted: %f %f %f | d squared = %f\n",
- ME,r.x,r.y,r.z,dmin);
+ printf("%s atom inserted (%d/%d): %f %f %f\n",
+ ME,(iap->cnt_steps+1)*iap->ins_atoms,
+ iap->ins_steps*iap->ins_atoms,r.x,r.y,r.z);
+ printf(" -> d2 = %f/%f\n",dmin,iap->cr*iap->cr);
cnt+=1;
}
}
if(csp->type&CHSATTR_AVGRST) {
if(csp->avgrst)
- mdrun->sattr|=CHSATTR_AVGRST;
+ mdrun->sattr|=SATTR_AVGRST;
else
- mdrun->sattr&=(~(CHSATTR_AVGRST));
+ mdrun->sattr&=(~(SATTR_AVGRST));
}
if(csp->type&CHSATTR_RSTEPS) {
mdrun->relax_steps=csp->rsteps;
tau=mdrun->timestep;
/* check whether relaxation steps are necessary */
- if(!((check_pressure(moldyn,mdrun)==FALSE)|\
- (check_temperature(moldyn,mdrun)==FALSE))) {
+ if((check_pressure(moldyn,mdrun)==TRUE)&\
+ (check_temperature(moldyn,mdrun)==TRUE)) {
/* be verbose */
stage_print("\n###########################\n");
/* stage specific stuff */
switch(stage->type) {
+ case STAGE_DISPLACE_ATOM:
+ stage_print(" -> displace atom\n\n");
+ displace_atom(moldyn,mdrun);
+ change_stage=TRUE;
+ break;
case STAGE_INSERT_ATOMS:
stage_print(" -> insert atoms\n\n");
iap=stage->params;
}
insert_atoms(moldyn,mdrun);
iap->cnt_steps+=1;
- break;
+ break;
case STAGE_CONTINUE:
stage_print(" -> continue\n\n");
if(stage->executed==TRUE) {
}
cp=stage->params;
steps=cp->runs;
- break;
+ break;
case STAGE_ANNEAL:
stage_print(" -> anneal\n\n");
ap=stage->params;
change_stage=TRUE;
break;
}
- set_temperature(moldyn,moldyn->t_ref+ap->dt);
+ if(moldyn->t_ref+ap->dt>=0.0)
+ set_temperature(moldyn,
+ moldyn->t_ref+ap->dt);
ap->count+=1;
break;
case STAGE_CHAATTR:
/* sanity check! */
/* prepare simulation */
+
moldyn_init(&moldyn,argc,argv);
+
+ if(mdrun.continue_file)
+ moldyn_read_save_file(&moldyn,mdrun.continue_file);
+
if(set_int_alg(&moldyn,mdrun.intalgo)<0)
return -1;
+
+ /* potential */
set_cutoff(&moldyn,mdrun.cutoff);
if(set_potential(&moldyn,mdrun.potential)<0)
return -1;
ME,mdrun.potential);
return -1;
}
+
+ /* if it is a continue run, reset schedule and skip lattice init */
+ if(mdrun.continue_file) {
+ memset(&(moldyn.schedule),0,sizeof(t_moldyn_schedule));
+ goto addsched;
+ }
+
+ /* initial lattice and dimensions */
set_dim(&moldyn,mdrun.dim.x,mdrun.dim.y,mdrun.dim.z,mdrun.vis);
set_pbc(&moldyn,mdrun.pbcx,mdrun.pbcy,mdrun.pbcz);
switch(mdrun.lattice) {
return -1;
}
moldyn_bc_check(&moldyn);
+
+ /* temperature and pressure */
set_temperature(&moldyn,mdrun.temperature);
set_pressure(&moldyn,mdrun.pressure);
thermal_init(&moldyn,TRUE);
+
+addsched:
+ /* first schedule */
moldyn_add_schedule(&moldyn,mdrun.prerun,mdrun.timestep);
+
+ /* log */
moldyn_set_log_dir(&moldyn,mdrun.sdir);
moldyn_set_report(&moldyn,"CHANGE ME","CHANGE ME TOO");
if(mdrun.elog)