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