]> hackdaworld.org Git - physik/posic.git/blob - moldyn.c
7f9745ad29ef2f8e7dd52a52301a4f0a05ad0283
[physik/posic.git] / moldyn.c
1 /*
2  * moldyn.c - molecular dynamics library main file
3  *
4  * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
5  *
6  */
7
8 #define _GNU_SOURCE
9 #include <stdio.h>
10 #include <stdlib.h>
11 #include <string.h>
12 #include <sys/types.h>
13 #include <sys/stat.h>
14 #include <fcntl.h>
15 #include <unistd.h>
16 #include <math.h>
17
18 #include "moldyn.h"
19
20 #include "math/math.h"
21 #include "init/init.h"
22 #include "random/random.h"
23 #include "visual/visual.h"
24
25 int moldyn_usage(char **argv) {
26
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("-V <steps> <filebase> (rasmol file)\n");
34         printf("--- physics options ---\n");
35         printf("-T <temperature> [K] (%f)\n",MOLDYN_TEMP);
36         printf("-t <timestep tau> [s] (%.15f)\n",MOLDYN_TAU);
37         printf("-R <runs> (%d)\n",MOLDYN_RUNS);
38         printf("\n");
39
40         return 0;
41 }
42
43 int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv) {
44
45         int i;
46
47         memset(moldyn,0,sizeof(t_moldyn));
48
49         /* default values */
50         moldyn->t=MOLDYN_TEMP;
51         moldyn->tau=MOLDYN_TAU;
52         moldyn->time_steps=MOLDYN_RUNS;
53
54         /* parse argv */
55         for(i=1;i<argc;i++) {
56                 if(argv[i][0]=='-') {
57                         switch(argv[i][1]){
58                                 case 'E':
59                                         moldyn->ewrite=atoi(argv[++i]);
60                                         strncpy(moldyn->efb,argv[++i],64);
61                                         break;
62                                 case 'M':
63                                         moldyn->mwrite=atoi(argv[++i]);
64                                         strncpy(moldyn->mfb,argv[++i],64);
65                                         break;
66                                 case 'D':
67                                         moldyn->dwrite=atoi(argv[++i]);
68                                         strncpy(moldyn->dfb,argv[++i],64);
69                                         break;
70                                 case 'S':
71                                         moldyn->swrite=atoi(argv[++i]);
72                                         strncpy(moldyn->sfb,argv[++i],64);
73                                         break;
74                                 case 'V':
75                                         moldyn->vwrite=atoi(argv[++i]);
76                                         strncpy(moldyn->vfb,argv[++i],64);
77                                         break;
78                                 case 'T':
79                                         moldyn->t=atof(argv[++i]);
80                                         break;
81                                 case 't':
82                                         moldyn->tau=atof(argv[++i]);
83                                         break;
84                                 case 'R':
85                                         moldyn->time_steps=atoi(argv[++i]);
86                                         break;
87                                 default:
88                                         printf("unknown option %s\n",argv[i]);
89                                         moldyn_usage(argv);
90                                         return -1;
91                         }
92                 } else {
93                         moldyn_usage(argv);
94                         return -1;
95                 }
96         }
97
98         return 0;
99 }
100
101 int moldyn_log_init(t_moldyn *moldyn,void *v) {
102
103         moldyn->lvstat=0;
104         t_visual *vis;
105
106         vis=v;
107
108         if(moldyn->ewrite) {
109                 moldyn->efd=open(moldyn->efb,O_WRONLY|O_CREAT|O_TRUNC);
110                 if(moldyn->efd<0) {
111                         perror("[moldyn] efd open");
112                         return moldyn->efd;
113                 }
114                 dprintf(moldyn->efd,"# moldyn total energy logfile\n");
115                 moldyn->lvstat|=MOLDYN_LVSTAT_TOTAL_E;
116         }
117
118         if(moldyn->mwrite) {
119                 moldyn->mfd=open(moldyn->mfb,O_WRONLY|O_CREAT|O_TRUNC);
120                 if(moldyn->mfd<0) {
121                         perror("[moldyn] mfd open");
122                         return moldyn->mfd;
123                 }
124                 dprintf(moldyn->mfd,"# moldyn total momentum logfile\n");
125                 moldyn->lvstat|=MOLDYN_LVSTAT_TOTAL_M;
126         }
127
128         if(moldyn->swrite)
129                 moldyn->lvstat|=MOLDYN_LVSTAT_SAVE;
130
131         if(moldyn->dwrite) {
132                 moldyn->dfd=open(moldyn->dfb,O_WRONLY|O_CREAT|O_TRUNC);
133                 if(moldyn->dfd<0) {
134                         perror("[moldyn] dfd open");
135                         return moldyn->dfd;
136                 }
137                 write(moldyn->dfd,moldyn,sizeof(t_moldyn));
138                 moldyn->lvstat|=MOLDYN_LVSTAT_DUMP;
139         }
140
141         if((moldyn->vwrite)&&(vis)) {
142                 moldyn->visual=vis;
143                 visual_init(vis,moldyn->vfb);
144                 moldyn->lvstat|=MOLDYN_LVSTAT_VISUAL;
145         }
146
147         moldyn->lvstat|=MOLDYN_LVSTAT_INITIALIZED;
148
149         return 0;
150 }
151
152 int moldyn_shutdown(t_moldyn *moldyn) {
153
154         if(moldyn->efd) close(moldyn->efd);
155         if(moldyn->mfd) close(moldyn->efd);
156         if(moldyn->dfd) close(moldyn->efd);
157         if(moldyn->visual) visual_tini(moldyn->visual);
158
159         return 0;
160 }
161
162 int create_lattice(unsigned char type,int element,double mass,double lc,
163                    int a,int b,int c,t_atom **atom) {
164
165         int count;
166         int ret;
167         t_3dvec origin;
168
169         count=a*b*c;
170
171         if(type==FCC) count*=4;
172         if(type==DIAMOND) count*=8;
173
174         *atom=malloc(count*sizeof(t_atom));
175         if(*atom==NULL) {
176                 perror("malloc (atoms)");
177                 return -1;
178         }
179
180         v3_zero(&origin);
181
182         switch(type) {
183                 case FCC:
184                         ret=fcc_init(a,b,c,lc,*atom,&origin);
185                         break;
186                 case DIAMOND:
187                         ret=diamond_init(a,b,c,lc,*atom,&origin);
188                         break;
189                 default:
190                         printf("unknown lattice type (%02x)\n",type);
191                         return -1;
192         }
193
194         /* debug */
195         if(ret!=count) {
196                 printf("ok, there is something wrong ...\n");
197                 printf("calculated -> %d atoms\n",count);
198                 printf("created -> %d atoms\n",ret);
199                 return -1;
200         }
201
202         while(count) {
203                 (*atom)[count-1].element=element;
204                 (*atom)[count-1].mass=mass;
205                 count-=1;
206         }
207
208         return ret;
209 }
210
211 int destroy_lattice(t_atom *atom) {
212
213         if(atom) free(atom);
214
215         return 0;
216 }
217
218 int thermal_init(t_moldyn *moldyn,t_random *random,int count) {
219
220         /*
221          * - gaussian distribution of velocities
222          * - zero total momentum
223          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
224          */
225
226         int i;
227         double v,sigma;
228         t_3dvec p_total,delta;
229         t_atom *atom;
230
231         atom=moldyn->atom;
232
233         /* gaussian distribution of velocities */
234         v3_zero(&p_total);
235         for(i=0;i<count;i++) {
236                 sigma=sqrt(2.0*K_BOLTZMANN*moldyn->t/atom[i].mass);
237                 /* x direction */
238                 v=sigma*rand_get_gauss(random);
239                 atom[i].v.x=v;
240                 p_total.x+=atom[i].mass*v;
241                 /* y direction */
242                 v=sigma*rand_get_gauss(random);
243                 atom[i].v.y=v;
244                 p_total.y+=atom[i].mass*v;
245                 /* z direction */
246                 v=sigma*rand_get_gauss(random);
247                 atom[i].v.z=v;
248                 p_total.z+=atom[i].mass*v;
249         }
250
251         /* zero total momentum */
252         v3_scale(&p_total,&p_total,1.0/count);
253         for(i=0;i<count;i++) {
254                 v3_scale(&delta,&p_total,1.0/atom[i].mass);
255                 v3_sub(&(atom[i].v),&(atom[i].v),&delta);
256         }
257
258         /* velocity scaling */
259         scale_velocity(moldyn,count);
260
261         return 0;
262 }
263
264 int scale_velocity(t_moldyn *moldyn,int count) {
265
266         int i;
267         double e,c;
268         t_atom *atom;
269
270         atom=moldyn->atom;
271
272         /*
273          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
274          */
275         e=0.0;
276         for(i=0;i<count;i++)
277                 e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
278         c=sqrt((2.0*e)/(3.0*count*K_BOLTZMANN*moldyn->t));
279         for(i=0;i<count;i++)
280                 v3_scale(&(atom[i].v),&(atom[i].v),(1.0/c));
281
282         return 0;
283 }
284
285 double get_e_kin(t_atom *atom,int count) {
286
287         int i;
288         double e;
289
290         e=0.0;
291
292         for(i=0;i<count;i++) {
293                 e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
294         }
295
296         return e;
297 }
298
299 double get_e_pot(t_moldyn *moldyn) {
300
301         return(moldyn->potential(moldyn));
302 }
303
304 double get_total_energy(t_moldyn *moldyn) {
305
306         double e;
307
308         e=get_e_kin(moldyn->atom,moldyn->count);
309         e+=get_e_pot(moldyn);
310
311         return e;
312 }
313
314 t_3dvec get_total_p(t_atom *atom, int count) {
315
316         t_3dvec p,p_total;
317         int i;
318
319         v3_zero(&p_total);
320         for(i=0;i<count;i++) {
321                 v3_scale(&p,&(atom[i].v),atom[i].mass);
322                 v3_add(&p_total,&p_total,&p);
323         }
324
325         return p_total;
326 }
327
328 double estimate_time_step(t_moldyn *moldyn,double nn_dist,double t) {
329
330         double tau;
331
332         tau=0.05*nn_dist/(sqrt(3.0*K_BOLTZMANN*t/moldyn->atom[0].mass));
333         tau*=1.0E-9;
334         if(tau<moldyn->tau)
335                 printf("[moldyn] warning: time step  (%f > %.15f)\n",
336                        moldyn->tau,tau);
337
338         return tau;     
339 }
340
341
342 /*
343  *
344  * 'integration of newtons equation' - algorithms
345  *
346  */
347
348 /* start the integration */
349
350 int moldyn_integrate(t_moldyn *moldyn) {
351
352         int i;
353         unsigned int e,m,s,d,v;
354         t_3dvec p;
355
356         int fd;
357         char fb[128];
358
359         e=moldyn->ewrite;
360         m=moldyn->mwrite;
361         s=moldyn->swrite;
362         d=moldyn->dwrite;
363         v=moldyn->vwrite;
364
365         if(!(moldyn->lvstat&MOLDYN_LVSTAT_INITIALIZED)) {
366                 printf("[moldyn] warning, lv system not initialized\n");
367                 return -1;
368         }
369
370         /* calculate initial forces */
371         moldyn->force(moldyn);
372
373         for(i=0;i<moldyn->time_steps;i++) {
374                 /* integration step */
375                 moldyn->integrate(moldyn);
376
377                 /* check for log & visualiziation */
378                 if(e) {
379                         if(!(i%e))
380                                 dprintf(moldyn->efd,
381                                         "%.15f %.45f\n",i*moldyn->tau,
382                                         get_total_energy(moldyn));
383                 }
384                 if(m) {
385                         if(!(i%m)) {
386                                 p=get_total_p(moldyn->atom,moldyn->count);
387                                 dprintf(moldyn->mfd,
388                                         "%.15f %.45f\n",i*moldyn->tau,
389                                         v3_norm(&p));
390                         }
391                 }
392                 if(s) {
393                         if(!(i%s)) {
394                                 snprintf(fb,128,"%s-%f-%.15f.save",moldyn->sfb,
395                                          moldyn->t,i*moldyn->tau);
396                                 fd=open(fb,O_WRONLY|O_TRUNC|O_CREAT);
397                                 if(fd<0) perror("[moldyn] save fd open");
398                                 else {
399                                         write(fd,moldyn,sizeof(t_moldyn));
400                                         write(fd,moldyn->atom,
401                                               moldyn->count*sizeof(t_atom));
402                                 }
403                         }       
404                 }
405                 if(d) {
406                         if(!(i%d))
407                                 write(moldyn->dfd,moldyn->atom,
408                                       moldyn->count*sizeof(t_atom));
409
410                 }
411                 if(v) {
412                         if(!(i%v))
413                                 visual_atoms(moldyn->visual,i*moldyn->tau,
414                                              moldyn->atom,moldyn->count);
415                 }
416         }
417
418         return 0;
419 }
420
421 /* velocity verlet */
422
423 int velocity_verlet(t_moldyn *moldyn) {
424
425         int i,count;
426         double tau,tau_square;
427         t_3dvec delta;
428         t_atom *atom;
429
430         atom=moldyn->atom;
431         count=moldyn->count;
432         tau=moldyn->tau;
433
434         tau_square=tau*tau;
435
436         for(i=0;i<count;i++) {
437                 /* new positions */
438                 v3_scale(&delta,&(atom[i].v),tau);
439                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
440                 v3_scale(&delta,&(atom[i].f),0.5*tau_square/atom[i].mass);
441                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
442                 v3_per_bound(&(atom[i].r),&(moldyn->dim));
443
444                 /* velocities */
445                 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
446                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
447         }
448
449         /* forces depending on chosen potential */
450         moldyn->force(moldyn);
451
452         for(i=0;i<count;i++) {
453                 /* again velocities */
454                 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
455                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
456         }
457
458         return 0;
459 }
460
461
462 /*
463  *
464  * potentials & corresponding forces
465  * 
466  */
467
468 /* harmonic oscillator potential and force */
469
470 double potential_harmonic_oscillator(t_moldyn *moldyn) {
471
472         t_ho_params *params;
473         t_atom *atom;
474         int i,j;
475         int count;
476         t_3dvec distance;
477         double d,u;
478         double sc,equi_dist;
479
480         params=moldyn->pot_params;
481         atom=moldyn->atom;
482         sc=params->spring_constant;
483         equi_dist=params->equilibrium_distance;
484         count=moldyn->count;
485
486         u=0.0;
487         for(i=0;i<count;i++) {
488                 for(j=0;j<i;j++) {
489                         v3_sub(&distance,&(atom[i].r),&(atom[j].r));
490                         d=v3_norm(&distance);
491                         u+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
492                 }
493         }
494
495         return u;
496 }
497
498 int force_harmonic_oscillator(t_moldyn *moldyn) {
499
500         t_ho_params *params;
501         int i,j,count;
502         t_atom *atom;
503         t_3dvec distance;
504         t_3dvec force;
505         double d;
506         double sc,equi_dist;
507
508         atom=moldyn->atom;      
509         count=moldyn->count;
510         params=moldyn->pot_params;
511         sc=params->spring_constant;
512         equi_dist=params->equilibrium_distance;
513
514         for(i=0;i<count;i++) v3_zero(&(atom[i].f));
515
516         for(i=0;i<count;i++) {
517                 for(j=0;j<i;j++) {
518                         v3_sub(&distance,&(atom[i].r),&(atom[j].r));
519                         v3_per_bound(&distance,&(moldyn->dim));
520                         d=v3_norm(&distance);
521                         if(d<=moldyn->cutoff) {
522                                 v3_scale(&force,&distance,
523                                          -sc*(1.0-(equi_dist/d)));
524                                 v3_add(&(atom[i].f),&(atom[i].f),&force);
525                                 v3_sub(&(atom[j].f),&(atom[j].f),&force);
526                         }
527                 }
528         }
529
530         return 0;
531 }
532
533
534 /* lennard jones potential & force for one sort of atoms */
535  
536 double potential_lennard_jones(t_moldyn *moldyn) {
537
538         t_lj_params *params;
539         t_atom *atom;
540         int i,j;
541         int count;
542         t_3dvec distance;
543         double d,help;
544         double u;
545         double eps,sig6,sig12;
546
547         params=moldyn->pot_params;
548         atom=moldyn->atom;
549         count=moldyn->count;
550         eps=params->epsilon4;
551         sig6=params->sigma6;
552         sig12=params->sigma12;
553
554         u=0.0;
555         for(i=0;i<count;i++) {
556                 for(j=0;j<i;j++) {
557                         v3_sub(&distance,&(atom[i].r),&(atom[j].r));
558                         d=1.0/v3_absolute_square(&distance);    /* 1/r^2 */
559                         help=d*d;                               /* 1/r^4 */
560                         help*=d;                                /* 1/r^6 */
561                         d=help*help;                            /* 1/r^12 */
562                         u+=eps*(sig6*help-sig12*d);
563                 }
564         }
565         
566         return u;
567 }
568
569 int force_lennard_jones(t_moldyn *moldyn) {
570
571         t_lj_params *params;
572         int i,j,count;
573         t_atom *atom;
574         t_3dvec distance;
575         t_3dvec force;
576         double d,h1,h2;
577         double eps,sig6,sig12;
578
579         atom=moldyn->atom;      
580         count=moldyn->count;
581         params=moldyn->pot_params;
582         eps=params->epsilon4;
583         sig6=6*params->sigma6;
584         sig12=12*params->sigma12;
585
586         for(i=0;i<count;i++) v3_zero(&(atom[i].f));
587
588         for(i=0;i<count;i++) {
589                 for(j=0;j<i;j++) {
590                         v3_sub(&distance,&(atom[i].r),&(atom[j].r));
591                         v3_per_bound(&distance,&(moldyn->dim));
592                         d=v3_absolute_square(&distance);
593                         if(d<=moldyn->cutoff_square) {
594                                 h1=1.0/d;                       /* 1/r^2 */
595                                 d=h1*h1;                        /* 1/r^4 */
596                                 h2=d*d;                         /* 1/r^8 */
597                                 h1*=d;                          /* 1/r^6 */
598                                 h1*=h2;                         /* 1/r^14 */
599                                 h1*=sig12;
600                                 h2*=sig6;
601                                 /* actually there would be a '-',       *
602                                  * but f=-d/dr potential                */
603                                 d=h1+h2;
604                                 d*=eps;
605                                 v3_scale(&force,&distance,d);
606                                 v3_add(&(atom[i].f),&(atom[i].f),&force);
607                                 v3_sub(&(atom[j].f),&(atom[j].f),&force);
608                         }
609                 }
610         }
611
612         return 0;
613 }
614