MolLib/MD/MDForce.c

   1 /*
   2  *  MDForce.c
   3  *
   4  *  Created by Toshi Nagata on 2005/06/07.
   5  *  Copyright 2005 Toshi Nagata. All rights reserved.
   6  *
   7  This program is free software; you can redistribute it and/or modify
   8  it under the terms of the GNU General Public License as published by
   9  the Free Software Foundation version 2 of the License.
  10
  11  This program is distributed in the hope that it will be useful,
  12  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  GNU General Public License for more details.
  15  */
  16
  17 #include "MDForce.h"
  18 #include "MDGraphite.h"
  19 #include "MDEwald.h"
  20
  21 #include <stdlib.h>
  22 #include <string.h>
  23
  24 extern int do_custom_bond_callback(MDArena *arena, Double r, void *procObj, Double *energy, Double *force);
  25
  26 static void
  27 s_custom_bond_force(MDArena *arena, Double r, MDCustomBondPar *cp, Double *energy, Double *force)
  28 {
  29         float s;
  30         switch (cp->type) {
  31                 case kMorseType:
  32                         s = exp(-cp->u.morse.a * (r - cp->u.morse.r0));
  33                         *energy = cp->u.morse.D * (1.0 + s * (s - 2.0));
  34                         *force = -2.0 * cp->u.morse.a * cp->u.morse.D * s * (1.0 - s);
  35                         break;
  36                 case kTCLProcType:
  37                 /*      retval = do_custom_bond_callback(arena, r, cp->u.proc, energy, force);
  38                         if (retval != 0)
  39                                 md_panic(arena, NULL); */
  40                         break;
  41                 default:
  42                         *energy = *force = 0.0;
  43                         break;
  44         }
  45 }
  46
  47 static void
  48 s_calc_bond_force(MDArena *arena)
  49 {
  50         Atom *ap = arena->mol->atoms;
  51         Int nbonds = arena->mol->nbonds;
  52         Int *bonds = arena->mol->bonds;
  53         Int *bond_par_i = arena->bond_par_i;
  54         Int *custom_bond_par_i = arena->custom_bond_par_i;
  55         BondPar *bond_pars = arena->par->bondPars;
  56         Double *anbond_r0 = arena->anbond_r0;
  57         Double *energies = &arena->energies[kBondIndex];
  58         Vector *forces = &arena->forces[kBondIndex * arena->mol->natoms];
  59         int i;
  60         for (i = 0; i < nbonds; i++, bonds += 2, bond_par_i++) {
  61                 Vector r12;
  62                 Double k0, k1, w1, w2, r0;
  63                 Int idx;
  64                 BondPar *bp;
  65                 if (*bond_par_i < 0)
  66                         continue;  /*  Ignore this entry  */
  67         /*      if (bond_uniq != NULL && bond_uniq[i] >= 0)
  68                         continue;  *//*  Non-unique bond  */
  69                 if (ap[bonds[0]].mm_exclude || ap[bonds[1]].mm_exclude)
  70                         continue;  /*  Skip non-occupied atoms  */
  71                 VecSub(r12, ap[bonds[1]].r, ap[bonds[0]].r);
  72                 w1 = VecLength(r12);
  73                 if (custom_bond_par_i != NULL && (idx = custom_bond_par_i[i]) > 0 && idx <= arena->ncustom_bond_pars) {
  74                         Double energy, force;
  75                         s_custom_bond_force(arena, w1, arena->custom_bond_pars + (idx - 1), &energy, &force);
  76                         k0 = energy;
  77                         k1 = force / w1;
  78                         r0 = 0.0;  /*  To keep complier happy  */
  79                 } else {
  80                         bp = bond_pars + *bond_par_i;
  81                         r0 = bp->r0;
  82                         if (anbond_r0 != NULL) {
  83                                 if (anbond_r0[i] > 0.0)
  84                                         r0 += anbond_r0[i];
  85                         }
  86                         w2 = w1 - r0;
  87                         k0 = bp->k * w2 * w2;         /*  Energy  */
  88                         k1 = -2.0 * bp->k * w2 / w1;  /*  Force / r  */
  89                 }
  90                 VecScaleSelf(r12, k1);
  91                 *energies += k0;
  92                 VecDec(forces[bonds[0]], r12);
  93                 VecInc(forces[bonds[1]], r12);
  94                 if (arena->debug_result && arena->debug_output_level > 1) {
  95                         fprintf(arena->debug_result, "bond force %d-%d: r=%f, r0=%f, k1=%f, {%f %f %f}\n", bonds[0]+1, bonds[1]+1, w1, r0, k1 * INTERNAL2KCAL, r12.x * INTERNAL2KCAL, r12.y * INTERNAL2KCAL, r12.z * INTERNAL2KCAL);
  96                 }
  97         }
  98 }
  99
 100 static inline int
 101 s_calc_angle_force_one(const Vector *r1, const Vector *r2, const Vector *r3, const AnglePar *anp, Double *angle, Double *energy, Double *k1out, Vector *force1, Vector *force2, Vector *force3)
 102 {
 103         Vector r21, r23, f1, f2, v1;
 104         Double w1, w2, w3, cost, sint, t, k0, k1;
 105         VecSub(r21, *r1, *r2);
 106         VecSub(r23, *r3, *r2);
 107         w1 = 1.0 / VecLength(r21);
 108         w2 = 1.0 / VecLength(r23);
 109         VecScaleSelf(r21, w1);
 110         VecScaleSelf(r23, w2);
 111         cost = VecDot(r21, r23);
 112         if (cost > 0.999999 || cost < -0.999999) {
 113                 /*      printf("Cannot handle linear angle %d-%d-%d: skipped.\n", angles[0]+1, angles[1]+1, angles[2]+1); */
 114                 return -1;  /*  Cannot handle linear angle  */
 115         }
 116         sint = sqrt(1.0 - cost * cost);
 117         /*      if (sint < 1e-5 && sint > -1e-5)
 118          continue;  *//*  Cannot handle linear angle  */
 119         t = atan2(sint, cost) - anp->a0;
 120         k0 = anp->k * t * t;
 121
 122         if (sint < 0.1 && sint > -0.1) {
 123                 /* ---- Method 1 ---- */
 124                 /* This is slower than method 2, but it is safer when sint is close to 0 */
 125                 k1 = -2 * anp->k * t;
 126                 VecCross(v1, r21, r23);
 127                 VecCross(f1, r21, v1);
 128                 w3 = w1 * k1 / VecLength(f1);
 129                 if (!isfinite(w3))
 130                         return -1;
 131                 VecScaleSelf(f1, w3);
 132                 VecCross(f2, v1, r23);
 133                 w3 = w2 * k1 / VecLength(f2);
 134                 if (!isfinite(w3))
 135                         return -1;
 136                 VecScaleSelf(f2, w3);
 137         } else {
 138                 /* ---- Method 2 ---- */
 139                 k1 = -2 * anp->k * t / sint;
 140                 VecScale(f1, r21, cost);
 141                 VecDec(f1, r23);
 142                 w3 = w1 * k1;
 143                 VecScaleSelf(f1, w3);
 144                 VecScale(f2, r23, cost);
 145                 VecDec(f2, r21);
 146                 w3 = w2 * k1;
 147                 VecScaleSelf(f2, w3);
 148         }
 149
 150         *energy = k0;
 151         *force1 = f1;
 152         *force3 = f2;
 153         *angle = t + anp->a0;
 154         *k1out = k1;
 155         VecInc(f1, f2);
 156         VecScale(*force2, f1, -1);
 157         return 0;
 158 }
 159
 160 static void
 161 s_calc_angle_force(MDArena *arena)
 162 {
 163         Atom *ap = arena->mol->atoms;
 164         Int nangles = arena->mol->nangles;
 165         Int *angles = arena->mol->angles;
 166         Int *angle_par_i = arena->angle_par_i;
 167         AnglePar *angle_pars = arena->par->anglePars;
 168 /*      Int *angle_uniq = (arena->nsyms > 0 ? arena->sym_angle_uniq : NULL); */
 169         Double *energies = &arena->energies[kAngleIndex];
 170         Vector *forces = &arena->forces[kAngleIndex * arena->mol->natoms];
 171         int i;
 172         for (i = 0; i < nangles; i++, angles += 3, angle_par_i++) {
 173                 Atom *ap1, *ap2, *ap3;
 174                 AnglePar *anp;
 175                 Double en, ang, k1;
 176                 Vector f1, f2, f3;
 177                 if (*angle_par_i < 0)
 178                         continue;  /*  Ignore this entry  */
 179                 ap1 = ATOM_AT_INDEX(ap, angles[0]);
 180                 ap2 = ATOM_AT_INDEX(ap, angles[1]);
 181                 ap3 = ATOM_AT_INDEX(ap, angles[2]);
 182         /*      if (angle_uniq != NULL && angle_uniq[i] >= 0)
 183                         continue; */ /*  Non-unique angle  */
 184                 if (ap1->mm_exclude || ap2->mm_exclude || ap3->mm_exclude)
 185                         continue;  /*  Skip non-occupied atoms  */
 186                 if (arena->nalchem_flags > 0) {
 187                         if (angles[0] < arena->nalchem_flags && angles[1] < arena->nalchem_flags
 188                                 && angles[2] < arena->nalchem_flags
 189                                 && (arena->alchem_flags[angles[0]] | arena->alchem_flags[angles[1]] | arena->alchem_flags[angles[2]]) == 3)
 190                                 continue;  /*  Interaction between vanishing and appearing groups is ignored  */
 191                 }
 192                 anp = angle_pars + *angle_par_i;
 193                 if (s_calc_angle_force_one(&(ap1->r), &(ap2->r), &(ap3->r), angle_pars + *angle_par_i, &ang, &en, &k1, &f1, &f2, &f3) != 0)
 194                         continue;
 195                 *energies += en;
 196                 VecInc(forces[angles[0]], f1);
 197                 VecInc(forces[angles[1]], f2);
 198                 VecInc(forces[angles[2]], f3);
 199                 if (arena->debug_result && arena->debug_output_level > 1) {
 200                         fprintf(arena->debug_result, "angle force %d-%d-%d: a=%f, a0=%f, k1=%f, {%f %f %f}, {%f %f %f}\n", angles[0]+1, angles[1]+1, angles[2]+1, ang*180/PI, anp->a0*180/PI, k1 * INTERNAL2KCAL, -f2.x * INTERNAL2KCAL, -f2.y * INTERNAL2KCAL, -f2.z * INTERNAL2KCAL, f3.x * INTERNAL2KCAL, f3.y * INTERNAL2KCAL, f3.z * INTERNAL2KCAL);
 201                 }
 202         }
 203 }
 204
 205 static inline int
 206 s_calc_dihedral_force_one(const Vector *r1, const Vector *r2, const Vector *r3, const Vector *r4, const TorsionPar *tp, Double *phi, Double *energy, Double *k1out, Vector *force1, Vector *force2, Vector *force3, Vector *force4)
 207 {
 208         Vector r21, r32, r43;
 209         Vector v1, v2, v3;
 210         Double w1, w2, w3, k0, k1;
 211         Double cosphi, sinphi;
 212         Vector f1, f2, f3;
 213         Int n;
 214         VecSub(r21, *r1, *r2);
 215         VecSub(r32, *r2, *r3);
 216         VecSub(r43, *r3, *r4);
 217         VecCross(v1, r21, r32);
 218         VecCross(v2, r32, r43);
 219         VecCross(v3, r32, v1);
 220         w1 = VecLength(v1);
 221         w2 = VecLength(v2);
 222         w3 = VecLength(v3);
 223         if (w1 < 1e-5 || w2 < 1e-5 || w3 < 1e-5)
 224                 return -1;  /*  The dihedral cannot be defined  */
 225         w1 = 1.0/w1;
 226         w2 = 1.0/w2;
 227         w3 = 1.0/w3;
 228         VecScaleSelf(v1, w1);
 229         VecScaleSelf(v2, w2);
 230         VecScaleSelf(v3, w3);
 231
 232         /*  The dihedral angle value  */
 233         cosphi = VecDot(v1, v2);
 234         sinphi = VecDot(v3, v2);
 235         *phi = -atan2(sinphi, cosphi);
 236
 237         /*  Repeat for multiple dihedral terms  */
 238         k0 = k1 = 0.0;
 239         for (n = 0; n < tp->mult; n++) {
 240                 Double k = tp->k[n];
 241                 Double phi0 = tp->phi0[n];
 242                 Int period = tp->period[n];
 243                 if (period > 0) {
 244                         k0 += k * (1 + cos(period * *phi + phi0));
 245                         k1 -= period * k * sin(period * *phi + phi0);
 246                 } else {
 247                         /*  A simple quadratic term  */
 248                         phi0 = *phi - phi0;
 249                         if (phi0 < -PI)
 250                                 phi0 += 2 * PI;
 251                         else if (phi0 > PI)
 252                                 phi0 -= 2 * PI;
 253                         k0 += k * phi0 * phi0;
 254                         k1 += 2 * k * phi0;
 255                 }
 256         }
 257
 258         if (sinphi < -0.1 || sinphi > 0.1) {
 259                 /*  The sin form  */
 260                 Vector v4, v5, vw0;
 261                 VecScale(v4, v1, cosphi);
 262                 VecDec(v4, v2);
 263                 VecScaleSelf(v4, w1);
 264                 VecScale(v5, v2, cosphi);
 265                 VecDec(v5, v1);
 266                 VecScaleSelf(v5, w2);
 267                 k1 /= sinphi;
 268                 VecCross(f1, r32, v4);
 269                 VecCross(f3, v5, r32);
 270                 VecCross(f2, v4, r21);
 271                 VecCross(vw0, r43, v5);
 272                 VecInc(f2, vw0);
 273                 VecScaleSelf(f1, k1);
 274                 VecScaleSelf(f2, k1);
 275                 VecScaleSelf(f3, k1);
 276         } else {
 277                 /*  The cos form  */
 278                 Vector v6, v7, vw1, vw2;
 279                 VecScale(v6, v3, sinphi);
 280                 VecDec(v6, v2);
 281                 VecScaleSelf(v6, w3);
 282                 VecScale(v7, v2, sinphi);
 283                 VecDec(v7, v3);
 284                 VecScaleSelf(v7, w2);
 285                 k1 = -k1 / cosphi;
 286                 VecCross(vw1, v6, r32);
 287                 VecCross(f1, r32, vw1);
 288                 VecCross(f3, v7, r32);
 289                 VecCross(f2, r43, v7);
 290                 VecCross(vw1, r32, r21);
 291                 VecCross(vw2, v6, vw1);
 292                 VecInc(f2, vw2);
 293                 VecCross(vw1, r32, v6);
 294                 VecCross(vw2, r21, vw1);
 295                 VecInc(f2, vw2);
 296                 VecScaleSelf(f1, k1);
 297                 VecScaleSelf(f2, k1);
 298                 VecScaleSelf(f3, k1);
 299         }
 300         /*      if (dihedral_uniq != NULL)
 301          k0 *= -dihedral_uniq[i]; */
 302         *energy = k0;
 303         *k1out = k1;
 304         *force1 = f1;
 305         VecDec(f1, f2);
 306         VecScale(*force2, f1, -1);
 307         VecDec(f2, f3);
 308         VecScale(*force3, f2, -1);
 309         VecScale(*force4, f3, -1);
 310         return 0;
 311 }
 312
 313 static void
 314 s_calc_dihedral_force_sub(MDArena *arena, Atom *ap, Int ndihedrals, Int *dihedrals, Int *dihedral_par_i, TorsionPar *dihedral_pars, Int *dihedral_uniq, Double *energies, Vector *forces)
 315 {
 316         int i;
 317         if (arena->mol->nsyms == 0)
 318                 dihedral_uniq = NULL;  /*  Ignore the symmetry info  */
 319         for (i = 0; i < ndihedrals; i++, dihedrals += 4, dihedral_par_i++) {
 320                 Atom *ap1, *ap2, *ap3, *ap4;
 321                 Double phi, en, k1;
 322                 Vector f1, f2, f3, f4;
 323                 TorsionPar *tp;
 324
 325                 if (*dihedral_par_i < 0)
 326                         continue;  /*  Ignore this entry  */
 327         /*      if (dihedral_uniq != NULL && dihedral_uniq[i] >= 0)
 328                         continue;  *//*  Non-unique dihedral  */
 329                 ap1 = ATOM_AT_INDEX(ap, dihedrals[0]);
 330                 ap2 = ATOM_AT_INDEX(ap, dihedrals[1]);
 331                 ap3 = ATOM_AT_INDEX(ap, dihedrals[2]);
 332                 ap4 = ATOM_AT_INDEX(ap, dihedrals[3]);
 333                 if (ap1->mm_exclude || ap2->mm_exclude || ap3->mm_exclude || ap4->mm_exclude)
 334                         continue;  /*  Skip non-occupied atoms  */
 335                 if (arena->nalchem_flags > 0) {
 336                         if (dihedrals[0] < arena->nalchem_flags && dihedrals[1] < arena->nalchem_flags
 337                                 && dihedrals[2] < arena->nalchem_flags && dihedrals[3] < arena->nalchem_flags
 338                                 && (arena->alchem_flags[dihedrals[0]] | arena->alchem_flags[dihedrals[1]]
 339                                         | arena->alchem_flags[dihedrals[2]] | arena->alchem_flags[dihedrals[3]]) == 3)
 340                                 continue;  /*  Interaction between vanishing and appearing groups is ignored  */
 341                 }
 342                 tp = dihedral_pars + *dihedral_par_i;
 343                 if (s_calc_dihedral_force_one(&(ap1->r), &(ap2->r), &(ap3->r), &(ap4->r), tp, &phi, &en, &k1, &f1, &f2, &f3, &f4) != 0)
 344                         continue;
 345
 346                 /*
 347                 VecSub(r21, ap[dihedrals[0]].r, ap[dihedrals[1]].r);
 348                 VecSub(r32, ap[dihedrals[1]].r, ap[dihedrals[2]].r);
 349                 VecSub(r43, ap[dihedrals[2]].r, ap[dihedrals[3]].r);
 350                 VecCross(v1, r21, r32);
 351                 VecCross(v2, r32, r43);
 352                 VecCross(v3, r32, v1);
 353                 w1 = VecLength(v1);
 354                 w2 = VecLength(v2);
 355                 w3 = VecLength(v3);
 356                 if (w1 < 1e-5 || w2 < 1e-5 || w3 < 1e-5)
 357                         continue;  //  The dihedral cannot be defined
 358                 w1 = 1.0/w1;
 359                 w2 = 1.0/w2;
 360                 w3 = 1.0/w3;
 361                 VecScaleSelf(v1, w1);
 362                 VecScaleSelf(v2, w2);
 363                 VecScaleSelf(v3, w3);
 364
 365                 //  The dihedral angle value
 366                 cosphi = VecDot(v1, v2);
 367                 sinphi = VecDot(v3, v2);
 368                 phi = -atan2(sinphi, cosphi);
 369
 370                 //  Repeat for multiple dihedral terms
 371                 k0 = k1 = 0.0;
 372                 for (n = 0; n < tp->mult; n++) {
 373                         Double k = tp->k[n];
 374                         Double phi0 = tp->phi0[n];
 375                         Int period = tp->period[n];
 376                         if (period > 0) {
 377                                 k0 += k * (1 + cos(period * phi + phi0));
 378                                 k1 -= period * k * sin(period * phi + phi0);
 379                         } else {
 380                                 //  A simple quadratic term
 381                                 phi -= phi0;
 382                                 if (phi < -PI)
 383                                         phi += 2 * PI;
 384                                 else if (phi > PI)
 385                                         phi -= 2 * PI;
 386                                 k0 += k * phi * phi;
 387                                 k1 += 2 * k * phi;
 388                         }
 389                 }
 390
 391                 if (sinphi < -0.1 || sinphi > 0.1) {
 392                         //  The sin form
 393                         Vector v4, v5, vw0;
 394                         VecScale(v4, v1, cosphi);
 395                         VecDec(v4, v2);
 396                         VecScaleSelf(v4, w1);
 397                         VecScale(v5, v2, cosphi);
 398                         VecDec(v5, v1);
 399                         VecScaleSelf(v5, w2);
 400                         k1 /= sinphi;
 401                         VecCross(f1, r32, v4);
 402                         VecCross(f3, v5, r32);
 403                         VecCross(f2, v4, r21);
 404                         VecCross(vw0, r43, v5);
 405                         VecInc(f2, vw0);
 406                         VecScaleSelf(f1, k1);
 407                         VecScaleSelf(f2, k1);
 408                         VecScaleSelf(f3, k1);
 409                 } else {
 410                         //  The cos form
 411                         Vector v6, v7, vw1, vw2;
 412                         VecScale(v6, v3, sinphi);
 413                         VecDec(v6, v2);
 414                         VecScaleSelf(v6, w3);
 415                         VecScale(v7, v2, sinphi);
 416                         VecDec(v7, v3);
 417                         VecScaleSelf(v7, w2);
 418                         k1 = -k1 / cosphi;
 419                         VecCross(vw1, v6, r32);
 420                         VecCross(f1, r32, vw1);
 421                         VecCross(f3, v7, r32);
 422                         VecCross(f2, r43, v7);
 423                         VecCross(vw1, r32, r21);
 424                         VecCross(vw2, v6, vw1);
 425                         VecInc(f2, vw2);
 426                         VecCross(vw1, r32, v6);
 427                         VecCross(vw2, r21, vw1);
 428                         VecInc(f2, vw2);
 429                         VecScaleSelf(f1, k1);
 430                         VecScaleSelf(f2, k1);
 431                         VecScaleSelf(f3, k1);
 432                 }
 433                 *energies += k0;
 434                 VecInc(forces[dihedrals[0]], f1);
 435                 VecDec(f1, f2);
 436                 VecDec(forces[dihedrals[1]], f1);
 437                 VecDec(f2, f3);
 438                 VecDec(forces[dihedrals[2]], f2);
 439                 VecDec(forces[dihedrals[3]], f3);
 440 */
 441                 *energies += en;
 442                 VecInc(forces[dihedrals[0]], f1);
 443                 VecInc(forces[dihedrals[1]], f2);
 444                 VecInc(forces[dihedrals[2]], f3);
 445                 VecInc(forces[dihedrals[3]], f4);
 446                 if (arena->debug_result && arena->debug_output_level > 1) {
 447                         fprintf(arena->debug_result, "dihedral(improper) force %d-%d-%d-%d: phi=%f, k1=%f, {%f %f %f}, {%f %f %f}, {%f %f %f}\n", dihedrals[0]+1, dihedrals[1]+1, dihedrals[2]+1, dihedrals[3]+1, phi*180/PI, k1 * INTERNAL2KCAL, f1.x * INTERNAL2KCAL, f1.y * INTERNAL2KCAL, f1.z * INTERNAL2KCAL, f2.x * INTERNAL2KCAL, f2.y * INTERNAL2KCAL, f2.z * INTERNAL2KCAL, f3.x * INTERNAL2KCAL, f3.y * INTERNAL2KCAL, f3.z * INTERNAL2KCAL);
 448                 }
 449         }
 450 }
 451
 452 static void
 453 s_calc_dihedral_force(MDArena *arena)
 454 {
 455         Molecule *mol = arena->mol;
 456         Parameter *par = arena->par;
 457         Double *energies = &arena->energies[kDihedralIndex];
 458         Vector *forces = &arena->forces[kDihedralIndex * arena->mol->natoms];
 459         s_calc_dihedral_force_sub(arena, mol->atoms, mol->ndihedrals, mol->dihedrals, arena->dihedral_par_i, par->dihedralPars, NULL/*arena->sym_dihedral_uniq*/, energies, forces);
 460 }
 461
 462 static void
 463 s_calc_improper_force(MDArena *arena)
 464 {
 465         Molecule *mol = arena->mol;
 466         Parameter *par = arena->par;
 467         Double *energies = &arena->energies[kImproperIndex];
 468         Vector *forces = &arena->forces[kImproperIndex * arena->mol->natoms];
 469         s_calc_dihedral_force_sub(arena, mol->atoms, mol->nimpropers, mol->impropers, arena->improper_par_i, par->improperPars, NULL/*arena->sym_improper_uniq*/, energies, forces);
 470 }
 471
 472 /*  ==================================================================== */
 473 /*  md_check_verlet_list: only for debugging the verlet list generation  */
 474 /*  ==================================================================== */
 475
 476 typedef struct md_check_record {
 477         Int i, j, dx, dy, n;
 478         Double len;
 479 } md_check_record;
 480
 481 static int
 482 s_md_check_verlet_comparator(const void *ap, const void *bp)
 483 {
 484         Double d = ((const md_check_record *)ap)->len - ((const md_check_record *)bp)->len;
 485         return (d < 0 ? -1 : (d > 0 ? 1 : 0));
 486 }
 487
 488 void
 489 md_check_verlet_list(MDArena *arena)
 490 {
 491         int i, j, k;
 492         int dx, dy, nn;
 493         int ndx, ndy, ndz;
 494         Atom *api, *apj;
 495         Vector cell_offsets[27];
 496         XtalCell *cell = arena->mol->cell;
 497         md_check_record *cr;
 498         int ncr;
 499
 500         ndx = (arena->periodic_a ? 1 : 0);
 501         ndy = (arena->periodic_b ? 1 : 0);
 502         ndz = (arena->periodic_c ? 1 : 0);
 503         if (cell != NULL && (ndx != 0 || ndy != 0 || ndz != 0)) {
 504                 nn = 0;
 505                 for (i = -1; i <= 1; i++) {
 506                         for (j = -1; j <= 1; j++) {
 507                                 for (k = -1; k <= 1; k++) {
 508                                         VecZero(cell_offsets[nn]);
 509                                         VecScaleInc(cell_offsets[nn], cell->axes[0], i);
 510                                         VecScaleInc(cell_offsets[nn], cell->axes[1], j);
 511                                         VecScaleInc(cell_offsets[nn], cell->axes[2], k);
 512                                         nn++;
 513                                 }
 514                         }
 515                 }
 516         }
 517         /*  Debugging is done with x-y 2-dimensional system  */
 518         ncr = arena->mol->natoms * (arena->mol->natoms - 1) / 2 * 9;
 519         cr = (md_check_record *)malloc(sizeof(md_check_record) * ncr);
 520         k = 0;
 521         for (i = 0, api = arena->mol->atoms; i < arena->mol->natoms; i++, api = ATOM_NEXT(api)) {
 522                 for (j = i + 1, apj = ATOM_AT_INDEX(arena->mol->atoms, j); j < arena->mol->natoms; j++, apj = ATOM_NEXT(apj)) {
 523                         Vector dr;
 524                         VecSub(dr, apj->r, api->r);
 525                         for (dx = -1; dx <= 1; dx++) {
 526                                 for (dy = -1; dy <= 1; dy++) {
 527                                         Vector dr2 = dr;
 528                                         nn = dx * 9 + dy * 3 + 13;
 529                                         VecInc(dr2, cell_offsets[nn]);
 530                                         cr[k].i = i;
 531                                         cr[k].j = j;
 532                                         cr[k].dx = dx;
 533                                         cr[k].dy = dy;
 534                                         cr[k].n = -1;
 535                                         cr[k].len = VecLength(dr2);
 536                                         k++;
 537                                 }
 538                         }
 539                 }
 540         }
 541         for (k = 0; k < arena->nverlets; k++) {
 542                 i = arena->verlets[k].n1;
 543                 j = arena->verlets[k].n2;
 544                 if (i > j) {
 545                         j = i;
 546                         i = arena->verlets[k].n2;
 547                 }
 548                 dx = arena->verlets[k].symop.dx;
 549                 dy = arena->verlets[k].symop.dy;
 550                 nn = (i * arena->mol->natoms - i * (i + 1) / 2 + (j - i) - 1) * 9 + (dx + 1) * 3 + dy + 1;
 551                 if (cr[nn].i != i || cr[nn].j != j || cr[nn].dx != dx || cr[nn].dy != dy)
 552                         fprintf(arena->log_result, "Debug: internal inconsistency\n");
 553                 cr[nn].n = k;
 554         }
 555
 556         mergesort(cr, ncr, sizeof(md_check_record), s_md_check_verlet_comparator);
 557         for (i = 0; i < ncr; i++) {
 558                 Double len2;
 559                 len2 = -1;
 560                 if (cr[i].n >= 0)
 561                         len2 = arena->verlets[cr[i].n].length;
 562                 fprintf(arena->log_result, "Debug: %5d (i=%4d,j=%4d) [dx=%4d, dy=%4d] n=%4d %f %f\n", i, cr[i].i, cr[i].j, cr[i].dx, cr[i].dy, cr[i].n, cr[i].len, len2);
 563                 if (cr[i].len > arena->pairlist_distance)
 564                         break;
 565         }
 566         while (i < ncr) {
 567                 if (cr[i].n != -1) {
 568                         fprintf(arena->log_result, "Debug: %5d (i=%4d,j=%4d) [dx=%4d, dy=%4d] n=%4d %f %f\n", i, cr[i].i, cr[i].j, cr[i].dx, cr[i].dy, cr[i].n, cr[i].len, arena->verlets[cr[i].n].length);
 569                 }
 570                 i++;
 571         }
 572         fflush(arena->log_result);
 573         free(cr);
 574 }
 575
 576 /*  ==================================================================== */
 577
 578 /*  Count all lattice points within distance d from point v  */
 579 /*  Ref. "Algorithms for the Shortest and Closest Lattice Vector Problems"
 580     Guillaume Hanrot, Xavier Pujol, Damien Stehle
 581     Coding and Cryptology; Lecture Notes in Computer Science, 2011 vol. 6639/2011 pp. 159-190
 582     DOI: 10.1007/978-3-642-20901-7_10  */
 583 static int
 584 s_enum_neighbors(MDArena *arena, Vector v, Double d)
 585 {
 586         Vector bn[3];  /*  Base vectors for the cell axes  */
 587         Double bl[3];  /*  Square of the length of base vectors  */
 588         Double mu[3];  /*  Non-diagonal terms for Gram-Schmidt orthogonalization  */
 589         XtalCell *cell = arena->mol->cell;
 590         Int dim;       /*  Number of periodic axes  */
 591         Int ijk[3];    /*  Renumber the indices. For 1-dimensional case, axes[ijk[0]] is the periodic axis.
 592                                            For 2-dimensional case, axes[ijk[0]] and axes[ijk[1]] are the periodic axes.
 593                                            Only uses exchange of two elements, so that ijk[ijk[i]] == i for i=0,1,2  */
 594         Vector axes[3]; /*  Renumbered cell axes  */
 595         Double t[3];   /*  Projected vector v on bn[]  */
 596         Int count;     /*  Number of results  */
 597         Int x[3];      /*  Lattice point candidate  */
 598         Double r[3];
 599         Int i;
 600         Double w;
 601
 602         dim = (arena->periodic_a != 0) + (arena->periodic_b != 0) + (arena->periodic_c != 0);
 603         ijk[0] = 0; ijk[1] = 1; ijk[2] = 2;
 604         if (dim == 0) {
 605                 /*  Non-periodic case: check whether (0,0,0) is within d or not  */
 606                 if (VecLength(v) < d) {
 607                         x[0] = x[1] = x[2] = 0;
 608                         AssignArray(&arena->lattice_offsets, &arena->nlattice_offsets, sizeof(Int) * 3, 0, x);
 609                         return 1;
 610                 } else return 0;
 611         }
 612
 613         if (dim == 1) {
 614                 if (arena->periodic_b) {
 615                         ijk[0] = 1; ijk[1] = 0;
 616                 } else if (arena->periodic_c) {
 617                         ijk[0] = 2; ijk[2] = 0;
 618                 }
 619         } else if (dim == 2) {
 620                 if (!arena->periodic_a) {
 621                         ijk[2] = 0; ijk[0] = 2;
 622                 } else if (!arena->periodic_b) {
 623                         ijk[1] = 0; ijk[0] = 1;
 624                 }
 625         }
 626         axes[0] = cell->axes[ijk[0]];
 627         axes[1] = cell->axes[ijk[1]];
 628         axes[2] = cell->axes[ijk[2]];
 629
 630         /*  Gram-Schmidt orthogonalization  */
 631         bn[0] = axes[0];
 632         bl[0] = VecLength2(bn[0]);
 633         mu[0] = VecDot(axes[1], bn[0]) / bl[0];  /*  bl[0] should be non-zero  */
 634         if (dim == 1) {
 635                 VecZero(bn[1]);
 636                 VecZero(bn[2]);
 637                 bl[1] = bl[2] = 0.0;
 638                 mu[1] = mu[2] = 0.0;
 639         } else {
 640                 bn[1] = axes[1];
 641                 VecScaleInc(bn[1], bn[0], -mu[0]);
 642                 bl[1] = VecLength2(bn[1]);
 643                 mu[1] = VecDot(axes[2], bn[0]) / bl[0];
 644                 if (dim == 2 || bl[1] < 1e-10) {
 645                         VecZero(bn[2]);
 646                         bl[2] = mu[2] = 0.0;
 647                         if (bl[1] < 1e-10)
 648                                 dim = 1;
 649                         else dim = 2;
 650                 } else {
 651                         mu[2] = VecDot(axes[2], bn[1]) / bl[1];
 652                         bn[2] = axes[2];
 653                         VecScaleInc(bn[2], bn[0], -mu[1]);
 654                         VecScaleInc(bn[2], bn[1], -mu[2]);
 655                         bl[2] = VecLength2(bn[2]);
 656                         if (bl[2] < 1e-10)
 657                                 dim = 2;
 658                         else dim = 3;
 659                 }
 660         }
 661
 662         /*  Project the target vector  */
 663         t[0] = t[1] = t[2] = 0.0;
 664         t[0] = VecDot(v, bn[0]) / bl[0];
 665         if (dim >= 2) {
 666                 t[1] = VecDot(v, bn[1]) / bl[1];
 667                 if (dim >= 3) {
 668                         t[2] = VecDot(v, bn[2]) / bl[2];
 669                 }
 670         }
 671
 672         /*  Enumerate  */
 673         count = 0;
 674         x[0] = x[1] = x[2] = 0;
 675         r[0] = r[1] = r[2] = 0.0;
 676         i = dim - 1;
 677         w = -10000.0;
 678         x[i] = ceil(t[i] - d / sqrt(bl[i]));
 679         while (1) {
 680                 Int j;
 681                 Double w2;
 682                 if (w == -10000.0) {
 683                         w = 0.0;
 684                         for (j = i + 1; j < dim; j++) {
 685                                 w += x[j] * mu[j + i - 1];
 686                         }
 687                 }
 688                 w2 = x[i] - t[i] + w;
 689                 r[i] = w2 * w2 * bl[i];
 690                 w2 = 0.0;
 691                 for (j = i; j < dim; j++) {
 692                         w2 += r[j];
 693                 }
 694                 w2 = d * d - w2;
 695                 if (w2 >= 0.0) {
 696                         if (i == 0) {
 697                                 /*  Found  */
 698                                 Int xx[3];
 699                                 xx[0] = x[ijk[0]];
 700                                 xx[1] = x[ijk[1]];
 701                                 xx[2] = x[ijk[2]];
 702                                 AssignArray(&arena->lattice_offsets, &arena->nlattice_offsets, sizeof(Int) * 3, count, xx);
 703                                 count++;
 704                                 x[0]++;
 705                                 w = -10000.0;
 706                         } else {
 707                                 /*  Step down  */
 708                                 i--;
 709                                 w = 0.0;
 710                                 for (j = i + 1; j < dim; j++) {
 711                                         w += x[j] * mu[j + i - 1];
 712                                 }
 713                                 x[i] = ceil(t[i] - w - sqrt(w2 / bl[i]));
 714                         }
 715                 } else {
 716                         i++;
 717                         if (i >= dim)
 718                                 break;
 719                         x[i]++;
 720                         w = -10000.0;
 721                 }
 722         }
 723         return count;
 724 }
 725
 726 /*  Update the Verlet list (for non-bonding interaction)  */
 727 static void
 728 s_make_verlet_list(MDArena *arena)
 729 {
 730         int i, j, k, n, nn, natoms;
 731         int dx, dy, dz, ndx, ndy, ndz;
 732         Molecule *mol = arena->mol;
 733         Atom *atoms = mol->atoms;
 734         Atom *api, *apj;
 735         MDVerlet *vl = arena->verlets;
 736         Vector *vdr = arena->verlets_dr;
 737         Double limit, limit2;
 738         Byte use_sym;
 739         XtalCell *cell = mol->cell;
 740         Parameter *par = arena->par;
 741         Double vdw_cutoff;
 742
 743         natoms = mol->natoms;
 744         for (i = 0; i < natoms; i++)
 745                 VecZero(vdr[i]);
 746
 747         ndx = (arena->periodic_a ? 1 : 0);
 748         ndy = (arena->periodic_b ? 1 : 0);
 749         ndz = (arena->periodic_c ? 1 : 0);
 750         limit = arena->pairlist_distance;
 751         limit2 = limit * limit;
 752         n = 0;
 753         use_sym = (arena->mol->nsyms > 0);
 754         for (i = 0, api = atoms; i < natoms; i++, api++) {
 755                 MDExclusion *exinfo = arena->exinfo + i;
 756                 Int index4 = (exinfo + 1)->index0;
 757                 Int vdw_idx1 = arena->vdw_par_i[i];
 758                 Int vdw_idx, vdw_idx2;
 759                 arena->verlet_i[i] = n;
 760
 761                 if (api->anchor != NULL)
 762                         continue;  /*  Skip pi_anchors  */
 763
 764                 for (j = i, apj = atoms + j; j < natoms; j++, apj++) {
 765                         Vector rij;
 766                         Double lenij2;
 767                         Int *ip, index0;
 768                         int exflag = 0;
 769                         int mult = 1;
 770                 /*      int dxbase, dybase, dzbase; */
 771                         int count;
 772
 773                         /*  Pi anchors  */
 774                         if (apj->anchor != NULL)
 775                                 continue;
 776
 777                         /*  Fixed atoms  */
 778                         if (api->fix_force < 0 && apj->fix_force < 0)
 779                                 continue;
 780
 781                         /*  Non-occupied atoms  */
 782                         if (api->mm_exclude || apj->mm_exclude)
 783                                 continue;
 784
 785                         /*  If no symmetry expansion, then no interaction with self  */
 786                         if (ndx + ndy + ndz == 0 && i == j)
 787                                 continue;
 788
 789                         VecSub(rij, apj->r, api->r);
 790
 791                         /*  Calculate the cell offset for the nearest neighbor  */
 792                         if (ndx + ndy + ndz != 0 && apj->periodic_exclude == 0 && cell != NULL) {
 793                                 count = s_enum_neighbors(arena, rij, limit);
 794                         } else {
 795                                 static Int sZeros[3] = {0, 0, 0};
 796                                 AssignArray(&arena->lattice_offsets, &arena->nlattice_offsets, sizeof(Int) * 3, 0, sZeros);
 797                                 count = 1;
 798                         }
 799                         /*  Check the specific cutoff table for (i, j) pair  */
 800                         vdw_idx2 = arena->vdw_par_i[j];
 801                         if (vdw_idx1 < 0 || vdw_idx2 < 0)
 802                                 vdw_idx = par->nvdwPars * par->nvdwPars;  /*  A null record  */
 803                         else if (vdw_idx1 < vdw_idx2)
 804                                 vdw_idx = vdw_idx1 * par->nvdwPars + vdw_idx2;
 805                         else
 806                                 vdw_idx = vdw_idx2 * par->nvdwPars + vdw_idx1;
 807                         vdw_cutoff = arena->cutoff;
 808                         for (k = par->nvdwCutoffPars - 1; k >= 0; k--) {
 809                                 VdwCutoffPar *cr = par->vdwCutoffPars + k;
 810                                 if (((cr->type1 == kAtomTypeWildcard || cr->type1 == api->type) && (cr->type2 == kAtomTypeWildcard || cr->type2 == apj->type)) ||
 811                                         ((cr->type1 == kAtomTypeWildcard || cr->type1 == apj->type) && (cr->type2 == kAtomTypeWildcard || cr->type2 == api->type))) {
 812                                         vdw_cutoff = cr->cutoff;
 813                                         break;
 814                                 }
 815                                 if ((cr->type1 == i && cr->type2 == j) || (cr->type1 == j && cr->type2 == i)) {
 816                                         vdw_cutoff = cr->cutoff;
 817                                         break;
 818                                 }
 819                         }
 820                         if (vdw_cutoff < 0) {
 821                                 /*  A negative value of specific cutoff means "cutoff at r_eq * (-specific_cutoff)  */
 822                                 MDVdwCache *cp = &(arena->vdw_cache[vdw_idx]);
 823                                 Double r_eq = pow(cp->par.A / cp->par.B * 2.0, 1.0/6.0);
 824                                 vdw_cutoff = r_eq * (-vdw_cutoff);
 825                         }
 826
 827                         /*  Search for pairs  */
 828                         for (nn = 0; nn < count; nn++) {
 829                                 Vector rij0 = rij;
 830                                 dx = arena->lattice_offsets[nn * 3];
 831                                 dy = arena->lattice_offsets[nn * 3 + 1];
 832                                 dz = arena->lattice_offsets[nn * 3 + 2];
 833                                 if (dx == 0 && dy == 0 && dz == 0) {
 834                                         /*  Pair within the unit cell  */
 835                                         if (i == j)
 836                                                 continue;  /*  No interaction with self  */
 837                                         /*  Is this pair to be excluded?  */
 838                                         for (index0 = exinfo->index0, ip = arena->exlist + index0; index0 < index4; index0++, ip++) {
 839                                                 if (*ip == j)
 840                                                         break;
 841                                         }
 842                                         if (index0 < exinfo->index3)
 843                                                 continue;  /*  Special exclusion, 1-2, 1-3  */
 844                                         if (index0 < index4)
 845                                                 exflag = 1;  /*  1-4 interaction  */
 846                                 } else if (apj->periodic_exclude) {
 847                                         continue;
 848                                 } else {
 849                                         VecScaleInc(rij0, cell->axes[0], dx);
 850                                         VecScaleInc(rij0, cell->axes[1], dy);
 851                                         VecScaleInc(rij0, cell->axes[2], dz);
 852                                 /*      VecInc(rij0, cell_offsets[nn]); */
 853                                         exflag = 0;
 854                                 }
 855
 856                                 lenij2 = VecLength2(rij0);
 857                                 if (lenij2 <= limit2) {
 858                                         MDVerlet *vlp;
 859                                         if (n >= arena->max_nverlets) {
 860                                                 arena->max_nverlets += 32;
 861                                                 vl = (MDVerlet *)realloc(vl, sizeof(MDVerlet) * arena->max_nverlets);
 862                                                 if (vl == NULL)
 863                                                         md_panic(arena, "Low memory");
 864                                         }
 865                                         vlp = &vl[n];
 866                                         vlp->vdw_type = (exflag ? 1 : 0);
 867                                         vlp->mult = mult;
 868                                         vlp->symop.dx = dx;
 869                                         vlp->symop.dy = dy;
 870                                         vlp->symop.dz = dz;
 871                                         vlp->symop.sym = 0;
 872                                         vlp->symop.alive = (vlp->symop.dx != 0 || vlp->symop.dy != 0 || vlp->symop.dz != 0);
 873                                         vlp->index = vdw_idx;
 874                                         vlp->n1 = i;
 875                                         vlp->n2 = j;
 876                                         vlp->vdw_cutoff = vdw_cutoff;
 877                                         vlp->length = sqrt(lenij2);
 878                                         n++;
 879                                 } /* end if lenij2 <= limit */
 880                         } /* end loop nn */
 881                 } /* end loop j */
 882
 883         } /* end loop i */
 884         arena->verlet_i[natoms] = n;
 885         arena->nverlets = n;
 886         arena->verlets = vl;
 887         arena->last_verlet_step = arena->step;
 888
 889         if (arena->debug_result && arena->debug_output_level > 2) {
 890                 fprintf(arena->debug_result, "\n  Verlet list at step %d\n", arena->step);
 891                 fprintf(arena->debug_result, "  {atom1 atom2 vdw_type (=1 if 1-4 bonded) vdw_index cell_translation}\n");
 892                 for (i = 0; i < arena->nverlets; i++) {
 893                         fprintf(arena->debug_result, "{%d %d %d %d %06d}%c", vl[i].n1+1, vl[i].n2+1, vl[i].vdw_type, vl[i].index, (vl[i].symop.dx + 50) * 10000 + (vl[i].symop.dy + 50) * 100 + (vl[i].symop.dz + 50), (i % 4 == 3 ? '\n' : ' '));
 894                 }
 895                 if (i % 4 != 0)
 896                         fprintf(arena->debug_result, "\n");
 897         }
 898 }
 899
 900 /*  Calculate the nonbonded force  */
 901 /*  group_flags[] is used for calculation of pair-specific interaction between
 902     two groups of atoms  */
 903 static void
 904 s_calc_nonbonded_force_sub(MDArena *arena, Double *energies, Double *eenergies, Vector *forces, Vector *eforces, const MDGroupFlags *group_flags_1, const MDGroupFlags *group_flags_2)
 905 {
 906         int i;
 907         Vector *vdr;
 908         Double limit, elimit, slimit, vlimit, dielec_r;
 909         Double slimit6, slimit12, vlimit6, vlimit12;
 910         Double lambda, dlambda;
 911 /*      Double beta2;  */
 912         MDVerlet *vl;
 913 /*      MDVdwCache *vdw_cache = arena->vdw_cache; */
 914         Atom *atoms = arena->mol->atoms;
 915         XtalCell *cell = arena->mol->cell;
 916         Vector *eforces_corr;
 917         Double *eenergies_corr;
 918         Int do_ewald;
 919
 920         if (group_flags_1 == NULL && arena->use_ewald) {
 921                 eforces_corr = &arena->forces[kESCorrectionIndex * arena->mol->natoms];
 922                 eenergies_corr = &arena->energies[kESCorrectionIndex];
 923                 do_ewald = 1;
 924         } else {
 925                 do_ewald = 0;
 926                 eforces_corr = NULL;
 927                 eenergies_corr = NULL;
 928         }
 929
 930 #if MINIMIZE_CELL
 931         memset(arena->cell_forces, 0, sizeof(Double) * 12);
 932 #endif
 933
 934         /*  Check whether the Verlet list needs update  */
 935         if (arena->last_verlet_step >= 0) {
 936                 i = arena->mol->natoms - 1;
 937                 vdr = arena->verlets_dr + i;
 938                 if (arena->cutoff > arena->electro_cutoff)
 939                         limit = (arena->pairlist_distance - arena->cutoff) * 0.5;
 940                 else
 941                         limit = (arena->pairlist_distance - arena->electro_cutoff) * 0.5;
 942                 limit = limit * limit;
 943                 for ( ; i >= 0; i--, vdr--) {
 944                         if (VecLength2(*vdr) >= limit)
 945                                 break;
 946                 }
 947         } else i = 0;
 948         if (i >= 0)
 949                 s_make_verlet_list(arena);
 950
 951         /*  Calculate the non-bonded interaction for each pair in the Verlet list  */
 952         vlimit = arena->cutoff * arena->cutoff;
 953         elimit = arena->electro_cutoff * arena->electro_cutoff;
 954         slimit = arena->switch_distance * arena->switch_distance;
 955         if (slimit > vlimit)
 956                 slimit = 0.0;  /*  No switching  */
 957         if (slimit > 0.0) {
 958                 slimit6 = 1/(slimit * slimit * slimit);
 959                 slimit12 = slimit6 * slimit6;
 960         }
 961         vlimit6 = 1/(vlimit * vlimit * vlimit);
 962         vlimit12 = vlimit6 * vlimit6;
 963         dielec_r = COULOMBIC / arena->dielectric;
 964 /*      if (arena->use_ewald)
 965                 beta2 = arena->ewald_beta * arena->ewald_beta;
 966         else beta2 = 0.0; */
 967         for (i = arena->nverlets - 1, vl = arena->verlets + i; i >= 0; i--, vl--) {
 968                 Double A, B, k0, k1;
 969                 MDVdwCache *vp = arena->vdw_cache + vl->index;
 970                 Atom *ap1, *ap2;
 971                 Vector rij, fij, rij2, fij2;
 972                 Double r2, w2, w6, w12;
 973                 if (group_flags_1 != NULL && group_flags_2 != NULL) {
 974                         if (!(get_group_flag(group_flags_1, vl->n1) && get_group_flag(group_flags_2, vl->n2))
 975                         && !(get_group_flag(group_flags_1, vl->n2) && get_group_flag(group_flags_2, vl->n1)))
 976                                 continue;
 977                 }
 978
 979                 if (arena->nalchem_flags > 0) {
 980                         char c1, c2;
 981                         if (vl->n1 < arena->nalchem_flags)
 982                                 c1 = arena->alchem_flags[vl->n1];
 983                         else c1 = 0;
 984                         if (vl->n2 < arena->nalchem_flags)
 985                                 c2 = arena->alchem_flags[vl->n2];
 986                         else c2 = 0;
 987                         if ((c1 == 1 && c2 == 2) || (c1 == 2 && c2 == 1))
 988                                 continue;
 989                         if (c1 == 1 || c2 == 1) {
 990                                 lambda = (1.0 - arena->alchem_lambda);
 991                                 dlambda = -arena->alchem_dlambda;
 992                         } else if (c1 == 2 || c2 == 2) {
 993                                 lambda = arena->alchem_lambda;
 994                                 dlambda = arena->alchem_dlambda;
 995                         } else {
 996                                 lambda = 1.0;
 997                                 dlambda = 0.0;
 998                         }
 999                 } else {
1000                         lambda = 1.0;
1001                         dlambda = 0.0;
1002                 }
1003
1004                 if (vl->vdw_type == 1) {
1005                         A = vp->par.A14;
1006                         B = vp->par.B14;
1007                 /*      vofs = vp->vcut14; */
1008                 } else {
1009                         A = vp->par.A;
1010                         B = vp->par.B;
1011                 /*      vofs = vp->vcut; */
1012                 }
1013                 ap1 = &atoms[vl->n1];
1014                 ap2 = &atoms[vl->n2];
1015                 rij = ap2->r;
1016                 if (vl->symop.alive) {
1017                         VecScaleInc(rij, cell->axes[0], vl->symop.dx);
1018                         VecScaleInc(rij, cell->axes[1], vl->symop.dy);
1019                         VecScaleInc(rij, cell->axes[2], vl->symop.dz);
1020                 }
1021                 VecDec(rij, ap1->r);
1022                 r2 = VecLength2(rij);
1023                 if (r2 >= elimit)
1024                         continue;
1025                 limit = vl->vdw_cutoff * vl->vdw_cutoff;
1026                 if (r2 >= limit)
1027                         continue;
1028
1029                 fij.x = fij.y = fij.z = 0.0;
1030                 fij2 = fij;
1031
1032                 /*  Coulombic force  */
1033                 w12 = ap1->charge * ap2->charge * dielec_r;
1034                 if (w12 != 0.0) {
1035                         Double k2, k3, er;
1036                         if (do_ewald || !arena->use_xplor_shift) {
1037                                 w2 = 1.0 / sqrt(r2);
1038                                 w6 = w2 / r2;
1039                                 k0 = w12 * w2;
1040                                 k1 = w12 * w6;
1041                                 if (do_ewald) {
1042                                         k2 = arena->ewald_beta / w2;
1043                                         er = erf(k2);
1044                                         k3 = -k1 * (-er + 2 * k2 * exp(-k2 * k2) * PI2R);
1045                                         k2 = -k0 * er;
1046                                 } else {
1047                                 /*      vofs = -w12 / arena->electro_cutoff;
1048                                         k0 += vofs; */  /*  Discontinue potential at cutoff  */
1049                                 }
1050                         } else {
1051                                 w2 = 1.0 / sqrt(r2);
1052                                 k0 = r2 / elimit - 1.0;
1053                                 k0 = k0 * k0;
1054                                 k0 = w12 * k0 * w2;
1055                                 k1 = (3.0 * r2 / elimit - 2.0) / elimit - 1.0 / r2;
1056                                 k1 = -w12 * k1 * w2;
1057                         }
1058                         if (vl->vdw_type == 1) {
1059                                 k0 *= arena->scale14_elect;
1060                                 k1 *= arena->scale14_elect;
1061                         }
1062                         if (do_ewald) {
1063                                 VecScale(fij, rij, k3);
1064                                 *eenergies_corr += k2 / vl->mult;
1065                                 VecDec(eforces_corr[vl->n1], fij);
1066                                 VecInc(eforces_corr[vl->n2], fij);
1067                                 if (arena->debug_result && arena->debug_output_level > 1) {
1068                                         fprintf(arena->debug_result, "Electrostatic correction force %d-%d: r=%f, k0=%f, k1=%f, {%f %f %f}\n", vl->n1+1, vl->n2+1, sqrt(r2), k2/KCAL2INTERNAL, k3*sqrt(r2)/KCAL2INTERNAL, fij.x/KCAL2INTERNAL, fij.y/KCAL2INTERNAL, fij.z/KCAL2INTERNAL);
1069                                 }
1070                         }
1071                         VecScale(fij, rij, k1);
1072                         *eenergies += k0 / vl->mult;
1073                         if (eforces != NULL) {
1074                                 VecDec(eforces[vl->n1], fij);
1075                                 VecInc(eforces[vl->n2], fij);
1076                         }
1077                         fij2 = fij;
1078                         if (arena->debug_result && arena->debug_output_level > 1) {
1079                                 fprintf(arena->debug_result, "nonbonded(electrostatic) force %d-%d: r=%f, k0=%f, k1=%f, {%f %f %f}\n", vl->n1+1, vl->n2+1, sqrt(r2), k0/KCAL2INTERNAL, k1*sqrt(r2)/KCAL2INTERNAL, fij.x/KCAL2INTERNAL, fij.y/KCAL2INTERNAL, fij.z/KCAL2INTERNAL);
1080                         }
1081                 }
1082
1083                 /*  van der Waals force  */
1084                 /*  w2 = 1/(r**2), w6 = 1/(r**6), w12 = 1/(r**12)  */
1085                 w2 = 1.0 / r2;
1086                 w6 = w2 * w2 * w2;
1087                 w12 = w6 * w6;
1088                 if (slimit == 0 || r2 < slimit) {
1089                         k0 = A * w12 - B * w6;
1090                         k1 = 6 * w2 * (2.0 * A * w12 - B * w6);
1091                 } else {
1092                         /*  Between switch_distance and cutoff: Linear switching function  */
1093                         w2 = 1.0 / slimit;
1094                         w6 = w2 * w2 * w2;
1095                         w12 = w6 * w6;
1096                         k1 = (A * w12 - B * w6) / (arena->cutoff - arena->switch_distance);
1097                         w2 = sqrt(r2);
1098                         k0 = k1 * (arena->cutoff - w2);
1099                         k1 = -k1 / w2;
1100                 }
1101                 if (vl->vdw_type == 1) {
1102                         k0 *= arena->scale14_vdw;
1103                         k1 *= arena->scale14_vdw;
1104                 }
1105                 k0 /= vl->mult;
1106                 k1 *= lambda;
1107                 VecScale(fij, rij, k1);
1108                 *energies += k0 * lambda;
1109                 if (forces != NULL) {
1110                         VecDec(forces[vl->n1], fij);
1111                         VecInc(forces[vl->n2], fij);
1112                 }
1113                 VecInc(fij2, fij);
1114                 if (dlambda != 0.0)
1115                         arena->alchem_energy += k0 * dlambda;
1116
1117                 if (arena->debug_result && arena->debug_output_level > 1) {
1118                         fprintf(arena->debug_result, "nonbonded(vdw) force %d-%d: r=%f, k0=%f, k1=%f, {%f %f %f}\n", vl->n1+1, vl->n2+1, sqrt(r2), k0/KCAL2INTERNAL, k1*sqrt(r2)/KCAL2INTERNAL, fij.x/KCAL2INTERNAL, fij.y/KCAL2INTERNAL, fij.z/KCAL2INTERNAL);
1119                 }
1120
1121 #if MINIMIZE_CELL
1122                 if (arena->is_minimizing && arena->minimize_cell && (ap2->symop.alive || vl->symop.alive)) {
1123                         /*  Force for changing cell parameters  */
1124                         Vector rr;
1125                         Int j, k;
1126                         Transform da, tr, symtr, symtr2;
1127                         if (ap2->symop.alive)
1128                                 rij2 = ATOM_AT_INDEX(arena->mol->atoms, ap2->symbase)->r;
1129                         else rij2 = ap2->r;
1130                         TransformVec(&rr, cell->rtr, &rij2);
1131                         memmove(symtr, SYMMETRY_AT_INDEX(arena->mol->syms, ap2->symop.sym), sizeof(Transform));
1132                         symtr[9] += ap2->symop.dx + vl->symop.dx;
1133                         symtr[10] += ap2->symop.dy + vl->symop.dy;
1134                         symtr[11] += ap2->symop.dz + vl->symop.dz;
1135                         TransformMul(symtr2, symtr, cell->rtr);
1136                         TransformMul(symtr2, cell->tr, symtr2);
1137                         for (j = 0; j < 12; j++) {
1138                                 Vector rrr;
1139                                 if (arena->periodic_a == 0 && (j == 0 || j == 1 || j == 2 || j == 3 || j == 6 || j == 9))
1140                                         continue;
1141                                 if (arena->periodic_b == 0 && (j == 1 || j == 3 || j == 4 || j == 5 || j == 7 || j == 10))
1142                                         continue;
1143                                 if (arena->periodic_c == 0 && (j == 2 || j == 5 || j == 6 || j == 7 || j == 8 || j == 11))
1144                                         continue;
1145                                 memset(da, 0, sizeof(Transform));
1146                                 da[j] = 1.0;
1147                                 TransformMul(tr, symtr2, da);
1148                                 /*  The translation part should be adjusted, because the derivative matrix
1149                                         has the (3,3) element as 0.0 (instead of 1.0). Thus,
1150                                     | R v | | R' v' | = | R*R' R*v' |   (instead of R*v'+v)
1151                                     | 0 1 | | 0  0  |   | 0    0    |
1152                                     da * symtr does not have this problem, because the derivative matrix is
1153                                     multipled from the left.  */
1154                                 tr[9] -= symtr2[9];
1155                                 tr[10] -= symtr2[10];
1156                                 tr[11] -= symtr2[11];
1157                                 TransformMul(da, da, symtr);
1158                                 for (k = 0; k < 12; k++)
1159                                         da[k] -= tr[k];
1160                                 TransformVec(&rrr, da, &rr);
1161                                 arena->cell_forces[j] += VecDot(fij2, rrr);
1162                         }
1163                 }
1164 #endif
1165
1166         }
1167
1168         if (arena->use_ewald != 0) {
1169                 /*  Calculate correction terms for excluded atom pairs  */
1170                 Atom *api, *apj;
1171                 Int j, k;
1172                 Vector rij;
1173                 Double d, dd, w0, w1, w12;
1174                 for (i = 0, api = arena->mol->atoms; i < arena->mol->natoms; i++, api = ATOM_NEXT(api)) {
1175                         MDExclusion *exinfo = arena->exinfo + i;
1176                         for (k = exinfo->index0; k < exinfo->index3; k++) {
1177                                 j = arena->exlist[k];
1178                                 if (j <= i)
1179                                         continue;
1180                                 apj = ATOM_AT_INDEX(arena->mol->atoms, j);
1181                                 VecSub(rij, api->r, apj->r);
1182                                 d = VecLength(rij);
1183                                 w12 = -api->charge * apj->charge * COULOMBIC / arena->dielectric * 0.5;
1184                                 dd = arena->ewald_beta * d;
1185                                 w0 = w12 * erf(dd) / d;
1186                                 w1 = (2 * w0 - w12 * 4 * dd * PI2R * exp(-dd * dd)) / (d * d);
1187                                 VecScaleSelf(rij, w1);
1188                                 if (arena->debug_result && arena->debug_output_level > 1) {
1189                                         fprintf(arena->debug_result, "Electrostatic correction force (excluded) %d-%d: r=%f, k0=%f, k1=%f, {%f %f %f}\n", i+1, j+1, d, w0/KCAL2INTERNAL, w0*d/KCAL2INTERNAL, rij.x/KCAL2INTERNAL, rij.y/KCAL2INTERNAL, rij.z/KCAL2INTERNAL);
1190                                 }
1191                                 VecInc(eforces_corr[i], rij);
1192                                 VecDec(eforces_corr[j], rij);
1193                                 *eenergies_corr += w0;
1194                         }
1195                 }
1196         }
1197 }
1198
1199 static void
1200 s_calc_nonbonded_force(MDArena *arena)
1201 {
1202         Double *energies = &arena->energies[kVDWIndex];
1203         Double *eenergies = &arena->energies[kElectrostaticIndex];
1204         Vector *forces = &arena->forces[kVDWIndex * arena->mol->natoms];
1205         Vector *eforces = &arena->forces[kElectrostaticIndex * arena->mol->natoms];
1206         s_calc_nonbonded_force_sub(arena, energies, eenergies, forces, eforces, NULL, NULL);
1207 }
1208
1209 static void
1210 s_calc_auxiliary_force(MDArena *arena)
1211 {
1212         Double *energies = &arena->energies[kAuxiliaryIndex];
1213         Vector *forces = &arena->forces[kAuxiliaryIndex * arena->mol->natoms];
1214         Atom *atoms = arena->mol->atoms;
1215         int natoms = arena->mol->natoms;
1216         int i, j;
1217
1218         /*  Centric force - OBSOLETE */
1219 /*      if (arena->centric_potential_force > 0) {
1220                 Vector center = arena->mol->atoms[arena->centric_potential_center].r;
1221                 Double rin = arena->centric_potential_inner_limit;
1222                 Double rout = arena->centric_potential_outer_limit;
1223                 Double k = arena->centric_potential_force;
1224                 for (i = 0; i < natoms; i++) {
1225                         Vector r21;
1226                         Double w1, w2, k0, k1;
1227                         VecSub(r21, atoms[i].r, center);
1228                         w2 = VecLength2(r21);
1229                         w1 = sqrt(w2);
1230                         if (w1 <= rin) {
1231                                 k0 = k1 = 0.0;
1232                         } else if (rout > 0 && w1 > rout) {
1233                                 k0 = k * (rout - rin) * (w1 + w1 - rout + rin);
1234                                 k1 = -2.0 * k * (rout - rin) / w1;
1235                         } else {
1236                                 k0 = k * (w1 - rin) * (w1 - rin);
1237                                 k1 = -2.0 * k * (1 - rin / w1);
1238                         }
1239                         *energies += k0;
1240                         VecScaleSelf(r21, k1);
1241                         VecInc(forces[i], r21);
1242 #if DEBUG
1243                         if (arena->debug_result && arena->debug_output_level > 1) {
1244                                 fprintf(arena->debug_result, "auxiliary(centric) force %d: r=%f, k1=%f, {%f %f %f}\n", i, w1, k1*w1 * INTERNAL2KCAL, r21.x * INTERNAL2KCAL, r21.y * INTERNAL2KCAL, r21.z * INTERNAL2KCAL);
1245                         }
1246 #endif
1247                 }
1248         }
1249         */
1250
1251         /*  Spherical boundary conditions  */
1252         if (arena->spherical_bc_force > 0) {
1253                 Vector center = arena->spherical_bc_center;
1254                 Double rin = arena->spherical_bc_inner_limit;
1255                 Double rout = arena->spherical_bc_outer_limit;
1256                 Double k = arena->spherical_bc_force;
1257                 for (i = 0; i < natoms; i++) {
1258                         Vector r21;
1259                         Double w1, w2, k0, k1;
1260                         if (atoms[i].anchor != NULL)
1261                                 continue;
1262                         VecSub(r21, atoms[i].r, center);
1263                         w2 = VecLength2(r21);
1264                         w1 = sqrt(w2);
1265                         if (w1 <= rin) {
1266                                 k0 = k1 = 0.0;
1267                         } else if (rout > 0 && w1 > rout) {
1268                                 k0 = k * (rout - rin) * (w1 + w1 - rout + rin);
1269                                 k1 = -2.0 * k * (rout - rin) / w1;
1270                         } else {
1271                                 k0 = k * (w1 - rin) * (w1 - rin);
1272                                 k1 = -2.0 * k * (1 - rin / w1);
1273                         }
1274                         *energies += k0;
1275                         VecScaleSelf(r21, k1);
1276                         VecInc(forces[i], r21);
1277                         if (arena->debug_result && arena->debug_output_level > 1) {
1278                                 fprintf(arena->debug_result, "auxiliary(spherical BC) force %d: r=%f, k1=%f, {%f %f %f}\n", i, w1, k1*w1 * INTERNAL2KCAL, r21.x * INTERNAL2KCAL, r21.y * INTERNAL2KCAL, r21.z * INTERNAL2KCAL);
1279                         }
1280                 }
1281         }
1282
1283         /*  Box force  */
1284         if (arena->box_potential_force > 0) {
1285                 Double xsize = arena->box_potential_xsize;
1286                 Double ysize = arena->box_potential_ysize;
1287                 Double zsize = arena->box_potential_zsize;
1288                 Double k = arena->box_potential_force;
1289                 for (i = 0; i < natoms; i++) {
1290                         Vector r = atoms[i].r;
1291                         if (atoms[i].anchor != NULL)
1292                                 continue;
1293                         if (r.x > xsize)
1294                                 r.x -= xsize;
1295                         else if (r.x < -xsize)
1296                                 r.x += xsize;
1297                         else r.x = 0.0;
1298                         if (r.y > ysize)
1299                                 r.y -= ysize;
1300                         else if (r.y < -ysize)
1301                                 r.y += ysize;
1302                         else r.y = 0.0;
1303                         if (r.z > zsize)
1304                                 r.z -= zsize;
1305                         else if (r.z < -zsize)
1306                                 r.z += zsize;
1307                         else r.z = 0.0;
1308                         *energies += k * (r.x * r.x + r.y * r.y + r.z * r.z);
1309                         VecScaleInc(forces[i], r, -2);
1310                         if (arena->debug_result && arena->debug_output_level > 1) {
1311                                 fprintf(arena->debug_result, "auxiliary(box) force %d: {%f %f %f}\n", i, -2*r.x * INTERNAL2KCAL, -2*r.y * INTERNAL2KCAL, -2*r.z * INTERNAL2KCAL);
1312                         }
1313                 }
1314         }
1315
1316         /*  Fix atoms  */
1317         for (i = j = 0; i < natoms; i++) {
1318                 Atom *ap = &atoms[i];
1319                 Vector r21;
1320                 Double w1, w2, k0, k1;
1321                 if (ap->fix_force <= 0.0)
1322                         continue;
1323                 VecSub(r21, ap->r, ap->fix_pos);
1324                 w2 = VecLength2(r21);
1325                 w1 = sqrt(w2);
1326                 k0 = ap->fix_force * w2;
1327                 k1 = -2.0 * ap->fix_force * w1;
1328                 *energies += k0;
1329                 VecScaleSelf(r21, k1);
1330                 VecInc(forces[i], r21);
1331                 j++;
1332                 if (arena->debug_result && arena->debug_output_level > 1) {
1333                         fprintf(arena->debug_result, "auxiliary(fix) force %d: r=%f, k1=%f, {%f %f %f}\n", i, w1, k1*w1 * INTERNAL2KCAL, r21.x * INTERNAL2KCAL, r21.y * INTERNAL2KCAL, r21.z * INTERNAL2KCAL);
1334                 }
1335         }
1336
1337         /*  External forces (used for artificial deformation of the structure)  */
1338         if (arena->nexforces > 0) {
1339                 for (i = 0; i < arena->nexforces; i++) {
1340                         VecInc(forces[i], arena->exforces[i]);
1341                 }
1342         }
1343
1344         /*  Graphite  */
1345         if (arena->graphite != NULL && arena->use_graphite) {
1346                 graphite_force(arena->graphite, arena, energies, forces);
1347         }
1348 }
1349
1350 static void
1351 s_md_debug_output_positions(MDArena *arena)
1352 {
1353         int i;
1354         if (arena->debug_result == NULL || arena->debug_output_level == 0)
1355                 return;
1356         fprintf(arena->debug_result, "\n  Atom positions, velocities, and forces at step %d\n", arena->step);
1357         fprintf(arena->debug_result, "%5s %7s %7s %7s %7s %7s %7s %7s %7s %7s\n", "No.", "x", "y", "z", "vx", "vy", "vz", "fx", "fy", "fz");
1358         for (i = 0; i < arena->mol->natoms; i++) {
1359                 Atom *ap = &arena->mol->atoms[i];
1360                 fprintf(arena->debug_result, "%5d %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n", i+1, ap->r.x, ap->r.y, ap->r.z, ap->v.x, ap->v.y, ap->v.z, ap->f.x * INTERNAL2KCAL, ap->f.y * INTERNAL2KCAL, ap->f.z * INTERNAL2KCAL);
1361         }
1362 }
1363
1364 void
1365 calc_force(MDArena *arena)
1366 {
1367         Int i, j, natoms;
1368         Vector *ff, *fa;
1369         Molecule *mol = arena->mol;
1370         Int doSurface = 0, doEwald = 0;
1371         Atom *ap;
1372
1373         natoms = mol->natoms;
1374
1375         /*  Clear the energy/force storage  */
1376         for (i = 0; i < natoms; i++)
1377                 VecZero(mol->atoms[i].f);
1378
1379         memset(arena->energies, 0, sizeof(Double) * kSlowIndex);
1380         memset(arena->forces, 0, sizeof(Vector) * kSlowIndex * natoms);
1381         arena->total_energy = 0.0;
1382         arena->alchem_energy = 0.0;
1383
1384         if (arena->step == arena->start_step || arena->step % arena->surface_potential_freq == 0) {
1385                 doSurface = 1;
1386                 arena->energies[kSurfaceIndex] = 0.0;
1387                 memset(arena->forces + kSurfaceIndex * natoms, 0, sizeof(Vector) * natoms);
1388         }
1389
1390         if (arena->step == arena->start_step || arena->step % arena->surface_potential_freq == 0) {
1391                 doEwald = 1;
1392                 arena->energies[kPMEIndex] = 0.0;
1393                 memset(arena->forces + kPMEIndex * natoms, 0, sizeof(Vector) * natoms);
1394         }
1395
1396         s_calc_bond_force(arena);
1397         s_calc_angle_force(arena);
1398         s_calc_dihedral_force(arena);
1399         s_calc_improper_force(arena);
1400         s_calc_nonbonded_force(arena);
1401         s_calc_auxiliary_force(arena);
1402
1403         if (doEwald && arena->use_ewald != 0)
1404                 calc_ewald_force(arena);
1405
1406         if (doSurface && arena->probe_radius > 0.0) {
1407         /*      if (arena->sphere_points == 0)
1408                         calc_surface_force(arena);
1409                 else
1410                         calc_surface_force_2(arena); */
1411                 calc_surface_force(arena);
1412         }
1413
1414         /*  Distribute forces on pi-anchor atoms to their components  */
1415         for (i = 0; i < natoms; i++) {
1416                 Int k, n, *ip;
1417                 Double w;
1418                 ap = &(mol->atoms[i]);
1419                 if (ap->anchor == NULL)
1420                         continue;
1421                 ip = AtomConnectData(&ap->anchor->connect);
1422                 n = ap->anchor->connect.count;
1423                 for (k = 0; k < n; k++) {
1424                         ff = &arena->forces[ip[k]];
1425                         fa = &arena->forces[i];
1426                         w = ap->anchor->coeffs[k];
1427                         for (j = 0; j < kKineticIndex; j++) {
1428                                 VecScaleInc(*ff, *fa, w);
1429                                 ff += natoms;
1430                                 fa += natoms;
1431                         }
1432                 }
1433                 fa = &arena->forces[i];
1434                 for (j = 0; j < kKineticIndex; j++) {
1435                         VecZero(*fa);
1436                         fa += natoms;
1437                 }
1438         }
1439
1440         /*  Sum up all partial forces and energies  */
1441         arena->total_energy = 0.0;
1442         for (i = 0; i < kKineticIndex; i++)
1443                 arena->total_energy += arena->energies[i];
1444
1445         for (i = 0; i < natoms; i++) {
1446                 if (mol->atoms[i].mm_exclude)
1447                         continue;
1448                 fa = &mol->atoms[i].f;
1449                 ff = &arena->forces[i];
1450                 for (j = 0; j < kKineticIndex; j++) {
1451                         VecInc(*fa, *ff);
1452                         ff += natoms;
1453                 }
1454         }
1455
1456         /*  The total (internal) force must be zero  */
1457 /*      {
1458                 Vector f;
1459                 Double w;
1460                 VecZero(f);
1461                 for (i = 0; i < natoms; i++) {
1462                         Vector *fp = &(mol->atoms[i].f);
1463                         VecInc(f, *fp);
1464                 }
1465                 w = (natoms > 0 ? 1.0 / natoms : 0.0);
1466                 VecScaleSelf(f, w);
1467                 for (i = 0; i < natoms; i++) {
1468                         Vector *fp = &(mol->atoms[i].f);
1469                         VecDec(*fp, f);
1470                 }
1471         } */
1472
1473         s_md_debug_output_positions(arena);
1474 }
1475
1476 void
1477 calc_pair_interaction(MDArena *arena, const MDGroupFlags *group_flags_1, const MDGroupFlags *group_flags_2)
1478 {
1479         Vector *forces, *eforces;
1480         if (arena->pair_forces != NULL) {
1481                 forces = arena->pair_forces;
1482                 eforces = forces + arena->mol->natoms;
1483                 memset(forces, 0, sizeof(Vector) * arena->mol->natoms * 2);
1484         } else forces = eforces = NULL;
1485         arena->pair_energies[0] = arena->pair_energies[1] = 0.0;
1486         s_calc_nonbonded_force_sub(arena, &arena->pair_energies[0], &arena->pair_energies[1], forces, eforces, group_flags_1, group_flags_2);
1487 }