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_complete_params(t_albe_mult_params *p) {
  24
  25         printf("[moldyn] albe 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->S2mixed=p->Smixed*p->Smixed;
  29
  30         printf("[moldyn] albe mult parameter info:\n");
  31         printf("  S (A)  | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed);
  32         printf("  R (A)  | %f | %f | %f\n",p->R[0],p->R[1],p->Rmixed);
  33         printf("  A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV);
  34         printf("  B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV);
  35         printf("  lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1],
  36                                           p->lambda_m);
  37         printf("  mu     | %f | %f | %f\n",p->mu[0],p->mu[1],p->mu_m);
  38         printf("  gamma  | %f | %f\n",p->gamma[0],p->gamma[1]);
  39         printf("  c      | %f | %f\n",p->c[0],p->c[1]);
  40         printf("  d      | %f | %f\n",p->d[0],p->d[1]);
  41         printf("  h      | %f | %f\n",p->h[0],p->h[1]);
  42
  43         return 0;
  44 }
  45
  46 /* albe 1 body part */
  47 int albe_mult_1bp(t_moldyn *moldyn,t_atom *ai) {
  48
  49         int brand;
  50         t_albe_mult_params *params;
  51         t_albe_exchange *exchange;
  52
  53         brand=ai->brand;
  54         params=moldyn->pot_params;
  55         exchange=&(params->exchange);
  56
  57         /*
  58          * simple: point constant parameters only depending on atom i to
  59          *         their right values
  60          */
  61
  62         exchange->gamma_i=&(params->gamma[brand]);
  63         exchange->c_i=&(params->c[brand]);
  64         exchange->d_i=&(params->d[brand]);
  65         exchange->h_i=&(params->h[brand]);
  66
  67         exchange->ci2=params->c[brand]*params->c[brand];
  68         exchange->di2=params->d[brand]*params->d[brand];
  69         exchange->ci2di2=exchange->ci2/exchange->di2;
  70
  71         return 0;
  72 }
  73
  74 /* albe 3 body potential function (first ij loop) */
  75 int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
  76
  77         t_albe_mult_params *params;
  78         t_albe_exchange *exchange;
  79         unsigned char brand;
  80         double S2;
  81         t_3dvec dist_ij;
  82         double d_ij2,d_ij;
  83
  84         params=moldyn->pot_params;
  85         exchange=&(params->exchange);
  86
  87         /* reset zeta sum */
  88         exchange->zeta_ij=0.0;
  89
  90         /*
  91          * set ij depending values
  92          */
  93
  94         brand=ai->brand;
  95
  96         if(brand==aj->brand)
  97                 S2=params->S2[brand];
  98         else
  99                 S2=params->S2mixed;
 100
 101         /* dist_ij, d_ij2 */
 102         v3_sub(&dist_ij,&(aj->r),&(ai->r));
 103         if(bc) check_per_bound(moldyn,&dist_ij);
 104         d_ij2=v3_absolute_square(&dist_ij);
 105
 106         /* if d_ij2 > S2 => no force & potential energy contribution */
 107         if(d_ij2>S2) {
 108                 moldyn->run3bp=0;
 109                 return 0;
 110         }
 111
 112         /* d_ij */
 113         d_ij=sqrt(d_ij2);
 114
 115         /* store values */
 116         exchange->dist_ij=dist_ij;
 117         exchange->d_ij2=d_ij2;
 118         exchange->d_ij=d_ij;
 119
 120         /* reset k counter for first k loop */
 121         exchange->kcount=0;
 122
 123         return 0;
 124 }
 125
 126 /* albe 3 body potential function (first k loop) */
 127 int albe_mult_3bp_k1(t_moldyn *moldyn,
 128                         t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 129
 130         t_albe_mult_params *params;
 131         t_albe_exchange *exchange;
 132         unsigned char brand;
 133         double R,S,S2;
 134         t_3dvec dist_ij,dist_ik;
 135         double d_ik2,d_ik,d_ij;
 136         double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
 137         double f_c_ik,df_c_ik;
 138         int kcount;
 139
 140         params=moldyn->pot_params;
 141         exchange=&(params->exchange);
 142         kcount=exchange->kcount;
 143
 144         if(kcount>ALBE_MAXN) {
 145                 printf("FATAL: neighbours = %d\n",kcount);
 146                 printf("  -> %d %d %d\n",ai->tag,aj->tag,ak->tag);
 147         }
 148
 149         /* ik constants */
 150         brand=ai->brand;
 151         if(brand==ak->brand) {
 152                 R=params->R[brand];
 153                 S=params->S[brand];
 154                 S2=params->S2[brand];
 155         }
 156         else {
 157                 R=params->Rmixed;
 158                 S=params->Smixed;
 159                 S2=params->S2mixed;
 160         }
 161
 162         /* dist_ik, d_ik2 */
 163         v3_sub(&dist_ik,&(ak->r),&(ai->r));
 164         if(bc) check_per_bound(moldyn,&dist_ik);
 165         d_ik2=v3_absolute_square(&dist_ik);
 166
 167         /* store data for second k loop */
 168         exchange->dist_ik[kcount]=dist_ik;
 169         exchange->d_ik2[kcount]=d_ik2;
 170
 171         /* return if not within cutoff */
 172         if(d_ik2>S2) {
 173                 exchange->kcount++;
 174                 return 0;
 175         }
 176
 177         /* d_ik */
 178         d_ik=sqrt(d_ik2);
 179
 180         /* dist_ij, d_ij */
 181         dist_ij=exchange->dist_ij;
 182         d_ij=exchange->d_ij;
 183
 184         /* cos theta */
 185         cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
 186
 187         /* g_ijk */
 188         h_cos=*(exchange->h_i)-cos_theta;
 189         d2_h_cos2=exchange->di2+(h_cos*h_cos);
 190         frac=exchange->ci2/d2_h_cos2;
 191         g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
 192         dg=-2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2;
 193
 194         /* zeta sum += f_c_ik * g_ijk */
 195         if(d_ik<=R) {
 196                 exchange->zeta_ij+=g;
 197                 f_c_ik=1.0;
 198                 df_c_ik=0.0;
 199         }
 200         else {
 201                 s_r=S-R;
 202                 arg=M_PI*(d_ik-R)/s_r;
 203                 f_c_ik=0.5+0.5*cos(arg);
 204                 df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
 205                 exchange->zeta_ij+=f_c_ik*g;
 206         }
 207
 208         /* store even more data for second k loop */
 209         exchange->g[kcount]=g;
 210         exchange->dg[kcount]=dg;
 211         exchange->d_ik[kcount]=d_ik;
 212         exchange->cos_theta[kcount]=cos_theta;
 213         exchange->f_c_ik[kcount]=f_c_ik;
 214         exchange->df_c_ik[kcount]=df_c_ik;
 215
 216         /* increase k counter */
 217         exchange->kcount++;
 218
 219         return 0;
 220 }
 221
 222 int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 223
 224         t_albe_mult_params *params;
 225         t_albe_exchange *exchange;
 226         t_3dvec force;
 227         double f_a,df_a,b,db,f_c,df_c;
 228         double f_r,df_r;
 229         double scale;
 230         double mu,B;
 231         double lambda,A;
 232         double d_ij,r0;
 233         unsigned char brand;
 234         double S,R,s_r,arg;
 235
 236         params=moldyn->pot_params;
 237         exchange=&(params->exchange);
 238
 239         brand=aj->brand;
 240         if(brand==ai->brand) {
 241                 S=params->S[brand];
 242                 R=params->R[brand];
 243                 B=params->B[brand];
 244                 A=params->A[brand];
 245                 r0=params->r0[brand];
 246                 mu=params->mu[brand];
 247                 lambda=params->lambda[brand];
 248         }
 249         else {
 250                 S=params->Smixed;
 251                 R=params->Rmixed;
 252                 B=params->Bmixed;
 253                 A=params->Amixed;
 254                 r0=params->r0_mixed;
 255                 mu=params->mu_m;
 256                 lambda=params->lambda_m;
 257         }
 258
 259         d_ij=exchange->d_ij;
 260
 261         /* f_c, df_c */
 262         if(d_ij<R) {
 263                 f_c=1.0;
 264                 df_c=0.0;
 265         }
 266         else {
 267                 s_r=S-R;
 268                 arg=M_PI*(d_ij-R)/s_r;
 269                 f_c=0.5+0.5*cos(arg);
 270                 df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
 271         }
 272
 273         /* f_a, df_a */
 274         f_a=-B*exp(-mu*(d_ij-r0));
 275         df_a=mu*f_a/d_ij;
 276
 277         /* f_r, df_r */
 278         f_r=A*exp(-lambda*(d_ij-r0));
 279         df_r=lambda*f_r/d_ij;
 280
 281         /* b, db */
 282         if(exchange->zeta_ij==0.0) {
 283                 b=1.0;
 284                 db=0.0;
 285         }
 286         else {
 287                 b=1.0/sqrt(1.0+exchange->zeta_ij);
 288                 db=-0.5*b/(1.0+exchange->zeta_ij);
 289         }
 290
 291         /* force contribution */
 292         scale=-0.5*(f_c*(df_r+b*df_a)+df_c*(f_r+b*df_a));
 293         v3_scale(&force,&(exchange->dist_ij),scale);
 294         v3_add(&(ai->f),&(ai->f),&force);
 295         v3_sub(&(aj->f),&(aj->f),&force); // dri rij = - drj rij
 296
 297 #ifdef DEBUG
 298         if((ai==&(moldyn->atom[0]))|(aj==&(moldyn->atom[0]))) {
 299                 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,aj->tag);
 300                 printf("adding %f %f %f\n",force.x,force.y,force.z);
 301                 if(ai==&(moldyn->atom[0]))
 302                         printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
 303                 if(aj==&(moldyn->atom[0]))
 304                         printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
 305                 printf("energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
 306                                                     f_c,b,f_a,f_r);
 307         }
 308 #endif
 309
 310         /* virial */
 311         if(aj<ai)
 312                 virial_calc(ai,&force,&(exchange->dist_ij));
 313
 314         /* dzeta prefactor = - 0.5 f_c f_a db */
 315         exchange->pre_dzeta=-0.5*f_a*f_c*db;
 316
 317         /* energy contribution */
 318         moldyn->energy+=0.5*f_c*(f_r+b*f_a);
 319
 320         /* reset k counter for second k loop */
 321         exchange->kcount=0;
 322
 323         return 0;
 324 }
 325
 326 /* albe 3 body potential function (second k loop) */
 327 int albe_mult_3bp_k2(t_moldyn *moldyn,
 328                         t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 329
 330         t_albe_mult_params *params;
 331         t_albe_exchange *exchange;
 332         int kcount;
 333         t_3dvec dist_ik,dist_ij;
 334         double d_ik2,d_ik,d_ij2,d_ij;
 335         unsigned char brand;
 336         double S2;
 337         double g,dg,cos_theta;
 338         double pre_dzeta;
 339         double f_c_ik,df_c_ik;
 340         double dijdik_inv,fcdg,dfcg;
 341         t_3dvec dcosdri,dcosdrj,dcosdrk;
 342         t_3dvec force,tmp;
 343
 344         params=moldyn->pot_params;
 345         exchange=&(params->exchange);
 346         kcount=exchange->kcount;
 347
 348         if(kcount>ALBE_MAXN)
 349                 printf("FATAL: neighbours!\n");
 350
 351         /* d_ik2 */
 352         d_ik2=exchange->d_ik2[kcount];
 353
 354         brand=ak->brand;
 355         if(brand==ai->brand)
 356                 S2=params->S2[brand];
 357         else
 358                 S2=params->S2mixed;
 359
 360         /* return if d_ik > S */
 361         if(d_ik2>S2) {
 362                 exchange->kcount++;
 363                 return 0;
 364         }
 365
 366         /* prefactor dzeta */
 367         pre_dzeta=exchange->pre_dzeta;
 368
 369         /* dist_ik, d_ik */
 370         dist_ik=exchange->dist_ik[kcount];
 371         d_ik=exchange->d_ik[kcount];
 372
 373         /* f_c_ik, df_c_ik */
 374         f_c_ik=exchange->f_c_ik[kcount];
 375         df_c_ik=exchange->df_c_ik[kcount];
 376
 377         /* dist_ij, d_ij2, d_ij */
 378         dist_ij=exchange->dist_ij;
 379         d_ij2=exchange->d_ij2;
 380         d_ij=exchange->d_ij;
 381
 382         /* g, dg, cos_theta */
 383         g=exchange->g[kcount];
 384         dg=exchange->dg[kcount];
 385         cos_theta=exchange->cos_theta[kcount];
 386
 387         /* cos_theta derivatives wrt i,j,k */
 388         dijdik_inv=1.0/(d_ij*d_ik);
 389         v3_scale(&dcosdrj,&dist_ik,dijdik_inv);
 390         v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
 391         v3_add(&dcosdrj,&dcosdrj,&tmp);
 392         v3_scale(&dcosdrk,&dist_ij,dijdik_inv);
 393         v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
 394         v3_add(&dcosdrk,&dcosdrk,&tmp);
 395         v3_add(&dcosdri,&dcosdrj,&dcosdrk);
 396         v3_scale(&dcosdri,&dcosdri,-1.0);
 397
 398         /* f_c_ik * dg, df_c_ik * g */
 399         fcdg=f_c_ik*dg;
 400         dfcg=df_c_ik*g;
 401
 402         /* derivative wrt i */
 403         v3_scale(&force,&dist_ik,dfcg);
 404         v3_scale(&tmp,&dcosdri,fcdg);
 405         v3_add(&force,&force,&tmp);
 406         v3_scale(&force,&force,pre_dzeta);
 407
 408         /* force contribution */
 409         v3_add(&(ai->f),&(ai->f),&force);
 410
 411 #ifdef DEBUG
 412         if(ai==&(moldyn->atom[0])) {
 413                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
 414                 printf("adding %f %f %f\n",force.x,force.y,force.z);
 415                 printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
 416         }
 417 #endif
 418
 419         /* virial */
 420         //virial_calc(ai,&force,&dist_ij);
 421
 422         /* derivative wrt j */
 423         v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
 424
 425         /* force contribution */
 426         v3_add(&(aj->f),&(aj->f),&force);
 427
 428 #ifdef DEBUG
 429         if(aj==&(moldyn->atom[0])) {
 430                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
 431                 printf("adding %f %f %f\n",force.x,force.y,force.z);
 432                 printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
 433         }
 434 #endif
 435
 436         /* virial */
 437         //v3_scale(&force,&force,-1.0);
 438         if(aj<ai)
 439                 virial_calc(ai,&force,&dist_ij);
 440
 441         /* derivative wrt k */
 442         v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
 443         v3_scale(&tmp,&dcosdrk,fcdg);
 444         v3_add(&force,&force,&tmp);
 445         v3_scale(&force,&force,pre_dzeta);
 446
 447         /* force contribution */
 448         v3_add(&(ak->f),&(ak->f),&force);
 449
 450 #ifdef DEBUG
 451         if(ak==&(moldyn->atom[0])) {
 452                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
 453                 printf("adding %f %f %f\n",force.x,force.y,force.z);
 454                 printf("total k: %f %f %f\n",ak->f.x,ak->f.y,ak->f.z);
 455         }
 456 #endif
 457
 458         /* virial */
 459         //v3_scale(&force,&force,-1.0);
 460         if(aj<ai)
 461                 virial_calc(ai,&force,&dist_ik);
 462
 463         /* increase k counter */
 464         exchange->kcount++;
 465
 466         return 0;
 467
 468 }