2 * moldyn.c - molecular dynamics library main file
4 * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
12 #include <sys/types.h>
20 #include "math/math.h"
21 #include "init/init.h"
22 #include "random/random.h"
23 #include "visual/visual.h"
25 int moldyn_usage(char **argv) {
27 printf("\n%s usage:\n\n",argv[0]);
28 printf("--- general options ---\n");
29 printf("-E <steps> <file> (log total energy)\n");
30 printf("-M <steps> <file> (log total momentum)\n");
31 printf("-D <steps> <file> (dump total information)\n");
32 printf("-S <steps> <filebase> (single save file)\n");
33 printf("--- physics options ---\n");
34 printf("-T <temperature> [K] (%f)\n",MOLDYN_TEMP);
35 printf("-t <timestep tau> [s] (%f)\n",MOLDYN_TAU);
36 printf("-R <runs> (%d)\n",MOLDYN_RUNS);
42 int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv) {
46 memset(moldyn,0,sizeof(t_moldyn));
49 moldyn->t=MOLDYN_TEMP;
50 moldyn->tau=MOLDYN_TAU;
51 moldyn->time_steps=MOLDYN_RUNS;
58 moldyn->ewrite=atoi(argv[++i]);
59 strncpy(moldyn->efb,argv[++i],64);
62 moldyn->mwrite=atoi(argv[++i]);
63 strncpy(moldyn->mfb,argv[++i],64);
66 moldyn->dwrite=atoi(argv[++i]);
67 strncpy(moldyn->dfb,argv[++i],64);
70 moldyn->swrite=atoi(argv[++i]);
71 strncpy(moldyn->sfb,argv[++i],64);
75 moldyn->t=atof(argv[++i]);
77 moldyn->tau=atof(argv[++i]);
80 moldyn->time_steps=atoi(argv[++i]);
83 printf("unknown option %s\n",argv[i]);
96 int moldyn_log_init(t_moldyn *moldyn) {
101 moldyn->efd=open(moldyn->efb,O_WRONLY|O_CREAT|O_TRUNC);
103 perror("[moldyn] efd open");
106 dprintf(moldyn->efd,"# moldyn total energy logfile\n");
107 moldyn->lvstat|=MOLDYN_LVSTAT_TOTAL_E;
111 moldyn->mfd=open(moldyn->mfb,O_WRONLY|O_CREAT|O_TRUNC);
113 perror("[moldyn] mfd open");
116 dprintf(moldyn->mfd,"# moldyn total momentum logfile\n");
117 moldyn->lvstat|=MOLDYN_LVSTAT_TOTAL_M;
121 moldyn->lvstat|=MOLDYN_LVSTAT_SAVE;
124 moldyn->dfd=open(moldyn->dfb,O_WRONLY|O_CREAT|O_TRUNC);
126 perror("[moldyn] dfd open");
129 write(moldyn->dfd,moldyn,sizeof(t_moldyn));
130 moldyn->lvstat|=MOLDYN_LVSTAT_DUMP;
134 moldyn->lvstat|=MOLDYN_LVSTAT_VISUAL;
139 int create_lattice(unsigned char type,int element,double mass,double lc,
140 int a,int b,int c,t_atom **atom) {
148 if(type==FCC) count*=4;
149 if(type==DIAMOND) count*=8;
151 *atom=malloc(count*sizeof(t_atom));
153 perror("malloc (atoms)");
161 ret=fcc_init(a,b,c,lc,*atom,&origin);
164 ret=diamond_init(a,b,c,lc,*atom,&origin);
167 printf("unknown lattice type (%02x)\n",type);
173 printf("ok, there is something wrong ...\n");
174 printf("calculated -> %d atoms\n",count);
175 printf("created -> %d atoms\n",ret);
180 (*atom)[count-1].element=element;
181 (*atom)[count-1].mass=mass;
188 int destroy_lattice(t_atom *atom) {
195 int thermal_init(t_moldyn *moldyn,t_random *random,int count) {
198 * - gaussian distribution of velocities
199 * - zero total momentum
200 * - velocity scaling (E = 3/2 N k T), E: kinetic energy
205 t_3dvec p_total,delta;
210 /* gaussian distribution of velocities */
212 for(i=0;i<count;i++) {
213 sigma=sqrt(2.0*K_BOLTZMANN*moldyn->t/atom[i].mass);
215 v=sigma*rand_get_gauss(random);
217 p_total.x+=atom[i].mass*v;
219 v=sigma*rand_get_gauss(random);
221 p_total.y+=atom[i].mass*v;
223 v=sigma*rand_get_gauss(random);
225 p_total.z+=atom[i].mass*v;
228 /* zero total momentum */
229 v3_scale(&p_total,&p_total,1.0/count);
230 for(i=0;i<count;i++) {
231 v3_scale(&delta,&p_total,1.0/atom[i].mass);
232 v3_sub(&(atom[i].v),&(atom[i].v),&delta);
235 /* velocity scaling */
236 scale_velocity(moldyn,count);
241 int scale_velocity(t_moldyn *moldyn,int count) {
250 * - velocity scaling (E = 3/2 N k T), E: kinetic energy
254 e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
255 c=sqrt((2.0*e)/(3.0*count*K_BOLTZMANN*moldyn->t));
257 v3_scale(&(atom[i].v),&(atom[i].v),(1.0/c));
262 double get_e_kin(t_atom *atom,int count) {
269 for(i=0;i<count;i++) {
270 e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
276 double get_e_pot(t_moldyn *moldyn) {
278 return(moldyn->potential(moldyn));
281 double get_total_energy(t_moldyn *moldyn) {
285 e=get_e_kin(moldyn->atom,moldyn->count);
286 e+=get_e_pot(moldyn);
291 t_3dvec get_total_p(t_atom *atom, int count) {
297 for(i=0;i<count;i++) {
298 v3_scale(&p,&(atom[i].v),atom[i].mass);
299 v3_add(&p_total,&p_total,&p);
305 double estimate_time_step(t_moldyn *moldyn,double nn_dist,double t) {
309 tau=0.05*nn_dist/(sqrt(3.0*K_BOLTZMANN*t/moldyn->atom[0].mass));
312 printf("[moldyn] warning: time step (%f > %.15f)\n",
321 * 'integration of newtons equation' - algorithms
325 /* start the integration */
327 int moldyn_integrate(t_moldyn *moldyn) {
330 unsigned int e,m,s,d,v;
331 unsigned char lvstat;
343 if(!(lvstat&MOLDYN_LVSTAT_INITIALIZED)) {
344 printf("[moldyn] warning, lv system not initialized\n");
348 /* calculate initial forces */
349 moldyn->force(moldyn);
351 for(i=0;i<moldyn->time_steps;i++) {
352 /* integration step */
353 moldyn->integrate(moldyn);
355 /* check for log & visualiziation */
359 "%.15f %.45f\n",i*moldyn->tau,
360 get_total_energy(moldyn));
364 p=get_total_p(moldyn->atom,moldyn->count);
366 "%.15f %.45f\n",i*moldyn->tau,
372 snprintf(fb,128,"%s-%f-%.15f.save",moldyn->sfb,
373 moldyn->t,i*moldyn->tau);
374 fd=open(fb,O_WRONLY|O_TRUNC|O_CREAT);
375 if(fd<0) perror("[moldyn] save fd open");
377 write(fd,moldyn,sizeof(t_moldyn));
378 write(fd,moldyn->atom,
379 moldyn->count*sizeof(t_atom));
385 write(moldyn->dfd,moldyn->atom,
386 moldyn->count*sizeof(t_atom));
391 visual_atoms(moldyn->visual,i*moldyn->tau,
392 moldyn->atom,moldyn->count);
399 /* velocity verlet */
401 int velocity_verlet(t_moldyn *moldyn) {
404 double tau,tau_square;
414 for(i=0;i<count;i++) {
416 v3_scale(&delta,&(atom[i].v),tau);
417 v3_add(&(atom[i].r),&(atom[i].r),&delta);
418 v3_scale(&delta,&(atom[i].f),0.5*tau_square/atom[i].mass);
419 v3_add(&(atom[i].r),&(atom[i].r),&delta);
420 v3_per_bound(&(atom[i].r),&(moldyn->dim));
423 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
424 v3_add(&(atom[i].v),&(atom[i].v),&delta);
427 /* forces depending on chosen potential */
428 moldyn->force(moldyn);
430 for(i=0;i<count;i++) {
431 /* again velocities */
432 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
433 v3_add(&(atom[i].v),&(atom[i].v),&delta);
442 * potentials & corresponding forces
446 /* harmonic oscillator potential and force */
448 double potential_harmonic_oscillator(t_moldyn *moldyn) {
458 params=moldyn->pot_params;
460 sc=params->spring_constant;
461 equi_dist=params->equilibrium_distance;
465 for(i=0;i<count;i++) {
467 v3_sub(&distance,&(atom[i].r),&(atom[j].r));
468 d=v3_norm(&distance);
469 u+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
476 int force_harmonic_oscillator(t_moldyn *moldyn) {
488 params=moldyn->pot_params;
489 sc=params->spring_constant;
490 equi_dist=params->equilibrium_distance;
492 for(i=0;i<count;i++) v3_zero(&(atom[i].f));
494 for(i=0;i<count;i++) {
496 v3_sub(&distance,&(atom[i].r),&(atom[j].r));
497 v3_per_bound(&distance,&(moldyn->dim));
498 d=v3_norm(&distance);
499 if(d<=moldyn->cutoff) {
500 v3_scale(&force,&distance,
501 (-sc*(1.0-(equi_dist/d))));
502 v3_add(&(atom[i].f),&(atom[i].f),&force);
503 v3_sub(&(atom[j].f),&(atom[j].f),&force);
512 /* lennard jones potential & force for one sort of atoms */
514 double potential_lennard_jones(t_moldyn *moldyn) {
523 double eps,sig6,sig12;
525 params=moldyn->pot_params;
530 sig12=params->sigma12;
533 for(i=0;i<count;i++) {
535 v3_sub(&distance,&(atom[j].r),&(atom[i].r));
536 d=1.0/v3_absolute_square(&distance); /* 1/r^2 */
537 help=d*d; /* 1/r^4 */
539 d=help*help; /* 1/r^12 */
540 u+=eps*(sig12*d-sig6*help);
547 int force_lennard_jones(t_moldyn *moldyn) {
555 double eps,sig6,sig12;
559 params=moldyn->pot_params;
562 sig12=params->sigma12;
564 for(i=0;i<count;i++) v3_zero(&(atom[i].f));
566 for(i=0;i<count;i++) {
568 v3_sub(&distance,&(atom[j].r),&(atom[i].r));
569 v3_per_bound(&distance,&(moldyn->dim));
570 d=v3_absolute_square(&distance);
571 if(d<=moldyn->cutoff_square) {
572 h1=1.0/d; /* 1/r^2 */
581 v3_scale(&force,&distance,d);
582 v3_add(&(atom[j].f),&(atom[j].f),&force);
583 v3_sub(&(atom[i].f),&(atom[i].f),&force);