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_set_params(t_moldyn *moldyn,int element1,int element2) {
  24
  25         t_tersoff_mult_params *p;
  26
  27         // set cutoff before parameters (actually only necessary for some pots)
  28         if(moldyn->cutoff==0.0) {
  29                 printf("[tersoff] WARNING: no cutoff!\n");
  30                 return -1;
  31         }
  32
  33         /* alloc mem for potential parameters */
  34          moldyn->pot_params=malloc(sizeof(t_tersoff_mult_params));
  35          if(moldyn->pot_params==NULL) {
  36                 perror("[tersoff] pot params alloc");
  37                 return -1;
  38         }
  39
  40         /* these are now tersoff parameters */
  41         p=moldyn->pot_params;
  42
  43         // only 1 combination by now :p
  44         switch(element1) {
  45                 case SI:
  46                         /* type: silicon */
  47                         p->S[0]=TM_S_SI;
  48                         p->R[0]=TM_R_SI;
  49                         p->A[0]=TM_A_SI;
  50                         p->B[0]=TM_B_SI;
  51                         p->lambda[0]=TM_LAMBDA_SI;
  52                         p->mu[0]=TM_MU_SI;
  53                         p->beta[0]=TM_BETA_SI;
  54                         p->n[0]=TM_N_SI;
  55                         p->c[0]=TM_C_SI;
  56                         p->d[0]=TM_D_SI;
  57                         p->h[0]=TM_H_SI;
  58                         switch(element2) {
  59                                 case C:
  60                                         p->chi=TM_CHI_SIC;
  61                                         break;
  62                                 default:
  63                                         printf("[tersoff] WARNING: element2\n");
  64                                         return -1;
  65                         }
  66                         break;
  67                 default:
  68                         printf("[tersoff] WARNING: element1\n");
  69                         return -1;
  70         }
  71
  72         switch(element2) {
  73                 case C:
  74                         /* type carbon */
  75                         p->S[1]=TM_S_C;
  76                         p->R[1]=TM_R_C;
  77                         p->A[1]=TM_A_C;
  78                         p->B[1]=TM_B_C;
  79                         p->lambda[1]=TM_LAMBDA_C;
  80                         p->mu[1]=TM_MU_C;
  81                         p->beta[1]=TM_BETA_C;
  82                         p->n[1]=TM_N_C;
  83                         p->c[1]=TM_C_C;
  84                         p->d[1]=TM_D_C;
  85                         p->h[1]=TM_H_C;
  86                         break;
  87                 default:
  88                         printf("[tersoff] WARNING: element1\n");
  89                         return -1;
  90         }
  91
  92         printf("[tersoff] parameter completion\n");
  93         p->S2[0]=p->S[0]*p->S[0];
  94         p->S2[1]=p->S[1]*p->S[1];
  95         p->Smixed=sqrt(p->S[0]*p->S[1]);
  96         p->S2mixed=p->Smixed*p->Smixed;
  97         p->Rmixed=sqrt(p->R[0]*p->R[1]);
  98         p->Amixed=sqrt(p->A[0]*p->A[1]);
  99         p->Bmixed=sqrt(p->B[0]*p->B[1]);
 100         p->lambda_m=0.5*(p->lambda[0]+p->lambda[1]);
 101         p->mu_m=0.5*(p->mu[0]+p->mu[1]);
 102
 103         printf("[tersoff] mult parameter info:\n");
 104         printf("  S (A)  | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed);
 105         printf("  R (A)  | %f | %f | %f\n",p->R[0],p->R[1],p->Rmixed);
 106         printf("  A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV);
 107         printf("  B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV);
 108         printf("  lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1],
 109                                           p->lambda_m);
 110         printf("  mu     | %f | %f | %f\n",p->mu[0],p->mu[1],p->mu_m);
 111         printf("  beta   | %.10f | %.10f\n",p->beta[0],p->beta[1]);
 112         printf("  n      | %f | %f\n",p->n[0],p->n[1]);
 113         printf("  c      | %f | %f\n",p->c[0],p->c[1]);
 114         printf("  d      | %f | %f\n",p->d[0],p->d[1]);
 115         printf("  h      | %f | %f\n",p->h[0],p->h[1]);
 116         printf("  chi    | %f \n",p->chi);
 117
 118         return 0;
 119 }
 120
 121 /* tersoff 1 body part */
 122 int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) {
 123
 124         int brand;
 125         t_tersoff_mult_params *params;
 126         t_tersoff_exchange *exchange;
 127
 128         brand=ai->brand;
 129         params=moldyn->pot_params;
 130         exchange=&(params->exchange);
 131
 132         /*
 133          * simple: point constant parameters only depending on atom i to
 134          *         their right values
 135          */
 136
 137         exchange->beta_i=&(params->beta[brand]);
 138         exchange->n_i=&(params->n[brand]);
 139         exchange->c_i=&(params->c[brand]);
 140         exchange->d_i=&(params->d[brand]);
 141         exchange->h_i=&(params->h[brand]);
 142
 143         exchange->betaini=pow(*(exchange->beta_i),*(exchange->n_i));
 144         exchange->ci2=params->c[brand]*params->c[brand];
 145         exchange->di2=params->d[brand]*params->d[brand];
 146         exchange->ci2di2=exchange->ci2/exchange->di2;
 147
 148         return 0;
 149 }
 150
 151 /* tersoff 2 body part */
 152 int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 153
 154         t_tersoff_mult_params *params;
 155         t_3dvec dist_ij,force;
 156         double d_ij,d_ij2;
 157         double A,R,S,S2,lambda;
 158         double f_r,df_r;
 159         double f_c,df_c;
 160         int brand;
 161         double s_r;
 162         double arg;
 163         double energy;
 164
 165         printf("WARNING! - tersoff_mult_2bp is obsolete.\n");
 166         printf("WARNING! - repulsive part handled in 3bp/j2 routine.\n");
 167
 168         /* use newtons third law */
 169         if(ai<aj) return 0;
 170
 171         params=moldyn->pot_params;
 172         brand=aj->brand;
 173
 174         /* determine cutoff square */
 175         if(brand==ai->brand)
 176                 S2=params->S2[brand];
 177         else
 178                 S2=params->S2mixed;
 179
 180         /* dist_ij, d_ij2 */
 181         v3_sub(&dist_ij,&(aj->r),&(ai->r));
 182         if(bc) check_per_bound(moldyn,&dist_ij);
 183         d_ij2=v3_absolute_square(&dist_ij);
 184
 185         /* if d_ij2 > S2 => no force & potential energy contribution */
 186         if(d_ij2>S2) {
 187                 return 0;
 188         }
 189
 190         /* now we will need the distance */
 191         d_ij=sqrt(d_ij2);
 192
 193         /* more constants */
 194         if(brand==ai->brand) {
 195                 S=params->S[brand];
 196                 R=params->R[brand];
 197                 A=params->A[brand];
 198                 lambda=params->lambda[brand];
 199         }
 200         else {
 201                 S=params->Smixed;
 202                 R=params->Rmixed;
 203                 A=params->Amixed;
 204                 lambda=params->lambda_m;
 205         }
 206
 207         /* f_r_ij, df_r_ij */
 208         f_r=A*exp(-lambda*d_ij);
 209         df_r=lambda*f_r/d_ij;
 210
 211         /* f_c, df_c */
 212         if(d_ij<R) {
 213                 f_c=1.0;
 214                 df_c=0.0;
 215                 v3_scale(&force,&dist_ij,-df_r);
 216         }
 217         else {
 218                 s_r=S-R;
 219                 arg=M_PI*(d_ij-R)/s_r;
 220                 f_c=0.5+0.5*cos(arg);
 221                 df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
 222                 v3_scale(&force,&dist_ij,-df_c*f_r-df_r*f_c);
 223         }
 224
 225         /* add forces */
 226         v3_add(&(ai->f),&(ai->f),&force);
 227         v3_scale(&force,&force,-1.0); // reason: dri rij = - drj rij
 228         v3_add(&(aj->f),&(aj->f),&force);
 229
 230 #ifdef DEBUG
 231         if((ai==&(moldyn->atom[0]))|(aj==&(moldyn->atom[0]))) {
 232                 printf("force 2bp: [%d %d]\n",ai->tag,aj->tag);
 233                 printf("adding %f %f %f\n",force.x,force.y,force.z);
 234                 if(ai==&(moldyn->atom[0]))
 235                         printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
 236                 if(aj==&(moldyn->atom[0]))
 237                         printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
 238         }
 239 #endif
 240
 241         /* virial */
 242         virial_calc(aj,&force,&dist_ij);
 243
 244         /* energy 2bp contribution */
 245         energy=f_r*f_c;
 246         moldyn->energy+=energy;
 247         ai->e+=0.5*energy;
 248         aj->e+=0.5*energy;
 249
 250         return 0;
 251 }
 252
 253 /* tersoff 3 body potential function (first ij loop) */
 254 int tersoff_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 255
 256         t_tersoff_mult_params *params;
 257         t_tersoff_exchange *exchange;
 258         unsigned char brand;
 259         double S2;
 260         t_3dvec dist_ij;
 261         double d_ij2,d_ij;
 262
 263         params=moldyn->pot_params;
 264         exchange=&(params->exchange);
 265
 266         /* reset zeta sum */
 267         exchange->zeta_ij=0.0;
 268
 269         /*
 270          * set ij depending values
 271          */
 272
 273         brand=ai->brand;
 274
 275         if(brand==aj->brand)
 276                 S2=params->S2[brand];
 277         else
 278                 S2=params->S2mixed;
 279
 280         /* dist_ij, d_ij2 */
 281         v3_sub(&dist_ij,&(aj->r),&(ai->r));
 282         if(bc) check_per_bound(moldyn,&dist_ij);
 283         d_ij2=v3_absolute_square(&dist_ij);
 284
 285         /* if d_ij2 > S2 => no force & potential energy contribution */
 286         if(d_ij2>S2) {
 287                 moldyn->run3bp=0;
 288                 return 0;
 289         }
 290
 291         /* d_ij */
 292         d_ij=sqrt(d_ij2);
 293
 294         /* store values */
 295         exchange->dist_ij=dist_ij;
 296         exchange->d_ij2=d_ij2;
 297         exchange->d_ij=d_ij;
 298
 299         /* reset k counter for first k loop */
 300         exchange->kcount=0;
 301
 302         return 0;
 303 }
 304
 305 /* tersoff 3 body potential function (first k loop) */
 306 int tersoff_mult_3bp_k1(t_moldyn *moldyn,
 307                         t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 308
 309         t_tersoff_mult_params *params;
 310         t_tersoff_exchange *exchange;
 311         unsigned char brand;
 312         double R,S,S2;
 313         t_3dvec dist_ij,dist_ik;
 314         double d_ik2,d_ik,d_ij;
 315         double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
 316         double f_c_ik,df_c_ik;
 317         int kcount;
 318
 319         params=moldyn->pot_params;
 320         exchange=&(params->exchange);
 321         kcount=exchange->kcount;
 322
 323         if(kcount>TERSOFF_MAXN) {
 324                 printf("FATAL: neighbours = %d\n",kcount);
 325                 printf("  -> %d %d %d\n",ai->tag,aj->tag,ak->tag);
 326         }
 327
 328         /* ik constants */
 329         brand=ai->brand;
 330         if(brand==ak->brand) {
 331                 R=params->R[brand];
 332                 S=params->S[brand];
 333                 S2=params->S2[brand];
 334         }
 335         else {
 336                 R=params->Rmixed;
 337                 S=params->Smixed;
 338                 S2=params->S2mixed;
 339         }
 340
 341         /* dist_ik, d_ik2 */
 342         v3_sub(&dist_ik,&(ak->r),&(ai->r));
 343         if(bc) check_per_bound(moldyn,&dist_ik);
 344         d_ik2=v3_absolute_square(&dist_ik);
 345
 346         /* store data for second k loop */
 347         exchange->dist_ik[kcount]=dist_ik;
 348         exchange->d_ik2[kcount]=d_ik2;
 349
 350         /* return if not within cutoff */
 351         if(d_ik2>S2) {
 352                 exchange->kcount++;
 353                 return 0;
 354         }
 355
 356         /* d_ik */
 357         d_ik=sqrt(d_ik2);
 358
 359         /* dist_ij, d_ij */
 360         dist_ij=exchange->dist_ij;
 361         d_ij=exchange->d_ij;
 362
 363         /* cos theta */
 364         cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
 365
 366         /* g_ijk */
 367         h_cos=*(exchange->h_i)-cos_theta;
 368         d2_h_cos2=exchange->di2+(h_cos*h_cos);
 369         frac=exchange->ci2/d2_h_cos2;
 370         g=1.0+exchange->ci2di2-frac;
 371         dg=-2.0*frac*h_cos/d2_h_cos2;
 372
 373         /* zeta sum += f_c_ik * g_ijk */
 374         if(d_ik<=R) {
 375                 exchange->zeta_ij+=g;
 376                 f_c_ik=1.0;
 377                 df_c_ik=0.0;
 378         }
 379         else {
 380                 s_r=S-R;
 381                 arg=M_PI*(d_ik-R)/s_r;
 382                 f_c_ik=0.5+0.5*cos(arg);
 383                 df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
 384                 exchange->zeta_ij+=f_c_ik*g;
 385         }
 386
 387 #ifdef DEBUG
 388         if((ai==&(moldyn->atom[0]))|
 389            (aj==&(moldyn->atom[864]))|
 390            (ak==&(moldyn->atom[1003]))) {
 391                 printf(" -> %f %f %f\n",exchange->ci2di2,frac,h_cos);
 392         }
 393 #endif
 394
 395         /* store even more data for second k loop */
 396         exchange->g[kcount]=g;
 397         exchange->dg[kcount]=dg;
 398         exchange->d_ik[kcount]=d_ik;
 399         exchange->cos_theta[kcount]=cos_theta;
 400         exchange->f_c_ik[kcount]=f_c_ik;
 401         exchange->df_c_ik[kcount]=df_c_ik;
 402
 403         /* increase k counter */
 404         exchange->kcount++;
 405
 406         return 0;
 407 }
 408
 409 int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 410
 411         t_tersoff_mult_params *params;
 412         t_tersoff_exchange *exchange;
 413         t_3dvec force;
 414         double f_a,df_a,b,db,f_c,df_c;
 415         double f_r,df_r;
 416         double scale;
 417         double mu,B,chi;
 418         double lambda,A;
 419         double d_ij;
 420         unsigned char brand;
 421         double ni,tmp;
 422         double S,R,s_r,arg;
 423         double energy;
 424
 425         params=moldyn->pot_params;
 426         exchange=&(params->exchange);
 427
 428         brand=aj->brand;
 429         if(brand==ai->brand) {
 430                 S=params->S[brand];
 431                 R=params->R[brand];
 432                 B=params->B[brand];
 433                 A=params->A[brand];
 434                 mu=params->mu[brand];
 435                 lambda=params->lambda[brand];
 436                 chi=1.0;
 437         }
 438         else {
 439                 S=params->Smixed;
 440                 R=params->Rmixed;
 441                 B=params->Bmixed;
 442                 A=params->Amixed;
 443                 mu=params->mu_m;
 444                 lambda=params->lambda_m;
 445                 chi=params->chi;
 446         }
 447
 448         d_ij=exchange->d_ij;
 449
 450         /* f_c, df_c */
 451         if(d_ij<R) {
 452                 f_c=1.0;
 453                 df_c=0.0;
 454         }
 455         else {
 456                 s_r=S-R;
 457                 arg=M_PI*(d_ij-R)/s_r;
 458                 f_c=0.5+0.5*cos(arg);
 459                 df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
 460         }
 461
 462         /* f_a, df_a */
 463         f_a=-B*exp(-mu*d_ij);
 464         df_a=mu*f_a/d_ij;
 465
 466         /* f_r, df_r */
 467         f_r=A*exp(-lambda*d_ij);
 468         df_r=lambda*f_r/d_ij;
 469
 470         /* b, db */
 471         if(exchange->zeta_ij==0.0) {
 472                 b=chi;
 473                 db=0.0;
 474         }
 475         else {
 476                 ni=*(exchange->n_i);
 477                 tmp=exchange->betaini*pow(exchange->zeta_ij,ni-1.0);
 478                 b=(1.0+exchange->zeta_ij*tmp);
 479                 db=chi*pow(b,-1.0/(2.0*ni)-1.0);
 480                 b=db*b;
 481                 db*=-0.5*tmp;
 482         }
 483
 484         /* force contribution */
 485         scale=-0.5*(f_c*(df_r+b*df_a)+df_c*(f_r+b*df_a));
 486         v3_scale(&force,&(exchange->dist_ij),scale);
 487         v3_add(&(ai->f),&(ai->f),&force);
 488         v3_scale(&force,&force,-1.0); // dri rij = - drj rij
 489         v3_add(&(aj->f),&(aj->f),&force);
 490
 491 #ifdef DEBUG
 492         if((ai==&(moldyn->atom[0]))|(aj==&(moldyn->atom[0]))) {
 493                 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,aj->tag);
 494                 printf("adding %f %f %f\n",force.x,force.y,force.z);
 495                 if(ai==&(moldyn->atom[0]))
 496                         printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
 497                 if(aj==&(moldyn->atom[0]))
 498                         printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
 499                 printf("energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
 500                                                     f_c,b,f_a,f_r);
 501                 printf("        %f %f %f\n",exchange->zeta_ij,.0,.0);
 502         }
 503 #endif
 504
 505         /* virial */
 506         virial_calc(aj,&force,&(exchange->dist_ij));
 507
 508         /* dzeta prefactor = - 0.5 f_c f_a db */
 509         exchange->pre_dzeta=-0.5*f_a*f_c*db;
 510
 511         /* energy contribution */
 512         energy=0.5*f_c*(b*f_a+f_r);
 513         moldyn->energy+=energy;
 514         ai->e+=energy;
 515
 516         /* reset k counter for second k loop */
 517         exchange->kcount=0;
 518
 519         return 0;
 520 }
 521
 522 /* tersoff 3 body potential function (second k loop) */
 523 int tersoff_mult_3bp_k2(t_moldyn *moldyn,
 524                         t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 525
 526         t_tersoff_mult_params *params;
 527         t_tersoff_exchange *exchange;
 528         int kcount;
 529         t_3dvec dist_ik,dist_ij;
 530         double d_ik2,d_ik,d_ij2,d_ij;
 531         unsigned char brand;
 532         double S2;
 533         double g,dg,cos_theta;
 534         double pre_dzeta;
 535         double f_c_ik,df_c_ik;
 536         double dijdik_inv,fcdg,dfcg;
 537         t_3dvec dcosdri,dcosdrj,dcosdrk;
 538         t_3dvec force,tmp;
 539
 540         params=moldyn->pot_params;
 541         exchange=&(params->exchange);
 542         kcount=exchange->kcount;
 543
 544         if(kcount>TERSOFF_MAXN)
 545                 printf("FATAL: neighbours!\n");
 546
 547         /* d_ik2 */
 548         d_ik2=exchange->d_ik2[kcount];
 549
 550         brand=ak->brand;
 551         if(brand==ai->brand)
 552                 S2=params->S2[brand];
 553         else
 554                 S2=params->S2mixed;
 555
 556         /* return if d_ik > S */
 557         if(d_ik2>S2) {
 558                 exchange->kcount++;
 559                 return 0;
 560         }
 561
 562         /* prefactor dzeta */
 563         pre_dzeta=exchange->pre_dzeta;
 564
 565         /* dist_ik, d_ik */
 566         dist_ik=exchange->dist_ik[kcount];
 567         d_ik=exchange->d_ik[kcount];
 568
 569         /* f_c_ik, df_c_ik */
 570         f_c_ik=exchange->f_c_ik[kcount];
 571         df_c_ik=exchange->df_c_ik[kcount];
 572
 573         /* dist_ij, d_ij2, d_ij */
 574         dist_ij=exchange->dist_ij;
 575         d_ij2=exchange->d_ij2;
 576         d_ij=exchange->d_ij;
 577
 578         /* g, dg, cos_theta */
 579         g=exchange->g[kcount];
 580         dg=exchange->dg[kcount];
 581         cos_theta=exchange->cos_theta[kcount];
 582
 583         /* cos_theta derivatives wrt i,j,k */
 584         dijdik_inv=1.0/(d_ij*d_ik);
 585         v3_scale(&dcosdrj,&dist_ik,dijdik_inv);
 586         v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
 587         v3_add(&dcosdrj,&dcosdrj,&tmp);
 588         v3_scale(&dcosdrk,&dist_ij,dijdik_inv);
 589         v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
 590         v3_add(&dcosdrk,&dcosdrk,&tmp);
 591         v3_add(&dcosdri,&dcosdrj,&dcosdrk);
 592         v3_scale(&dcosdri,&dcosdri,-1.0);
 593
 594         /* f_c_ik * dg, df_c_ik * g */
 595         fcdg=f_c_ik*dg;
 596         dfcg=df_c_ik*g;
 597
 598         /* derivative wrt i */
 599         v3_scale(&force,&dist_ik,dfcg);
 600         v3_scale(&tmp,&dcosdri,fcdg);
 601         v3_add(&force,&force,&tmp);
 602         v3_scale(&force,&force,pre_dzeta);
 603
 604         /* force contribution */
 605         v3_add(&(ai->f),&(ai->f),&force);
 606
 607 #ifdef DEBUG
 608         if(ai==&(moldyn->atom[0])) {
 609                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
 610                 printf("adding %f %f %f\n",force.x,force.y,force.z);
 611                 printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
 612         }
 613 #endif
 614
 615         /* derivative wrt j */
 616         v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
 617
 618         /* force contribution */
 619         v3_add(&(aj->f),&(aj->f),&force);
 620
 621 #ifdef DEBUG
 622         if(aj==&(moldyn->atom[0])) {
 623                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
 624                 printf("adding %f %f %f\n",force.x,force.y,force.z);
 625                 printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
 626         }
 627 #endif
 628
 629         /* virial */
 630         v3_scale(&force,&force,-1.0);
 631         virial_calc(ai,&force,&dist_ij);
 632
 633         /* derivative wrt k */
 634         v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
 635         v3_scale(&tmp,&dcosdrk,fcdg);
 636         v3_add(&force,&force,&tmp);
 637         v3_scale(&force,&force,pre_dzeta);
 638
 639         /* force contribution */
 640         v3_add(&(ak->f),&(ak->f),&force);
 641
 642 #ifdef DEBUG
 643         if(ak==&(moldyn->atom[0])) {
 644                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
 645                 printf("adding %f %f %f\n",force.x,force.y,force.z);
 646                 printf("total k: %f %f %f\n",ak->f.x,ak->f.y,ak->f.z);
 647         }
 648 #endif
 649
 650         /* virial */
 651         v3_scale(&force,&force,-1.0);
 652         virial_calc(ai,&force,&dist_ik);
 653
 654         /* increase k counter */
 655         exchange->kcount++;
 656
 657         return 0;
 658
 659 }
 660
 661 int tersoff_mult_check_2b_bond(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 662
 663         t_tersoff_mult_params *params;
 664         t_3dvec dist;
 665         double d;
 666         u8 brand;
 667
 668         v3_sub(&dist,&(aj->r),&(ai->r));
 669         if(bc) check_per_bound(moldyn,&dist);
 670         d=v3_absolute_square(&dist);
 671
 672         params=moldyn->pot_params;
 673         brand=ai->brand;
 674
 675         if(brand==aj->brand) {
 676                 if(d<=params->S2[brand])
 677                         return TRUE;
 678         }
 679         else {
 680                 if(d<=params->S2mixed)
 681                         return TRUE;
 682         }
 683
 684         return FALSE;
 685 }
 686