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