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         /* alloc mem for potential parameters */
  28         moldyn->pot_params=malloc(sizeof(t_albe_mult_params));
  29         if(moldyn->pot_params==NULL) {
  30                 perror("[albe] pot params alloc");
  31                 return -1;
  32         }
  33
  34         /* these are now albe parameters */
  35         p=moldyn->pot_params;
  36
  37         // only 1 combination by now :p
  38         switch(element1) {
  39                 case SI:
  40                         /* type: silicon */
  41                         p->S[0]=ALBE_S_SI;
  42                         p->R[0]=ALBE_R_SI;
  43                         p->A[0]=ALBE_A_SI;
  44                         p->B[0]=ALBE_B_SI;
  45                         p->r0[0]=ALBE_R0_SI;
  46                         p->lambda[0]=ALBE_LAMBDA_SI;
  47                         p->mu[0]=ALBE_MU_SI;
  48                         p->gamma[0]=ALBE_GAMMA_SI;
  49                         p->c[0]=ALBE_C_SI;
  50                         p->d[0]=ALBE_D_SI;
  51                         p->h[0]=ALBE_H_SI;
  52                         switch(element2) {
  53                                 case C:
  54                                         /* type: carbon */
  55                                         p->S[1]=ALBE_S_C;
  56                                         p->R[1]=ALBE_R_C;
  57                                         p->A[1]=ALBE_A_C;
  58                                         p->B[1]=ALBE_B_C;
  59                                         p->r0[1]=ALBE_R0_C;
  60                                         p->lambda[1]=ALBE_LAMBDA_C;
  61                                         p->mu[1]=ALBE_MU_C;
  62                                         p->gamma[1]=ALBE_GAMMA_C;
  63                                         p->c[1]=ALBE_C_C;
  64                                         p->d[1]=ALBE_D_C;
  65                                         p->h[1]=ALBE_H_C;
  66                                         /* mixed type: silicon carbide */
  67                                         p->Smixed=ALBE_S_SIC;
  68                                         p->Rmixed=ALBE_R_SIC;
  69                                         p->Amixed=ALBE_A_SIC;
  70                                         p->Bmixed=ALBE_B_SIC;
  71                                         p->r0_mixed=ALBE_R0_SIC;
  72                                         p->lambda_m=ALBE_LAMBDA_SIC;
  73                                         p->mu_m=ALBE_MU_SIC;
  74                                         p->gamma_m=ALBE_GAMMA_SIC;
  75                                         p->c_mixed=ALBE_C_SIC;
  76                                         p->d_mixed=ALBE_D_SIC;
  77                                         p->h_mixed=ALBE_H_SIC;
  78                                         break;
  79                                 default:
  80                                         printf("[albe] WARNING: element2\n");
  81                                         return -1;
  82                         }
  83                         break;
  84                 default:
  85                         printf("[albe] WARNING: element1\n");
  86                         return -1;
  87         }
  88
  89         printf("[albe] parameter completion\n");
  90         p->S2[0]=p->S[0]*p->S[0];
  91         p->S2[1]=p->S[1]*p->S[1];
  92         p->S2mixed=p->Smixed*p->Smixed;
  93
  94         printf("[albe] mult parameter info:\n");
  95         printf("  S (A)  | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed);
  96         printf("  R (A)  | %f | %f | %f\n",p->R[0],p->R[1],p->Rmixed);
  97         printf("  A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV);
  98         printf("  B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV);
  99         printf("  lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1],
 100                                           p->lambda_m);
 101         printf("  mu     | %f | %f | %f\n",p->mu[0],p->mu[1],p->mu_m);
 102         printf("  gamma  | %f | %f\n",p->gamma[0],p->gamma[1]);
 103         printf("  c      | %f | %f\n",p->c[0],p->c[1]);
 104         printf("  d      | %f | %f\n",p->d[0],p->d[1]);
 105         printf("  h      | %f | %f\n",p->h[0],p->h[1]);
 106
 107         return 0;
 108 }
 109
 110 /* albe 3 body potential function (first ij loop) */
 111 int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 112
 113         t_albe_mult_params *params;
 114         t_albe_exchange *exchange;
 115         unsigned char brand;
 116         double S2;
 117         t_3dvec dist_ij;
 118         double d_ij2,d_ij;
 119
 120         params=moldyn->pot_params;
 121         exchange=&(params->exchange);
 122
 123         /* reset zeta sum */
 124         exchange->zeta_ij=0.0;
 125
 126         /*
 127          * set ij depending values
 128          */
 129
 130         brand=ai->brand;
 131         if(brand==aj->brand) {
 132                 S2=params->S2[brand];
 133         }
 134         else {
 135                 S2=params->S2mixed;
 136         }
 137
 138         /* dist_ij, d_ij2 */
 139         v3_sub(&dist_ij,&(aj->r),&(ai->r));
 140         if(bc) check_per_bound(moldyn,&dist_ij);
 141         d_ij2=v3_absolute_square(&dist_ij);
 142
 143         /* if d_ij2 > S2 => no force & potential energy contribution */
 144         if(d_ij2>S2) {
 145                 moldyn->run3bp=0;
 146                 return 0;
 147         }
 148
 149         /* d_ij */
 150         d_ij=sqrt(d_ij2);
 151
 152         /* store values */
 153         exchange->dist_ij=dist_ij;
 154         exchange->d_ij2=d_ij2;
 155         exchange->d_ij=d_ij;
 156
 157         /* reset k counter for first k loop */
 158         exchange->kcount=0;
 159
 160         return 0;
 161 }
 162
 163 /* albe 3 body potential function (first k loop) */
 164 int albe_mult_3bp_k1(t_moldyn *moldyn,
 165                      t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 166
 167         t_albe_mult_params *params;
 168         t_albe_exchange *exchange;
 169         unsigned char brand;
 170         double R,S,S2;
 171         t_3dvec dist_ij,dist_ik;
 172         double d_ik2,d_ik,d_ij;
 173         double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
 174         double f_c_ik,df_c_ik;
 175         int kcount;
 176
 177         params=moldyn->pot_params;
 178         exchange=&(params->exchange);
 179         kcount=exchange->kcount;
 180
 181         if(kcount>ALBE_MAXN) {
 182                 printf("FATAL: neighbours = %d\n",kcount);
 183                 printf("  -> %d %d %d\n",ai->tag,aj->tag,ak->tag);
 184         }
 185
 186         /* ik constants */
 187         brand=ai->brand;
 188         if(brand==ak->brand) {
 189                 R=params->R[brand];
 190                 S=params->S[brand];
 191                 S2=params->S2[brand];
 192                 /* albe needs i,k depending c,d,h and gamma values */
 193                 exchange->gamma_i=&(params->gamma[brand]);
 194                 exchange->c_i=&(params->c[brand]);
 195                 exchange->d_i=&(params->d[brand]);
 196                 exchange->h_i=&(params->h[brand]);
 197         }
 198         else {
 199                 R=params->Rmixed;
 200                 S=params->Smixed;
 201                 S2=params->S2mixed;
 202                 /* albe needs i,k depending c,d,h and gamma values */
 203                 exchange->gamma_i=&(params->gamma_m);
 204                 exchange->c_i=&(params->c_mixed);
 205                 exchange->d_i=&(params->d_mixed);
 206                 exchange->h_i=&(params->h_mixed);
 207         }
 208         exchange->ci2=*(exchange->c_i)**(exchange->c_i);
 209         exchange->di2=*(exchange->d_i)**(exchange->d_i);
 210         exchange->ci2di2=exchange->ci2/exchange->di2;
 211
 212         /* dist_ik, d_ik2 */
 213         v3_sub(&dist_ik,&(ak->r),&(ai->r));
 214         if(bc) check_per_bound(moldyn,&dist_ik);
 215         d_ik2=v3_absolute_square(&dist_ik);
 216
 217         /* store data for second k loop */
 218         exchange->dist_ik[kcount]=dist_ik;
 219         exchange->d_ik2[kcount]=d_ik2;
 220
 221         /* return if not within cutoff */
 222         if(d_ik2>S2) {
 223                 exchange->kcount++;
 224                 return 0;
 225         }
 226
 227         /* d_ik */
 228         d_ik=sqrt(d_ik2);
 229
 230         /* dist_ij, d_ij */
 231         dist_ij=exchange->dist_ij;
 232         d_ij=exchange->d_ij;
 233
 234         /* cos theta */
 235         cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
 236
 237         /* g_ijk */
 238         h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism
 239         d2_h_cos2=exchange->di2+(h_cos*h_cos);
 240         frac=exchange->ci2/d2_h_cos2;
 241         g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
 242         dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
 243
 244         /* zeta sum += f_c_ik * g_ijk */
 245         if(d_ik<=R) {
 246                 exchange->zeta_ij+=g;
 247                 f_c_ik=1.0;
 248                 df_c_ik=0.0;
 249         }
 250         else {
 251                 s_r=S-R;
 252                 arg=M_PI*(d_ik-R)/s_r;
 253                 f_c_ik=0.5+0.5*cos(arg);
 254                 df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
 255                 exchange->zeta_ij+=f_c_ik*g;
 256         }
 257
 258         /* store even more data for second k loop */
 259         exchange->g[kcount]=g;
 260         exchange->dg[kcount]=dg;
 261         exchange->d_ik[kcount]=d_ik;
 262         exchange->cos_theta[kcount]=cos_theta;
 263         exchange->f_c_ik[kcount]=f_c_ik;
 264         exchange->df_c_ik[kcount]=df_c_ik;
 265
 266         /* increase k counter */
 267         exchange->kcount++;
 268
 269         return 0;
 270 }
 271
 272 int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 273
 274         t_albe_mult_params *params;
 275         t_albe_exchange *exchange;
 276         t_3dvec force;
 277         double f_a,df_a,b,db,f_c,df_c;
 278         double f_r,df_r;
 279         double scale;
 280         double mu,B;
 281         double lambda,A;
 282         double d_ij,r0;
 283         unsigned char brand;
 284         double S,R,s_r,arg;
 285         double energy;
 286
 287         params=moldyn->pot_params;
 288         exchange=&(params->exchange);
 289
 290         brand=aj->brand;
 291         if(brand==ai->brand) {
 292                 S=params->S[brand];
 293                 R=params->R[brand];
 294                 B=params->B[brand];
 295                 A=params->A[brand];
 296                 r0=params->r0[brand];
 297                 mu=params->mu[brand];
 298                 lambda=params->lambda[brand];
 299         }
 300         else {
 301                 S=params->Smixed;
 302                 R=params->Rmixed;
 303                 B=params->Bmixed;
 304                 A=params->Amixed;
 305                 r0=params->r0_mixed;
 306                 mu=params->mu_m;
 307                 lambda=params->lambda_m;
 308         }
 309
 310         d_ij=exchange->d_ij;
 311
 312         /* f_c, df_c */
 313         if(d_ij<R) {
 314                 f_c=1.0;
 315                 df_c=0.0;
 316         }
 317         else {
 318                 s_r=S-R;
 319                 arg=M_PI*(d_ij-R)/s_r;
 320                 f_c=0.5+0.5*cos(arg);
 321                 df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
 322         }
 323
 324         /* f_a, df_a */
 325         f_a=-B*exp(-mu*(d_ij-r0));
 326         df_a=mu*f_a/d_ij;
 327
 328         /* f_r, df_r */
 329         f_r=A*exp(-lambda*(d_ij-r0));
 330         df_r=lambda*f_r/d_ij;
 331
 332         /* b, db */
 333         if(exchange->zeta_ij==0.0) {
 334                 b=1.0;
 335                 db=0.0;
 336         }
 337         else {
 338                 b=1.0/sqrt(1.0+exchange->zeta_ij);
 339                 db=-0.5*b/(1.0+exchange->zeta_ij);
 340         }
 341
 342         /* force contribution for atom i */
 343         scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
 344         v3_scale(&force,&(exchange->dist_ij),scale);
 345         v3_add(&(ai->f),&(ai->f),&force);
 346
 347         /* force contribution for atom j */
 348         v3_scale(&force,&force,-1.0); // dri rij = - drj rij
 349         v3_add(&(aj->f),&(aj->f),&force);
 350
 351         /* virial */
 352         virial_calc(aj,&force,&(exchange->dist_ij));
 353
 354 #ifdef DEBUG
 355 if(moldyn->time>DSTART&&moldyn->time<DEND) {
 356         if((ai==&(moldyn->atom[DATOM]))|(aj==&(moldyn->atom[DATOM]))) {
 357                 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,aj->tag);
 358                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
 359                 if(ai==&(moldyn->atom[0]))
 360                         printf("  total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
 361                 if(aj==&(moldyn->atom[0]))
 362                         printf("  total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
 363                 printf("  energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
 364                                                     f_c,b,f_a,f_r);
 365                 printf("          %f %f %f\n",exchange->zeta_ij,.0,.0);
 366         }
 367 }
 368 #endif
 369
 370         /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
 371         exchange->pre_dzeta=0.5*f_a*f_c*db;
 372
 373         /* energy contribution */
 374         energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
 375         moldyn->energy+=energy;
 376         ai->e+=energy;
 377
 378         /* reset k counter for second k loop */
 379         exchange->kcount=0;
 380
 381         return 0;
 382 }
 383
 384 /* albe 3 body potential function (second k loop) */
 385 int albe_mult_3bp_k2(t_moldyn *moldyn,
 386                      t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 387
 388         t_albe_mult_params *params;
 389         t_albe_exchange *exchange;
 390         int kcount;
 391         t_3dvec dist_ik,dist_ij;
 392         double d_ik2,d_ik,d_ij2,d_ij;
 393         unsigned char brand;
 394         double S2;
 395         double g,dg,cos_theta;
 396         double pre_dzeta;
 397         double f_c_ik,df_c_ik;
 398         double dijdik_inv,fcdg,dfcg;
 399         t_3dvec dcosdrj,dcosdrk;
 400         t_3dvec force,tmp;
 401
 402         params=moldyn->pot_params;
 403         exchange=&(params->exchange);
 404         kcount=exchange->kcount;
 405
 406         if(kcount>ALBE_MAXN)
 407                 printf("FATAL: neighbours!\n");
 408
 409         /* d_ik2 */
 410         d_ik2=exchange->d_ik2[kcount];
 411
 412         brand=ak->brand;
 413         if(brand==ai->brand)
 414                 S2=params->S2[brand];
 415         else
 416                 S2=params->S2mixed;
 417
 418         /* return if d_ik > S */
 419         if(d_ik2>S2) {
 420                 exchange->kcount++;
 421                 return 0;
 422         }
 423
 424         /* prefactor dzeta */
 425         pre_dzeta=exchange->pre_dzeta;
 426
 427         /* dist_ik, d_ik */
 428         dist_ik=exchange->dist_ik[kcount];
 429         d_ik=exchange->d_ik[kcount];
 430
 431         /* f_c_ik, df_c_ik */
 432         f_c_ik=exchange->f_c_ik[kcount];
 433         df_c_ik=exchange->df_c_ik[kcount];
 434
 435         /* dist_ij, d_ij2, d_ij */
 436         dist_ij=exchange->dist_ij;
 437         d_ij2=exchange->d_ij2;
 438         d_ij=exchange->d_ij;
 439
 440         /* g, dg, cos_theta */
 441         g=exchange->g[kcount];
 442         dg=exchange->dg[kcount];
 443         cos_theta=exchange->cos_theta[kcount];
 444
 445         /* cos_theta derivatives wrt j,k */
 446         dijdik_inv=1.0/(d_ij*d_ik);
 447         v3_scale(&dcosdrj,&dist_ik,dijdik_inv);         // j
 448         v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
 449         v3_add(&dcosdrj,&dcosdrj,&tmp);
 450         v3_scale(&dcosdrk,&dist_ij,dijdik_inv);         // k
 451         v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
 452         v3_add(&dcosdrk,&dcosdrk,&tmp);
 453
 454         /* f_c_ik * dg, df_c_ik * g */
 455         fcdg=f_c_ik*dg;
 456         dfcg=df_c_ik*g;
 457
 458         /* derivative wrt j */
 459         v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
 460
 461         /* force contribution */
 462         v3_add(&(aj->f),&(aj->f),&force);
 463
 464 #ifdef DEBUG
 465 if(moldyn->time>DSTART&&moldyn->time<DEND) {
 466         if(aj==&(moldyn->atom[DATOM])) {
 467                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
 468                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
 469                 printf("  total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
 470                 printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
 471                 printf("    d ij ik = %f %f\n",d_ij,d_ik);
 472         }
 473 }
 474 #endif
 475
 476         /* force contribution to atom i */
 477         v3_scale(&force,&force,-1.0);
 478         v3_add(&(ai->f),&(ai->f),&force);
 479
 480         /* virial */
 481         virial_calc(ai,&force,&dist_ij);
 482
 483         /* derivative wrt k */
 484         v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
 485         v3_scale(&tmp,&dcosdrk,fcdg);
 486         v3_add(&force,&force,&tmp);
 487         v3_scale(&force,&force,pre_dzeta);
 488
 489         /* force contribution */
 490         v3_add(&(ak->f),&(ak->f),&force);
 491
 492 #ifdef DEBUG
 493 if(moldyn->time>DSTART&&moldyn->time<DEND) {
 494         if(ak==&(moldyn->atom[DATOM])) {
 495                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
 496                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
 497                 printf("  total k: %f %f %f\n",ak->f.x,ak->f.y,ak->f.z);
 498                 printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
 499                 printf("    d ij ik = %f %f\n",d_ij,d_ik);
 500         }
 501 }
 502 #endif
 503
 504         /* force contribution to atom i */
 505         v3_scale(&force,&force,-1.0);
 506         v3_add(&(ai->f),&(ai->f),&force);
 507
 508         /* virial */
 509         virial_calc(ai,&force,&dist_ik);
 510
 511         /* increase k counter */
 512         exchange->kcount++;
 513
 514         return 0;
 515
 516 }
 517
 518 int albe_mult_check_2b_bond(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,u8 bc) {
 519
 520         t_albe_mult_params *params;
 521         t_3dvec dist;
 522         double d;
 523         u8 brand;
 524
 525         v3_sub(&dist,&(jtom->r),&(itom->r));
 526         if(bc) check_per_bound(moldyn,&dist);
 527         d=v3_absolute_square(&dist);
 528
 529         params=moldyn->pot_params;
 530         brand=itom->brand;
 531
 532         if(brand==jtom->brand) {
 533                 if(d<=params->S2[brand])
 534                         return TRUE;
 535         }
 536         else {
 537                 if(d<=params->S2mixed)
 538                         return TRUE;
 539         }
 540
 541         return FALSE;
 542 }