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_complete_params(t_albe_mult_params *p) {
25 printf("[moldyn] albe parameter completion\n");
26 p->S2[0]=p->S[0]*p->S[0];
27 p->S2[1]=p->S[1]*p->S[1];
28 p->S2mixed=p->Smixed*p->Smixed;
30 printf("[moldyn] albe mult parameter info:\n");
31 printf(" S (A) | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed);
32 printf(" R (A) | %f | %f | %f\n",p->R[0],p->R[1],p->Rmixed);
33 printf(" A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV);
34 printf(" B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV);
35 printf(" lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1],
37 printf(" mu | %f | %f | %f\n",p->mu[0],p->mu[1],p->mu_m);
38 printf(" gamma | %f | %f\n",p->gamma[0],p->gamma[1]);
39 printf(" c | %f | %f\n",p->c[0],p->c[1]);
40 printf(" d | %f | %f\n",p->d[0],p->d[1]);
41 printf(" h | %f | %f\n",p->h[0],p->h[1]);
46 /* albe 3 body potential function (first ij loop) */
47 int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
49 t_albe_mult_params *params;
50 t_albe_exchange *exchange;
56 params=moldyn->pot_params;
57 exchange=&(params->exchange);
60 exchange->zeta_ij=0.0;
63 * set ij depending values
67 if(brand==aj->brand) {
75 v3_sub(&dist_ij,&(aj->r),&(ai->r));
76 if(bc) check_per_bound(moldyn,&dist_ij);
77 d_ij2=v3_absolute_square(&dist_ij);
79 /* if d_ij2 > S2 => no force & potential energy contribution */
84 //if(ai==&(moldyn->atom[5832]))
85 // printf(" %d: %d (%f)\n",ai->tag,aj->tag,moldyn->time);
91 exchange->dist_ij=dist_ij;
92 exchange->d_ij2=d_ij2;
95 /* reset k counter for first k loop */
101 /* albe 3 body potential function (first k loop) */
102 int albe_mult_3bp_k1(t_moldyn *moldyn,
103 t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
105 t_albe_mult_params *params;
106 t_albe_exchange *exchange;
109 t_3dvec dist_ij,dist_ik;
110 double d_ik2,d_ik,d_ij;
111 double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
112 double f_c_ik,df_c_ik;
115 params=moldyn->pot_params;
116 exchange=&(params->exchange);
117 kcount=exchange->kcount;
119 if(kcount>ALBE_MAXN) {
120 printf("FATAL: neighbours = %d\n",kcount);
121 printf(" -> %d %d %d\n",ai->tag,aj->tag,ak->tag);
126 if(brand==ak->brand) {
129 S2=params->S2[brand];
130 /* albe needs i,k depending c,d,h and gamma values */
131 exchange->gamma_i=&(params->gamma[brand]);
132 exchange->c_i=&(params->c[brand]);
133 exchange->d_i=&(params->d[brand]);
134 exchange->h_i=&(params->h[brand]);
140 /* albe needs i,k depending c,d,h and gamma values */
141 exchange->gamma_i=&(params->gamma_m);
142 exchange->c_i=&(params->c_mixed);
143 exchange->d_i=&(params->d_mixed);
144 exchange->h_i=&(params->h_mixed);
146 exchange->ci2=*(exchange->c_i)**(exchange->c_i);
147 exchange->di2=*(exchange->d_i)**(exchange->d_i);
148 exchange->ci2di2=exchange->ci2/exchange->di2;
151 v3_sub(&dist_ik,&(ak->r),&(ai->r));
152 if(bc) check_per_bound(moldyn,&dist_ik);
153 d_ik2=v3_absolute_square(&dist_ik);
155 /* store data for second k loop */
156 exchange->dist_ik[kcount]=dist_ik;
157 exchange->d_ik2[kcount]=d_ik2;
159 /* return if not within cutoff */
169 dist_ij=exchange->dist_ij;
173 cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
176 h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism
177 d2_h_cos2=exchange->di2+(h_cos*h_cos);
178 frac=exchange->ci2/d2_h_cos2;
179 g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
180 dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
182 /* zeta sum += f_c_ik * g_ijk */
184 exchange->zeta_ij+=g;
190 arg=M_PI*(d_ik-R)/s_r;
191 f_c_ik=0.5+0.5*cos(arg);
192 df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
193 exchange->zeta_ij+=f_c_ik*g;
196 /* store even more data for second k loop */
197 exchange->g[kcount]=g;
198 exchange->dg[kcount]=dg;
199 exchange->d_ik[kcount]=d_ik;
200 exchange->cos_theta[kcount]=cos_theta;
201 exchange->f_c_ik[kcount]=f_c_ik;
202 exchange->df_c_ik[kcount]=df_c_ik;
204 /* increase k counter */
210 int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
212 t_albe_mult_params *params;
213 t_albe_exchange *exchange;
215 double f_a,df_a,b,db,f_c,df_c;
225 params=moldyn->pot_params;
226 exchange=&(params->exchange);
229 if(brand==ai->brand) {
234 r0=params->r0[brand];
235 mu=params->mu[brand];
236 lambda=params->lambda[brand];
245 lambda=params->lambda_m;
257 arg=M_PI*(d_ij-R)/s_r;
258 f_c=0.5+0.5*cos(arg);
259 df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
263 f_a=-B*exp(-mu*(d_ij-r0));
267 f_r=A*exp(-lambda*(d_ij-r0));
268 df_r=lambda*f_r/d_ij;
271 if(exchange->zeta_ij==0.0) {
276 b=1.0/sqrt(1.0+exchange->zeta_ij);
277 db=-0.5*b/(1.0+exchange->zeta_ij);
280 /* force contribution */
281 scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a));
282 v3_scale(&force,&(exchange->dist_ij),scale);
283 v3_add(&(ai->f),&(ai->f),&force);
284 v3_sub(&(aj->f),&(aj->f),&force); // dri rij = - drj rij
287 if(moldyn->time>DSTART&&moldyn->time<DEND) {
288 if((ai==&(moldyn->atom[5832]))|(aj==&(moldyn->atom[5832]))) {
289 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,aj->tag);
290 printf("adding %f %f %f\n",force.x,force.y,force.z);
291 if(ai==&(moldyn->atom[0]))
292 printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
293 if(aj==&(moldyn->atom[0]))
294 printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
295 printf("energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
297 printf(" %f %f %f\n",exchange->zeta_ij,.0,.0);
303 virial_calc(ai,&force,&(exchange->dist_ij));
305 /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
306 exchange->pre_dzeta=0.5*f_a*f_c*db;
308 /* energy contribution */
309 energy=0.5*f_c*(f_r-b*f_a);
310 moldyn->energy+=energy;
313 /* reset k counter for second k loop */
319 /* albe 3 body potential function (second k loop) */
320 int albe_mult_3bp_k2(t_moldyn *moldyn,
321 t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
323 t_albe_mult_params *params;
324 t_albe_exchange *exchange;
326 t_3dvec dist_ik,dist_ij;
327 double d_ik2,d_ik,d_ij2,d_ij;
330 double g,dg,cos_theta;
332 double f_c_ik,df_c_ik;
333 double dijdik_inv,fcdg,dfcg;
334 t_3dvec dcosdri,dcosdrj,dcosdrk;
337 params=moldyn->pot_params;
338 exchange=&(params->exchange);
339 kcount=exchange->kcount;
342 printf("FATAL: neighbours!\n");
345 d_ik2=exchange->d_ik2[kcount];
349 S2=params->S2[brand];
353 /* return if d_ik > S */
359 /* prefactor dzeta */
360 pre_dzeta=exchange->pre_dzeta;
363 dist_ik=exchange->dist_ik[kcount];
364 d_ik=exchange->d_ik[kcount];
366 /* f_c_ik, df_c_ik */
367 f_c_ik=exchange->f_c_ik[kcount];
368 df_c_ik=exchange->df_c_ik[kcount];
370 /* dist_ij, d_ij2, d_ij */
371 dist_ij=exchange->dist_ij;
372 d_ij2=exchange->d_ij2;
375 /* g, dg, cos_theta */
376 g=exchange->g[kcount];
377 dg=exchange->dg[kcount];
378 cos_theta=exchange->cos_theta[kcount];
380 /* cos_theta derivatives wrt i,j,k */
381 dijdik_inv=1.0/(d_ij*d_ik);
382 v3_scale(&dcosdrj,&dist_ik,dijdik_inv); // j
383 v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
384 v3_add(&dcosdrj,&dcosdrj,&tmp);
385 v3_scale(&dcosdrk,&dist_ij,dijdik_inv); // k
386 v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
387 v3_add(&dcosdrk,&dcosdrk,&tmp);
388 v3_add(&dcosdri,&dcosdrj,&dcosdrk); // i
389 v3_scale(&dcosdri,&dcosdri,-1.0);
391 /* f_c_ik * dg, df_c_ik * g */
395 /* derivative wrt i */
396 v3_scale(&force,&dist_ik,dfcg);
397 v3_scale(&tmp,&dcosdri,fcdg);
398 v3_add(&force,&force,&tmp);
399 v3_scale(&force,&force,pre_dzeta);
401 /* force contribution */
402 v3_add(&(ai->f),&(ai->f),&force);
405 if(moldyn->time>DSTART&&moldyn->time<DEND) {
406 if(ai==&(moldyn->atom[5832])) {
407 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
408 printf("adding %f %f %f\n",force.x,force.y,force.z);
409 printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
415 //virial_calc(ai,&force,&dist_ij);
417 /* derivative wrt j */
418 v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
420 /* force contribution */
421 v3_add(&(aj->f),&(aj->f),&force);
424 if(moldyn->time>DSTART&&moldyn->time<DEND) {
425 if(aj==&(moldyn->atom[5832])) {
426 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
427 printf("adding %f %f %f\n",force.x,force.y,force.z);
428 printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
434 v3_scale(&force,&force,-1.0);
435 virial_calc(ai,&force,&dist_ij);
437 /* derivative wrt k */
438 v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
439 v3_scale(&tmp,&dcosdrk,fcdg);
440 v3_add(&force,&force,&tmp);
441 v3_scale(&force,&force,pre_dzeta);
443 /* force contribution */
444 v3_add(&(ak->f),&(ak->f),&force);
447 if(moldyn->time>DSTART&&moldyn->time<DEND) {
448 if(ak==&(moldyn->atom[5832])) {
449 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
450 printf("adding %f %f %f\n",force.x,force.y,force.z);
451 printf("total k: %f %f %f\n",ak->f.x,ak->f.y,ak->f.z);
457 v3_scale(&force,&force,-1.0);
458 virial_calc(ai,&force,&dist_ik);
460 /* increase k counter */