2 * albe_orig.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"
20 #include "albe_orig.h"
22 /* create mixed terms from parameters and set them */
23 int albe_orig_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
25 t_albe_orig_mult_params *p;
27 // set cutoff before parameters (actually only necessary for some pots)
28 if(moldyn->cutoff==0.0) {
29 printf("[albe orig] 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 orig] 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_m=p->c_mixed*p->c_mixed;
89 p->d_mixed=ALBE_D_SIC;
90 p->d2_m=p->d_mixed*p->d_mixed;
91 p->c2d2_m=p->c2_m/p->d2_m;
92 p->h_mixed=ALBE_H_SIC;
95 printf("[albe orig] WARNING: element2");
101 printf("[albe orig] WARNING: element1\n");
105 printf("[albe orig] parameter completion\n");
106 p->S2[0]=p->S[0]*p->S[0];
107 p->S2[1]=p->S[1]*p->S[1];
108 p->S2mixed=p->Smixed*p->Smixed;
110 printf("[albe orig] mult parameter info:\n");
111 printf(" S (A) | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed);
112 printf(" R (A) | %f | %f | %f\n",p->R[0],p->R[1],p->Rmixed);
113 printf(" A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV);
114 printf(" B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV);
115 printf(" lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1],
117 printf(" mu | %f | %f | %f\n",p->mu[0],p->mu[1],p->mu_m);
118 printf(" gamma | %f | %f | %f\n",p->gamma[0],p->gamma[1],p->gamma_m);
119 printf(" c | %f | %f | %f\n",p->c[0],p->c[1],p->c_mixed);
120 printf(" d | %f | %f | %f\n",p->d[0],p->d[1],p->d_mixed);
121 printf(" h | %f | %f | %f\n",p->h[0],p->h[1],p->h_mixed);
126 /* albe 3 body potential function (first ij loop) */
127 int albe_orig_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
129 t_albe_orig_mult_params *params;
130 t_albe_orig_exchange *exchange;
136 params=moldyn->pot_params;
137 exchange=&(params->exchange);
140 exchange->zeta_ij=0.0;
143 * set ij depending values
147 if(brand==aj->brand) {
148 S2=params->S2[brand];
155 v3_sub(&dist_ij,&(aj->r),&(ai->r));
156 if(bc) check_per_bound(moldyn,&dist_ij);
157 d_ij2=v3_absolute_square(&dist_ij);
159 /* if d_ij2 > S2 => no force & potential energy contribution */
169 exchange->dist_ij=dist_ij;
170 exchange->d_ij2=d_ij2;
173 /* reset k counter for first k loop */
179 /* albe 3 body potential function (first k loop) */
180 int albe_orig_mult_3bp_k1(t_moldyn *moldyn,
181 t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
183 t_albe_orig_mult_params *params;
184 t_albe_orig_exchange *exchange;
187 t_3dvec dist_ij,dist_ik;
188 double d_ik2,d_ik,d_ij;
189 double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
190 double f_c_ik,df_c_ik;
193 /* continue if aj equals ak */
197 params=moldyn->pot_params;
198 exchange=&(params->exchange);
199 kcount=exchange->kcount;
201 if(kcount>ALBE_ORIG_MAXN) {
202 printf("FATAL: neighbours = %d\n",kcount);
203 printf(" -> %d %d %d\n",ai->tag,aj->tag,ak->tag);
208 if(brand==ak->brand) {
211 S2=params->S2[brand];
212 /* albe needs i,k depending c,d,h and gamma values */
213 exchange->gamma_i=&(params->gamma[brand]);
214 exchange->c_i=&(params->c[brand]);
215 exchange->d_i=&(params->d[brand]);
216 exchange->h_i=&(params->h[brand]);
217 exchange->ci2=&(params->c2[brand]);
218 exchange->di2=&(params->d2[brand]);
219 exchange->ci2di2=&(params->c2d2[brand]);
225 /* albe needs i,k depending c,d,h and gamma values */
226 exchange->gamma_i=&(params->gamma_m);
227 exchange->c_i=&(params->c_mixed);
228 exchange->d_i=&(params->d_mixed);
229 exchange->h_i=&(params->h_mixed);
230 exchange->ci2=&(params->c2_m);
231 exchange->di2=&(params->d2_m);
232 exchange->ci2di2=&(params->c2d2_m);
236 v3_sub(&dist_ik,&(ak->r),&(ai->r));
237 if(bc) check_per_bound(moldyn,&dist_ik);
238 d_ik2=v3_absolute_square(&dist_ik);
240 /* store data for second k loop */
241 exchange->dist_ik[kcount]=dist_ik;
242 exchange->d_ik2[kcount]=d_ik2;
244 /* return if not within cutoff */
254 dist_ij=exchange->dist_ij;
258 cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
261 h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism
262 d2_h_cos2=*(exchange->di2)+(h_cos*h_cos);
263 frac=*(exchange->ci2)/d2_h_cos2;
264 g=*(exchange->gamma_i)*(1.0+*(exchange->ci2di2)-frac);
265 dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
267 /* zeta sum += f_c_ik * g_ijk */
269 exchange->zeta_ij+=g;
275 arg=M_PI*(d_ik-R)/s_r;
276 f_c_ik=0.5+0.5*cos(arg);
277 df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
278 exchange->zeta_ij+=f_c_ik*g;
281 /* store even more data for second k loop */
282 exchange->g[kcount]=g;
283 exchange->dg[kcount]=dg;
284 exchange->d_ik[kcount]=d_ik;
285 exchange->cos_theta[kcount]=cos_theta;
286 exchange->f_c_ik[kcount]=f_c_ik;
287 exchange->df_c_ik[kcount]=df_c_ik;
289 /* increase k counter */
295 int albe_orig_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
297 t_albe_orig_mult_params *params;
298 t_albe_orig_exchange *exchange;
300 double f_a,df_a,b,db,f_c,df_c;
310 params=moldyn->pot_params;
311 exchange=&(params->exchange);
314 if(brand==ai->brand) {
319 r0=params->r0[brand];
320 mu=params->mu[brand];
321 lambda=params->lambda[brand];
330 lambda=params->lambda_m;
342 arg=M_PI*(d_ij-R)/s_r;
343 f_c=0.5+0.5*cos(arg);
344 df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
348 f_a=-B*exp(-mu*(d_ij-r0));
352 f_r=A*exp(-lambda*(d_ij-r0));
353 df_r=lambda*f_r/d_ij;
356 if(exchange->zeta_ij==0.0) {
361 b=1.0/sqrt(1.0+exchange->zeta_ij);
362 db=-0.5*b/(1.0+exchange->zeta_ij);
365 /* force contribution for atom i */
366 scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
367 v3_scale(&force,&(exchange->dist_ij),scale);
368 v3_add(&(ai->f),&(ai->f),&force);
370 /* force contribution for atom j */
371 v3_scale(&force,&force,-1.0); // dri rij = - drj rij
372 v3_add(&(aj->f),&(aj->f),&force);
375 virial_calc(ai,&force,&(exchange->dist_ij));
378 if(moldyn->time>DSTART&&moldyn->time<DEND) {
379 if((ai==&(moldyn->atom[DATOM]))|(aj==&(moldyn->atom[DATOM]))) {
380 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,aj->tag);
381 printf(" adding %f %f %f\n",force.x,force.y,force.z);
382 if(ai==&(moldyn->atom[0]))
383 printf(" total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
384 if(aj==&(moldyn->atom[0]))
385 printf(" total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
386 printf(" energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
388 printf(" %f %f %f\n",exchange->zeta_ij,.0,.0);
393 /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
394 exchange->pre_dzeta=0.5*f_a*f_c*db;
396 /* energy contribution */
397 energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
398 moldyn->energy+=energy;
401 /* reset k counter for second k loop */
407 /* albe 3 body potential function (second k loop) */
408 int albe_orig_mult_3bp_k2(t_moldyn *moldyn,
409 t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
411 t_albe_orig_mult_params *params;
412 t_albe_orig_exchange *exchange;
414 t_3dvec dist_ik,dist_ij;
415 double d_ik2,d_ik,d_ij2,d_ij;
418 double g,dg,cos_theta;
420 double f_c_ik,df_c_ik;
421 double dijdik_inv,fcdg,dfcg;
422 t_3dvec dcosdrj,dcosdrk;
425 /* continue if aj equals ak */
429 params=moldyn->pot_params;
430 exchange=&(params->exchange);
431 kcount=exchange->kcount;
433 if(kcount>ALBE_ORIG_MAXN)
434 printf("FATAL: neighbours!\n");
437 d_ik2=exchange->d_ik2[kcount];
441 S2=params->S2[brand];
445 /* return if d_ik > S */
451 /* prefactor dzeta */
452 pre_dzeta=exchange->pre_dzeta;
455 dist_ik=exchange->dist_ik[kcount];
456 d_ik=exchange->d_ik[kcount];
458 /* f_c_ik, df_c_ik */
459 f_c_ik=exchange->f_c_ik[kcount];
460 df_c_ik=exchange->df_c_ik[kcount];
462 /* dist_ij, d_ij2, d_ij */
463 dist_ij=exchange->dist_ij;
464 d_ij2=exchange->d_ij2;
467 /* g, dg, cos_theta */
468 g=exchange->g[kcount];
469 dg=exchange->dg[kcount];
470 cos_theta=exchange->cos_theta[kcount];
472 /* cos_theta derivatives wrt j,k */
473 dijdik_inv=1.0/(d_ij*d_ik);
474 v3_scale(&dcosdrj,&dist_ik,dijdik_inv); // j
475 v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
476 v3_add(&dcosdrj,&dcosdrj,&tmp);
477 v3_scale(&dcosdrk,&dist_ij,dijdik_inv); // k
478 v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
479 v3_add(&dcosdrk,&dcosdrk,&tmp);
481 /* f_c_ik * dg, df_c_ik * g */
485 /* derivative wrt j */
486 v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
488 /* force contribution */
489 v3_add(&(aj->f),&(aj->f),&force);
492 if(moldyn->time>DSTART&&moldyn->time<DEND) {
493 if(aj==&(moldyn->atom[DATOM])) {
494 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
495 printf(" adding %f %f %f\n",force.x,force.y,force.z);
496 printf(" total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
497 printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
498 printf(" d ij ik = %f %f\n",d_ij,d_ik);
504 virial_calc(ai,&force,&dist_ij);
506 /* force contribution to atom i */
507 v3_scale(&force,&force,-1.0);
508 v3_add(&(ai->f),&(ai->f),&force);
510 /* derivative wrt k */
511 v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
512 v3_scale(&tmp,&dcosdrk,fcdg);
513 v3_add(&force,&force,&tmp);
514 v3_scale(&force,&force,pre_dzeta);
516 /* force contribution */
517 v3_add(&(ak->f),&(ak->f),&force);
520 if(moldyn->time>DSTART&&moldyn->time<DEND) {
521 if(ak==&(moldyn->atom[DATOM])) {
522 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
523 printf(" adding %f %f %f\n",force.x,force.y,force.z);
524 printf(" total k: %f %f %f\n",ak->f.x,ak->f.y,ak->f.z);
525 printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
526 printf(" d ij ik = %f %f\n",d_ij,d_ik);
532 virial_calc(ai,&force,&dist_ik);
534 /* force contribution to atom i */
535 v3_scale(&force,&force,-1.0);
536 v3_add(&(ai->f),&(ai->f),&force);
538 /* increase k counter */
545 int albe_orig_mult_check_2b_bond(t_moldyn *moldyn,
546 t_atom *itom,t_atom *jtom,u8 bc) {
548 t_albe_orig_mult_params *params;
553 v3_sub(&dist,&(jtom->r),&(itom->r));
554 if(bc) check_per_bound(moldyn,&dist);
555 d=v3_absolute_square(&dist);
557 params=moldyn->pot_params;
560 if(brand==jtom->brand) {
561 if(d<=params->S2[brand])
565 if(d<=params->S2mixed)