2 * tersoff.c - tersoff potential
4 * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
12 #include <sys/types.h>
18 #include "../moldyn.h"
19 #include "../math/math.h"
22 /* create mixed terms from parameters and set them */
23 int tersoff_mult_complete_params(t_tersoff_mult_params *p) {
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]);
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],
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);
54 /* tersoff 1 body part */
55 int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) {
58 t_tersoff_mult_params *params;
59 t_tersoff_exchange *exchange;
62 params=moldyn->pot_params;
63 exchange=&(params->exchange);
66 * simple: point constant parameters only depending on atom i to
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]);
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;
84 /* tersoff 2 body part */
85 int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
87 t_tersoff_mult_params *params;
88 t_3dvec dist_ij,force;
90 double A,R,S,S2,lambda;
98 printf("WARNING! - tersoff_mult_2bp is obsolete.\n");
99 printf("WARNING! - repulsive part handled in 3bp/j2 routine.\n");
101 /* use newtons third law */
104 params=moldyn->pot_params;
107 /* determine cutoff square */
109 S2=params->S2[brand];
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);
118 /* if d_ij2 > S2 => no force & potential energy contribution */
123 /* now we will need the distance */
127 if(brand==ai->brand) {
131 lambda=params->lambda[brand];
137 lambda=params->lambda_m;
140 /* f_r_ij, df_r_ij */
141 f_r=A*exp(-lambda*d_ij);
142 df_r=lambda*f_r/d_ij;
148 v3_scale(&force,&dist_ij,-df_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);
159 v3_add(&(ai->f),&(ai->f),&force);
160 v3_sub(&(aj->f),&(aj->f),&force); // reason: dri rij = - drj rij
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);
174 virial_calc(ai,&force,&dist_ij);
176 /* energy 2bp contribution */
178 moldyn->energy+=energy;
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) {
188 t_tersoff_mult_params *params;
189 t_tersoff_exchange *exchange;
195 params=moldyn->pot_params;
196 exchange=&(params->exchange);
199 exchange->zeta_ij=0.0;
202 * set ij depending values
208 S2=params->S2[brand];
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);
217 /* if d_ij2 > S2 => no force & potential energy contribution */
227 exchange->dist_ij=dist_ij;
228 exchange->d_ij2=d_ij2;
231 /* reset k counter for first k loop */
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) {
241 t_tersoff_mult_params *params;
242 t_tersoff_exchange *exchange;
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;
251 params=moldyn->pot_params;
252 exchange=&(params->exchange);
253 kcount=exchange->kcount;
255 if(kcount>TERSOFF_MAXN) {
256 printf("FATAL: neighbours = %d\n",kcount);
257 printf(" -> %d %d %d\n",ai->tag,aj->tag,ak->tag);
262 if(brand==ak->brand) {
265 S2=params->S2[brand];
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);
278 /* store data for second k loop */
279 exchange->dist_ik[kcount]=dist_ik;
280 exchange->d_ik2[kcount]=d_ik2;
282 /* return if not within cutoff */
292 dist_ij=exchange->dist_ij;
296 cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
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;
305 /* zeta sum += f_c_ik * g_ijk */
307 exchange->zeta_ij+=g;
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;
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);
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;
335 /* increase k counter */
341 int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
343 t_tersoff_mult_params *params;
344 t_tersoff_exchange *exchange;
346 double f_a,df_a,b,db,f_c,df_c;
357 params=moldyn->pot_params;
358 exchange=&(params->exchange);
361 if(brand==ai->brand) {
366 mu=params->mu[brand];
367 lambda=params->lambda[brand];
376 lambda=params->lambda_m;
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));
395 f_a=-B*exp(-mu*d_ij);
399 f_r=A*exp(-lambda*d_ij);
400 df_r=lambda*f_r/d_ij;
403 if(exchange->zeta_ij==0.0) {
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);
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
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),
432 printf(" %f %f %f\n",exchange->zeta_ij,.0,.0);
438 virial_calc(ai,&force,&(exchange->dist_ij));
440 /* dzeta prefactor = - 0.5 f_c f_a db */
441 exchange->pre_dzeta=-0.5*f_a*f_c*db;
443 /* energy contribution */
444 energy=0.5*f_c*(b*f_a+f_r);
445 moldyn->energy+=energy;
448 /* reset k counter for second k loop */
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) {
458 t_tersoff_mult_params *params;
459 t_tersoff_exchange *exchange;
461 t_3dvec dist_ik,dist_ij;
462 double d_ik2,d_ik,d_ij2,d_ij;
465 double g,dg,cos_theta;
467 double f_c_ik,df_c_ik;
468 double dijdik_inv,fcdg,dfcg;
469 t_3dvec dcosdri,dcosdrj,dcosdrk;
472 params=moldyn->pot_params;
473 exchange=&(params->exchange);
474 kcount=exchange->kcount;
476 if(kcount>TERSOFF_MAXN)
477 printf("FATAL: neighbours!\n");
480 d_ik2=exchange->d_ik2[kcount];
484 S2=params->S2[brand];
488 /* return if d_ik > S */
494 /* prefactor dzeta */
495 pre_dzeta=exchange->pre_dzeta;
498 dist_ik=exchange->dist_ik[kcount];
499 d_ik=exchange->d_ik[kcount];
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];
505 /* dist_ij, d_ij2, d_ij */
506 dist_ij=exchange->dist_ij;
507 d_ij2=exchange->d_ij2;
510 /* g, dg, cos_theta */
511 g=exchange->g[kcount];
512 dg=exchange->dg[kcount];
513 cos_theta=exchange->cos_theta[kcount];
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);
526 /* f_c_ik * dg, df_c_ik * g */
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);
536 /* force contribution */
537 v3_add(&(ai->f),&(ai->f),&force);
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);
548 //virial_calc(ai,&force,&dist_ij);
550 /* derivative wrt j */
551 v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
553 /* force contribution */
554 v3_add(&(aj->f),&(aj->f),&force);
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);
565 //v3_scale(&force,&force,-1.0);
567 virial_calc(ai,&force,&dist_ij);
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);
575 /* force contribution */
576 v3_add(&(ak->f),&(ak->f),&force);
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);
587 //v3_scale(&force,&force,-1.0);
589 virial_calc(ai,&force,&dist_ik);
591 /* increase k counter */