]> www.hackdaworld.org Git - physik/posic.git/blob - moldyn.c
9bd216d7fc66d8e171459e60d32159a7eb5340ba
[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 <sys/time.h>
17 #include <time.h>
18 #include <math.h>
19
20 #include "moldyn.h"
21 #include "report/report.h"
22
23 /* potential includes */
24 #include "potentials/harmonic_oscillator.h"
25 #include "potentials/lennard_jones.h"
26 #include "potentials/albe.h"
27 #ifdef TERSOFF_ORIG
28 #include "potentials/tersoff_orig.h"
29 #else
30 #include "potentials/tersoff.h"
31 #endif
32
33 /*
34  * the moldyn functions
35  */
36
37 int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
38
39         printf("[moldyn] init\n");
40
41         memset(moldyn,0,sizeof(t_moldyn));
42
43         moldyn->argc=argc;
44         moldyn->args=argv;
45
46         rand_init(&(moldyn->random),NULL,1);
47         moldyn->random.status|=RAND_STAT_VERBOSE;
48
49         return 0;
50 }
51
52 int moldyn_shutdown(t_moldyn *moldyn) {
53
54         printf("[moldyn] shutdown\n");
55
56         moldyn_log_shutdown(moldyn);
57         link_cell_shutdown(moldyn);
58         rand_close(&(moldyn->random));
59         free(moldyn->atom);
60
61         return 0;
62 }
63
64 int set_int_alg(t_moldyn *moldyn,u8 algo) {
65
66         printf("[moldyn] integration algorithm: ");
67
68         switch(algo) {
69                 case MOLDYN_INTEGRATE_VERLET:
70                         moldyn->integrate=velocity_verlet;
71                         printf("velocity verlet\n");
72                         break;
73                 default:
74                         printf("unknown integration algorithm: %02x\n",algo);
75                         printf("unknown\n");
76                         return -1;
77         }
78
79         return 0;
80 }
81
82 int set_cutoff(t_moldyn *moldyn,double cutoff) {
83
84         moldyn->cutoff=cutoff;
85         moldyn->cutoff_square=cutoff*cutoff;
86
87         printf("[moldyn] cutoff [A]: %f\n",moldyn->cutoff);
88
89         return 0;
90 }
91
92 int set_temperature(t_moldyn *moldyn,double t_ref) {
93
94         moldyn->t_ref=t_ref;
95
96         printf("[moldyn] temperature [K]: %f\n",moldyn->t_ref);
97
98         return 0;
99 }
100
101 int set_pressure(t_moldyn *moldyn,double p_ref) {
102
103         moldyn->p_ref=p_ref;
104
105         printf("[moldyn] pressure [bar]: %f\n",moldyn->p_ref/BAR);
106
107         return 0;
108 }
109
110 int set_p_scale(t_moldyn *moldyn,u8 ptype,double ptc) {
111
112         moldyn->pt_scale&=(~(P_SCALE_MASK));
113         moldyn->pt_scale|=ptype;
114         moldyn->p_tc=ptc;
115
116         printf("[moldyn] p scaling:\n");
117
118         printf("  p: %s",ptype?"yes":"no ");
119         if(ptype)
120                 printf(" | type: %02x | factor: %f",ptype,ptc);
121         printf("\n");
122
123         return 0;
124 }
125
126 int set_t_scale(t_moldyn *moldyn,u8 ttype,double ttc) {
127
128         moldyn->pt_scale&=(~(T_SCALE_MASK));
129         moldyn->pt_scale|=ttype;
130         moldyn->t_tc=ttc;
131
132         printf("[moldyn] t scaling:\n");
133
134         printf("  t: %s",ttype?"yes":"no ");
135         if(ttype)
136                 printf(" | type: %02x | factor: %f",ttype,ttc);
137         printf("\n");
138
139         return 0;
140 }
141
142 int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc) {
143
144         moldyn->pt_scale=(ptype|ttype);
145         moldyn->t_tc=ttc;
146         moldyn->p_tc=ptc;
147
148         printf("[moldyn] p/t scaling:\n");
149
150         printf("  p: %s",ptype?"yes":"no ");
151         if(ptype)
152                 printf(" | type: %02x | factor: %f",ptype,ptc);
153         printf("\n");
154
155         printf("  t: %s",ttype?"yes":"no ");
156         if(ttype)
157                 printf(" | type: %02x | factor: %f",ttype,ttc);
158         printf("\n");
159
160         return 0;
161 }
162
163 int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize) {
164
165         moldyn->dim.x=x;
166         moldyn->dim.y=y;
167         moldyn->dim.z=z;
168
169         moldyn->volume=x*y*z;
170
171         if(visualize) {
172                 moldyn->vis.dim.x=x;
173                 moldyn->vis.dim.y=y;
174                 moldyn->vis.dim.z=z;
175         }
176
177         printf("[moldyn] dimensions in A and A^3 respectively:\n");
178         printf("  x: %f\n",moldyn->dim.x);
179         printf("  y: %f\n",moldyn->dim.y);
180         printf("  z: %f\n",moldyn->dim.z);
181         printf("  volume: %f\n",moldyn->volume);
182         printf("  visualize simulation box: %s\n",visualize?"yes":"no");
183
184         return 0;
185 }
186
187 int set_nn_dist(t_moldyn *moldyn,double dist) {
188
189         moldyn->nnd=dist;
190
191         return 0;
192 }
193
194 int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z) {
195
196         printf("[moldyn] periodic boundary conditions:\n");
197
198         if(x)
199                 moldyn->status|=MOLDYN_STAT_PBX;
200
201         if(y)
202                 moldyn->status|=MOLDYN_STAT_PBY;
203
204         if(z)
205                 moldyn->status|=MOLDYN_STAT_PBZ;
206
207         printf("  x: %s\n",x?"yes":"no");
208         printf("  y: %s\n",y?"yes":"no");
209         printf("  z: %s\n",z?"yes":"no");
210
211         return 0;
212 }
213
214 int set_potential(t_moldyn *moldyn,u8 type) {
215
216         switch(type) {
217                 case MOLDYN_POTENTIAL_TM:
218                         moldyn->func1b=tersoff_mult_1bp;
219                         moldyn->func3b_j1=tersoff_mult_3bp_j1;
220                         moldyn->func3b_k1=tersoff_mult_3bp_k1;
221                         moldyn->func3b_j2=tersoff_mult_3bp_j2;
222                         moldyn->func3b_k2=tersoff_mult_3bp_k2;
223                         moldyn->check_2b_bond=tersoff_mult_check_2b_bond;
224                         break;
225                 case MOLDYN_POTENTIAL_AM:
226                         moldyn->func3b_j1=albe_mult_3bp_j1;
227                         moldyn->func3b_k1=albe_mult_3bp_k1;
228                         moldyn->func3b_j2=albe_mult_3bp_j2;
229                         moldyn->func3b_k2=albe_mult_3bp_k2;
230                         moldyn->check_2b_bond=albe_mult_check_2b_bond;
231                         break;
232                 case MOLDYN_POTENTIAL_HO:
233                         moldyn->func2b=harmonic_oscillator;
234                         moldyn->check_2b_bond=harmonic_oscillator_check_2b_bond;
235                         break;
236                 case MOLDYN_POTENTIAL_LJ:
237                         moldyn->func2b=lennard_jones;
238                         moldyn->check_2b_bond=lennard_jones_check_2b_bond;
239                         break;
240                 default:
241                         printf("[moldyn] set potential: unknown type %02x\n",
242                                type);
243                         return -1;
244         }
245
246         return 0;
247 }
248
249 int set_avg_skip(t_moldyn *moldyn,int skip) {
250
251         printf("[moldyn] skip %d steps before starting average calc\n",skip);
252         moldyn->avg_skip=skip;
253
254         return 0;
255 }
256
257 int moldyn_set_log_dir(t_moldyn *moldyn,char *dir) {
258
259         strncpy(moldyn->vlsdir,dir,127);
260
261         return 0;
262 }
263
264 int moldyn_set_report(t_moldyn *moldyn,char *author,char *title) {
265
266         strncpy(moldyn->rauthor,author,63);
267         strncpy(moldyn->rtitle,title,63);
268
269         return 0;
270 }
271         
272 int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) {
273
274         char filename[128];
275         int ret;
276
277         printf("[moldyn] set log: ");
278
279         switch(type) {
280                 case LOG_TOTAL_ENERGY:
281                         moldyn->ewrite=timer;
282                         snprintf(filename,127,"%s/energy",moldyn->vlsdir);
283                         moldyn->efd=open(filename,
284                                          O_WRONLY|O_CREAT|O_EXCL,
285                                          S_IRUSR|S_IWUSR);
286                         if(moldyn->efd<0) {
287                                 perror("[moldyn] energy log fd open");
288                                 return moldyn->efd;
289                         }
290                         dprintf(moldyn->efd,"# total energy log file\n");
291                         printf("total energy (%d)\n",timer);
292                         break;
293                 case LOG_TOTAL_MOMENTUM:
294                         moldyn->mwrite=timer;
295                         snprintf(filename,127,"%s/momentum",moldyn->vlsdir);
296                         moldyn->mfd=open(filename,
297                                          O_WRONLY|O_CREAT|O_EXCL,
298                                          S_IRUSR|S_IWUSR);
299                         if(moldyn->mfd<0) {
300                                 perror("[moldyn] momentum log fd open");
301                                 return moldyn->mfd;
302                         }
303                         dprintf(moldyn->efd,"# total momentum log file\n");
304                         printf("total momentum (%d)\n",timer);
305                         break;
306                 case LOG_PRESSURE:
307                         moldyn->pwrite=timer;
308                         snprintf(filename,127,"%s/pressure",moldyn->vlsdir);
309                         moldyn->pfd=open(filename,
310                                          O_WRONLY|O_CREAT|O_EXCL,
311                                          S_IRUSR|S_IWUSR);
312                         if(moldyn->pfd<0) {
313                                 perror("[moldyn] pressure log file\n");
314                                 return moldyn->pfd;
315                         }
316                         dprintf(moldyn->pfd,"# pressure log file\n");
317                         printf("pressure (%d)\n",timer);
318                         break;
319                 case LOG_TEMPERATURE:
320                         moldyn->twrite=timer;
321                         snprintf(filename,127,"%s/temperature",moldyn->vlsdir);
322                         moldyn->tfd=open(filename,
323                                          O_WRONLY|O_CREAT|O_EXCL,
324                                          S_IRUSR|S_IWUSR);
325                         if(moldyn->tfd<0) {
326                                 perror("[moldyn] temperature log file\n");
327                                 return moldyn->tfd;
328                         }
329                         dprintf(moldyn->tfd,"# temperature log file\n");
330                         printf("temperature (%d)\n",timer);
331                         break;
332                 case LOG_VOLUME:
333                         moldyn->vwrite=timer;
334                         snprintf(filename,127,"%s/volume",moldyn->vlsdir);
335                         moldyn->vfd=open(filename,
336                                          O_WRONLY|O_CREAT|O_EXCL,
337                                          S_IRUSR|S_IWUSR);
338                         if(moldyn->vfd<0) {
339                                 perror("[moldyn] volume log file\n");
340                                 return moldyn->vfd;
341                         }
342                         dprintf(moldyn->vfd,"# volume log file\n");
343                         printf("volume (%d)\n",timer);
344                         break;
345                 case SAVE_STEP:
346                         moldyn->swrite=timer;
347                         printf("save file (%d)\n",timer);
348                         break;
349                 case VISUAL_STEP:
350                         moldyn->awrite=timer;
351                         ret=visual_init(moldyn,moldyn->vlsdir);
352                         if(ret<0) {
353                                 printf("[moldyn] visual init failure\n");
354                                 return ret;
355                         }
356                         printf("visual file (%d)\n",timer);
357                         break;
358                 case CREATE_REPORT:
359                         snprintf(filename,127,"%s/report.tex",moldyn->vlsdir);
360                         moldyn->rfd=open(filename,
361                                          O_WRONLY|O_CREAT|O_EXCL,
362                                          S_IRUSR|S_IWUSR);
363                         if(moldyn->rfd<0) {
364                                 perror("[moldyn] report fd open");      
365                                 return moldyn->rfd;
366                         }
367                         printf("report -> ");
368                         if(moldyn->efd) {
369                                 snprintf(filename,127,"%s/e_plot.scr",
370                                          moldyn->vlsdir);
371                                 moldyn->epfd=open(filename,
372                                                  O_WRONLY|O_CREAT|O_EXCL,
373                                                  S_IRUSR|S_IWUSR);
374                                 if(moldyn->epfd<0) {
375                                         perror("[moldyn] energy plot fd open");
376                                         return moldyn->epfd;
377                                 }
378                                 dprintf(moldyn->epfd,e_plot_script);
379                                 close(moldyn->epfd);
380                                 printf("energy ");
381                         }
382                         if(moldyn->pfd) {
383                                 snprintf(filename,127,"%s/pressure_plot.scr",
384                                          moldyn->vlsdir);
385                                 moldyn->ppfd=open(filename,
386                                                   O_WRONLY|O_CREAT|O_EXCL,
387                                                   S_IRUSR|S_IWUSR);
388                                 if(moldyn->ppfd<0) {
389                                         perror("[moldyn] p plot fd open");
390                                         return moldyn->ppfd;
391                                 }
392                                 dprintf(moldyn->ppfd,pressure_plot_script);
393                                 close(moldyn->ppfd);
394                                 printf("pressure ");
395                         }
396                         if(moldyn->tfd) {
397                                 snprintf(filename,127,"%s/temperature_plot.scr",
398                                          moldyn->vlsdir);
399                                 moldyn->tpfd=open(filename,
400                                                   O_WRONLY|O_CREAT|O_EXCL,
401                                                   S_IRUSR|S_IWUSR);
402                                 if(moldyn->tpfd<0) {
403                                         perror("[moldyn] t plot fd open");
404                                         return moldyn->tpfd;
405                                 }
406                                 dprintf(moldyn->tpfd,temperature_plot_script);
407                                 close(moldyn->tpfd);
408                                 printf("temperature ");
409                         }
410                         dprintf(moldyn->rfd,report_start,
411                                 moldyn->rauthor,moldyn->rtitle);
412                         printf("\n");
413                         break;
414                 default:
415                         printf("unknown log type: %02x\n",type);
416                         return -1;
417         }
418
419         return 0;
420 }
421
422 int moldyn_log_shutdown(t_moldyn *moldyn) {
423
424         char sc[256];
425
426         printf("[moldyn] log shutdown\n");
427         if(moldyn->efd) {
428                 close(moldyn->efd);
429                 if(moldyn->rfd) {
430                         dprintf(moldyn->rfd,report_energy);
431                         snprintf(sc,255,"cd %s && gnuplot e_plot.scr",
432                                  moldyn->vlsdir);
433                         system(sc);
434                 }
435         }
436         if(moldyn->mfd) close(moldyn->mfd);
437         if(moldyn->pfd) {
438                 close(moldyn->pfd);
439                 if(moldyn->rfd)
440                         dprintf(moldyn->rfd,report_pressure);
441                         snprintf(sc,255,"cd %s && gnuplot pressure_plot.scr",
442                                  moldyn->vlsdir);
443                         system(sc);
444         }
445         if(moldyn->tfd) {
446                 close(moldyn->tfd);
447                 if(moldyn->rfd)
448                         dprintf(moldyn->rfd,report_temperature);
449                         snprintf(sc,255,"cd %s && gnuplot temperature_plot.scr",
450                                  moldyn->vlsdir);
451                         system(sc);
452         }
453         if(moldyn->rfd) {
454                 dprintf(moldyn->rfd,report_end);
455                 close(moldyn->rfd);
456                 snprintf(sc,255,"cd %s && pdflatex report >/dev/null 2>&1",
457                          moldyn->vlsdir);
458                 system(sc);
459                 snprintf(sc,255,"cd %s && pdflatex report >/dev/null 2>&1",
460                          moldyn->vlsdir);
461                 system(sc);
462                 snprintf(sc,255,"cd %s && dvipdf report >/dev/null 2>&1",
463                          moldyn->vlsdir);
464                 system(sc);
465         }
466
467         return 0;
468 }
469
470 /*
471  * creating lattice functions
472  */
473
474 int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
475                    u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin) {
476
477         int new,count;
478         int ret;
479         t_3dvec orig;
480         void *ptr;
481         t_atom *atom;
482         char name[16];
483
484         new=a*b*c;
485         count=moldyn->count;
486
487         /* how many atoms do we expect */
488         if(type==NONE) {
489                 new*=1;
490                 printf("[moldyn] WARNING: create 'none' lattice called");
491         }
492         if(type==CUBIC) new*=1;
493         if(type==FCC) new*=4;
494         if(type==DIAMOND) new*=8;
495
496         /* allocate space for atoms */
497         ptr=realloc(moldyn->atom,(count+new)*sizeof(t_atom));
498         if(!ptr) {
499                 perror("[moldyn] realloc (create lattice)");
500                 return -1;
501         }
502         moldyn->atom=ptr;
503         atom=&(moldyn->atom[count]);
504
505         /* no atoms on the boundaries (only reason: it looks better!) */
506         if(!origin) {
507                 orig.x=0.5*lc;
508                 orig.y=0.5*lc;
509                 orig.z=0.5*lc;
510         }
511         else {
512                 orig.x=origin->x;
513                 orig.y=origin->y;
514                 orig.z=origin->z;
515         }
516
517         switch(type) {
518                 case CUBIC:
519                         set_nn_dist(moldyn,lc);
520                         ret=cubic_init(a,b,c,lc,atom,&orig);
521                         strcpy(name,"cubic");
522                         break;
523                 case FCC:
524                         if(!origin)
525                                 v3_scale(&orig,&orig,0.5);
526                         set_nn_dist(moldyn,0.5*sqrt(2.0)*lc);
527                         ret=fcc_init(a,b,c,lc,atom,&orig);
528                         strcpy(name,"fcc");
529                         break;
530                 case DIAMOND:
531                         if(!origin)
532                                 v3_scale(&orig,&orig,0.25);
533                         set_nn_dist(moldyn,0.25*sqrt(3.0)*lc);
534                         ret=diamond_init(a,b,c,lc,atom,&orig);
535                         strcpy(name,"diamond");
536                         break;
537                 default:
538                         printf("unknown lattice type (%02x)\n",type);
539                         return -1;
540         }
541
542         /* debug */
543         if(ret!=new) {
544                 printf("[moldyn] creating lattice failed\n");
545                 printf("  amount of atoms\n");
546                 printf("  - expected: %d\n",new);
547                 printf("  - created: %d\n",ret);
548                 return -1;
549         }
550
551         moldyn->count+=new;
552         printf("[moldyn] created %s lattice with %d atoms\n",name,new);
553
554         for(ret=0;ret<new;ret++) {
555                 atom[ret].element=element;
556                 atom[ret].mass=mass;
557                 atom[ret].attr=attr;
558                 atom[ret].brand=brand;
559                 atom[ret].tag=count+ret;
560                 check_per_bound(moldyn,&(atom[ret].r));
561                 atom[ret].r_0=atom[ret].r;
562         }
563
564         /* update total system mass */
565         total_mass_calc(moldyn);
566
567         return ret;
568 }
569
570 int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
571              t_3dvec *r,t_3dvec *v) {
572
573         t_atom *atom;
574         void *ptr;
575         int count;
576         
577         atom=moldyn->atom;
578         count=(moldyn->count)++;        // asshole style!
579
580         ptr=realloc(atom,(count+1)*sizeof(t_atom));
581         if(!ptr) {
582                 perror("[moldyn] realloc (add atom)");
583                 return -1;
584         }
585         moldyn->atom=ptr;
586
587         atom=moldyn->atom;
588
589         /* initialize new atom */
590         memset(&(atom[count]),0,sizeof(t_atom));
591         atom[count].r=*r;
592         atom[count].v=*v;
593         atom[count].element=element;
594         atom[count].mass=mass;
595         atom[count].brand=brand;
596         atom[count].tag=count;
597         atom[count].attr=attr;
598         check_per_bound(moldyn,&(atom[count].r));
599         atom[count].r_0=atom[count].r;
600
601         /* update total system mass */
602         total_mass_calc(moldyn);
603
604         return 0;
605 }
606
607 int del_atom(t_moldyn *moldyn,int tag) {
608
609         t_atom *new,*old;
610         int cnt;
611
612         old=moldyn->atom;
613
614         new=(t_atom *)malloc((moldyn->count-1)*sizeof(t_atom));
615         if(!new) {
616                 perror("[moldyn]malloc (del atom)");
617                 return -1;
618         }
619
620         for(cnt=0;cnt<tag;cnt++)
621                 new[cnt]=old[cnt];
622         
623         for(cnt=tag+1;cnt<moldyn->count;cnt++) {
624                 new[cnt-1]=old[cnt];
625                 new[cnt-1].tag=cnt-1;
626         }
627
628         moldyn->count-=1;
629         moldyn->atom=new;
630
631         free(old);
632
633         return 0;
634 }
635
636 /* cubic init */
637 int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
638
639         int count;
640         t_3dvec r;
641         int i,j,k;
642         t_3dvec o;
643
644         count=0;
645         if(origin)
646                 v3_copy(&o,origin);
647         else
648                 v3_zero(&o);
649
650         r.x=o.x;
651         for(i=0;i<a;i++) {
652                 r.y=o.y;
653                 for(j=0;j<b;j++) {
654                         r.z=o.z;
655                         for(k=0;k<c;k++) {
656                                 v3_copy(&(atom[count].r),&r);
657                                 count+=1;
658                                 r.z+=lc;
659                         }
660                         r.y+=lc;
661                 }
662                 r.x+=lc;
663         }
664
665         for(i=0;i<count;i++) {
666                 atom[i].r.x-=(a*lc)/2.0;
667                 atom[i].r.y-=(b*lc)/2.0;
668                 atom[i].r.z-=(c*lc)/2.0;
669         }
670
671         return count;
672 }
673
674 /* fcc lattice init */
675 int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
676
677         int count;
678         int i,j,k,l;
679         t_3dvec o,r,n;
680         t_3dvec basis[3];
681
682         count=0;
683         if(origin)
684                 v3_copy(&o,origin);
685         else
686                 v3_zero(&o);
687
688         /* construct the basis */
689         memset(basis,0,3*sizeof(t_3dvec));
690         basis[0].x=0.5*lc;
691         basis[0].y=0.5*lc;
692         basis[1].x=0.5*lc;
693         basis[1].z=0.5*lc;
694         basis[2].y=0.5*lc;
695         basis[2].z=0.5*lc;
696
697         /* fill up the room */
698         r.x=o.x;
699         for(i=0;i<a;i++) {
700                 r.y=o.y;
701                 for(j=0;j<b;j++) {
702                         r.z=o.z;
703                         for(k=0;k<c;k++) {
704                                 /* first atom */
705                                 v3_copy(&(atom[count].r),&r);
706                                 count+=1;
707                                 r.z+=lc;
708                                 /* the three face centered atoms */
709                                 for(l=0;l<3;l++) {
710                                         v3_add(&n,&r,&basis[l]);
711                                         v3_copy(&(atom[count].r),&n);
712                                         count+=1;
713                                 }
714                         }
715                         r.y+=lc;
716                 }
717                 r.x+=lc;
718         }
719                                 
720         /* coordinate transformation */
721         for(i=0;i<count;i++) {
722                 atom[i].r.x-=(a*lc)/2.0;
723                 atom[i].r.y-=(b*lc)/2.0;
724                 atom[i].r.z-=(c*lc)/2.0;
725         }
726
727         return count;
728 }
729
730 int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
731
732         int count;
733         t_3dvec o;
734
735         count=fcc_init(a,b,c,lc,atom,origin);
736
737         o.x=0.25*lc;
738         o.y=0.25*lc;
739         o.z=0.25*lc;
740
741         if(origin) v3_add(&o,&o,origin);
742
743         count+=fcc_init(a,b,c,lc,&atom[count],&o);
744
745         return count;
746 }
747
748 int destroy_atoms(t_moldyn *moldyn) {
749
750         if(moldyn->atom) free(moldyn->atom);
751
752         return 0;
753 }
754
755 int thermal_init(t_moldyn *moldyn,u8 equi_init) {
756
757         /*
758          * - gaussian distribution of velocities
759          * - zero total momentum
760          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
761          */
762
763         int i;
764         double v,sigma;
765         t_3dvec p_total,delta;
766         t_atom *atom;
767         t_random *random;
768
769         atom=moldyn->atom;
770         random=&(moldyn->random);
771
772         printf("[moldyn] thermal init (equi init: %s)\n",equi_init?"yes":"no");
773
774         /* gaussian distribution of velocities */
775         v3_zero(&p_total);
776         for(i=0;i<moldyn->count;i++) {
777                 sigma=sqrt(2.0*K_BOLTZMANN*moldyn->t_ref/atom[i].mass);
778                 /* x direction */
779                 v=sigma*rand_get_gauss(random);
780                 atom[i].v.x=v;
781                 p_total.x+=atom[i].mass*v;
782                 /* y direction */
783                 v=sigma*rand_get_gauss(random);
784                 atom[i].v.y=v;
785                 p_total.y+=atom[i].mass*v;
786                 /* z direction */
787                 v=sigma*rand_get_gauss(random);
788                 atom[i].v.z=v;
789                 p_total.z+=atom[i].mass*v;
790         }
791
792         /* zero total momentum */
793         v3_scale(&p_total,&p_total,1.0/moldyn->count);
794         for(i=0;i<moldyn->count;i++) {
795                 v3_scale(&delta,&p_total,1.0/atom[i].mass);
796                 v3_sub(&(atom[i].v),&(atom[i].v),&delta);
797         }
798
799         /* velocity scaling */
800         scale_velocity(moldyn,equi_init);
801
802         return 0;
803 }
804
805 double total_mass_calc(t_moldyn *moldyn) {
806
807         int i;
808
809         moldyn->mass=0.0;
810
811         for(i=0;i<moldyn->count;i++)
812                 moldyn->mass+=moldyn->atom[i].mass;
813
814         return moldyn->mass;
815 }
816
817 double temperature_calc(t_moldyn *moldyn) {
818
819         /* assume up to date kinetic energy, which is 3/2 N k_B T */
820
821         moldyn->t=(2.0*moldyn->ekin)/(3.0*K_BOLTZMANN*moldyn->count);
822
823         return moldyn->t;
824 }
825
826 double get_temperature(t_moldyn *moldyn) {
827
828         return moldyn->t;
829 }
830
831 int scale_velocity(t_moldyn *moldyn,u8 equi_init) {
832
833         int i;
834         double e,scale;
835         t_atom *atom;
836         int count;
837
838         atom=moldyn->atom;
839
840         /*
841          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
842          */
843
844         /* get kinetic energy / temperature & count involved atoms */
845         e=0.0;
846         count=0;
847         for(i=0;i<moldyn->count;i++) {
848                 if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB)) {
849                         e+=atom[i].mass*v3_absolute_square(&(atom[i].v));
850                         count+=1;
851                 }
852         }
853         e*=0.5;
854         if(count!=0) moldyn->t=e/(1.5*count*K_BOLTZMANN);
855         else return 0;  /* no atoms involved in scaling! */
856         
857         /* (temporary) hack for e,t = 0 */
858         if(e==0.0) {
859         moldyn->t=0.0;
860                 if(moldyn->t_ref!=0.0) {
861                         thermal_init(moldyn,equi_init);
862                         return 0;
863                 }
864                 else
865                         return 0; /* no scaling needed */
866         }
867
868
869         /* get scaling factor */
870         scale=moldyn->t_ref/moldyn->t;
871         if(equi_init&TRUE)
872                 scale*=2.0;
873         else
874                 if(moldyn->pt_scale&T_SCALE_BERENDSEN)
875                         scale=1.0+(scale-1.0)/moldyn->t_tc;
876         scale=sqrt(scale);
877
878         /* velocity scaling */
879         for(i=0;i<moldyn->count;i++) {
880                 if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB))
881                         v3_scale(&(atom[i].v),&(atom[i].v),scale);
882         }
883
884         return 0;
885 }
886
887 double ideal_gas_law_pressure(t_moldyn *moldyn) {
888
889         double p;
890
891         p=moldyn->count*moldyn->t*K_BOLTZMANN/moldyn->volume;
892
893         return p;
894 }
895
896 double virial_sum(t_moldyn *moldyn) {
897
898         int i;
899         t_virial *virial;
900
901         /* virial (sum over atom virials) */
902         moldyn->virial=0.0;
903         moldyn->vir.xx=0.0;
904         moldyn->vir.yy=0.0;
905         moldyn->vir.zz=0.0;
906         moldyn->vir.xy=0.0;
907         moldyn->vir.xz=0.0;
908         moldyn->vir.yz=0.0;
909         for(i=0;i<moldyn->count;i++) {
910                 virial=&(moldyn->atom[i].virial);
911                 moldyn->virial+=(virial->xx+virial->yy+virial->zz);
912                 moldyn->vir.xx+=virial->xx;
913                 moldyn->vir.yy+=virial->yy;
914                 moldyn->vir.zz+=virial->zz;
915                 moldyn->vir.xy+=virial->xy;
916                 moldyn->vir.xz+=virial->xz;
917                 moldyn->vir.yz+=virial->yz;
918         }
919
920         /* global virial (absolute coordinates) */
921         virial=&(moldyn->gvir);
922         moldyn->gv=virial->xx+virial->yy+virial->zz;
923
924         return moldyn->virial;
925 }
926
927 double pressure_calc(t_moldyn *moldyn) {
928
929         /*
930          * PV = NkT + <W>
931          * with W = 1/3 sum_i f_i r_i (- skipped!)
932          * virial = sum_i f_i r_i
933          * 
934          * => P = (2 Ekin + virial) / (3V)
935          */
936
937         /* assume up to date virial & up to date kinetic energy */
938
939         /* pressure (atom virials) */
940         moldyn->p=2.0*moldyn->ekin+moldyn->virial;
941         moldyn->p/=(3.0*moldyn->volume);
942
943         /* pressure (absolute coordinates) */
944         moldyn->gp=2.0*moldyn->ekin+moldyn->gv;
945         moldyn->gp/=(3.0*moldyn->volume);
946
947         return moldyn->p;
948 }
949
950 int average_reset(t_moldyn *moldyn) {
951
952         printf("[moldyn] average reset\n");
953
954         /* update skip value */
955         moldyn->avg_skip=moldyn->total_steps;
956
957         /* kinetic energy */
958         moldyn->k_sum=0.0;
959         moldyn->k2_sum=0.0;
960         
961         /* potential energy */
962         moldyn->v_sum=0.0;
963         moldyn->v2_sum=0.0;
964
965         /* temperature */
966         moldyn->t_sum=0.0;
967
968         /* virial */
969         moldyn->virial_sum=0.0;
970         moldyn->gv_sum=0.0;
971
972         /* pressure */
973         moldyn->p_sum=0.0;
974         moldyn->gp_sum=0.0;
975         moldyn->tp_sum=0.0;
976
977         return 0;
978 }
979
980 int average_and_fluctuation_calc(t_moldyn *moldyn) {
981
982         int denom;
983
984         if(moldyn->total_steps<moldyn->avg_skip)
985                 return 0;
986
987         denom=moldyn->total_steps+1-moldyn->avg_skip;
988
989         /* assume up to date energies, temperature, pressure etc */
990
991         /* kinetic energy */
992         moldyn->k_sum+=moldyn->ekin;
993         moldyn->k2_sum+=(moldyn->ekin*moldyn->ekin);
994         moldyn->k_avg=moldyn->k_sum/denom;
995         moldyn->k2_avg=moldyn->k2_sum/denom;
996         moldyn->dk2_avg=moldyn->k2_avg-(moldyn->k_avg*moldyn->k_avg);
997
998         /* potential energy */
999         moldyn->v_sum+=moldyn->energy;
1000         moldyn->v2_sum+=(moldyn->energy*moldyn->energy);
1001         moldyn->v_avg=moldyn->v_sum/denom;
1002         moldyn->v2_avg=moldyn->v2_sum/denom;
1003         moldyn->dv2_avg=moldyn->v2_avg-(moldyn->v_avg*moldyn->v_avg);
1004
1005         /* temperature */
1006         moldyn->t_sum+=moldyn->t;
1007         moldyn->t_avg=moldyn->t_sum/denom;
1008
1009         /* virial */
1010         moldyn->virial_sum+=moldyn->virial;
1011         moldyn->virial_avg=moldyn->virial_sum/denom;
1012         moldyn->gv_sum+=moldyn->gv;
1013         moldyn->gv_avg=moldyn->gv_sum/denom;
1014
1015         /* pressure */
1016         moldyn->p_sum+=moldyn->p;
1017         moldyn->p_avg=moldyn->p_sum/denom;
1018         moldyn->gp_sum+=moldyn->gp;
1019         moldyn->gp_avg=moldyn->gp_sum/denom;
1020         moldyn->tp_sum+=moldyn->tp;
1021         moldyn->tp_avg=moldyn->tp_sum/denom;
1022
1023         return 0;
1024 }
1025
1026 int get_heat_capacity(t_moldyn *moldyn) {
1027
1028         double temp2,ighc;
1029
1030         /* averages needed for heat capacity calc */
1031         if(moldyn->total_steps<moldyn->avg_skip)
1032                 return 0;
1033
1034         /* (temperature average)^2 */
1035         temp2=moldyn->t_avg*moldyn->t_avg;
1036         printf("[moldyn] specific heat capacity for T=%f K [J/(kg K)]\n",
1037                moldyn->t_avg);
1038
1039         /* ideal gas contribution */
1040         ighc=3.0*moldyn->count*K_BOLTZMANN/2.0;
1041         printf("  ideal gas contribution: %f\n",
1042                ighc/moldyn->mass*KILOGRAM/JOULE);
1043
1044         /* specific heat for nvt ensemble */
1045         moldyn->c_v_nvt=moldyn->dv2_avg/(K_BOLTZMANN*temp2)+ighc;
1046         moldyn->c_v_nvt/=moldyn->mass;
1047
1048         /* specific heat for nve ensemble */
1049         moldyn->c_v_nve=ighc/(1.0-(moldyn->dv2_avg/(ighc*K_BOLTZMANN*temp2)));
1050         moldyn->c_v_nve/=moldyn->mass;
1051
1052         printf("  NVE: %f\n",moldyn->c_v_nve*KILOGRAM/JOULE);
1053         printf("  NVT: %f\n",moldyn->c_v_nvt*KILOGRAM/JOULE);
1054 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)));
1055
1056         return 0;
1057 }
1058
1059 double thermodynamic_pressure_calc(t_moldyn *moldyn) {
1060
1061         t_3dvec dim;
1062         //t_3dvec *tp;
1063         double h,dv;
1064         double y0,y1;
1065         double su,sd;
1066         t_atom *store;
1067
1068         /*
1069          * dU = - p dV
1070          *
1071          * => p = - dU/dV
1072          *
1073          */
1074
1075         /* store atomic configuration + dimension */
1076         store=malloc(moldyn->count*sizeof(t_atom));
1077         if(store==NULL) {
1078                 printf("[moldyn] allocating store mem failed\n");
1079                 return -1;
1080         }
1081         memcpy(store,moldyn->atom,moldyn->count*sizeof(t_atom));
1082         dim=moldyn->dim;
1083
1084         /* x1, y1 */
1085         sd=0.00001;
1086         h=(1.0-sd)*(1.0-sd)*(1.0-sd);
1087         su=pow(2.0-h,ONE_THIRD)-1.0;
1088         dv=(1.0-h)*moldyn->volume;
1089
1090         /* scale up dimension and atom positions */
1091         scale_dim(moldyn,SCALE_UP,su,TRUE,TRUE,TRUE);
1092         scale_atoms(moldyn,SCALE_UP,su,TRUE,TRUE,TRUE);
1093         link_cell_shutdown(moldyn);
1094         link_cell_init(moldyn,QUIET);
1095         potential_force_calc(moldyn);
1096         y1=moldyn->energy;
1097
1098         /* restore atomic configuration + dim */
1099         memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
1100         moldyn->dim=dim;
1101
1102         /* scale down dimension and atom positions */
1103         scale_dim(moldyn,SCALE_DOWN,sd,TRUE,TRUE,TRUE);
1104         scale_atoms(moldyn,SCALE_DOWN,sd,TRUE,TRUE,TRUE);
1105         link_cell_shutdown(moldyn);
1106         link_cell_init(moldyn,QUIET);
1107         potential_force_calc(moldyn);
1108         y0=moldyn->energy;
1109         
1110         /* calculate pressure */
1111         moldyn->tp=-(y1-y0)/(2.0*dv);
1112
1113         /* restore atomic configuration */
1114         memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
1115         moldyn->dim=dim;
1116         link_cell_shutdown(moldyn);
1117         link_cell_init(moldyn,QUIET);
1118         //potential_force_calc(moldyn);
1119
1120         /* free store buffer */
1121         if(store)
1122                 free(store);
1123
1124         return moldyn->tp;
1125 }
1126
1127 double get_pressure(t_moldyn *moldyn) {
1128
1129         return moldyn->p;
1130
1131 }
1132
1133 int scale_dim(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z) {
1134
1135         t_3dvec *dim;
1136
1137         dim=&(moldyn->dim);
1138
1139         if(dir==SCALE_UP)
1140                 scale=1.0+scale;
1141
1142         if(dir==SCALE_DOWN)
1143                 scale=1.0-scale;
1144
1145         if(x) dim->x*=scale;
1146         if(y) dim->y*=scale;
1147         if(z) dim->z*=scale;
1148
1149         return 0;
1150 }
1151
1152 int scale_atoms(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z) {
1153
1154         int i;
1155         t_3dvec *r;
1156
1157         if(dir==SCALE_UP)
1158                 scale=1.0+scale;
1159
1160         if(dir==SCALE_DOWN)
1161                 scale=1.0-scale;
1162
1163         for(i=0;i<moldyn->count;i++) {
1164                 r=&(moldyn->atom[i].r);
1165                 if(x) r->x*=scale;
1166                 if(y) r->y*=scale;
1167                 if(z) r->z*=scale;
1168         }
1169
1170         return 0;
1171 }
1172
1173 int scale_volume(t_moldyn *moldyn) {
1174
1175         t_3dvec *dim,*vdim;
1176         double scale;
1177         t_linkcell *lc;
1178
1179         vdim=&(moldyn->vis.dim);
1180         dim=&(moldyn->dim);
1181         lc=&(moldyn->lc);
1182
1183         /* scaling factor */
1184         if(moldyn->pt_scale&P_SCALE_BERENDSEN) {
1185                 scale=1.0-(moldyn->p_ref-moldyn->p)*moldyn->p_tc;
1186                 scale=pow(scale,ONE_THIRD);
1187         }
1188         else {
1189                 scale=pow(moldyn->p/moldyn->p_ref,ONE_THIRD);
1190         }
1191
1192         /* scale the atoms and dimensions */
1193         scale_atoms(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
1194         scale_dim(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
1195
1196         /* visualize dimensions */
1197         if(vdim->x!=0) {
1198                 vdim->x=dim->x;
1199                 vdim->y=dim->y;
1200                 vdim->z=dim->z;
1201         }
1202
1203         /* recalculate scaled volume */
1204         moldyn->volume=dim->x*dim->y*dim->z;
1205
1206         /* adjust/reinit linkcell */
1207         if(((int)(dim->x/moldyn->cutoff)!=lc->nx)||
1208            ((int)(dim->y/moldyn->cutoff)!=lc->ny)||
1209            ((int)(dim->z/moldyn->cutoff)!=lc->nx)) {
1210                 link_cell_shutdown(moldyn);
1211                 link_cell_init(moldyn,QUIET);
1212         } else {
1213                 lc->x*=scale;
1214                 lc->y*=scale;
1215                 lc->z*=scale;
1216         }
1217
1218         return 0;
1219
1220 }
1221
1222 double e_kin_calc(t_moldyn *moldyn) {
1223
1224         int i;
1225         t_atom *atom;
1226
1227         atom=moldyn->atom;
1228         moldyn->ekin=0.0;
1229
1230         for(i=0;i<moldyn->count;i++) {
1231                 atom[i].ekin=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
1232                 moldyn->ekin+=atom[i].ekin;
1233         }
1234
1235         return moldyn->ekin;
1236 }
1237
1238 double get_total_energy(t_moldyn *moldyn) {
1239
1240         return(moldyn->ekin+moldyn->energy);
1241 }
1242
1243 t_3dvec get_total_p(t_moldyn *moldyn) {
1244
1245         t_3dvec p,p_total;
1246         int i;
1247         t_atom *atom;
1248
1249         atom=moldyn->atom;
1250
1251         v3_zero(&p_total);
1252         for(i=0;i<moldyn->count;i++) {
1253                 v3_scale(&p,&(atom[i].v),atom[i].mass);
1254                 v3_add(&p_total,&p_total,&p);
1255         }
1256
1257         return p_total;
1258 }
1259
1260 double estimate_time_step(t_moldyn *moldyn,double nn_dist) {
1261
1262         double tau;
1263
1264         /* nn_dist is the nearest neighbour distance */
1265
1266         tau=(0.05*nn_dist*moldyn->atom[0].mass)/sqrt(3.0*K_BOLTZMANN*moldyn->t);
1267
1268         return tau;     
1269 }
1270
1271 /*
1272  * numerical tricks
1273  */
1274
1275 /* linked list / cell method */
1276
1277 int link_cell_init(t_moldyn *moldyn,u8 vol) {
1278
1279         t_linkcell *lc;
1280         int i;
1281
1282         lc=&(moldyn->lc);
1283
1284         /* partitioning the md cell */
1285         lc->nx=moldyn->dim.x/moldyn->cutoff;
1286         lc->x=moldyn->dim.x/lc->nx;
1287         lc->ny=moldyn->dim.y/moldyn->cutoff;
1288         lc->y=moldyn->dim.y/lc->ny;
1289         lc->nz=moldyn->dim.z/moldyn->cutoff;
1290         lc->z=moldyn->dim.z/lc->nz;
1291         lc->cells=lc->nx*lc->ny*lc->nz;
1292
1293 #ifdef STATIC_LISTS
1294         lc->subcell=malloc(lc->cells*sizeof(int*));
1295 #else
1296         lc->subcell=malloc(lc->cells*sizeof(t_list));
1297 #endif
1298
1299         if(lc->subcell==NULL) {
1300                 perror("[moldyn] cell init (malloc)");
1301                 return -1;
1302         }
1303
1304         if(lc->cells<27)
1305                 printf("[moldyn] FATAL: less then 27 subcells!\n");
1306
1307         if(vol) {
1308 #ifdef STATIC_LISTS
1309                 printf("[moldyn] initializing 'static' linked cells (%d)\n",
1310                        lc->cells);
1311 #else
1312                 printf("[moldyn] initializing 'dynamic' linked cells (%d)\n",
1313                        lc->cells);
1314 #endif
1315                 printf("  x: %d x %f A\n",lc->nx,lc->x);
1316                 printf("  y: %d x %f A\n",lc->ny,lc->y);
1317                 printf("  z: %d x %f A\n",lc->nz,lc->z);
1318         }
1319
1320 #ifdef STATIC_LISTS
1321         /* list init */
1322         for(i=0;i<lc->cells;i++) {
1323                 lc->subcell[i]=malloc((MAX_ATOMS_PER_LIST+1)*sizeof(int));
1324                 if(lc->subcell[i]==NULL) {
1325                         perror("[moldyn] list init (malloc)");
1326                         return -1;
1327                 }
1328                 /*
1329                 if(i==0)
1330                         printf(" ---> %d malloc %p (%p)\n",
1331                                i,lc->subcell[0],lc->subcell);
1332                 */
1333         }
1334 #else
1335         for(i=0;i<lc->cells;i++)
1336                 list_init_f(&(lc->subcell[i]));
1337 #endif
1338
1339         /* update the list */
1340         link_cell_update(moldyn);
1341
1342         return 0;
1343 }
1344
1345 int link_cell_update(t_moldyn *moldyn) {
1346
1347         int count,i,j,k;
1348         int nx,ny;
1349         t_atom *atom;
1350         t_linkcell *lc;
1351 #ifdef STATIC_LISTS
1352         int p;
1353 #endif
1354
1355         atom=moldyn->atom;
1356         lc=&(moldyn->lc);
1357
1358         nx=lc->nx;
1359         ny=lc->ny;
1360
1361         for(i=0;i<lc->cells;i++)
1362 #ifdef STATIC_LISTS
1363                 memset(lc->subcell[i],0,(MAX_ATOMS_PER_LIST+1)*sizeof(int));
1364 #else
1365                 list_destroy_f(&(lc->subcell[i]));
1366 #endif
1367
1368         for(count=0;count<moldyn->count;count++) {
1369                 i=((atom[count].r.x+(moldyn->dim.x/2))/lc->x);
1370                 j=((atom[count].r.y+(moldyn->dim.y/2))/lc->y);
1371                 k=((atom[count].r.z+(moldyn->dim.z/2))/lc->z);
1372         
1373 #ifdef STATIC_LISTS
1374                 p=0;
1375                 while(lc->subcell[i+j*nx+k*nx*ny][p]!=0)
1376                         p++;
1377
1378                 if(p>=MAX_ATOMS_PER_LIST) {
1379                         printf("[moldyn] FATAL: amount of atoms too high!\n");
1380                         return -1;
1381                 }
1382
1383                 lc->subcell[i+j*nx+k*nx*ny][p]=count;
1384 #else
1385                 list_add_immediate_f(&(lc->subcell[i+j*nx+k*nx*ny]),
1386                                      &(atom[count]));
1387                 /*
1388                 if(j==0&&k==0)
1389                         printf(" ---> %d %d malloc %p (%p)\n",
1390                                i,count,lc->subcell[i].current,lc->subcell);
1391                 */
1392 #endif
1393         }
1394
1395         return 0;
1396 }
1397
1398 int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,
1399 #ifdef STATIC_LISTS
1400                               int **cell
1401 #else
1402                               t_list *cell
1403 #endif
1404                              ) {
1405
1406         t_linkcell *lc;
1407         int a;
1408         int count1,count2;
1409         int ci,cj,ck;
1410         int nx,ny,nz;
1411         int x,y,z;
1412         u8 bx,by,bz;
1413
1414         lc=&(moldyn->lc);
1415         nx=lc->nx;
1416         ny=lc->ny;
1417         nz=lc->nz;
1418         count1=1;
1419         count2=27;
1420         a=nx*ny;
1421
1422         if(i>=nx||j>=ny||k>=nz)
1423                 printf("[moldyn] WARNING: lcni %d/%d %d/%d %d/%d\n",
1424                        i,nx,j,ny,k,nz);
1425
1426         cell[0]=lc->subcell[i+j*nx+k*a];
1427         for(ci=-1;ci<=1;ci++) {
1428                 bx=0;
1429                 x=i+ci;
1430                 if((x<0)||(x>=nx)) {
1431                         x=(x+nx)%nx;
1432                         bx=1;
1433                 }
1434                 for(cj=-1;cj<=1;cj++) {
1435                         by=0;
1436                         y=j+cj;
1437                         if((y<0)||(y>=ny)) {
1438                                 y=(y+ny)%ny;
1439                                 by=1;
1440                         }
1441                         for(ck=-1;ck<=1;ck++) {
1442                                 bz=0;
1443                                 z=k+ck;
1444                                 if((z<0)||(z>=nz)) {
1445                                         z=(z+nz)%nz;
1446                                         bz=1;
1447                                 }
1448                                 if(!(ci|cj|ck)) continue;
1449                                 if(bx|by|bz) {
1450                                         cell[--count2]=lc->subcell[x+y*nx+z*a];
1451                                 }
1452                                 else {
1453                                         cell[count1++]=lc->subcell[x+y*nx+z*a];
1454                                 }
1455                         }
1456                 }
1457         }
1458
1459         lc->dnlc=count1;
1460
1461         return count1;
1462 }
1463
1464 int link_cell_shutdown(t_moldyn *moldyn) {
1465
1466         int i;
1467         t_linkcell *lc;
1468
1469         lc=&(moldyn->lc);
1470
1471         for(i=0;i<lc->cells;i++) {
1472 #ifdef STATIC_LISTS
1473                 free(lc->subcell[i]);
1474 #else
1475                 //printf(" ---> %d free %p\n",i,lc->subcell[i].start);
1476                 list_destroy_f(&(lc->subcell[i]));
1477 #endif
1478         }
1479
1480         free(lc->subcell);
1481
1482         return 0;
1483 }
1484
1485 int moldyn_add_schedule(t_moldyn *moldyn,int runs,double tau) {
1486
1487         int count;
1488         void *ptr;
1489         t_moldyn_schedule *schedule;
1490
1491         schedule=&(moldyn->schedule);
1492         count=++(schedule->total_sched);
1493
1494         ptr=realloc(schedule->runs,count*sizeof(int));
1495         if(!ptr) {
1496                 perror("[moldyn] realloc (runs)");
1497                 return -1;
1498         }
1499         schedule->runs=ptr;
1500         schedule->runs[count-1]=runs;
1501
1502         ptr=realloc(schedule->tau,count*sizeof(double));
1503         if(!ptr) {
1504                 perror("[moldyn] realloc (tau)");
1505                 return -1;
1506         }
1507         schedule->tau=ptr;
1508         schedule->tau[count-1]=tau;
1509
1510         printf("[moldyn] schedule added:\n");
1511         printf("  number: %d | runs: %d | tau: %f\n",count-1,runs,tau);
1512                                        
1513
1514         return 0;
1515 }
1516
1517 int moldyn_set_schedule_hook(t_moldyn *moldyn,set_hook hook,void *hook_params) {
1518
1519         moldyn->schedule.hook=hook;
1520         moldyn->schedule.hook_params=hook_params;
1521         
1522         return 0;
1523 }
1524
1525 /*
1526  *
1527  * 'integration of newtons equation' - algorithms
1528  *
1529  */
1530
1531 /* start the integration */
1532
1533 int moldyn_integrate(t_moldyn *moldyn) {
1534
1535         int i;
1536         unsigned int e,m,s,v,p,t,a;
1537         t_3dvec momentum;
1538         t_moldyn_schedule *sched;
1539         t_atom *atom;
1540         int fd;
1541         char dir[128];
1542         double ds;
1543         double energy_scale;
1544         struct timeval t1,t2;
1545         //double tp;
1546
1547         sched=&(moldyn->schedule);
1548         atom=moldyn->atom;
1549
1550         /* initialize linked cell method */
1551         link_cell_init(moldyn,VERBOSE);
1552
1553         /* logging & visualization */
1554         e=moldyn->ewrite;
1555         m=moldyn->mwrite;
1556         s=moldyn->swrite;
1557         v=moldyn->vwrite;
1558         a=moldyn->awrite;
1559         p=moldyn->pwrite;
1560         t=moldyn->twrite;
1561
1562         /* sqaure of some variables */
1563         moldyn->tau_square=moldyn->tau*moldyn->tau;
1564
1565         /* get current time */
1566         gettimeofday(&t1,NULL);
1567
1568         /* calculate initial forces */
1569         potential_force_calc(moldyn);
1570 #ifdef DEBUG
1571 //return 0;
1572 #endif
1573
1574         /* some stupid checks before we actually start calculating bullshit */
1575         if(moldyn->cutoff>0.5*moldyn->dim.x)
1576                 printf("[moldyn] WARNING: cutoff > 0.5 x dim.x\n");
1577         if(moldyn->cutoff>0.5*moldyn->dim.y)
1578                 printf("[moldyn] WARNING: cutoff > 0.5 x dim.y\n");
1579         if(moldyn->cutoff>0.5*moldyn->dim.z)
1580                 printf("[moldyn] WARNING: cutoff > 0.5 x dim.z\n");
1581         if(moldyn->count) {
1582                 ds=0.5*atom[0].f.x*moldyn->tau_square/atom[0].mass;
1583                 if(ds>0.05*moldyn->nnd)
1584                 printf("[moldyn] WARNING: forces too high / tau too small!\n");
1585         }
1586
1587         /* zero absolute time */
1588         // should have right values!
1589         //moldyn->time=0.0;
1590         //moldyn->total_steps=0;
1591
1592         /* debugging, ignore */
1593         moldyn->debug=0;
1594
1595         /* tell the world */
1596         printf("[moldyn] integration start, go get a coffee ...\n");
1597
1598         /* executing the schedule */
1599         sched->count=0;
1600         while(sched->count<sched->total_sched) {
1601
1602                 /* setting amount of runs and finite time step size */
1603                 moldyn->tau=sched->tau[sched->count];
1604                 moldyn->tau_square=moldyn->tau*moldyn->tau;
1605                 moldyn->time_steps=sched->runs[sched->count];
1606
1607                 /* energy scaling factor (might change!) */
1608                 energy_scale=moldyn->count*EV;
1609
1610         /* integration according to schedule */
1611
1612         for(i=0;i<moldyn->time_steps;i++) {
1613
1614                 /* integration step */
1615                 moldyn->integrate(moldyn);
1616
1617                 /* calculate kinetic energy, temperature and pressure */
1618                 e_kin_calc(moldyn);
1619                 temperature_calc(moldyn);
1620                 virial_sum(moldyn);
1621                 pressure_calc(moldyn);
1622                 /*
1623                 thermodynamic_pressure_calc(moldyn);
1624                 printf("\n\nDEBUG: numeric pressure calc: %f\n\n",
1625                        moldyn->tp/BAR);
1626                 */
1627
1628                 /* calculate fluctuations + averages */
1629                 average_and_fluctuation_calc(moldyn);
1630
1631                 /* p/t scaling */
1632                 if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT))
1633                         scale_velocity(moldyn,FALSE);
1634                 if(moldyn->pt_scale&(P_SCALE_BERENDSEN|P_SCALE_DIRECT))
1635                         scale_volume(moldyn);
1636
1637                 /* check for log & visualization */
1638                 if(e) {
1639                         if(!(moldyn->total_steps%e))
1640                                 dprintf(moldyn->efd,
1641                                         "%f %f %f %f\n",
1642                                         moldyn->time,moldyn->ekin/energy_scale,
1643                                         moldyn->energy/energy_scale,
1644                                         get_total_energy(moldyn)/energy_scale);
1645                 }
1646                 if(m) {
1647                         if(!(moldyn->total_steps%m)) {
1648                                 momentum=get_total_p(moldyn);
1649                                 dprintf(moldyn->mfd,
1650                                         "%f %f %f %f %f\n",moldyn->time,
1651                                         momentum.x,momentum.y,momentum.z,
1652                                         v3_norm(&momentum));
1653                         }
1654                 }
1655                 if(p) {
1656                         if(!(moldyn->total_steps%p)) {
1657                                 dprintf(moldyn->pfd,
1658                                         "%f %f %f %f %f %f %f\n",moldyn->time,
1659                                          moldyn->p/BAR,moldyn->p_avg/BAR,
1660                                          moldyn->gp/BAR,moldyn->gp_avg/BAR,
1661                                          moldyn->tp/BAR,moldyn->tp_avg/BAR);
1662                         }
1663                 }
1664                 if(t) {
1665                         if(!(moldyn->total_steps%t)) {
1666                                 dprintf(moldyn->tfd,
1667                                         "%f %f %f\n",
1668                                         moldyn->time,moldyn->t,moldyn->t_avg);
1669                         }
1670                 }
1671                 if(v) {
1672                         if(!(moldyn->total_steps%v)) {
1673                                 dprintf(moldyn->vfd,
1674                                         "%f %f\n",moldyn->time,moldyn->volume);
1675                         }
1676                 }
1677                 if(s) {
1678                         if(!(moldyn->total_steps%s)) {
1679                                 snprintf(dir,128,"%s/s-%07.f.save",
1680                                          moldyn->vlsdir,moldyn->time);
1681                                 fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT,
1682                                         S_IRUSR|S_IWUSR);
1683                                 if(fd<0) perror("[moldyn] save fd open");
1684                                 else {
1685                                         write(fd,moldyn,sizeof(t_moldyn));
1686                                         write(fd,moldyn->atom,
1687                                               moldyn->count*sizeof(t_atom));
1688                                 }
1689                                 close(fd);
1690                         }       
1691                 }
1692                 if(a) {
1693                         if(!(moldyn->total_steps%a)) {
1694                                 visual_atoms(moldyn);
1695                         }
1696                 }
1697
1698                 /* display progress */
1699                 //if(!(moldyn->total_steps%10)) {
1700                         /* get current time */
1701                         gettimeofday(&t2,NULL);
1702
1703 printf("\rsched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)",
1704        sched->count,i,moldyn->total_steps,
1705        moldyn->t,moldyn->t_avg,
1706        moldyn->p/BAR,moldyn->p_avg/BAR,
1707        //moldyn->p/BAR,(moldyn->p-2.0*moldyn->ekin/(3.0*moldyn->volume))/BAR,
1708        moldyn->volume,
1709        (int)(t2.tv_sec-t1.tv_sec));
1710
1711                         fflush(stdout);
1712
1713                         /* copy over time */
1714                         t1=t2;
1715                 //}
1716
1717                 /* increase absolute time */
1718                 moldyn->time+=moldyn->tau;
1719                 moldyn->total_steps+=1;
1720
1721         }
1722
1723                 /* check for hooks */
1724                 if(sched->hook) {
1725                         printf("\n ## schedule hook %d start ##\n",
1726                                sched->count);
1727                         sched->hook(moldyn,sched->hook_params);
1728                         printf(" ## schedule hook end ##\n");
1729                 }
1730
1731                 /* increase the schedule counter */
1732                 sched->count+=1;
1733
1734         }
1735
1736         return 0;
1737 }
1738
1739 /* velocity verlet */
1740
1741 int velocity_verlet(t_moldyn *moldyn) {
1742
1743         int i,count;
1744         double tau,tau_square,h;
1745         t_3dvec delta;
1746         t_atom *atom;
1747
1748         atom=moldyn->atom;
1749         count=moldyn->count;
1750         tau=moldyn->tau;
1751         tau_square=moldyn->tau_square;
1752
1753         for(i=0;i<count;i++) {
1754                 /* check whether fixed atom */
1755                 if(atom[i].attr&ATOM_ATTR_FP)
1756                         continue;
1757                 /* new positions */
1758                 h=0.5/atom[i].mass;
1759                 v3_scale(&delta,&(atom[i].v),tau);
1760                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
1761                 v3_scale(&delta,&(atom[i].f),h*tau_square);
1762                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
1763                 check_per_bound(moldyn,&(atom[i].r));
1764
1765                 /* velocities [actually v(t+tau/2)] */
1766                 v3_scale(&delta,&(atom[i].f),h*tau);
1767                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
1768         }
1769
1770         /* criticial check */
1771         moldyn_bc_check(moldyn);
1772
1773         /* neighbour list update */
1774         link_cell_update(moldyn);
1775
1776         /* forces depending on chosen potential */
1777         potential_force_calc(moldyn);
1778
1779         for(i=0;i<count;i++) {
1780                 /* check whether fixed atom */
1781                 if(atom[i].attr&ATOM_ATTR_FP)
1782                         continue;
1783                 /* again velocities [actually v(t+tau)] */
1784                 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
1785                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
1786         }
1787
1788         return 0;
1789 }
1790
1791
1792 /*
1793  *
1794  * potentials & corresponding forces & virial routine
1795  * 
1796  */
1797
1798 /* generic potential and force calculation */
1799
1800 int potential_force_calc(t_moldyn *moldyn) {
1801
1802         int i,j,k,count;
1803         t_atom *itom,*jtom,*ktom;
1804         t_virial *virial;
1805         t_linkcell *lc;
1806 #ifdef STATIC_LISTS
1807         int *neighbour_i[27];
1808         int p,q;
1809         t_atom *atom;
1810 #else
1811         t_list neighbour_i[27];
1812         t_list neighbour_i2[27];
1813         t_list *this,*that;
1814 #endif
1815         u8 bc_ij,bc_ik;
1816         int dnlc;
1817
1818         count=moldyn->count;
1819         itom=moldyn->atom;
1820         lc=&(moldyn->lc);
1821 #ifdef STATIC_LISTS
1822         atom=moldyn->atom;
1823 #endif
1824
1825         /* reset energy */
1826         moldyn->energy=0.0;
1827
1828         /* reset global virial */
1829         memset(&(moldyn->gvir),0,sizeof(t_virial));
1830
1831         /* reset force, site energy and virial of every atom */
1832         for(i=0;i<count;i++) {
1833
1834                 /* reset force */
1835                 v3_zero(&(itom[i].f));
1836
1837                 /* reset virial */
1838                 virial=(&(itom[i].virial));
1839                 virial->xx=0.0;
1840                 virial->yy=0.0;
1841                 virial->zz=0.0;
1842                 virial->xy=0.0;
1843                 virial->xz=0.0;
1844                 virial->yz=0.0;
1845         
1846                 /* reset site energy */
1847                 itom[i].e=0.0;
1848
1849         }
1850
1851         /* get energy, force and virial of every atom */
1852
1853         /* first (and only) loop over atoms i */
1854         for(i=0;i<count;i++) {
1855
1856                 /* single particle potential/force */
1857                 if(itom[i].attr&ATOM_ATTR_1BP)
1858                         if(moldyn->func1b)
1859                                 moldyn->func1b(moldyn,&(itom[i]));
1860
1861                 if(!(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)))
1862                         continue;
1863
1864                 /* 2 body pair potential/force */
1865         
1866                 link_cell_neighbour_index(moldyn,
1867                                           (itom[i].r.x+moldyn->dim.x/2)/lc->x,
1868                                           (itom[i].r.y+moldyn->dim.y/2)/lc->y,
1869                                           (itom[i].r.z+moldyn->dim.z/2)/lc->z,
1870                                           neighbour_i);
1871
1872                 dnlc=lc->dnlc;
1873
1874                 /* first loop over atoms j */
1875                 if(moldyn->func2b) {
1876                         for(j=0;j<27;j++) {
1877
1878                                 bc_ij=(j<dnlc)?0:1;
1879 #ifdef STATIC_LISTS
1880                                 p=0;
1881
1882                                 while(neighbour_i[j][p]!=0) {
1883
1884                                         jtom=&(atom[neighbour_i[j][p]]);
1885                                         p++;
1886
1887                                         if(jtom==&(itom[i]))
1888                                                 continue;
1889
1890                                         if((jtom->attr&ATOM_ATTR_2BP)&
1891                                            (itom[i].attr&ATOM_ATTR_2BP)) {
1892                                                 moldyn->func2b(moldyn,
1893                                                                &(itom[i]),
1894                                                                jtom,
1895                                                                bc_ij);
1896                                         }
1897                                 }
1898 #else
1899                                 this=&(neighbour_i[j]);
1900                                 list_reset_f(this);
1901
1902                                 if(this->start==NULL)
1903                                         continue;
1904
1905                                 do {
1906                                         jtom=this->current->data;
1907
1908                                         if(jtom==&(itom[i]))
1909                                                 continue;
1910
1911                                         if((jtom->attr&ATOM_ATTR_2BP)&
1912                                            (itom[i].attr&ATOM_ATTR_2BP)) {
1913                                                 moldyn->func2b(moldyn,
1914                                                                &(itom[i]),
1915                                                                jtom,
1916                                                                bc_ij);
1917                                         }
1918                                 } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
1919 #endif
1920
1921                         }
1922                 }
1923
1924                 /* 3 body potential/force */
1925
1926                 if(!(itom[i].attr&ATOM_ATTR_3BP))
1927                         continue;
1928
1929                 /* copy the neighbour lists */
1930 #ifdef STATIC_LISTS
1931                 /* no copy needed for static lists */
1932 #else
1933                 memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
1934 #endif
1935
1936                 /* second loop over atoms j */
1937                 for(j=0;j<27;j++) {
1938
1939                         bc_ij=(j<dnlc)?0:1;
1940 #ifdef STATIC_LISTS
1941                         p=0;
1942
1943                         while(neighbour_i[j][p]!=0) {
1944
1945                                 jtom=&(atom[neighbour_i[j][p]]);
1946                                 p++;
1947 #else
1948                         this=&(neighbour_i[j]);
1949                         list_reset_f(this);
1950
1951                         if(this->start==NULL)
1952                                 continue;
1953
1954                         do {
1955
1956                                 jtom=this->current->data;
1957 #endif
1958
1959                                 if(jtom==&(itom[i]))
1960                                         continue;
1961
1962                                 if(!(jtom->attr&ATOM_ATTR_3BP))
1963                                         continue;
1964
1965                                 /* reset 3bp run */
1966                                 moldyn->run3bp=1;
1967
1968                                 if(moldyn->func3b_j1)
1969                                         moldyn->func3b_j1(moldyn,
1970                                                           &(itom[i]),
1971                                                           jtom,
1972                                                           bc_ij);
1973
1974                                 /* in first j loop, 3bp run can be skipped */
1975                                 if(!(moldyn->run3bp))
1976                                         continue;
1977                         
1978                                 /* first loop over atoms k */
1979                                 if(moldyn->func3b_k1) {
1980
1981                                 for(k=0;k<27;k++) {
1982
1983                                         bc_ik=(k<dnlc)?0:1;
1984 #ifdef STATIC_LISTS
1985                                         q=0;
1986
1987                                         while(neighbour_i[j][q]!=0) {
1988
1989                                                 ktom=&(atom[neighbour_i[k][q]]);
1990                                                 q++;
1991 #else
1992                                         that=&(neighbour_i2[k]);
1993                                         list_reset_f(that);
1994                                         
1995                                         if(that->start==NULL)
1996                                                 continue;
1997
1998                                         do {
1999                                                 ktom=that->current->data;
2000 #endif
2001
2002                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
2003                                                         continue;
2004
2005                                                 if(ktom==jtom)
2006                                                         continue;
2007
2008                                                 if(ktom==&(itom[i]))
2009                                                         continue;
2010
2011                                                 moldyn->func3b_k1(moldyn,
2012                                                                   &(itom[i]),
2013                                                                   jtom,
2014                                                                   ktom,
2015                                                                   bc_ik|bc_ij);
2016 #ifdef STATIC_LISTS
2017                                         }
2018 #else
2019                                         } while(list_next_f(that)!=\
2020                                                 L_NO_NEXT_ELEMENT);
2021 #endif
2022
2023                                 }
2024
2025                                 }
2026
2027                                 if(moldyn->func3b_j2)
2028                                         moldyn->func3b_j2(moldyn,
2029                                                           &(itom[i]),
2030                                                           jtom,
2031                                                           bc_ij);
2032
2033                                 /* second loop over atoms k */
2034                                 if(moldyn->func3b_k2) {
2035
2036                                 for(k=0;k<27;k++) {
2037
2038                                         bc_ik=(k<dnlc)?0:1;
2039 #ifdef STATIC_LISTS
2040                                         q=0;
2041
2042                                         while(neighbour_i[j][q]!=0) {
2043
2044                                                 ktom=&(atom[neighbour_i[k][q]]);
2045                                                 q++;
2046 #else
2047                                         that=&(neighbour_i2[k]);
2048                                         list_reset_f(that);
2049                                         
2050                                         if(that->start==NULL)
2051                                                 continue;
2052
2053                                         do {
2054                                                 ktom=that->current->data;
2055 #endif
2056
2057                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
2058                                                         continue;
2059
2060                                                 if(ktom==jtom)
2061                                                         continue;
2062
2063                                                 if(ktom==&(itom[i]))
2064                                                         continue;
2065
2066                                                 moldyn->func3b_k2(moldyn,
2067                                                                   &(itom[i]),
2068                                                                   jtom,
2069                                                                   ktom,
2070                                                                   bc_ik|bc_ij);
2071
2072 #ifdef STATIC_LISTS
2073                                         }
2074 #else
2075                                         } while(list_next_f(that)!=\
2076                                                 L_NO_NEXT_ELEMENT);
2077 #endif
2078
2079                                 }
2080                                 
2081                                 }
2082
2083                                 /* 2bp post function */
2084                                 if(moldyn->func3b_j3) {
2085                                         moldyn->func3b_j3(moldyn,
2086                                                           &(itom[i]),
2087                                                           jtom,bc_ij);
2088                                 }
2089 #ifdef STATIC_LISTS
2090                         }
2091 #else
2092                         } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
2093 #endif
2094                 
2095                 }
2096                 
2097 #ifdef DEBUG
2098         //printf("\n\n");
2099 #endif
2100 #ifdef VDEBUG
2101         printf("\n\n");
2102 #endif
2103
2104         }
2105
2106 #ifdef DEBUG
2107         //printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
2108         if(moldyn->time>DSTART&&moldyn->time<DEND) {
2109                 printf("force:\n");
2110                 printf("  x: %0.40f\n",moldyn->atom[DATOM].f.x);
2111                 printf("  y: %0.40f\n",moldyn->atom[DATOM].f.y);
2112                 printf("  z: %0.40f\n",moldyn->atom[DATOM].f.z);
2113         }
2114 #endif
2115
2116         /* some postprocessing */
2117         for(i=0;i<count;i++) {
2118                 /* calculate global virial */
2119                 moldyn->gvir.xx+=itom[i].r.x*itom[i].f.x;
2120                 moldyn->gvir.yy+=itom[i].r.y*itom[i].f.y;
2121                 moldyn->gvir.zz+=itom[i].r.z*itom[i].f.z;
2122                 moldyn->gvir.xy+=itom[i].r.y*itom[i].f.x;
2123                 moldyn->gvir.xz+=itom[i].r.z*itom[i].f.x;
2124                 moldyn->gvir.yz+=itom[i].r.z*itom[i].f.y;
2125
2126                 /* check forces regarding the given timestep */
2127                 if(v3_norm(&(itom[i].f))>\
2128                    0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square)
2129                         printf("[moldyn] WARNING: pfc (high force: atom %d)\n",
2130                                i);
2131         }
2132
2133         return 0;
2134 }
2135
2136 /*
2137  * virial calculation
2138  */
2139
2140 //inline int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) {
2141 int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) {
2142
2143         a->virial.xx+=f->x*d->x;
2144         a->virial.yy+=f->y*d->y;
2145         a->virial.zz+=f->z*d->z;
2146         a->virial.xy+=f->x*d->y;
2147         a->virial.xz+=f->x*d->z;
2148         a->virial.yz+=f->y*d->z;
2149
2150         return 0;
2151 }
2152
2153 /*
2154  * periodic boundary checking
2155  */
2156
2157 //inline int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
2158 int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
2159         
2160         double x,y,z;
2161         t_3dvec *dim;
2162
2163         dim=&(moldyn->dim);
2164
2165         x=dim->x/2;
2166         y=dim->y/2;
2167         z=dim->z/2;
2168
2169         if(moldyn->status&MOLDYN_STAT_PBX) {
2170                 if(a->x>=x) a->x-=dim->x;
2171                 else if(-a->x>x) a->x+=dim->x;
2172         }
2173         if(moldyn->status&MOLDYN_STAT_PBY) {
2174                 if(a->y>=y) a->y-=dim->y;
2175                 else if(-a->y>y) a->y+=dim->y;
2176         }
2177         if(moldyn->status&MOLDYN_STAT_PBZ) {
2178                 if(a->z>=z) a->z-=dim->z;
2179                 else if(-a->z>z) a->z+=dim->z;
2180         }
2181
2182         return 0;
2183 }
2184         
2185 /*
2186  * debugging / critical check functions
2187  */
2188
2189 int moldyn_bc_check(t_moldyn *moldyn) {
2190
2191         t_atom *atom;
2192         t_3dvec *dim;
2193         int i;
2194         double x;
2195         u8 byte;
2196         int j,k;
2197
2198         atom=moldyn->atom;
2199         dim=&(moldyn->dim);
2200         x=dim->x/2;
2201
2202         for(i=0;i<moldyn->count;i++) {
2203                 if(atom[i].r.x>=dim->x/2||-atom[i].r.x>dim->x/2) {
2204                         printf("FATAL: atom %d: x: %.20f (%.20f)\n",
2205                                i,atom[i].r.x,dim->x/2);
2206                         printf("diagnostic:\n");
2207                         printf("-----------\natom.r.x:\n");
2208                         for(j=0;j<8;j++) {
2209                                 memcpy(&byte,(u8 *)(&(atom[i].r.x))+j,1);
2210                                 for(k=0;k<8;k++)
2211                                         printf("%d%c",
2212                                         ((byte)&(1<<k))?1:0,
2213                                         (k==7)?'\n':'|');
2214                         }
2215                         printf("---------------\nx=dim.x/2:\n");
2216                         for(j=0;j<8;j++) {
2217                                 memcpy(&byte,(u8 *)(&x)+j,1);
2218                                 for(k=0;k<8;k++)
2219                                         printf("%d%c",
2220                                         ((byte)&(1<<k))?1:0,
2221                                         (k==7)?'\n':'|');
2222                         }
2223                         if(atom[i].r.x==x) printf("the same!\n");
2224                         else printf("different!\n");
2225                 }
2226                 if(atom[i].r.y>=dim->y/2||-atom[i].r.y>dim->y/2)
2227                         printf("FATAL: atom %d: y: %.20f (%.20f)\n",
2228                                i,atom[i].r.y,dim->y/2);
2229                 if(atom[i].r.z>=dim->z/2||-atom[i].r.z>dim->z/2)
2230                         printf("FATAL: atom %d: z: %.20f (%.20f)\n",
2231                                i,atom[i].r.z,dim->z/2);
2232         }
2233
2234         return 0;
2235 }
2236
2237 /*
2238  * restore function
2239  */
2240
2241 int moldyn_read_save_file(t_moldyn *moldyn,char *file) {
2242
2243         int fd;
2244         int cnt,size;
2245         int fsize;
2246         int corr;
2247
2248         fd=open(file,O_RDONLY);
2249         if(fd<0) {
2250                 perror("[moldyn] load save file open");
2251                 return fd;
2252         }
2253
2254         fsize=lseek(fd,0,SEEK_END);
2255         lseek(fd,0,SEEK_SET);
2256
2257         size=sizeof(t_moldyn);
2258
2259         while(size) {
2260                 cnt=read(fd,moldyn,size);
2261                 if(cnt<0) {
2262                         perror("[moldyn] load save file read (moldyn)");
2263                         return cnt;
2264                 }
2265                 size-=cnt;
2266         }
2267
2268         size=moldyn->count*sizeof(t_atom);
2269
2270         /* correcting possible atom data offset */
2271         corr=0;
2272         if(fsize!=sizeof(t_moldyn)+size) {
2273                 corr=fsize-sizeof(t_moldyn)-size;
2274                 printf("[moldyn] WARNING: lsf (illegal file size)\n");
2275                 printf("  moifying offset:\n");
2276                 printf("  - current pos: %d\n",sizeof(t_moldyn));
2277                 printf("  - atom size: %d\n",size);
2278                 printf("  - file size: %d\n",fsize);
2279                 printf("  => correction: %d\n",corr);
2280                 lseek(fd,corr,SEEK_CUR);
2281         }
2282
2283         moldyn->atom=(t_atom *)malloc(size);
2284         if(moldyn->atom==NULL) {
2285                 perror("[moldyn] load save file malloc (atoms)");
2286                 return -1;
2287         }
2288
2289         while(size) {
2290                 cnt=read(fd,moldyn->atom,size);
2291                 if(cnt<0) {
2292                         perror("[moldyn] load save file read (atoms)");
2293                         return cnt;
2294                 }
2295                 size-=cnt;
2296         }
2297
2298         // hooks etc ...
2299
2300         return 0;
2301 }
2302
2303 int moldyn_free_save_file(t_moldyn *moldyn) {
2304
2305         free(moldyn->atom);
2306
2307         return 0;
2308 }
2309
2310 int moldyn_load(t_moldyn *moldyn) {
2311
2312         // later ...
2313
2314         return 0;
2315 }
2316
2317 /*
2318  * function to find/callback all combinations of 2 body bonds
2319  */
2320
2321 int process_2b_bonds(t_moldyn *moldyn,void *data,
2322                      int (*process)(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
2323                                     void *data,u8 bc)) {
2324
2325         t_linkcell *lc;
2326 #ifdef STATIC_LISTS
2327         int *neighbour[27];
2328         int p;
2329 #else
2330         t_list neighbour[27];
2331 #endif
2332         u8 bc;
2333         t_atom *itom,*jtom;
2334         int i,j;
2335         t_list *this;
2336
2337         lc=&(moldyn->lc);
2338         itom=moldyn->atom;
2339         
2340         for(i=0;i<moldyn->count;i++) {
2341                 /* neighbour indexing */
2342                 link_cell_neighbour_index(moldyn,
2343                                           (itom[i].r.x+moldyn->dim.x/2)/lc->x,
2344                                           (itom[i].r.y+moldyn->dim.y/2)/lc->x,
2345                                           (itom[i].r.z+moldyn->dim.z/2)/lc->x,
2346                                           neighbour);
2347
2348                 for(j=0;j<27;j++) {
2349
2350                         bc=(j<lc->dnlc)?0:1;
2351
2352 #ifdef STATIC_LISTS
2353                         p=0;
2354
2355                         while(neighbour[j][p]!=0) {
2356
2357                                 jtom=&(moldyn->atom[neighbour[j][p]]);
2358                                 p++;
2359 #else
2360                         this=&(neighbour[j]);
2361                         list_reset_f(this);
2362
2363                         if(this->start==NULL)
2364                                 continue;
2365
2366                         do {
2367
2368                                 jtom=this->current->data;
2369 #endif
2370
2371                                 /* process bond */
2372                                 process(moldyn,&(itom[i]),jtom,data,bc);
2373
2374 #ifdef STATIC_LISTS
2375                         }
2376 #else
2377                         } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
2378 #endif
2379                 }
2380         }
2381
2382         return 0;
2383
2384 }
2385
2386 /*
2387  * post processing functions
2388  */
2389
2390 int get_line(int fd,char *line,int max) {
2391
2392         int count,ret;
2393
2394         count=0;
2395
2396         while(1) {
2397                 if(count==max) return count;
2398                 ret=read(fd,line+count,1);
2399                 if(ret<=0) return ret;
2400                 if(line[count]=='\n') {
2401                         memset(line+count,0,max-count-1);
2402                         //line[count]='\0';
2403                         return count+1;
2404                 }
2405                 count+=1;
2406         }
2407 }
2408
2409 int pair_correlation_init(t_moldyn *moldyn,double dr) {
2410
2411         
2412         return 0;
2413 }
2414
2415 int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) {
2416
2417         int i;
2418         t_atom *atom;
2419         t_3dvec dist;
2420         double d2;
2421         int a_cnt;
2422         int b_cnt;
2423
2424         atom=moldyn->atom;
2425         dc[0]=0;
2426         dc[1]=0;
2427         dc[2]=0;
2428         a_cnt=0;
2429         b_cnt=0;
2430
2431         for(i=0;i<moldyn->count;i++) {
2432
2433                 v3_sub(&dist,&(atom[i].r),&(atom[i].r_0));
2434                 check_per_bound(moldyn,&dist);
2435                 d2=v3_absolute_square(&dist);
2436
2437                 if(atom[i].brand) {
2438                         b_cnt+=1;
2439                         dc[1]+=d2;
2440                 }
2441                 else {
2442                         a_cnt+=1;
2443                         dc[0]+=d2;
2444                 }
2445
2446                 dc[2]+=d2;
2447         }
2448
2449         dc[0]*=(1.0/(6.0*moldyn->time*a_cnt));
2450         dc[1]*=(1.0/(6.0*moldyn->time*b_cnt));
2451         dc[2]*=(1.0/(6.0*moldyn->time*moldyn->count));
2452                 
2453         return 0;
2454 }
2455
2456 int bonding_analyze(t_moldyn *moldyn,double *cnt) {
2457
2458         return 0;
2459 }
2460
2461 int calculate_pair_correlation_process(t_moldyn *moldyn,t_atom *itom,
2462                                        t_atom *jtom,void *data,u8 bc) {
2463
2464         t_3dvec dist;
2465         double d;
2466         int s;
2467         t_pcc *pcc;
2468
2469         /* only count pairs once,
2470          * skip same atoms */
2471         if(itom->tag>=jtom->tag)
2472                 return 0;
2473
2474         /*
2475          * pair correlation calc
2476          */
2477
2478         /* get pcc data */
2479         pcc=data;
2480
2481         /* distance */
2482         v3_sub(&dist,&(jtom->r),&(itom->r));
2483         if(bc) check_per_bound(moldyn,&dist);
2484         d=v3_absolute_square(&dist);
2485
2486         /* ignore if greater cutoff */
2487         if(d>moldyn->cutoff_square)
2488                 return 0;
2489
2490         /* fill the slots */
2491         d=sqrt(d);
2492         s=(int)(d/pcc->dr);
2493
2494         /* should never happen but it does 8) -
2495          * related to -ffloat-store problem! */
2496         if(s>=pcc->o1) {
2497                 printf("[moldyn] WARNING: pcc (%d/%d)",
2498                        s,pcc->o1);
2499                 printf("\n");
2500                 s=pcc->o1-1;
2501         }
2502
2503         if(itom->brand!=jtom->brand) {
2504                 /* mixed */
2505                 pcc->stat[s]+=1;
2506         }
2507         else {
2508                 /* type a - type a bonds */
2509                 if(itom->brand==0)
2510                         pcc->stat[s+pcc->o1]+=1;
2511                 else
2512                 /* type b - type b bonds */
2513                         pcc->stat[s+pcc->o2]+=1;
2514         }
2515
2516         return 0;
2517 }
2518
2519 int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
2520
2521         t_pcc pcc;
2522         double norm;
2523         int i;
2524
2525         pcc.dr=dr;
2526         pcc.o1=moldyn->cutoff/dr;
2527         pcc.o2=2*pcc.o1;
2528
2529         if(pcc.o1*dr<=moldyn->cutoff)
2530                 printf("[moldyn] WARNING: pcc (low #slots)\n");
2531
2532         printf("[moldyn] pair correlation calc info:\n");
2533         printf("  time: %f\n",moldyn->time);
2534         printf("  count: %d\n",moldyn->count);
2535         printf("  cutoff: %f\n",moldyn->cutoff);
2536         printf("  temperature: cur=%f avg=%f\n",moldyn->t,moldyn->t_avg);
2537
2538         if(ptr!=NULL) {
2539                 pcc.stat=(double *)ptr;
2540         }
2541         else {
2542                 pcc.stat=(double *)malloc(3*pcc.o1*sizeof(double));
2543                 if(pcc.stat==NULL) {
2544                         perror("[moldyn] pair correlation malloc");
2545                         return -1;
2546                 }
2547         }
2548
2549         memset(pcc.stat,0,3*pcc.o1*sizeof(double));
2550
2551         /* process */
2552         process_2b_bonds(moldyn,&pcc,calculate_pair_correlation_process);
2553
2554         /* normalization */
2555         for(i=1;i<pcc.o1;i++) {
2556                  // normalization: 4 pi r^2 dr
2557                  // here: not double counting pairs -> 2 pi r r dr
2558                  // ... and actually it's a constant times r^2
2559                 norm=i*i*dr*dr;
2560                 pcc.stat[i]/=norm;
2561                 pcc.stat[pcc.o1+i]/=norm;
2562                 pcc.stat[pcc.o2+i]/=norm;
2563         }
2564         /* */
2565
2566         if(ptr==NULL) {
2567                 /* todo: store/print pair correlation function */
2568                 free(pcc.stat);
2569         }
2570
2571         return 0;
2572 }
2573
2574 int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
2575                          void *data,u8 bc) {
2576
2577         t_ba *ba;
2578         t_3dvec dist;
2579         double d;
2580
2581         if(itom->tag>=jtom->tag)
2582                 return 0;
2583
2584         /* distance */
2585         v3_sub(&dist,&(jtom->r),&(itom->r));
2586         if(bc) check_per_bound(moldyn,&dist);
2587         d=v3_absolute_square(&dist);
2588
2589         /* ignore if greater or equal cutoff */
2590         if(d>moldyn->cutoff_square)
2591                 return 0;
2592
2593         /* check for potential bond */
2594         if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE)
2595                 return 0;
2596
2597         /* now count this bonding ... */
2598         ba=data;
2599
2600         /* increase total bond counter
2601          * ... double counting!
2602          */
2603         ba->tcnt+=2;
2604
2605         if(itom->brand==0)
2606                 ba->acnt[jtom->tag]+=1;
2607         else
2608                 ba->bcnt[jtom->tag]+=1;
2609         
2610         if(jtom->brand==0)
2611                 ba->acnt[itom->tag]+=1;
2612         else
2613                 ba->bcnt[itom->tag]+=1;
2614
2615         return 0;
2616 }
2617
2618 int bond_analyze(t_moldyn *moldyn,double *quality) {
2619
2620         // by now: # bonds of type 'a-4b' and 'b-4a' / # bonds total
2621
2622         int qcnt;
2623         int ccnt,cset;
2624         t_ba ba;
2625         int i;
2626         t_atom *atom;
2627
2628         ba.acnt=malloc(moldyn->count*sizeof(int));
2629         if(ba.acnt==NULL) {
2630                 perror("[moldyn] bond analyze malloc (a)");
2631                 return -1;
2632         }
2633         memset(ba.acnt,0,moldyn->count*sizeof(int));
2634
2635         ba.bcnt=malloc(moldyn->count*sizeof(int));
2636         if(ba.bcnt==NULL) {
2637                 perror("[moldyn] bond analyze malloc (b)");
2638                 return -1;
2639         }
2640         memset(ba.bcnt,0,moldyn->count*sizeof(int));
2641
2642         ba.tcnt=0;
2643         qcnt=0;
2644         ccnt=0;
2645         cset=0;
2646
2647         atom=moldyn->atom;
2648
2649         process_2b_bonds(moldyn,&ba,bond_analyze_process);
2650
2651         for(i=0;i<moldyn->count;i++) {
2652                 if(atom[i].brand==0) {
2653                         if((ba.acnt[i]==0)&(ba.bcnt[i]==4))
2654                                 qcnt+=4;
2655                 }
2656                 else {
2657                         if((ba.acnt[i]==4)&(ba.bcnt[i]==0)) {
2658                                 qcnt+=4;
2659                                 ccnt+=1;
2660                         }
2661                         cset+=1;
2662                 }
2663         }
2664
2665         printf("[moldyn] bond analyze: c_cnt=%d | set=%d\n",ccnt,cset);
2666         printf("[moldyn] bond analyze: q_cnt=%d | tot=%d\n",qcnt,ba.tcnt);
2667
2668         if(quality) {
2669                 quality[0]=1.0*ccnt/cset;
2670                 quality[1]=1.0*qcnt/ba.tcnt;
2671         }
2672         else {
2673                 printf("[moldyn] bond analyze: c_bnd_q=%f\n",1.0*qcnt/ba.tcnt);
2674                 printf("[moldyn] bond analyze:   tot_q=%f\n",1.0*qcnt/ba.tcnt);
2675         }
2676
2677         return 0;
2678 }
2679
2680 /*
2681  * visualization code
2682  */
2683
2684 int visual_init(t_moldyn *moldyn,char *filebase) {
2685
2686         strncpy(moldyn->vis.fb,filebase,128);
2687
2688         return 0;
2689 }
2690
2691 int visual_bonds_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
2692                          void *data,u8 bc) {
2693
2694         t_vb *vb;
2695
2696         vb=data;
2697
2698         if(itom->tag>=jtom->tag)
2699                 return 0;
2700         
2701         if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE)
2702                 return 0;
2703
2704         if((itom->attr&ATOM_ATTR_VB)|(jtom->attr&ATOM_ATTR_VB))
2705                 dprintf(vb->fd,"# [B] %f %f %f %f %f %f\n",
2706                         itom->r.x,itom->r.y,itom->r.z,
2707                         jtom->r.x,jtom->r.y,jtom->r.z);
2708
2709         return 0;
2710 }
2711
2712 int visual_atoms(t_moldyn *moldyn) {
2713
2714         int i;
2715         char file[128+64];
2716         t_3dvec dim;
2717         double help;
2718         t_visual *v;
2719         t_atom *atom;
2720         t_vb vb;
2721
2722         v=&(moldyn->vis);
2723         dim.x=v->dim.x;
2724         dim.y=v->dim.y;
2725         dim.z=v->dim.z;
2726         atom=moldyn->atom;
2727
2728         help=(dim.x+dim.y);
2729
2730         sprintf(file,"%s/atomic_conf_%07.f.xyz",v->fb,moldyn->time);
2731         vb.fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
2732         if(vb.fd<0) {
2733                 perror("open visual save file fd");
2734                 return -1;
2735         }
2736
2737         /* write the actual data file */
2738
2739         // povray header
2740         dprintf(vb.fd,"# [P] %d %07.f <%f,%f,%f>\n",
2741                 moldyn->count,moldyn->time,help/40.0,help/40.0,-0.8*help);
2742
2743         // atomic configuration
2744         for(i=0;i<moldyn->count;i++)
2745                 // atom type, positions, color and kinetic energy
2746                 dprintf(vb.fd,"%s %f %f %f %s %f\n",pse_name[atom[i].element],
2747                                                     atom[i].r.x,
2748                                                     atom[i].r.y,
2749                                                     atom[i].r.z,
2750                                                     pse_col[atom[i].element],
2751                                                     atom[i].ekin);
2752         
2753         // bonds between atoms
2754         process_2b_bonds(moldyn,&vb,visual_bonds_process);
2755         
2756         // boundaries
2757         if(dim.x) {
2758                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2759                         -dim.x/2,-dim.y/2,-dim.z/2,
2760                         dim.x/2,-dim.y/2,-dim.z/2);
2761                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2762                         -dim.x/2,-dim.y/2,-dim.z/2,
2763                         -dim.x/2,dim.y/2,-dim.z/2);
2764                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2765                         dim.x/2,dim.y/2,-dim.z/2,
2766                         dim.x/2,-dim.y/2,-dim.z/2);
2767                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2768                         -dim.x/2,dim.y/2,-dim.z/2,
2769                         dim.x/2,dim.y/2,-dim.z/2);
2770
2771                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2772                         -dim.x/2,-dim.y/2,dim.z/2,
2773                         dim.x/2,-dim.y/2,dim.z/2);
2774                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2775                         -dim.x/2,-dim.y/2,dim.z/2,
2776                         -dim.x/2,dim.y/2,dim.z/2);
2777                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2778                         dim.x/2,dim.y/2,dim.z/2,
2779                         dim.x/2,-dim.y/2,dim.z/2);
2780                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2781                         -dim.x/2,dim.y/2,dim.z/2,
2782                         dim.x/2,dim.y/2,dim.z/2);
2783
2784                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2785                         -dim.x/2,-dim.y/2,dim.z/2,
2786                         -dim.x/2,-dim.y/2,-dim.z/2);
2787                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2788                         -dim.x/2,dim.y/2,dim.z/2,
2789                         -dim.x/2,dim.y/2,-dim.z/2);
2790                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2791                         dim.x/2,-dim.y/2,dim.z/2,
2792                         dim.x/2,-dim.y/2,-dim.z/2);
2793                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2794                         dim.x/2,dim.y/2,dim.z/2,
2795                         dim.x/2,dim.y/2,-dim.z/2);
2796         }
2797
2798         close(vb.fd);
2799
2800         return 0;
2801 }
2802