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 // set cutoff before parameters (actually only necessary for some pots)
28 if(moldyn->cutoff==0.0) {
29 printf("[tersoff] WARNING: no cutoff!\n");
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");
40 /* these are now tersoff parameters */
43 // only 1 combination by now :p
51 p->lambda[0]=TM_LAMBDA_SI;
53 p->beta[0]=TM_BETA_SI;
63 printf("[tersoff] WARNING: element2\n");
68 printf("[tersoff] WARNING: element1\n");
79 p->lambda[1]=TM_LAMBDA_C;
88 printf("[tersoff] WARNING: element2\n");
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];
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],
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);
129 /* tersoff 2 body part */
130 int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
132 t_tersoff_mult_params *params;
133 t_3dvec dist_ij,force;
135 double A,R,S,S2,lambda;
143 printf("WARNING! - tersoff_mult_2bp is obsolete.\n");
144 printf("WARNING! - repulsive part handled in 3bp/j2 routine.\n");
146 /* use newtons third law */
149 params=moldyn->pot_params;
152 /* determine cutoff square */
154 S2=params->S2[brand];
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);
163 /* if d_ij2 > S2 => no force & potential energy contribution */
168 /* now we will need the distance */
172 if(brand==ai->brand) {
176 lambda=params->lambda[brand];
182 lambda=params->lambda_m;
185 /* f_r_ij, df_r_ij */
186 f_r=A*exp(-lambda*d_ij);
187 df_r=lambda*f_r/d_ij;
193 v3_scale(&force,&dist_ij,-df_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);
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);
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);
220 virial_calc(ai,&force,&dist_ij);
222 /* energy 2bp contribution */
224 moldyn->energy+=energy;
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) {
234 t_tersoff_mult_params *params;
235 t_tersoff_exchange *exchange;
241 params=moldyn->pot_params;
242 exchange=&(params->exchange);
245 exchange->zeta_ij=0.0;
248 * set ij depending values
254 S2=params->S2[brand];
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);
263 /* if d_ij2 > S2 => no force & potential energy contribution */
273 exchange->dist_ij=dist_ij;
274 exchange->d_ij2=d_ij2;
277 /* reset k counter for first k loop */
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) {
287 t_tersoff_mult_params *params;
288 t_tersoff_exchange *exchange;
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;
297 params=moldyn->pot_params;
298 exchange=&(params->exchange);
299 kcount=exchange->kcount;
301 if(kcount>TERSOFF_MAXN) {
302 printf("FATAL: neighbours = %d\n",kcount);
303 printf(" -> %d %d %d\n",ai->tag,aj->tag,ak->tag);
308 if(brand==ak->brand) {
311 S2=params->S2[brand];
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);
324 /* store data for second k loop */
325 exchange->dist_ik[kcount]=dist_ik;
326 exchange->d_ik2[kcount]=d_ik2;
328 /* return if not within cutoff */
338 dist_ij=exchange->dist_ij;
342 cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
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;
351 /* zeta sum += f_c_ik * g_ijk */
353 exchange->zeta_ij+=g;
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;
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);
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;
378 /* increase k counter */
384 int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
386 t_tersoff_mult_params *params;
387 t_tersoff_exchange *exchange;
389 double f_a,df_a,b,db,f_c,df_c;
400 params=moldyn->pot_params;
401 exchange=&(params->exchange);
404 if(brand==aj->brand) {
409 mu=params->mu[brand];
410 lambda=params->lambda[brand];
419 lambda=params->lambda_m;
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));
438 f_a=-B*exp(-mu*d_ij);
442 f_r=A*exp(-lambda*d_ij);
443 df_r=lambda*f_r/d_ij;
446 if(exchange->zeta_ij==0.0) {
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);
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);
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),
476 printf(" %f %f %f\n",exchange->zeta_ij,.0,.0);
481 virial_calc(ai,&force,&(exchange->dist_ij));
483 /* dzeta prefactor = - 0.5 f_c f_a db */
484 exchange->pre_dzeta=-0.5*f_a*f_c*db;
486 /* energy contribution */
487 energy=0.5*f_c*(b*f_a+f_r);
488 moldyn->energy+=energy;
491 /* reset k counter for second k loop */
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) {
501 t_tersoff_mult_params *params;
502 t_tersoff_exchange *exchange;
504 t_3dvec dist_ik,dist_ij;
505 double d_ik2,d_ik,d_ij2,d_ij;
508 double g,dg,cos_theta;
510 double f_c_ik,df_c_ik;
511 double dijdik_inv,fcdg,dfcg;
512 t_3dvec dcosdrj,dcosdrk;
515 params=moldyn->pot_params;
516 exchange=&(params->exchange);
517 kcount=exchange->kcount;
519 if(kcount>TERSOFF_MAXN)
520 printf("FATAL: neighbours!\n");
523 d_ik2=exchange->d_ik2[kcount];
527 S2=params->S2[brand];
531 /* return if d_ik > S */
537 /* prefactor dzeta */
538 pre_dzeta=exchange->pre_dzeta;
541 dist_ik=exchange->dist_ik[kcount];
542 d_ik=exchange->d_ik[kcount];
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];
548 /* dist_ij, d_ij2, d_ij */
549 dist_ij=exchange->dist_ij;
550 d_ij2=exchange->d_ij2;
553 /* g, dg, cos_theta */
554 g=exchange->g[kcount];
555 dg=exchange->dg[kcount];
556 cos_theta=exchange->cos_theta[kcount];
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);
567 /* f_c_ik * dg, df_c_ik * g */
571 /* derivative wrt j */
572 v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
574 /* force contribution */
575 v3_add(&(aj->f),&(aj->f),&force);
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);
588 virial_calc(ai,&force,&dist_ij);
590 /* force contribution to atom i */
591 v3_scale(&force,&force,-1.0);
592 v3_add(&(ai->f),&(ai->f),&force);
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);
600 /* force contribution */
601 v3_add(&(ak->f),&(ak->f),&force);
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);
614 virial_calc(ai,&force,&dist_ik);
616 /* force contribution to atom i */
617 v3_scale(&force,&force,-1.0);
618 v3_add(&(ai->f),&(ai->f),&force);
620 /* increase k counter */
627 int tersoff_mult_check_2b_bond(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
629 t_tersoff_mult_params *params;
634 v3_sub(&dist,&(aj->r),&(ai->r));
635 if(bc) check_per_bound(moldyn,&dist);
636 d=v3_absolute_square(&dist);
638 params=moldyn->pot_params;
641 if(brand==aj->brand) {
642 if(d<=params->S2[brand])
646 if(d<=params->S2mixed)