2 * albe.c - albe 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 albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
25 t_albe_mult_params *p;
27 // set cutoff before parameters (actually only necessary for some pots)
28 if(moldyn->cutoff==0.0) {
29 printf("[albe] WARNING: no cutoff!\n");
33 /* alloc mem for potential parameters */
34 moldyn->pot_params=malloc(sizeof(t_albe_mult_params));
35 if(moldyn->pot_params==NULL) {
36 perror("[albe] pot params alloc");
40 /* these are now albe parameters */
43 // only 1 combination by now :p
52 p->lambda[0]=ALBE_LAMBDA_SI;
54 p->gamma[0]=ALBE_GAMMA_SI;
56 p->c2[0]=p->c[0]*p->c[0];
58 p->d2[0]=p->d[0]*p->d[0];
59 p->c2d2[0]=p->c2[0]/p->d2[0];
69 p->lambda[1]=ALBE_LAMBDA_C;
71 p->gamma[1]=ALBE_GAMMA_C;
73 p->c2[1]=p->c[1]*p->c[1];
75 p->d2[1]=p->d[1]*p->d[1];
76 p->c2d2[1]=p->c2[1]/p->d2[1];
78 /* mixed type: silicon carbide */
83 p->r0_mixed=ALBE_R0_SIC;
84 p->lambda_m=ALBE_LAMBDA_SIC;
86 p->gamma_m=ALBE_GAMMA_SIC;
87 p->c_mixed=ALBE_C_SIC;
88 p->c2_mixed=p->c_mixed*p->c_mixed;
89 p->d_mixed=ALBE_D_SIC;
90 p->d2_mixed=p->d_mixed*p->d_mixed;
91 p->c2d2_m=p->c2_mixed/p->d2_mixed;
92 p->h_mixed=ALBE_H_SIC;
95 printf("[albe] WARNING: element2\n");
100 printf("[albe] WARNING: element1\n");
104 printf("[albe] parameter completion\n");
105 p->S2[0]=p->S[0]*p->S[0];
106 p->S2[1]=p->S[1]*p->S[1];
107 p->S2mixed=p->Smixed*p->Smixed;
108 p->c2[0]=p->c[0]*p->c[0];
109 p->c2[1]=p->c[1]*p->c[1];
110 p->c2_mixed=p->c_mixed*p->c_mixed;
111 p->d2[0]=p->d[0]*p->d[0];
112 p->d2[1]=p->d[1]*p->d[1];
113 p->d2_mixed=p->d_mixed*p->d_mixed;
114 p->c2d2[0]=p->c2[0]/p->d2[0];
115 p->c2d2[1]=p->c2[1]/p->d2[1];
116 p->c2d2_m=p->c2_mixed/p->d2_mixed;
118 printf("[albe] mult parameter info:\n");
119 printf(" S (A) | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed);
120 printf(" R (A) | %f | %f | %f\n",p->R[0],p->R[1],p->Rmixed);
121 printf(" A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV);
122 printf(" B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV);
123 printf(" lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1],
125 printf(" mu | %f | %f | %f\n",p->mu[0],p->mu[1],p->mu_m);
126 printf(" gamma | %f | %f | %f\n",p->gamma[0],p->gamma[1],p->gamma_m);
127 printf(" c | %f | %f | %f\n",p->c[0],p->c[1],p->c_mixed);
128 printf(" d | %f | %f | %f\n",p->d[0],p->d[1],p->d_mixed);
129 printf(" c2 | %f | %f | %f\n",p->c2[0],p->c2[1],p->c2_mixed);
130 printf(" d2 | %f | %f | %f\n",p->d2[0],p->d2[1],p->d2_mixed);
131 printf(" c2d2 | %f | %f | %f\n",p->c2d2[0],p->c2d2[1],p->c2d2_m);
132 printf(" h | %f | %f | %f\n",p->h[0],p->h[1],p->h_mixed);
138 int albe_mult_i0(t_moldyn *moldyn,t_atom *ai) {
140 t_albe_mult_params *params;
141 t_albe_exchange *exchange;
145 params=moldyn->pot_params;
146 exchange=&(params->exchange);
148 /* zero exchange values */
149 memset(exchange->zeta,0,ALBE_MAXN*sizeof(double));
150 for(i=0;i<ALBE_MAXN;i++)
151 memset(exchange->dzeta[i],0,ALBE_MAXN*sizeof(t_3dvec));
158 /* first j loop within first i loop */
159 int albe_mult_i0_j0(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
161 t_albe_mult_params *params;
162 t_albe_exchange *exchange;
164 double S2,S,R,d2,d,s_r,arg;
169 params=moldyn->pot_params;
170 exchange=&(params->exchange);
175 /* set ij depending values */
177 if(brand==aj->brand) {
178 S2=params->S2[brand];
185 v3_sub(&dist,&(aj->r),&(ai->r));
186 if(bc) check_per_bound(moldyn,&dist);
187 exchange->dist[j]=dist;
188 d2=v3_absolute_square(&dist);
191 /* if d_ij2 > S2 => no force & potential energy contribution */
200 /* more ij depending values */
201 if(brand==aj->brand) {
204 /* albe needs i,(j/k) depending c,d,h and gamma values */
205 exchange->gamma_[j]=&(params->gamma[brand]);
206 exchange->c_[j]=&(params->c[brand]);
207 exchange->d_[j]=&(params->d[brand]);
208 exchange->h_[j]=&(params->h[brand]);
209 exchange->c2_[j]=&(params->c2[brand]);
210 exchange->d2_[j]=&(params->d2[brand]);
211 exchange->c2d2_[j]=&(params->c2d2[brand]);
216 /* albe needs i,(j/k) depending c,d,h and gamma values */
217 exchange->gamma_[j]=&(params->gamma_m);
218 exchange->c_[j]=&(params->c_mixed);
219 exchange->d_[j]=&(params->d_mixed);
220 exchange->h_[j]=&(params->h_mixed);
221 exchange->c2_[j]=&(params->c2_mixed);
222 exchange->d2_[j]=&(params->d2_mixed);
223 exchange->c2d2_[j]=&(params->c2d2_m);
232 exchange->f_c[j]=1.0;
233 exchange->df_c[j]=0.0;
238 exchange->f_c[j]=0.5+0.5*cos(arg);
239 exchange->df_c[j]=0.5*sin(arg)*(M_PI/(s_r*d));
242 /* reset k counter */
248 /* first k loop within first j loop within first i loop */
249 int albe_mult_i0_j0_k0(t_moldyn *moldyn,
250 t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
252 t_albe_mult_params *params;
253 t_albe_exchange *exchange;
257 double dj,dk,djdk_inv,cos_theta;
258 double gj,dgj,h_cos_j,d2_h_cos2_j,frac_j;
259 double gk,dgk,h_cos_k,d2_h_cos2_k,frac_k;
260 t_3dvec dcosdrj,dcosdrk,tmp;
261 t_3dvec *dzjj,*dzkk,*dzjk,*dzkj;
263 params=moldyn->pot_params;
264 exchange=&(params->exchange);
271 /* k<j & check whether to run k */
275 printf("FATAL: too many neighbours! (%d)\n",k);
276 printf(" atom i:%d | j:%d | k:%d\n",ai->tag,aj->tag,ak->tag);
278 if((k>=j)|(exchange->skip[k])) {
284 distj=exchange->dist[j];
285 distk=exchange->dist[k];
288 djdk_inv=1.0/(dj*dk);
291 cos_theta=v3_scalar_product(&distj,&distk)*djdk_inv;
293 /* g(cos(theta)) ij and ik values */
294 h_cos_j=*(exchange->h_[j])+cos_theta; // + in albe formalism
295 d2_h_cos2_j=*exchange->d2_[j]+(h_cos_j*h_cos_j);
296 frac_j=*exchange->c2_[j]/d2_h_cos2_j;
297 gj=1.0+*exchange->c2d2_[j]-frac_j;
298 gj*=*(exchange->gamma_[j]);
299 dgj=*(exchange->gamma_[j])*2.0*frac_j*h_cos_j/d2_h_cos2_j; // + in albe
300 if(ak->brand==aj->brand) {
305 h_cos_k=*(exchange->h_[k])+cos_theta;
306 d2_h_cos2_k=*exchange->d2_[k]+(h_cos_k*h_cos_k);
307 frac_k=*exchange->c2_[k]/d2_h_cos2_k;
308 gk=1.0+*exchange->c2d2_[k]-frac_k;
309 gk*=*(exchange->gamma_[k]);
310 dgk=*(exchange->gamma_[k])*2.0*frac_k*h_cos_k/d2_h_cos2_k;
314 if(ai==&(moldyn->atom[DATOM]))
315 printf("zeta_ij: %f %f %f %f\n",f_c_ik*g,f_c_ik,g,d_ik);
318 /* store even more data for second k loop */
319 exchange->g[kcount]=g;
320 exchange->dg[kcount]=dg;
321 exchange->d_ik[kcount]=d_ik;
322 exchange->cos_theta[kcount]=cos_theta;
323 exchange->f_c_ik[kcount]=f_c_ik;
324 exchange->df_c_ik[kcount]=df_c_ik;
326 /* increase k counter */
332 /* first j loop within first i loop */
333 int albe_mult_i0_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
335 t_albe_mult_params *params;
336 t_albe_exchange *exchange;
338 params=moldyn->pot_params;
339 exchange=&(params->exchange);
341 /* increase j counter */
347 /* second j loop within first i loop */
348 int albe_mult_i0_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
350 t_albe_mult_params *params;
351 t_albe_exchange *exchange;
354 double d,f_a,df_a,f_r,df_r,f_c,df_c,b,db;
355 double A,B,mu,lambda,r0;
361 params=moldyn->pot_params;
362 exchange=&(params->exchange);
367 /* skip if j not within cutoff */
368 if(exchange->skip[j]) {
376 dist=&(exchange->dist[j]);
377 f_c=exchange->f_c[j];
378 df_c=exchange->df_c[j];
380 /* determine parameters to calculate fa, dfa, fr, dfr */
382 if(brand==ai->brand) {
385 r0=params->r0[brand];
386 mu=params->mu[brand];
387 lambda=params->lambda[brand];
394 lambda=params->lambda_m;
398 f_a=-B*exp(-mu*(d-r0));
402 f_r=A*exp(-lambda*(d-r0));
406 b=1.0/sqrt(1.0+exchange->zeta[j]);
407 db=-0.5*b/(1.0+exchange->zeta[j]);
409 /* energy contribution */
410 energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
411 moldyn->energy+=energy;
414 /* force contribution for atom i due to ij bond */
415 scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
416 v3_scale(&force,dist,scale);
417 v3_add(&(ai->f),&(ai->f),&force);
421 printf("force: %.15f %.15f %.15f | %d %d (ij) %.15f\n",force.x,force.y,force.z,ai->tag,aj->tag,exchange->zeta[j]);
422 printf(" t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z);
426 /* force contribution for atom j due to ij bond */
427 v3_scale(&force,&force,-1.0); // dri rij = - drj rij
428 v3_add(&(aj->f),&(aj->f),&force);
431 virial_calc(ai,&force,&(exchange->dist_ij));
434 if((ai==&(moldyn->atom[DATOM]))|(aj==&(moldyn->atom[DATOM]))) {
435 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,aj->tag);
436 printf(" adding %f %f %f\n",force.x,force.y,force.z);
437 if(ai==&(moldyn->atom[DATOM]))
438 printf(" total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
439 if(aj==&(moldyn->atom[DATOM]))
440 printf(" total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
441 printf(" energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
443 printf(" %f %f %f\n",exchange->zeta_ij,.0,.0);
448 virial_calc(ai,&force,dist);
450 /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
451 exchange->pre_dzeta=0.5*f_a*f_c*db;
453 /* force contribution (drj derivative) */
454 v3_scale(&force,&(exchange->dzeta[j][j]),exchange->pre_dzeta);
455 v3_add(&(aj->f),&(aj->f),&force);
459 printf("force: %.15f %.15f %.15f | %d %d (j der)\n",force.x,force.y,force.z,aj->tag,ai->tag);
460 printf(" t: %.15f %.15f %.15f\n",aj->f.x,aj->f.y,aj->f.z);
465 virial_calc(ai,&force,dist);
467 v3_scale(&force,&force,-1.0);
468 v3_add(&(ai->f),&(ai->f),&force);
472 printf("force: %.15f %.15f %.15f | %d %d (i contr j der)\n",force.x,force.y,force.z,ai->tag,aj->tag);
473 printf(" t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z);
477 /* reset k counter for second k loop */
483 /* second k loop within second j loop within first i loop */
484 int albe_mult_i0_j2_k0(t_moldyn *moldyn,
485 t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
487 t_albe_mult_params *params;
488 t_albe_exchange *exchange;
493 params=moldyn->pot_params;
494 exchange=&(params->exchange);
501 /* prefactor dzeta */
502 pre_dzeta=exchange->pre_dzeta;
505 dist_ik=exchange->dist_ik[kcount];
506 d_ik=exchange->d_ik[kcount];
508 /* f_c_ik, df_c_ik */
509 f_c_ik=exchange->f_c_ik[kcount];
510 df_c_ik=exchange->df_c_ik[kcount];
512 /* dist_ij, d_ij2, d_ij */
513 dist_ij=exchange->dist_ij;
514 d_ij2=exchange->d_ij2;
517 /* g, dg, cos_theta */
518 g=exchange->g[kcount];
519 dg=exchange->dg[kcount];
520 cos_theta=exchange->cos_theta[kcount];
522 /* cos_theta derivatives wrt j,k */
523 dijdik_inv=1.0/(d_ij*d_ik);
524 v3_scale(&dcosdrj,&dist_ik,dijdik_inv); // j
525 v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
526 v3_add(&dcosdrj,&dcosdrj,&tmp);
527 v3_scale(&dcosdrk,&dist_ij,dijdik_inv); // k
528 v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
529 v3_add(&dcosdrk,&dcosdrk,&tmp);
531 /* f_c_ik * dg, df_c_ik * g */
535 /* derivative wrt j */
536 v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
538 /* force contribution */
539 v3_add(&(aj->f),&(aj->f),&force);
542 if(aj==&(moldyn->atom[DATOM])) {
543 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
544 printf(" adding %f %f %f\n",force.x,force.y,force.z);
545 printf(" total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
546 printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
547 printf(" d ij ik = %f %f\n",d_ij,d_ik);
552 virial_calc(ai,&force,&dist_ij);
554 /* force contribution to atom i */
555 v3_scale(&force,&force,-1.0);
556 v3_add(&(ai->f),&(ai->f),&force);
558 /* derivative wrt k */
559 v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
560 v3_scale(&tmp,&dcosdrk,fcdg);
561 v3_add(&force,&force,&tmp);
562 v3_scale(&force,&force,pre_dzeta);
564 v3_scale(&force,&force,-1.0);
565 v3_add(&(ai->f),&(ai->f),&force);
568 if(ak==&(moldyn->atom[DATOM])) {
569 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
570 printf(" adding %f %f %f\n",force.x,force.y,force.z);
571 printf(" total k: %f %f %f\n",ak->f.x,ak->f.y,ak->f.z);
572 printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
573 printf(" d ij ik = %f %f\n",d_ij,d_ik);
578 virial_calc(ai,&force,&dist_ik);
580 /* force contribution to atom i */
581 v3_scale(&force,&force,-1.0);
582 v3_add(&(ai->f),&(ai->f),&force);
584 /* increase k counter */
590 int albe_mult_i0_j3(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
592 t_albe_mult_params *params;
593 t_albe_exchange *exchange;
595 params=moldyn->pot_params;
596 exchange=&(params->exchange);
598 /* increase j counter */
604 int albe_mult_check_2b_bond(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,u8 bc) {
606 t_albe_mult_params *params;
611 v3_sub(&dist,&(jtom->r),&(itom->r));
612 if(bc) check_per_bound(moldyn,&dist);
613 d=v3_absolute_square(&dist);
615 params=moldyn->pot_params;
618 if(brand==jtom->brand) {
619 if(d<=params->S2[brand])
623 if(d<=params->S2mixed)