potentials/tersoff.c

   1 /*
   2  * tersoff.c - tersoff 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 "tersoff.h"
  21
  22 /* create mixed terms from parameters and set them */
  23 int tersoff_mult_complete_params(t_tersoff_mult_params *p) {
  24
  25         printf("[moldyn] tersoff 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->Smixed=sqrt(p->S[0]*p->S[1]);
  29         p->S2mixed=p->Smixed*p->Smixed;
  30         p->Rmixed=sqrt(p->R[0]*p->R[1]);
  31         p->Amixed=sqrt(p->A[0]*p->A[1]);
  32         p->Bmixed=sqrt(p->B[0]*p->B[1]);
  33         p->lambda_m=0.5*(p->lambda[0]+p->lambda[1]);
  34         p->mu_m=0.5*(p->mu[0]+p->mu[1]);
  35
  36         printf("[moldyn] tersoff mult parameter info:\n");
  37         printf("  S (A)  | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed);
  38         printf("  R (A)  | %f | %f | %f\n",p->R[0],p->R[1],p->Rmixed);
  39         printf("  A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV);
  40         printf("  B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV);
  41         printf("  lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1],
  42                                           p->lambda_m);
  43         printf("  mu     | %f | %f | %f\n",p->mu[0],p->mu[1],p->mu_m);
  44         printf("  beta   | %.10f | %.10f\n",p->beta[0],p->beta[1]);
  45         printf("  n      | %f | %f\n",p->n[0],p->n[1]);
  46         printf("  c      | %f | %f\n",p->c[0],p->c[1]);
  47         printf("  d      | %f | %f\n",p->d[0],p->d[1]);
  48         printf("  h      | %f | %f\n",p->h[0],p->h[1]);
  49         printf("  chi    | %f \n",p->chi);
  50
  51         return 0;
  52 }
  53
  54 /* tersoff 1 body part */
  55 int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) {
  56
  57         int brand;
  58         t_tersoff_mult_params *params;
  59         t_tersoff_exchange *exchange;
  60
  61         brand=ai->brand;
  62         params=moldyn->pot_params;
  63         exchange=&(params->exchange);
  64
  65         /*
  66          * simple: point constant parameters only depending on atom i to
  67          *         their right values
  68          */
  69
  70         exchange->beta_i=&(params->beta[brand]);
  71         exchange->n_i=&(params->n[brand]);
  72         exchange->c_i=&(params->c[brand]);
  73         exchange->d_i=&(params->d[brand]);
  74         exchange->h_i=&(params->h[brand]);
  75
  76         exchange->betaini=pow(*(exchange->beta_i),*(exchange->n_i));
  77         exchange->ci2=params->c[brand]*params->c[brand];
  78         exchange->di2=params->d[brand]*params->d[brand];
  79         exchange->ci2di2=exchange->ci2/exchange->di2;
  80
  81         return 0;
  82 }
  83
  84 /* tersoff 2 body part */
  85 int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
  86
  87         t_tersoff_mult_params *params;
  88         t_3dvec dist_ij,force;
  89         double d_ij,d_ij2;
  90         double A,R,S,S2,lambda;
  91         double f_r,df_r;
  92         double f_c,df_c;
  93         int brand;
  94         double s_r;
  95         double arg;
  96         double energy;
  97
  98         printf("WARNING! - tersoff_mult_2bp is obsolete.\n");
  99         printf("WARNING! - repulsive part handled in 3bp/j2 routine.\n");
 100
 101         /* use newtons third law */
 102         if(ai<aj) return 0;
 103
 104         params=moldyn->pot_params;
 105         brand=aj->brand;
 106
 107         /* determine cutoff square */
 108         if(brand==ai->brand)
 109                 S2=params->S2[brand];
 110         else
 111                 S2=params->S2mixed;
 112
 113         /* dist_ij, d_ij2 */
 114         v3_sub(&dist_ij,&(aj->r),&(ai->r));
 115         if(bc) check_per_bound(moldyn,&dist_ij);
 116         d_ij2=v3_absolute_square(&dist_ij);
 117
 118         /* if d_ij2 > S2 => no force & potential energy contribution */
 119         if(d_ij2>S2) {
 120                 return 0;
 121         }
 122
 123         /* now we will need the distance */
 124         d_ij=sqrt(d_ij2);
 125
 126         /* more constants */
 127         if(brand==ai->brand) {
 128                 S=params->S[brand];
 129                 R=params->R[brand];
 130                 A=params->A[brand];
 131                 lambda=params->lambda[brand];
 132         }
 133         else {
 134                 S=params->Smixed;
 135                 R=params->Rmixed;
 136                 A=params->Amixed;
 137                 lambda=params->lambda_m;
 138         }
 139
 140         /* f_r_ij, df_r_ij */
 141         f_r=A*exp(-lambda*d_ij);
 142         df_r=lambda*f_r/d_ij;
 143
 144         /* f_c, df_c */
 145         if(d_ij<R) {
 146                 f_c=1.0;
 147                 df_c=0.0;
 148                 v3_scale(&force,&dist_ij,-df_r);
 149         }
 150         else {
 151                 s_r=S-R;
 152                 arg=M_PI*(d_ij-R)/s_r;
 153                 f_c=0.5+0.5*cos(arg);
 154                 df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
 155                 v3_scale(&force,&dist_ij,-df_c*f_r-df_r*f_c);
 156         }
 157
 158         /* add forces */
 159         v3_add(&(ai->f),&(ai->f),&force);
 160         v3_sub(&(aj->f),&(aj->f),&force); // reason: dri rij = - drj rij
 161
 162 #ifdef DEBUG
 163         if((ai==&(moldyn->atom[0]))|(aj==&(moldyn->atom[0]))) {
 164                 printf("force 2bp: [%d %d]\n",ai->tag,aj->tag);
 165                 printf("adding %f %f %f\n",force.x,force.y,force.z);
 166                 if(ai==&(moldyn->atom[0]))
 167                         printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
 168                 if(aj==&(moldyn->atom[0]))
 169                         printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
 170         }
 171 #endif
 172
 173         /* virial */
 174         virial_calc(ai,&force,&dist_ij);
 175
 176         /* energy 2bp contribution */
 177         energy=f_r*f_c;
 178         moldyn->energy+=energy;
 179         ai->e+=0.5*energy;
 180         aj->e+=0.5*energy;
 181
 182         return 0;
 183 }
 184
 185 /* tersoff 3 body potential function (first ij loop) */
 186 int tersoff_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 187
 188         t_tersoff_mult_params *params;
 189         t_tersoff_exchange *exchange;
 190         unsigned char brand;
 191         double S2;
 192         t_3dvec dist_ij;
 193         double d_ij2,d_ij;
 194
 195         params=moldyn->pot_params;
 196         exchange=&(params->exchange);
 197
 198         /* reset zeta sum */
 199         exchange->zeta_ij=0.0;
 200
 201         /*
 202          * set ij depending values
 203          */
 204
 205         brand=ai->brand;
 206
 207         if(brand==aj->brand)
 208                 S2=params->S2[brand];
 209         else
 210                 S2=params->S2mixed;
 211
 212         /* dist_ij, d_ij2 */
 213         v3_sub(&dist_ij,&(aj->r),&(ai->r));
 214         if(bc) check_per_bound(moldyn,&dist_ij);
 215         d_ij2=v3_absolute_square(&dist_ij);
 216
 217         /* if d_ij2 > S2 => no force & potential energy contribution */
 218         if(d_ij2>S2) {
 219                 moldyn->run3bp=0;
 220                 return 0;
 221         }
 222
 223         /* d_ij */
 224         d_ij=sqrt(d_ij2);
 225
 226         /* store values */
 227         exchange->dist_ij=dist_ij;
 228         exchange->d_ij2=d_ij2;
 229         exchange->d_ij=d_ij;
 230
 231         /* reset k counter for first k loop */
 232         exchange->kcount=0;
 233
 234         return 0;
 235 }
 236
 237 /* tersoff 3 body potential function (first k loop) */
 238 int tersoff_mult_3bp_k1(t_moldyn *moldyn,
 239                         t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 240
 241         t_tersoff_mult_params *params;
 242         t_tersoff_exchange *exchange;
 243         unsigned char brand;
 244         double R,S,S2;
 245         t_3dvec dist_ij,dist_ik;
 246         double d_ik2,d_ik,d_ij;
 247         double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
 248         double f_c_ik,df_c_ik;
 249         int kcount;
 250
 251         params=moldyn->pot_params;
 252         exchange=&(params->exchange);
 253         kcount=exchange->kcount;
 254
 255         if(kcount>TERSOFF_MAXN) {
 256                 printf("FATAL: neighbours = %d\n",kcount);
 257                 printf("  -> %d %d %d\n",ai->tag,aj->tag,ak->tag);
 258         }
 259
 260         /* ik constants */
 261         brand=ai->brand;
 262         if(brand==ak->brand) {
 263                 R=params->R[brand];
 264                 S=params->S[brand];
 265                 S2=params->S2[brand];
 266         }
 267         else {
 268                 R=params->Rmixed;
 269                 S=params->Smixed;
 270                 S2=params->S2mixed;
 271         }
 272
 273         /* dist_ik, d_ik2 */
 274         v3_sub(&dist_ik,&(ak->r),&(ai->r));
 275         if(bc) check_per_bound(moldyn,&dist_ik);
 276         d_ik2=v3_absolute_square(&dist_ik);
 277
 278         /* store data for second k loop */
 279         exchange->dist_ik[kcount]=dist_ik;
 280         exchange->d_ik2[kcount]=d_ik2;
 281
 282         /* return if not within cutoff */
 283         if(d_ik2>S2) {
 284                 exchange->kcount++;
 285                 return 0;
 286         }
 287
 288         /* d_ik */
 289         d_ik=sqrt(d_ik2);
 290
 291         /* dist_ij, d_ij */
 292         dist_ij=exchange->dist_ij;
 293         d_ij=exchange->d_ij;
 294
 295         /* cos theta */
 296         cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
 297
 298         /* g_ijk */
 299         h_cos=*(exchange->h_i)-cos_theta;
 300         d2_h_cos2=exchange->di2+(h_cos*h_cos);
 301         frac=exchange->ci2/d2_h_cos2;
 302         g=1.0+exchange->ci2di2-frac;
 303         dg=-2.0*frac*h_cos/d2_h_cos2;
 304
 305         /* zeta sum += f_c_ik * g_ijk */
 306         if(d_ik<=R) {
 307                 exchange->zeta_ij+=g;
 308                 f_c_ik=1.0;
 309                 df_c_ik=0.0;
 310         }
 311         else {
 312                 s_r=S-R;
 313                 arg=M_PI*(d_ik-R)/s_r;
 314                 f_c_ik=0.5+0.5*cos(arg);
 315                 df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
 316                 exchange->zeta_ij+=f_c_ik*g;
 317         }
 318
 319 #ifdef DEBUG
 320         if((ai==&(moldyn->atom[0]))|
 321            (aj==&(moldyn->atom[864]))|
 322            (ak==&(moldyn->atom[1003]))) {
 323                 printf(" -> %f %f %f\n",exchange->ci2di2,frac,h_cos);
 324         }
 325 #endif
 326
 327         /* store even more data for second k loop */
 328         exchange->g[kcount]=g;
 329         exchange->dg[kcount]=dg;
 330         exchange->d_ik[kcount]=d_ik;
 331         exchange->cos_theta[kcount]=cos_theta;
 332         exchange->f_c_ik[kcount]=f_c_ik;
 333         exchange->df_c_ik[kcount]=df_c_ik;
 334
 335         /* increase k counter */
 336         exchange->kcount++;
 337
 338         return 0;
 339 }
 340
 341 int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 342
 343         t_tersoff_mult_params *params;
 344         t_tersoff_exchange *exchange;
 345         t_3dvec force;
 346         double f_a,df_a,b,db,f_c,df_c;
 347         double f_r,df_r;
 348         double scale;
 349         double mu,B,chi;
 350         double lambda,A;
 351         double d_ij;
 352         unsigned char brand;
 353         double ni,tmp;
 354         double S,R,s_r,arg;
 355         double energy;
 356
 357         params=moldyn->pot_params;
 358         exchange=&(params->exchange);
 359
 360         brand=aj->brand;
 361         if(brand==ai->brand) {
 362                 S=params->S[brand];
 363                 R=params->R[brand];
 364                 B=params->B[brand];
 365                 A=params->A[brand];
 366                 mu=params->mu[brand];
 367                 lambda=params->lambda[brand];
 368                 chi=1.0;
 369         }
 370         else {
 371                 S=params->Smixed;
 372                 R=params->Rmixed;
 373                 B=params->Bmixed;
 374                 A=params->Amixed;
 375                 mu=params->mu_m;
 376                 lambda=params->lambda_m;
 377                 chi=params->chi;
 378         }
 379
 380         d_ij=exchange->d_ij;
 381
 382         /* f_c, df_c */
 383         if(d_ij<R) {
 384                 f_c=1.0;
 385                 df_c=0.0;
 386         }
 387         else {
 388                 s_r=S-R;
 389                 arg=M_PI*(d_ij-R)/s_r;
 390                 f_c=0.5+0.5*cos(arg);
 391                 df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
 392         }
 393
 394         /* f_a, df_a */
 395         f_a=-B*exp(-mu*d_ij);
 396         df_a=mu*f_a/d_ij;
 397
 398         /* f_r, df_r */
 399         f_r=A*exp(-lambda*d_ij);
 400         df_r=lambda*f_r/d_ij;
 401
 402         /* b, db */
 403         if(exchange->zeta_ij==0.0) {
 404                 b=chi;
 405                 db=0.0;
 406         }
 407         else {
 408                 ni=*(exchange->n_i);
 409                 tmp=exchange->betaini*pow(exchange->zeta_ij,ni-1.0);
 410                 b=(1.0+exchange->zeta_ij*tmp);
 411                 db=chi*pow(b,-1.0/(2.0*ni)-1.0);
 412                 b=db*b;
 413                 db*=-0.5*tmp;
 414         }
 415
 416         /* force contribution */
 417         scale=-0.5*(f_c*(df_r+b*df_a)+df_c*(f_r+b*df_a));
 418         v3_scale(&force,&(exchange->dist_ij),scale);
 419         v3_add(&(ai->f),&(ai->f),&force);
 420         v3_sub(&(aj->f),&(aj->f),&force); // dri rij = - drj rij
 421
 422 #ifdef DEBUG
 423         if((ai==&(moldyn->atom[0]))|(aj==&(moldyn->atom[0]))) {
 424                 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,aj->tag);
 425                 printf("adding %f %f %f\n",force.x,force.y,force.z);
 426                 if(ai==&(moldyn->atom[0]))
 427                         printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
 428                 if(aj==&(moldyn->atom[0]))
 429                         printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
 430                 printf("energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
 431                                                     f_c,b,f_a,f_r);
 432                 printf("        %f %f %f\n",exchange->zeta_ij,.0,.0);
 433         }
 434 #endif
 435
 436         /* virial */
 437         if(aj<ai)
 438                 virial_calc(ai,&force,&(exchange->dist_ij));
 439
 440         /* dzeta prefactor = - 0.5 f_c f_a db */
 441         exchange->pre_dzeta=-0.5*f_a*f_c*db;
 442
 443         /* energy contribution */
 444         energy=0.5*f_c*(b*f_a+f_r);
 445         moldyn->energy+=energy;
 446         ai->e+=energy;
 447
 448         /* reset k counter for second k loop */
 449         exchange->kcount=0;
 450
 451         return 0;
 452 }
 453
 454 /* tersoff 3 body potential function (second k loop) */
 455 int tersoff_mult_3bp_k2(t_moldyn *moldyn,
 456                         t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 457
 458         t_tersoff_mult_params *params;
 459         t_tersoff_exchange *exchange;
 460         int kcount;
 461         t_3dvec dist_ik,dist_ij;
 462         double d_ik2,d_ik,d_ij2,d_ij;
 463         unsigned char brand;
 464         double S2;
 465         double g,dg,cos_theta;
 466         double pre_dzeta;
 467         double f_c_ik,df_c_ik;
 468         double dijdik_inv,fcdg,dfcg;
 469         t_3dvec dcosdri,dcosdrj,dcosdrk;
 470         t_3dvec force,tmp;
 471
 472         params=moldyn->pot_params;
 473         exchange=&(params->exchange);
 474         kcount=exchange->kcount;
 475
 476         if(kcount>TERSOFF_MAXN)
 477                 printf("FATAL: neighbours!\n");
 478
 479         /* d_ik2 */
 480         d_ik2=exchange->d_ik2[kcount];
 481
 482         brand=ak->brand;
 483         if(brand==ai->brand)
 484                 S2=params->S2[brand];
 485         else
 486                 S2=params->S2mixed;
 487
 488         /* return if d_ik > S */
 489         if(d_ik2>S2) {
 490                 exchange->kcount++;
 491                 return 0;
 492         }
 493
 494         /* prefactor dzeta */
 495         pre_dzeta=exchange->pre_dzeta;
 496
 497         /* dist_ik, d_ik */
 498         dist_ik=exchange->dist_ik[kcount];
 499         d_ik=exchange->d_ik[kcount];
 500
 501         /* f_c_ik, df_c_ik */
 502         f_c_ik=exchange->f_c_ik[kcount];
 503         df_c_ik=exchange->df_c_ik[kcount];
 504
 505         /* dist_ij, d_ij2, d_ij */
 506         dist_ij=exchange->dist_ij;
 507         d_ij2=exchange->d_ij2;
 508         d_ij=exchange->d_ij;
 509
 510         /* g, dg, cos_theta */
 511         g=exchange->g[kcount];
 512         dg=exchange->dg[kcount];
 513         cos_theta=exchange->cos_theta[kcount];
 514
 515         /* cos_theta derivatives wrt i,j,k */
 516         dijdik_inv=1.0/(d_ij*d_ik);
 517         v3_scale(&dcosdrj,&dist_ik,dijdik_inv);
 518         v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
 519         v3_add(&dcosdrj,&dcosdrj,&tmp);
 520         v3_scale(&dcosdrk,&dist_ij,dijdik_inv);
 521         v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
 522         v3_add(&dcosdrk,&dcosdrk,&tmp);
 523         v3_add(&dcosdri,&dcosdrj,&dcosdrk);
 524         v3_scale(&dcosdri,&dcosdri,-1.0);
 525
 526         /* f_c_ik * dg, df_c_ik * g */
 527         fcdg=f_c_ik*dg;
 528         dfcg=df_c_ik*g;
 529
 530         /* derivative wrt i */
 531         v3_scale(&force,&dist_ik,dfcg);
 532         v3_scale(&tmp,&dcosdri,fcdg);
 533         v3_add(&force,&force,&tmp);
 534         v3_scale(&force,&force,pre_dzeta);
 535
 536         /* force contribution */
 537         v3_add(&(ai->f),&(ai->f),&force);
 538
 539 #ifdef DEBUG
 540         if(ai==&(moldyn->atom[0])) {
 541                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
 542                 printf("adding %f %f %f\n",force.x,force.y,force.z);
 543                 printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
 544         }
 545 #endif
 546
 547         /* virial */
 548         //virial_calc(ai,&force,&dist_ij);
 549
 550         /* derivative wrt j */
 551         v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
 552
 553         /* force contribution */
 554         v3_add(&(aj->f),&(aj->f),&force);
 555
 556 #ifdef DEBUG
 557         if(aj==&(moldyn->atom[0])) {
 558                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
 559                 printf("adding %f %f %f\n",force.x,force.y,force.z);
 560                 printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
 561         }
 562 #endif
 563
 564         /* virial */
 565         //v3_scale(&force,&force,-1.0);
 566         if(aj<ai)
 567                 virial_calc(ai,&force,&dist_ij);
 568
 569         /* derivative wrt k */
 570         v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
 571         v3_scale(&tmp,&dcosdrk,fcdg);
 572         v3_add(&force,&force,&tmp);
 573         v3_scale(&force,&force,pre_dzeta);
 574
 575         /* force contribution */
 576         v3_add(&(ak->f),&(ak->f),&force);
 577
 578 #ifdef DEBUG
 579         if(ak==&(moldyn->atom[0])) {
 580                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
 581                 printf("adding %f %f %f\n",force.x,force.y,force.z);
 582                 printf("total k: %f %f %f\n",ak->f.x,ak->f.y,ak->f.z);
 583         }
 584 #endif
 585
 586         /* virial */
 587         //v3_scale(&force,&force,-1.0);
 588         if(aj<ai)
 589                 virial_calc(ai,&force,&dist_ik);
 590
 591         /* increase k counter */
 592         exchange->kcount++;
 593
 594         return 0;
 595
 596 }