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