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