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