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->pot1b_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_tersoff_exchange *exchange;
  89         t_3dvec dist_ij,force;
  90         double d_ij,d_ij2;
  91         double A,B,R,S,S2,lambda,mu;
  92         double f_r,df_r;
  93         double f_c,df_c;
  94         int brand;
  95         double s_r;
  96         double arg;
  97
  98         params=moldyn->pot2b_params;
  99         brand=aj->brand;
 100         exchange=&(params->exchange);
 101
 102         /* clear 3bp and 2bp post run */
 103         exchange->run3bp=0;
 104         exchange->run2bp_post=0;
 105
 106         /* reset S > r > R mark */
 107         exchange->d_ij_between_rs=0;
 108
 109         /*
 110          * calc of 2bp contribution of V_ij and dV_ij/ji
 111          *
 112          * for Vij and dV_ij we need:
 113          * - f_c_ij, df_c_ij
 114          * - f_r_ij, df_r_ij
 115          *
 116          * for dV_ji we need:
 117          * - f_c_ji = f_c_ij, df_c_ji = df_c_ij
 118          * - f_r_ji = f_r_ij; df_r_ji = df_r_ij
 119          *
 120          */
 121
 122         /* constants */
 123         if(brand==ai->brand) {
 124                 S=params->S[brand];
 125                 S2=params->S2[brand];
 126                 R=params->R[brand];
 127                 A=params->A[brand];
 128                 B=params->B[brand];
 129                 lambda=params->lambda[brand];
 130                 mu=params->mu[brand];
 131                 exchange->chi=1.0;
 132         }
 133         else {
 134                 S=params->Smixed;
 135                 S2=params->S2mixed;
 136                 R=params->Rmixed;
 137                 A=params->Amixed;
 138                 B=params->Bmixed;
 139                 lambda=params->lambda_m;
 140                 mu=params->mu_m;
 141                 params->exchange.chi=params->chi;
 142         }
 143
 144         /* dist_ij, d_ij */
 145         v3_sub(&dist_ij,&(aj->r),&(ai->r));
 146         if(bc) check_per_bound(moldyn,&dist_ij);
 147         d_ij2=v3_absolute_square(&dist_ij);
 148
 149         /* if d_ij2 > S2 => no force & potential energy contribution */
 150         if(d_ij2>S2)
 151                 return 0;
 152
 153         /* now we will need the distance */
 154         //d_ij=v3_norm(&dist_ij);
 155         d_ij=sqrt(d_ij2);
 156
 157         /* save for use in 3bp */
 158         exchange->d_ij=d_ij;
 159         exchange->d_ij2=d_ij2;
 160         exchange->dist_ij=dist_ij;
 161
 162         /* more constants */
 163         exchange->beta_j=&(params->beta[brand]);
 164         exchange->n_j=&(params->n[brand]);
 165         exchange->c_j=&(params->c[brand]);
 166         exchange->d_j=&(params->d[brand]);
 167         exchange->h_j=&(params->h[brand]);
 168         if(brand==ai->brand) {
 169                 exchange->betajnj=exchange->betaini;
 170                 exchange->cj2=exchange->ci2;
 171                 exchange->dj2=exchange->di2;
 172                 exchange->cj2dj2=exchange->ci2di2;
 173         }
 174         else {
 175                 exchange->betajnj=pow(*(exchange->beta_j),*(exchange->n_j));
 176                 exchange->cj2=params->c[brand]*params->c[brand];
 177                 exchange->dj2=params->d[brand]*params->d[brand];
 178                 exchange->cj2dj2=exchange->cj2/exchange->dj2;
 179         }
 180
 181         /* f_r_ij = f_r_ji, df_r_ij = df_r_ji */
 182         f_r=A*exp(-lambda*d_ij);
 183         df_r=lambda*f_r/d_ij;
 184
 185         /* f_a, df_a calc (again, same for ij and ji) | save for later use! */
 186         exchange->f_a=-B*exp(-mu*d_ij);
 187         exchange->df_a=mu*exchange->f_a/d_ij;
 188
 189         /* f_c, df_c calc (again, same for ij and ji) */
 190         if(d_ij<R) {
 191                 /* f_c = 1, df_c = 0 */
 192                 f_c=1.0;
 193                 df_c=0.0;
 194                 /* two body contribution (ij, ji) */
 195                 v3_scale(&force,&dist_ij,-df_r);
 196         }
 197         else {
 198                 s_r=S-R;
 199                 arg=M_PI*(d_ij-R)/s_r;
 200                 f_c=0.5+0.5*cos(arg);
 201                 df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
 202                 /* two body contribution (ij, ji) */
 203                 v3_scale(&force,&dist_ij,-df_c*f_r-df_r*f_c);
 204                 /* tell 3bp that S > r > R */
 205                 exchange->d_ij_between_rs=1;
 206         }
 207
 208         /* add forces of 2bp (ij, ji) contribution
 209          * dVij = dVji and we sum up both: no 1/2) */
 210         v3_add(&(ai->f),&(ai->f),&force);
 211
 212         /* virial */
 213         ai->virial.xx-=force.x*dist_ij.x;
 214         ai->virial.yy-=force.y*dist_ij.y;
 215         ai->virial.zz-=force.z*dist_ij.z;
 216         ai->virial.xy-=force.x*dist_ij.y;
 217         ai->virial.xz-=force.x*dist_ij.z;
 218         ai->virial.yz-=force.y*dist_ij.z;
 219
 220 #ifdef DEBUG
 221 if(ai==&(moldyn->atom[0])) {
 222         printf("dVij, dVji (2bp) contrib:\n");
 223         printf("%f | %f\n",force.x,ai->f.x);
 224         printf("%f | %f\n",force.y,ai->f.y);
 225         printf("%f | %f\n",force.z,ai->f.z);
 226 }
 227 #endif
 228 #ifdef VDEBUG
 229 if(ai==&(moldyn->atom[0])) {
 230         printf("dVij, dVji (2bp) contrib:\n");
 231         printf("%f | %f\n",force.x*dist_ij.x,ai->virial.xx);
 232         printf("%f | %f\n",force.y*dist_ij.y,ai->virial.yy);
 233         printf("%f | %f\n",force.z*dist_ij.z,ai->virial.zz);
 234 }
 235 #endif
 236
 237         /* energy 2bp contribution (ij, ji) is 0.5 f_r f_c ... */
 238         moldyn->energy+=(0.5*f_r*f_c);
 239
 240         /* save for use in 3bp */
 241         exchange->f_c=f_c;
 242         exchange->df_c=df_c;
 243
 244         /* enable the run of 3bp function and 2bp post processing */
 245         exchange->run3bp=1;
 246         exchange->run2bp_post=1;
 247
 248         /* reset 3bp sums */
 249         exchange->zeta_ij=0.0;
 250         exchange->zeta_ji=0.0;
 251         v3_zero(&(exchange->dzeta_ij));
 252         v3_zero(&(exchange->dzeta_ji));
 253
 254         return 0;
 255 }
 256
 257 /* tersoff 2 body post part */
 258
 259 int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 260
 261         /*
 262          * here we have to allow for the 3bp sums
 263          *
 264          * that is:
 265          * - zeta_ij, dzeta_ij
 266          * - zeta_ji, dzeta_ji
 267          *
 268          * to compute the 3bp contribution to:
 269          * - Vij, dVij
 270          * - dVji
 271          *
 272          */
 273
 274         t_tersoff_mult_params *params;
 275         t_tersoff_exchange *exchange;
 276
 277         t_3dvec force,temp;
 278         t_3dvec *dist_ij;
 279         double b,db,tmp;
 280         double f_c,df_c,f_a,df_a;
 281         double chi,ni,betaini,nj,betajnj;
 282         double zeta;
 283
 284         params=moldyn->pot2b_params;
 285         exchange=&(params->exchange);
 286
 287         /* we do not run if f_c_ij was detected to be 0! */
 288         if(!(exchange->run2bp_post))
 289                 return 0;
 290
 291         f_c=exchange->f_c;
 292         df_c=exchange->df_c;
 293         f_a=exchange->f_a;
 294         df_a=exchange->df_a;
 295         betaini=exchange->betaini;
 296         betajnj=exchange->betajnj;
 297         ni=*(exchange->n_i);
 298         nj=*(exchange->n_j);
 299         chi=exchange->chi;
 300         dist_ij=&(exchange->dist_ij);
 301
 302         /* Vij and dVij */
 303         zeta=exchange->zeta_ij;
 304         if(zeta==0.0) {
 305                 moldyn->debug++;                /* just for debugging ... */
 306                 b=chi;
 307                 v3_scale(&force,dist_ij,df_a*b*f_c);
 308         }
 309         else {
 310                 tmp=betaini*pow(zeta,ni-1.0);           /* beta^n * zeta^n-1 */
 311                 b=(1+zeta*tmp);                         /* 1 + beta^n zeta^n */
 312                 db=chi*pow(b,-1.0/(2*ni)-1);            /* x(...)^(-1/2n - 1) */
 313                 b=db*b;                                 /* b_ij */
 314                 db*=-0.5*tmp;                           /* db_ij */
 315                 v3_scale(&force,&(exchange->dzeta_ij),f_a*db);
 316                 v3_scale(&temp,dist_ij,df_a*b);
 317                 v3_add(&force,&force,&temp);
 318                 v3_scale(&force,&force,f_c);
 319         }
 320         v3_scale(&temp,dist_ij,df_c*b*f_a);
 321         v3_add(&force,&force,&temp);
 322         v3_scale(&force,&force,-0.5);
 323
 324         /* add force */
 325         v3_add(&(ai->f),&(ai->f),&force);
 326
 327         /* virial */
 328         ai->virial.xx-=force.x*dist_ij->x;
 329         ai->virial.yy-=force.y*dist_ij->y;
 330         ai->virial.zz-=force.z*dist_ij->z;
 331         ai->virial.xy-=force.x*dist_ij->y;
 332         ai->virial.xz-=force.x*dist_ij->z;
 333         ai->virial.yz-=force.y*dist_ij->z;
 334
 335 #ifdef DEBUG
 336 if(ai==&(moldyn->atom[0])) {
 337         printf("dVij (3bp) contrib:\n");
 338         printf("%f | %f\n",force.x,ai->f.x);
 339         printf("%f | %f\n",force.y,ai->f.y);
 340         printf("%f | %f\n",force.z,ai->f.z);
 341 }
 342 #endif
 343 #ifdef VDEBUG
 344 if(ai==&(moldyn->atom[0])) {
 345         printf("dVij (3bp) contrib:\n");
 346         printf("%f | %f\n",force.x*dist_ij->x,ai->virial.xx);
 347         printf("%f | %f\n",force.y*dist_ij->y,ai->virial.yy);
 348         printf("%f | %f\n",force.z*dist_ij->z,ai->virial.zz);
 349 }
 350 #endif
 351
 352         /* add energy of 3bp sum */
 353         moldyn->energy+=(0.5*f_c*b*f_a);
 354
 355         /* dVji */
 356         zeta=exchange->zeta_ji;
 357         if(zeta==0.0) {
 358                 moldyn->debug++;
 359                 b=chi;
 360                 v3_scale(&force,dist_ij,df_a*b*f_c);
 361         }
 362         else {
 363                 tmp=betajnj*pow(zeta,nj-1.0);           /* beta^n * zeta^n-1 */
 364                 b=(1+zeta*tmp);                         /* 1 + beta^n zeta^n */
 365                 db=chi*pow(b,-1.0/(2*nj)-1);            /* x(...)^(-1/2n - 1) */
 366                 b=db*b;                                 /* b_ij */
 367                 db*=-0.5*tmp;                           /* db_ij */
 368                 v3_scale(&force,&(exchange->dzeta_ji),f_a*db);
 369                 v3_scale(&temp,dist_ij,df_a*b);
 370                 v3_add(&force,&force,&temp);
 371                 v3_scale(&force,&force,f_c);
 372         }
 373         v3_scale(&temp,dist_ij,df_c*b*f_a);
 374         v3_add(&force,&force,&temp);
 375         v3_scale(&force,&force,-0.5);
 376
 377         /* add force */
 378         v3_add(&(ai->f),&(ai->f),&force);
 379
 380         /* virial - plus sign, as dist_ij = - dist_ji - (really??) */
 381 // TEST ... with a minus instead
 382         ai->virial.xx-=force.x*dist_ij->x;
 383         ai->virial.yy-=force.y*dist_ij->y;
 384         ai->virial.zz-=force.z*dist_ij->z;
 385         ai->virial.xy-=force.x*dist_ij->y;
 386         ai->virial.xz-=force.x*dist_ij->z;
 387         ai->virial.yz-=force.y*dist_ij->z;
 388
 389 #ifdef DEBUG
 390 if(ai==&(moldyn->atom[0])) {
 391         printf("dVji (3bp) contrib:\n");
 392         printf("%f | %f\n",force.x,ai->f.x);
 393         printf("%f | %f\n",force.y,ai->f.y);
 394         printf("%f | %f\n",force.z,ai->f.z);
 395 }
 396 #endif
 397 #ifdef VDEBUG
 398 if(ai==&(moldyn->atom[0])) {
 399         printf("dVji (3bp) contrib:\n");
 400         printf("%f | %f\n",force.x*dist_ij->x,ai->virial.xx);
 401         printf("%f | %f\n",force.y*dist_ij->y,ai->virial.yy);
 402         printf("%f | %f\n",force.z*dist_ij->z,ai->virial.zz);
 403 }
 404 #endif
 405
 406         return 0;
 407 }
 408
 409 /* tersoff 3 body part */
 410
 411 int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 412
 413         t_tersoff_mult_params *params;
 414         t_tersoff_exchange *exchange;
 415         t_3dvec dist_ij,dist_ik,dist_jk;
 416         t_3dvec temp1,temp2;
 417         t_3dvec *dzeta;
 418         double R,S,S2,s_r;
 419         double B,mu;
 420         double d_ij,d_ik,d_jk,d_ij2,d_ik2,d_jk2;
 421         double rr,dd;
 422         double f_c,df_c;
 423         double f_c_ik,df_c_ik,arg;
 424         double f_c_jk;
 425         double n,c,d,h;
 426         double c2,d2,c2d2;
 427         double cos_theta,d_costheta1,d_costheta2;
 428         double h_cos,d2_h_cos2;
 429         double frac,g,zeta,chi;
 430         double tmp;
 431         int brand;
 432
 433         params=moldyn->pot3b_params;
 434         exchange=&(params->exchange);
 435
 436         if(!(exchange->run3bp))
 437                 return 0;
 438
 439         /*
 440          * calc of 3bp contribution of V_ij and dV_ij/ji/jk &
 441          * 2bp contribution of dV_jk
 442          *
 443          * for Vij and dV_ij we still need:
 444          * - b_ij, db_ij (zeta_ij)
 445          *   - f_c_ik, df_c_ik, constants_i, cos_theta_ijk, d_costheta_ijk
 446          *
 447          * for dV_ji we still need:
 448          * - b_ji, db_ji (zeta_ji)
 449          *   - f_c_jk, d_c_jk, constants_j, cos_theta_jik, d_costheta_jik
 450          *
 451          * for dV_jk we need:
 452          * - f_c_jk
 453          * - f_a_jk
 454          * - db_jk (zeta_jk)
 455          *   - f_c_ji, df_c_ji, constants_j, cos_theta_jki, d_costheta_jki
 456          *
 457          */
 458
 459         /*
 460          * get exchange data
 461          */
 462
 463         /* dist_ij, d_ij - this is < S_ij ! */
 464         dist_ij=exchange->dist_ij;
 465         d_ij=exchange->d_ij;
 466         d_ij2=exchange->d_ij2;
 467
 468         /* f_c_ij, df_c_ij (same for ji) */
 469         f_c=exchange->f_c;
 470         df_c=exchange->df_c;
 471
 472         /*
 473          * calculate unknown values now ...
 474          */
 475
 476         /* V_ij and dV_ij stuff (in b_ij there is f_c_ik) */
 477
 478         /* dist_ik, d_ik */
 479         v3_sub(&dist_ik,&(ak->r),&(ai->r));
 480         if(bc) check_per_bound(moldyn,&dist_ik);
 481         d_ik2=v3_absolute_square(&dist_ik);
 482
 483         /* ik constants */
 484         brand=ai->brand;
 485         if(brand==ak->brand) {
 486                 R=params->R[brand];
 487                 S=params->S[brand];
 488                 S2=params->S2[brand];
 489         }
 490         else {
 491                 R=params->Rmixed;
 492                 S=params->Smixed;
 493                 S2=params->S2mixed;
 494         }
 495
 496         /* zeta_ij/dzeta_ij contribution only for d_ik < S */
 497         if(d_ik2<S2) {
 498
 499                 /* now we need d_ik */
 500                 d_ik=sqrt(d_ik2);
 501
 502                 /* get constants_i from exchange data */
 503                 n=*(exchange->n_i);
 504                 c=*(exchange->c_i);
 505                 d=*(exchange->d_i);
 506                 h=*(exchange->h_i);
 507                 c2=exchange->ci2;
 508                 d2=exchange->di2;
 509                 c2d2=exchange->ci2di2;
 510
 511                 /* cosine of theta_ijk by scalaproduct */
 512                 rr=v3_scalar_product(&dist_ij,&dist_ik);
 513                 dd=d_ij*d_ik;
 514                 cos_theta=rr/dd;
 515
 516                 /* d_costheta */
 517                 tmp=1.0/dd;
 518                 d_costheta1=cos_theta/d_ij2-tmp;
 519                 d_costheta2=cos_theta/d_ik2-tmp;
 520
 521                 /* some usefull values */
 522                 h_cos=(h-cos_theta);
 523                 d2_h_cos2=d2+(h_cos*h_cos);
 524                 frac=c2/(d2_h_cos2);
 525
 526                 /* g(cos_theta) */
 527                 g=1.0+c2d2-frac;
 528
 529                 /* d_costheta_ij and dg(cos_theta) - needed in any case! */
 530                 v3_scale(&temp1,&dist_ij,d_costheta1);
 531                 v3_scale(&temp2,&dist_ik,d_costheta2);
 532                 v3_add(&temp1,&temp1,&temp2);
 533                 v3_scale(&temp1,&temp1,-2.0*frac*h_cos/d2_h_cos2); /* dg */
 534
 535                 /* f_c_ik & df_c_ik + {d,}zeta contribution */
 536                 dzeta=&(exchange->dzeta_ij);
 537                 if(d_ik<R) {
 538                         /* {d,}f_c_ik */
 539                         // => f_c_ik=1.0;
 540                         // => df_c_ik=0.0; of course we do not set this!
 541
 542                         /* zeta_ij */
 543                         exchange->zeta_ij+=g;
 544
 545                         /* dzeta_ij */
 546                         v3_add(dzeta,dzeta,&temp1);
 547                 }
 548                 else {
 549                         /* {d,}f_c_ik */
 550                         s_r=S-R;
 551                         arg=M_PI*(d_ik-R)/s_r;
 552                         f_c_ik=0.5+0.5*cos(arg);
 553                         df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
 554
 555                         /* zeta_ij */
 556                         exchange->zeta_ij+=f_c_ik*g;
 557
 558                         /* dzeta_ij */
 559                         v3_scale(&temp1,&temp1,f_c_ik);
 560                         v3_scale(&temp2,&dist_ik,g*df_c_ik);
 561                         v3_add(&temp1,&temp1,&temp2);
 562                         v3_add(dzeta,dzeta,&temp1);
 563                 }
 564         }
 565
 566         /* dV_ji stuff (in b_ji there is f_c_jk) + dV_jk stuff! */
 567
 568         /* dist_jk, d_jk */
 569         v3_sub(&dist_jk,&(ak->r),&(aj->r));
 570         if(bc) check_per_bound(moldyn,&dist_jk);
 571         d_jk2=v3_absolute_square(&dist_jk);
 572
 573         /* jk constants */
 574         brand=aj->brand;
 575         if(brand==ak->brand) {
 576                 R=params->R[brand];
 577                 S=params->S[brand];
 578                 S2=params->S2[brand];
 579                 B=params->B[brand];
 580                 mu=params->mu[brand];
 581                 chi=1.0;
 582         }
 583         else {
 584                 R=params->Rmixed;
 585                 S=params->Smixed;
 586                 S2=params->S2mixed;
 587                 B=params->Bmixed;
 588                 mu=params->mu_m;
 589                 chi=params->chi;
 590         }
 591
 592         /* zeta_ji/dzeta_ji contribution only for d_jk < S_jk */
 593         if(d_jk2<S2) {
 594
 595                 /* now we need d_ik */
 596                 d_jk=sqrt(d_jk2);
 597
 598                 /* constants_j from exchange data */
 599                 n=*(exchange->n_j);
 600                 c=*(exchange->c_j);
 601                 d=*(exchange->d_j);
 602                 h=*(exchange->h_j);
 603                 c2=exchange->cj2;
 604                 d2=exchange->dj2;
 605                 c2d2=exchange->cj2dj2;
 606
 607                 /* cosine of theta_jik by scalaproduct */
 608                 rr=-v3_scalar_product(&dist_ij,&dist_jk); /* -1, as ij -> ji */
 609                 dd=d_ij*d_jk;
 610                 cos_theta=rr/dd;
 611
 612                 /* d_costheta */
 613                 d_costheta1=1.0/dd;
 614                 d_costheta2=cos_theta/d_ij2;
 615
 616                 /* some usefull values */
 617                 h_cos=(h-cos_theta);
 618                 d2_h_cos2=d2+(h_cos*h_cos);
 619                 frac=c2/(d2_h_cos2);
 620
 621                 /* g(cos_theta) */
 622                 g=1.0+c2d2-frac;
 623
 624                 /* d_costheta_jik and dg(cos_theta) - needed in any case! */
 625                 v3_scale(&temp1,&dist_jk,d_costheta1);
 626                 v3_scale(&temp2,&dist_ij,-d_costheta2); /* ji -> ij => -1 */
 627                 //v3_add(&temp1,&temp1,&temp2);
 628                 v3_sub(&temp1,&temp1,&temp2); /* there is a minus! */
 629                 v3_scale(&temp1,&temp1,-2.0*frac*h_cos/d2_h_cos2); /* dg */
 630
 631                 /* store dg in temp2 and use it for dVjk later */
 632                 v3_copy(&temp2,&temp1);
 633
 634                 /* f_c_jk + {d,}zeta contribution (df_c_jk = 0) */
 635                 dzeta=&(exchange->dzeta_ji);
 636                 if(d_jk<R) {
 637                         /* f_c_jk */
 638                         f_c_jk=1.0;
 639
 640                         /* zeta_ji */
 641                         exchange->zeta_ji+=g;
 642
 643                         /* dzeta_ji */
 644                         v3_add(dzeta,dzeta,&temp1);
 645                 }
 646                 else {
 647                         /* f_c_jk */
 648                         s_r=S-R;
 649                         arg=M_PI*(d_jk-R)/s_r;
 650                         f_c_jk=0.5+0.5*cos(arg);
 651
 652                         /* zeta_ji */
 653                         exchange->zeta_ji+=f_c_jk*g;
 654
 655                         /* dzeta_ji */
 656                         v3_scale(&temp1,&temp1,f_c_jk);
 657                         v3_add(dzeta,dzeta,&temp1);
 658                 }
 659
 660                 /* dV_jk stuff | add force contribution on atom i immediately */
 661                 if(exchange->d_ij_between_rs) {
 662                         zeta=f_c*g;
 663                         v3_scale(&temp1,&temp2,f_c);
 664                         v3_scale(&temp2,&dist_ij,df_c*g);
 665                         v3_add(&temp2,&temp2,&temp1); /* -> dzeta_jk in temp2 */
 666                 }
 667                 else {
 668                         zeta=g;
 669                         // dzeta_jk is simply dg, which is stored in temp2
 670                 }
 671                 /* betajnj * zeta_jk ^ nj-1 */
 672                 tmp=exchange->betajnj*pow(zeta,(n-1.0));
 673                 tmp=-chi/2.0*pow((1+tmp*zeta),(-1.0/(2.0*n)-1))*tmp;
 674                 v3_scale(&temp2,&temp2,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5);
 675                 v3_add(&(ai->f),&(ai->f),&temp2); /* -1 skipped in f_a calc ^ */
 676                                                   /* scaled with 0.5 ^ */
 677
 678                 /* virial */
 679                 ai->virial.xx-=temp2.x*dist_jk.x;
 680                 ai->virial.yy-=temp2.y*dist_jk.y;
 681                 ai->virial.zz-=temp2.z*dist_jk.z;
 682                 ai->virial.xy-=temp2.x*dist_jk.y;
 683                 ai->virial.xz-=temp2.x*dist_jk.z;
 684                 ai->virial.yz-=temp2.y*dist_jk.z;
 685
 686 #ifdef DEBUG
 687 if(ai==&(moldyn->atom[0])) {
 688         printf("dVjk (3bp) contrib:\n");
 689         printf("%f | %f\n",temp2.x,ai->f.x);
 690         printf("%f | %f\n",temp2.y,ai->f.y);
 691         printf("%f | %f\n",temp2.z,ai->f.z);
 692 }
 693 #endif
 694 #ifdef VDEBUG
 695 if(ai==&(moldyn->atom[0])) {
 696         printf("dVjk (3bp) contrib:\n");
 697         printf("%f | %f\n",temp2.x*dist_jk.x,ai->virial.xx);
 698         printf("%f | %f\n",temp2.y*dist_jk.y,ai->virial.yy);
 699         printf("%f | %f\n",temp2.z*dist_jk.z,ai->virial.zz);
 700 }
 701 #endif
 702
 703         }
 704
 705         return 0;
 706 }
 707