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