2 * mdrun.c - main code to run a md simulation
4 * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
11 #include "potentials/harmonic_oscillator.h"
12 #include "potentials/lennard_jones.h"
13 #include "potentials/albe.h"
15 #include "potentials/tersoff_orig.h"
17 #include "potentials/tersoff.h"
26 int mdrun_usage(void) {
28 printf("%s usage:\n",ME);
33 int mdrun_parse_argv(t_mdrun *mdrun,int argc,char **argv) {
40 printf("%s unknown switch: %s\n",ME,argv[i]);
46 strncpy(mdrun->cfile,argv[++i],128);
49 strncpy(mdrun->sdir,argv[++i],128);
55 printf("%s unknown option: %s\n",ME,argv[i]);
65 int mdrun_parse_config(t_mdrun *mdrun) {
74 /* open config file */
75 fd=open(mdrun->cfile,O_RDONLY);
77 snprintf(error,128,"%s open cfile %s",ME,mdrun->cfile);
82 /* read, parse, set */
87 ret=get_line(fd,line,128);
91 // ignore # lines and \n
92 if((line[0]=='#')|(ret==1))
99 wptr=strtok(NULL," \t");
101 wptr=strtok(line," \t");
104 strncpy(word[wcnt],wptr,32);
108 if(!strncmp(word[0],"potential",9)) {
109 if(!strncmp(word[1],"albe",4))
110 mdrun->potential=MOLDYN_POTENTIAL_AM;
111 if(!strncmp(word[1],"tersoff",7))
112 mdrun->potential=MOLDYN_POTENTIAL_TM;
113 if(!strncmp(word[1],"ho",2))
114 mdrun->potential=MOLDYN_POTENTIAL_HO;
115 if(!strncmp(word[1],"lj",2))
116 mdrun->potential=MOLDYN_POTENTIAL_LJ;
118 else if(!strncmp(word[0],"cutoff",6))
119 mdrun->cutoff=atof(word[1]);
120 else if(!strncmp(word[0],"intalgo",7)) {
121 if(!strncmp(word[1],"verlet",5))
122 mdrun->intalgo=MOLDYN_INTEGRATE_VERLET;
124 else if(!strncmp(word[0],"timestep",8))
125 mdrun->timestep=atof(word[1]);
126 else if(!strncmp(word[0],"volume",6)) {
127 mdrun->dim.x=atof(word[1]);
128 mdrun->dim.y=atof(word[2]);
129 mdrun->dim.z=atof(word[3]);
130 if(strncmp(word[4],"0",1))
133 else if(!strncmp(word[0],"pbc",3)) {
134 if(strncmp(word[1],"0",1))
138 if(strncmp(word[2],"0",1))
142 if(strncmp(word[3],"0",1))
147 else if(!strncmp(word[0],"temperature",11))
148 mdrun->temperature=atof(word[1]);
149 else if(!strncmp(word[0],"pressure",8))
150 mdrun->pressure=atof(word[1]);
151 else if(!strncmp(word[0],"lattice",7)) {
152 if(!strncmp(word[1],"zincblende",10))
153 mdrun->lattice=ZINCBLENDE;
154 if(!strncmp(word[1],"fcc",3))
156 if(!strncmp(word[1],"diamond",7))
157 mdrun->lattice=DIAMOND;
159 else if(!strncmp(word[0],"element1",8)) {
160 mdrun->element1=atoi(word[1]);
161 switch(mdrun->element1) {
169 printf("%s unknown element1: %s|%d\n",
170 ME,word[1],mdrun->element1);
174 else if(!strncmp(word[0],"element2",8)) {
175 mdrun->element2=atoi(word[1]);
176 switch(mdrun->element2) {
184 printf("%s unknown element2: %s|%d\n",
185 ME,word[1],mdrun->element2);
189 else if(!strncmp(word[0],"filllc",6)) {
190 mdrun->lx=atoi(word[1]);
191 mdrun->ly=atoi(word[2]);
192 mdrun->lz=atoi(word[3]);
194 else if(!strncmp(word[0],"aattrib",7)) {
195 // for aatrib line we need a special stage
196 // containing one schedule of 0 loops ...
197 printf("%s aattrib: %s\n",ME,word[1]);
199 else if(!strncmp(word[0],"eattrib",7)) {
200 mdrun->p_tau=atof(word[1]);
201 mdrun->t_tau=atof(word[2]);
203 else if(!strncmp(word[0],"prerun",6))
204 mdrun->prerun=atoi(word[1]);
205 else if(!strncmp(word[0],"avgskip",7))
206 mdrun->avgskip=atoi(word[1]);
207 else if(!strncmp(word[0],"elog",4))
208 mdrun->elog=atoi(word[1]);
209 else if(!strncmp(word[0],"tlog",4))
210 mdrun->tlog=atoi(word[1]);
211 else if(!strncmp(word[0],"plog",4))
212 mdrun->plog=atoi(word[1]);
213 else if(!strncmp(word[0],"vlog",4))
214 mdrun->vlog=atoi(word[1]);
215 else if(!strncmp(word[0],"save",4))
216 mdrun->save=atoi(word[1]);
217 else if(!strncmp(word[0],"visualize",9))
218 mdrun->visualize=atoi(word[1]);
219 else if(!strncmp(word[0],"stage",5)) {
220 // for every stage line, add a stage
221 printf("%s stage: %s\n",ME,word[1]);
224 printf("%s unknown command %s, skipped!\n",ME,wptr);
232 int main(int argc,char **argv) {
238 /* clear mdrun struct */
239 memset(&mdrun,0,sizeof(t_mdrun));
241 /* parse arguments */
242 if(mdrun_parse_argv(&mdrun,argc,argv)<0)
245 /* parse config file */
246 mdrun_parse_config(&mdrun);
250 /* prepare simulation */
251 moldyn_init(&moldyn,argc,argv);
252 if(set_int_alg(&moldyn,mdrun.intalgo)<0)
254 set_cutoff(&moldyn,mdrun.cutoff);
255 if(set_potential(&moldyn,mdrun.potential)<0)
257 switch(mdrun.potential) {
258 case MOLDYN_POTENTIAL_AM:
259 albe_mult_set_params(&moldyn,
263 case MOLDYN_POTENTIAL_TM:
264 tersoff_mult_set_params(&moldyn,
268 case MOLDYN_POTENTIAL_HO:
269 harmonic_oscillator_set_params(&moldyn,mdrun.element1);
271 case MOLDYN_POTENTIAL_LJ:
272 lennard_jones_set_params(&moldyn,mdrun.element1);
275 printf("%s unknown potential: %02x\n",
279 set_dim(&moldyn,mdrun.dim.x,mdrun.dim.y,mdrun.dim.z,mdrun.vis);
280 set_pbc(&moldyn,mdrun.pbcx,mdrun.pbcy,mdrun.pbcz);
281 switch(mdrun.lattice) {
283 create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element1,
284 mdrun.m1,0,0,mdrun.lx,
285 mdrun.ly,mdrun.lz,NULL);
288 create_lattice(&moldyn,DIAMOND,mdrun.lc,mdrun.element1,
289 mdrun.m1,0,0,mdrun.lx,
290 mdrun.ly,mdrun.lz,NULL);
293 o.x=0.5*0.25*mdrun.lc; o.y=o.x; o.z=o.x;
294 create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element1,
295 mdrun.m1,0,0,mdrun.lx,
296 mdrun.ly,mdrun.lz,&o);
297 o.x+=0.25*mdrun.lc; o.y=o.x; o.z=o.x;
298 create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element2,
299 mdrun.m2,0,0,mdrun.lx,
300 mdrun.ly,mdrun.lz,&o);
303 printf("%s unknown lattice type: %02x\n",
307 moldyn_bc_check(&moldyn);
308 set_temperature(&moldyn,mdrun.temperature);
309 set_pressure(&moldyn,mdrun.pressure);
310 thermal_init(&moldyn,TRUE);
311 moldyn_add_schedule(&moldyn,mdrun.prerun,mdrun.timestep);
312 moldyn_set_log_dir(&moldyn,mdrun.sdir);
313 moldyn_set_report(&moldyn,"CHANGE ME","CHANGE ME TOO");
315 moldyn_set_log(&moldyn,LOG_TOTAL_ENERGY,mdrun.elog);
317 moldyn_set_log(&moldyn,LOG_TEMPERATURE,mdrun.tlog);
319 moldyn_set_log(&moldyn,LOG_PRESSURE,mdrun.plog);
321 moldyn_set_log(&moldyn,LOG_VOLUME,mdrun.vlog);
323 moldyn_set_log(&moldyn,VISUAL_STEP,mdrun.visualize);
325 moldyn_set_log(&moldyn,SAVE_STEP,mdrun.save);
326 moldyn_set_log(&moldyn,CREATE_REPORT,0);
327 set_avg_skip(&moldyn,mdrun.avgskip);
329 /* prepare the hook function */
331 /* run the simulation */
332 moldyn_integrate(&moldyn);
335 moldyn_shutdown(&moldyn);