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