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         p->c2[0]=p->c[0]*p->c[0];
 100         p->c2[1]=p->c[1]*p->c[1];
 101         p->c2_mixed=p->c_mixed*p->c_mixed;
 102         p->d2[0]=p->d[0]*p->d[0];
 103         p->d2[1]=p->d[1]*p->d[1];
 104         p->d2_mixed=p->d_mixed*p->d_mixed;
 105         p->c2d2[0]=p->c2[0]/p->d2[0];
 106         p->c2d2[1]=p->c2[1]/p->d2[1];
 107         p->c2d2_m=p->c2_mixed/p->d2_mixed;
 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\n",p->gamma[0],p->gamma[1]);
 118         printf("  c      | %f | %f\n",p->c[0],p->c[1]);
 119         printf("  d      | %f | %f\n",p->d[0],p->d[1]);
 120         printf("  h      | %f | %f\n",p->h[0],p->h[1]);
 121
 122         return 0;
 123 }
 124
 125 /* albe 3 body potential function (first ij loop) */
 126 int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 127
 128         t_albe_mult_params *params;
 129         t_albe_exchange *exchange;
 130         unsigned char brand;
 131         double S2;
 132         t_3dvec dist_ij;
 133         double d_ij2,d_ij;
 134
 135         params=moldyn->pot_params;
 136         exchange=&(params->exchange);
 137
 138         /* reset zeta sum */
 139         exchange->zeta_ij=0.0;
 140
 141         /*
 142          * set ij depending values
 143          */
 144
 145         brand=ai->brand;
 146         if(brand==aj->brand) {
 147                 S2=params->S2[brand];
 148         }
 149         else {
 150                 S2=params->S2mixed;
 151         }
 152
 153         /* dist_ij, d_ij2 */
 154         v3_sub(&dist_ij,&(aj->r),&(ai->r));
 155         if(bc) check_per_bound(moldyn,&dist_ij);
 156         d_ij2=v3_absolute_square(&dist_ij);
 157
 158         /* if d_ij2 > S2 => no force & potential energy contribution */
 159         if(d_ij2>S2) {
 160                 moldyn->run3bp=0;
 161                 return 0;
 162         }
 163
 164         /* d_ij */
 165         d_ij=sqrt(d_ij2);
 166
 167         /* store values */
 168         exchange->dist_ij=dist_ij;
 169         exchange->d_ij2=d_ij2;
 170         exchange->d_ij=d_ij;
 171
 172         /* reset k counter for first k loop */
 173         exchange->kcount=0;
 174
 175         return 0;
 176 }
 177
 178 /* albe 3 body potential function (first k loop) */
 179 int albe_mult_3bp_k1(t_moldyn *moldyn,
 180                      t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 181
 182         t_albe_mult_params *params;
 183         t_albe_exchange *exchange;
 184         unsigned char brand;
 185         double R,S,S2;
 186         t_3dvec dist_ij,dist_ik;
 187         double d_ik2,d_ik,d_ij;
 188         double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
 189         double f_c_ik,df_c_ik;
 190         int kcount;
 191
 192         params=moldyn->pot_params;
 193         exchange=&(params->exchange);
 194         kcount=exchange->kcount;
 195
 196         if(kcount>ALBE_MAXN) {
 197                 printf("FATAL: neighbours = %d\n",kcount);
 198                 printf("  -> %d %d %d\n",ai->tag,aj->tag,ak->tag);
 199         }
 200
 201         /* ik constants */
 202         brand=ai->brand;
 203         if(brand==ak->brand) {
 204                 R=params->R[brand];
 205                 S=params->S[brand];
 206                 S2=params->S2[brand];
 207                 /* albe needs i,k depending c,d,h and gamma values */
 208                 exchange->gamma_i=&(params->gamma[brand]);
 209                 exchange->c_i=&(params->c[brand]);
 210                 exchange->d_i=&(params->d[brand]);
 211                 exchange->h_i=&(params->h[brand]);
 212         }
 213         else {
 214                 R=params->Rmixed;
 215                 S=params->Smixed;
 216                 S2=params->S2mixed;
 217                 /* albe needs i,k depending c,d,h and gamma values */
 218                 exchange->gamma_i=&(params->gamma_m);
 219                 exchange->c_i=&(params->c_mixed);
 220                 exchange->d_i=&(params->d_mixed);
 221                 exchange->h_i=&(params->h_mixed);
 222         }
 223         exchange->ci2=*(exchange->c_i)**(exchange->c_i);
 224         exchange->di2=*(exchange->d_i)**(exchange->d_i);
 225         exchange->ci2di2=exchange->ci2/exchange->di2;
 226
 227         /* dist_ik, d_ik2 */
 228         v3_sub(&dist_ik,&(ak->r),&(ai->r));
 229         if(bc) check_per_bound(moldyn,&dist_ik);
 230         d_ik2=v3_absolute_square(&dist_ik);
 231
 232         /* store data for second k loop */
 233         exchange->dist_ik[kcount]=dist_ik;
 234         exchange->d_ik2[kcount]=d_ik2;
 235
 236         /* return if not within cutoff */
 237         if(d_ik2>S2) {
 238                 exchange->kcount++;
 239                 return 0;
 240         }
 241
 242         /* d_ik */
 243         d_ik=sqrt(d_ik2);
 244
 245         /* dist_ij, d_ij */
 246         dist_ij=exchange->dist_ij;
 247         d_ij=exchange->d_ij;
 248
 249         /* cos theta */
 250         cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
 251
 252         /* g_ijk */
 253         h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism
 254         d2_h_cos2=exchange->di2+(h_cos*h_cos);
 255         frac=exchange->ci2/d2_h_cos2;
 256         g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
 257         dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
 258
 259         /* zeta sum += f_c_ik * g_ijk */
 260         if(d_ik<=R) {
 261                 exchange->zeta_ij+=g;
 262                 f_c_ik=1.0;
 263                 df_c_ik=0.0;
 264         }
 265         else {
 266                 s_r=S-R;
 267                 arg=M_PI*(d_ik-R)/s_r;
 268                 f_c_ik=0.5+0.5*cos(arg);
 269                 df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
 270                 exchange->zeta_ij+=f_c_ik*g;
 271         }
 272
 273         /* store even more data for second k loop */
 274         exchange->g[kcount]=g;
 275         exchange->dg[kcount]=dg;
 276         exchange->d_ik[kcount]=d_ik;
 277         exchange->cos_theta[kcount]=cos_theta;
 278         exchange->f_c_ik[kcount]=f_c_ik;
 279         exchange->df_c_ik[kcount]=df_c_ik;
 280
 281         /* increase k counter */
 282         exchange->kcount++;
 283
 284         return 0;
 285 }
 286
 287 int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 288
 289         t_albe_mult_params *params;
 290         t_albe_exchange *exchange;
 291         t_3dvec force;
 292         double f_a,df_a,b,db,f_c,df_c;
 293         double f_r,df_r;
 294         double scale;
 295         double mu,B;
 296         double lambda,A;
 297         double d_ij,r0;
 298         unsigned char brand;
 299         double S,R,s_r,arg;
 300         double energy;
 301
 302         params=moldyn->pot_params;
 303         exchange=&(params->exchange);
 304
 305         brand=aj->brand;
 306         if(brand==ai->brand) {
 307                 S=params->S[brand];
 308                 R=params->R[brand];
 309                 B=params->B[brand];
 310                 A=params->A[brand];
 311                 r0=params->r0[brand];
 312                 mu=params->mu[brand];
 313                 lambda=params->lambda[brand];
 314         }
 315         else {
 316                 S=params->Smixed;
 317                 R=params->Rmixed;
 318                 B=params->Bmixed;
 319                 A=params->Amixed;
 320                 r0=params->r0_mixed;
 321                 mu=params->mu_m;
 322                 lambda=params->lambda_m;
 323         }
 324
 325         d_ij=exchange->d_ij;
 326
 327         /* f_c, df_c */
 328         if(d_ij<R) {
 329                 f_c=1.0;
 330                 df_c=0.0;
 331         }
 332         else {
 333                 s_r=S-R;
 334                 arg=M_PI*(d_ij-R)/s_r;
 335                 f_c=0.5+0.5*cos(arg);
 336                 df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
 337         }
 338
 339         /* f_a, df_a */
 340         f_a=-B*exp(-mu*(d_ij-r0));
 341         df_a=mu*f_a/d_ij;
 342
 343         /* f_r, df_r */
 344         f_r=A*exp(-lambda*(d_ij-r0));
 345         df_r=lambda*f_r/d_ij;
 346
 347         /* b, db */
 348         if(exchange->zeta_ij==0.0) {
 349                 b=1.0;
 350                 db=0.0;
 351         }
 352         else {
 353                 b=1.0/sqrt(1.0+exchange->zeta_ij);
 354                 db=-0.5*b/(1.0+exchange->zeta_ij);
 355         }
 356
 357         /* force contribution for atom i */
 358         scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
 359         v3_scale(&force,&(exchange->dist_ij),scale);
 360         v3_add(&(ai->f),&(ai->f),&force);
 361
 362         /* force contribution for atom j */
 363         v3_scale(&force,&force,-1.0); // dri rij = - drj rij
 364         v3_add(&(aj->f),&(aj->f),&force);
 365
 366         /* virial */
 367         virial_calc(ai,&force,&(exchange->dist_ij));
 368
 369 #ifdef DEBUG
 370 if(moldyn->time>DSTART&&moldyn->time<DEND) {
 371         if((ai==&(moldyn->atom[DATOM]))|(aj==&(moldyn->atom[DATOM]))) {
 372                 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,aj->tag);
 373                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
 374                 if(ai==&(moldyn->atom[0]))
 375                         printf("  total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
 376                 if(aj==&(moldyn->atom[0]))
 377                         printf("  total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
 378                 printf("  energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
 379                                                     f_c,b,f_a,f_r);
 380                 printf("          %f %f %f\n",exchange->zeta_ij,.0,.0);
 381         }
 382 }
 383 #endif
 384
 385         /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
 386         exchange->pre_dzeta=0.5*f_a*f_c*db;
 387
 388         /* energy contribution */
 389         energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
 390         moldyn->energy+=energy;
 391         ai->e+=energy;
 392
 393         /* reset k counter for second k loop */
 394         exchange->kcount=0;
 395
 396         return 0;
 397 }
 398
 399 /* albe 3 body potential function (second k loop) */
 400 int albe_mult_3bp_k2(t_moldyn *moldyn,
 401                      t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 402
 403         t_albe_mult_params *params;
 404         t_albe_exchange *exchange;
 405         int kcount;
 406         t_3dvec dist_ik,dist_ij;
 407         double d_ik2,d_ik,d_ij2,d_ij;
 408         unsigned char brand;
 409         double S2;
 410         double g,dg,cos_theta;
 411         double pre_dzeta;
 412         double f_c_ik,df_c_ik;
 413         double dijdik_inv,fcdg,dfcg;
 414         t_3dvec dcosdrj,dcosdrk;
 415         t_3dvec force,tmp;
 416
 417         params=moldyn->pot_params;
 418         exchange=&(params->exchange);
 419         kcount=exchange->kcount;
 420
 421         if(kcount>ALBE_MAXN)
 422                 printf("FATAL: neighbours!\n");
 423
 424         /* d_ik2 */
 425         d_ik2=exchange->d_ik2[kcount];
 426
 427         brand=ak->brand;
 428         if(brand==ai->brand)
 429                 S2=params->S2[brand];
 430         else
 431                 S2=params->S2mixed;
 432
 433         /* return if d_ik > S */
 434         if(d_ik2>S2) {
 435                 exchange->kcount++;
 436                 return 0;
 437         }
 438
 439         /* prefactor dzeta */
 440         pre_dzeta=exchange->pre_dzeta;
 441
 442         /* dist_ik, d_ik */
 443         dist_ik=exchange->dist_ik[kcount];
 444         d_ik=exchange->d_ik[kcount];
 445
 446         /* f_c_ik, df_c_ik */
 447         f_c_ik=exchange->f_c_ik[kcount];
 448         df_c_ik=exchange->df_c_ik[kcount];
 449
 450         /* dist_ij, d_ij2, d_ij */
 451         dist_ij=exchange->dist_ij;
 452         d_ij2=exchange->d_ij2;
 453         d_ij=exchange->d_ij;
 454
 455         /* g, dg, cos_theta */
 456         g=exchange->g[kcount];
 457         dg=exchange->dg[kcount];
 458         cos_theta=exchange->cos_theta[kcount];
 459
 460         /* cos_theta derivatives wrt j,k */
 461         dijdik_inv=1.0/(d_ij*d_ik);
 462         v3_scale(&dcosdrj,&dist_ik,dijdik_inv);         // j
 463         v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
 464         v3_add(&dcosdrj,&dcosdrj,&tmp);
 465         v3_scale(&dcosdrk,&dist_ij,dijdik_inv);         // k
 466         v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
 467         v3_add(&dcosdrk,&dcosdrk,&tmp);
 468
 469         /* f_c_ik * dg, df_c_ik * g */
 470         fcdg=f_c_ik*dg;
 471         dfcg=df_c_ik*g;
 472
 473         /* derivative wrt j */
 474         v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
 475
 476         /* force contribution */
 477         v3_add(&(aj->f),&(aj->f),&force);
 478
 479 #ifdef DEBUG
 480 if(moldyn->time>DSTART&&moldyn->time<DEND) {
 481         if(aj==&(moldyn->atom[DATOM])) {
 482                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
 483                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
 484                 printf("  total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
 485                 printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
 486                 printf("    d ij ik = %f %f\n",d_ij,d_ik);
 487         }
 488 }
 489 #endif
 490
 491         /* virial */
 492         virial_calc(ai,&force,&dist_ij);
 493
 494         /* force contribution to atom i */
 495         v3_scale(&force,&force,-1.0);
 496         v3_add(&(ai->f),&(ai->f),&force);
 497
 498         /* derivative wrt k */
 499         v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
 500         v3_scale(&tmp,&dcosdrk,fcdg);
 501         v3_add(&force,&force,&tmp);
 502         v3_scale(&force,&force,pre_dzeta);
 503
 504         /* force contribution */
 505         v3_add(&(ak->f),&(ak->f),&force);
 506
 507 #ifdef DEBUG
 508 if(moldyn->time>DSTART&&moldyn->time<DEND) {
 509         if(ak==&(moldyn->atom[DATOM])) {
 510                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
 511                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
 512                 printf("  total k: %f %f %f\n",ak->f.x,ak->f.y,ak->f.z);
 513                 printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
 514                 printf("    d ij ik = %f %f\n",d_ij,d_ik);
 515         }
 516 }
 517 #endif
 518
 519         /* virial */
 520         virial_calc(ai,&force,&dist_ik);
 521
 522         /* force contribution to atom i */
 523         v3_scale(&force,&force,-1.0);
 524         v3_add(&(ai->f),&(ai->f),&force);
 525
 526         /* increase k counter */
 527         exchange->kcount++;
 528
 529         return 0;
 530
 531 }
 532
 533 int albe_mult_check_2b_bond(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,u8 bc) {
 534
 535         t_albe_mult_params *params;
 536         t_3dvec dist;
 537         double d;
 538         u8 brand;
 539
 540         v3_sub(&dist,&(jtom->r),&(itom->r));
 541         if(bc) check_per_bound(moldyn,&dist);
 542         d=v3_absolute_square(&dist);
 543
 544         params=moldyn->pot_params;
 545         brand=itom->brand;
 546
 547         if(brand==jtom->brand) {
 548                 if(d<=params->S2[brand])
 549                         return TRUE;
 550         }
 551         else {
 552                 if(d<=params->S2mixed)
 553                         return TRUE;
 554         }
 555
 556         return FALSE;
 557 }