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