]> hackdaworld.org Git - physik/posic.git/blob - moldyn.c
1b52a5628b3e3fe4df32ec4a76483795a1491c15
[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 #include "list/list.h"
25
26 int moldyn_usage(char **argv) {
27
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("-C <cutoff radius> [m] (%.15f)\n",MOLDYN_CUTOFF);
39         printf("-R <runs> (%d)\n",MOLDYN_RUNS);
40         printf(" -- integration algo --\n");
41         printf("  -I <number> (%d)\n",MOLDYN_INTEGRATE_DEFAULT);
42         printf("     0: velocity verlet\n");
43         printf(" -- potential --\n");
44         printf("  -P <number> <param1 param2 ...>\n");
45         printf("     0: harmonic oscillator\n");
46         printf("        param1: spring constant\n");
47         printf("        param2: equilibrium distance\n");
48         printf("     1: lennard jones\n");
49         printf("        param1: epsilon\n");
50         printf("        param2: sigma\n");
51         printf("\n");
52
53         return 0;
54 }
55
56 int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv) {
57
58         int i;
59         t_ho_params hop;
60         t_lj_params ljp;
61         double s,e;
62
63         memset(moldyn,0,sizeof(t_moldyn));
64
65         /* default values */
66         moldyn->t=MOLDYN_TEMP;
67         moldyn->tau=MOLDYN_TAU;
68         moldyn->time_steps=MOLDYN_RUNS;
69         moldyn->integrate=velocity_verlet;
70         moldyn->potential_force_function=lennard_jones;
71
72         /* parse argv */
73         for(i=1;i<argc;i++) {
74                 if(argv[i][0]=='-') {
75                         switch(argv[i][1]){
76                                 case 'E':
77                                         moldyn->ewrite=atoi(argv[++i]);
78                                         strncpy(moldyn->efb,argv[++i],64);
79                                         break;
80                                 case 'M':
81                                         moldyn->mwrite=atoi(argv[++i]);
82                                         strncpy(moldyn->mfb,argv[++i],64);
83                                         break;
84                                 case 'D':
85                                         moldyn->dwrite=atoi(argv[++i]);
86                                         strncpy(moldyn->dfb,argv[++i],64);
87                                         break;
88                                 case 'S':
89                                         moldyn->swrite=atoi(argv[++i]);
90                                         strncpy(moldyn->sfb,argv[++i],64);
91                                         break;
92                                 case 'V':
93                                         moldyn->vwrite=atoi(argv[++i]);
94                                         strncpy(moldyn->vfb,argv[++i],64);
95                                         break;
96                                 case 'T':
97                                         moldyn->t=atof(argv[++i]);
98                                         break;
99                                 case 't':
100                                         moldyn->tau=atof(argv[++i]);
101                                         break;
102                                 case 'C':
103                                         moldyn->cutoff=atof(argv[++i]);
104                                         break;
105                                 case 'R':
106                                         moldyn->time_steps=atoi(argv[++i]);
107                                         break;
108                                 case 'I':
109         /* integration algorithm */
110         switch(atoi(argv[++i])) {
111                 case MOLDYN_INTEGRATE_VERLET:
112                         moldyn->integrate=velocity_verlet;
113                         break;
114                 default:
115                         printf("unknown integration algo %s\n",argv[i]);
116                         moldyn_usage(argv);
117                         return -1;
118         }
119
120                                 case 'P':
121         /* potential + params */
122         switch(atoi(argv[++i])) {
123                 case MOLDYN_POTENTIAL_HO:
124                         hop.spring_constant=atof(argv[++i]);
125                         hop.equilibrium_distance=atof(argv[++i]);
126                         moldyn->pot_params=malloc(sizeof(t_ho_params));
127                         memcpy(moldyn->pot_params,&hop,sizeof(t_ho_params));
128                         moldyn->potential_force_function=harmonic_oscillator;
129                         break;
130                 case MOLDYN_POTENTIAL_LJ:
131                         e=atof(argv[++i]);
132                         s=atof(argv[++i]);
133                         ljp.epsilon4=4*e;
134                         ljp.sigma6=s*s*s*s*s*s;
135                         ljp.sigma12=ljp.sigma6*ljp.sigma6;
136                         moldyn->pot_params=malloc(sizeof(t_lj_params));
137                         memcpy(moldyn->pot_params,&ljp,sizeof(t_lj_params));
138                         moldyn->potential_force_function=lennard_jones;
139                         break;
140                 default:
141                         printf("unknown potential %s\n",argv[i]);
142                         moldyn_usage(argv);
143                         return -1;
144         }
145
146                                 default:
147                                         printf("unknown option %s\n",argv[i]);
148                                         moldyn_usage(argv);
149                                         return -1;
150                         }
151                 } else {
152                         moldyn_usage(argv);
153                         return -1;
154                 }
155         }
156
157         return 0;
158 }
159
160 int moldyn_log_init(t_moldyn *moldyn) {
161
162         moldyn->lvstat=0;
163         t_visual *vis;
164
165         vis=&(moldyn->vis);
166
167         if(moldyn->ewrite) {
168                 moldyn->efd=open(moldyn->efb,O_WRONLY|O_CREAT|O_TRUNC);
169                 if(moldyn->efd<0) {
170                         perror("[moldyn] efd open");
171                         return moldyn->efd;
172                 }
173                 dprintf(moldyn->efd,"# moldyn total energy logfile\n");
174                 moldyn->lvstat|=MOLDYN_LVSTAT_TOTAL_E;
175         }
176
177         if(moldyn->mwrite) {
178                 moldyn->mfd=open(moldyn->mfb,O_WRONLY|O_CREAT|O_TRUNC);
179                 if(moldyn->mfd<0) {
180                         perror("[moldyn] mfd open");
181                         return moldyn->mfd;
182                 }
183                 dprintf(moldyn->mfd,"# moldyn total momentum logfile\n");
184                 moldyn->lvstat|=MOLDYN_LVSTAT_TOTAL_M;
185         }
186
187         if(moldyn->swrite)
188                 moldyn->lvstat|=MOLDYN_LVSTAT_SAVE;
189
190         if(moldyn->dwrite) {
191                 moldyn->dfd=open(moldyn->dfb,O_WRONLY|O_CREAT|O_TRUNC);
192                 if(moldyn->dfd<0) {
193                         perror("[moldyn] dfd open");
194                         return moldyn->dfd;
195                 }
196                 write(moldyn->dfd,moldyn,sizeof(t_moldyn));
197                 moldyn->lvstat|=MOLDYN_LVSTAT_DUMP;
198         }
199
200         if((moldyn->vwrite)&&(vis)) {
201                 moldyn->visual=vis;
202                 visual_init(vis,moldyn->vfb);
203                 moldyn->lvstat|=MOLDYN_LVSTAT_VISUAL;
204         }
205
206         moldyn->lvstat|=MOLDYN_LVSTAT_INITIALIZED;
207
208         return 0;
209 }
210
211 int moldyn_log_shutdown(t_moldyn *moldyn) {
212
213         if(moldyn->efd) close(moldyn->efd);
214         if(moldyn->mfd) close(moldyn->efd);
215         if(moldyn->dfd) close(moldyn->efd);
216         if(moldyn->visual) visual_tini(moldyn->visual);
217
218         return 0;
219 }
220
221 int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
222
223         int ret;
224
225         ret=moldyn_parse_argv(moldyn,argc,argv);
226         if(ret<0) return ret;
227
228         ret=moldyn_log_init(moldyn);
229         if(ret<0) return ret;
230
231         rand_init(&(moldyn->random),NULL,1);
232         moldyn->random.status|=RAND_STAT_VERBOSE;
233
234         moldyn->status=0;
235
236         return 0;
237 }
238
239 int moldyn_shutdown(t_moldyn *moldyn) {
240
241         moldyn_log_shutdown(moldyn);
242         rand_close(&(moldyn->random));
243         free(moldyn->atom);
244
245         return 0;
246 }
247
248 int create_lattice(unsigned char type,int element,double mass,double lc,
249                    int a,int b,int c,t_atom **atom) {
250
251         int count;
252         int ret;
253         t_3dvec origin;
254
255         count=a*b*c;
256
257         if(type==FCC) count*=4;
258         if(type==DIAMOND) count*=8;
259
260         *atom=malloc(count*sizeof(t_atom));
261         if(*atom==NULL) {
262                 perror("malloc (atoms)");
263                 return -1;
264         }
265
266         v3_zero(&origin);
267
268         switch(type) {
269                 case FCC:
270                         ret=fcc_init(a,b,c,lc,*atom,&origin);
271                         break;
272                 case DIAMOND:
273                         ret=diamond_init(a,b,c,lc,*atom,&origin);
274                         break;
275                 default:
276                         printf("unknown lattice type (%02x)\n",type);
277                         return -1;
278         }
279
280         /* debug */
281         if(ret!=count) {
282                 printf("ok, there is something wrong ...\n");
283                 printf("calculated -> %d atoms\n",count);
284                 printf("created -> %d atoms\n",ret);
285                 return -1;
286         }
287
288         while(count) {
289                 (*atom)[count-1].element=element;
290                 (*atom)[count-1].mass=mass;
291                 count-=1;
292         }
293
294         return ret;
295 }
296
297 int destroy_lattice(t_atom *atom) {
298
299         if(atom) free(atom);
300
301         return 0;
302 }
303
304 int thermal_init(t_moldyn *moldyn) {
305
306         /*
307          * - gaussian distribution of velocities
308          * - zero total momentum
309          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
310          */
311
312         int i;
313         double v,sigma;
314         t_3dvec p_total,delta;
315         t_atom *atom;
316         t_random *random;
317
318         atom=moldyn->atom;
319         random=&(moldyn->random);
320
321         /* gaussian distribution of velocities */
322         v3_zero(&p_total);
323         for(i=0;i<moldyn->count;i++) {
324                 sigma=sqrt(2.0*K_BOLTZMANN*moldyn->t/atom[i].mass);
325                 /* x direction */
326                 v=sigma*rand_get_gauss(random);
327                 atom[i].v.x=v;
328                 p_total.x+=atom[i].mass*v;
329                 /* y direction */
330                 v=sigma*rand_get_gauss(random);
331                 atom[i].v.y=v;
332                 p_total.y+=atom[i].mass*v;
333                 /* z direction */
334                 v=sigma*rand_get_gauss(random);
335                 atom[i].v.z=v;
336                 p_total.z+=atom[i].mass*v;
337         }
338
339         /* zero total momentum */
340         v3_scale(&p_total,&p_total,1.0/moldyn->count);
341         for(i=0;i<moldyn->count;i++) {
342                 v3_scale(&delta,&p_total,1.0/atom[i].mass);
343                 v3_sub(&(atom[i].v),&(atom[i].v),&delta);
344         }
345
346         /* velocity scaling */
347         scale_velocity(moldyn);
348
349         return 0;
350 }
351
352 int scale_velocity(t_moldyn *moldyn) {
353
354         int i;
355         double e,c;
356         t_atom *atom;
357
358         atom=moldyn->atom;
359
360         /*
361          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
362          */
363         e=0.0;
364         for(i=0;i<moldyn->count;i++)
365                 e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
366         c=sqrt((2.0*e)/(3.0*moldyn->count*K_BOLTZMANN*moldyn->t));
367         for(i=0;i<moldyn->count;i++)
368                 v3_scale(&(atom[i].v),&(atom[i].v),(1.0/c));
369
370         return 0;
371 }
372
373 double get_e_kin(t_atom *atom,int count) {
374
375         int i;
376         double e;
377
378         e=0.0;
379
380         for(i=0;i<count;i++) {
381                 e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
382         }
383
384         return e;
385 }
386
387 double get_e_pot(t_moldyn *moldyn) {
388
389         return moldyn->energy;
390 }
391
392 double get_total_energy(t_moldyn *moldyn) {
393
394         double e;
395
396         e=get_e_kin(moldyn->atom,moldyn->count);
397         e+=get_e_pot(moldyn);
398
399         return e;
400 }
401
402 t_3dvec get_total_p(t_atom *atom, int count) {
403
404         t_3dvec p,p_total;
405         int i;
406
407         v3_zero(&p_total);
408         for(i=0;i<count;i++) {
409                 v3_scale(&p,&(atom[i].v),atom[i].mass);
410                 v3_add(&p_total,&p_total,&p);
411         }
412
413         return p_total;
414 }
415
416 double estimate_time_step(t_moldyn *moldyn,double nn_dist,double t) {
417
418         double tau;
419
420         tau=0.05*nn_dist/(sqrt(3.0*K_BOLTZMANN*t/moldyn->atom[0].mass));
421         tau*=1.0E-9;
422         if(tau<moldyn->tau)
423                 printf("[moldyn] warning: time step  (%f > %.15f)\n",
424                        moldyn->tau,tau);
425
426         return tau;     
427 }
428
429 /*
430  * numerical tricks
431  */
432
433 /* linked list / cell method */
434
435 int link_cell_init(t_moldyn *moldyn) {
436
437         t_linkcell *lc;
438         int i;
439
440         lc=&(moldyn->lc);
441
442         /* list log fd */
443         lc->listfd=open("/dev/null",O_WRONLY);
444
445         /* partitioning the md cell */
446         lc->nx=moldyn->dim.x/moldyn->cutoff;
447         lc->x=moldyn->dim.x/lc->nx;
448         lc->ny=moldyn->dim.y/moldyn->cutoff;
449         lc->y=moldyn->dim.y/lc->ny;
450         lc->nz=moldyn->dim.z/moldyn->cutoff;
451         lc->z=moldyn->dim.z/lc->nz;
452
453         lc->cells=lc->nx*lc->ny*lc->nz;
454         lc->subcell=malloc(lc->cells*sizeof(t_list));
455
456         printf("initializing linked cells (%d)\n",lc->cells);
457
458         for(i=0;i<lc->cells;i++)
459                 //list_init(&(lc->subcell[i]),1);
460                 list_init(&(lc->subcell[i]));
461
462         link_cell_update(moldyn);
463         
464         return 0;
465 }
466
467 int link_cell_update(t_moldyn *moldyn) {
468
469         int count,i,j,k;
470         int nx,ny,nz;
471         t_atom *atom;
472         t_linkcell *lc;
473
474         atom=moldyn->atom;
475         lc=&(moldyn->lc);
476
477         nx=lc->nx;
478         ny=lc->ny;
479         nz=lc->nz;
480
481         for(i=0;i<lc->cells;i++)
482                 list_destroy(&(moldyn->lc.subcell[i]));
483         
484         for(count=0;count<moldyn->count;count++) {
485                 i=(atom[count].r.x+(moldyn->dim.x/2))/lc->x;
486                 j=(atom[count].r.y+(moldyn->dim.y/2))/lc->y;
487                 k=(atom[count].r.z+(moldyn->dim.z/2))/lc->z;
488                 list_add_immediate_ptr(&(moldyn->lc.subcell[i+j*nx+k*nx*ny]),
489                                        &(atom[count]));
490         }
491
492         return 0;
493 }
494
495 int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) {
496
497         t_linkcell *lc;
498         int a;
499         int count1,count2;
500         int ci,cj,ck;
501         int nx,ny,nz;
502         int x,y,z;
503         unsigned char bx,by,bz;
504
505         lc=&(moldyn->lc);
506         nx=lc->nx;
507         ny=lc->ny;
508         nz=lc->nz;
509         count1=1;
510         count2=27;
511         a=nx*ny;
512
513
514         cell[0]=lc->subcell[i+j*nx+k*a];
515         for(ci=-1;ci<=1;ci++) {
516                 bx=0;
517                 x=i+ci;
518                 if((x<0)||(x>=nx)) {
519                         x=(x+nx)%nx;
520                         bx=1;
521                 }
522                 for(cj=-1;cj<=1;cj++) {
523                         by=0;
524                         y=j+cj;
525                         if((y<0)||(y>=ny)) {
526                                 y=(y+ny)%ny;
527                                 by=1;
528                         }
529                         for(ck=-1;ck<=1;ck++) {
530                                 bz=0;
531                                 z=k+ck;
532                                 if((z<0)||(z>=nz)) {
533                                         z=(z+nz)%nz;
534                                         bz=1;
535                                 }
536                                 if(!(ci|cj|ck)) continue;
537                                 if(bx|by|bz) {
538                                         cell[--count2]=lc->subcell[x+y*nx+z*a];
539                                 }
540                                 else {
541                                         cell[count1++]=lc->subcell[x+y*nx+z*a];
542                                 }
543                         }
544                 }
545         }
546
547         return count2;
548 }
549
550 int link_cell_shutdown(t_moldyn *moldyn) {
551
552         int i;
553         t_linkcell *lc;
554
555         lc=&(moldyn->lc);
556
557         for(i=0;i<lc->nx*lc->ny*lc->nz;i++)
558                 list_shutdown(&(moldyn->lc.subcell[i]));
559
560         if(lc->listfd) close(lc->listfd);
561
562         return 0;
563 }
564
565 /*
566  *
567  * 'integration of newtons equation' - algorithms
568  *
569  */
570
571 /* start the integration */
572
573 int moldyn_integrate(t_moldyn *moldyn) {
574
575         int i;
576         unsigned int e,m,s,d,v;
577         t_3dvec p;
578
579         int fd;
580         char fb[128];
581
582         /* initialize linked cell method */
583         link_cell_init(moldyn);
584
585         /* logging & visualization */
586         e=moldyn->ewrite;
587         m=moldyn->mwrite;
588         s=moldyn->swrite;
589         d=moldyn->dwrite;
590         v=moldyn->vwrite;
591
592         if(!(moldyn->lvstat&MOLDYN_LVSTAT_INITIALIZED)) {
593                 printf("[moldyn] warning, lv system not initialized\n");
594                 return -1;
595         }
596
597         /* sqaure of some variables */
598         moldyn->tau_square=moldyn->tau*moldyn->tau;
599         moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff;
600
601         /* calculate initial forces */
602         moldyn->potential_force_function(moldyn);
603
604         for(i=0;i<moldyn->time_steps;i++) {
605
606                 /* integration step */
607                 moldyn->integrate(moldyn);
608
609                 /* check for log & visualization */
610                 if(e) {
611                         if(!(i%e))
612                                 dprintf(moldyn->efd,
613                                         "%.15f %.45f\n",i*moldyn->tau,
614                                         get_total_energy(moldyn));
615                 }
616                 if(m) {
617                         if(!(i%m)) {
618                                 p=get_total_p(moldyn->atom,moldyn->count);
619                                 dprintf(moldyn->mfd,
620                                         "%.15f %.45f\n",i*moldyn->tau,
621                                         v3_norm(&p));
622                         }
623                 }
624                 if(s) {
625                         if(!(i%s)) {
626                                 snprintf(fb,128,"%s-%f-%.15f.save",moldyn->sfb,
627                                          moldyn->t,i*moldyn->tau);
628                                 fd=open(fb,O_WRONLY|O_TRUNC|O_CREAT);
629                                 if(fd<0) perror("[moldyn] save fd open");
630                                 else {
631                                         write(fd,moldyn,sizeof(t_moldyn));
632                                         write(fd,moldyn->atom,
633                                               moldyn->count*sizeof(t_atom));
634                                 }
635                         }       
636                 }
637                 if(d) {
638                         if(!(i%d))
639                                 write(moldyn->dfd,moldyn->atom,
640                                       moldyn->count*sizeof(t_atom));
641
642                 }
643                 if(v) {
644                         if(!(i%v)) {
645                                 visual_atoms(moldyn->visual,i*moldyn->tau,
646                                              moldyn->atom,moldyn->count);
647                                 printf("\rsteps: %d",i);
648                                 fflush(stdout);
649                         }
650                 }
651         }
652
653         return 0;
654 }
655
656 /* velocity verlet */
657
658 int velocity_verlet(t_moldyn *moldyn) {
659
660         int i,count;
661         double tau,tau_square;
662         t_3dvec delta;
663         t_atom *atom;
664
665         atom=moldyn->atom;
666         count=moldyn->count;
667         tau=moldyn->tau;
668         tau_square=moldyn->tau_square;
669
670         for(i=0;i<count;i++) {
671                 /* new positions */
672                 v3_scale(&delta,&(atom[i].v),tau);
673                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
674                 v3_scale(&delta,&(atom[i].f),0.5*tau_square/atom[i].mass);
675                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
676                 v3_per_bound(&(atom[i].r),&(moldyn->dim));
677
678                 /* velocities */
679                 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
680                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
681         }
682
683         /* neighbour list update */
684 printf("list update ...\n");
685         link_cell_update(moldyn);
686 printf("done\n");
687
688         /* forces depending on chosen potential */
689 printf("calc potential/force ...\n");
690         moldyn->potential_force_function(moldyn);
691 printf("done\n");
692
693         for(i=0;i<count;i++) {
694                 /* again velocities */
695                 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
696                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
697         }
698
699         return 0;
700 }
701
702
703 /*
704  *
705  * potentials & corresponding forces
706  * 
707  */
708
709 /* harmonic oscillator potential and force */
710
711 int harmonic_oscillator(t_moldyn *moldyn) {
712
713         t_ho_params *params;
714         t_atom *atom,*btom;
715         t_linkcell *lc;
716         t_list *this,neighbour[27];
717         int i,j,c;
718         int count;
719         t_3dvec force,distance;
720         double d,u;
721         double sc,equi_dist;
722         int ni,nj,nk;
723
724         params=moldyn->pot_params;
725         atom=moldyn->atom;
726         lc=&(moldyn->lc);
727         sc=params->spring_constant;
728         equi_dist=params->equilibrium_distance;
729         count=moldyn->count;
730
731         /* reset energy counter */
732         u=0.0;
733
734         for(i=0;i<count;i++) {
735                 /* reset force */
736                 v3_zero(&(atom[i].f));
737
738                 /* determine cell + neighbours */
739                 ni=(atom[i].r.x+(moldyn->dim.x/2))/lc->x;
740                 nj=(atom[i].r.y+(moldyn->dim.y/2))/lc->y;
741                 nk=(atom[i].r.z+(moldyn->dim.z/2))/lc->z;
742                 c=link_cell_neighbour_index(moldyn,ni,nj,nk,neighbour);
743
744                 /*
745                  * processing cell of atom i
746                  * => no need to check for empty list (1 element at minimum)
747                  */
748                 this=&(neighbour[0]);
749                 list_reset(this);
750                 do {
751                         btom=this->current->data;
752                         if(btom==&(atom[i]))
753                                 continue;
754                         v3_sub(&distance,&(atom[i].r),&(btom->r));
755                         d=v3_norm(&distance);
756                         if(d<=moldyn->cutoff) {
757                                 u+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
758                                 v3_scale(&force,&distance,
759                                          -sc*(1.0-(equi_dist/d)));
760                                 v3_add(&(atom[i].f),&(atom[i].f),&force);
761                         }
762                 } while(list_next(this)!=L_NO_NEXT_ELEMENT);
763
764                 /*
765                  * direct neighbour cells
766                  * => no boundary condition check necessary
767                  */
768                 for(j=1;j<c;j++) {
769                         this=&(neighbour[j]);
770                         list_reset(this); /* there might not be a single atom */
771                         if(this->start!=NULL) {
772
773                         do {
774                                 btom=this->current->data;
775                                 v3_sub(&distance,&(atom[i].r),&(btom->r));
776                                 d=v3_norm(&distance);
777                                 if(d<=moldyn->cutoff) {
778                                         u+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
779                                         v3_scale(&force,&distance,
780                                                  -sc*(1.0-(equi_dist/d)));
781                                         v3_add(&(atom[i].f),&(atom[i].f),
782                                                &force);
783                                 }
784                         } while(list_next(this)!=L_NO_NEXT_ELEMENT);
785
786                         }
787                 }
788
789                 /*
790                  * indirect neighbour cells
791                  * => check boundary conditions
792                  */
793                 for(j=c;j<27;j++) {
794                         this=&(neighbour[j]);
795                         list_reset(this); /* check boundary conditions */
796                         if(this->start!=NULL) {
797
798                         do {
799                                 btom=this->current->data;
800                                 v3_sub(&distance,&(atom[i].r),&(btom->r));
801                                 v3_per_bound(&distance,&(moldyn->dim));
802                                 d=v3_norm(&distance);
803                                 if(d<=moldyn->cutoff) {
804                                         u+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
805                                         v3_scale(&force,&distance,
806                                                  -sc*(1.0-(equi_dist/d)));
807                                         v3_add(&(atom[i].f),&(atom[i].f),
808                                                &force);
809                                 }
810                         } while(list_next(this)!=L_NO_NEXT_ELEMENT);
811
812                         }
813                 }
814         }
815
816         moldyn->energy=0.5*u;
817
818         return 0;
819 }
820
821 /* lennard jones potential & force for one sort of atoms */
822  
823 int lennard_jones(t_moldyn *moldyn) {
824
825         t_lj_params *params;
826         t_atom *atom,*btom;
827         t_linkcell *lc;
828         t_list *this,neighbour[27];
829         int i,j,c;
830         int count;
831         t_3dvec force,distance;
832         double d,h1,h2,u;
833         double eps,sig6,sig12;
834         int ni,nj,nk;
835
836         params=moldyn->pot_params;
837         atom=moldyn->atom;
838         lc=&(moldyn->lc);
839         count=moldyn->count;
840         eps=params->epsilon4;
841         sig6=params->sigma6;
842         sig12=params->sigma12;
843
844         /* reset energy counter */
845         u=0.0;
846
847         for(i=0;i<count;i++) {
848                 /* reset force */
849                 v3_zero(&(atom[i].f));
850
851                 /* determine cell + neighbours */
852                 ni=(atom[i].r.x+(moldyn->dim.x/2))/lc->x;
853                 nj=(atom[i].r.y+(moldyn->dim.y/2))/lc->y;
854                 nk=(atom[i].r.z+(moldyn->dim.z/2))/lc->z;
855                 c=link_cell_neighbour_index(moldyn,ni,nj,nk,neighbour);
856
857                 /* processing cell of atom i */
858                 this=&(neighbour[0]);
859                 list_reset(this); /* list has 1 element at minimum */
860                 do {
861                         btom=this->current->data;
862                         if(btom==&(atom[i]))
863                                 continue;
864                         v3_sub(&distance,&(atom[i].r),&(btom->r));
865                         d=v3_absolute_square(&distance);        /* 1/r^2 */
866                         if(d<=moldyn->cutoff_square) {
867                                 d=1.0/d;        /* 1/r^2 */
868                                 h2=d*d;         /* 1/r^4 */
869                                 h2*=d;          /* 1/r^6 */
870                                 h1=h2*h2;       /* 1/r^12 */
871                                 u+=eps*(sig12*h1-sig6*h2);
872                                 h2*=d;          /* 1/r^8 */
873                                 h1*=d;          /* 1/r^14 */
874                                 h2*=6*sig6;
875                                 h1*=12*sig12;
876                                 d=+h1-h2;
877                                 d*=eps;
878                                 v3_scale(&force,&distance,d);
879                                 v3_add(&(atom[i].f),&(atom[i].f),&force);
880                         }
881                 } while(list_next(this)!=L_NO_NEXT_ELEMENT);
882
883                 /* neighbours not doing boundary condition overflow */
884                 for(j=1;j<c;j++) {
885                         this=&(neighbour[j]);
886                         list_reset(this); /* there might not be a single atom */
887                         if(this->start!=NULL) {
888
889                         do {
890                                 btom=this->current->data;
891                                 v3_sub(&distance,&(atom[i].r),&(btom->r));
892                                 d=v3_absolute_square(&distance); /* r^2 */
893                                 if(d<=moldyn->cutoff_square) {
894                                         d=1.0/d;        /* 1/r^2 */
895                                         h2=d*d;         /* 1/r^4 */
896                                         h2*=d;          /* 1/r^6 */
897                                         h1=h2*h2;       /* 1/r^12 */
898                                         u+=eps*(sig12*h1-sig6*h2);
899                                         h2*=d;          /* 1/r^8 */
900                                         h1*=d;          /* 1/r^14 */
901                                         h2*=6*sig6;
902                                         h1*=12*sig12;
903                                         d=+h1-h2;
904                                         d*=eps;
905                                         v3_scale(&force,&distance,d);
906                                         v3_add(&(atom[i].f),&(atom[i].f),
907                                                &force);
908                                 }
909                         } while(list_next(this)!=L_NO_NEXT_ELEMENT);
910                                 
911                         }
912                 }
913
914                 /* neighbours due to boundary conditions */
915                 for(j=c;j<27;j++) {
916                         this=&(neighbour[j]);
917                         list_reset(this); /* check boundary conditions */
918                         if(this->start!=NULL) {
919
920                         do {
921                                 btom=this->current->data;
922                                 v3_sub(&distance,&(atom[i].r),&(btom->r));
923                                 v3_per_bound(&distance,&(moldyn->dim));
924                                 d=v3_absolute_square(&distance); /* r^2 */
925                                 if(d<=moldyn->cutoff_square) {
926                                         d=1.0/d;        /* 1/r^2 */
927                                         h2=d*d;         /* 1/r^4 */
928                                         h2*=d;          /* 1/r^6 */
929                                         h1=h2*h2;       /* 1/r^12 */
930                                         u+=eps*(sig12*h1-sig6*h2);
931                                         h2*=d;          /* 1/r^8 */
932                                         h1*=d;          /* 1/r^14 */
933                                         h2*=6*sig6;
934                                         h1*=12*sig12;
935                                         d=+h1-h2;
936                                         d*=eps;
937                                         v3_scale(&force,&distance,d);
938                                         v3_add(&(atom[i].f),&(atom[i].f),
939                                                &force);
940                                 }
941                         } while(list_next(this)!=L_NO_NEXT_ELEMENT);
942
943                         }
944                 }
945         }
946
947         moldyn->energy=0.5*u;
948         
949         return 0;
950 }
951
952 /* tersoff potential & force for 2 sorts of atoms */
953
954 int tersoff(t_moldyn *moldyn) {
955
956         t_tersoff_params *params;
957         t_atom *atom,*btom,*ktom;
958         t_linkcell *lc;
959         t_list *this,*thisk,neighbour[27],neighbourk[27];
960         int i,j,k,c,ck;
961         int count;
962         double u;
963         int ni,nj,nk;
964         int ki,kj,kk;
965         
966
967         params=moldyn->pot_params;
968         atom=moldyn->atom;
969         lc=&(moldyn->lc);
970         count=moldyn->count;
971         
972         /* reset energy counter */
973         u=0.0;
974
975         for(i=0;i<count;i++) {
976                 /* reset force */
977                 v3_zero(&(atom[i].f));
978
979                 /* determin cell neighbours */
980                 ni=(atom[i].r.x+(moldyn->dim.x/2))/lc->x;
981                 nj=(atom[i].r.y+(moldyn->dim.y/2))/lc->y;
982                 nk=(atom[i].r.z+(moldyn->dim.z/2))/lc->z;
983                 c=link_cell_neighbour_index(moldyn,ni,nj,nk,neighbour);
984
985                 /*
986                  * processing cell of atom i
987                  * => no need to check for empty list (1 element at minimum)
988                  */
989                 this=&(neighbour[0]);
990                 list_reset(this);
991                 do {
992                         btom=this->current->data;
993                         if(btom==&(atom[i]))
994                                 continue;
995
996                         /* 2 body stuff */
997
998                         v3_sub(&dist_ij,btom,&(atom[i]));
999                         d_ij=v3_norm(&dist_ij);
1000                         if(d_ij<=S) {
1001
1002                                 S=;
1003                                 R=;
1004                                 A=;
1005                                 lambda=;
1006                                 B=;
1007                                 mu=;
1008
1009                                 if(d_ij<=R) {
1010                                         f_c=1.0;
1011                                         df_c=0.0;
1012                                 }
1013                                 else {
1014                                         s_r=S-R;
1015                                         arg1=PI*(d_ij-R)/s_r;
1016                                         f_c=0.5+0.5*cos(arg1);
1017                                         df_c=-0.5*sin(arg1)*(PI/(s_r*d_ij));
1018                                         f_r=A*exp(-lambda*d_ij);
1019                                         f_a=-B*exp(-mu*d_ij);
1020                                 }
1021                         }
1022                         else 
1023                                 continue;               
1024
1025                         
1026                         /* end 2 body stuff */ 
1027
1028                         /* determine cell neighbours of btom */
1029                         ki=(btom->r.x+(moldyn->dim.x/2))/lc->x;
1030                         kj=(btom->r.y+(moldyn->dim.y/2))/lc->y;
1031                         kk=(btom->r.z+(moldyn->dim.z/2))/lc->z;
1032                         ck=link_cell_neighbour_index(moldyn,ki,kj,kk,
1033                                                      neighbourk);
1034
1035                         /* cell of btom */
1036                         thisk=&(neighbourk[0]);
1037                         list_reset(thisk);
1038                         do {
1039                                 ktom=thisk->current->data;
1040                                 if(ktom==btom)
1041                                         continue;
1042                                 if(ktom==&(atom[i]))
1043                                         continue;
1044                                 
1045                                 /* 3 body stuff (1) */
1046                                 
1047                                 theta_ijk=;
1048                                 sin_theta=;
1049                                 cos_theta=;
1050                                 hi_cos=;
1051                                 hi_cos_square=;
1052
1053
1054                                 /* end 3 body stuff (1) */
1055
1056
1057                         } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
1058
1059                         /* direct neighbours of btom cell */
1060                         for(k=1;k<ck;k++) {
1061                                 thisk=&(neighbourk[k]);
1062                                 list_reset(thisk);
1063                                 if(thisk->start!=NULL) {
1064
1065                                 do {
1066                                         ktom=thisk->current->data;
1067                                         if(ktom==&(atom[i]))
1068                                                 continue;
1069
1070                                 /* 3 body stuff (2) */
1071
1072                                 } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
1073
1074                                 }
1075                         }
1076
1077                         /* indirect neighbours of btom cell */
1078                         for(k=ck;k<27;k++) {
1079                                 thisk=&(neighbourk[k]);
1080                                 list_reset(thisk);
1081                                 if(thisk->start!=NULL) {
1082
1083                                 do {
1084                                         ktom=thisk->current->data;
1085                                         if(ktom==&(atom[i]))
1086                                                 continue;
1087
1088                                 /* 3 body stuff */
1089
1090                                 } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
1091
1092                                 }
1093                         }
1094
1095
1096                 } while(list_next(this)!=L_NO_NEXT_ELEMENT);
1097
1098                 /*
1099                  * direct neighbour cells of atom i
1100                  */
1101                 for(j=1;j<c;j++) {
1102                         this=&(neighbour[j]);
1103                         list_reset(this);
1104                         if(this->start!=NULL) {
1105
1106                         do {
1107                                 btom=this->current->data;
1108
1109                                 /* 2 body stuff */
1110
1111
1112                         /* determine cell neighbours of btom */
1113                         ki=(btom->r.x+(moldyn->dim.x/2))/lc->x;
1114                         kj=(btom->r.y+(moldyn->dim.y/2))/lc->y;
1115                         kk=(btom->r.z+(moldyn->dim.z/2))/lc->z;
1116                         ck=link_cell_neighbour_index(moldyn,ki,kj,kk,
1117                                                      neighbourk);
1118
1119                         /* cell of btom */
1120                         thisk=&(neighbourk[0]);
1121                         list_reset(thisk);
1122                         do {
1123                                 ktom=thisk->current->data;
1124                                 if(ktom==btom)
1125                                         continue;
1126                                 if(ktom==&(atom[i]))
1127                                         continue;
1128                                 
1129                                 /* 3 body stuff (1) */
1130
1131                         } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
1132
1133                         /* direct neighbours of btom cell */
1134                         for(k=1;k<ck;k++) {
1135                                 thisk=&(neighbourk[k]);
1136                                 list_reset(thisk);
1137                                 if(thisk->start!=NULL) {
1138
1139                                 do {
1140                                         ktom=thisk->current->data;
1141                                         if(ktom==&(atom[i]))
1142                                                 continue;
1143
1144                                 /* 3 body stuff (2) */
1145
1146                                 } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
1147
1148                                 }
1149                         }
1150
1151                         /* indirect neighbours of btom cell */
1152                         for(k=ck;k<27;k++) {
1153                                 thisk=&(neighbourk[k]);
1154                                 list_reset(thisk);
1155                                 if(thisk->start!=NULL) {
1156
1157                                 do {
1158                                         ktom=thisk->current->data;
1159                                         if(ktom==&(atom[i]))
1160                                                 continue;
1161
1162                                 /* 3 body stuff (3) */
1163
1164                                 } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
1165
1166                                 }
1167                         }
1168
1169
1170                         } while(list_next(this)!=L_NO_NEXT_ELEMENT);
1171
1172                         }
1173                 }
1174
1175                 /*
1176                  * indirect neighbour cells of atom i
1177                  */
1178                 for(j=c;j<27;j++) {
1179                         this=&(neighbour[j]);
1180                         list_reset(this);
1181                         if(this->start!=NULL) {
1182
1183                         do {
1184                                 btom=this->current->data;
1185
1186                                 /* 2 body stuff */
1187
1188
1189                         /* determine cell neighbours of btom */
1190                         ki=(btom->r.x+(moldyn->dim.x/2))/lc->x;
1191                         kj=(btom->r.y+(moldyn->dim.y/2))/lc->y;
1192                         kk=(btom->r.z+(moldyn->dim.z/2))/lc->z;
1193                         ck=link_cell_neighbour_index(moldyn,ki,kj,kk,
1194                                                      neighbourk);
1195
1196                         /* cell of btom */
1197                         thisk=&(neighbourk[0]);
1198                         list_reset(thisk);
1199                         do {
1200                                 ktom=thisk->current->data;
1201                                 if(ktom==btom)
1202                                         continue;
1203                                 if(ktom==&(atom[i]))
1204                                         continue;
1205                                 
1206                                 /* 3 body stuff (1) */
1207
1208                         } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
1209
1210                         /* direct neighbours of btom cell */
1211                         for(k=1;k<ck;k++) {
1212                                 thisk=&(neighbourk[k]);
1213                                 list_reset(thisk);
1214                                 if(thisk->start!=NULL) {
1215
1216                                 do {
1217                                         ktom=thisk->current->data;
1218                                         if(ktom==&(atom[i]))
1219                                                 continue;
1220
1221                                 /* 3 body stuff (2) */
1222
1223                                 } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
1224
1225                                 }
1226                         }
1227
1228                         /* indirect neighbours of btom cell */
1229                         for(k=ck;k<27;k++) {
1230                                 thisk=&(neighbourk[k]);
1231                                 list_reset(thisk);
1232                                 if(thisk->start!=NULL) {
1233
1234                                 do {
1235                                         ktom=thisk->current->data;
1236                                         if(ktom==&(atom[i]))
1237                                                 continue;
1238
1239                                 /* 3 body stuff (3) */
1240
1241                                 } while(list_next(thisk)!=L_NO_NEXT_ELEMENT);
1242
1243                                 }
1244                         }
1245
1246
1247                         } while(list_next(this)!=L_NO_NEXT_ELEMENT);
1248
1249                         }
1250                 }
1251                 
1252         }
1253
1254         moldyn->energy=0.5*u;
1255
1256         return 0;
1257 }
1258