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