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;
66 p->lambda[1]=ALBE_LAMBDA_C;
68 p->gamma[1]=ALBE_GAMMA_C;
72 /* mixed type: silicon carbide */
77 p->r0_mixed=ALBE_R0_SIC;
78 p->lambda_m=ALBE_LAMBDA_SIC;
80 p->gamma_m=ALBE_GAMMA_SIC;
81 p->c_mixed=ALBE_C_SIC;
82 p->d_mixed=ALBE_D_SIC;
83 p->h_mixed=ALBE_H_SIC;
86 printf("[albe] WARNING: element2\n");
91 printf("[albe] WARNING: element1\n");
95 printf("[albe] parameter completion\n");
96 p->S2[0]=p->S[0]*p->S[0];
97 p->S2[1]=p->S[1]*p->S[1];
98 p->S2mixed=p->Smixed*p->Smixed;
100 printf("[albe] mult parameter info:\n");
101 printf(" S (A) | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed);
102 printf(" R (A) | %f | %f | %f\n",p->R[0],p->R[1],p->Rmixed);
103 printf(" A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV);
104 printf(" B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV);
105 printf(" lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1],
107 printf(" mu | %f | %f | %f\n",p->mu[0],p->mu[1],p->mu_m);
108 printf(" gamma | %f | %f\n",p->gamma[0],p->gamma[1]);
109 printf(" c | %f | %f\n",p->c[0],p->c[1]);
110 printf(" d | %f | %f\n",p->d[0],p->d[1]);
111 printf(" h | %f | %f\n",p->h[0],p->h[1]);
117 int albe_mult_i0(t_moldyn *moldyn,t_atom *ai) {
120 t_albe_mult_params *params;
121 t_albe_exchange *exchange;
123 params=moldyn->pot_params;
124 exchange=&(params->exchange);
126 /* zero exchange values */
127 memset(exchange->dist,0,ALBE_MAXN*sizeof(t_3dvec));
128 memset(exchange->d2,0,ALBE_MAXN*sizeof(double));
129 memset(exchange->d,0,ALBE_MAXN*sizeof(double));
130 memset(exchange->zeta,0,ALBE_MAXN*sizeof(double));
131 for(i=0;i<ALBE_MAXN;i++)
132 memset(exchange->dzeta[i],0,ALBE_MAXN*sizeof(double));
133 memset(exchange->skip,0,ALBE_MAXN*sizeof(u8));
140 /* first j loop within first i loop */
141 int albe_mult_i0_j0(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
143 t_albe_mult_params *params;
144 t_albe_exchange *exchange;
145 double S2,S,R,d2,d,s_r,arg;
150 params=moldyn->pot_params;
151 exchange=&(params->exchange);
156 /* set ij depending values */
158 if(brand==aj->brand) {
159 S2=params->S2[brand];
166 v3_sub(&dist,&(aj->r),&(ai->r));
167 exchange->dist[j]=dist;
168 if(bc) check_per_bound(moldyn,&dist);
169 d2=v3_absolute_square(&dist);
172 /* if d_ij2 > S2 => no force & potential energy contribution */
180 /* more ij depending values */
181 if(brand==aj->brand) {
184 /* set albe needs i,(j/k) depending c,d,h and gamma values */
185 exchange->gamma_[j]=&(params->gamma[brand]);
186 exchange->c_[j]=&(params->c[brand]);
187 exchange->d_[j]=&(params->d[brand]);
188 exchange->h_[j]=&(params->h[brand]);
193 /* albe needs i,(j/k) depending c,d,h and gamma values */
194 exchange->gamma_[j]=&(params->gamma_m);
195 exchange->c_[j]=&(params->c_mixed);
196 exchange->d_[j]=&(params->d_mixed);
197 exchange->h_[j]=&(params->h_mixed);
199 exchange->c2_[j]=*exchange->c_[j]**exchange->c_[j];
200 exchange->d2_[j]=*exchange->d_[j]**exchange->d_[j];
201 exchange->c2d2_[j]=exchange->c2_[j]/exchange->d2_[j];
204 d=sqrt(exchange->d2[j]);
209 exchange->f_c[j]=1.0;
210 exchange->df_c[j]=0.0;
215 exchange->f_c[j]=0.5+0.5*cos(arg);
216 exchange->df_c[j]=0.5*sin(arg)*(M_PI/(s_r*d));
219 /* reset k counter */
225 /* first k loop within first j loop within first i loop */
226 int albe_mult_i0_j0_k0(t_moldyn *moldyn,
227 t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
229 t_albe_mult_params *params;
230 t_albe_exchange *exchange;
233 double dj,dk,djdk_inv,cos_theta;
234 double gj,dgj,h_cos_j,d2_h_cos2_j,frac_j;
235 double gk,dgk,h_cos_k,d2_h_cos2_k,frac_k;
236 t_3dvec dcosdrj,dcosdrk,tmp,**dzeta;
238 params=moldyn->pot_params;
239 exchange=&(params->exchange);
241 /* k<j & check whether to run k */
245 printf("FATAL: too many neighbours! (%d)\n",k);
246 printf(" atom i:%d | j:%d | k:%d\n",ai->tag,aj->tag,ak->tag);
248 if((k>=j)|(exchange->skip[k])) {
254 distj=exchange->dist[j];
255 distk=exchange->dist[k];
258 djdk_inv=1.0/(dj*dk);
261 cos_theta=v3_scalar_product(&distj,&distk)*djdk_inv;
263 /* g(cos(theta)) - in albe: ik-depending values! */
264 h_cos_j=*(exchange->h_[j])+cos_theta; // + in albe formalism
265 d2_h_cos2_j=exchange->d2_[j]+(h_cos_j*h_cos_j);
266 frac_j=exchange->c2_[j]/d2_h_cos2_j;
267 gj=1.0+exchange->c2d2_[j]-frac_j;
268 gj*=*(exchange->gamma_[j]);
269 dgj=*(exchange->gamma_[j])*2.0*frac_j*h_cos_j/d2_h_cos2_j; // + in albe
270 if(ak->brand==aj->brand) {
275 h_cos_k=*(exchange->h_[k])+cos_theta;
276 d2_h_cos2_k=exchange->d2_[k]+(h_cos_k*h_cos_k);
277 frac_k=exchange->c2_[k]/d2_h_cos2_k;
278 gk=1.0+exchange->c2d2_[k]-frac_k;
279 gk*=*(exchange->gamma_[k]);
280 dgk=*(exchange->gamma_[k])*2.0*frac_k*h_cos_k/d2_h_cos2_k;
284 exchange->zeta[j]+=(exchange->f_c[k]*gk);
285 exchange->zeta[k]+=(exchange->f_c[j]*gj);
287 /* cos theta derivatives */
288 v3_scale(&dcosdrj,&distk,djdk_inv); // j
289 v3_scale(&tmp,&distj,-cos_theta/exchange->d2[j]);
290 v3_add(&dcosdrj,&dcosdrj,&tmp);
291 v3_scale(&dcosdrk,&distj,djdk_inv); // k
292 v3_scale(&tmp,&distk,-cos_theta/exchange->d2[k]);
293 v3_add(&dcosdrk,&dcosdrk,&tmp);
295 /* zeta derivatives */
296 dzeta=exchange->dzeta;
297 v3_scale(&tmp,&dcosdrj,exchange->f_c[k]*dgk);
298 v3_add(&dzeta[j][j],&dzeta[j][j],&tmp); // j j
299 v3_scale(&tmp,&dcosdrk,exchange->f_c[j]*dgj);
300 v3_add(&dzeta[k][k],&dzeta[k][k],&tmp); // k k
301 v3_scale(&tmp,&distk,exchange->df_c[k]*gk/dk);
302 v3_add(&dzeta[j][k],&dzeta[j][k],&tmp);
303 v3_scale(&tmp,&dcosdrk,exchange->f_c[k]*dgk);
304 v3_add(&dzeta[j][k],&dzeta[j][k],&tmp); // j k
305 v3_scale(&tmp,&distj,exchange->df_c[j]*gj/dj);
306 v3_add(&dzeta[k][j],&dzeta[k][j],&tmp);
307 v3_scale(&tmp,&dcosdrj,exchange->f_c[j]*dgj);
308 v3_add(&dzeta[k][j],&dzeta[k][j],&tmp); // k j
310 /* increase k counter */
316 /* first j loop within first i loop */
317 int albe_mult_i0_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
319 t_albe_mult_params *params;
320 t_albe_exchange *exchange;
322 params=moldyn->pot_params;
323 exchange=&(params->exchange);
325 /* increase j counter */
331 /* second j loop within first i loop */
332 int albe_mult_i0_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
334 t_albe_mult_params *params;
335 t_albe_exchange *exchange;
338 double d,f_a,df_a,f_r,df_r,f_c,df_c,b,db;
339 double A,B,mu,lambda,r0;
345 params=moldyn->pot_params;
346 exchange=&(params->exchange);
351 /* skip if j not within cutoff */
352 if(exchange->skip[j]) {
360 dist=&(exchange->dist[j]);
361 f_c=exchange->f_c[j];
362 df_c=exchange->df_c[j];
364 /* determine parameters to calculate fa, dfa, fr, dfr */
366 if(brand==ai->brand) {
369 r0=params->r0[brand];
370 mu=params->mu[brand];
371 lambda=params->lambda[brand];
378 lambda=params->lambda_m;
382 f_a=-B*exp(-mu*(d-r0));
386 f_r=A*exp(-lambda*(d-r0));
390 b=1.0/sqrt(1.0+exchange->zeta[j]);
391 db=-0.5*b/(1.0+exchange->zeta[j]);
393 /* energy contribution */
394 energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
395 moldyn->energy+=energy;
398 /* force contribution for atom i due to ij bond */
399 scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
400 v3_scale(&force,dist,scale);
401 v3_add(&(ai->f),&(ai->f),&force);
403 /* force contribution for atom j due to ij bond */
404 v3_scale(&force,&force,-1.0); // dri rij = - drj rij
405 v3_add(&(aj->f),&(aj->f),&force);
408 virial_calc(aj,&force,dist);
410 /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
411 exchange->pre_dzeta=0.5*f_a*exchange->f_c[j]*db;
413 /* force contribution (drj derivative) */
414 v3_scale(&force,&(exchange->dzeta[j][j]),exchange->pre_dzeta);
415 v3_add(&(aj->f),&(aj->f),&force);
416 v3_scale(&force,&force,-1.0);
417 v3_add(&(ai->f),&(ai->f),&force);
420 virial_calc(ai,&force,dist);
422 /* reset k counter for second k loop */
428 /* second k loop within second j loop within first i loop */
429 int albe_mult_i0_j2_k0(t_moldyn *moldyn,
430 t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
432 t_albe_mult_params *params;
433 t_albe_exchange *exchange;
438 params=moldyn->pot_params;
439 exchange=&(params->exchange);
441 /* k!=j & check whether to run k */
444 if((k==j)|(exchange->skip[k])) {
449 /* force contribution (drk derivative) */
450 v3_scale(&force,&(exchange->dzeta[j][k]),exchange->pre_dzeta);
451 v3_add(&(ak->f),&(ak->f),&force);
452 v3_scale(&force,&force,-1.0);
453 v3_add(&(ai->f),&(ai->f),&force);
456 virial_calc(ai,&force,&(exchange->dist[k]));
458 /* increase k counter */
464 int albe_mult_i0_j3(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
466 t_albe_mult_params *params;
467 t_albe_exchange *exchange;
469 params=moldyn->pot_params;
470 exchange=&(params->exchange);
472 /* increase j counter */
478 int albe_mult_check_2b_bond(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,u8 bc) {
480 t_albe_mult_params *params;
485 v3_sub(&dist,&(jtom->r),&(itom->r));
486 if(bc) check_per_bound(moldyn,&dist);
487 d=v3_absolute_square(&dist);
489 params=moldyn->pot_params;
492 if(brand==jtom->brand) {
493 if(d<=params->S2[brand])
497 if(d<=params->S2mixed)