potentials/albe.c

   1 /*
   2  * albe.c - albe potential
   3  *
   4  * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
   5  *
   6  */
   7
   8 #define _GNU_SOURCE
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <string.h>
  12 #include <sys/types.h>
  13 #include <sys/stat.h>
  14 #include <fcntl.h>
  15 #include <unistd.h>
  16 #include <math.h>
  17
  18 #include "../moldyn.h"
  19 #include "../math/math.h"
  20 #include "albe.h"
  21
  22 /* create mixed terms from parameters and set them */
  23 int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
  24
  25         t_albe_mult_params *p;
  26
  27         // set cutoff before parameters (actually only necessary for some pots)
  28         if(moldyn->cutoff==0.0) {
  29                 printf("[albe] WARNING: no cutoff!\n");
  30                 return -1;
  31         }
  32
  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");
  37                 return -1;
  38         }
  39
  40         /* these are now albe parameters */
  41         p=moldyn->pot_params;
  42
  43         // only 1 combination by now :p
  44         switch(element1) {
  45                 case SI:
  46                         /* type: silicon */
  47                         p->S[0]=ALBE_S_SI;
  48                         p->R[0]=ALBE_R_SI;
  49                         p->A[0]=ALBE_A_SI;
  50                         p->B[0]=ALBE_B_SI;
  51                         p->r0[0]=ALBE_R0_SI;
  52                         p->lambda[0]=ALBE_LAMBDA_SI;
  53                         p->mu[0]=ALBE_MU_SI;
  54                         p->gamma[0]=ALBE_GAMMA_SI;
  55                         p->c[0]=ALBE_C_SI;
  56                         p->c2[0]=p->c[0]*p->c[0];
  57                         p->d[0]=ALBE_D_SI;
  58                         p->d2[0]=p->d[0]*p->d[0];
  59                         p->c2d2[0]=p->c2[0]/p->d2[0];
  60                         p->h[0]=ALBE_H_SI;
  61                         switch(element2) {
  62                                 case C:
  63                                         /* type: carbon */
  64                                         p->S[1]=ALBE_S_C;
  65                                         p->R[1]=ALBE_R_C;
  66                                         p->A[1]=ALBE_A_C;
  67                                         p->B[1]=ALBE_B_C;
  68                                         p->r0[1]=ALBE_R0_C;
  69                                         p->lambda[1]=ALBE_LAMBDA_C;
  70                                         p->mu[1]=ALBE_MU_C;
  71                                         p->gamma[1]=ALBE_GAMMA_C;
  72                                         p->c[1]=ALBE_C_C;
  73                                         p->c2[1]=p->c[1]*p->c[1];
  74                                         p->d[1]=ALBE_D_C;
  75                                         p->d2[1]=p->d[1]*p->d[1];
  76                                         p->c2d2[1]=p->c2[1]/p->d2[1];
  77                                         p->h[1]=ALBE_H_C;
  78                                         /* mixed type: silicon carbide */
  79                                         p->Smixed=ALBE_S_SIC;
  80                                         p->Rmixed=ALBE_R_SIC;
  81                                         p->Amixed=ALBE_A_SIC;
  82                                         p->Bmixed=ALBE_B_SIC;
  83                                         p->r0_mixed=ALBE_R0_SIC;
  84                                         p->lambda_m=ALBE_LAMBDA_SIC;
  85                                         p->mu_m=ALBE_MU_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;
  93                                         break;
  94                                 default:
  95                                         printf("[albe] WARNING: element2\n");
  96                                         return -1;
  97                         }
  98                         break;
  99                 default:
 100                         printf("[albe] WARNING: element1\n");
 101                         return -1;
 102         }
 103
 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
 109         printf("[albe] mult parameter info:\n");
 110         printf("  S (A)  | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed);
 111         printf("  R (A)  | %f | %f | %f\n",p->R[0],p->R[1],p->Rmixed);
 112         printf("  A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV);
 113         printf("  B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV);
 114         printf("  lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1],
 115                                           p->lambda_m);
 116         printf("  mu     | %f | %f | %f\n",p->mu[0],p->mu[1],p->mu_m);
 117         printf("  gamma  | %f | %f | %f\n",p->gamma[0],p->gamma[1],p->gamma_m);
 118         printf("  c      | %f | %f | %f\n",p->c[0],p->c[1],p->c_mixed);
 119         printf("  d      | %f | %f | %f\n",p->d[0],p->d[1],p->d_mixed);
 120         printf("  c2     | %f | %f | %f\n",p->c2[0],p->c2[1],p->c2_mixed);
 121         printf("  d2     | %f | %f | %f\n",p->d2[0],p->d2[1],p->d2_mixed);
 122         printf("  c2d2   | %f | %f | %f\n",p->c2d2[0],p->c2d2[1],p->c2d2_m);
 123         printf("  h      | %f | %f | %f\n",p->h[0],p->h[1],p->h_mixed);
 124
 125         return 0;
 126 }
 127
 128 /* first i loop */
 129 int albe_mult_i0(t_moldyn *moldyn,t_atom *ai) {
 130
 131         t_albe_mult_params *params;
 132         t_albe_exchange *exchange;
 133
 134         int i;
 135
 136         params=moldyn->pot_params;
 137         exchange=&(params->exchange);
 138
 139         /* zero exchange values */
 140         memset(exchange->zeta,0,ALBE_MAXN*sizeof(double));
 141         for(i=0;i<ALBE_MAXN;i++)
 142                 memset(exchange->dzeta[i],0,ALBE_MAXN*sizeof(t_3dvec));
 143         exchange->jcnt=0;
 144         exchange->j2cnt=0;
 145
 146         return 0;
 147 }
 148
 149 /* first j loop within first i loop */
 150 int albe_mult_i0_j0(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 151
 152         t_albe_mult_params *params;
 153         t_albe_exchange *exchange;
 154
 155         double S2,S,R,d2,d,s_r,arg;
 156         t_3dvec dist;
 157         int j;
 158         u8 brand;
 159
 160         params=moldyn->pot_params;
 161         exchange=&(params->exchange);
 162
 163         /* get j counter */
 164         j=exchange->jcnt;
 165
 166         /* set ij depending values */
 167         brand=ai->brand;
 168         if(brand==aj->brand) {
 169                 S2=params->S2[brand];
 170         }
 171         else {
 172                 S2=params->S2mixed;
 173         }
 174
 175         /* dist_ij, d_ij2 */
 176         v3_sub(&dist,&(aj->r),&(ai->r));
 177         if(bc) check_per_bound(moldyn,&dist);
 178         exchange->dist[j]=dist;
 179         d2=v3_absolute_square(&dist);
 180         exchange->d2[j]=d2;
 181
 182         /* if d_ij2 > S2 => no force & potential energy contribution */
 183         if(d2>S2) {
 184                 moldyn->run3bp=0;
 185                 exchange->skip[j]=1;
 186                 exchange->jcnt+=1;
 187                 return 0;
 188         }
 189         exchange->skip[j]=0;
 190
 191         /* more ij depending values */
 192         if(brand==aj->brand) {
 193                 R=params->R[brand];
 194                 S=params->S[brand];
 195                 /* albe needs i,(j/k) depending c,d,h and gamma values */
 196                 exchange->gamma_[j]=&(params->gamma[brand]);
 197                 exchange->c_[j]=&(params->c[brand]);
 198                 exchange->d_[j]=&(params->d[brand]);
 199                 exchange->h_[j]=&(params->h[brand]);
 200                 exchange->c2_[j]=&(params->c2[brand]);
 201                 exchange->d2_[j]=&(params->d2[brand]);
 202                 exchange->c2d2_[j]=&(params->c2d2[brand]);
 203         }
 204         else {
 205                 R=params->Rmixed;
 206                 S=params->Smixed;
 207                 /* albe needs i,(j/k) depending c,d,h and gamma values */
 208                 exchange->gamma_[j]=&(params->gamma_m);
 209                 exchange->c_[j]=&(params->c_mixed);
 210                 exchange->d_[j]=&(params->d_mixed);
 211                 exchange->h_[j]=&(params->h_mixed);
 212                 exchange->c2_[j]=&(params->c2_mixed);
 213                 exchange->d2_[j]=&(params->d2_mixed);
 214                 exchange->c2d2_[j]=&(params->c2d2_m);
 215         }
 216
 217         /* d_ij */
 218         d=sqrt(d2);
 219         exchange->d[j]=d;
 220
 221         /* f_c, df_c */
 222         if(d<R) {
 223                 exchange->f_c[j]=1.0;
 224                 exchange->df_c[j]=0.0;
 225         }
 226         else {
 227                 s_r=S-R;
 228                 arg=M_PI*(d-R)/s_r;
 229                 exchange->f_c[j]=0.5+0.5*cos(arg);
 230                 exchange->df_c[j]=0.5*sin(arg)*(M_PI/(s_r*d));
 231         }
 232
 233         /* reset k counter */
 234         exchange->kcnt=0;
 235
 236         return 0;
 237 }
 238
 239 /* first k loop within first j loop within first i loop */
 240 int albe_mult_i0_j0_k0(t_moldyn *moldyn,
 241                        t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 242
 243         t_albe_mult_params *params;
 244         t_albe_exchange *exchange;
 245
 246         int j,k;
 247         t_3dvec distj,distk;
 248         double dj,dk,djdk_inv,cos_theta;
 249         double gj,dgj,h_cos_j,d2_h_cos2_j,frac_j;
 250         double gk,dgk,h_cos_k,d2_h_cos2_k,frac_k;
 251         t_3dvec dcosdrj,dcosdrk,tmp;
 252         t_3dvec *dzjj,*dzkk,*dzjk,*dzkj;
 253
 254         params=moldyn->pot_params;
 255         exchange=&(params->exchange);
 256
 257         if(aj==ak) {
 258                 exchange->kcnt+=1;
 259                 return 0;
 260         }
 261
 262         /* k<j & check whether to run k */
 263         j=exchange->jcnt;
 264         k=exchange->kcnt;
 265         if(k>=ALBE_MAXN) {
 266                 printf("FATAL: too many neighbours! (%d)\n",k);
 267                 printf("  atom i:%d | j:%d | k:%d\n",ai->tag,aj->tag,ak->tag);
 268         }
 269         if((k>=j)|(exchange->skip[k])) {
 270                 exchange->kcnt+=1;
 271                 return 0;
 272         }
 273
 274         /* distances */
 275         distj=exchange->dist[j];
 276         distk=exchange->dist[k];
 277         dj=exchange->d[j];
 278         dk=exchange->d[k];
 279         djdk_inv=1.0/(dj*dk);
 280
 281         /* cos theta */
 282         cos_theta=v3_scalar_product(&distj,&distk)*djdk_inv;
 283
 284         /* g(cos(theta)) ij and ik values */
 285         h_cos_j=*(exchange->h_[j])+cos_theta; // + in albe formalism
 286         d2_h_cos2_j=*exchange->d2_[j]+(h_cos_j*h_cos_j);
 287         frac_j=*exchange->c2_[j]/d2_h_cos2_j;
 288         gj=1.0+*exchange->c2d2_[j]-frac_j;
 289         gj*=*(exchange->gamma_[j]);
 290         dgj=*(exchange->gamma_[j])*2.0*frac_j*h_cos_j/d2_h_cos2_j; // + in albe
 291         if(ak->brand==aj->brand) {
 292                 gk=gj;
 293                 dgk=dgj;
 294         }
 295         else {
 296                 h_cos_k=*(exchange->h_[k])+cos_theta;
 297                 d2_h_cos2_k=*exchange->d2_[k]+(h_cos_k*h_cos_k);
 298                 frac_k=*exchange->c2_[k]/d2_h_cos2_k;
 299                 gk=1.0+*exchange->c2d2_[k]-frac_k;
 300                 gk*=*(exchange->gamma_[k]);
 301                 dgk=*(exchange->gamma_[k])*2.0*frac_k*h_cos_k/d2_h_cos2_k;
 302         }
 303
 304         /* zeta - for albe: ik depending g function */
 305 //if(ai->tag==0) {
 306 //      printf("------> %.15f %.15f\n",dj,dk);
 307 //      printf("------> %.15f %.15f\n",dj,dk);
 308 //}
 309
 310         exchange->zeta[j]+=(exchange->f_c[k]*gk);
 311         exchange->zeta[k]+=(exchange->f_c[j]*gj);
 312
 313         /* cos theta derivatives */
 314         v3_scale(&dcosdrj,&distk,djdk_inv);             // j
 315         v3_scale(&tmp,&distj,-cos_theta/exchange->d2[j]);
 316         v3_add(&dcosdrj,&dcosdrj,&tmp);
 317         v3_scale(&dcosdrk,&distj,djdk_inv);             // k
 318         v3_scale(&tmp,&distk,-cos_theta/exchange->d2[k]);
 319         v3_add(&dcosdrk,&dcosdrk,&tmp);
 320
 321         /* zeta derivatives */
 322         dzjj=&(exchange->dzeta[j][j]);
 323         dzkk=&(exchange->dzeta[k][k]);
 324         dzjk=&(exchange->dzeta[j][k]);
 325         dzkj=&(exchange->dzeta[k][j]);
 326         v3_scale(&tmp,&dcosdrj,exchange->f_c[k]*dgk);
 327         v3_add(dzjj,dzjj,&tmp);                         // j j
 328         v3_scale(&tmp,&dcosdrk,exchange->f_c[j]*dgj);
 329         v3_add(dzkk,dzkk,&tmp);                         // k k
 330         v3_scale(&tmp,&distk,-exchange->df_c[k]*gk);    // dk rik = - di rik
 331         v3_add(dzjk,dzjk,&tmp);
 332         v3_scale(&tmp,&dcosdrk,exchange->f_c[k]*dgk);
 333         v3_add(dzjk,dzjk,&tmp);                         // j k
 334         v3_scale(&tmp,&distj,-exchange->df_c[j]*gj);    // dj rij = - di rij
 335         v3_add(dzkj,dzkj,&tmp);
 336         v3_scale(&tmp,&dcosdrj,exchange->f_c[j]*dgj);
 337         v3_add(dzkj,dzkj,&tmp);                         // k j
 338
 339         /* increase k counter */
 340         exchange->kcnt+=1;
 341
 342         return 0;
 343 }
 344
 345 /* first j loop within first i loop */
 346 int albe_mult_i0_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 347
 348         t_albe_mult_params *params;
 349         t_albe_exchange *exchange;
 350
 351         params=moldyn->pot_params;
 352         exchange=&(params->exchange);
 353
 354         /* increase j counter */
 355         exchange->jcnt+=1;
 356
 357         return 0;
 358 }
 359
 360 /* second j loop within first i loop */
 361 int albe_mult_i0_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 362
 363         t_albe_mult_params *params;
 364         t_albe_exchange *exchange;
 365
 366         int j;
 367         double d,f_a,df_a,f_r,df_r,f_c,df_c,b,db;
 368         double A,B,mu,lambda,r0;
 369         double energy;
 370         t_3dvec *dist,force;
 371         double scale;
 372         u8 brand;
 373
 374         params=moldyn->pot_params;
 375         exchange=&(params->exchange);
 376
 377         /* get j counter */
 378         j=exchange->j2cnt;
 379
 380         /* skip if j not within cutoff */
 381         if(exchange->skip[j]) {
 382                 moldyn->run3bp=0;
 383                 exchange->j2cnt+=1;
 384                 return 0;
 385         }
 386
 387         /* distance */
 388         d=exchange->d[j];
 389         dist=&(exchange->dist[j]);
 390         f_c=exchange->f_c[j];
 391         df_c=exchange->df_c[j];
 392
 393         /* determine parameters to calculate fa, dfa, fr, dfr */
 394         brand=aj->brand;
 395         if(brand==ai->brand) {
 396                 B=params->B[brand];
 397                 A=params->A[brand];
 398                 r0=params->r0[brand];
 399                 mu=params->mu[brand];
 400                 lambda=params->lambda[brand];
 401         }
 402         else {
 403                 B=params->Bmixed;
 404                 A=params->Amixed;
 405                 r0=params->r0_mixed;
 406                 mu=params->mu_m;
 407                 lambda=params->lambda_m;
 408         }
 409
 410         /* f_a, df_a */
 411         f_a=-B*exp(-mu*(d-r0));
 412         df_a=mu*f_a/d;
 413
 414         /* f_r, df_r */
 415         f_r=A*exp(-lambda*(d-r0));
 416         df_r=lambda*f_r/d;
 417
 418         /* b, db */
 419         b=1.0/sqrt(1.0+exchange->zeta[j]);
 420         db=-0.5*b/(1.0+exchange->zeta[j]);
 421
 422         /* energy contribution */
 423         energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
 424         moldyn->energy+=energy;
 425         ai->e+=energy;
 426
 427         /* force contribution for atom i due to ij bond */
 428         scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
 429         v3_scale(&force,dist,scale);
 430         v3_add(&(ai->f),&(ai->f),&force);
 431
 432 #ifdef NDEBUG
 433 if(ai->tag==0) {
 434 printf("force: %.15f %.15f %.15f | %d %d (ij) %.15f\n",force.x,force.y,force.z,ai->tag,aj->tag,exchange->zeta[j]);
 435 printf("    t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z);
 436 }
 437 #endif
 438
 439         /* force contribution for atom j due to ij bond */
 440         v3_scale(&force,&force,-1.0); // dri rij = - drj rij
 441         v3_add(&(aj->f),&(aj->f),&force);
 442
 443 #ifdef NDEBUG
 444 if(aj->tag==0) {
 445 printf("force: %.15f %.15f %.15f | %d %d (ji) %.15f\n",force.x,force.y,force.z,aj->tag,ai->tag,exchange->zeta[j]);
 446 printf("    t: %.15f %.15f %.15f\n",aj->f.x,aj->f.y,aj->f.z);
 447 }
 448 #endif
 449
 450         /* virial */
 451         virial_calc(aj,&force,dist);
 452
 453         /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
 454         exchange->pre_dzeta=0.5*f_a*f_c*db;
 455
 456         /* force contribution (drj derivative) */
 457         v3_scale(&force,&(exchange->dzeta[j][j]),exchange->pre_dzeta);
 458         v3_add(&(aj->f),&(aj->f),&force);
 459
 460 #ifdef NDEBUG
 461 if(aj->tag==0) {
 462 printf("force: %.15f %.15f %.15f | %d %d (j der)\n",force.x,force.y,force.z,aj->tag,ai->tag);
 463 printf("    t: %.15f %.15f %.15f\n",aj->f.x,aj->f.y,aj->f.z);
 464 }
 465 #endif
 466
 467         v3_scale(&force,&force,-1.0);
 468         v3_add(&(ai->f),&(ai->f),&force);
 469
 470 #ifdef NDEBUG
 471 if(ai->tag==0) {
 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);
 474 }
 475 #endif
 476
 477         /* virial */
 478         virial_calc(ai,&force,dist);
 479
 480         /* reset k counter for second k loop */
 481         exchange->kcnt=0;
 482
 483         return 0;
 484 }
 485
 486 /* second k loop within second j loop within first i loop */
 487 int albe_mult_i0_j2_k0(t_moldyn *moldyn,
 488                        t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 489
 490         t_albe_mult_params *params;
 491         t_albe_exchange *exchange;
 492
 493         int j,k;
 494         t_3dvec force;
 495
 496         params=moldyn->pot_params;
 497         exchange=&(params->exchange);
 498
 499         if(aj==ak) {
 500                 exchange->kcnt+=1;
 501                 return 0;
 502         }
 503
 504         /* k!=j & check whether to run k */
 505         k=exchange->kcnt;
 506         j=exchange->j2cnt;
 507         if((k==j)|(exchange->skip[k])) {
 508                 exchange->kcnt+=1;
 509                 return 0;
 510         }
 511
 512         /* force contribution (drk derivative) */
 513         v3_scale(&force,&(exchange->dzeta[j][k]),exchange->pre_dzeta);
 514         v3_add(&(ak->f),&(ak->f),&force);
 515
 516 #ifdef NDEBUG
 517 if(ak->tag==0) {
 518 printf("force: %.15f %.15f %.15f | %d %d %d (k der)\n",force.x,force.y,force.z,ai->tag,aj->tag,ak->tag);
 519 printf("    t: %.15f %.15f %.15f\n",ak->f.x,ak->f.y,ak->f.z);
 520 }
 521 #endif
 522
 523         v3_scale(&force,&force,-1.0);
 524         v3_add(&(ai->f),&(ai->f),&force);
 525
 526 #ifdef NDEBUG
 527 if(ai->tag==0) {
 528 printf("force: %.15f %.15f %.15f | %d %d %d -- %d(i contr k der)\n",force.x,force.y,force.z,ai->tag,aj->tag,ak->tag,k);
 529 printf("    t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z);
 530 printf("  ## %f\n",exchange->d[k]);
 531 }
 532 #endif
 533
 534         /* virial */
 535         virial_calc(ai,&force,&(exchange->dist[k]));
 536
 537         /* increase k counter */
 538         exchange->kcnt+=1;
 539
 540         return 0;
 541 }
 542
 543 int albe_mult_i0_j3(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 544
 545         t_albe_mult_params *params;
 546         t_albe_exchange *exchange;
 547
 548         params=moldyn->pot_params;
 549         exchange=&(params->exchange);
 550
 551         /* increase j counter */
 552         exchange->j2cnt+=1;
 553
 554         return 0;
 555 }
 556
 557 int albe_mult_check_2b_bond(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,u8 bc) {
 558
 559         t_albe_mult_params *params;
 560         t_3dvec dist;
 561         double d;
 562         u8 brand;
 563
 564         v3_sub(&dist,&(jtom->r),&(itom->r));
 565         if(bc) check_per_bound(moldyn,&dist);
 566         d=v3_absolute_square(&dist);
 567
 568         params=moldyn->pot_params;
 569         brand=itom->brand;
 570
 571         if(brand==jtom->brand) {
 572                 if(d<=params->S2[brand])
 573                         return TRUE;
 574         }
 575         else {
 576                 if(d<=params->S2mixed)
 577                         return TRUE;
 578         }
 579
 580         return FALSE;
 581 }