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"
24 #include "list/list.h"
26 int moldyn_usage(char **argv) {
28 printf("\n%s usage:\n\n",argv[0]);
29 printf("--- general options ---\n");
30 printf("-E <steps> <file> (log total energy)\n");
31 printf("-M <steps> <file> (log total momentum)\n");
32 printf("-D <steps> <file> (dump total information)\n");
33 printf("-S <steps> <filebase> (single save file)\n");
34 printf("-V <steps> <filebase> (rasmol file)\n");
35 printf("--- physics options ---\n");
36 printf("-T <temperature> [K] (%f)\n",MOLDYN_TEMP);
37 printf("-t <timestep tau> [s] (%.15f)\n",MOLDYN_TAU);
38 printf("-R <runs> (%d)\n",MOLDYN_RUNS);
39 printf(" -- integration algo --\n");
40 printf(" -I <number> (%d)\n",MOLDYN_INTEGRATE_DEFAULT);
41 printf(" 0: velocity verlet\n");
42 printf(" -- potential --\n");
43 printf(" -P <number> <param1 param2 ...>\n");
44 printf(" 0: harmonic oscillator\n");
45 printf(" param1: spring constant\n");
46 printf(" param2: equilibrium distance\n");
47 printf(" 1: lennard jones\n");
48 printf(" param1: epsilon\n");
49 printf(" param2: sigma\n");
55 int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv) {
62 memset(moldyn,0,sizeof(t_moldyn));
65 moldyn->t=MOLDYN_TEMP;
66 moldyn->tau=MOLDYN_TAU;
67 moldyn->time_steps=MOLDYN_RUNS;
68 moldyn->integrate=velocity_verlet;
69 moldyn->potential_force_function=lennard_jones;
76 moldyn->ewrite=atoi(argv[++i]);
77 strncpy(moldyn->efb,argv[++i],64);
80 moldyn->mwrite=atoi(argv[++i]);
81 strncpy(moldyn->mfb,argv[++i],64);
84 moldyn->dwrite=atoi(argv[++i]);
85 strncpy(moldyn->dfb,argv[++i],64);
88 moldyn->swrite=atoi(argv[++i]);
89 strncpy(moldyn->sfb,argv[++i],64);
92 moldyn->vwrite=atoi(argv[++i]);
93 strncpy(moldyn->vfb,argv[++i],64);
96 moldyn->t=atof(argv[++i]);
99 moldyn->tau=atof(argv[++i]);
102 moldyn->time_steps=atoi(argv[++i]);
105 /* integration algorithm */
106 switch(atoi(argv[++i])) {
107 case MOLDYN_INTEGRATE_VERLET:
108 moldyn->integrate=velocity_verlet;
111 printf("unknown integration algo %s\n",argv[i]);
117 /* potential + params */
118 switch(atoi(argv[++i])) {
119 case MOLDYN_POTENTIAL_HO:
120 hop.spring_constant=atof(argv[++i]);
121 hop.equilibrium_distance=atof(argv[++i]);
122 moldyn->pot_params=malloc(sizeof(t_ho_params));
123 memcpy(moldyn->pot_params,&hop,sizeof(t_ho_params));
124 moldyn->potential_force_function=harmonic_oscillator;
126 case MOLDYN_POTENTIAL_LJ:
130 ljp.sigma6=s*s*s*s*s*s;
131 ljp.sigma12=ljp.sigma6*ljp.sigma6;
132 moldyn->pot_params=malloc(sizeof(t_lj_params));
133 memcpy(moldyn->pot_params,&ljp,sizeof(t_lj_params));
134 moldyn->potential_force_function=lennard_jones;
137 printf("unknown potential %s\n",argv[i]);
143 printf("unknown option %s\n",argv[i]);
156 int moldyn_log_init(t_moldyn *moldyn,void *v) {
164 moldyn->efd=open(moldyn->efb,O_WRONLY|O_CREAT|O_TRUNC);
166 perror("[moldyn] efd open");
169 dprintf(moldyn->efd,"# moldyn total energy logfile\n");
170 moldyn->lvstat|=MOLDYN_LVSTAT_TOTAL_E;
174 moldyn->mfd=open(moldyn->mfb,O_WRONLY|O_CREAT|O_TRUNC);
176 perror("[moldyn] mfd open");
179 dprintf(moldyn->mfd,"# moldyn total momentum logfile\n");
180 moldyn->lvstat|=MOLDYN_LVSTAT_TOTAL_M;
184 moldyn->lvstat|=MOLDYN_LVSTAT_SAVE;
187 moldyn->dfd=open(moldyn->dfb,O_WRONLY|O_CREAT|O_TRUNC);
189 perror("[moldyn] dfd open");
192 write(moldyn->dfd,moldyn,sizeof(t_moldyn));
193 moldyn->lvstat|=MOLDYN_LVSTAT_DUMP;
196 if((moldyn->vwrite)&&(vis)) {
198 visual_init(vis,moldyn->vfb);
199 moldyn->lvstat|=MOLDYN_LVSTAT_VISUAL;
202 moldyn->lvstat|=MOLDYN_LVSTAT_INITIALIZED;
207 int moldyn_shutdown(t_moldyn *moldyn) {
209 if(moldyn->efd) close(moldyn->efd);
210 if(moldyn->mfd) close(moldyn->efd);
211 if(moldyn->dfd) close(moldyn->efd);
212 if(moldyn->visual) visual_tini(moldyn->visual);
217 int create_lattice(unsigned char type,int element,double mass,double lc,
218 int a,int b,int c,t_atom **atom) {
226 if(type==FCC) count*=4;
227 if(type==DIAMOND) count*=8;
229 *atom=malloc(count*sizeof(t_atom));
231 perror("malloc (atoms)");
239 ret=fcc_init(a,b,c,lc,*atom,&origin);
242 ret=diamond_init(a,b,c,lc,*atom,&origin);
245 printf("unknown lattice type (%02x)\n",type);
251 printf("ok, there is something wrong ...\n");
252 printf("calculated -> %d atoms\n",count);
253 printf("created -> %d atoms\n",ret);
258 (*atom)[count-1].element=element;
259 (*atom)[count-1].mass=mass;
266 int destroy_lattice(t_atom *atom) {
273 int thermal_init(t_moldyn *moldyn,t_random *random) {
276 * - gaussian distribution of velocities
277 * - zero total momentum
278 * - velocity scaling (E = 3/2 N k T), E: kinetic energy
283 t_3dvec p_total,delta;
288 /* gaussian distribution of velocities */
290 for(i=0;i<moldyn->count;i++) {
291 sigma=sqrt(2.0*K_BOLTZMANN*moldyn->t/atom[i].mass);
293 v=sigma*rand_get_gauss(random);
295 p_total.x+=atom[i].mass*v;
297 v=sigma*rand_get_gauss(random);
299 p_total.y+=atom[i].mass*v;
301 v=sigma*rand_get_gauss(random);
303 p_total.z+=atom[i].mass*v;
306 /* zero total momentum */
307 v3_scale(&p_total,&p_total,1.0/moldyn->count);
308 for(i=0;i<moldyn->count;i++) {
309 v3_scale(&delta,&p_total,1.0/atom[i].mass);
310 v3_sub(&(atom[i].v),&(atom[i].v),&delta);
313 /* velocity scaling */
314 scale_velocity(moldyn);
319 int scale_velocity(t_moldyn *moldyn) {
328 * - velocity scaling (E = 3/2 N k T), E: kinetic energy
331 for(i=0;i<moldyn->count;i++)
332 e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
333 c=sqrt((2.0*e)/(3.0*moldyn->count*K_BOLTZMANN*moldyn->t));
334 for(i=0;i<moldyn->count;i++)
335 v3_scale(&(atom[i].v),&(atom[i].v),(1.0/c));
340 double get_e_kin(t_atom *atom,int count) {
347 for(i=0;i<count;i++) {
348 e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
354 double get_e_pot(t_moldyn *moldyn) {
356 return moldyn->energy;
359 double get_total_energy(t_moldyn *moldyn) {
363 e=get_e_kin(moldyn->atom,moldyn->count);
364 e+=get_e_pot(moldyn);
369 t_3dvec get_total_p(t_atom *atom, int count) {
375 for(i=0;i<count;i++) {
376 v3_scale(&p,&(atom[i].v),atom[i].mass);
377 v3_add(&p_total,&p_total,&p);
383 double estimate_time_step(t_moldyn *moldyn,double nn_dist,double t) {
387 tau=0.05*nn_dist/(sqrt(3.0*K_BOLTZMANN*t/moldyn->atom[0].mass));
390 printf("[moldyn] warning: time step (%f > %.15f)\n",
400 /* linked list / cell method */
402 int link_cell_init(t_moldyn *moldyn) {
408 /* partitioning the md cell */
409 lc->nx=moldyn->dim.x/moldyn->cutoff;
410 lc->x=moldyn->dim.x/lc->nx;
411 lc->ny=moldyn->dim.y/moldyn->cutoff;
412 lc->y=moldyn->dim.y/lc->ny;
413 lc->nz=moldyn->dim.z/moldyn->cutoff;
414 lc->z=moldyn->dim.z/lc->nz;
416 lc->subcell=malloc(lc->nx*lc->ny*lc->nz*sizeof(t_list));
418 link_cell_update(moldyn);
423 int link_cell_update(t_moldyn *moldyn) {
437 for(i=0;i<nx*ny*nz;i++)
438 list_destroy(&(moldyn->lc.subcell[i]));
440 for(count=0;count<moldyn->count;count++) {
441 i=atom[count].r.x/lc->x;
442 j=atom[count].r.y/lc->y;
443 k=atom[count].r.z/lc->z;
444 list_add_immediate_ptr(&(moldyn->lc.subcell[i+j*nx+k*nx*ny]),
451 int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) {
459 unsigned char bx,by,bz;
469 cell[0]=lc->subcell[i+j*nx+k*a];
470 for(ci=-1;ci<=1;ci++) {
477 for(cj=-1;cj<=1;cj++) {
484 for(ck=-1;ck<=1;ck++) {
491 if(!(x|y|z)) continue;
493 cell[--count2]=lc->subcell[x+y*nx+z*a];
496 cell[count1++]=lc->subcell[x+y*nx+z*a];
505 int link_cell_shutdown(t_moldyn *moldyn) {
512 for(i=0;i<lc->nx*lc->ny*lc->nz;i++)
513 list_shutdown(&(moldyn->lc.subcell[i]));
520 * 'integration of newtons equation' - algorithms
524 /* start the integration */
526 int moldyn_integrate(t_moldyn *moldyn) {
529 unsigned int e,m,s,d,v;
535 /* logging & visualization */
542 if(!(moldyn->lvstat&MOLDYN_LVSTAT_INITIALIZED)) {
543 printf("[moldyn] warning, lv system not initialized\n");
547 /* sqaure of some variables */
548 moldyn->tau_square=moldyn->tau*moldyn->tau;
549 moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff;
551 /* create the neighbour list */
552 link_cell_update(moldyn);
554 /* calculate initial forces */
555 moldyn->potential_force_function(moldyn);
557 for(i=0;i<moldyn->time_steps;i++) {
561 /* neighbour list update */
562 link_cell_update(moldyn);
564 /* integration step */
565 moldyn->integrate(moldyn);
567 /* check for log & visualization */
571 "%.15f %.45f\n",i*moldyn->tau,
572 get_total_energy(moldyn));
576 p=get_total_p(moldyn->atom,moldyn->count);
578 "%.15f %.45f\n",i*moldyn->tau,
584 snprintf(fb,128,"%s-%f-%.15f.save",moldyn->sfb,
585 moldyn->t,i*moldyn->tau);
586 fd=open(fb,O_WRONLY|O_TRUNC|O_CREAT);
587 if(fd<0) perror("[moldyn] save fd open");
589 write(fd,moldyn,sizeof(t_moldyn));
590 write(fd,moldyn->atom,
591 moldyn->count*sizeof(t_atom));
597 write(moldyn->dfd,moldyn->atom,
598 moldyn->count*sizeof(t_atom));
603 visual_atoms(moldyn->visual,i*moldyn->tau,
604 moldyn->atom,moldyn->count);
611 /* velocity verlet */
613 int velocity_verlet(t_moldyn *moldyn) {
616 double tau,tau_square;
623 tau_square=moldyn->tau_square;
625 for(i=0;i<count;i++) {
627 v3_scale(&delta,&(atom[i].v),tau);
628 v3_add(&(atom[i].r),&(atom[i].r),&delta);
629 v3_scale(&delta,&(atom[i].f),0.5*tau_square/atom[i].mass);
630 v3_add(&(atom[i].r),&(atom[i].r),&delta);
631 v3_per_bound(&(atom[i].r),&(moldyn->dim));
634 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
635 v3_add(&(atom[i].v),&(atom[i].v),&delta);
638 /* forces depending on chosen potential */
639 moldyn->potential_force_function(moldyn);
641 for(i=0;i<count;i++) {
642 /* again velocities */
643 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
644 v3_add(&(atom[i].v),&(atom[i].v),&delta);
653 * potentials & corresponding forces
657 /* harmonic oscillator potential and force */
659 int harmonic_oscillator(t_moldyn *moldyn) {
664 t_list *this,neighbour[27];
667 t_3dvec force,distance;
672 params=moldyn->pot_params;
675 sc=params->spring_constant;
676 equi_dist=params->equilibrium_distance;
680 for(i=0;i<count;i++) {
681 /* determine cell + neighbours */
682 ni=atom[i].r.x/lc->x;
683 nj=atom[i].r.y/lc->y;
684 nk=atom[i].r.z/lc->z;
685 c=link_cell_neighbour_index(moldyn,ni,nj,nk,neighbour);
687 /* processing cell of atom i */
688 this=&(neighbour[0]);
689 list_reset(this); /* list has 1 element at minimum */
691 btom=this->current->data;
694 v3_sub(&distance,&(atom[i].r),&(btom->r));
695 d=v3_norm(&distance);
696 u+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
697 v3_scale(&force,&distance,-sc*(1.0-(equi_dist/d)));
698 v3_add(&(atom[i].f),&(atom[i].f),&force);
699 } while(list_next(this)!=L_NO_NEXT_ELEMENT);
701 /* neighbours not doing boundary condition overflow */
703 this=&(neighbour[j]);
704 list_reset(this); /* there might not be a single atom */
705 if(this->start!=NULL) {
708 btom=this->current->data;
709 v3_sub(&distance,&(atom[i].r),&(btom->r));
710 d=v3_norm(&distance);
711 if(d<=moldyn->cutoff) {
712 u+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
713 v3_scale(&force,&distance,
714 -sc*(1.0-(equi_dist/d)));
715 v3_add(&(atom[i].f),&(atom[i].f),
718 } while(list_next(this)!=L_NO_NEXT_ELEMENT);
723 /* neighbours due to boundary conditions */
725 this=&(neighbour[j]);
726 list_reset(this); /* check boundary conditions */
727 if(this->start!=NULL) {
730 btom=this->current->data;
731 v3_sub(&distance,&(atom[i].r),&(btom->r));
732 v3_per_bound(&distance,&(moldyn->dim));
733 d=v3_norm(&distance);
734 if(d<=moldyn->cutoff) {
735 u+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
736 v3_scale(&force,&distance,
737 -sc*(1.0-(equi_dist/d)));
738 v3_add(&(atom[i].f),&(atom[i].f),
741 } while(list_next(this)!=L_NO_NEXT_ELEMENT);
752 /* lennard jones potential & force for one sort of atoms */
754 int lennard_jones(t_moldyn *moldyn) {
759 t_list *this,neighbour[27];
762 t_3dvec force,distance;
764 double eps,sig6,sig12;
767 params=moldyn->pot_params;
771 eps=params->epsilon4;
773 sig12=params->sigma12;
776 for(i=0;i<count;i++) {
777 /* determine cell + neighbours */
778 ni=atom[i].r.x/lc->x;
779 nj=atom[i].r.y/lc->y;
780 nk=atom[i].r.z/lc->z;
781 c=link_cell_neighbour_index(moldyn,ni,nj,nk,neighbour);
783 /* processing cell of atom i */
784 this=&(neighbour[0]);
785 list_reset(this); /* list has 1 element at minimum */
787 btom=this->current->data;
790 v3_sub(&distance,&(atom[i].r),&(btom->r));
791 d=1.0/v3_absolute_square(&distance); /* 1/r^2 */
794 h1=h2*h2; /* 1/r^12 */
795 u+=eps*(sig12*h1-sig6*h2);
802 v3_scale(&force,&distance,d);
803 v3_add(&(atom[i].f),&(atom[i].f),&force);
804 } while(list_next(this)!=L_NO_NEXT_ELEMENT);
806 /* neighbours not doing boundary condition overflow */
808 this=&(neighbour[j]);
809 list_reset(this); /* there might not be a single atom */
810 if(this->start!=NULL) {
813 btom=this->current->data;
814 v3_sub(&distance,&(atom[i].r),&(btom->r));
815 d=v3_absolute_square(&distance); /* r^2 */
816 if(d<=moldyn->cutoff_square) {
820 h1=h2*h2; /* 1/r^12 */
821 u+=eps*(sig12*h1-sig6*h2);
828 v3_scale(&force,&distance,d);
829 v3_add(&(atom[i].f),&(atom[i].f),
832 } while(list_next(this)!=L_NO_NEXT_ELEMENT);
837 /* neighbours due to boundary conditions */
839 this=&(neighbour[j]);
840 list_reset(this); /* check boundary conditions */
841 if(this->start!=NULL) {
844 btom=this->current->data;
845 v3_sub(&distance,&(atom[i].r),&(btom->r));
846 v3_per_bound(&distance,&(moldyn->dim));
847 d=v3_absolute_square(&distance); /* r^2 */
848 if(d<=moldyn->cutoff_square) {
852 h1=h2*h2; /* 1/r^12 */
853 u+=eps*(sig12*h1-sig6*h2);
860 v3_scale(&force,&distance,d);
861 v3_add(&(atom[i].f),&(atom[i].f),
864 } while(list_next(this)!=L_NO_NEXT_ELEMENT);