]> www.hackdaworld.org Git - physik/posic.git/blob - potentials/albe.c
b00e173757e0b4bb3d538b110c610cb8bb46d891
[physik/posic.git] / potentials / albe.c
1 /*
2  * albe.c - albe potential
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 <math.h>
17
18 #include "../moldyn.h"
19 #include "../math/math.h"
20 #include "albe.h"
21
22 /* create mixed terms from parameters and set them */
23 int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
24
25         t_albe_mult_params *p;
26
27         // set cutoff before parameters (actually only necessary for some pots)
28         if(moldyn->cutoff==0.0) {
29                 printf("[albe] WARNING: no cutoff!\n");
30                 return -1;
31         }
32
33         /* alloc mem for potential parameters */
34         moldyn->pot_params=malloc(sizeof(t_albe_mult_params));
35         if(moldyn->pot_params==NULL) {
36                 perror("[albe] pot params alloc");
37                 return -1;
38         }
39
40         /* these are now albe parameters */
41         p=moldyn->pot_params;
42
43         // only 1 combination by now :p
44         switch(element1) {
45                 case SI:
46                         /* type: silicon */
47                         p->S[0]=ALBE_S_SI;
48                         p->R[0]=ALBE_R_SI;
49                         p->A[0]=ALBE_A_SI;
50                         p->B[0]=ALBE_B_SI;
51                         p->r0[0]=ALBE_R0_SI;
52                         p->lambda[0]=ALBE_LAMBDA_SI;
53                         p->mu[0]=ALBE_MU_SI;
54                         p->gamma[0]=ALBE_GAMMA_SI;
55                         p->c[0]=ALBE_C_SI;
56                         p->d[0]=ALBE_D_SI;
57                         p->h[0]=ALBE_H_SI;
58                         switch(element2) {
59                                 case C:
60                                         /* type: carbon */
61                                         p->S[1]=ALBE_S_C;
62                                         p->R[1]=ALBE_R_C;
63                                         p->A[1]=ALBE_A_C;
64                                         p->B[1]=ALBE_B_C;
65                                         p->r0[1]=ALBE_R0_C;
66                                         p->lambda[1]=ALBE_LAMBDA_C;
67                                         p->mu[1]=ALBE_MU_C;
68                                         p->gamma[1]=ALBE_GAMMA_C;
69                                         p->c[1]=ALBE_C_C;
70                                         p->d[1]=ALBE_D_C;
71                                         p->h[1]=ALBE_H_C;
72                                         /* mixed type: silicon carbide */
73                                         p->Smixed=ALBE_S_SIC;
74                                         p->Rmixed=ALBE_R_SIC;
75                                         p->Amixed=ALBE_A_SIC;
76                                         p->Bmixed=ALBE_B_SIC;
77                                         p->r0_mixed=ALBE_R0_SIC;
78                                         p->lambda_m=ALBE_LAMBDA_SIC;
79                                         p->mu_m=ALBE_MU_SIC;
80                                         p->gamma_m=ALBE_GAMMA_SIC;
81                                         p->c_mixed=ALBE_C_SIC;
82                                         p->d_mixed=ALBE_D_SIC;
83                                         p->h_mixed=ALBE_H_SIC;
84                                         break;
85                                 default:
86                                         printf("[albe] WARNING: element2\n");
87                                         return -1;
88                         }
89                         break;
90                 default:
91                         printf("[albe] WARNING: element1\n");
92                         return -1;
93         }
94
95         printf("[albe] parameter completion\n");
96         p->S2[0]=p->S[0]*p->S[0];
97         p->S2[1]=p->S[1]*p->S[1];
98         p->S2mixed=p->Smixed*p->Smixed;
99         p->c2[0]=p->c[0]*p->c[0];
100         p->c2[1]=p->c[1]*p->c[1];
101         p->c2_mixed=p->c_mixed*p->c_mixed;
102         p->d2[0]=p->d[0]*p->d[0];
103         p->d2[1]=p->d[1]*p->d[1];
104         p->d2_mixed=p->d_mixed*p->d_mixed;
105         p->c2d2[0]=p->c2[0]/p->d2[0];
106         p->c2d2[1]=p->c2[1]/p->d2[1];
107         p->c2d2_m=p->c2_mixed/p->d2_mixed;
108
109         printf("[albe] mult parameter info:\n");
110         printf("  S (A)  | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed);
111         printf("  R (A)  | %f | %f | %f\n",p->R[0],p->R[1],p->Rmixed);
112         printf("  A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV);
113         printf("  B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV);
114         printf("  lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1],
115                                           p->lambda_m);
116         printf("  mu     | %f | %f | %f\n",p->mu[0],p->mu[1],p->mu_m);
117         printf("  gamma  | %f | %f\n",p->gamma[0],p->gamma[1]);
118         printf("  c      | %f | %f\n",p->c[0],p->c[1]);
119         printf("  d      | %f | %f\n",p->d[0],p->d[1]);
120         printf("  h      | %f | %f\n",p->h[0],p->h[1]);
121
122         return 0;
123 }
124
125 /* albe 3 body potential function (first ij loop) */
126 int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
127
128         t_albe_mult_params *params;
129         t_albe_exchange *exchange;
130         unsigned char brand;
131         double S2;
132         t_3dvec dist_ij;
133         double d_ij2,d_ij;
134
135         params=moldyn->pot_params;
136         exchange=&(params->exchange);
137
138         /* reset zeta sum */
139         exchange->zeta_ij=0.0;
140
141         /*
142          * set ij depending values
143          */
144
145         brand=ai->brand;
146         if(brand==aj->brand) {
147                 S2=params->S2[brand];
148         }
149         else {
150                 S2=params->S2mixed;
151         }
152
153         /* dist_ij, d_ij2 */
154         v3_sub(&dist_ij,&(aj->r),&(ai->r));
155         if(bc) check_per_bound(moldyn,&dist_ij);
156         d_ij2=v3_absolute_square(&dist_ij);
157
158         /* if d_ij2 > S2 => no force & potential energy contribution */
159         if(d_ij2>S2) {
160                 moldyn->run3bp=0;
161                 return 0;
162         }
163
164         /* d_ij */
165         d_ij=sqrt(d_ij2);
166
167         /* store values */
168         exchange->dist_ij=dist_ij;
169         exchange->d_ij2=d_ij2;
170         exchange->d_ij=d_ij;
171
172         /* reset k counter for first k loop */
173         exchange->kcount=0;
174                 
175         return 0;
176 }
177
178 /* albe 3 body potential function (first k loop) */
179 #ifdef PTHREADS
180 void *albe_mult_3bp_k1(void *ptr) {
181 #else
182 int albe_mult_3bp_k1(t_moldyn *moldyn,
183                      t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
184 #endif
185
186         t_albe_mult_params *params;
187         t_albe_exchange *exchange;
188         unsigned char brand;
189         double R,S,S2;
190         t_3dvec dist_ij,dist_ik;
191         double d_ik2,d_ik,d_ij;
192         double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
193         double f_c_ik,df_c_ik;
194         int kcount;
195 #ifdef PTHREADS
196         t_kdata *kdata;
197         t_moldyn *moldyn;
198         t_atom *ai,*aj,*ak;
199         u8 bc;
200
201         kdata=ptr;
202         moldyn=kdata->moldyn;
203         ai=kdata->ai;
204         aj=kdata->aj;
205         ak=kdata->ak;
206         bc=kdata->bc;
207 #endif
208
209         params=moldyn->pot_params;
210         exchange=&(params->exchange);
211         kcount=exchange->kcount;
212
213         if(kcount>ALBE_MAXN) {
214                 printf("FATAL: neighbours = %d\n",kcount);
215                 printf("  -> %d %d %d\n",ai->tag,aj->tag,ak->tag);
216         }
217
218         /* ik constants */
219         brand=ai->brand;
220         if(brand==ak->brand) {
221                 R=params->R[brand];
222                 S=params->S[brand];
223                 S2=params->S2[brand];
224                 /* albe needs i,k depending c,d,h and gamma values */
225                 exchange->gamma_i=&(params->gamma[brand]);
226                 exchange->c_i=&(params->c[brand]);
227                 exchange->d_i=&(params->d[brand]);
228                 exchange->h_i=&(params->h[brand]);
229         }
230         else {
231                 R=params->Rmixed;
232                 S=params->Smixed;
233                 S2=params->S2mixed;
234                 /* albe needs i,k depending c,d,h and gamma values */
235                 exchange->gamma_i=&(params->gamma_m);
236                 exchange->c_i=&(params->c_mixed);
237                 exchange->d_i=&(params->d_mixed);
238                 exchange->h_i=&(params->h_mixed);
239         }
240         exchange->ci2=*(exchange->c_i)**(exchange->c_i);
241         exchange->di2=*(exchange->d_i)**(exchange->d_i);
242         exchange->ci2di2=exchange->ci2/exchange->di2;
243
244         /* dist_ik, d_ik2 */
245         v3_sub(&dist_ik,&(ak->r),&(ai->r));
246         if(bc) check_per_bound(moldyn,&dist_ik);
247         d_ik2=v3_absolute_square(&dist_ik);
248
249         /* store data for second k loop */
250         exchange->dist_ik[kcount]=dist_ik;
251         exchange->d_ik2[kcount]=d_ik2;
252
253         /* return if not within cutoff */
254         if(d_ik2>S2) {
255                 exchange->kcount++;
256                 return 0;
257         }
258
259         /* d_ik */
260         d_ik=sqrt(d_ik2);
261
262         /* dist_ij, d_ij */
263         dist_ij=exchange->dist_ij;
264         d_ij=exchange->d_ij;
265
266         /* cos theta */
267         cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
268
269         /* g_ijk */
270         h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism
271         d2_h_cos2=exchange->di2+(h_cos*h_cos);
272         frac=exchange->ci2/d2_h_cos2;
273         g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
274         dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
275
276         /* zeta sum += f_c_ik * g_ijk */
277         if(d_ik<=R) {
278                 exchange->zeta_ij+=g;
279                 f_c_ik=1.0;
280                 df_c_ik=0.0;
281         }
282         else {
283                 s_r=S-R;
284                 arg=M_PI*(d_ik-R)/s_r;
285                 f_c_ik=0.5+0.5*cos(arg);
286                 df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
287                 exchange->zeta_ij+=f_c_ik*g;
288         }
289
290         /* store even more data for second k loop */
291         exchange->g[kcount]=g;
292         exchange->dg[kcount]=dg;
293         exchange->d_ik[kcount]=d_ik;
294         exchange->cos_theta[kcount]=cos_theta;
295         exchange->f_c_ik[kcount]=f_c_ik;
296         exchange->df_c_ik[kcount]=df_c_ik;
297
298         /* increase k counter */
299         exchange->kcount++;
300
301 #ifdef PTHREADS
302 }
303 #else
304         return 0;
305 }
306 #endif
307
308 int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
309
310         t_albe_mult_params *params;
311         t_albe_exchange *exchange;
312         t_3dvec force;
313         double f_a,df_a,b,db,f_c,df_c;
314         double f_r,df_r;
315         double scale;
316         double mu,B;
317         double lambda,A;
318         double d_ij,r0;
319         unsigned char brand;
320         double S,R,s_r,arg;
321         double energy;
322
323         params=moldyn->pot_params;
324         exchange=&(params->exchange);
325
326         brand=aj->brand;
327         if(brand==ai->brand) {
328                 S=params->S[brand];
329                 R=params->R[brand];
330                 B=params->B[brand];
331                 A=params->A[brand];
332                 r0=params->r0[brand];
333                 mu=params->mu[brand];
334                 lambda=params->lambda[brand];
335         }
336         else {
337                 S=params->Smixed;
338                 R=params->Rmixed;
339                 B=params->Bmixed;
340                 A=params->Amixed;
341                 r0=params->r0_mixed;
342                 mu=params->mu_m;
343                 lambda=params->lambda_m;
344         }
345
346         d_ij=exchange->d_ij;
347
348         /* f_c, df_c */
349         if(d_ij<R) {
350                 f_c=1.0;
351                 df_c=0.0;
352         }
353         else {
354                 s_r=S-R;
355                 arg=M_PI*(d_ij-R)/s_r;
356                 f_c=0.5+0.5*cos(arg);
357                 df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
358         }
359
360         /* f_a, df_a */
361         f_a=-B*exp(-mu*(d_ij-r0));
362         df_a=mu*f_a/d_ij;
363
364         /* f_r, df_r */
365         f_r=A*exp(-lambda*(d_ij-r0));
366         df_r=lambda*f_r/d_ij;
367
368         /* b, db */
369         if(exchange->zeta_ij==0.0) {
370                 b=1.0;
371                 db=0.0;
372         }
373         else {
374                 b=1.0/sqrt(1.0+exchange->zeta_ij);
375                 db=-0.5*b/(1.0+exchange->zeta_ij);
376         }
377
378         /* force contribution for atom i */
379         scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
380         v3_scale(&force,&(exchange->dist_ij),scale);
381         v3_add(&(ai->f),&(ai->f),&force);
382
383         /* force contribution for atom j */
384         v3_scale(&force,&force,-1.0); // dri rij = - drj rij
385         v3_add(&(aj->f),&(aj->f),&force);
386
387         /* virial */
388         virial_calc(ai,&force,&(exchange->dist_ij));
389
390 #ifdef DEBUG
391 if(moldyn->time>DSTART&&moldyn->time<DEND) {
392         if((ai==&(moldyn->atom[DATOM]))|(aj==&(moldyn->atom[DATOM]))) {
393                 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,aj->tag);
394                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
395                 if(ai==&(moldyn->atom[0]))
396                         printf("  total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
397                 if(aj==&(moldyn->atom[0]))
398                         printf("  total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
399                 printf("  energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
400                                                     f_c,b,f_a,f_r);
401                 printf("          %f %f %f\n",exchange->zeta_ij,.0,.0);
402         }
403 }
404 #endif
405
406         /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
407         exchange->pre_dzeta=0.5*f_a*f_c*db;
408
409         /* energy contribution */
410         energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
411         moldyn->energy+=energy;
412         ai->e+=energy;
413
414         /* reset k counter for second k loop */
415         exchange->kcount=0;
416                 
417         return 0;
418 }
419
420 /* albe 3 body potential function (second k loop) */
421 int albe_mult_3bp_k2(t_moldyn *moldyn,
422                      t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
423
424         t_albe_mult_params *params;
425         t_albe_exchange *exchange;
426         int kcount;
427         t_3dvec dist_ik,dist_ij;
428         double d_ik2,d_ik,d_ij2,d_ij;
429         unsigned char brand;
430         double S2;
431         double g,dg,cos_theta;
432         double pre_dzeta;
433         double f_c_ik,df_c_ik;
434         double dijdik_inv,fcdg,dfcg;
435         t_3dvec dcosdrj,dcosdrk;
436         t_3dvec force,tmp;
437
438         params=moldyn->pot_params;
439         exchange=&(params->exchange);
440         kcount=exchange->kcount;
441
442         if(kcount>ALBE_MAXN)
443                 printf("FATAL: neighbours!\n");
444
445         /* d_ik2 */
446         d_ik2=exchange->d_ik2[kcount];
447
448         brand=ak->brand;
449         if(brand==ai->brand)
450                 S2=params->S2[brand];
451         else
452                 S2=params->S2mixed;
453
454         /* return if d_ik > S */
455         if(d_ik2>S2) {
456                 exchange->kcount++;
457                 return 0;
458         }
459
460         /* prefactor dzeta */
461         pre_dzeta=exchange->pre_dzeta;
462
463         /* dist_ik, d_ik */
464         dist_ik=exchange->dist_ik[kcount];
465         d_ik=exchange->d_ik[kcount];
466
467         /* f_c_ik, df_c_ik */
468         f_c_ik=exchange->f_c_ik[kcount];
469         df_c_ik=exchange->df_c_ik[kcount];
470
471         /* dist_ij, d_ij2, d_ij */
472         dist_ij=exchange->dist_ij;
473         d_ij2=exchange->d_ij2;
474         d_ij=exchange->d_ij;
475
476         /* g, dg, cos_theta */
477         g=exchange->g[kcount];
478         dg=exchange->dg[kcount];
479         cos_theta=exchange->cos_theta[kcount];
480
481         /* cos_theta derivatives wrt j,k */
482         dijdik_inv=1.0/(d_ij*d_ik);
483         v3_scale(&dcosdrj,&dist_ik,dijdik_inv);         // j
484         v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
485         v3_add(&dcosdrj,&dcosdrj,&tmp);
486         v3_scale(&dcosdrk,&dist_ij,dijdik_inv);         // k
487         v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
488         v3_add(&dcosdrk,&dcosdrk,&tmp);
489
490         /* f_c_ik * dg, df_c_ik * g */
491         fcdg=f_c_ik*dg;
492         dfcg=df_c_ik*g;
493
494         /* derivative wrt j */
495         v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
496
497         /* force contribution */
498         v3_add(&(aj->f),&(aj->f),&force);
499
500 #ifdef DEBUG
501 if(moldyn->time>DSTART&&moldyn->time<DEND) {
502         if(aj==&(moldyn->atom[DATOM])) {
503                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
504                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
505                 printf("  total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
506                 printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
507                 printf("    d ij ik = %f %f\n",d_ij,d_ik);
508         }
509 }
510 #endif
511
512         /* virial */
513         virial_calc(ai,&force,&dist_ij);
514
515         /* force contribution to atom i */
516         v3_scale(&force,&force,-1.0);
517         v3_add(&(ai->f),&(ai->f),&force);
518
519         /* derivative wrt k */
520         v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
521         v3_scale(&tmp,&dcosdrk,fcdg);
522         v3_add(&force,&force,&tmp);
523         v3_scale(&force,&force,pre_dzeta);
524
525         /* force contribution */
526         v3_add(&(ak->f),&(ak->f),&force);
527
528 #ifdef DEBUG
529 if(moldyn->time>DSTART&&moldyn->time<DEND) {
530         if(ak==&(moldyn->atom[DATOM])) {
531                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
532                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
533                 printf("  total k: %f %f %f\n",ak->f.x,ak->f.y,ak->f.z);
534                 printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
535                 printf("    d ij ik = %f %f\n",d_ij,d_ik);
536         }
537 }
538 #endif
539
540         /* virial */
541         virial_calc(ai,&force,&dist_ik);
542         
543         /* force contribution to atom i */
544         v3_scale(&force,&force,-1.0);
545         v3_add(&(ai->f),&(ai->f),&force);
546
547         /* increase k counter */
548         exchange->kcount++;     
549
550         return 0;
551
552 }
553
554 int albe_mult_check_2b_bond(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,u8 bc) {
555
556         t_albe_mult_params *params;
557         t_3dvec dist;
558         double d;
559         u8 brand;
560
561         v3_sub(&dist,&(jtom->r),&(itom->r));
562         if(bc) check_per_bound(moldyn,&dist);
563         d=v3_absolute_square(&dist);
564
565         params=moldyn->pot_params;
566         brand=itom->brand;
567
568         if(brand==jtom->brand) {
569                 if(d<=params->S2[brand])
570                         return TRUE;
571         }
572         else {
573                 if(d<=params->S2mixed)
574                         return TRUE;
575         }
576
577         return FALSE;
578 }