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