]> hackdaworld.org Git - physik/posic.git/blob - potentials/albe_fast.c
added count column to bond_analyze code
[physik/posic.git] / potentials / albe_fast.c
1 /*
2  * test: albe_new.c
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 #ifdef PARALLEL
21 #include <omp.h>
22 #endif
23
24 #include "../moldyn.h"
25 #include "../math/math.h"
26 #include "albe.h"
27
28 /*
29  * virial calculation
30  */
31
32 #define albe_v_calc(a,f,d)      a->virial.xx+=f->x*d->x; \
33                                 a->virial.yy+=f->y*d->y; \
34                                 a->virial.zz+=f->z*d->z; \
35                                 a->virial.xy+=f->x*d->y; \
36                                 a->virial.xz+=f->x*d->z; \
37                                 a->virial.yz+=f->y*d->z
38
39 int albe_potential_force_calc(t_moldyn *moldyn) {
40
41         int i,j,k,count;
42         t_atom *itom,*jtom,*ktom;
43         t_virial *virial;
44         t_linkcell *lc;
45 #ifdef STATIC_LISTS
46         int *neighbour_i[27];
47         int p,q;
48 #elif LOWMEM_LISTS
49         int neighbour_i[27];
50         int p,q;
51 #else
52         t_list neighbour_i[27];
53         t_list neighbour_i2[27];
54         t_list *this,*that;
55 #endif
56         u8 bc_ij,bc_ik;
57         int dnlc;
58
59         // needed to work
60         t_atom *ai;
61
62         // optimized
63         t_albe_mult_params *params;
64         t_albe_exchange *exchange;
65         t_3dvec dist_ij;
66         double d_ij2;
67         double d_ij;
68         u8 brand_i;
69         double S2;
70         int kcount;
71         double zeta_ij;
72         double pre_dzeta;
73
74         // more ...
75         double Rk,Sk,Sk2;
76         t_3dvec dist_ik;
77         double d_ik2,d_ik;
78         double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
79         double f_c_ik,df_c_ik;
80
81         t_3dvec force;
82         double f_a,df_a,b,db,f_c,df_c;
83         double f_r,df_r;
84         double scale;
85         double mu,B;
86         double lambda,A;
87         double r0;
88         double S,R;
89         double energy;
90
91         double dijdik_inv,fcdg,dfcg;
92         t_3dvec dcosdrj,dcosdrk;
93         t_3dvec tmp;
94
95         // even more ...
96         double gamma_i;
97         double c_i;
98         double d_i;
99         double h_i;
100         double ci2;
101         double di2;
102         double ci2di2;
103
104         count=moldyn->count;
105         itom=moldyn->atom;
106         lc=&(moldyn->lc);
107
108         // optimized
109         params=moldyn->pot_params;
110         exchange=&(params->exchange);
111
112
113         /* reset energy */
114         moldyn->energy=0.0;
115
116         /* reset global virial */
117         memset(&(moldyn->gvir),0,sizeof(t_virial));
118
119         /* reset force, site energy and virial of every atom */
120 #ifdef PARALLEL
121         #pragma omp parallel for private(virial)
122 #endif
123         for(i=0;i<count;i++) {
124
125                 /* reset force */
126                 v3_zero(&(itom[i].f));
127
128                 /* reset virial */
129                 virial=(&(itom[i].virial));
130                 virial->xx=0.0;
131                 virial->yy=0.0;
132                 virial->zz=0.0;
133                 virial->xy=0.0;
134                 virial->xz=0.0;
135                 virial->yz=0.0;
136         
137                 /* reset site energy */
138                 itom[i].e=0.0;
139
140         }
141
142         /* get energy, force and virial of every atom */
143
144         /* first (and only) loop over atoms i */
145         for(i=0;i<count;i++) {
146
147                 if(!(itom[i].attr&ATOM_ATTR_3BP))
148                         continue;
149
150                 link_cell_neighbour_index(moldyn,
151                                           (itom[i].r.x+moldyn->dim.x/2)/lc->x,
152                                           (itom[i].r.y+moldyn->dim.y/2)/lc->y,
153                                           (itom[i].r.z+moldyn->dim.z/2)/lc->z,
154                                           neighbour_i);
155
156                 dnlc=lc->dnlc;
157
158                 /* copy the neighbour lists */
159 #ifdef STATIC_LISTS
160 #elif LOWMEM_LISTS
161 #else
162                 memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
163 #endif
164
165                 ai=&(itom[i]);
166                 brand_i=ai->brand;
167
168                 /* loop over atoms j */
169                 for(j=0;j<27;j++) {
170
171                         bc_ij=(j<dnlc)?0:1;
172 #ifdef STATIC_LISTS
173                         p=0;
174
175                         while(neighbour_i[j][p]!=-1) {
176
177                                 jtom=&(itom[neighbour_i[j][p]]);
178                                 p++;
179 #elif LOWMEM_LISTS
180                         p=neighbour_i[j];
181
182                         while(p!=-1) {
183
184                                 jtom=&(itom[p]);
185                                 p=lc->subcell->list[p];
186 #else
187                         this=&(neighbour_i[j]);
188                         list_reset_f(this);
189
190                         if(this->start==NULL)
191                                 continue;
192
193                         do {
194
195                                 jtom=this->current->data;
196 #endif
197
198                                 if(jtom==&(itom[i]))
199                                         continue;
200
201                                 if(!(jtom->attr&ATOM_ATTR_3BP))
202                                         continue;
203
204                                 /* reset 3bp run */
205                                 moldyn->run3bp=1;
206
207
208 /* j1 func here ... */
209 /* albe 3 body potential function (first ij loop) */
210
211         /* reset zeta sum */
212         zeta_ij=0.0;
213
214         /*
215          * set ij depending values
216          */
217
218         if(brand_i==jtom->brand) {
219                 S2=params->S2[brand_i];
220         }
221         else {
222                 S2=params->S2mixed;
223         }
224
225         /* dist_ij, d_ij2 */
226         v3_sub(&dist_ij,&(jtom->r),&(ai->r));
227         if(bc_ij) check_per_bound(moldyn,&dist_ij);
228         d_ij2=v3_absolute_square(&dist_ij);
229
230         /* if d_ij2 > S2 => no force & potential energy contribution */
231         if(d_ij2>S2)
232                 continue;
233
234         /* d_ij */
235         d_ij=sqrt(d_ij2);
236
237         /* reset k counter for first k loop */
238         kcount=0;
239
240                                 /* first loop over atoms k */
241                                 for(k=0;k<27;k++) {
242
243                                         bc_ik=(k<dnlc)?0:1;
244 #ifdef STATIC_LISTS
245                                         q=0;
246
247                                         while(neighbour_i[k][q]!=-1) {
248
249                                                 ktom=&(itom[neighbour_i[k][q]]);
250                                                 q++;
251 #elif LOWMEM_LISTS
252                                         q=neighbour_i[k];
253
254                                         while(q!=-1) {
255
256                                                 ktom=&(itom[q]);
257                                                 q=lc->subcell->list[q];
258 #else
259                                         that=&(neighbour_i2[k]);
260                                         list_reset_f(that);
261                                         
262                                         if(that->start==NULL)
263                                                 continue;
264
265                                         do {
266                                                 ktom=that->current->data;
267 #endif
268
269                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
270                                                         continue;
271
272                                                 if(ktom==jtom)
273                                                         continue;
274
275                                                 if(ktom==&(itom[i]))
276                                                         continue;
277
278
279 /* k1 func here ... */
280 /* albe 3 body potential function (first k loop) */
281
282         if(kcount>ALBE_MAXN) {
283                 printf("FATAL: neighbours = %d\n",kcount);
284                 printf("  -> %d %d %d\n",ai->tag,jtom->tag,ktom->tag);
285         }
286
287         /* ik constants */
288         if(brand_i==ktom->brand) {
289                 Rk=params->R[brand_i];
290                 Sk=params->S[brand_i];
291                 Sk2=params->S2[brand_i];
292                 /* albe needs i,k depending c,d,h and gamma values */
293                 gamma_i=params->gamma[brand_i];
294                 c_i=params->c[brand_i];
295                 d_i=params->d[brand_i];
296                 h_i=params->h[brand_i];
297                 ci2=params->c2[brand_i];
298                 di2=params->d2[brand_i];
299                 ci2di2=params->c2d2[brand_i];
300         }
301         else {
302                 Rk=params->Rmixed;
303                 Sk=params->Smixed;
304                 Sk2=params->S2mixed;
305                 /* albe needs i,k depending c,d,h and gamma values */
306                 gamma_i=params->gamma_m;
307                 c_i=params->c_mixed;
308                 d_i=params->d_mixed;
309                 h_i=params->h_mixed;
310                 ci2=params->c2_mixed;
311                 di2=params->d2_mixed;
312                 ci2di2=params->c2d2_m;
313         }
314
315         /* dist_ik, d_ik2 */
316         v3_sub(&dist_ik,&(ktom->r),&(ai->r));
317         if(bc_ik) check_per_bound(moldyn,&dist_ik);
318         d_ik2=v3_absolute_square(&dist_ik);
319
320         /* store data for second k loop */
321         exchange->dist_ik[kcount]=dist_ik;
322         exchange->d_ik2[kcount]=d_ik2;
323
324         /* return if not within cutoff */
325         if(d_ik2>Sk2) {
326                 kcount++;
327                 continue;
328         }
329
330         /* d_ik */
331         d_ik=sqrt(d_ik2);
332
333         /* cos theta */
334         cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
335
336         /* g_ijk 
337         h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism
338         d2_h_cos2=exchange->di2+(h_cos*h_cos);
339         frac=exchange->ci2/d2_h_cos2;
340         g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
341         dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
342         */
343
344         h_cos=h_i+cos_theta; // + in albe formalism
345         d2_h_cos2=di2+(h_cos*h_cos);
346         frac=ci2/d2_h_cos2;
347         g=gamma_i*(1.0+ci2di2-frac);
348         dg=2.0*frac*gamma_i*h_cos/d2_h_cos2; // + in albe f..
349
350         /* zeta sum += f_c_ik * g_ijk */
351         if(d_ik<=Rk) {
352                 zeta_ij+=g;
353                 f_c_ik=1.0;
354                 df_c_ik=0.0;
355         }
356         else {
357                 s_r=Sk-Rk;
358                 arg=M_PI*(d_ik-Rk)/s_r;
359                 f_c_ik=0.5+0.5*cos(arg);
360                 df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
361                 zeta_ij+=f_c_ik*g;
362         }
363
364         /* store even more data for second k loop */
365         exchange->g[kcount]=g;
366         exchange->dg[kcount]=dg;
367         exchange->d_ik[kcount]=d_ik;
368         exchange->cos_theta[kcount]=cos_theta;
369         exchange->f_c_ik[kcount]=f_c_ik;
370         exchange->df_c_ik[kcount]=df_c_ik;
371
372         /* increase k counter */
373         kcount++;
374
375 #ifdef STATIC_LISTS
376                                         }
377 #elif LOWMEM_LISTS
378                                         }
379 #else
380                                         } while(list_next_f(that)!=\
381                                                 L_NO_NEXT_ELEMENT);
382 #endif
383
384                                 }
385
386
387 /* j2 func here ... */
388
389
390         if(brand_i==jtom->brand) {
391                 S=params->S[brand_i];
392                 R=params->R[brand_i];
393                 B=params->B[brand_i];
394                 A=params->A[brand_i];
395                 r0=params->r0[brand_i];
396                 mu=params->mu[brand_i];
397                 lambda=params->lambda[brand_i];
398         }
399         else {
400                 S=params->Smixed;
401                 R=params->Rmixed;
402                 B=params->Bmixed;
403                 A=params->Amixed;
404                 r0=params->r0_mixed;
405                 mu=params->mu_m;
406                 lambda=params->lambda_m;
407         }
408
409         /* f_c, df_c */
410         if(d_ij<R) {
411                 f_c=1.0;
412                 df_c=0.0;
413         }
414         else {
415                 s_r=S-R;
416                 arg=M_PI*(d_ij-R)/s_r;
417                 f_c=0.5+0.5*cos(arg);
418                 df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
419         }
420
421         /* f_a, df_a */
422         f_a=-B*exp(-mu*(d_ij-r0));
423         df_a=mu*f_a/d_ij;
424
425         /* f_r, df_r */
426         f_r=A*exp(-lambda*(d_ij-r0));
427         df_r=lambda*f_r/d_ij;
428
429         /* b, db */
430         if(zeta_ij==0.0) {
431                 b=1.0;
432                 db=0.0;
433         }
434         else {
435                 b=1.0/sqrt(1.0+zeta_ij);
436                 db=-0.5*b/(1.0+zeta_ij);
437         }
438
439         /* force contribution for atom i */
440         scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
441         v3_scale(&force,&(dist_ij),scale);
442         v3_add(&(ai->f),&(ai->f),&force);
443
444         /* force contribution for atom j */
445         v3_scale(&force,&force,-1.0); // dri rij = - drj rij
446         v3_add(&(jtom->f),&(jtom->f),&force);
447
448         /* virial */
449         virial_calc(ai,&force,&(dist_ij));
450
451 #ifdef DEBUG
452 if(moldyn->time>DSTART&&moldyn->time<DEND) {
453         if((ai==&(moldyn->atom[DATOM]))|(jtom==&(moldyn->atom[DATOM]))) {
454                 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,jtom->tag);
455                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
456                 if(ai==&(moldyn->atom[0]))
457                         printf("  total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
458                 if(jtom==&(moldyn->atom[0]))
459                         printf("  total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
460                 printf("  energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
461                                                     f_c,b,f_a,f_r);
462                 printf("          %f %f %f\n",zeta_ij,.0,.0);
463         }
464 }
465 #endif
466
467         /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
468         pre_dzeta=0.5*f_a*f_c*db;
469
470         /* energy contribution */
471         energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
472         moldyn->energy+=energy;
473         ai->e+=energy;
474
475         /* reset k counter for second k loop */
476         kcount=0;
477                 
478
479                                 /* second loop over atoms k */
480                                 for(k=0;k<27;k++) {
481
482                                         bc_ik=(k<dnlc)?0:1;
483 #ifdef STATIC_LISTS
484                                         q=0;
485
486                                         while(neighbour_i[k][q]!=-1) {
487
488                                                 ktom=&(itom[neighbour_i[k][q]]);
489                                                 q++;
490 #elif LOWMEM_LISTS
491                                         q=neighbour_i[k];
492
493                                         while(q!=-1) {
494
495                                                 ktom=&(itom[q]);
496                                                 q=lc->subcell->list[q];
497 #else
498                                         that=&(neighbour_i2[k]);
499                                         list_reset_f(that);
500                                         
501                                         if(that->start==NULL)
502                                                 continue;
503
504                                         do {
505                                                 ktom=that->current->data;
506 #endif
507
508                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
509                                                         continue;
510
511                                                 if(ktom==jtom)
512                                                         continue;
513
514                                                 if(ktom==&(itom[i]))
515                                                         continue;
516
517
518 /* k2 func here ... */
519 /* albe 3 body potential function (second k loop) */
520
521         if(kcount>ALBE_MAXN)
522                 printf("FATAL: neighbours!\n");
523
524         /* d_ik2 */
525         d_ik2=exchange->d_ik2[kcount];
526
527         if(brand_i==ktom->brand)
528                 Sk2=params->S2[brand_i];
529         else
530                 Sk2=params->S2mixed;
531
532         /* return if d_ik > S */
533         if(d_ik2>Sk2) {
534                 kcount++;
535                 continue;
536         }
537
538         /* dist_ik, d_ik */
539         dist_ik=exchange->dist_ik[kcount];
540         d_ik=exchange->d_ik[kcount];
541
542         /* f_c_ik, df_c_ik */
543         f_c_ik=exchange->f_c_ik[kcount];
544         df_c_ik=exchange->df_c_ik[kcount];
545
546         /* g, dg, cos_theta */
547         g=exchange->g[kcount];
548         dg=exchange->dg[kcount];
549         cos_theta=exchange->cos_theta[kcount];
550
551         /* cos_theta derivatives wrt j,k */
552         dijdik_inv=1.0/(d_ij*d_ik);
553         v3_scale(&dcosdrj,&dist_ik,dijdik_inv);         // j
554         v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
555         v3_add(&dcosdrj,&dcosdrj,&tmp);
556         v3_scale(&dcosdrk,&dist_ij,dijdik_inv);         // k
557         v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
558         v3_add(&dcosdrk,&dcosdrk,&tmp);
559
560         /* f_c_ik * dg, df_c_ik * g */
561         fcdg=f_c_ik*dg;
562         dfcg=df_c_ik*g;
563
564         /* derivative wrt j */
565         v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
566
567         /* force contribution */
568         v3_add(&(jtom->f),&(jtom->f),&force);
569
570 #ifdef DEBUG
571 if(moldyn->time>DSTART&&moldyn->time<DEND) {
572         if(jtom==&(moldyn->atom[DATOM])) {
573                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
574                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
575                 printf("  total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
576                 printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
577                 printf("    d ij ik = %f %f\n",d_ij,d_ik);
578         }
579 }
580 #endif
581
582         /* virial */
583         virial_calc(ai,&force,&dist_ij);
584
585         /* force contribution to atom i */
586         v3_scale(&force,&force,-1.0);
587         v3_add(&(ai->f),&(ai->f),&force);
588
589         /* derivative wrt k */
590         v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
591         v3_scale(&tmp,&dcosdrk,fcdg);
592         v3_add(&force,&force,&tmp);
593         v3_scale(&force,&force,pre_dzeta);
594
595         /* force contribution */
596         v3_add(&(ktom->f),&(ktom->f),&force);
597
598 #ifdef DEBUG
599 if(moldyn->time>DSTART&&moldyn->time<DEND) {
600         if(ktom==&(moldyn->atom[DATOM])) {
601                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
602                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
603                 printf("  total k: %f %f %f\n",ktom->f.x,ktom->f.y,ktom->f.z);
604                 printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
605                 printf("    d ij ik = %f %f\n",d_ij,d_ik);
606         }
607 }
608 #endif
609
610         /* virial */
611         virial_calc(ai,&force,&dist_ik);
612         
613         /* force contribution to atom i */
614         v3_scale(&force,&force,-1.0);
615         v3_add(&(ai->f),&(ai->f),&force);
616
617         /* increase k counter */
618         kcount++;       
619
620
621
622 #ifdef STATIC_LISTS
623                                         }
624 #elif LOWMEM_LISTS
625                                         }
626 #else
627                                         } while(list_next_f(that)!=\
628                                                 L_NO_NEXT_ELEMENT);
629 #endif
630
631                                 }
632                                 
633 #ifdef STATIC_LISTS
634                         }
635 #elif LOWMEM_LISTS
636                         }
637 #else
638                         } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
639 #endif
640                 
641                 }
642                 
643 #ifdef DEBUG
644         //printf("\n\n");
645 #endif
646 #ifdef VDEBUG
647         printf("\n\n");
648 #endif
649
650         }
651
652 #ifdef DEBUG
653         //printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
654         if(moldyn->time>DSTART&&moldyn->time<DEND) {
655                 printf("force:\n");
656                 printf("  x: %0.40f\n",moldyn->atom[DATOM].f.x);
657                 printf("  y: %0.40f\n",moldyn->atom[DATOM].f.y);
658                 printf("  z: %0.40f\n",moldyn->atom[DATOM].f.z);
659         }
660 #endif
661
662         /* some postprocessing */
663 #ifdef PARALLEL
664         #pragma omp parallel for
665 #endif
666         for(i=0;i<count;i++) {
667                 /* calculate global virial */
668                 moldyn->gvir.xx+=itom[i].r.x*itom[i].f.x;
669                 moldyn->gvir.yy+=itom[i].r.y*itom[i].f.y;
670                 moldyn->gvir.zz+=itom[i].r.z*itom[i].f.z;
671                 moldyn->gvir.xy+=itom[i].r.y*itom[i].f.x;
672                 moldyn->gvir.xz+=itom[i].r.z*itom[i].f.x;
673                 moldyn->gvir.yz+=itom[i].r.z*itom[i].f.y;
674
675                 /* check forces regarding the given timestep */
676                 if(v3_norm(&(itom[i].f))>\
677                    0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square)
678                         printf("[moldyn] WARNING: pfc (high force: atom %d)\n",
679                                i);
680         }
681
682         return 0;
683 }
684