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]].occupancy == 0.0 || ap[bonds[1]].occupancy == 0.0)
  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 void
 100 s_calc_angle_force(MDArena *arena)
 101 {
 102         Atom *ap = arena->mol->atoms;
 103         Int nangles = arena->mol->nangles;
 104         Int *angles = arena->mol->angles;
 105         Int *angle_par_i = arena->angle_par_i;
 106         AnglePar *angle_pars = arena->par->anglePars;
 107 /*      Int *angle_uniq = (arena->nsyms > 0 ? arena->sym_angle_uniq : NULL); */
 108         Double *energies = &arena->energies[kAngleIndex];
 109         Vector *forces = &arena->forces[kAngleIndex * arena->mol->natoms];
 110         int i;
 111         for (i = 0; i < nangles; i++, angles += 3, angle_par_i++) {
 112                 Vector r21, r23, f1, f2, v1;
 113                 Double k0, k1, w1, w2, w3;
 114                 Double cost, sint, t;
 115                 AnglePar *anp;
 116                 if (*angle_par_i < 0)
 117                         continue;  /*  Ignore this entry  */
 118         /*      if (angle_uniq != NULL && angle_uniq[i] >= 0)
 119                         continue; */ /*  Non-unique angle  */
 120                 if (ap[angles[0]].occupancy == 0.0 || ap[angles[1]].occupancy == 0.0 || ap[angles[2]].occupancy == 0.0)
 121                         continue;  /*  Skip non-occupied atoms  */
 122                 anp = angle_pars + *angle_par_i;
 123                 VecSub(r21, ap[angles[0]].r, ap[angles[1]].r);
 124                 VecSub(r23, ap[angles[2]].r, ap[angles[1]].r);
 125                 w1 = 1.0 / VecLength(r21);
 126                 w2 = 1.0 / VecLength(r23);
 127                 VecScaleSelf(r21, w1);
 128                 VecScaleSelf(r23, w2);
 129                 cost = VecDot(r21, r23);
 130                 if (cost > 0.999999 || cost < -0.999999) {
 131                 /*      printf("Cannot handle linear angle %d-%d-%d: skipped.\n", angles[0]+1, angles[1]+1, angles[2]+1); */
 132                         continue;  /*  Cannot handle linear angle  */
 133                 }
 134                 sint = sqrt(1.0 - cost * cost);
 135         /*      if (sint < 1e-5 && sint > -1e-5)
 136                         continue;  *//*  Cannot handle linear angle  */
 137                 t = atan2(sint, cost) - anp->a0;
 138                 k0 = anp->k * t * t;
 139
 140                 if (sint < 0.1 && sint > -0.1) {
 141                         /* ---- Method 1 ---- */
 142                         /* This is slower than method 2, but it is safer when sint is close to 0 */
 143                         k1 = -2 * anp->k * t;
 144                         VecCross(v1, r21, r23);
 145                         VecCross(f1, r21, v1);
 146                         w3 = w1 * k1 / VecLength(f1);
 147                         if (!isfinite(w3))
 148                                 continue;
 149                         VecScaleSelf(f1, w3);
 150                         VecCross(f2, v1, r23);
 151                         w3 = w2 * k1 / VecLength(f2);
 152                         if (!isfinite(w3))
 153                                 continue;
 154                         VecScaleSelf(f2, w3);
 155                 } else {
 156                         /* ---- Method 2 ---- */
 157                         k1 = -2 * anp->k * t / sint;
 158                         VecScale(f1, r21, cost);
 159                         VecDec(f1, r23);
 160                         w3 = w1 * k1;
 161                         VecScaleSelf(f1, w3);
 162                         VecScale(f2, r23, cost);
 163                         VecDec(f2, r21);
 164                         w3 = w2 * k1;
 165                         VecScaleSelf(f2, w3);
 166                 }
 167
 168                 *energies += k0;
 169                 VecInc(forces[angles[0]], f1);
 170                 VecInc(forces[angles[2]], f2);
 171                 VecInc(f1, f2);
 172                 VecDec(forces[angles[1]], f1);
 173
 174                 if (arena->debug_result && arena->debug_output_level > 1) {
 175                         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, (t+anp->a0)*180/PI, anp->a0*180/PI, k1, f1.x, f1.y, f1.z, f2.x, f2.y, f2.z);
 176                 }
 177
 178         }
 179 }
 180
 181 static void
 182 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)
 183 {
 184         int i;
 185         if (arena->mol->nsyms == 0)
 186                 dihedral_uniq = NULL;  /*  Ignore the symmetry info  */
 187         for (i = 0; i < ndihedrals; i++, dihedrals += 4, dihedral_par_i++) {
 188                 Vector r21, r32, r43;
 189                 Vector v1, v2, v3;
 190                 Double w1, w2, w3, k0, k1;
 191                 Double cosphi, sinphi, phi;
 192                 Vector f1, f2, f3;
 193                 TorsionPar *tp;
 194                 int n;
 195
 196                 if (*dihedral_par_i < 0)
 197                         continue;  /*  Ignore this entry  */
 198         /*      if (dihedral_uniq != NULL && dihedral_uniq[i] >= 0)
 199                         continue;  *//*  Non-unique dihedral  */
 200                 if (ap[dihedrals[0]].occupancy == 0.0 || ap[dihedrals[1]].occupancy == 0.0 || ap[dihedrals[2]].occupancy == 0.0 || ap[dihedrals[3]].occupancy == 0.0)
 201                         continue;  /*  Skip non-occupied atoms  */
 202                 else tp = dihedral_pars + *dihedral_par_i;
 203                 VecSub(r21, ap[dihedrals[0]].r, ap[dihedrals[1]].r);
 204                 VecSub(r32, ap[dihedrals[1]].r, ap[dihedrals[2]].r);
 205                 VecSub(r43, ap[dihedrals[2]].r, ap[dihedrals[3]].r);
 206                 VecCross(v1, r21, r32);
 207                 VecCross(v2, r32, r43);
 208                 VecCross(v3, r32, v1);
 209                 w1 = VecLength(v1);
 210                 w2 = VecLength(v2);
 211                 w3 = VecLength(v3);
 212                 if (w1 < 1e-5 || w2 < 1e-5 || w3 < 1e-5)
 213                         continue;  /*  The dihedral cannot be defined  */
 214                 w1 = 1.0/w1;
 215                 w2 = 1.0/w2;
 216                 w3 = 1.0/w3;
 217                 VecScaleSelf(v1, w1);
 218                 VecScaleSelf(v2, w2);
 219                 VecScaleSelf(v3, w3);
 220
 221                 /*  The dihedral angle value  */
 222                 cosphi = VecDot(v1, v2);
 223                 sinphi = VecDot(v3, v2);
 224                 phi = -atan2(sinphi, cosphi);
 225
 226                 /*  Repeat for multiple dihedral terms  */
 227                 k0 = k1 = 0.0;
 228                 for (n = 0; n < tp->mult; n++) {
 229                         Double k = tp->k[n];
 230                         Double phi0 = tp->phi0[n];
 231                         Int period = tp->period[n];
 232                         if (period > 0) {
 233                                 k0 += k * (1 + cos(period * phi + phi0));
 234                                 k1 -= period * k * sin(period * phi + phi0);
 235                         } else {
 236                                 /*  A simple quadratic term  */
 237                                 phi -= phi0;
 238                                 if (phi < -PI)
 239                                         phi += 2 * PI;
 240                                 else if (phi > PI)
 241                                         phi -= 2 * PI;
 242                                 k0 += k * phi * phi;
 243                                 k1 += 2 * k * phi;
 244                         }
 245                 }
 246
 247                 if (sinphi < -0.1 || sinphi > 0.1) {
 248                         /*  The sin form  */
 249                         Vector v4, v5, vw0;
 250                         VecScale(v4, v1, cosphi);
 251                         VecDec(v4, v2);
 252                         VecScaleSelf(v4, w1);
 253                         VecScale(v5, v2, cosphi);
 254                         VecDec(v5, v1);
 255                         VecScaleSelf(v5, w2);
 256                         k1 /= sinphi;
 257                         VecCross(f1, r32, v4);
 258                         VecCross(f3, v5, r32);
 259                         VecCross(f2, v4, r21);
 260                         VecCross(vw0, r43, v5);
 261                         VecInc(f2, vw0);
 262                         VecScaleSelf(f1, k1);
 263                         VecScaleSelf(f2, k1);
 264                         VecScaleSelf(f3, k1);
 265                 } else {
 266                         /*  The cos form  */
 267                         Vector v6, v7, vw1, vw2;
 268                         VecScale(v6, v3, sinphi);
 269                         VecDec(v6, v2);
 270                         VecScaleSelf(v6, w3);
 271                         VecScale(v7, v2, sinphi);
 272                         VecDec(v7, v3);
 273                         VecScaleSelf(v7, w2);
 274                         k1 = -k1 / cosphi;
 275                         VecCross(vw1, v6, r32);
 276                         VecCross(f1, r32, vw1);
 277                         VecCross(f3, v7, r32);
 278                         VecCross(f2, r43, v7);
 279                         VecCross(vw1, r32, r21);
 280                         VecCross(vw2, v6, vw1);
 281                         VecInc(f2, vw2);
 282                         VecCross(vw1, r32, v6);
 283                         VecCross(vw2, r21, vw1);
 284                         VecInc(f2, vw2);
 285                         VecScaleSelf(f1, k1);
 286                         VecScaleSelf(f2, k1);
 287                         VecScaleSelf(f3, k1);
 288                 }
 289         /*      if (dihedral_uniq != NULL)
 290                         k0 *= -dihedral_uniq[i]; */
 291                 *energies += k0;
 292                 VecInc(forces[dihedrals[0]], f1);
 293                 VecDec(f1, f2);
 294                 VecDec(forces[dihedrals[1]], f1);
 295                 VecDec(f2, f3);
 296                 VecDec(forces[dihedrals[2]], f2);
 297                 VecDec(forces[dihedrals[3]], f3);
 298
 299                 if (arena->debug_result && arena->debug_output_level > 1) {
 300                         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);
 301                 }
 302         }
 303 }
 304
 305 static void
 306 s_calc_dihedral_force(MDArena *arena)
 307 {
 308         Molecule *mol = arena->mol;
 309         Parameter *par = arena->par;
 310         Double *energies = &arena->energies[kDihedralIndex];
 311         Vector *forces = &arena->forces[kDihedralIndex * arena->mol->natoms];
 312         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);
 313 }
 314
 315 static void
 316 s_calc_improper_force(MDArena *arena)
 317 {
 318         Molecule *mol = arena->mol;
 319         Parameter *par = arena->par;
 320         Double *energies = &arena->energies[kImproperIndex];
 321         Vector *forces = &arena->forces[kImproperIndex * arena->mol->natoms];
 322         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);
 323 }
 324
 325 /*  ==================================================================== */
 326 /*  md_check_verlet_list: only for debugging the verlet list generation  */
 327 /*  ==================================================================== */
 328
 329 typedef struct md_check_record {
 330         Int i, j, dx, dy, n;
 331         Double len;
 332 } md_check_record;
 333
 334 static int
 335 s_md_check_verlet_comparator(const void *ap, const void *bp)
 336 {
 337         Double d = ((const md_check_record *)ap)->len - ((const md_check_record *)bp)->len;
 338         return (d < 0 ? -1 : (d > 0 ? 1 : 0));
 339 }
 340
 341 void
 342 md_check_verlet_list(MDArena *arena)
 343 {
 344         int i, j, k;
 345         int dx, dy, dz, nn;
 346         int ndx, ndy, ndz;
 347         Atom *api, *apj;
 348         Vector cell_offsets[27];
 349         XtalCell *cell = arena->mol->cell;
 350         md_check_record *cr;
 351         int ncr;
 352
 353         ndx = (arena->periodic_a ? 1 : 0);
 354         ndy = (arena->periodic_b ? 1 : 0);
 355         ndz = (arena->periodic_c ? 1 : 0);
 356         if (cell != NULL && (ndx != 0 || ndy != 0 || ndz != 0)) {
 357                 nn = 0;
 358                 for (i = -1; i <= 1; i++) {
 359                         for (j = -1; j <= 1; j++) {
 360                                 for (k = -1; k <= 1; k++) {
 361                                         VecZero(cell_offsets[nn]);
 362                                         VecScaleInc(cell_offsets[nn], cell->axes[0], i);
 363                                         VecScaleInc(cell_offsets[nn], cell->axes[1], j);
 364                                         VecScaleInc(cell_offsets[nn], cell->axes[2], k);
 365                                         nn++;
 366                                 }
 367                         }
 368                 }
 369         }
 370         /*  Debugging is done with x-y 2-dimensional system  */
 371         ncr = arena->mol->natoms * (arena->mol->natoms - 1) / 2 * 9;
 372         cr = (md_check_record *)malloc(sizeof(md_check_record) * ncr);
 373         k = 0;
 374         for (i = 0, api = arena->mol->atoms; i < arena->mol->natoms; i++, api = ATOM_NEXT(api)) {
 375                 for (j = i + 1, apj = ATOM_AT_INDEX(arena->mol->atoms, j); j < arena->mol->natoms; j++, apj = ATOM_NEXT(apj)) {
 376                         Vector dr;
 377                         VecSub(dr, apj->r, api->r);
 378                         for (dx = -1; dx <= 1; dx++) {
 379                                 for (dy = -1; dy <= 1; dy++) {
 380                                         Vector dr2 = dr;
 381                                         nn = dx * 9 + dy * 3 + 13;
 382                                         VecInc(dr2, cell_offsets[nn]);
 383                                         cr[k].i = i;
 384                                         cr[k].j = j;
 385                                         cr[k].dx = dx;
 386                                         cr[k].dy = dy;
 387                                         cr[k].n = -1;
 388                                         cr[k].len = VecLength(dr2);
 389                                         k++;
 390                                 }
 391                         }
 392                 }
 393         }
 394         for (k = 0; k < arena->nverlets; k++) {
 395                 i = arena->verlets[k].n1;
 396                 j = arena->verlets[k].n2;
 397                 if (i > j) {
 398                         j = i;
 399                         i = arena->verlets[k].n2;
 400                 }
 401                 dx = arena->verlets[k].symop.dx;
 402                 dy = arena->verlets[k].symop.dy;
 403                 nn = (i * arena->mol->natoms - i * (i + 1) / 2 + (j - i) - 1) * 9 + (dx + 1) * 3 + dy + 1;
 404                 if (cr[nn].i != i || cr[nn].j != j || cr[nn].dx != dx || cr[nn].dy != dy)
 405                         fprintf(arena->log_result, "Debug: internal inconsistency\n");
 406                 cr[nn].n = k;
 407         }
 408
 409         mergesort(cr, ncr, sizeof(md_check_record), s_md_check_verlet_comparator);
 410         for (i = 0; i < ncr; i++) {
 411                 Double len2;
 412                 len2 = -1;
 413                 if (cr[i].n >= 0)
 414                         len2 = arena->verlets[cr[i].n].length;
 415                 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);
 416                 if (cr[i].len > arena->pairlist_distance)
 417                         break;
 418         }
 419         while (i < ncr) {
 420                 if (cr[i].n != -1) {
 421                         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);
 422                 }
 423                 i++;
 424         }
 425         fflush(arena->log_result);
 426         free(cr);
 427 }
 428
 429 /*  ==================================================================== */
 430
 431
 432 /*  Update the Verlet list (for non-bonding interaction)  */
 433 static void
 434 s_make_verlet_list(MDArena *arena)
 435 {
 436         int i, j, k, n, nn, natoms;
 437         int dx, dy, dz, ndx, ndy, ndz;
 438         Molecule *mol = arena->mol;
 439         Atom *atoms = mol->atoms;
 440         Atom *api, *apj;
 441         MDVerlet *vl = arena->verlets;
 442         Vector *vdr = arena->verlets_dr;
 443         Double limit;
 444         Byte use_sym;
 445         XtalCell *cell = mol->cell;
 446         Vector cell_offsets[27];
 447         Parameter *par = arena->par;
 448         Double vdw_cutoff;
 449
 450         natoms = mol->natoms;
 451         for (i = 0; i < natoms; i++)
 452                 VecZero(vdr[i]);
 453
 454         ndx = (arena->periodic_a ? 1 : 0);
 455         ndy = (arena->periodic_b ? 1 : 0);
 456         ndz = (arena->periodic_c ? 1 : 0);
 457
 458         if (cell != NULL && (ndx != 0 || ndy != 0 || ndz != 0)) {
 459                 nn = 0;
 460                 for (i = -1; i <= 1; i++) {
 461                         for (j = -1; j <= 1; j++) {
 462                                 for (k = -1; k <= 1; k++) {
 463                                         VecZero(cell_offsets[nn]);
 464                                         VecScaleInc(cell_offsets[nn], cell->axes[0], i);
 465                                         VecScaleInc(cell_offsets[nn], cell->axes[1], j);
 466                                         VecScaleInc(cell_offsets[nn], cell->axes[2], k);
 467                                         nn++;
 468                                 }
 469                         }
 470                 }
 471         }
 472
 473         limit = arena->pairlist_distance * arena->pairlist_distance;
 474         n = 0;
 475         use_sym = (arena->mol->nsyms > 0);
 476         for (i = 0, api = atoms; i < natoms; i++, api++) {
 477                 MDExclusion *exinfo = arena->exinfo + i;
 478                 Int index4 = (exinfo + 1)->index0;
 479                 Int vdw_idx1 = arena->vdw_par_i[i];
 480                 Int vdw_idx, vdw_idx2;
 481                 arena->verlet_i[i] = n;
 482
 483                 for (j = i + 1, apj = atoms + j; j < natoms; j++, apj++) {
 484                         Vector rij;
 485                         Double lenij2;
 486                         Int *ip, index0;
 487                         int exflag = 0;
 488                         int mult = 1;
 489                         int dxbase, dybase, dzbase;
 490
 491                         /*  Fixed atoms  */
 492                         if (api->fix_force < 0 && apj->fix_force < 0)
 493                                 continue;
 494
 495                         /*  Non-occupied atoms  */
 496                         if (api->occupancy == 0.0 || apj->occupancy == 0.0)
 497                                 continue;
 498
 499                         VecSub(rij, apj->r, api->r);
 500
 501                         /*  Calculate the cell offset for the nearest neighbor  */
 502                         /*  NOTE: the offset is calculated independently for each axis. This may result
 503                             in unexpected choice when the angles between the axes are far from 90 deg */
 504                         if (apj->periodic_exclude == 0 && cell != NULL) {
 505                                 TransformPtr rtp = cell->rtr;
 506                                 TransformPtr tp = cell->tr;
 507                                 Double w;
 508                                 if (arena->periodic_a) {
 509                                         w = rtp[0] * rij.x + rtp[1] * rij.y + rtp[2] * rij.z;
 510                                         dxbase = floor(-w + 0.5);
 511                                 } else dxbase = 0;
 512                                 if (arena->periodic_b) {
 513                                         w = rtp[3] * rij.x + rtp[4] * rij.y + rtp[5] * rij.z;
 514                                         dybase = floor(-w + 0.5);
 515                                 } else dybase = 0;
 516                                 if (arena->periodic_c) {
 517                                         w = rtp[6] * rij.x + rtp[7] * rij.y + rtp[8] * rij.z;
 518                                         dzbase = floor(-w + 0.5);
 519                                 } else dzbase = 0;
 520                                 rij.x += tp[0] * dxbase + tp[1] * dybase + tp[2] * dzbase;
 521                                 rij.y += tp[3] * dxbase + tp[4] * dybase + tp[5] * dzbase;
 522                                 rij.z += tp[6] * dxbase + tp[7] * dybase + tp[8] * dzbase;
 523                         } else dxbase = dybase = dzbase = 0;
 524
 525                         /*  Non unique atom pair  */
 526                 /*      if (use_sym && api->symop.alive && apj->symop.alive)
 527                                 continue; */
 528
 529                         /*  Check the specific cutoff table for (i, j) pair  */
 530                         vdw_idx2 = arena->vdw_par_i[j];
 531                         if (vdw_idx1 < 0 || vdw_idx2 < 0)
 532                                 vdw_idx = par->nvdwPars * par->nvdwPars;  /*  A null record  */
 533                         else if (vdw_idx1 < vdw_idx2)
 534                                 vdw_idx = vdw_idx1 * par->nvdwPars + vdw_idx2;
 535                         else
 536                                 vdw_idx = vdw_idx2 * par->nvdwPars + vdw_idx1;
 537                         vdw_cutoff = arena->cutoff;
 538                         for (k = par->nvdwCutoffPars - 1; k >= 0; k--) {
 539                                 VdwCutoffPar *cr = par->vdwCutoffPars + k;
 540                                 if (cr->type == 0) {
 541                                         if (((cr->n1 < 0 || cr->n1 == api->type) && (cr->n2 < 0 || cr->n2 == apj->type)) ||
 542                                                 ((cr->n1 < 0 || cr->n1 == apj->type) && (cr->n2 < 0 || cr->n2 == api->type))) {
 543                                                 vdw_cutoff = cr->cutoff;
 544                                                 break;
 545                                         }
 546                                 } else {
 547                                         if ((cr->n1 <= i && i <= cr->n2 && cr->n3 <= j && j <= cr->n4)||
 548                                                 (cr->n3 <= i && i <= cr->n4 && cr->n1 <= j && j <= cr->n2)) {
 549                                                 vdw_cutoff = cr->cutoff;
 550                                                 break;
 551                                         }
 552                                 }
 553                         }
 554                         if (vdw_cutoff < 0) {
 555                                 /*  A negative value of specific cutoff means "cutoff at r_eq * (-specific_cutoff)  */
 556                                 MDVdwCache *cp = &(arena->vdw_cache[vdw_idx]);
 557                                 Double r_eq = pow(cp->par.A / cp->par.B * 2.0, 1.0/6.0);
 558                                 vdw_cutoff = r_eq * (-vdw_cutoff);
 559                         }
 560
 561                         /*  Search for pairs, taking periodicity into account  */
 562                         for (dx = -ndx; dx <= ndx; dx++) {
 563                                 for (dy = -ndy; dy <= ndy; dy++) {
 564                                         for (dz = -ndz; dz <= ndz; dz++) {
 565                                                 Vector rij0 = rij;
 566                                                 nn = dx * 9 + dy * 3 + dz + 13;
 567                                                 if (nn == 13) {
 568                                                         /*  Pair within the unit cell  */
 569                                                         /*  Is this pair to be excluded?  */
 570                                                         for (index0 = exinfo->index0, ip = arena->exlist + index0; index0 < index4; index0++, ip++) {
 571                                                                 if (*ip == j)
 572                                                                         break;
 573                                                         }
 574                                                         if (index0 < exinfo->index3)
 575                                                                 continue;  /*  Special exclusion, 1-2, 1-3  */
 576                                                         if (index0 < index4)
 577                                                                 exflag = 1;  /*  1-4 interaction  */
 578                                                 } else if (apj->periodic_exclude) {
 579                                                         continue;
 580                                                 } else {
 581                                                         VecInc(rij0, cell_offsets[nn]);
 582                                                         exflag = 0;
 583                                                 }
 584
 585                                                 lenij2 = VecLength2(rij0);
 586                                                 if (lenij2 <= limit) {
 587                                                         MDVerlet *vlp;
 588                                                         if (n >= arena->max_nverlets) {
 589                                                                 arena->max_nverlets += 32;
 590                                                                 vl = (MDVerlet *)realloc(vl, sizeof(MDVerlet) * arena->max_nverlets);
 591                                                                 if (vl == NULL)
 592                                                                         md_panic(arena, "Low memory");
 593                                                         }
 594                                                         vlp = &vl[n];
 595                                                         vlp->vdw_type = (exflag ? 1 : 0);
 596                                                         vlp->mult = mult;
 597                                                         vlp->symop.dx = dx + dxbase;
 598                                                         vlp->symop.dy = dy + dybase;
 599                                                         vlp->symop.dz = dz + dzbase;
 600                                                         vlp->symop.sym = 0;
 601                                                         vlp->symop.alive = (vlp->symop.dx != 0 || vlp->symop.dy != 0 || vlp->symop.dz != 0);
 602                                                         vlp->index = vdw_idx;
 603                                                         vlp->n1 = i;
 604                                                         vlp->n2 = j;
 605                                                         vlp->vdw_cutoff = vdw_cutoff;
 606                                                         vlp->length = sqrt(lenij2);
 607                                                         n++;
 608                                                 } /* end if lenij2 <= limit */
 609                                         } /* end loop dz */
 610                                 } /* end loop dy */
 611                         } /* end loop dx */
 612                 } /* end loop j */
 613
 614         } /* end loop i */
 615         arena->verlet_i[natoms] = n;
 616         arena->nverlets = n;
 617         arena->verlets = vl;
 618         arena->last_verlet_step = arena->step;
 619
 620         if (arena->debug_result && arena->debug_output_level > 2) {
 621                 fprintf(arena->debug_result, "\n  Verlet list at step %d\n", arena->step);
 622                 fprintf(arena->debug_result, "  {atom1 atom2 vdw_type (=1 if 1-4 bonded) vdw_index}\n");
 623                 for (i = 0; i < arena->nverlets; i++) {
 624                         fprintf(arena->debug_result, "{%d %d %d %d}%c", vl[i].n1+1, vl[i].n2+1, vl[i].vdw_type, vl[i].index, (i % 4 == 3 ? '\n' : ' '));
 625                 }
 626                 if (i % 4 != 0)
 627                         fprintf(arena->debug_result, "\n");
 628         }
 629 }
 630
 631 /*  Calculate the nonbonded force  */
 632 /*  group_flags[] is used for calculation of pair-specific interaction between
 633     two groups of atoms  */
 634 static void
 635 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)
 636 {
 637         int i;
 638         Vector *vdr;
 639         Double limit, elimit, dielec_r;
 640         Double lambda, dlambda;
 641
 642         MDVerlet *vl;
 643 /*      MDVdwCache *vdw_cache = arena->vdw_cache; */
 644         Atom *atoms = arena->mol->atoms;
 645         XtalCell *cell = arena->mol->cell;
 646
 647         /*  Check whether the Verlet list needs update  */
 648         if (arena->last_verlet_step >= 0) {
 649                 i = arena->mol->natoms - 1;
 650                 vdr = arena->verlets_dr + i;
 651                 if (arena->cutoff > arena->electro_cutoff)
 652                         limit = (arena->pairlist_distance - arena->cutoff) * 0.5;
 653                 else
 654                         limit = (arena->pairlist_distance - arena->electro_cutoff) * 0.5;
 655                 limit = limit * limit;
 656                 for ( ; i >= 0; i--, vdr--) {
 657                         if (VecLength2(*vdr) >= limit)
 658                                 break;
 659                 }
 660         } else i = 0;
 661         if (i >= 0)
 662                 s_make_verlet_list(arena);
 663
 664         /*  Calculate the non-bonded interaction for each pair in the Verlet list  */
 665 /*      limit = arena->cutoff * arena->cutoff; */
 666         elimit = arena->electro_cutoff * arena->electro_cutoff;
 667         dielec_r = COULOMBIC / arena->dielectric;
 668         for (i = arena->nverlets - 1, vl = arena->verlets + i; i >= 0; i--, vl--) {
 669                 Double A, B, vofs, k0, k1;
 670                 MDVdwCache *vp = arena->vdw_cache + vl->index;
 671                 Atom *ap1, *ap2;
 672                 Vector rij, fij;
 673                 Double r2, w2, w6, w12;
 674                 if (group_flags_1 != NULL && group_flags_2 != NULL) {
 675                         if (!(get_group_flag(group_flags_1, vl->n1) && get_group_flag(group_flags_2, vl->n2))
 676                         && !(get_group_flag(group_flags_1, vl->n2) && get_group_flag(group_flags_2, vl->n1)))
 677                                 continue;
 678                 }
 679
 680                 if (arena->nalchem_flags > 0) {
 681                         char c1, c2;
 682                         if (vl->n1 < arena->nalchem_flags)
 683                                 c1 = arena->alchem_flags[vl->n1];
 684                         else c1 = 0;
 685                         if (vl->n2 < arena->nalchem_flags)
 686                                 c2 = arena->alchem_flags[vl->n2];
 687                         else c2 = 0;
 688                         if ((c1 == 1 && c2 == 2) || (c1 == 2 && c2 == 1))
 689                                 continue;
 690                         if (c1 == 1 || c2 == 1) {
 691                                 lambda = (1.0 - arena->alchem_lambda);
 692                                 dlambda = -arena->alchem_dlambda;
 693                         } else if (c1 == 2 || c2 == 2) {
 694                                 lambda = arena->alchem_lambda;
 695                                 dlambda = arena->alchem_dlambda;
 696                         } else {
 697                                 lambda = 1.0;
 698                                 dlambda = 0.0;
 699                         }
 700                 } else {
 701                         lambda = 1.0;
 702                         dlambda = 0.0;
 703                 }
 704
 705                 if (vl->vdw_type == 1) {
 706                         A = vp->par.A14;
 707                         B = vp->par.B14;
 708                 /*      vofs = vp->vcut14; */
 709                 } else {
 710                         A = vp->par.A;
 711                         B = vp->par.B;
 712                 /*      vofs = vp->vcut; */
 713                 }
 714                 ap1 = &atoms[vl->n1];
 715                 ap2 = &atoms[vl->n2];
 716                 rij = ap2->r;
 717                 if (vl->symop.alive) {
 718                         VecScaleInc(rij, cell->axes[0], vl->symop.dx);
 719                         VecScaleInc(rij, cell->axes[1], vl->symop.dy);
 720                         VecScaleInc(rij, cell->axes[2], vl->symop.dz);
 721                 }
 722                 VecDec(rij, ap1->r);
 723                 r2 = VecLength2(rij);
 724                 if (r2 >= elimit)
 725                         continue;
 726                 limit = vl->vdw_cutoff * vl->vdw_cutoff;
 727                 if (r2 >= limit)
 728                         continue;
 729
 730                 fij.x = fij.y = fij.z = 0.0;
 731
 732                 /*  Coulombic force  */
 733                 w12 = ap1->charge * ap2->charge * dielec_r;
 734                 if (w12 != 0.0) {
 735                         if (arena->use_xplor_shift) {
 736                                 w2 = 1.0 / sqrt(r2);
 737                                 k0 = r2 / elimit - 1.0;
 738                                 k0 = k0 * k0;
 739                                 k0 = w12 * k0 * w2;
 740                                 k1 = (3.0 * r2 / elimit - 2.0) / elimit - 1.0 / r2;
 741                                 k1 = -w12 * k1 * w2;
 742                         } else {
 743                                 w2 = 1.0 / sqrt(r2);
 744                                 w6 = w2 / r2;
 745                                 k0 = w12 * w2;
 746                                 k1 = w12 * w6;
 747                                 vofs = -w12 / arena->electro_cutoff;
 748                                 k0 += vofs;
 749                         }
 750                         if (vl->vdw_type == 1) {
 751                                 k0 *= arena->scale14_elect;
 752                                 k1 *= arena->scale14_elect;
 753                         }
 754                         VecScale(fij, rij, k1);
 755                         *eenergies += k0 / vl->mult;
 756                         if (eforces != NULL) {
 757                                 VecDec(eforces[vl->n1], fij);
 758                                 VecInc(eforces[vl->n2], fij);
 759                         }
 760                         if (arena->debug_result && arena->debug_output_level > 1) {
 761                                 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);
 762                         }
 763                 }
 764
 765                 /*  van der Waals force  */
 766                 w2 = 1.0 / r2;
 767                 w6 = w2 * w2 * w2;
 768                 w12 = w6 * w6;
 769                 k0 = A * w12 - B * w6;
 770                 k1 = 6 * w2 * (2.0 * A * w12 - B * w6);
 771                 w2 = 1.0 / limit;
 772                 w6 = w2 * w2 * w2;
 773                 w12 = w6 * w6;
 774                 vofs = -A * w12 + B * w6;
 775                 k0 += vofs;
 776                 if (vl->vdw_type == 1) {
 777                         k0 *= arena->scale14_vdw;
 778                         k1 *= arena->scale14_vdw;
 779                 }
 780                 k0 /= vl->mult;
 781                 k1 *= lambda;
 782                 VecScale(fij, rij, k1);
 783                 *energies += k0 * lambda;
 784                 if (forces != NULL) {
 785                         VecDec(forces[vl->n1], fij);
 786                         VecInc(forces[vl->n2], fij);
 787                 }
 788                 if (dlambda != 0.0)
 789                         arena->alchem_energy += k0 * dlambda;
 790
 791                 if (arena->debug_result && arena->debug_output_level > 1) {
 792                         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);
 793                 }
 794         }
 795 }
 796
 797 static void
 798 s_calc_nonbonded_force(MDArena *arena)
 799 {
 800         Double *energies = &arena->energies[kVDWIndex];
 801         Double *eenergies = &arena->energies[kElectrostaticIndex];
 802         Vector *forces = &arena->forces[kVDWIndex * arena->mol->natoms];
 803         Vector *eforces = &arena->forces[kElectrostaticIndex * arena->mol->natoms];
 804         s_calc_nonbonded_force_sub(arena, energies, eenergies, forces, eforces, NULL, NULL);
 805 }
 806
 807 static void
 808 s_calc_auxiliary_force(MDArena *arena)
 809 {
 810         Double *energies = &arena->energies[kAuxiliaryIndex];
 811         Vector *forces = &arena->forces[kAuxiliaryIndex * arena->mol->natoms];
 812         Atom *atoms = arena->mol->atoms;
 813         int natoms = arena->mol->natoms;
 814         int i, j;
 815
 816         /*  Centric force - OBSOLETE */
 817 /*      if (arena->centric_potential_force > 0) {
 818                 Vector center = arena->mol->atoms[arena->centric_potential_center].r;
 819                 Double rin = arena->centric_potential_inner_limit;
 820                 Double rout = arena->centric_potential_outer_limit;
 821                 Double k = arena->centric_potential_force;
 822                 for (i = 0; i < natoms; i++) {
 823                         Vector r21;
 824                         Double w1, w2, k0, k1;
 825                         VecSub(r21, atoms[i].r, center);
 826                         w2 = VecLength2(r21);
 827                         w1 = sqrt(w2);
 828                         if (w1 <= rin) {
 829                                 k0 = k1 = 0.0;
 830                         } else if (rout > 0 && w1 > rout) {
 831                                 k0 = k * (rout - rin) * (w1 + w1 - rout + rin);
 832                                 k1 = -2.0 * k * (rout - rin) / w1;
 833                         } else {
 834                                 k0 = k * (w1 - rin) * (w1 - rin);
 835                                 k1 = -2.0 * k * (1 - rin / w1);
 836                         }
 837                         *energies += k0;
 838                         VecScaleSelf(r21, k1);
 839                         VecInc(forces[i], r21);
 840 #if DEBUG
 841                         if (arena->debug_result && arena->debug_output_level > 1) {
 842                                 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);
 843                         }
 844 #endif
 845                 }
 846         }
 847         */
 848
 849         /*  Spherical boundary conditions  */
 850         if (arena->spherical_bc_force > 0) {
 851                 Vector center = arena->spherical_bc_center;
 852                 Double rin = arena->spherical_bc_inner_limit;
 853                 Double rout = arena->spherical_bc_outer_limit;
 854                 Double k = arena->spherical_bc_force;
 855                 for (i = 0; i < natoms; i++) {
 856                         Vector r21;
 857                         Double w1, w2, k0, k1;
 858                         VecSub(r21, atoms[i].r, center);
 859                         w2 = VecLength2(r21);
 860                         w1 = sqrt(w2);
 861                         if (w1 <= rin) {
 862                                 k0 = k1 = 0.0;
 863                         } else if (rout > 0 && w1 > rout) {
 864                                 k0 = k * (rout - rin) * (w1 + w1 - rout + rin);
 865                                 k1 = -2.0 * k * (rout - rin) / w1;
 866                         } else {
 867                                 k0 = k * (w1 - rin) * (w1 - rin);
 868                                 k1 = -2.0 * k * (1 - rin / w1);
 869                         }
 870                         *energies += k0;
 871                         VecScaleSelf(r21, k1);
 872                         VecInc(forces[i], r21);
 873                         if (arena->debug_result && arena->debug_output_level > 1) {
 874                                 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);
 875                         }
 876                 }
 877         }
 878
 879         /*  Box force  */
 880         if (arena->box_potential_force > 0) {
 881                 Double xsize = arena->box_potential_xsize;
 882                 Double ysize = arena->box_potential_ysize;
 883                 Double zsize = arena->box_potential_zsize;
 884                 Double k = arena->box_potential_force;
 885                 for (i = 0; i < natoms; i++) {
 886                         Vector r = atoms[i].r;
 887                         if (r.x > xsize)
 888                                 r.x -= xsize;
 889                         else if (r.x < -xsize)
 890                                 r.x += xsize;
 891                         else r.x = 0.0;
 892                         if (r.y > ysize)
 893                                 r.y -= ysize;
 894                         else if (r.y < -ysize)
 895                                 r.y += ysize;
 896                         else r.y = 0.0;
 897                         if (r.z > zsize)
 898                                 r.z -= zsize;
 899                         else if (r.z < -zsize)
 900                                 r.z += zsize;
 901                         else r.z = 0.0;
 902                         *energies += k * (r.x * r.x + r.y * r.y + r.z * r.z);
 903                         VecScaleInc(forces[i], r, -2);
 904                         if (arena->debug_result && arena->debug_output_level > 1) {
 905                                 fprintf(arena->debug_result, "auxiliary(box) force %d: {%f %f %f}\n", i, -2*r.x, -2*r.y, -2*r.z);
 906                         }
 907                 }
 908         }
 909
 910         /*  Fix atoms  */
 911         for (i = j = 0; i < natoms; i++) {
 912                 Atom *ap = &atoms[i];
 913                 Vector r21;
 914                 Double w1, w2, k0, k1;
 915                 if (ap->fix_force <= 0.0)
 916                         continue;
 917                 VecSub(r21, ap->r, ap->fix_pos);
 918                 w2 = VecLength2(r21);
 919                 w1 = sqrt(w2);
 920                 k0 = ap->fix_force * w2;
 921                 k1 = -2.0 * ap->fix_force * w1;
 922                 *energies += k0;
 923                 VecScaleSelf(r21, k1);
 924                 VecInc(forces[i], r21);
 925                 j++;
 926                 if (arena->debug_result && arena->debug_output_level > 1) {
 927                         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);
 928                 }
 929         }
 930
 931         /*  Graphite  */
 932         if (arena->graphite != NULL && arena->use_graphite) {
 933                 graphite_force(arena->graphite, arena, energies, forces);
 934         }
 935 }
 936
 937 static void
 938 s_md_debug_output_positions(MDArena *arena)
 939 {
 940         int i;
 941         if (arena->debug_result == NULL || arena->debug_output_level == 0)
 942                 return;
 943         fprintf(arena->debug_result, "\n  Atom positions, velocities, and forces at step %d\n", arena->step);
 944         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");
 945         for (i = 0; i < arena->mol->natoms; i++) {
 946                 Atom *ap = &arena->mol->atoms[i];
 947                 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);
 948         }
 949 }
 950
 951 void
 952 calc_force(MDArena *arena)
 953 {
 954         Int i, j, natoms;
 955         Vector *ff, *fa;
 956         Molecule *mol = arena->mol;
 957         Int doSurface = 0;
 958
 959         natoms = mol->natoms;
 960
 961         /*  Clear the energy/force storage  */
 962         for (i = 0; i < natoms; i++)
 963                 VecZero(mol->atoms[i].f);
 964
 965         memset(arena->energies, 0, sizeof(Double) * kSlowIndex);
 966         memset(arena->forces, 0, sizeof(Vector) * kSlowIndex * natoms);
 967         arena->total_energy = 0.0;
 968         arena->alchem_energy = 0.0;
 969
 970         if (arena->step == arena->start_step || arena->step % arena->surface_potential_freq == 0) {
 971                 doSurface = 1;
 972                 arena->energies[kSurfaceIndex] = 0.0;
 973                 memset(arena->forces + kSurfaceIndex * natoms, 0, sizeof(Vector) * natoms);
 974         }
 975
 976         s_calc_bond_force(arena);
 977         s_calc_angle_force(arena);
 978         s_calc_dihedral_force(arena);
 979         s_calc_improper_force(arena);
 980         s_calc_nonbonded_force(arena);
 981         s_calc_auxiliary_force(arena);
 982
 983         if (doSurface && arena->probe_radius > 0.0) {
 984         /*      if (arena->sphere_points == 0)
 985                         calc_surface_force(arena);
 986                 else
 987                         calc_surface_force_2(arena); */
 988                 calc_surface_force(arena);
 989         }
 990
 991         /*  Sum up all partial forces and energies  */
 992         arena->total_energy = 0.0;
 993         for (i = 0; i < kKineticIndex; i++)
 994                 arena->total_energy += arena->energies[i];
 995
 996         for (i = 0; i < natoms; i++) {
 997                 fa = &mol->atoms[i].f;
 998                 ff = &arena->forces[i];
 999                 for (j = 0; j < kKineticIndex; j++) {
1000                         VecInc(*fa, *ff);
1001                         ff += natoms;
1002                 }
1003         }
1004
1005         /*  The total (internal) force must be zero  */
1006 /*      {
1007                 Vector f;
1008                 Double w;
1009                 VecZero(f);
1010                 for (i = 0; i < natoms; i++) {
1011                         Vector *fp = &(mol->atoms[i].f);
1012                         VecInc(f, *fp);
1013                 }
1014                 w = (natoms > 0 ? 1.0 / natoms : 0.0);
1015                 VecScaleSelf(f, w);
1016                 for (i = 0; i < natoms; i++) {
1017                         Vector *fp = &(mol->atoms[i].f);
1018                         VecDec(*fp, f);
1019                 }
1020         } */
1021
1022         s_md_debug_output_positions(arena);
1023 }
1024
1025 void
1026 calc_pair_interaction(MDArena *arena, const MDGroupFlags *group_flags_1, const MDGroupFlags *group_flags_2)
1027 {
1028         Vector *forces, *eforces;
1029         if (arena->pair_forces != NULL) {
1030                 forces = arena->pair_forces;
1031                 eforces = forces + arena->mol->natoms;
1032                 memset(forces, 0, sizeof(Vector) * arena->mol->natoms * 2);
1033         } else forces = eforces = NULL;
1034         arena->pair_energies[0] = arena->pair_energies[1] = 0.0;
1035         s_calc_nonbonded_force_sub(arena, &arena->pair_energies[0], &arena->pair_energies[1], forces, eforces, group_flags_1, group_flags_2);
1036 }