potentials/tersoff_orig.c

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