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