potentials/albe.c

   1 /*
   2  * albe.c - albe 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 "albe.h"
  21
  22 /* create mixed terms from parameters and set them */
  23 int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
  24
  25         t_albe_mult_params *p;
  26
  27         // set cutoff before parameters (actually only necessary for some pots)
  28         if(moldyn->cutoff==0.0) {
  29                 printf("[albe] WARNING: no cutoff!\n");
  30                 return -1;
  31         }
  32
  33         /* alloc mem for potential parameters */
  34         moldyn->pot_params=malloc(sizeof(t_albe_mult_params));
  35         if(moldyn->pot_params==NULL) {
  36                 perror("[albe] pot params alloc");
  37                 return -1;
  38         }
  39
  40         /* these are now albe 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]=ALBE_S_SI;
  48                         p->R[0]=ALBE_R_SI;
  49                         p->A[0]=ALBE_A_SI;
  50                         p->B[0]=ALBE_B_SI;
  51                         p->r0[0]=ALBE_R0_SI;
  52                         p->lambda[0]=ALBE_LAMBDA_SI;
  53                         p->mu[0]=ALBE_MU_SI;
  54                         p->gamma[0]=ALBE_GAMMA_SI;
  55                         p->c[0]=ALBE_C_SI;
  56                         p->c2[0]=p->c[0]*p->c[0];
  57                         p->d[0]=ALBE_D_SI;
  58                         p->d2[0]=p->d[0]*p->d[0];
  59                         p->c2d2[0]=p->c2[0]/p->d2[0];
  60                         p->h[0]=ALBE_H_SI;
  61                         switch(element2) {
  62                                 case C:
  63                                         /* type: carbon */
  64                                         p->S[1]=ALBE_S_C;
  65                                         p->R[1]=ALBE_R_C;
  66                                         p->A[1]=ALBE_A_C;
  67                                         p->B[1]=ALBE_B_C;
  68                                         p->r0[1]=ALBE_R0_C;
  69                                         p->lambda[1]=ALBE_LAMBDA_C;
  70                                         p->mu[1]=ALBE_MU_C;
  71                                         p->gamma[1]=ALBE_GAMMA_C;
  72                                         p->c[1]=ALBE_C_C;
  73                                         p->c2[1]=p->c[1]*p->c[1];
  74                                         p->d[1]=ALBE_D_C;
  75                                         p->d2[1]=p->d[1]*p->d[1];
  76                                         p->c2d2[1]=p->c2[1]/p->d2[1];
  77                                         p->h[1]=ALBE_H_C;
  78                                         /* mixed type: silicon carbide */
  79                                         p->Smixed=ALBE_S_SIC;
  80                                         p->Rmixed=ALBE_R_SIC;
  81                                         p->Amixed=ALBE_A_SIC;
  82                                         p->Bmixed=ALBE_B_SIC;
  83                                         p->r0_mixed=ALBE_R0_SIC;
  84                                         p->lambda_m=ALBE_LAMBDA_SIC;
  85                                         p->mu_m=ALBE_MU_SIC;
  86                                         p->gamma_m=ALBE_GAMMA_SIC;
  87                                         p->c_mixed=ALBE_C_SIC;
  88                                         p->c2_mixed=p->c_mixed*p->c_mixed;
  89                                         p->d_mixed=ALBE_D_SIC;
  90                                         p->d2_mixed=p->d_mixed*p->d_mixed;
  91                                         p->c2d2_m=p->c2_mixed/p->d2_mixed;
  92                                         p->h_mixed=ALBE_H_SIC;
  93                                         break;
  94                                 default:
  95                                         printf("[albe] WARNING: element2\n");
  96                                         return -1;
  97                         }
  98                         break;
  99                 default:
 100                         printf("[albe] WARNING: element1\n");
 101                         return -1;
 102         }
 103
 104         printf("[albe] parameter completion\n");
 105         p->S2[0]=p->S[0]*p->S[0];
 106         p->S2[1]=p->S[1]*p->S[1];
 107         p->S2mixed=p->Smixed*p->Smixed;
 108         p->c2[0]=p->c[0]*p->c[0];
 109         p->c2[1]=p->c[1]*p->c[1];
 110         p->c2_mixed=p->c_mixed*p->c_mixed;
 111         p->d2[0]=p->d[0]*p->d[0];
 112         p->d2[1]=p->d[1]*p->d[1];
 113         p->d2_mixed=p->d_mixed*p->d_mixed;
 114         p->c2d2[0]=p->c2[0]/p->d2[0];
 115         p->c2d2[1]=p->c2[1]/p->d2[1];
 116         p->c2d2_m=p->c2_mixed/p->d2_mixed;
 117
 118         printf("[albe] mult parameter info:\n");
 119         printf("  S (A)  | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed);
 120         printf("  R (A)  | %f | %f | %f\n",p->R[0],p->R[1],p->Rmixed);
 121         printf("  A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV);
 122         printf("  B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV);
 123         printf("  lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1],
 124                                           p->lambda_m);
 125         printf("  mu     | %f | %f | %f\n",p->mu[0],p->mu[1],p->mu_m);
 126         printf("  gamma  | %f | %f | %f\n",p->gamma[0],p->gamma[1],p->gamma_m);
 127         printf("  c      | %f | %f | %f\n",p->c[0],p->c[1],p->c_mixed);
 128         printf("  d      | %f | %f | %f\n",p->d[0],p->d[1],p->d_mixed);
 129         printf("  c2     | %f | %f | %f\n",p->c2[0],p->c2[1],p->c2_mixed);
 130         printf("  d2     | %f | %f | %f\n",p->d2[0],p->d2[1],p->d2_mixed);
 131         printf("  c2d2   | %f | %f | %f\n",p->c2d2[0],p->c2d2[1],p->c2d2_m);
 132         printf("  h      | %f | %f | %f\n",p->h[0],p->h[1],p->h_mixed);
 133
 134         return 0;
 135 }
 136
 137 /* first i loop */
 138 int albe_mult_i0(t_moldyn *moldyn,t_atom *ai) {
 139
 140         t_albe_mult_params *params;
 141         t_albe_exchange *exchange;
 142
 143         int i;
 144
 145         params=moldyn->pot_params;
 146         exchange=&(params->exchange);
 147
 148         /* zero exchange values */
 149         memset(exchange->zeta,0,ALBE_MAXN*sizeof(double));
 150         for(i=0;i<ALBE_MAXN;i++)
 151                 memset(exchange->dzeta[i],0,ALBE_MAXN*sizeof(t_3dvec));
 152         exchange->jcnt=0;
 153         exchange->j2cnt=0;
 154
 155         return 0;
 156 }
 157
 158 /* first j loop within first i loop */
 159 int albe_mult_i0_j0(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 160
 161         t_albe_mult_params *params;
 162         t_albe_exchange *exchange;
 163
 164         double S2,S,R,d2,d,s_r,arg;
 165         t_3dvec dist;
 166         int j;
 167         u8 brand;
 168
 169         params=moldyn->pot_params;
 170         exchange=&(params->exchange);
 171
 172         /* get j counter */
 173         j=exchange->jcnt;
 174
 175         /* set ij depending values */
 176         brand=ai->brand;
 177         if(brand==aj->brand) {
 178                 S2=params->S2[brand];
 179         }
 180         else {
 181                 S2=params->S2mixed;
 182         }
 183
 184         /* dist_ij, d_ij2 */
 185         v3_sub(&dist,&(aj->r),&(ai->r));
 186         if(bc) check_per_bound(moldyn,&dist);
 187         exchange->dist[j]=dist;
 188         d2=v3_absolute_square(&dist);
 189         exchange->d2[j]=d2;
 190
 191         /* if d_ij2 > S2 => no force & potential energy contribution */
 192         if(d2>S2) {
 193                 moldyn->run3bp=0;
 194                 exchange->skip[j]=1;
 195                 exchange->jcnt+=1;
 196                 return 0;
 197         }
 198         exchange->skip[j]=0;
 199
 200         /* more ij depending values */
 201         if(brand==aj->brand) {
 202                 R=params->R[brand];
 203                 S=params->S[brand];
 204                 /* albe needs i,(j/k) depending c,d,h and gamma values */
 205                 exchange->gamma_[j]=&(params->gamma[brand]);
 206                 exchange->c_[j]=&(params->c[brand]);
 207                 exchange->d_[j]=&(params->d[brand]);
 208                 exchange->h_[j]=&(params->h[brand]);
 209                 exchange->c2_[j]=&(params->c2[brand]);
 210                 exchange->d2_[j]=&(params->d2[brand]);
 211                 exchange->c2d2_[j]=&(params->c2d2[brand]);
 212         }
 213         else {
 214                 R=params->Rmixed;
 215                 S=params->Smixed;
 216                 /* albe needs i,(j/k) depending c,d,h and gamma values */
 217                 exchange->gamma_[j]=&(params->gamma_m);
 218                 exchange->c_[j]=&(params->c_mixed);
 219                 exchange->d_[j]=&(params->d_mixed);
 220                 exchange->h_[j]=&(params->h_mixed);
 221                 exchange->c2_[j]=&(params->c2_mixed);
 222                 exchange->d2_[j]=&(params->d2_mixed);
 223                 exchange->c2d2_[j]=&(params->c2d2_m);
 224         }
 225
 226         /* d_ij */
 227         d=sqrt(d2);
 228         exchange->d[j]=d;
 229
 230         /* f_c, df_c */
 231         if(d<R) {
 232                 exchange->f_c[j]=1.0;
 233                 exchange->df_c[j]=0.0;
 234         }
 235         else {
 236                 s_r=S-R;
 237                 arg=M_PI*(d-R)/s_r;
 238                 exchange->f_c[j]=0.5+0.5*cos(arg);
 239                 exchange->df_c[j]=0.5*sin(arg)*(M_PI/(s_r*d));
 240         }
 241
 242         /* reset k counter */
 243         exchange->kcnt=0;
 244
 245         return 0;
 246 }
 247
 248 /* first k loop within first j loop within first i loop */
 249 int albe_mult_i0_j0_k0(t_moldyn *moldyn,
 250                        t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 251
 252         t_albe_mult_params *params;
 253         t_albe_exchange *exchange;
 254
 255         int j,k;
 256         t_3dvec distj,distk;
 257         double dj,dk,djdk_inv,cos_theta;
 258         double gj,dgj,h_cos_j,d2_h_cos2_j,frac_j;
 259         double gk,dgk,h_cos_k,d2_h_cos2_k,frac_k;
 260         t_3dvec dcosdrj,dcosdrk,tmp;
 261         t_3dvec *dzjj,*dzkk,*dzjk,*dzkj;
 262
 263         params=moldyn->pot_params;
 264         exchange=&(params->exchange);
 265
 266         if(aj==ak) {
 267                 exchange->kcnt+=1;
 268                 return 0;
 269         }
 270
 271         /* k<j & check whether to run k */
 272         j=exchange->jcnt;
 273         k=exchange->kcnt;
 274         if(k>=ALBE_MAXN) {
 275                 printf("FATAL: too many neighbours! (%d)\n",k);
 276                 printf("  atom i:%d | j:%d | k:%d\n",ai->tag,aj->tag,ak->tag);
 277         }
 278         if((k>=j)|(exchange->skip[k])) {
 279                 exchange->kcnt+=1;
 280                 return 0;
 281         }
 282
 283         /* distances */
 284         distj=exchange->dist[j];
 285         distk=exchange->dist[k];
 286         dj=exchange->d[j];
 287         dk=exchange->d[k];
 288         djdk_inv=1.0/(dj*dk);
 289
 290         /* cos theta */
 291         cos_theta=v3_scalar_product(&distj,&distk)*djdk_inv;
 292
 293         /* g(cos(theta)) ij and ik values */
 294         h_cos_j=*(exchange->h_[j])+cos_theta; // + in albe formalism
 295         d2_h_cos2_j=*exchange->d2_[j]+(h_cos_j*h_cos_j);
 296         frac_j=*exchange->c2_[j]/d2_h_cos2_j;
 297         gj=1.0+*exchange->c2d2_[j]-frac_j;
 298         gj*=*(exchange->gamma_[j]);
 299         dgj=*(exchange->gamma_[j])*2.0*frac_j*h_cos_j/d2_h_cos2_j; // + in albe
 300         if(ak->brand==aj->brand) {
 301                 gk=gj;
 302                 dgk=dgj;
 303         }
 304         else {
 305                 h_cos_k=*(exchange->h_[k])+cos_theta;
 306                 d2_h_cos2_k=*exchange->d2_[k]+(h_cos_k*h_cos_k);
 307                 frac_k=*exchange->c2_[k]/d2_h_cos2_k;
 308                 gk=1.0+*exchange->c2d2_[k]-frac_k;
 309                 gk*=*(exchange->gamma_[k]);
 310                 dgk=*(exchange->gamma_[k])*2.0*frac_k*h_cos_k/d2_h_cos2_k;
 311         }
 312
 313 #ifdef DEBUG
 314         if(ai==&(moldyn->atom[DATOM]))
 315                 printf("zeta_ij: %f %f %f %f\n",f_c_ik*g,f_c_ik,g,d_ik);
 316 #endif
 317
 318         /* store even more data for second k loop */
 319         exchange->g[kcount]=g;
 320         exchange->dg[kcount]=dg;
 321         exchange->d_ik[kcount]=d_ik;
 322         exchange->cos_theta[kcount]=cos_theta;
 323         exchange->f_c_ik[kcount]=f_c_ik;
 324         exchange->df_c_ik[kcount]=df_c_ik;
 325
 326         /* increase k counter */
 327         exchange->kcnt+=1;
 328
 329         return 0;
 330 }
 331
 332 /* first j loop within first i loop */
 333 int albe_mult_i0_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 334
 335         t_albe_mult_params *params;
 336         t_albe_exchange *exchange;
 337
 338         params=moldyn->pot_params;
 339         exchange=&(params->exchange);
 340
 341         /* increase j counter */
 342         exchange->jcnt+=1;
 343
 344         return 0;
 345 }
 346
 347 /* second j loop within first i loop */
 348 int albe_mult_i0_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 349
 350         t_albe_mult_params *params;
 351         t_albe_exchange *exchange;
 352
 353         int j;
 354         double d,f_a,df_a,f_r,df_r,f_c,df_c,b,db;
 355         double A,B,mu,lambda,r0;
 356         double energy;
 357         t_3dvec *dist,force;
 358         double scale;
 359         u8 brand;
 360
 361         params=moldyn->pot_params;
 362         exchange=&(params->exchange);
 363
 364         /* get j counter */
 365         j=exchange->j2cnt;
 366
 367         /* skip if j not within cutoff */
 368         if(exchange->skip[j]) {
 369                 moldyn->run3bp=0;
 370                 exchange->j2cnt+=1;
 371                 return 0;
 372         }
 373
 374         /* distance */
 375         d=exchange->d[j];
 376         dist=&(exchange->dist[j]);
 377         f_c=exchange->f_c[j];
 378         df_c=exchange->df_c[j];
 379
 380         /* determine parameters to calculate fa, dfa, fr, dfr */
 381         brand=aj->brand;
 382         if(brand==ai->brand) {
 383                 B=params->B[brand];
 384                 A=params->A[brand];
 385                 r0=params->r0[brand];
 386                 mu=params->mu[brand];
 387                 lambda=params->lambda[brand];
 388         }
 389         else {
 390                 B=params->Bmixed;
 391                 A=params->Amixed;
 392                 r0=params->r0_mixed;
 393                 mu=params->mu_m;
 394                 lambda=params->lambda_m;
 395         }
 396
 397         /* f_a, df_a */
 398         f_a=-B*exp(-mu*(d-r0));
 399         df_a=mu*f_a/d;
 400
 401         /* f_r, df_r */
 402         f_r=A*exp(-lambda*(d-r0));
 403         df_r=lambda*f_r/d;
 404
 405         /* b, db */
 406         b=1.0/sqrt(1.0+exchange->zeta[j]);
 407         db=-0.5*b/(1.0+exchange->zeta[j]);
 408
 409         /* energy contribution */
 410         energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
 411         moldyn->energy+=energy;
 412         ai->e+=energy;
 413
 414         /* force contribution for atom i due to ij bond */
 415         scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
 416         v3_scale(&force,dist,scale);
 417         v3_add(&(ai->f),&(ai->f),&force);
 418
 419 #ifdef NDEBUG
 420 if(ai->tag==0) {
 421 printf("force: %.15f %.15f %.15f | %d %d (ij) %.15f\n",force.x,force.y,force.z,ai->tag,aj->tag,exchange->zeta[j]);
 422 printf("    t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z);
 423 }
 424 #endif
 425
 426         /* force contribution for atom j due to ij bond */
 427         v3_scale(&force,&force,-1.0); // dri rij = - drj rij
 428         v3_add(&(aj->f),&(aj->f),&force);
 429
 430         /* virial */
 431         virial_calc(ai,&force,&(exchange->dist_ij));
 432
 433 #ifdef DEBUG
 434         if((ai==&(moldyn->atom[DATOM]))|(aj==&(moldyn->atom[DATOM]))) {
 435                 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,aj->tag);
 436                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
 437                 if(ai==&(moldyn->atom[DATOM]))
 438                         printf("  total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
 439                 if(aj==&(moldyn->atom[DATOM]))
 440                         printf("  total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
 441                 printf("  energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
 442                                                     f_c,b,f_a,f_r);
 443                 printf("          %f %f %f\n",exchange->zeta_ij,.0,.0);
 444         }
 445 #endif
 446
 447         /* virial */
 448         virial_calc(ai,&force,dist);
 449
 450         /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
 451         exchange->pre_dzeta=0.5*f_a*f_c*db;
 452
 453         /* force contribution (drj derivative) */
 454         v3_scale(&force,&(exchange->dzeta[j][j]),exchange->pre_dzeta);
 455         v3_add(&(aj->f),&(aj->f),&force);
 456
 457 #ifdef NDEBUG
 458 if(aj->tag==0) {
 459 printf("force: %.15f %.15f %.15f | %d %d (j der)\n",force.x,force.y,force.z,aj->tag,ai->tag);
 460 printf("    t: %.15f %.15f %.15f\n",aj->f.x,aj->f.y,aj->f.z);
 461 }
 462 #endif
 463
 464         /* virial */
 465         virial_calc(ai,&force,dist);
 466
 467         v3_scale(&force,&force,-1.0);
 468         v3_add(&(ai->f),&(ai->f),&force);
 469
 470 #ifdef NDEBUG
 471 if(ai->tag==0) {
 472 printf("force: %.15f %.15f %.15f | %d %d (i contr j der)\n",force.x,force.y,force.z,ai->tag,aj->tag);
 473 printf("    t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z);
 474 }
 475 #endif
 476
 477         /* reset k counter for second k loop */
 478         exchange->kcnt=0;
 479
 480         return 0;
 481 }
 482
 483 /* second k loop within second j loop within first i loop */
 484 int albe_mult_i0_j2_k0(t_moldyn *moldyn,
 485                        t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 486
 487         t_albe_mult_params *params;
 488         t_albe_exchange *exchange;
 489
 490         int j,k;
 491         t_3dvec force;
 492
 493         params=moldyn->pot_params;
 494         exchange=&(params->exchange);
 495
 496         if(aj==ak) {
 497                 exchange->kcnt+=1;
 498                 return 0;
 499         }
 500
 501         /* prefactor dzeta */
 502         pre_dzeta=exchange->pre_dzeta;
 503
 504         /* dist_ik, d_ik */
 505         dist_ik=exchange->dist_ik[kcount];
 506         d_ik=exchange->d_ik[kcount];
 507
 508         /* f_c_ik, df_c_ik */
 509         f_c_ik=exchange->f_c_ik[kcount];
 510         df_c_ik=exchange->df_c_ik[kcount];
 511
 512         /* dist_ij, d_ij2, d_ij */
 513         dist_ij=exchange->dist_ij;
 514         d_ij2=exchange->d_ij2;
 515         d_ij=exchange->d_ij;
 516
 517         /* g, dg, cos_theta */
 518         g=exchange->g[kcount];
 519         dg=exchange->dg[kcount];
 520         cos_theta=exchange->cos_theta[kcount];
 521
 522         /* cos_theta derivatives wrt j,k */
 523         dijdik_inv=1.0/(d_ij*d_ik);
 524         v3_scale(&dcosdrj,&dist_ik,dijdik_inv);         // j
 525         v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
 526         v3_add(&dcosdrj,&dcosdrj,&tmp);
 527         v3_scale(&dcosdrk,&dist_ij,dijdik_inv);         // k
 528         v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
 529         v3_add(&dcosdrk,&dcosdrk,&tmp);
 530
 531         /* f_c_ik * dg, df_c_ik * g */
 532         fcdg=f_c_ik*dg;
 533         dfcg=df_c_ik*g;
 534
 535         /* derivative wrt j */
 536         v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
 537
 538         /* force contribution */
 539         v3_add(&(aj->f),&(aj->f),&force);
 540
 541 #ifdef DEBUG
 542         if(aj==&(moldyn->atom[DATOM])) {
 543                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
 544                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
 545                 printf("  total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
 546                 printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
 547                 printf("    d ij ik = %f %f\n",d_ij,d_ik);
 548         }
 549 #endif
 550
 551         /* virial */
 552         virial_calc(ai,&force,&dist_ij);
 553
 554         /* force contribution to atom i */
 555         v3_scale(&force,&force,-1.0);
 556         v3_add(&(ai->f),&(ai->f),&force);
 557
 558         /* derivative wrt k */
 559         v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
 560         v3_scale(&tmp,&dcosdrk,fcdg);
 561         v3_add(&force,&force,&tmp);
 562         v3_scale(&force,&force,pre_dzeta);
 563
 564         v3_scale(&force,&force,-1.0);
 565         v3_add(&(ai->f),&(ai->f),&force);
 566
 567 #ifdef DEBUG
 568         if(ak==&(moldyn->atom[DATOM])) {
 569                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
 570                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
 571                 printf("  total k: %f %f %f\n",ak->f.x,ak->f.y,ak->f.z);
 572                 printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
 573                 printf("    d ij ik = %f %f\n",d_ij,d_ik);
 574         }
 575 #endif
 576
 577         /* virial */
 578         virial_calc(ai,&force,&dist_ik);
 579
 580         /* force contribution to atom i */
 581         v3_scale(&force,&force,-1.0);
 582         v3_add(&(ai->f),&(ai->f),&force);
 583
 584         /* increase k counter */
 585         exchange->kcnt+=1;
 586
 587         return 0;
 588 }
 589
 590 int albe_mult_i0_j3(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 591
 592         t_albe_mult_params *params;
 593         t_albe_exchange *exchange;
 594
 595         params=moldyn->pot_params;
 596         exchange=&(params->exchange);
 597
 598         /* increase j counter */
 599         exchange->j2cnt+=1;
 600
 601         return 0;
 602 }
 603
 604 int albe_mult_check_2b_bond(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,u8 bc) {
 605
 606         t_albe_mult_params *params;
 607         t_3dvec dist;
 608         double d;
 609         u8 brand;
 610
 611         v3_sub(&dist,&(jtom->r),&(itom->r));
 612         if(bc) check_per_bound(moldyn,&dist);
 613         d=v3_absolute_square(&dist);
 614
 615         params=moldyn->pot_params;
 616         brand=itom->brand;
 617
 618         if(brand==jtom->brand) {
 619                 if(d<=params->S2[brand])
 620                         return TRUE;
 621         }
 622         else {
 623                 if(d<=params->S2mixed)
 624                         return TRUE;
 625         }
 626
 627         return FALSE;
 628 }