+/*
+ * mdrun.c - main code to run a md simulation
+ *
+ * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
+ *
+ */
+
+#include "mdrun.h"
+
+/* potential */
+#include "potentials/harmonic_oscillator.h"
+#include "potentials/lennard_jones.h"
+#include "potentials/albe.h"
+#ifdef TERSOFF_ORIG
+#include "potentials/tersoff_orig.h"
+#else
+#include "potentials/tersoff.h"
+#endif
+
+#define ME "[mdrun]"
+
+/*
+ * parse function
+ */
+
+int mdrun_usage(void) {
+
+ printf("%s usage:\n",ME);
+
+ return 0;
+}
+
+int mdrun_parse_argv(t_mdrun *mdrun,int argc,char **argv) {
+
+ int i;
+
+ for(i=1;i<argc;i++) {
+
+ if(argv[i][0]!='-') {
+ printf("%s unknown switch: %s\n",ME,argv[i]);
+ return -1;
+ }
+
+ switch(argv[i][1]) {
+ case 'c':
+ strncpy(mdrun->cfile,argv[++i],128);
+ break;
+ case 's':
+ strncpy(mdrun->sdir,argv[++i],128);
+ break;
+ case 'h':
+ mdrun_usage();
+ break;
+ default:
+ printf("%s unknown option: %s\n",ME,argv[i]);
+ mdrun_usage();
+ return -1;
+ }
+
+ }
+
+ return 0;
+}
+
+int mdrun_parse_config(t_mdrun *mdrun) {
+
+ int fd,ret;
+ char error[128];
+ char line[128];
+ char *wptr;
+ char word[16][32];
+ int wcnt;
+
+ /* open config file */
+ fd=open(mdrun->cfile,O_RDONLY);
+ if(fd<0) {
+ snprintf(error,128,"%s open cfile %s",ME,mdrun->cfile);
+ perror(error);
+ return fd;
+ }
+
+ /* read, parse, set */
+ ret=1;
+ while(ret>0) {
+
+ /* read */
+ ret=get_line(fd,line,128);
+
+ /* parse */
+
+ // ignore # lines and \n
+ if((line[0]=='#')|(ret==1))
+ continue;
+
+ // get command + args
+ wcnt=0;
+ while(1) {
+ if(wcnt)
+ wptr=strtok(NULL," \t");
+ else
+ wptr=strtok(line," \t");
+ if(wptr==NULL)
+ break;
+ strncpy(word[wcnt],wptr,32);
+ wcnt+=1;
+ }
+
+ if(!strncmp(word[0],"potential",9)) {
+ if(!strncmp(word[1],"albe",4))
+ mdrun->potential=MOLDYN_POTENTIAL_AM;
+ if(!strncmp(word[1],"tersoff",7))
+ mdrun->potential=MOLDYN_POTENTIAL_TM;
+ if(!strncmp(word[1],"ho",2))
+ mdrun->potential=MOLDYN_POTENTIAL_HO;
+ if(!strncmp(word[1],"lj",2))
+ mdrun->potential=MOLDYN_POTENTIAL_LJ;
+ }
+ else if(!strncmp(word[0],"cutoff",6))
+ mdrun->cutoff=atof(word[1]);
+ else if(!strncmp(word[0],"intalgo",7)) {
+ if(!strncmp(word[1],"verlet",5))
+ mdrun->intalgo=MOLDYN_INTEGRATE_VERLET;
+ }
+ else if(!strncmp(word[0],"timestep",8))
+ mdrun->timestep=atof(word[1]);
+ else if(!strncmp(word[0],"volume",6)) {
+ mdrun->dim.x=atof(word[1]);
+ mdrun->dim.y=atof(word[2]);
+ mdrun->dim.z=atof(word[3]);
+ if(strncmp(word[4],"0",1))
+ mdrun->vis=TRUE;
+ }
+ else if(!strncmp(word[0],"pbc",3)) {
+ if(strncmp(word[1],"0",1))
+ mdrun->pbcx=TRUE;
+ else
+ mdrun->pbcx=FALSE;
+ if(strncmp(word[2],"0",1))
+ mdrun->pbcy=TRUE;
+ else
+ mdrun->pbcy=FALSE;
+ if(strncmp(word[3],"0",1))
+ mdrun->pbcz=TRUE;
+ else
+ mdrun->pbcz=FALSE;
+ }
+ else if(!strncmp(word[0],"temperature",11))
+ mdrun->temperature=atof(word[1]);
+ else if(!strncmp(word[0],"pressure",8))
+ mdrun->pressure=atof(word[1]);
+ else if(!strncmp(word[0],"lattice",7)) {
+ if(!strncmp(word[1],"zincblende",10))
+ mdrun->lattice=ZINCBLENDE;
+ if(!strncmp(word[1],"fcc",3))
+ mdrun->lattice=FCC;
+ if(!strncmp(word[1],"diamond",7))
+ mdrun->lattice=DIAMOND;
+ }
+ else if(!strncmp(word[0],"element1",8)) {
+ mdrun->element1=atoi(word[1]);
+ switch(mdrun->element1) {
+ case SI:
+ mdrun->m1=M_SI;
+ break;
+ case C:
+ mdrun->m1=M_C;
+ break;
+ default:
+ printf("%s unknown element1: %s|%d\n",
+ ME,word[1],mdrun->element1);
+ return -1;
+ }
+ }
+ else if(!strncmp(word[0],"element2",8)) {
+ mdrun->element2=atoi(word[1]);
+ switch(mdrun->element2) {
+ case SI:
+ mdrun->m2=M_SI;
+ break;
+ case C:
+ mdrun->m2=M_C;
+ break;
+ default:
+ printf("%s unknown element2: %s|%d\n",
+ ME,word[1],mdrun->element2);
+ return -1;
+ }
+ }
+ else if(!strncmp(word[0],"filllc",6)) {
+ mdrun->lx=atoi(word[1]);
+ mdrun->ly=atoi(word[2]);
+ mdrun->lz=atoi(word[3]);
+ }
+ else if(!strncmp(word[0],"aattrib",7)) {
+ // for aatrib line we need a special stage
+ // containing one schedule of 0 loops ...
+ printf("%s aattrib: %s\n",ME,word[1]);
+ }
+ else if(!strncmp(word[0],"eattrib",7)) {
+ mdrun->p_tau=atof(word[1]);
+ mdrun->t_tau=atof(word[2]);
+ }
+ else if(!strncmp(word[0],"prerun",6))
+ mdrun->prerun=atoi(word[1]);
+ else if(!strncmp(word[0],"avgskip",7))
+ mdrun->avgskip=atoi(word[1]);
+ else if(!strncmp(word[0],"elog",4))
+ mdrun->elog=atoi(word[1]);
+ else if(!strncmp(word[0],"tlog",4))
+ mdrun->tlog=atoi(word[1]);
+ else if(!strncmp(word[0],"plog",4))
+ mdrun->plog=atoi(word[1]);
+ else if(!strncmp(word[0],"vlog",4))
+ mdrun->vlog=atoi(word[1]);
+ else if(!strncmp(word[0],"save",4))
+ mdrun->save=atoi(word[1]);
+ else if(!strncmp(word[0],"visualize",9))
+ mdrun->visualize=atoi(word[1]);
+ else if(!strncmp(word[0],"stage",5)) {
+ // for every stage line, add a stage
+ printf("%s stage: %s\n",ME,word[1]);
+ }
+ else {
+ printf("%s unknown command %s, skipped!\n",ME,wptr);
+ continue;
+ }
+ }
+
+ return 0;
+}
+
+int main(int argc,char **argv) {
+
+ t_mdrun mdrun;
+ t_moldyn moldyn;
+ t_3dvec o;
+
+ /* clear mdrun struct */
+ memset(&mdrun,0,sizeof(t_mdrun));
+
+ /* parse arguments */
+ if(mdrun_parse_argv(&mdrun,argc,argv)<0)
+ return -1;
+
+ /* parse config file */
+ mdrun_parse_config(&mdrun);
+
+ /* sanity check! */
+
+ /* prepare simulation */
+ moldyn_init(&moldyn,argc,argv);
+ if(set_int_alg(&moldyn,mdrun.intalgo)<0)
+ return -1;
+ set_cutoff(&moldyn,mdrun.cutoff);
+ if(set_potential(&moldyn,mdrun.potential)<0)
+ return -1;
+ switch(mdrun.potential) {
+ case MOLDYN_POTENTIAL_AM:
+ albe_mult_set_params(&moldyn,
+ mdrun.element1,
+ mdrun.element2);
+ break;
+ case MOLDYN_POTENTIAL_TM:
+ tersoff_mult_set_params(&moldyn,
+ mdrun.element1,
+ mdrun.element2);
+ break;
+ case MOLDYN_POTENTIAL_HO:
+ harmonic_oscillator_set_params(&moldyn,mdrun.element1);
+ break;
+ case MOLDYN_POTENTIAL_LJ:
+ lennard_jones_set_params(&moldyn,mdrun.element1);
+ break;
+ default:
+ printf("%s unknown potential: %02x\n",
+ ME,mdrun.potential);
+ return -1;
+ }
+ 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) {
+ case FCC:
+ create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element1,
+ mdrun.m1,0,0,mdrun.lx,
+ mdrun.ly,mdrun.lz,NULL);
+ break;
+ case DIAMOND:
+ create_lattice(&moldyn,DIAMOND,mdrun.lc,mdrun.element1,
+ mdrun.m1,0,0,mdrun.lx,
+ mdrun.ly,mdrun.lz,NULL);
+ break;
+ case ZINCBLENDE:
+ o.x=0.5*0.25*mdrun.lc; o.y=o.x; o.z=o.x;
+ create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element1,
+ mdrun.m1,0,0,mdrun.lx,
+ mdrun.ly,mdrun.lz,&o);
+ o.x+=0.25*mdrun.lc; o.y=o.x; o.z=o.x;
+ create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element2,
+ mdrun.m2,0,0,mdrun.lx,
+ mdrun.ly,mdrun.lz,&o);
+ break;
+ default:
+ printf("%s unknown lattice type: %02x\n",
+ ME,mdrun.lattice);
+ return -1;
+ }
+ moldyn_bc_check(&moldyn);
+ set_temperature(&moldyn,mdrun.temperature);
+ set_pressure(&moldyn,mdrun.pressure);
+ thermal_init(&moldyn,TRUE);
+ moldyn_add_schedule(&moldyn,mdrun.prerun,mdrun.timestep);
+ moldyn_set_log_dir(&moldyn,mdrun.sdir);
+ moldyn_set_report(&moldyn,"CHANGE ME","CHANGE ME TOO");
+ if(mdrun.elog)
+ moldyn_set_log(&moldyn,LOG_TOTAL_ENERGY,mdrun.elog);
+ if(mdrun.tlog)
+ moldyn_set_log(&moldyn,LOG_TEMPERATURE,mdrun.tlog);
+ if(mdrun.plog)
+ moldyn_set_log(&moldyn,LOG_PRESSURE,mdrun.plog);
+ if(mdrun.vlog)
+ moldyn_set_log(&moldyn,LOG_VOLUME,mdrun.vlog);
+ if(mdrun.visualize)
+ moldyn_set_log(&moldyn,VISUAL_STEP,mdrun.visualize);
+ if(mdrun.save)
+ moldyn_set_log(&moldyn,SAVE_STEP,mdrun.save);
+ moldyn_set_log(&moldyn,CREATE_REPORT,0);
+ set_avg_skip(&moldyn,mdrun.avgskip);
+
+ /* prepare the hook function */
+
+ /* run the simulation */
+ moldyn_integrate(&moldyn);
+
+ /* shutdown */
+ moldyn_shutdown(&moldyn);
+
+ return 0;
+}