]> www.hackdaworld.org Git - physik/posic.git/blob - moldyn.c
's/for/while/' <- runtime schedule adding possible now
[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 #include "report/report.h"
20
21 int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
22
23         printf("[moldyn] init\n");
24
25         memset(moldyn,0,sizeof(t_moldyn));
26
27         rand_init(&(moldyn->random),NULL,1);
28         moldyn->random.status|=RAND_STAT_VERBOSE;
29
30         return 0;
31 }
32
33 int moldyn_shutdown(t_moldyn *moldyn) {
34
35         printf("[moldyn] shutdown\n");
36
37         moldyn_log_shutdown(moldyn);
38         link_cell_shutdown(moldyn);
39         rand_close(&(moldyn->random));
40         free(moldyn->atom);
41
42         return 0;
43 }
44
45 int set_int_alg(t_moldyn *moldyn,u8 algo) {
46
47         printf("[moldyn] integration algorithm: ");
48
49         switch(algo) {
50                 case MOLDYN_INTEGRATE_VERLET:
51                         moldyn->integrate=velocity_verlet;
52                         printf("velocity verlet\n");
53                         break;
54                 default:
55                         printf("unknown integration algorithm: %02x\n",algo);
56                         printf("unknown\n");
57                         return -1;
58         }
59
60         return 0;
61 }
62
63 int set_cutoff(t_moldyn *moldyn,double cutoff) {
64
65         moldyn->cutoff=cutoff;
66
67         printf("[moldyn] cutoff [A]: %f\n",moldyn->cutoff);
68
69         return 0;
70 }
71
72 int set_temperature(t_moldyn *moldyn,double t_ref) {
73
74         moldyn->t_ref=t_ref;
75
76         printf("[moldyn] temperature [K]: %f\n",moldyn->t_ref);
77
78         return 0;
79 }
80
81 int set_pressure(t_moldyn *moldyn,double p_ref) {
82
83         moldyn->p_ref=p_ref;
84
85         printf("[moldyn] pressure [bar]: %f\n",moldyn->p_ref/BAR);
86
87         return 0;
88 }
89
90 int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc) {
91
92         moldyn->pt_scale=(ptype|ttype);
93         moldyn->t_tc=ttc;
94         moldyn->p_tc=ptc;
95
96         printf("[moldyn] p/t scaling:\n");
97
98         printf("  p: %s",ptype?"yes":"no ");
99         if(ptype)
100                 printf(" | type: %02x | factor: %f",ptype,ptc);
101         printf("\n");
102
103         printf("  t: %s",ttype?"yes":"no ");
104         if(ttype)
105                 printf(" | type: %02x | factor: %f",ttype,ttc);
106         printf("\n");
107
108         return 0;
109 }
110
111 int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize) {
112
113         moldyn->dim.x=x;
114         moldyn->dim.y=y;
115         moldyn->dim.z=z;
116
117         moldyn->volume=x*y*z;
118
119         if(visualize) {
120                 moldyn->vis.dim.x=x;
121                 moldyn->vis.dim.y=y;
122                 moldyn->vis.dim.z=z;
123         }
124
125         moldyn->dv=0.000001*moldyn->volume;
126
127         printf("[moldyn] dimensions in A and A^3 respectively:\n");
128         printf("  x: %f\n",moldyn->dim.x);
129         printf("  y: %f\n",moldyn->dim.y);
130         printf("  z: %f\n",moldyn->dim.z);
131         printf("  volume: %f\n",moldyn->volume);
132         printf("  visualize simulation box: %s\n",visualize?"yes":"no");
133         printf("  delta volume (pressure calc): %f\n",moldyn->dv);
134
135         return 0;
136 }
137
138 int set_nn_dist(t_moldyn *moldyn,double dist) {
139
140         moldyn->nnd=dist;
141
142         return 0;
143 }
144
145 int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z) {
146
147         printf("[moldyn] periodic boundary conditions:\n");
148
149         if(x)
150                 moldyn->status|=MOLDYN_STAT_PBX;
151
152         if(y)
153                 moldyn->status|=MOLDYN_STAT_PBY;
154
155         if(z)
156                 moldyn->status|=MOLDYN_STAT_PBZ;
157
158         printf("  x: %s\n",x?"yes":"no");
159         printf("  y: %s\n",y?"yes":"no");
160         printf("  z: %s\n",z?"yes":"no");
161
162         return 0;
163 }
164
165 int set_potential1b(t_moldyn *moldyn,pf_func1b func) {
166
167         moldyn->func1b=func;
168
169         return 0;
170 }
171
172 int set_potential2b(t_moldyn *moldyn,pf_func2b func) {
173
174         moldyn->func2b=func;
175
176         return 0;
177 }
178
179 int set_potential3b_j1(t_moldyn *moldyn,pf_func2b func) {
180
181         moldyn->func3b_j1=func;
182
183         return 0;
184 }
185
186 int set_potential3b_j2(t_moldyn *moldyn,pf_func2b func) {
187
188         moldyn->func3b_j2=func;
189
190         return 0;
191 }
192
193 int set_potential3b_j3(t_moldyn *moldyn,pf_func2b func) {
194
195         moldyn->func3b_j3=func;
196
197         return 0;
198 }
199
200 int set_potential3b_k1(t_moldyn *moldyn,pf_func3b func) {
201
202         moldyn->func3b_k1=func;
203
204         return 0;
205 }
206
207 int set_potential3b_k2(t_moldyn *moldyn,pf_func3b func) {
208
209         moldyn->func3b_k2=func;
210
211         return 0;
212 }
213
214 int set_potential_params(t_moldyn *moldyn,void *params) {
215
216         moldyn->pot_params=params;
217
218         return 0;
219 }
220
221 int set_avg_skip(t_moldyn *moldyn,int skip) {
222
223         printf("[moldyn] skip %d steps before starting average calc\n",skip);
224         moldyn->avg_skip=skip;
225
226         return 0;
227 }
228
229 int moldyn_set_log_dir(t_moldyn *moldyn,char *dir) {
230
231         strncpy(moldyn->vlsdir,dir,127);
232
233         return 0;
234 }
235
236 int moldyn_set_report(t_moldyn *moldyn,char *author,char *title) {
237
238         strncpy(moldyn->rauthor,author,63);
239         strncpy(moldyn->rtitle,title,63);
240
241         return 0;
242 }
243         
244 int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) {
245
246         char filename[128];
247         int ret;
248
249         printf("[moldyn] set log: ");
250
251         switch(type) {
252                 case LOG_TOTAL_ENERGY:
253                         moldyn->ewrite=timer;
254                         snprintf(filename,127,"%s/energy",moldyn->vlsdir);
255                         moldyn->efd=open(filename,
256                                          O_WRONLY|O_CREAT|O_EXCL,
257                                          S_IRUSR|S_IWUSR);
258                         if(moldyn->efd<0) {
259                                 perror("[moldyn] energy log fd open");
260                                 return moldyn->efd;
261                         }
262                         dprintf(moldyn->efd,"# total energy log file\n");
263                         printf("total energy (%d)\n",timer);
264                         break;
265                 case LOG_TOTAL_MOMENTUM:
266                         moldyn->mwrite=timer;
267                         snprintf(filename,127,"%s/momentum",moldyn->vlsdir);
268                         moldyn->mfd=open(filename,
269                                          O_WRONLY|O_CREAT|O_EXCL,
270                                          S_IRUSR|S_IWUSR);
271                         if(moldyn->mfd<0) {
272                                 perror("[moldyn] momentum log fd open");
273                                 return moldyn->mfd;
274                         }
275                         dprintf(moldyn->efd,"# total momentum log file\n");
276                         printf("total momentum (%d)\n",timer);
277                         break;
278                 case LOG_PRESSURE:
279                         moldyn->pwrite=timer;
280                         snprintf(filename,127,"%s/pressure",moldyn->vlsdir);
281                         moldyn->pfd=open(filename,
282                                          O_WRONLY|O_CREAT|O_EXCL,
283                                          S_IRUSR|S_IWUSR);
284                         if(moldyn->pfd<0) {
285                                 perror("[moldyn] pressure log file\n");
286                                 return moldyn->pfd;
287                         }
288                         dprintf(moldyn->pfd,"# pressure log file\n");
289                         printf("pressure (%d)\n",timer);
290                         break;
291                 case LOG_TEMPERATURE:
292                         moldyn->twrite=timer;
293                         snprintf(filename,127,"%s/temperature",moldyn->vlsdir);
294                         moldyn->tfd=open(filename,
295                                          O_WRONLY|O_CREAT|O_EXCL,
296                                          S_IRUSR|S_IWUSR);
297                         if(moldyn->tfd<0) {
298                                 perror("[moldyn] temperature log file\n");
299                                 return moldyn->tfd;
300                         }
301                         dprintf(moldyn->tfd,"# temperature log file\n");
302                         printf("temperature (%d)\n",timer);
303                         break;
304                 case SAVE_STEP:
305                         moldyn->swrite=timer;
306                         printf("save file (%d)\n",timer);
307                         break;
308                 case VISUAL_STEP:
309                         moldyn->vwrite=timer;
310                         ret=visual_init(&(moldyn->vis),moldyn->vlsdir);
311                         if(ret<0) {
312                                 printf("[moldyn] visual init failure\n");
313                                 return ret;
314                         }
315                         printf("visual file (%d)\n",timer);
316                         break;
317                 case CREATE_REPORT:
318                         snprintf(filename,127,"%s/report.tex",moldyn->vlsdir);
319                         moldyn->rfd=open(filename,
320                                          O_WRONLY|O_CREAT|O_EXCL,
321                                          S_IRUSR|S_IWUSR);
322                         if(moldyn->rfd<0) {
323                                 perror("[moldyn] report fd open");      
324                                 return moldyn->rfd;
325                         }
326                         printf("report -> ");
327                         if(moldyn->efd) {
328                                 snprintf(filename,127,"%s/e_plot.scr",
329                                          moldyn->vlsdir);
330                                 moldyn->epfd=open(filename,
331                                                  O_WRONLY|O_CREAT|O_EXCL,
332                                                  S_IRUSR|S_IWUSR);
333                                 if(moldyn->epfd<0) {
334                                         perror("[moldyn] energy plot fd open");
335                                         return moldyn->epfd;
336                                 }
337                                 dprintf(moldyn->epfd,e_plot_script);
338                                 close(moldyn->epfd);
339                                 printf("energy ");
340                         }
341                         if(moldyn->pfd) {
342                                 snprintf(filename,127,"%s/pressure_plot.scr",
343                                          moldyn->vlsdir);
344                                 moldyn->ppfd=open(filename,
345                                                   O_WRONLY|O_CREAT|O_EXCL,
346                                                   S_IRUSR|S_IWUSR);
347                                 if(moldyn->ppfd<0) {
348                                         perror("[moldyn] p plot fd open");
349                                         return moldyn->ppfd;
350                                 }
351                                 dprintf(moldyn->ppfd,pressure_plot_script);
352                                 close(moldyn->ppfd);
353                                 printf("pressure ");
354                         }
355                         if(moldyn->tfd) {
356                                 snprintf(filename,127,"%s/temperature_plot.scr",
357                                          moldyn->vlsdir);
358                                 moldyn->tpfd=open(filename,
359                                                   O_WRONLY|O_CREAT|O_EXCL,
360                                                   S_IRUSR|S_IWUSR);
361                                 if(moldyn->tpfd<0) {
362                                         perror("[moldyn] t plot fd open");
363                                         return moldyn->tpfd;
364                                 }
365                                 dprintf(moldyn->tpfd,temperature_plot_script);
366                                 close(moldyn->tpfd);
367                                 printf("temperature ");
368                         }
369                         dprintf(moldyn->rfd,report_start,
370                                 moldyn->rauthor,moldyn->rtitle);
371                         printf("\n");
372                         break;
373                 default:
374                         printf("unknown log type: %02x\n",type);
375                         return -1;
376         }
377
378         return 0;
379 }
380
381 int moldyn_log_shutdown(t_moldyn *moldyn) {
382
383         char sc[256];
384
385         printf("[moldyn] log shutdown\n");
386         if(moldyn->efd) {
387                 close(moldyn->efd);
388                 if(moldyn->rfd) {
389                         dprintf(moldyn->rfd,report_energy);
390                         snprintf(sc,255,"cd %s && gnuplot e_plot.scr",
391                                  moldyn->vlsdir);
392                         system(sc);
393                 }
394         }
395         if(moldyn->mfd) close(moldyn->mfd);
396         if(moldyn->pfd) {
397                 close(moldyn->pfd);
398                 if(moldyn->rfd)
399                         dprintf(moldyn->rfd,report_pressure);
400                         snprintf(sc,255,"cd %s && gnuplot pressure_plot.scr",
401                                  moldyn->vlsdir);
402                         system(sc);
403         }
404         if(moldyn->tfd) {
405                 close(moldyn->tfd);
406                 if(moldyn->rfd)
407                         dprintf(moldyn->rfd,report_temperature);
408                         snprintf(sc,255,"cd %s && gnuplot temperature_plot.scr",
409                                  moldyn->vlsdir);
410                         system(sc);
411         }
412         if(moldyn->rfd) {
413                 dprintf(moldyn->rfd,report_end);
414                 close(moldyn->rfd);
415                 snprintf(sc,255,"cd %s && pdflatex report >/dev/null 2>&1",
416                          moldyn->vlsdir);
417                 system(sc);
418                 snprintf(sc,255,"cd %s && pdflatex report >/dev/null 2>&1",
419                          moldyn->vlsdir);
420                 system(sc);
421                 snprintf(sc,255,"cd %s && dvipdf report >/dev/null 2>&1",
422                          moldyn->vlsdir);
423                 system(sc);
424         }
425         if(&(moldyn->vis)) visual_tini(&(moldyn->vis));
426
427         return 0;
428 }
429
430 /*
431  * creating lattice functions
432  */
433
434 int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
435                    u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin) {
436
437         int new,count;
438         int ret;
439         t_3dvec orig;
440         void *ptr;
441         t_atom *atom;
442
443         new=a*b*c;
444         count=moldyn->count;
445
446         /* how many atoms do we expect */
447         if(type==CUBIC) new*=1;
448         if(type==FCC) new*=4;
449         if(type==DIAMOND) new*=8;
450
451         /* allocate space for atoms */
452         ptr=realloc(moldyn->atom,(count+new)*sizeof(t_atom));
453         if(!ptr) {
454                 perror("[moldyn] realloc (create lattice)");
455                 return -1;
456         }
457         moldyn->atom=ptr;
458         atom=&(moldyn->atom[count]);
459
460         /* no atoms on the boundaries (only reason: it looks better!) */
461         if(!origin) {
462                 orig.x=0.5*lc;
463                 orig.y=0.5*lc;
464                 orig.z=0.5*lc;
465         }
466         else {
467                 orig.x=origin->x;
468                 orig.y=origin->y;
469                 orig.z=origin->z;
470         }
471
472         switch(type) {
473                 case CUBIC:
474                         set_nn_dist(moldyn,lc);
475                         ret=cubic_init(a,b,c,lc,atom,&orig);
476                         break;
477                 case FCC:
478                         if(!origin)
479                                 v3_scale(&orig,&orig,0.5);
480                         set_nn_dist(moldyn,0.5*sqrt(2.0)*lc);
481                         ret=fcc_init(a,b,c,lc,atom,&orig);
482                         break;
483                 case DIAMOND:
484                         if(!origin)
485                                 v3_scale(&orig,&orig,0.25);
486                         set_nn_dist(moldyn,0.25*sqrt(3.0)*lc);
487                         ret=diamond_init(a,b,c,lc,atom,&orig);
488                         break;
489                 default:
490                         printf("unknown lattice type (%02x)\n",type);
491                         return -1;
492         }
493
494         /* debug */
495         if(ret!=new) {
496                 printf("[moldyn] creating lattice failed\n");
497                 printf("  amount of atoms\n");
498                 printf("  - expected: %d\n",new);
499                 printf("  - created: %d\n",ret);
500                 return -1;
501         }
502
503         moldyn->count+=new;
504         printf("[moldyn] created lattice with %d atoms\n",new);
505
506         for(ret=0;ret<new;ret++) {
507                 atom[ret].element=element;
508                 atom[ret].mass=mass;
509                 atom[ret].attr=attr;
510                 atom[ret].brand=brand;
511                 atom[ret].tag=count+ret;
512                 check_per_bound(moldyn,&(atom[ret].r));
513         }
514
515         /* update total system mass */
516         total_mass_calc(moldyn);
517
518         return ret;
519 }
520
521 /* cubic init */
522 int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
523
524         int count;
525         t_3dvec r;
526         int i,j,k;
527         t_3dvec o;
528
529         count=0;
530         if(origin)
531                 v3_copy(&o,origin);
532         else
533                 v3_zero(&o);
534
535         r.x=o.x;
536         for(i=0;i<a;i++) {
537                 r.y=o.y;
538                 for(j=0;j<b;j++) {
539                         r.z=o.z;
540                         for(k=0;k<c;k++) {
541                                 v3_copy(&(atom[count].r),&r);
542                                 count+=1;
543                                 r.z+=lc;
544                         }
545                         r.y+=lc;
546                 }
547                 r.x+=lc;
548         }
549
550         for(i=0;i<count;i++) {
551                 atom[i].r.x-=(a*lc)/2.0;
552                 atom[i].r.y-=(b*lc)/2.0;
553                 atom[i].r.z-=(c*lc)/2.0;
554         }
555
556         return count;
557 }
558
559 /* fcc lattice init */
560 int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
561
562         int count;
563         int i,j,k,l;
564         t_3dvec o,r,n;
565         t_3dvec basis[3];
566
567         count=0;
568         if(origin)
569                 v3_copy(&o,origin);
570         else
571                 v3_zero(&o);
572
573         /* construct the basis */
574         memset(basis,0,3*sizeof(t_3dvec));
575         basis[0].x=0.5*lc;
576         basis[0].y=0.5*lc;
577         basis[1].x=0.5*lc;
578         basis[1].z=0.5*lc;
579         basis[2].y=0.5*lc;
580         basis[2].z=0.5*lc;
581
582         /* fill up the room */
583         r.x=o.x;
584         for(i=0;i<a;i++) {
585                 r.y=o.y;
586                 for(j=0;j<b;j++) {
587                         r.z=o.z;
588                         for(k=0;k<c;k++) {
589                                 /* first atom */
590                                 v3_copy(&(atom[count].r),&r);
591                                 count+=1;
592                                 r.z+=lc;
593                                 /* the three face centered atoms */
594                                 for(l=0;l<3;l++) {
595                                         v3_add(&n,&r,&basis[l]);
596                                         v3_copy(&(atom[count].r),&n);
597                                         count+=1;
598                                 }
599                         }
600                         r.y+=lc;
601                 }
602                 r.x+=lc;
603         }
604                                 
605         /* coordinate transformation */
606         for(i=0;i<count;i++) {
607                 atom[i].r.x-=(a*lc)/2.0;
608                 atom[i].r.y-=(b*lc)/2.0;
609                 atom[i].r.z-=(c*lc)/2.0;
610         }
611
612         return count;
613 }
614
615 int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
616
617         int count;
618         t_3dvec o;
619
620         count=fcc_init(a,b,c,lc,atom,origin);
621
622         o.x=0.25*lc;
623         o.y=0.25*lc;
624         o.z=0.25*lc;
625
626         if(origin) v3_add(&o,&o,origin);
627
628         count+=fcc_init(a,b,c,lc,&atom[count],&o);
629
630         return count;
631 }
632
633 int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
634              t_3dvec *r,t_3dvec *v) {
635
636         t_atom *atom;
637         void *ptr;
638         int count;
639         
640         atom=moldyn->atom;
641         count=(moldyn->count)++;
642
643         ptr=realloc(atom,(count+1)*sizeof(t_atom));
644         if(!ptr) {
645                 perror("[moldyn] realloc (add atom)");
646                 return -1;
647         }
648         moldyn->atom=ptr;
649
650         atom=moldyn->atom;
651         atom[count].r=*r;
652         atom[count].v=*v;
653         atom[count].element=element;
654         atom[count].mass=mass;
655         atom[count].brand=brand;
656         atom[count].tag=count;
657         atom[count].attr=attr;
658
659         /* update total system mass */
660         total_mass_calc(moldyn);
661
662         return 0;
663 }
664
665 int destroy_atoms(t_moldyn *moldyn) {
666
667         if(moldyn->atom) free(moldyn->atom);
668
669         return 0;
670 }
671
672 int thermal_init(t_moldyn *moldyn,u8 equi_init) {
673
674         /*
675          * - gaussian distribution of velocities
676          * - zero total momentum
677          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
678          */
679
680         int i;
681         double v,sigma;
682         t_3dvec p_total,delta;
683         t_atom *atom;
684         t_random *random;
685
686         atom=moldyn->atom;
687         random=&(moldyn->random);
688
689         printf("[moldyn] thermal init (equi init: %s)\n",equi_init?"yes":"no");
690
691         /* gaussian distribution of velocities */
692         v3_zero(&p_total);
693         for(i=0;i<moldyn->count;i++) {
694                 sigma=sqrt(2.0*K_BOLTZMANN*moldyn->t_ref/atom[i].mass);
695                 /* x direction */
696                 v=sigma*rand_get_gauss(random);
697                 atom[i].v.x=v;
698                 p_total.x+=atom[i].mass*v;
699                 /* y direction */
700                 v=sigma*rand_get_gauss(random);
701                 atom[i].v.y=v;
702                 p_total.y+=atom[i].mass*v;
703                 /* z direction */
704                 v=sigma*rand_get_gauss(random);
705                 atom[i].v.z=v;
706                 p_total.z+=atom[i].mass*v;
707         }
708
709         /* zero total momentum */
710         v3_scale(&p_total,&p_total,1.0/moldyn->count);
711         for(i=0;i<moldyn->count;i++) {
712                 v3_scale(&delta,&p_total,1.0/atom[i].mass);
713                 v3_sub(&(atom[i].v),&(atom[i].v),&delta);
714         }
715
716         /* velocity scaling */
717         scale_velocity(moldyn,equi_init);
718
719         return 0;
720 }
721
722 double total_mass_calc(t_moldyn *moldyn) {
723
724         int i;
725
726         moldyn->mass=0.0;
727
728         for(i=0;i<moldyn->count;i++)
729                 moldyn->mass+=moldyn->atom[i].mass;
730
731         return moldyn->mass;
732 }
733
734 double temperature_calc(t_moldyn *moldyn) {
735
736         /* assume up to date kinetic energy, which is 3/2 N k_B T */
737
738         moldyn->t=(2.0*moldyn->ekin)/(3.0*K_BOLTZMANN*moldyn->count);
739
740         return moldyn->t;
741 }
742
743 double get_temperature(t_moldyn *moldyn) {
744
745         return moldyn->t;
746 }
747
748 int scale_velocity(t_moldyn *moldyn,u8 equi_init) {
749
750         int i;
751         double e,scale;
752         t_atom *atom;
753         int count;
754
755         atom=moldyn->atom;
756
757         /*
758          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
759          */
760
761         /* get kinetic energy / temperature & count involved atoms */
762         e=0.0;
763         count=0;
764         for(i=0;i<moldyn->count;i++) {
765                 if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB)) {
766                         e+=atom[i].mass*v3_absolute_square(&(atom[i].v));
767                         count+=1;
768                 }
769         }
770         e*=0.5;
771         if(count!=0) moldyn->t=e/(1.5*count*K_BOLTZMANN);
772         else return 0;  /* no atoms involved in scaling! */
773         
774         /* (temporary) hack for e,t = 0 */
775         if(e==0.0) {
776         moldyn->t=0.0;
777                 if(moldyn->t_ref!=0.0) {
778                         thermal_init(moldyn,equi_init);
779                         return 0;
780                 }
781                 else
782                         return 0; /* no scaling needed */
783         }
784
785
786         /* get scaling factor */
787         scale=moldyn->t_ref/moldyn->t;
788         if(equi_init&TRUE)
789                 scale*=2.0;
790         else
791                 if(moldyn->pt_scale&T_SCALE_BERENDSEN)
792                         scale=1.0+(scale-1.0)/moldyn->t_tc;
793         scale=sqrt(scale);
794
795         /* velocity scaling */
796         for(i=0;i<moldyn->count;i++) {
797                 if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB))
798                         v3_scale(&(atom[i].v),&(atom[i].v),scale);
799         }
800
801         return 0;
802 }
803
804 double ideal_gas_law_pressure(t_moldyn *moldyn) {
805
806         double p;
807
808         p=moldyn->count*moldyn->t*K_BOLTZMANN/moldyn->volume;
809
810         return p;
811 }
812
813 double virial_sum(t_moldyn *moldyn) {
814
815         int i;
816         double v;
817         t_virial *virial;
818
819         /* virial (sum over atom virials) */
820         v=0.0;
821         for(i=0;i<moldyn->count;i++) {
822                 virial=&(moldyn->atom[i].virial);
823                 v+=(virial->xx+virial->yy+virial->zz);
824         }
825         moldyn->virial=v;
826
827         /* global virial (absolute coordinates) */
828         virial=&(moldyn->gvir);
829         moldyn->gv=virial->xx+virial->yy+virial->zz;
830
831         return moldyn->virial;
832 }
833
834 double pressure_calc(t_moldyn *moldyn) {
835
836         /*
837          * PV = NkT + <W>
838          * with W = 1/3 sum_i f_i r_i (- skipped!)
839          * virial = sum_i f_i r_i
840          * 
841          * => P = (2 Ekin + virial) / (3V)
842          */
843
844         /* assume up to date virial & up to date kinetic energy */
845
846         /* pressure (atom virials) */
847         moldyn->p=2.0*moldyn->ekin+moldyn->virial;
848         moldyn->p/=(3.0*moldyn->volume);
849
850         /* pressure (absolute coordinates) */
851         moldyn->gp=2.0*moldyn->ekin+moldyn->gv;
852         moldyn->gp/=(3.0*moldyn->volume);
853
854         return moldyn->p;
855 }
856
857 int average_and_fluctuation_calc(t_moldyn *moldyn) {
858
859         if(moldyn->total_steps<moldyn->avg_skip)
860                 return 0;
861
862         int denom=moldyn->total_steps+1-moldyn->avg_skip;
863
864         /* assume up to date energies, temperature, pressure etc */
865
866         /* kinetic energy */
867         moldyn->k_sum+=moldyn->ekin;
868         moldyn->k2_sum+=(moldyn->ekin*moldyn->ekin);
869         moldyn->k_avg=moldyn->k_sum/denom;
870         moldyn->k2_avg=moldyn->k2_sum/denom;
871         moldyn->dk2_avg=moldyn->k2_avg-(moldyn->k_avg*moldyn->k_avg);
872
873         /* potential energy */
874         moldyn->v_sum+=moldyn->energy;
875         moldyn->v2_sum+=(moldyn->energy*moldyn->energy);
876         moldyn->v_avg=moldyn->v_sum/denom;
877         moldyn->v2_avg=moldyn->v2_sum/denom;
878         moldyn->dv2_avg=moldyn->v2_avg-(moldyn->v_avg*moldyn->v_avg);
879
880         /* temperature */
881         moldyn->t_sum+=moldyn->t;
882         moldyn->t_avg=moldyn->t_sum/denom;
883
884         /* virial */
885         moldyn->virial_sum+=moldyn->virial;
886         moldyn->virial_avg=moldyn->virial_sum/denom;
887         moldyn->gv_sum+=moldyn->gv;
888         moldyn->gv_avg=moldyn->gv_sum/denom;
889
890         /* pressure */
891         moldyn->p_sum+=moldyn->p;
892         moldyn->p_avg=moldyn->p_sum/denom;
893         moldyn->gp_sum+=moldyn->gp;
894         moldyn->gp_avg=moldyn->gp_sum/denom;
895
896         return 0;
897 }
898
899 int get_heat_capacity(t_moldyn *moldyn) {
900
901         double temp2,ighc;
902
903         /* averages needed for heat capacity calc */
904         if(moldyn->total_steps<moldyn->avg_skip)
905                 return 0;
906
907         /* (temperature average)^2 */
908         temp2=moldyn->t_avg*moldyn->t_avg;
909         printf("[moldyn] specific heat capacity for T=%f K [J/(kg K)]\n",
910                moldyn->t_avg);
911
912         /* ideal gas contribution */
913         ighc=3.0*moldyn->count*K_BOLTZMANN/2.0;
914         printf("  ideal gas contribution: %f\n",
915                ighc/moldyn->mass*KILOGRAM/JOULE);
916
917         /* specific heat for nvt ensemble */
918         moldyn->c_v_nvt=moldyn->dv2_avg/(K_BOLTZMANN*temp2)+ighc;
919         moldyn->c_v_nvt/=moldyn->mass;
920
921         /* specific heat for nve ensemble */
922         moldyn->c_v_nve=ighc/(1.0-(moldyn->dv2_avg/(ighc*K_BOLTZMANN*temp2)));
923         moldyn->c_v_nve/=moldyn->mass;
924
925         printf("  NVE: %f\n",moldyn->c_v_nve*KILOGRAM/JOULE);
926         printf("  NVT: %f\n",moldyn->c_v_nvt*KILOGRAM/JOULE);
927 printf("  --> <dV2> sim: %f experimental: %f\n",moldyn->dv2_avg,1.5*moldyn->count*K_B2*moldyn->t_avg*moldyn->t_avg*(1.0-1.5*moldyn->count*K_BOLTZMANN/(700*moldyn->mass*JOULE/KILOGRAM)));
928
929         return 0;
930 }
931
932 double thermodynamic_pressure_calc(t_moldyn *moldyn) {
933
934         t_3dvec dim,*tp;
935         double u_up,u_down,dv;
936         double scale,p;
937         t_atom *store;
938
939         /*
940          * dU = - p dV
941          *
942          * => p = - dU/dV
943          *
944          */
945
946         scale=0.00001;
947         dv=8*scale*scale*scale*moldyn->volume;
948
949         store=malloc(moldyn->count*sizeof(t_atom));
950         if(store==NULL) {
951                 printf("[moldyn] allocating store mem failed\n");
952                 return -1;
953         }
954
955         /* save unscaled potential energy + atom/dim configuration */
956         memcpy(store,moldyn->atom,moldyn->count*sizeof(t_atom));
957         dim=moldyn->dim;
958
959         /* scale up dimension and atom positions */
960         scale_dim(moldyn,SCALE_UP,scale,TRUE,TRUE,TRUE);
961         scale_atoms(moldyn,SCALE_UP,scale,TRUE,TRUE,TRUE);
962         link_cell_shutdown(moldyn);
963         link_cell_init(moldyn,QUIET);
964         potential_force_calc(moldyn);
965         u_up=moldyn->energy;
966
967         /* restore atomic configuration + dim */
968         memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
969         moldyn->dim=dim;
970
971         /* scale down dimension and atom positions */
972         scale_dim(moldyn,SCALE_DOWN,scale,TRUE,TRUE,TRUE);
973         scale_atoms(moldyn,SCALE_DOWN,scale,TRUE,TRUE,TRUE);
974         link_cell_shutdown(moldyn);
975         link_cell_init(moldyn,QUIET);
976         potential_force_calc(moldyn);
977         u_down=moldyn->energy;
978         
979         /* calculate pressure */
980         p=-(u_up-u_down)/dv;
981 printf("-------> %.10f %.10f %f\n",u_up/EV/moldyn->count,u_down/EV/moldyn->count,p/BAR);
982
983         /* restore atomic configuration + dim */
984         memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
985         moldyn->dim=dim;
986
987         /* restore energy */
988         potential_force_calc(moldyn);
989
990         link_cell_shutdown(moldyn);
991         link_cell_init(moldyn,QUIET);
992
993         return p;
994 }
995
996 double get_pressure(t_moldyn *moldyn) {
997
998         return moldyn->p;
999
1000 }
1001
1002 int scale_dim(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z) {
1003
1004         t_3dvec *dim;
1005
1006         dim=&(moldyn->dim);
1007
1008         if(dir==SCALE_UP)
1009                 scale=1.0+scale;
1010
1011         if(dir==SCALE_DOWN)
1012                 scale=1.0-scale;
1013
1014         if(x) dim->x*=scale;
1015         if(y) dim->y*=scale;
1016         if(z) dim->z*=scale;
1017
1018         return 0;
1019 }
1020
1021 int scale_atoms(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z) {
1022
1023         int i;
1024         t_3dvec *r;
1025
1026         if(dir==SCALE_UP)
1027                 scale=1.0+scale;
1028
1029         if(dir==SCALE_DOWN)
1030                 scale=1.0-scale;
1031
1032         for(i=0;i<moldyn->count;i++) {
1033                 r=&(moldyn->atom[i].r);
1034                 if(x) r->x*=scale;
1035                 if(y) r->y*=scale;
1036                 if(z) r->z*=scale;
1037         }
1038
1039         return 0;
1040 }
1041
1042 int scale_volume(t_moldyn *moldyn) {
1043
1044         t_3dvec *dim,*vdim;
1045         double scale;
1046         t_linkcell *lc;
1047
1048         vdim=&(moldyn->vis.dim);
1049         dim=&(moldyn->dim);
1050         lc=&(moldyn->lc);
1051
1052         /* scaling factor */
1053         if(moldyn->pt_scale&P_SCALE_BERENDSEN) {
1054                 scale=1.0-(moldyn->p_ref-moldyn->p)/moldyn->p_tc;
1055                 scale=pow(scale,ONE_THIRD);
1056         }
1057         else {
1058                 scale=pow(moldyn->p/moldyn->p_ref,ONE_THIRD);
1059         }
1060 moldyn->debug=scale;
1061
1062         /* scale the atoms and dimensions */
1063         scale_atoms(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
1064         scale_dim(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
1065
1066         /* visualize dimensions */
1067         if(vdim->x!=0) {
1068                 vdim->x=dim->x;
1069                 vdim->y=dim->y;
1070                 vdim->z=dim->z;
1071         }
1072
1073         /* recalculate scaled volume */
1074         moldyn->volume=dim->x*dim->y*dim->z;
1075
1076         /* adjust/reinit linkcell */
1077         if(((int)(dim->x/moldyn->cutoff)!=lc->nx)||
1078            ((int)(dim->y/moldyn->cutoff)!=lc->ny)||
1079            ((int)(dim->z/moldyn->cutoff)!=lc->nx)) {
1080                 link_cell_shutdown(moldyn);
1081                 link_cell_init(moldyn,QUIET);
1082         } else {
1083                 lc->x*=scale;
1084                 lc->y*=scale;
1085                 lc->z*=scale;
1086         }
1087
1088         return 0;
1089
1090 }
1091
1092 double e_kin_calc(t_moldyn *moldyn) {
1093
1094         int i;
1095         t_atom *atom;
1096
1097         atom=moldyn->atom;
1098         moldyn->ekin=0.0;
1099
1100         for(i=0;i<moldyn->count;i++)
1101                 moldyn->ekin+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
1102
1103         return moldyn->ekin;
1104 }
1105
1106 double get_total_energy(t_moldyn *moldyn) {
1107
1108         return(moldyn->ekin+moldyn->energy);
1109 }
1110
1111 t_3dvec get_total_p(t_moldyn *moldyn) {
1112
1113         t_3dvec p,p_total;
1114         int i;
1115         t_atom *atom;
1116
1117         atom=moldyn->atom;
1118
1119         v3_zero(&p_total);
1120         for(i=0;i<moldyn->count;i++) {
1121                 v3_scale(&p,&(atom[i].v),atom[i].mass);
1122                 v3_add(&p_total,&p_total,&p);
1123         }
1124
1125         return p_total;
1126 }
1127
1128 double estimate_time_step(t_moldyn *moldyn,double nn_dist) {
1129
1130         double tau;
1131
1132         /* nn_dist is the nearest neighbour distance */
1133
1134         tau=(0.05*nn_dist*moldyn->atom[0].mass)/sqrt(3.0*K_BOLTZMANN*moldyn->t);
1135
1136         return tau;     
1137 }
1138
1139 /*
1140  * numerical tricks
1141  */
1142
1143 /* linked list / cell method */
1144
1145 int link_cell_init(t_moldyn *moldyn,u8 vol) {
1146
1147         t_linkcell *lc;
1148         int i;
1149
1150         lc=&(moldyn->lc);
1151
1152         /* partitioning the md cell */
1153         lc->nx=moldyn->dim.x/moldyn->cutoff;
1154         lc->x=moldyn->dim.x/lc->nx;
1155         lc->ny=moldyn->dim.y/moldyn->cutoff;
1156         lc->y=moldyn->dim.y/lc->ny;
1157         lc->nz=moldyn->dim.z/moldyn->cutoff;
1158         lc->z=moldyn->dim.z/lc->nz;
1159
1160         lc->cells=lc->nx*lc->ny*lc->nz;
1161         lc->subcell=malloc(lc->cells*sizeof(t_list));
1162
1163         if(lc->cells<27)
1164                 printf("[moldyn] FATAL: less then 27 subcells!\n");
1165
1166         if(vol) {
1167                 printf("[moldyn] initializing linked cells (%d)\n",lc->cells);
1168                 printf("  x: %d x %f A\n",lc->nx,lc->x);
1169                 printf("  y: %d x %f A\n",lc->ny,lc->y);
1170                 printf("  z: %d x %f A\n",lc->nz,lc->z);
1171         }
1172
1173         for(i=0;i<lc->cells;i++)
1174                 list_init_f(&(lc->subcell[i]));
1175
1176         link_cell_update(moldyn);
1177         
1178         return 0;
1179 }
1180
1181 int link_cell_update(t_moldyn *moldyn) {
1182
1183         int count,i,j,k;
1184         int nx,ny;
1185         t_atom *atom;
1186         t_linkcell *lc;
1187         double x,y,z;
1188
1189         atom=moldyn->atom;
1190         lc=&(moldyn->lc);
1191
1192         nx=lc->nx;
1193         ny=lc->ny;
1194
1195         x=moldyn->dim.x/2;
1196         y=moldyn->dim.y/2;
1197         z=moldyn->dim.z/2;
1198
1199         for(i=0;i<lc->cells;i++)
1200                 list_destroy_f(&(lc->subcell[i]));
1201         
1202         for(count=0;count<moldyn->count;count++) {
1203                 i=((atom[count].r.x+(moldyn->dim.x/2))/lc->x);
1204                 j=((atom[count].r.y+(moldyn->dim.y/2))/lc->y);
1205                 k=((atom[count].r.z+(moldyn->dim.z/2))/lc->z);
1206                 list_add_immediate_f(&(lc->subcell[i+j*nx+k*nx*ny]),
1207                                      &(atom[count]));
1208         }
1209
1210         return 0;
1211 }
1212
1213 int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) {
1214
1215         t_linkcell *lc;
1216         int a;
1217         int count1,count2;
1218         int ci,cj,ck;
1219         int nx,ny,nz;
1220         int x,y,z;
1221         u8 bx,by,bz;
1222
1223         lc=&(moldyn->lc);
1224         nx=lc->nx;
1225         ny=lc->ny;
1226         nz=lc->nz;
1227         count1=1;
1228         count2=27;
1229         a=nx*ny;
1230
1231         cell[0]=lc->subcell[i+j*nx+k*a];
1232         for(ci=-1;ci<=1;ci++) {
1233                 bx=0;
1234                 x=i+ci;
1235                 if((x<0)||(x>=nx)) {
1236                         x=(x+nx)%nx;
1237                         bx=1;
1238                 }
1239                 for(cj=-1;cj<=1;cj++) {
1240                         by=0;
1241                         y=j+cj;
1242                         if((y<0)||(y>=ny)) {
1243                                 y=(y+ny)%ny;
1244                                 by=1;
1245                         }
1246                         for(ck=-1;ck<=1;ck++) {
1247                                 bz=0;
1248                                 z=k+ck;
1249                                 if((z<0)||(z>=nz)) {
1250                                         z=(z+nz)%nz;
1251                                         bz=1;
1252                                 }
1253                                 if(!(ci|cj|ck)) continue;
1254                                 if(bx|by|bz) {
1255                                         cell[--count2]=lc->subcell[x+y*nx+z*a];
1256                                 }
1257                                 else {
1258                                         cell[count1++]=lc->subcell[x+y*nx+z*a];
1259                                 }
1260                         }
1261                 }
1262         }
1263
1264         lc->dnlc=count1;
1265
1266         return count1;
1267 }
1268
1269 int link_cell_shutdown(t_moldyn *moldyn) {
1270
1271         int i;
1272         t_linkcell *lc;
1273
1274         lc=&(moldyn->lc);
1275
1276         for(i=0;i<lc->nx*lc->ny*lc->nz;i++)
1277                 list_destroy_f(&(moldyn->lc.subcell[i]));
1278
1279         free(lc->subcell);
1280
1281         return 0;
1282 }
1283
1284 int moldyn_add_schedule(t_moldyn *moldyn,int runs,double tau) {
1285
1286         int count;
1287         void *ptr;
1288         t_moldyn_schedule *schedule;
1289
1290         schedule=&(moldyn->schedule);
1291         count=++(schedule->total_sched);
1292
1293         ptr=realloc(schedule->runs,count*sizeof(int));
1294         if(!ptr) {
1295                 perror("[moldyn] realloc (runs)");
1296                 return -1;
1297         }
1298         schedule->runs=ptr;
1299         schedule->runs[count-1]=runs;
1300
1301         ptr=realloc(schedule->tau,count*sizeof(double));
1302         if(!ptr) {
1303                 perror("[moldyn] realloc (tau)");
1304                 return -1;
1305         }
1306         schedule->tau=ptr;
1307         schedule->tau[count-1]=tau;
1308
1309         printf("[moldyn] schedule added:\n");
1310         printf("  number: %d | runs: %d | tau: %f\n",count-1,runs,tau);
1311                                        
1312
1313         return 0;
1314 }
1315
1316 int moldyn_set_schedule_hook(t_moldyn *moldyn,set_hook hook,void *hook_params) {
1317
1318         moldyn->schedule.hook=hook;
1319         moldyn->schedule.hook_params=hook_params;
1320         
1321         return 0;
1322 }
1323
1324 /*
1325  *
1326  * 'integration of newtons equation' - algorithms
1327  *
1328  */
1329
1330 /* start the integration */
1331
1332 int moldyn_integrate(t_moldyn *moldyn) {
1333
1334         int i;
1335         unsigned int e,m,s,v,p,t;
1336         t_3dvec momentum;
1337         t_moldyn_schedule *sched;
1338         t_atom *atom;
1339         int fd;
1340         char dir[128];
1341         double ds;
1342         double energy_scale;
1343         //double tp;
1344
1345         sched=&(moldyn->schedule);
1346         atom=moldyn->atom;
1347
1348         /* initialize linked cell method */
1349         link_cell_init(moldyn,VERBOSE);
1350
1351         /* logging & visualization */
1352         e=moldyn->ewrite;
1353         m=moldyn->mwrite;
1354         s=moldyn->swrite;
1355         v=moldyn->vwrite;
1356         p=moldyn->pwrite;
1357         t=moldyn->twrite;
1358
1359         /* sqaure of some variables */
1360         moldyn->tau_square=moldyn->tau*moldyn->tau;
1361         moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff;
1362
1363         /* energy scaling factor */
1364         energy_scale=moldyn->count*EV;
1365
1366         /* calculate initial forces */
1367         potential_force_calc(moldyn);
1368 #ifdef DEBUG
1369 return 0;
1370 #endif
1371
1372         /* some stupid checks before we actually start calculating bullshit */
1373         if(moldyn->cutoff>0.5*moldyn->dim.x)
1374                 printf("[moldyn] warning: cutoff > 0.5 x dim.x\n");
1375         if(moldyn->cutoff>0.5*moldyn->dim.y)
1376                 printf("[moldyn] warning: cutoff > 0.5 x dim.y\n");
1377         if(moldyn->cutoff>0.5*moldyn->dim.z)
1378                 printf("[moldyn] warning: cutoff > 0.5 x dim.z\n");
1379         ds=0.5*atom[0].f.x*moldyn->tau_square/atom[0].mass;
1380         if(ds>0.05*moldyn->nnd)
1381                 printf("[moldyn] warning: forces too high / tau too small!\n");
1382
1383         /* zero absolute time */
1384         moldyn->time=0.0;
1385         moldyn->total_steps=0;
1386
1387         /* debugging, ignore */
1388         moldyn->debug=0;
1389
1390         /* tell the world */
1391         printf("[moldyn] integration start, go get a coffee ...\n");
1392
1393         /* executing the schedule */
1394         sched->count=0;
1395         while(sched->count<sched->total_sched) {
1396
1397                 /* setting amount of runs and finite time step size */
1398                 moldyn->tau=sched->tau[sched->count];
1399                 moldyn->tau_square=moldyn->tau*moldyn->tau;
1400                 moldyn->time_steps=sched->runs[sched->count];
1401
1402         /* integration according to schedule */
1403
1404         for(i=0;i<moldyn->time_steps;i++) {
1405
1406                 /* integration step */
1407                 moldyn->integrate(moldyn);
1408
1409                 /* calculate kinetic energy, temperature and pressure */
1410                 e_kin_calc(moldyn);
1411                 temperature_calc(moldyn);
1412                 virial_sum(moldyn);
1413                 pressure_calc(moldyn);
1414                 average_and_fluctuation_calc(moldyn);
1415
1416                 /* p/t scaling */
1417                 if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT))
1418                         scale_velocity(moldyn,FALSE);
1419                 if(moldyn->pt_scale&(P_SCALE_BERENDSEN|P_SCALE_DIRECT))
1420                         scale_volume(moldyn);
1421
1422                 /* check for log & visualization */
1423                 if(e) {
1424                         if(!(i%e))
1425                                 dprintf(moldyn->efd,
1426                                         "%f %f %f %f\n",
1427                                         moldyn->time,moldyn->ekin/energy_scale,
1428                                         moldyn->energy/energy_scale,
1429                                         get_total_energy(moldyn)/energy_scale);
1430                 }
1431                 if(m) {
1432                         if(!(i%m)) {
1433                                 momentum=get_total_p(moldyn);
1434                                 dprintf(moldyn->mfd,
1435                                         "%f %f %f %f %f\n",moldyn->time,
1436                                         momentum.x,momentum.y,momentum.z,
1437                                         v3_norm(&momentum));
1438                         }
1439                 }
1440                 if(p) {
1441                         if(!(i%p)) {
1442                                 dprintf(moldyn->pfd,
1443                                         "%f %f %f %f %f\n",moldyn->time,
1444                                          moldyn->p/BAR,moldyn->p_avg/BAR,
1445                                          moldyn->gp/BAR,moldyn->gp_avg/BAR);
1446                         }
1447                 }
1448                 if(t) {
1449                         if(!(i%t)) {
1450                                 dprintf(moldyn->tfd,
1451                                         "%f %f %f\n",
1452                                         moldyn->time,moldyn->t,moldyn->t_avg);
1453                         }
1454                 }
1455                 if(s) {
1456                         if(!(i%s)) {
1457                                 snprintf(dir,128,"%s/s-%07.f.save",
1458                                          moldyn->vlsdir,moldyn->time);
1459                                 fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT);
1460                                 if(fd<0) perror("[moldyn] save fd open");
1461                                 else {
1462                                         write(fd,moldyn,sizeof(t_moldyn));
1463                                         write(fd,moldyn->atom,
1464                                               moldyn->count*sizeof(t_atom));
1465                                 }
1466                                 close(fd);
1467                         }       
1468                 }
1469                 if(v) {
1470                         if(!(i%v)) {
1471                                 visual_atoms(&(moldyn->vis),moldyn->time,
1472                                              moldyn->atom,moldyn->count);
1473                         }
1474                 }
1475
1476                 /* display progress */
1477                 if(!(i%10)) {
1478         printf("\rsched:%d, steps:%d, T:%3.1f/%3.1f P:%4.1f/%4.1f V:%6.1f",
1479                sched->count,i,
1480                moldyn->t,moldyn->t_avg,
1481                moldyn->p_avg/BAR,moldyn->p/BAR,
1482                moldyn->volume);
1483         fflush(stdout);
1484                 }
1485
1486                 /* increase absolute time */
1487                 moldyn->time+=moldyn->tau;
1488                 moldyn->total_steps+=1;
1489
1490         }
1491
1492                 /* check for hooks */
1493                 if(sched->count+1<sched->total_sched) {
1494                         if(sched->hook) {
1495                                 printf("\n ## schedule hook %d/%d start ##\n",
1496                                        sched->count+1,sched->total_sched-1);
1497                                 sched->hook(moldyn,sched->hook_params);
1498                                 printf(" ## schedule hook end ##\n");
1499                         }
1500                 }
1501
1502                 /* increase the schedule counter */
1503                 sched->count+=1;
1504
1505         }
1506
1507         return 0;
1508 }
1509
1510 /* velocity verlet */
1511
1512 int velocity_verlet(t_moldyn *moldyn) {
1513
1514         int i,count;
1515         double tau,tau_square,h;
1516         t_3dvec delta;
1517         t_atom *atom;
1518
1519         atom=moldyn->atom;
1520         count=moldyn->count;
1521         tau=moldyn->tau;
1522         tau_square=moldyn->tau_square;
1523
1524         for(i=0;i<count;i++) {
1525                 /* new positions */
1526                 h=0.5/atom[i].mass;
1527                 v3_scale(&delta,&(atom[i].v),tau);
1528                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
1529                 v3_scale(&delta,&(atom[i].f),h*tau_square);
1530                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
1531                 check_per_bound(moldyn,&(atom[i].r));
1532
1533                 /* velocities [actually v(t+tau/2)] */
1534                 v3_scale(&delta,&(atom[i].f),h*tau);
1535                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
1536         }
1537
1538         /* neighbour list update */
1539         link_cell_update(moldyn);
1540
1541         /* forces depending on chosen potential */
1542         potential_force_calc(moldyn);
1543
1544         for(i=0;i<count;i++) {
1545                 /* again velocities [actually v(t+tau)] */
1546                 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
1547                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
1548         }
1549
1550         return 0;
1551 }
1552
1553
1554 /*
1555  *
1556  * potentials & corresponding forces & virial routine
1557  * 
1558  */
1559
1560 /* generic potential and force calculation */
1561
1562 int potential_force_calc(t_moldyn *moldyn) {
1563
1564         int i,j,k,count;
1565         t_atom *itom,*jtom,*ktom;
1566         t_virial *virial;
1567         t_linkcell *lc;
1568         t_list neighbour_i[27];
1569         t_list neighbour_i2[27];
1570         t_list *this,*that;
1571         u8 bc_ij,bc_ik;
1572         int dnlc;
1573
1574         count=moldyn->count;
1575         itom=moldyn->atom;
1576         lc=&(moldyn->lc);
1577
1578         /* reset energy */
1579         moldyn->energy=0.0;
1580
1581         /* reset global virial */
1582         memset(&(moldyn->gvir),0,sizeof(t_virial));
1583
1584         /* reset force, site energy and virial of every atom */
1585         for(i=0;i<count;i++) {
1586
1587                 /* reset force */
1588                 v3_zero(&(itom[i].f));
1589
1590                 /* reset virial */
1591                 virial=(&(itom[i].virial));
1592                 virial->xx=0.0;
1593                 virial->yy=0.0;
1594                 virial->zz=0.0;
1595                 virial->xy=0.0;
1596                 virial->xz=0.0;
1597                 virial->yz=0.0;
1598         
1599                 /* reset site energy */
1600                 itom[i].e=0.0;
1601
1602         }
1603
1604         /* get energy, force and virial of every atom */
1605
1606         /* first (and only) loop over atoms i */
1607         for(i=0;i<count;i++) {
1608
1609                 /* single particle potential/force */
1610                 if(itom[i].attr&ATOM_ATTR_1BP)
1611                         if(moldyn->func1b)
1612                                 moldyn->func1b(moldyn,&(itom[i]));
1613
1614                 if(!(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)))
1615                         continue;
1616
1617                 /* 2 body pair potential/force */
1618         
1619                 link_cell_neighbour_index(moldyn,
1620                                           (itom[i].r.x+moldyn->dim.x/2)/lc->x,
1621                                           (itom[i].r.y+moldyn->dim.y/2)/lc->y,
1622                                           (itom[i].r.z+moldyn->dim.z/2)/lc->z,
1623                                           neighbour_i);
1624
1625                 dnlc=lc->dnlc;
1626
1627                 /* first loop over atoms j */
1628                 if(moldyn->func2b) {
1629                         for(j=0;j<27;j++) {
1630
1631                                 this=&(neighbour_i[j]);
1632                                 list_reset_f(this);
1633
1634                                 if(this->start==NULL)
1635                                         continue;
1636
1637                                 bc_ij=(j<dnlc)?0:1;
1638
1639                                 do {
1640                                         jtom=this->current->data;
1641
1642                                         if(jtom==&(itom[i]))
1643                                                 continue;
1644
1645                                         if((jtom->attr&ATOM_ATTR_2BP)&
1646                                            (itom[i].attr&ATOM_ATTR_2BP)) {
1647                                                 moldyn->func2b(moldyn,
1648                                                                &(itom[i]),
1649                                                                jtom,
1650                                                                bc_ij);
1651                                         }
1652                                 } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
1653
1654                         }
1655                 }
1656
1657                 /* 3 body potential/force */
1658
1659                 if(!(itom[i].attr&ATOM_ATTR_3BP))
1660                         continue;
1661
1662                 /* copy the neighbour lists */
1663                 memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
1664
1665                 /* second loop over atoms j */
1666                 for(j=0;j<27;j++) {
1667
1668                         this=&(neighbour_i[j]);
1669                         list_reset_f(this);
1670
1671                         if(this->start==NULL)
1672                                 continue;
1673
1674                         bc_ij=(j<dnlc)?0:1;
1675
1676                         do {
1677                                 jtom=this->current->data;
1678
1679                                 if(jtom==&(itom[i]))
1680                                         continue;
1681
1682                                 if(!(jtom->attr&ATOM_ATTR_3BP))
1683                                         continue;
1684
1685                                 /* reset 3bp run */
1686                                 moldyn->run3bp=1;
1687
1688                                 if(moldyn->func3b_j1)
1689                                         moldyn->func3b_j1(moldyn,
1690                                                           &(itom[i]),
1691                                                           jtom,
1692                                                           bc_ij);
1693
1694                                 /* in first j loop, 3bp run can be skipped */
1695                                 if(!(moldyn->run3bp))
1696                                         continue;
1697                         
1698                                 /* first loop over atoms k */
1699                                 if(moldyn->func3b_k1) {
1700
1701                                 for(k=0;k<27;k++) {
1702
1703                                         that=&(neighbour_i2[k]);
1704                                         list_reset_f(that);
1705                                         
1706                                         if(that->start==NULL)
1707                                                 continue;
1708
1709                                         bc_ik=(k<dnlc)?0:1;
1710
1711                                         do {
1712
1713                                                 ktom=that->current->data;
1714
1715                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
1716                                                         continue;
1717
1718                                                 if(ktom==jtom)
1719                                                         continue;
1720
1721                                                 if(ktom==&(itom[i]))
1722                                                         continue;
1723
1724                                                 moldyn->func3b_k1(moldyn,
1725                                                                   &(itom[i]),
1726                                                                   jtom,
1727                                                                   ktom,
1728                                                                   bc_ik|bc_ij);
1729
1730                                         } while(list_next_f(that)!=\
1731                                                 L_NO_NEXT_ELEMENT);
1732
1733                                 }
1734
1735                                 }
1736
1737                                 if(moldyn->func3b_j2)
1738                                         moldyn->func3b_j2(moldyn,
1739                                                           &(itom[i]),
1740                                                           jtom,
1741                                                           bc_ij);
1742
1743                                 /* second loop over atoms k */
1744                                 if(moldyn->func3b_k2) {
1745
1746                                 for(k=0;k<27;k++) {
1747
1748                                         that=&(neighbour_i2[k]);
1749                                         list_reset_f(that);
1750                                         
1751                                         if(that->start==NULL)
1752                                                 continue;
1753
1754                                         bc_ik=(k<dnlc)?0:1;
1755
1756                                         do {
1757
1758                                                 ktom=that->current->data;
1759
1760                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
1761                                                         continue;
1762
1763                                                 if(ktom==jtom)
1764                                                         continue;
1765
1766                                                 if(ktom==&(itom[i]))
1767                                                         continue;
1768
1769                                                 moldyn->func3b_k2(moldyn,
1770                                                                   &(itom[i]),
1771                                                                   jtom,
1772                                                                   ktom,
1773                                                                   bc_ik|bc_ij);
1774
1775                                         } while(list_next_f(that)!=\
1776                                                 L_NO_NEXT_ELEMENT);
1777
1778                                 }
1779                                 
1780                                 }
1781
1782                                 /* 2bp post function */
1783                                 if(moldyn->func3b_j3) {
1784                                         moldyn->func3b_j3(moldyn,
1785                                                           &(itom[i]),
1786                                                           jtom,bc_ij);
1787                                 }
1788                                         
1789                         } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
1790                 
1791                 }
1792                 
1793 #ifdef DEBUG
1794         //printf("\n\n");
1795 #endif
1796 #ifdef VDEBUG
1797         printf("\n\n");
1798 #endif
1799
1800         }
1801
1802 #ifdef DEBUG
1803         printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
1804 #endif
1805
1806         /* calculate global virial */
1807         for(i=0;i<count;i++) {
1808                 moldyn->gvir.xx+=moldyn->atom[i].r.x*moldyn->atom[i].f.x;
1809                 moldyn->gvir.yy+=moldyn->atom[i].r.y*moldyn->atom[i].f.y;
1810                 moldyn->gvir.zz+=moldyn->atom[i].r.z*moldyn->atom[i].f.z;
1811                 moldyn->gvir.xy+=moldyn->atom[i].r.y*moldyn->atom[i].f.x;
1812                 moldyn->gvir.xz+=moldyn->atom[i].r.z*moldyn->atom[i].f.x;
1813                 moldyn->gvir.yz+=moldyn->atom[i].r.z*moldyn->atom[i].f.y;
1814         }
1815
1816         return 0;
1817 }
1818
1819 /*
1820  * virial calculation
1821  */
1822
1823 //inline int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) {
1824 int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) {
1825
1826         a->virial.xx+=f->x*d->x;
1827         a->virial.yy+=f->y*d->y;
1828         a->virial.zz+=f->z*d->z;
1829         a->virial.xy+=f->x*d->y;
1830         a->virial.xz+=f->x*d->z;
1831         a->virial.yz+=f->y*d->z;
1832
1833         return 0;
1834 }
1835
1836 /*
1837  * periodic boundary checking
1838  */
1839
1840 //inline int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
1841 int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
1842         
1843         double x,y,z;
1844         t_3dvec *dim;
1845
1846         dim=&(moldyn->dim);
1847
1848         x=dim->x/2;
1849         y=dim->y/2;
1850         z=dim->z/2;
1851
1852         if(moldyn->status&MOLDYN_STAT_PBX) {
1853                 if(a->x>=x) a->x-=dim->x;
1854                 else if(-a->x>x) a->x+=dim->x;
1855         }
1856         if(moldyn->status&MOLDYN_STAT_PBY) {
1857                 if(a->y>=y) a->y-=dim->y;
1858                 else if(-a->y>y) a->y+=dim->y;
1859         }
1860         if(moldyn->status&MOLDYN_STAT_PBZ) {
1861                 if(a->z>=z) a->z-=dim->z;
1862                 else if(-a->z>z) a->z+=dim->z;
1863         }
1864
1865         return 0;
1866 }
1867         
1868 /*
1869  * debugging / critical check functions
1870  */
1871
1872 int moldyn_bc_check(t_moldyn *moldyn) {
1873
1874         t_atom *atom;
1875         t_3dvec *dim;
1876         int i;
1877         double x;
1878         u8 byte;
1879         int j,k;
1880
1881         atom=moldyn->atom;
1882         dim=&(moldyn->dim);
1883         x=dim->x/2;
1884
1885         for(i=0;i<moldyn->count;i++) {
1886                 if(atom[i].r.x>=dim->x/2||-atom[i].r.x>dim->x/2) {
1887                         printf("FATAL: atom %d: x: %.20f (%.20f)\n",
1888                                i,atom[i].r.x,dim->x/2);
1889                         printf("diagnostic:\n");
1890                         printf("-----------\natom.r.x:\n");
1891                         for(j=0;j<8;j++) {
1892                                 memcpy(&byte,(u8 *)(&(atom[i].r.x))+j,1);
1893                                 for(k=0;k<8;k++)
1894                                         printf("%d%c",
1895                                         ((byte)&(1<<k))?1:0,
1896                                         (k==7)?'\n':'|');
1897                         }
1898                         printf("---------------\nx=dim.x/2:\n");
1899                         for(j=0;j<8;j++) {
1900                                 memcpy(&byte,(u8 *)(&x)+j,1);
1901                                 for(k=0;k<8;k++)
1902                                         printf("%d%c",
1903                                         ((byte)&(1<<k))?1:0,
1904                                         (k==7)?'\n':'|');
1905                         }
1906                         if(atom[i].r.x==x) printf("the same!\n");
1907                         else printf("different!\n");
1908                 }
1909                 if(atom[i].r.y>=dim->y/2||-atom[i].r.y>dim->y/2)
1910                         printf("FATAL: atom %d: y: %.20f (%.20f)\n",
1911                                i,atom[i].r.y,dim->y/2);
1912                 if(atom[i].r.z>=dim->z/2||-atom[i].r.z>dim->z/2)
1913                         printf("FATAL: atom %d: z: %.20f (%.20f)\n",
1914                                i,atom[i].r.z,dim->z/2);
1915         }
1916
1917         return 0;
1918 }
1919
1920 /*
1921  * post processing functions
1922  */
1923
1924 int get_line(int fd,char *line,int max) {
1925
1926         int count,ret;
1927
1928         count=0;
1929
1930         while(1) {
1931                 if(count==max) return count;
1932                 ret=read(fd,line+count,1);
1933                 if(ret<=0) return ret;
1934                 if(line[count]=='\n') {
1935                         line[count]='\0';
1936                         return count+1;
1937                 }
1938                 count+=1;
1939         }
1940 }
1941