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