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->d[0]=ALBE_D_SI;
  57                         p->h[0]=ALBE_H_SI;
  58                         switch(element2) {
  59                                 case C:
  60                                         /* type: carbon */
  61                                         p->S[1]=ALBE_S_C;
  62                                         p->R[1]=ALBE_R_C;
  63                                         p->A[1]=ALBE_A_C;
  64                                         p->B[1]=ALBE_B_C;
  65                                         p->r0[1]=ALBE_R0_C;
  66                                         p->lambda[1]=ALBE_LAMBDA_C;
  67                                         p->mu[1]=ALBE_MU_C;
  68                                         p->gamma[1]=ALBE_GAMMA_C;
  69                                         p->c[1]=ALBE_C_C;
  70                                         p->d[1]=ALBE_D_C;
  71                                         p->h[1]=ALBE_H_C;
  72                                         /* mixed type: silicon carbide */
  73                                         p->Smixed=ALBE_S_SIC;
  74                                         p->Rmixed=ALBE_R_SIC;
  75                                         p->Amixed=ALBE_A_SIC;
  76                                         p->Bmixed=ALBE_B_SIC;
  77                                         p->r0_mixed=ALBE_R0_SIC;
  78                                         p->lambda_m=ALBE_LAMBDA_SIC;
  79                                         p->mu_m=ALBE_MU_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;
  84                                         break;
  85                                 default:
  86                                         printf("[albe] WARNING: element2\n");
  87                                         return -1;
  88                         }
  89                         break;
  90                 default:
  91                         printf("[albe] WARNING: element1\n");
  92                         return -1;
  93         }
  94
  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;
  99
 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],
 106                                           p->lambda_m);
 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]);
 112
 113         return 0;
 114 }
 115
 116 /* first i loop */
 117 int albe_mult_i0(t_moldyn *moldyn,t_atom *ai) {
 118
 119         int i;
 120         t_albe_mult_params *params;
 121         t_albe_exchange *exchange;
 122
 123         params=moldyn->pot_params;
 124         exchange=&(params->exchange);
 125
 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));
 134         exchange->jcnt=0;
 135         exchange->j2cnt=0;
 136
 137         return 0;
 138 }
 139
 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) {
 142
 143         t_albe_mult_params *params;
 144         t_albe_exchange *exchange;
 145         double S2,S,R,d2,d,s_r,arg;
 146         t_3dvec dist;
 147         int j;
 148         u8 brand;
 149
 150         params=moldyn->pot_params;
 151         exchange=&(params->exchange);
 152
 153         /* get j counter */
 154         j=exchange->jcnt;
 155
 156         /* set ij depending values */
 157         brand=ai->brand;
 158         if(brand==aj->brand) {
 159                 S2=params->S2[brand];
 160         }
 161         else {
 162                 S2=params->S2mixed;
 163         }
 164
 165         /* dist_ij, d_ij2 */
 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);
 170         exchange->d2[j]=d2;
 171
 172         /* if d_ij2 > S2 => no force & potential energy contribution */
 173         if(d2>S2) {
 174                 moldyn->run3bp=0;
 175                 exchange->skip[j]=1;
 176                 exchange->jcnt+=1;
 177                 return 0;
 178         }
 179
 180         /* more ij depending values */
 181         if(brand==aj->brand) {
 182                 R=params->R[brand];
 183                 S=params->S[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]);
 189         }
 190         else {
 191                 R=params->Rmixed;
 192                 S=params->Smixed;
 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);
 198         }
 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];
 202
 203         /* d_ij */
 204         d=sqrt(exchange->d2[j]);
 205         exchange->d[j]=d;
 206
 207         /* f_c, df_c */
 208         if(d<R) {
 209                 exchange->f_c[j]=1.0;
 210                 exchange->df_c[j]=0.0;
 211         }
 212         else {
 213                 s_r=S-R;
 214                 arg=M_PI*(d-R)/s_r;
 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));
 217         }
 218
 219         /* reset k counter */
 220         exchange->kcnt=0;
 221
 222         return 0;
 223 }
 224
 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) {
 228
 229         t_albe_mult_params *params;
 230         t_albe_exchange *exchange;
 231         int j,k;
 232         t_3dvec distj,distk;
 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;
 237
 238         params=moldyn->pot_params;
 239         exchange=&(params->exchange);
 240
 241         /* k<j & check whether to run k */
 242         j=exchange->jcnt;
 243         k=exchange->kcnt;
 244         if(k>=ALBE_MAXN) {
 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);
 247         }
 248         if((k>=j)|(exchange->skip[k])) {
 249                 exchange->kcnt+=1;
 250                 return 0;
 251         }
 252
 253         /* distances */
 254         distj=exchange->dist[j];
 255         distk=exchange->dist[k];
 256         dj=exchange->d[j];
 257         dk=exchange->d[k];
 258         djdk_inv=1.0/(dj*dk);
 259
 260         /* cos theta */
 261         cos_theta=v3_scalar_product(&distj,&distk)*djdk_inv;
 262
 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) {
 271                 gk=gj;
 272                 dgk=dgj;
 273         }
 274         else {
 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;
 281         }
 282
 283         /* zeta */
 284         exchange->zeta[j]+=(exchange->f_c[k]*gk);
 285         exchange->zeta[k]+=(exchange->f_c[j]*gj);
 286
 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);
 294
 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
 309
 310         /* increase k counter */
 311         exchange->kcnt+=1;
 312
 313         return 0;
 314 }
 315
 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) {
 318
 319         t_albe_mult_params *params;
 320         t_albe_exchange *exchange;
 321
 322         params=moldyn->pot_params;
 323         exchange=&(params->exchange);
 324
 325         /* increase j counter */
 326         exchange->jcnt+=1;
 327
 328         return 0;
 329 }
 330
 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) {
 333
 334         t_albe_mult_params *params;
 335         t_albe_exchange *exchange;
 336
 337         int j;
 338         double d,f_a,df_a,f_r,df_r,f_c,df_c,b,db;
 339         double A,B,mu,lambda,r0;
 340         double energy;
 341         t_3dvec *dist,force;
 342         double scale;
 343         u8 brand;
 344
 345         params=moldyn->pot_params;
 346         exchange=&(params->exchange);
 347
 348         /* get j counter */
 349         j=exchange->j2cnt;
 350
 351         /* skip if j not within cutoff */
 352         if(exchange->skip[j]) {
 353                 moldyn->run3bp=0;
 354                 exchange->j2cnt+=1;
 355                 return 0;
 356         }
 357
 358         /* distance */
 359         d=exchange->d[j];
 360         dist=&(exchange->dist[j]);
 361         f_c=exchange->f_c[j];
 362         df_c=exchange->df_c[j];
 363
 364         /* determine parameters to calculate fa, dfa, fr, dfr */
 365         brand=aj->brand;
 366         if(brand==ai->brand) {
 367                 B=params->B[brand];
 368                 A=params->A[brand];
 369                 r0=params->r0[brand];
 370                 mu=params->mu[brand];
 371                 lambda=params->lambda[brand];
 372         }
 373         else {
 374                 B=params->Bmixed;
 375                 A=params->Amixed;
 376                 r0=params->r0_mixed;
 377                 mu=params->mu_m;
 378                 lambda=params->lambda_m;
 379         }
 380
 381         /* f_a, df_a */
 382         f_a=-B*exp(-mu*(d-r0));
 383         df_a=mu*f_a/d;
 384
 385         /* f_r, df_r */
 386         f_r=A*exp(-lambda*(d-r0));
 387         df_r=lambda*f_r/d;
 388
 389         /* b, db */
 390         b=1.0/sqrt(1.0+exchange->zeta[j]);
 391         db=-0.5*b/(1.0+exchange->zeta[j]);
 392
 393         /* energy contribution */
 394         energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
 395         moldyn->energy+=energy;
 396         ai->e+=energy;
 397
 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);
 402
 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);
 406
 407         /* virial */
 408         virial_calc(aj,&force,dist);
 409
 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;
 412
 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);
 418
 419         /* virial */
 420         virial_calc(ai,&force,dist);
 421
 422         /* reset k counter for second k loop */
 423         exchange->kcnt=0;
 424
 425         return 0;
 426 }
 427
 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) {
 431
 432         t_albe_mult_params *params;
 433         t_albe_exchange *exchange;
 434
 435         int j,k;
 436         t_3dvec force;
 437
 438         params=moldyn->pot_params;
 439         exchange=&(params->exchange);
 440
 441         /* k!=j & check whether to run k */
 442         k=exchange->kcnt;
 443         j=exchange->j2cnt;
 444         if((k==j)|(exchange->skip[k])) {
 445                 exchange->kcnt+=1;
 446                 return 0;
 447         }
 448
 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);
 454
 455         /* virial */
 456         virial_calc(ai,&force,&(exchange->dist[k]));
 457
 458         /* increase k counter */
 459         exchange->kcnt+=1;
 460
 461         return 0;
 462 }
 463
 464 int albe_mult_i0_j3(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 465
 466         t_albe_mult_params *params;
 467         t_albe_exchange *exchange;
 468
 469         params=moldyn->pot_params;
 470         exchange=&(params->exchange);
 471
 472         /* increase j counter */
 473         exchange->j2cnt+=1;
 474
 475         return 0;
 476 }
 477
 478 int albe_mult_check_2b_bond(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,u8 bc) {
 479
 480         t_albe_mult_params *params;
 481         t_3dvec dist;
 482         double d;
 483         u8 brand;
 484
 485         v3_sub(&dist,&(jtom->r),&(itom->r));
 486         if(bc) check_per_bound(moldyn,&dist);
 487         d=v3_absolute_square(&dist);
 488
 489         params=moldyn->pot_params;
 490         brand=itom->brand;
 491
 492         if(brand==jtom->brand) {
 493                 if(d<=params->S2[brand])
 494                         return TRUE;
 495         }
 496         else {
 497                 if(d<=params->S2mixed)
 498                         return TRUE;
 499         }
 500
 501         return FALSE;
 502 }