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_set_params(t_moldyn *moldyn,int element1,int element2) {
25 t_tersoff_mult_params *p;
27 /* alloc mem for potential parameters */
28 moldyn->pot_params=malloc(sizeof(t_tersoff_mult_params));
29 if(moldyn->pot_params==NULL) {
30 perror("[tersoff] pot params alloc");
34 /* these are now tersoff parameters */
37 // only 1 combination by now :p
45 p->lambda[0]=TM_LAMBDA_SI;
47 p->beta[0]=TM_BETA_SI;
57 printf("[tersoff] WARNING: element2\n");
62 printf("[tersoff] WARNING: element1\n");
73 p->lambda[1]=TM_LAMBDA_C;
82 printf("[tersoff] WARNING: element1\n");
86 printf("[tersoff] parameter completion\n");
87 p->S2[0]=p->S[0]*p->S[0];
88 p->S2[1]=p->S[1]*p->S[1];
89 p->Smixed=sqrt(p->S[0]*p->S[1]);
90 p->S2mixed=p->Smixed*p->Smixed;
91 p->Rmixed=sqrt(p->R[0]*p->R[1]);
92 p->Amixed=sqrt(p->A[0]*p->A[1]);
93 p->Bmixed=sqrt(p->B[0]*p->B[1]);
94 p->lambda_m=0.5*(p->lambda[0]+p->lambda[1]);
95 p->mu_m=0.5*(p->mu[0]+p->mu[1]);
97 printf("[tersoff] mult parameter info:\n");
98 printf(" S (A) | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed);
99 printf(" R (A) | %f | %f | %f\n",p->R[0],p->R[1],p->Rmixed);
100 printf(" A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV);
101 printf(" B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV);
102 printf(" lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1],
104 printf(" mu | %f | %f | %f\n",p->mu[0],p->mu[1],p->mu_m);
105 printf(" beta | %.10f | %.10f\n",p->beta[0],p->beta[1]);
106 printf(" n | %f | %f\n",p->n[0],p->n[1]);
107 printf(" c | %f | %f\n",p->c[0],p->c[1]);
108 printf(" d | %f | %f\n",p->d[0],p->d[1]);
109 printf(" h | %f | %f\n",p->h[0],p->h[1]);
110 printf(" chi | %f \n",p->chi);
115 /* tersoff 1 body part */
116 int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) {
119 t_tersoff_mult_params *params;
120 t_tersoff_exchange *exchange;
123 params=moldyn->pot_params;
124 exchange=&(params->exchange);
127 * simple: point constant parameters only depending on atom i to
131 exchange->beta_i=&(params->beta[brand]);
132 exchange->n_i=&(params->n[brand]);
133 exchange->c_i=&(params->c[brand]);
134 exchange->d_i=&(params->d[brand]);
135 exchange->h_i=&(params->h[brand]);
137 exchange->betaini=pow(*(exchange->beta_i),*(exchange->n_i));
138 exchange->ci2=params->c[brand]*params->c[brand];
139 exchange->di2=params->d[brand]*params->d[brand];
140 exchange->ci2di2=exchange->ci2/exchange->di2;
145 /* tersoff 2 body part */
146 int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
148 t_tersoff_mult_params *params;
149 t_3dvec dist_ij,force;
151 double A,R,S,S2,lambda;
159 printf("WARNING! - tersoff_mult_2bp is obsolete.\n");
160 printf("WARNING! - repulsive part handled in 3bp/j2 routine.\n");
162 /* use newtons third law */
165 params=moldyn->pot_params;
168 /* determine cutoff square */
170 S2=params->S2[brand];
175 v3_sub(&dist_ij,&(aj->r),&(ai->r));
176 if(bc) check_per_bound(moldyn,&dist_ij);
177 d_ij2=v3_absolute_square(&dist_ij);
179 /* if d_ij2 > S2 => no force & potential energy contribution */
184 /* now we will need the distance */
188 if(brand==ai->brand) {
192 lambda=params->lambda[brand];
198 lambda=params->lambda_m;
201 /* f_r_ij, df_r_ij */
202 f_r=A*exp(-lambda*d_ij);
203 df_r=lambda*f_r/d_ij;
209 v3_scale(&force,&dist_ij,-df_r);
213 arg=M_PI*(d_ij-R)/s_r;
214 f_c=0.5+0.5*cos(arg);
215 df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
216 v3_scale(&force,&dist_ij,-df_c*f_r-df_r*f_c);
220 v3_add(&(ai->f),&(ai->f),&force);
221 v3_sub(&(aj->f),&(aj->f),&force); // reason: dri rij = - drj rij
224 if((ai==&(moldyn->atom[0]))|(aj==&(moldyn->atom[0]))) {
225 printf("force 2bp: [%d %d]\n",ai->tag,aj->tag);
226 printf("adding %f %f %f\n",force.x,force.y,force.z);
227 if(ai==&(moldyn->atom[0]))
228 printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
229 if(aj==&(moldyn->atom[0]))
230 printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
235 virial_calc(ai,&force,&dist_ij);
237 /* energy 2bp contribution */
239 moldyn->energy+=energy;
246 /* tersoff 3 body potential function (first ij loop) */
247 int tersoff_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
249 t_tersoff_mult_params *params;
250 t_tersoff_exchange *exchange;
256 params=moldyn->pot_params;
257 exchange=&(params->exchange);
260 exchange->zeta_ij=0.0;
263 * set ij depending values
269 S2=params->S2[brand];
274 v3_sub(&dist_ij,&(aj->r),&(ai->r));
275 if(bc) check_per_bound(moldyn,&dist_ij);
276 d_ij2=v3_absolute_square(&dist_ij);
278 /* if d_ij2 > S2 => no force & potential energy contribution */
288 exchange->dist_ij=dist_ij;
289 exchange->d_ij2=d_ij2;
292 /* reset k counter for first k loop */
298 /* tersoff 3 body potential function (first k loop) */
299 int tersoff_mult_3bp_k1(t_moldyn *moldyn,
300 t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
302 t_tersoff_mult_params *params;
303 t_tersoff_exchange *exchange;
306 t_3dvec dist_ij,dist_ik;
307 double d_ik2,d_ik,d_ij;
308 double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
309 double f_c_ik,df_c_ik;
312 params=moldyn->pot_params;
313 exchange=&(params->exchange);
314 kcount=exchange->kcount;
316 if(kcount>TERSOFF_MAXN) {
317 printf("FATAL: neighbours = %d\n",kcount);
318 printf(" -> %d %d %d\n",ai->tag,aj->tag,ak->tag);
323 if(brand==ak->brand) {
326 S2=params->S2[brand];
335 v3_sub(&dist_ik,&(ak->r),&(ai->r));
336 if(bc) check_per_bound(moldyn,&dist_ik);
337 d_ik2=v3_absolute_square(&dist_ik);
339 /* store data for second k loop */
340 exchange->dist_ik[kcount]=dist_ik;
341 exchange->d_ik2[kcount]=d_ik2;
343 /* return if not within cutoff */
353 dist_ij=exchange->dist_ij;
357 cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
360 h_cos=*(exchange->h_i)-cos_theta;
361 d2_h_cos2=exchange->di2+(h_cos*h_cos);
362 frac=exchange->ci2/d2_h_cos2;
363 g=1.0+exchange->ci2di2-frac;
364 dg=-2.0*frac*h_cos/d2_h_cos2;
366 /* zeta sum += f_c_ik * g_ijk */
368 exchange->zeta_ij+=g;
374 arg=M_PI*(d_ik-R)/s_r;
375 f_c_ik=0.5+0.5*cos(arg);
376 df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
377 exchange->zeta_ij+=f_c_ik*g;
381 if((ai==&(moldyn->atom[0]))|
382 (aj==&(moldyn->atom[864]))|
383 (ak==&(moldyn->atom[1003]))) {
384 printf(" -> %f %f %f\n",exchange->ci2di2,frac,h_cos);
388 /* store even more data for second k loop */
389 exchange->g[kcount]=g;
390 exchange->dg[kcount]=dg;
391 exchange->d_ik[kcount]=d_ik;
392 exchange->cos_theta[kcount]=cos_theta;
393 exchange->f_c_ik[kcount]=f_c_ik;
394 exchange->df_c_ik[kcount]=df_c_ik;
396 /* increase k counter */
402 int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
404 t_tersoff_mult_params *params;
405 t_tersoff_exchange *exchange;
407 double f_a,df_a,b,db,f_c,df_c;
418 params=moldyn->pot_params;
419 exchange=&(params->exchange);
422 if(brand==ai->brand) {
427 mu=params->mu[brand];
428 lambda=params->lambda[brand];
437 lambda=params->lambda_m;
450 arg=M_PI*(d_ij-R)/s_r;
451 f_c=0.5+0.5*cos(arg);
452 df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
456 f_a=-B*exp(-mu*d_ij);
460 f_r=A*exp(-lambda*d_ij);
461 df_r=lambda*f_r/d_ij;
464 if(exchange->zeta_ij==0.0) {
470 tmp=exchange->betaini*pow(exchange->zeta_ij,ni-1.0);
471 b=(1.0+exchange->zeta_ij*tmp);
472 db=chi*pow(b,-1.0/(2.0*ni)-1.0);
477 /* force contribution */
478 scale=-0.5*(f_c*(df_r+b*df_a)+df_c*(f_r+b*df_a));
479 v3_scale(&force,&(exchange->dist_ij),scale);
480 v3_add(&(ai->f),&(ai->f),&force);
481 v3_sub(&(aj->f),&(aj->f),&force); // dri rij = - drj rij
484 if((ai==&(moldyn->atom[0]))|(aj==&(moldyn->atom[0]))) {
485 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,aj->tag);
486 printf("adding %f %f %f\n",force.x,force.y,force.z);
487 if(ai==&(moldyn->atom[0]))
488 printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
489 if(aj==&(moldyn->atom[0]))
490 printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
491 printf("energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
493 printf(" %f %f %f\n",exchange->zeta_ij,.0,.0);
499 virial_calc(ai,&force,&(exchange->dist_ij));
501 /* dzeta prefactor = - 0.5 f_c f_a db */
502 exchange->pre_dzeta=-0.5*f_a*f_c*db;
504 /* energy contribution */
505 energy=0.5*f_c*(b*f_a+f_r);
506 moldyn->energy+=energy;
509 /* reset k counter for second k loop */
515 /* tersoff 3 body potential function (second k loop) */
516 int tersoff_mult_3bp_k2(t_moldyn *moldyn,
517 t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
519 t_tersoff_mult_params *params;
520 t_tersoff_exchange *exchange;
522 t_3dvec dist_ik,dist_ij;
523 double d_ik2,d_ik,d_ij2,d_ij;
526 double g,dg,cos_theta;
528 double f_c_ik,df_c_ik;
529 double dijdik_inv,fcdg,dfcg;
530 t_3dvec dcosdri,dcosdrj,dcosdrk;
533 params=moldyn->pot_params;
534 exchange=&(params->exchange);
535 kcount=exchange->kcount;
537 if(kcount>TERSOFF_MAXN)
538 printf("FATAL: neighbours!\n");
541 d_ik2=exchange->d_ik2[kcount];
545 S2=params->S2[brand];
549 /* return if d_ik > S */
555 /* prefactor dzeta */
556 pre_dzeta=exchange->pre_dzeta;
559 dist_ik=exchange->dist_ik[kcount];
560 d_ik=exchange->d_ik[kcount];
562 /* f_c_ik, df_c_ik */
563 f_c_ik=exchange->f_c_ik[kcount];
564 df_c_ik=exchange->df_c_ik[kcount];
566 /* dist_ij, d_ij2, d_ij */
567 dist_ij=exchange->dist_ij;
568 d_ij2=exchange->d_ij2;
571 /* g, dg, cos_theta */
572 g=exchange->g[kcount];
573 dg=exchange->dg[kcount];
574 cos_theta=exchange->cos_theta[kcount];
576 /* cos_theta derivatives wrt i,j,k */
577 dijdik_inv=1.0/(d_ij*d_ik);
578 v3_scale(&dcosdrj,&dist_ik,dijdik_inv);
579 v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
580 v3_add(&dcosdrj,&dcosdrj,&tmp);
581 v3_scale(&dcosdrk,&dist_ij,dijdik_inv);
582 v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
583 v3_add(&dcosdrk,&dcosdrk,&tmp);
584 v3_add(&dcosdri,&dcosdrj,&dcosdrk);
585 v3_scale(&dcosdri,&dcosdri,-1.0);
587 /* f_c_ik * dg, df_c_ik * g */
591 /* derivative wrt i */
592 v3_scale(&force,&dist_ik,dfcg);
593 v3_scale(&tmp,&dcosdri,fcdg);
594 v3_add(&force,&force,&tmp);
595 v3_scale(&force,&force,pre_dzeta);
597 /* force contribution */
598 v3_add(&(ai->f),&(ai->f),&force);
601 if(ai==&(moldyn->atom[0])) {
602 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
603 printf("adding %f %f %f\n",force.x,force.y,force.z);
604 printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
609 //virial_calc(ai,&force,&dist_ij);
611 /* derivative wrt j */
612 v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
614 /* force contribution */
615 v3_add(&(aj->f),&(aj->f),&force);
618 if(aj==&(moldyn->atom[0])) {
619 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
620 printf("adding %f %f %f\n",force.x,force.y,force.z);
621 printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
626 //v3_scale(&force,&force,-1.0);
628 virial_calc(ai,&force,&dist_ij);
630 /* derivative wrt k */
631 v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
632 v3_scale(&tmp,&dcosdrk,fcdg);
633 v3_add(&force,&force,&tmp);
634 v3_scale(&force,&force,pre_dzeta);
636 /* force contribution */
637 v3_add(&(ak->f),&(ak->f),&force);
640 if(ak==&(moldyn->atom[0])) {
641 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
642 printf("adding %f %f %f\n",force.x,force.y,force.z);
643 printf("total k: %f %f %f\n",ak->f.x,ak->f.y,ak->f.z);
648 //v3_scale(&force,&force,-1.0);
650 virial_calc(ai,&force,&dist_ik);
652 /* increase k counter */