]> www.hackdaworld.org Git - physik/posic.git/blob - potentials/albe.c
new albe potential + new force calc routine (old potentials need to be
[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->c2[0]=p->c[0]*p->c[0];
57                         p->d[0]=ALBE_D_SI;
58                         p->d2[0]=p->d[0]*p->d[0];
59                         p->c2d2[0]=p->c2[0]/p->d2[0];
60                         p->h[0]=ALBE_H_SI;
61                         switch(element2) {
62                                 case C:
63                                         /* type: carbon */
64                                         p->S[1]=ALBE_S_C;
65                                         p->R[1]=ALBE_R_C;
66                                         p->A[1]=ALBE_A_C;
67                                         p->B[1]=ALBE_B_C;
68                                         p->r0[1]=ALBE_R0_C;
69                                         p->lambda[1]=ALBE_LAMBDA_C;
70                                         p->mu[1]=ALBE_MU_C;
71                                         p->gamma[1]=ALBE_GAMMA_C;
72                                         p->c[1]=ALBE_C_C;
73                                         p->c2[1]=p->c[1]*p->c[1];
74                                         p->d[1]=ALBE_D_C;
75                                         p->d2[1]=p->d[1]*p->d[1];
76                                         p->c2d2[1]=p->c2[1]/p->d2[1];
77                                         p->h[1]=ALBE_H_C;
78                                         /* mixed type: silicon carbide */
79                                         p->Smixed=ALBE_S_SIC;
80                                         p->Rmixed=ALBE_R_SIC;
81                                         p->Amixed=ALBE_A_SIC;
82                                         p->Bmixed=ALBE_B_SIC;
83                                         p->r0_mixed=ALBE_R0_SIC;
84                                         p->lambda_m=ALBE_LAMBDA_SIC;
85                                         p->mu_m=ALBE_MU_SIC;
86                                         p->gamma_m=ALBE_GAMMA_SIC;
87                                         p->c_mixed=ALBE_C_SIC;
88                                         p->c2_mixed=p->c_mixed*p->c_mixed;
89                                         p->d_mixed=ALBE_D_SIC;
90                                         p->d2_mixed=p->d_mixed*p->d_mixed;
91                                         p->c2d2_m=p->c2_mixed/p->d2_mixed;
92                                         p->h_mixed=ALBE_H_SIC;
93                                         break;
94                                 default:
95                                         printf("[albe] WARNING: element2\n");
96                                         return -1;
97                         }
98                         break;
99                 default:
100                         printf("[albe] WARNING: element1\n");
101                         return -1;
102         }
103
104         printf("[albe] parameter completion\n");
105         p->S2[0]=p->S[0]*p->S[0];
106         p->S2[1]=p->S[1]*p->S[1];
107         p->S2mixed=p->Smixed*p->Smixed;
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 | %f\n",p->gamma[0],p->gamma[1],p->gamma_m);
118         printf("  c      | %f | %f | %f\n",p->c[0],p->c[1],p->c_mixed);
119         printf("  d      | %f | %f | %f\n",p->d[0],p->d[1],p->d_mixed);
120         printf("  c2     | %f | %f | %f\n",p->c2[0],p->c2[1],p->c2_mixed);
121         printf("  d2     | %f | %f | %f\n",p->d2[0],p->d2[1],p->d2_mixed);
122         printf("  c2d2   | %f | %f | %f\n",p->c2d2[0],p->c2d2[1],p->c2d2_m);
123         printf("  h      | %f | %f | %f\n",p->h[0],p->h[1],p->h_mixed);
124
125         return 0;
126 }
127
128 /* first i loop */
129 int albe_mult_i0(t_moldyn *moldyn,t_atom *ai) {
130
131         t_albe_mult_params *params;
132         t_albe_exchange *exchange;
133         
134         int i;
135
136         params=moldyn->pot_params;
137         exchange=&(params->exchange);
138
139         /* zero exchange values */
140         memset(exchange->zeta,0,ALBE_MAXN*sizeof(double));
141         for(i=0;i<ALBE_MAXN;i++)
142                 memset(exchange->dzeta[i],0,ALBE_MAXN*sizeof(t_3dvec));
143         exchange->jcnt=0;
144         exchange->j2cnt=0;
145
146         return 0;
147 }
148
149 /* first j loop within first i loop */
150 int albe_mult_i0_j0(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
151
152         t_albe_mult_params *params;
153         t_albe_exchange *exchange;
154
155         double S2,S,R,d2,d,s_r,arg;
156         t_3dvec dist;
157         int j;
158         u8 brand;
159         
160         params=moldyn->pot_params;
161         exchange=&(params->exchange);
162
163         /* get j counter */
164         j=exchange->jcnt;
165
166         /* set ij depending values */
167         brand=ai->brand;
168         if(brand==aj->brand) {
169                 S2=params->S2[brand];
170         }
171         else {
172                 S2=params->S2mixed;
173         }
174
175         /* dist_ij, d_ij2 */
176         v3_sub(&dist,&(aj->r),&(ai->r));
177         if(bc) check_per_bound(moldyn,&dist);
178         exchange->dist[j]=dist;
179         d2=v3_absolute_square(&dist);
180         exchange->d2[j]=d2;
181
182         /* if d_ij2 > S2 => no force & potential energy contribution */
183         if(d2>S2) {
184                 moldyn->run3bp=0;
185                 exchange->skip[j]=1;
186                 exchange->jcnt+=1;
187                 return 0;
188         }
189         exchange->skip[j]=0;
190
191         /* more ij depending values */
192         if(brand==aj->brand) {
193                 R=params->R[brand];
194                 S=params->S[brand];
195                 /* albe needs i,(j/k) depending c,d,h and gamma values */
196                 exchange->gamma_[j]=&(params->gamma[brand]);
197                 exchange->c_[j]=&(params->c[brand]);
198                 exchange->d_[j]=&(params->d[brand]);
199                 exchange->h_[j]=&(params->h[brand]);
200                 exchange->c2_[j]=&(params->c2[brand]);
201                 exchange->d2_[j]=&(params->d2[brand]);
202                 exchange->c2d2_[j]=&(params->c2d2[brand]);
203         }
204         else {
205                 R=params->Rmixed;
206                 S=params->Smixed;
207                 /* albe needs i,(j/k) depending c,d,h and gamma values */
208                 exchange->gamma_[j]=&(params->gamma_m);
209                 exchange->c_[j]=&(params->c_mixed);
210                 exchange->d_[j]=&(params->d_mixed);
211                 exchange->h_[j]=&(params->h_mixed);
212                 exchange->c2_[j]=&(params->c2_mixed);
213                 exchange->d2_[j]=&(params->d2_mixed);
214                 exchange->c2d2_[j]=&(params->c2d2_m);
215         }
216
217         /* d_ij */
218         d=sqrt(d2);
219         exchange->d[j]=d;
220         
221         /* f_c, df_c */
222         if(d<R) {
223                 exchange->f_c[j]=1.0;
224                 exchange->df_c[j]=0.0;
225         }
226         else {
227                 s_r=S-R;
228                 arg=M_PI*(d-R)/s_r;
229                 exchange->f_c[j]=0.5+0.5*cos(arg);
230                 exchange->df_c[j]=0.5*sin(arg)*(M_PI/(s_r*d));
231         }
232
233         /* reset k counter */
234         exchange->kcnt=0;
235
236         return 0;
237 }
238
239 /* first k loop within first j loop within first i loop */
240 int albe_mult_i0_j0_k0(t_moldyn *moldyn,
241                        t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
242
243         t_albe_mult_params *params;
244         t_albe_exchange *exchange;
245
246         int j,k;
247         t_3dvec distj,distk;
248         double dj,dk,djdk_inv,cos_theta;
249         double gj,dgj,h_cos_j,d2_h_cos2_j,frac_j;
250         double gk,dgk,h_cos_k,d2_h_cos2_k,frac_k;
251         t_3dvec dcosdrj,dcosdrk,tmp;
252         t_3dvec *dzjj,*dzkk,*dzjk,*dzkj;
253
254         params=moldyn->pot_params;
255         exchange=&(params->exchange);
256
257         if(aj==ak) {
258                 exchange->kcnt+=1;
259                 return 0;
260         }
261
262         /* k<j & check whether to run k */
263         j=exchange->jcnt;
264         k=exchange->kcnt;
265         if(k>=ALBE_MAXN) {
266                 printf("FATAL: too many neighbours! (%d)\n",k);
267                 printf("  atom i:%d | j:%d | k:%d\n",ai->tag,aj->tag,ak->tag);
268         }
269         if((k>=j)|(exchange->skip[k])) {
270                 exchange->kcnt+=1;
271                 return 0;
272         }
273
274         /* distances */
275         distj=exchange->dist[j];
276         distk=exchange->dist[k];
277         dj=exchange->d[j];
278         dk=exchange->d[k];
279         djdk_inv=1.0/(dj*dk);
280
281         /* cos theta */
282         cos_theta=v3_scalar_product(&distj,&distk)*djdk_inv;
283
284         /* g(cos(theta)) ij and ik values */
285         h_cos_j=*(exchange->h_[j])+cos_theta; // + in albe formalism
286         d2_h_cos2_j=*exchange->d2_[j]+(h_cos_j*h_cos_j);
287         frac_j=*exchange->c2_[j]/d2_h_cos2_j;
288         gj=1.0+*exchange->c2d2_[j]-frac_j;
289         gj*=*(exchange->gamma_[j]);
290         dgj=*(exchange->gamma_[j])*2.0*frac_j*h_cos_j/d2_h_cos2_j; // + in albe 
291         if(ak->brand==aj->brand) {
292                 gk=gj;
293                 dgk=dgj;
294         }
295         else {
296                 h_cos_k=*(exchange->h_[k])+cos_theta;
297                 d2_h_cos2_k=*exchange->d2_[k]+(h_cos_k*h_cos_k);
298                 frac_k=*exchange->c2_[k]/d2_h_cos2_k;
299                 gk=1.0+*exchange->c2d2_[k]-frac_k;
300                 gk*=*(exchange->gamma_[k]);
301                 dgk=*(exchange->gamma_[k])*2.0*frac_k*h_cos_k/d2_h_cos2_k;
302         }
303
304         /* zeta - for albe: ik depending g function */
305 //if(ai->tag==0) {
306 //      printf("------> %.15f %.15f\n",dj,dk);
307 //      printf("------> %.15f %.15f\n",dj,dk);
308 //}
309
310         exchange->zeta[j]+=(exchange->f_c[k]*gk);
311         exchange->zeta[k]+=(exchange->f_c[j]*gj);
312
313         /* cos theta derivatives */
314         v3_scale(&dcosdrj,&distk,djdk_inv);             // j
315         v3_scale(&tmp,&distj,-cos_theta/exchange->d2[j]);
316         v3_add(&dcosdrj,&dcosdrj,&tmp);
317         v3_scale(&dcosdrk,&distj,djdk_inv);             // k
318         v3_scale(&tmp,&distk,-cos_theta/exchange->d2[k]);
319         v3_add(&dcosdrk,&dcosdrk,&tmp);
320
321         /* zeta derivatives */
322         dzjj=&(exchange->dzeta[j][j]);
323         dzkk=&(exchange->dzeta[k][k]);
324         dzjk=&(exchange->dzeta[j][k]);
325         dzkj=&(exchange->dzeta[k][j]);
326         v3_scale(&tmp,&dcosdrj,exchange->f_c[k]*dgk);
327         v3_add(dzjj,dzjj,&tmp);                         // j j
328         v3_scale(&tmp,&dcosdrk,exchange->f_c[j]*dgj);
329         v3_add(dzkk,dzkk,&tmp);                         // k k
330         v3_scale(&tmp,&distk,-exchange->df_c[k]*gk);    // dk rik = - di rik
331         v3_add(dzjk,dzjk,&tmp);
332         v3_scale(&tmp,&dcosdrk,exchange->f_c[k]*dgk);
333         v3_add(dzjk,dzjk,&tmp);                         // j k
334         v3_scale(&tmp,&distj,-exchange->df_c[j]*gj);    // dj rij = - di rij
335         v3_add(dzkj,dzkj,&tmp);
336         v3_scale(&tmp,&dcosdrj,exchange->f_c[j]*dgj);
337         v3_add(dzkj,dzkj,&tmp);                         // k j
338
339         /* increase k counter */
340         exchange->kcnt+=1;
341                 
342         return 0;
343 }
344
345 /* first j loop within first i loop */
346 int albe_mult_i0_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
347
348         t_albe_mult_params *params;
349         t_albe_exchange *exchange;
350
351         params=moldyn->pot_params;
352         exchange=&(params->exchange);
353
354         /* increase j counter */
355         exchange->jcnt+=1;
356
357         return 0;
358 }
359
360 /* second j loop within first i loop */
361 int albe_mult_i0_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
362
363         t_albe_mult_params *params;
364         t_albe_exchange *exchange;
365
366         int j;
367         double d,f_a,df_a,f_r,df_r,f_c,df_c,b,db;
368         double A,B,mu,lambda,r0;
369         double energy;
370         t_3dvec *dist,force;
371         double scale;
372         u8 brand;
373
374         params=moldyn->pot_params;
375         exchange=&(params->exchange);
376
377         /* get j counter */
378         j=exchange->j2cnt;
379
380         /* skip if j not within cutoff */
381         if(exchange->skip[j]) {
382                 moldyn->run3bp=0;
383                 exchange->j2cnt+=1;
384                 return 0;
385         }
386
387         /* distance */
388         d=exchange->d[j];
389         dist=&(exchange->dist[j]);
390         f_c=exchange->f_c[j];
391         df_c=exchange->df_c[j];
392
393         /* determine parameters to calculate fa, dfa, fr, dfr */
394         brand=aj->brand;
395         if(brand==ai->brand) {
396                 B=params->B[brand];
397                 A=params->A[brand];
398                 r0=params->r0[brand];
399                 mu=params->mu[brand];
400                 lambda=params->lambda[brand];
401         }
402         else {
403                 B=params->Bmixed;
404                 A=params->Amixed;
405                 r0=params->r0_mixed;
406                 mu=params->mu_m;
407                 lambda=params->lambda_m;
408         }
409
410         /* f_a, df_a */
411         f_a=-B*exp(-mu*(d-r0));
412         df_a=mu*f_a/d;
413
414         /* f_r, df_r */
415         f_r=A*exp(-lambda*(d-r0));
416         df_r=lambda*f_r/d;
417
418         /* b, db */
419         b=1.0/sqrt(1.0+exchange->zeta[j]);
420         db=-0.5*b/(1.0+exchange->zeta[j]);
421
422         /* energy contribution */
423         energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
424         moldyn->energy+=energy;
425         ai->e+=energy;
426
427         /* force contribution for atom i due to ij bond */
428         scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
429         v3_scale(&force,dist,scale);
430         v3_add(&(ai->f),&(ai->f),&force);
431
432 #ifdef NDEBUG
433 if(ai->tag==0) {
434 printf("force: %.15f %.15f %.15f | %d %d (ij) %.15f\n",force.x,force.y,force.z,ai->tag,aj->tag,exchange->zeta[j]);
435 printf("    t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z);
436 }
437 #endif
438
439         /* force contribution for atom j due to ij bond */
440         v3_scale(&force,&force,-1.0); // dri rij = - drj rij
441         v3_add(&(aj->f),&(aj->f),&force);
442
443 #ifdef NDEBUG
444 if(aj->tag==0) {
445 printf("force: %.15f %.15f %.15f | %d %d (ji) %.15f\n",force.x,force.y,force.z,aj->tag,ai->tag,exchange->zeta[j]);
446 printf("    t: %.15f %.15f %.15f\n",aj->f.x,aj->f.y,aj->f.z);
447 }
448 #endif
449
450         /* virial */
451         virial_calc(aj,&force,dist);
452
453         /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
454         exchange->pre_dzeta=0.5*f_a*f_c*db;
455
456         /* force contribution (drj derivative) */
457         v3_scale(&force,&(exchange->dzeta[j][j]),exchange->pre_dzeta);
458         v3_add(&(aj->f),&(aj->f),&force);
459
460 #ifdef NDEBUG
461 if(aj->tag==0) {
462 printf("force: %.15f %.15f %.15f | %d %d (j der)\n",force.x,force.y,force.z,aj->tag,ai->tag);
463 printf("    t: %.15f %.15f %.15f\n",aj->f.x,aj->f.y,aj->f.z);
464 }
465 #endif
466
467         v3_scale(&force,&force,-1.0);
468         v3_add(&(ai->f),&(ai->f),&force);
469
470 #ifdef NDEBUG
471 if(ai->tag==0) {
472 printf("force: %.15f %.15f %.15f | %d %d (i contr j der)\n",force.x,force.y,force.z,ai->tag,aj->tag);
473 printf("    t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z);
474 }
475 #endif
476
477         /* virial */
478         virial_calc(ai,&force,dist);
479
480         /* reset k counter for second k loop */
481         exchange->kcnt=0;
482                 
483         return 0;
484 }
485
486 /* second k loop within second j loop within first i loop */
487 int albe_mult_i0_j2_k0(t_moldyn *moldyn,
488                        t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
489
490         t_albe_mult_params *params;
491         t_albe_exchange *exchange;
492
493         int j,k;
494         t_3dvec force;
495
496         params=moldyn->pot_params;
497         exchange=&(params->exchange);
498
499         if(aj==ak) {
500                 exchange->kcnt+=1;
501                 return 0;
502         }
503
504         /* k!=j & check whether to run k */
505         k=exchange->kcnt;
506         j=exchange->j2cnt;
507         if((k==j)|(exchange->skip[k])) {
508                 exchange->kcnt+=1;
509                 return 0;
510         }
511         
512         /* force contribution (drk derivative) */
513         v3_scale(&force,&(exchange->dzeta[j][k]),exchange->pre_dzeta);
514         v3_add(&(ak->f),&(ak->f),&force);
515
516 #ifdef NDEBUG
517 if(ak->tag==0) {
518 printf("force: %.15f %.15f %.15f | %d %d %d (k der)\n",force.x,force.y,force.z,ai->tag,aj->tag,ak->tag);
519 printf("    t: %.15f %.15f %.15f\n",ak->f.x,ak->f.y,ak->f.z);
520 }
521 #endif
522
523         v3_scale(&force,&force,-1.0);
524         v3_add(&(ai->f),&(ai->f),&force);
525
526 #ifdef NDEBUG
527 if(ai->tag==0) {
528 printf("force: %.15f %.15f %.15f | %d %d %d -- %d(i contr k der)\n",force.x,force.y,force.z,ai->tag,aj->tag,ak->tag,k);
529 printf("    t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z);
530 printf("  ## %f\n",exchange->d[k]);
531 }
532 #endif
533
534         /* virial */
535         virial_calc(ai,&force,&(exchange->dist[k]));
536
537         /* increase k counter */
538         exchange->kcnt+=1;
539
540         return 0;
541 }
542
543 int albe_mult_i0_j3(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
544
545         t_albe_mult_params *params;
546         t_albe_exchange *exchange;
547
548         params=moldyn->pot_params;
549         exchange=&(params->exchange);
550
551         /* increase j counter */
552         exchange->j2cnt+=1;
553
554         return 0;
555 }
556
557 int albe_mult_check_2b_bond(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,u8 bc) {
558
559         t_albe_mult_params *params;
560         t_3dvec dist;
561         double d;
562         u8 brand;
563
564         v3_sub(&dist,&(jtom->r),&(itom->r));
565         if(bc) check_per_bound(moldyn,&dist);
566         d=v3_absolute_square(&dist);
567
568         params=moldyn->pot_params;
569         brand=itom->brand;
570
571         if(brand==jtom->brand) {
572                 if(d<=params->S2[brand])
573                         return TRUE;
574         }
575         else {
576                 if(d<=params->S2mixed)
577                         return TRUE;
578         }
579
580         return FALSE;
581 }