]> hackdaworld.org Git - physik/posic.git/blob - potentials/albe.c
changed state machine behaviour
[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_complete_params(t_albe_mult_params *p) {
24
25         printf("[moldyn] albe parameter completion\n");
26         p->S2[0]=p->S[0]*p->S[0];
27         p->S2[1]=p->S[1]*p->S[1];
28         p->S2mixed=p->Smixed*p->Smixed;
29
30         printf("[moldyn] albe mult parameter info:\n");
31         printf("  S (A)  | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed);
32         printf("  R (A)  | %f | %f | %f\n",p->R[0],p->R[1],p->Rmixed);
33         printf("  A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV);
34         printf("  B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV);
35         printf("  lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1],
36                                           p->lambda_m);
37         printf("  mu     | %f | %f | %f\n",p->mu[0],p->mu[1],p->mu_m);
38         printf("  gamma  | %f | %f\n",p->gamma[0],p->gamma[1]);
39         printf("  c      | %f | %f\n",p->c[0],p->c[1]);
40         printf("  d      | %f | %f\n",p->d[0],p->d[1]);
41         printf("  h      | %f | %f\n",p->h[0],p->h[1]);
42
43         return 0;
44 }
45
46 /* albe 3 body potential function (first ij loop) */
47 int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
48
49         t_albe_mult_params *params;
50         t_albe_exchange *exchange;
51         unsigned char brand;
52         double S2;
53         t_3dvec dist_ij;
54         double d_ij2,d_ij;
55
56         params=moldyn->pot_params;
57         exchange=&(params->exchange);
58
59         /* reset zeta sum */
60         exchange->zeta_ij=0.0;
61
62         /*
63          * set ij depending values
64          */
65
66         brand=ai->brand;
67         if(brand==aj->brand) {
68                 S2=params->S2[brand];
69         }
70         else {
71                 S2=params->S2mixed;
72         }
73
74         /* dist_ij, d_ij2 */
75         v3_sub(&dist_ij,&(aj->r),&(ai->r));
76         if(bc) check_per_bound(moldyn,&dist_ij);
77         d_ij2=v3_absolute_square(&dist_ij);
78
79         /* if d_ij2 > S2 => no force & potential energy contribution */
80         if(d_ij2>S2) {
81                 moldyn->run3bp=0;
82                 return 0;
83         }
84
85         /* d_ij */
86         d_ij=sqrt(d_ij2);
87
88         /* store values */
89         exchange->dist_ij=dist_ij;
90         exchange->d_ij2=d_ij2;
91         exchange->d_ij=d_ij;
92
93         /* reset k counter for first k loop */
94         exchange->kcount=0;
95                 
96         return 0;
97 }
98
99 /* albe 3 body potential function (first k loop) */
100 int albe_mult_3bp_k1(t_moldyn *moldyn,
101                      t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
102
103         t_albe_mult_params *params;
104         t_albe_exchange *exchange;
105         unsigned char brand;
106         double R,S,S2;
107         t_3dvec dist_ij,dist_ik;
108         double d_ik2,d_ik,d_ij;
109         double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
110         double f_c_ik,df_c_ik;
111         int kcount;
112
113         params=moldyn->pot_params;
114         exchange=&(params->exchange);
115         kcount=exchange->kcount;
116
117         if(kcount>ALBE_MAXN) {
118                 printf("FATAL: neighbours = %d\n",kcount);
119                 printf("  -> %d %d %d\n",ai->tag,aj->tag,ak->tag);
120         }
121
122         /* ik constants */
123         brand=ai->brand;
124         if(brand==ak->brand) {
125                 R=params->R[brand];
126                 S=params->S[brand];
127                 S2=params->S2[brand];
128                 /* albe needs i,k depending c,d,h and gamma values */
129                 exchange->gamma_i=&(params->gamma[brand]);
130                 exchange->c_i=&(params->c[brand]);
131                 exchange->d_i=&(params->d[brand]);
132                 exchange->h_i=&(params->h[brand]);
133         }
134         else {
135                 R=params->Rmixed;
136                 S=params->Smixed;
137                 S2=params->S2mixed;
138                 /* albe needs i,k depending c,d,h and gamma values */
139                 exchange->gamma_i=&(params->gamma_m);
140                 exchange->c_i=&(params->c_mixed);
141                 exchange->d_i=&(params->d_mixed);
142                 exchange->h_i=&(params->h_mixed);
143         }
144         exchange->ci2=*(exchange->c_i)**(exchange->c_i);
145         exchange->di2=*(exchange->d_i)**(exchange->d_i);
146         exchange->ci2di2=exchange->ci2/exchange->di2;
147
148         /* dist_ik, d_ik2 */
149         v3_sub(&dist_ik,&(ak->r),&(ai->r));
150         if(bc) check_per_bound(moldyn,&dist_ik);
151         d_ik2=v3_absolute_square(&dist_ik);
152
153         /* store data for second k loop */
154         exchange->dist_ik[kcount]=dist_ik;
155         exchange->d_ik2[kcount]=d_ik2;
156
157         /* return if not within cutoff */
158         if(d_ik2>S2) {
159                 exchange->kcount++;
160                 return 0;
161         }
162
163         /* d_ik */
164         d_ik=sqrt(d_ik2);
165
166         /* dist_ij, d_ij */
167         dist_ij=exchange->dist_ij;
168         d_ij=exchange->d_ij;
169
170         /* cos theta */
171         cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
172
173         /* g_ijk */
174         h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism
175         d2_h_cos2=exchange->di2+(h_cos*h_cos);
176         frac=exchange->ci2/d2_h_cos2;
177         g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
178         dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
179
180         /* zeta sum += f_c_ik * g_ijk */
181         if(d_ik<=R) {
182                 exchange->zeta_ij+=g;
183                 f_c_ik=1.0;
184                 df_c_ik=0.0;
185         }
186         else {
187                 s_r=S-R;
188                 arg=M_PI*(d_ik-R)/s_r;
189                 f_c_ik=0.5+0.5*cos(arg);
190                 df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
191                 exchange->zeta_ij+=f_c_ik*g;
192         }
193
194         /* store even more data for second k loop */
195         exchange->g[kcount]=g;
196         exchange->dg[kcount]=dg;
197         exchange->d_ik[kcount]=d_ik;
198         exchange->cos_theta[kcount]=cos_theta;
199         exchange->f_c_ik[kcount]=f_c_ik;
200         exchange->df_c_ik[kcount]=df_c_ik;
201
202         /* increase k counter */
203         exchange->kcount++;
204
205         return 0;
206 }
207
208 int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
209
210         t_albe_mult_params *params;
211         t_albe_exchange *exchange;
212         t_3dvec force;
213         double f_a,df_a,b,db,f_c,df_c;
214         double f_r,df_r;
215         double scale;
216         double mu,B;
217         double lambda,A;
218         double d_ij,r0;
219         unsigned char brand;
220         double S,R,s_r,arg;
221         double energy;
222
223         params=moldyn->pot_params;
224         exchange=&(params->exchange);
225
226         brand=aj->brand;
227         if(brand==ai->brand) {
228                 S=params->S[brand];
229                 R=params->R[brand];
230                 B=params->B[brand];
231                 A=params->A[brand];
232                 r0=params->r0[brand];
233                 mu=params->mu[brand];
234                 lambda=params->lambda[brand];
235         }
236         else {
237                 S=params->Smixed;
238                 R=params->Rmixed;
239                 B=params->Bmixed;
240                 A=params->Amixed;
241                 r0=params->r0_mixed;
242                 mu=params->mu_m;
243                 lambda=params->lambda_m;
244         }
245
246         d_ij=exchange->d_ij;
247
248         /* f_c, df_c */
249         if(d_ij<R) {
250                 f_c=1.0;
251                 df_c=0.0;
252         }
253         else {
254                 s_r=S-R;
255                 arg=M_PI*(d_ij-R)/s_r;
256                 f_c=0.5+0.5*cos(arg);
257                 df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
258         }
259
260         /* f_a, df_a */
261         f_a=-B*exp(-mu*(d_ij-r0));
262         df_a=mu*f_a/d_ij;
263
264         /* f_r, df_r */
265         f_r=A*exp(-lambda*(d_ij-r0));
266         df_r=lambda*f_r/d_ij;
267
268         /* b, db */
269         if(exchange->zeta_ij==0.0) {
270                 b=1.0;
271                 db=0.0;
272         }
273         else {
274                 b=1.0/sqrt(1.0+exchange->zeta_ij);
275                 db=-0.5*b/(1.0+exchange->zeta_ij);
276         }
277
278         /* force contribution */
279         scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a));
280         v3_scale(&force,&(exchange->dist_ij),scale);
281         v3_add(&(ai->f),&(ai->f),&force);
282         v3_sub(&(aj->f),&(aj->f),&force); // dri rij = - drj rij
283
284 #ifdef DEBUG
285 if(moldyn->time>DSTART&&moldyn->time<DEND) {
286         if((ai==&(moldyn->atom[DATOM]))|(aj==&(moldyn->atom[DATOM]))) {
287                 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,aj->tag);
288                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
289                 if(ai==&(moldyn->atom[0]))
290                         printf("  total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
291                 if(aj==&(moldyn->atom[0]))
292                         printf("  total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
293                 printf("  energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
294                                                     f_c,b,f_a,f_r);
295                 printf("          %f %f %f\n",exchange->zeta_ij,.0,.0);
296         }
297 }
298 #endif
299
300         /* virial */
301         virial_calc(ai,&force,&(exchange->dist_ij));
302
303         /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
304         exchange->pre_dzeta=0.5*f_a*f_c*db;
305
306         /* energy contribution */
307         energy=0.5*f_c*(f_r-b*f_a);
308         moldyn->energy+=energy;
309         ai->e+=energy;
310
311         /* reset k counter for second k loop */
312         exchange->kcount=0;
313                 
314         return 0;
315 }
316
317 /* albe 3 body potential function (second k loop) */
318 int albe_mult_3bp_k2(t_moldyn *moldyn,
319                      t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
320
321         t_albe_mult_params *params;
322         t_albe_exchange *exchange;
323         int kcount;
324         t_3dvec dist_ik,dist_ij;
325         double d_ik2,d_ik,d_ij2,d_ij;
326         unsigned char brand;
327         double S2;
328         double g,dg,cos_theta;
329         double pre_dzeta;
330         double f_c_ik,df_c_ik;
331         double dijdik_inv,fcdg,dfcg;
332         t_3dvec dcosdri,dcosdrj,dcosdrk;
333         t_3dvec force,tmp;
334
335         params=moldyn->pot_params;
336         exchange=&(params->exchange);
337         kcount=exchange->kcount;
338
339         if(kcount>ALBE_MAXN)
340                 printf("FATAL: neighbours!\n");
341
342         /* d_ik2 */
343         d_ik2=exchange->d_ik2[kcount];
344
345         brand=ak->brand;
346         if(brand==ai->brand)
347                 S2=params->S2[brand];
348         else
349                 S2=params->S2mixed;
350
351         /* return if d_ik > S */
352         if(d_ik2>S2) {
353                 exchange->kcount++;
354                 return 0;
355         }
356
357         /* prefactor dzeta */
358         pre_dzeta=exchange->pre_dzeta;
359
360         /* dist_ik, d_ik */
361         dist_ik=exchange->dist_ik[kcount];
362         d_ik=exchange->d_ik[kcount];
363
364         /* f_c_ik, df_c_ik */
365         f_c_ik=exchange->f_c_ik[kcount];
366         df_c_ik=exchange->df_c_ik[kcount];
367
368         /* dist_ij, d_ij2, d_ij */
369         dist_ij=exchange->dist_ij;
370         d_ij2=exchange->d_ij2;
371         d_ij=exchange->d_ij;
372
373         /* g, dg, cos_theta */
374         g=exchange->g[kcount];
375         dg=exchange->dg[kcount];
376         cos_theta=exchange->cos_theta[kcount];
377
378         /* cos_theta derivatives wrt i,j,k */
379         dijdik_inv=1.0/(d_ij*d_ik);
380         v3_scale(&dcosdrj,&dist_ik,dijdik_inv);         // j
381         v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
382         v3_add(&dcosdrj,&dcosdrj,&tmp);
383         v3_scale(&dcosdrk,&dist_ij,dijdik_inv);         // k
384         v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
385         v3_add(&dcosdrk,&dcosdrk,&tmp);
386         v3_add(&dcosdri,&dcosdrj,&dcosdrk);             // i
387         v3_scale(&dcosdri,&dcosdri,-1.0);
388
389         /* f_c_ik * dg, df_c_ik * g */
390         fcdg=f_c_ik*dg;
391         dfcg=df_c_ik*g;
392
393         /* derivative wrt i */
394         v3_scale(&force,&dist_ik,dfcg);
395         v3_scale(&tmp,&dcosdri,fcdg);
396         v3_add(&force,&force,&tmp);
397         v3_scale(&force,&force,pre_dzeta);
398
399         /* force contribution */
400         v3_add(&(ai->f),&(ai->f),&force);
401         
402 #ifdef DEBUG
403 if(moldyn->time>DSTART&&moldyn->time<DEND) {
404         if(ai==&(moldyn->atom[DATOM])) {
405                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
406                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
407                 printf("  total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
408                 printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
409                 printf("    d ij ik = %f %f\n",d_ij,d_ik);
410         }
411 }
412 #endif
413
414         /* virial */
415         //virial_calc(ai,&force,&dist_ij);
416
417         /* derivative wrt j */
418         v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
419
420         /* force contribution */
421         v3_add(&(aj->f),&(aj->f),&force);
422
423 #ifdef DEBUG
424 if(moldyn->time>DSTART&&moldyn->time<DEND) {
425         if(aj==&(moldyn->atom[DATOM])) {
426                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
427                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
428                 printf("  total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
429                 printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
430                 printf("    d ij ik = %f %f\n",d_ij,d_ik);
431         }
432 }
433 #endif
434
435         /* virial */
436         v3_scale(&force,&force,-1.0);
437         virial_calc(ai,&force,&dist_ij);
438
439         /* derivative wrt k */
440         v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
441         v3_scale(&tmp,&dcosdrk,fcdg);
442         v3_add(&force,&force,&tmp);
443         v3_scale(&force,&force,pre_dzeta);
444
445         /* force contribution */
446         v3_add(&(ak->f),&(ak->f),&force);
447
448 #ifdef DEBUG
449 if(moldyn->time>DSTART&&moldyn->time<DEND) {
450         if(ak==&(moldyn->atom[DATOM])) {
451                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
452                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
453                 printf("  total k: %f %f %f\n",ak->f.x,ak->f.y,ak->f.z);
454                 printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
455                 printf("    d ij ik = %f %f\n",d_ij,d_ik);
456         }
457 }
458 #endif
459
460         /* virial */
461         v3_scale(&force,&force,-1.0);
462         virial_calc(ai,&force,&dist_ik);
463         
464         /* increase k counter */
465         exchange->kcount++;     
466
467         return 0;
468
469 }