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