Bug fixes to triclinic boxes for reference AmoebaMultipoleForce

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceMultipoleForce.cpp

 
-/* Portions copyright (c) 2006-2014 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2015 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -5036,7 +5036,7 @@ void AmoebaReferencePmeMultipoleForce::computeAmoebaBsplines(const vector<Multip
 
             RealOpenMM w  = position[0]*_recipBoxVectors[0][jj]+position[1]*_recipBoxVectors[1][jj]+position[2]*_recipBoxVectors[2][jj];
             RealOpenMM fr = _pmeGridDimensions[jj]*(w-(int)(w+0.5)+0.5);
-            int ifr       = static_cast<int>(fr);
+            int ifr       = static_cast<int>(floor(fr));
             w             = fr - ifr;
             igrid[jj]     = ifr - AMOEBA_PME_ORDER + 1;
             igrid[jj]    += igrid[jj] < 0 ? _pmeGridDimensions[jj] : 0;
@@ -5197,14 +5197,14 @@ void AmoebaReferencePmeMultipoleForce::transformPotentialToCartesianCoordinates(
     for (int i = 0; i < 3; i++) {
         for (int j = 0; j < 6; j++) {
             b[i][j] = a[index1[i]][index1[j]]*a[index2[i]][index2[j]];
-            if (index1[i] != index2[i])
+            if (index1[j] != index2[j])
                 b[i][j] *= 2;
         }
     }
     for (int i = 3; i < 6; i++) {
         for (int j = 0; j < 6; j++) {
             b[i][j] = a[index1[i]][index1[j]]*a[index2[i]][index2[j]];
-            if (index1[i] != index2[i])
+            if (index1[j] != index2[j])
                 b[i][j] += a[index1[i]][index2[j]]*a[index2[i]][index1[j]];
         }
     }
@@ -5410,7 +5410,7 @@ void AmoebaReferencePmeMultipoleForce::computeFixedPotentialFromGrid()
     }
 }
 
-t_complex AmoebaReferencePmeMultipoleForce::computeInducedDipoleGridValue(const int2& particleGridIndices, const RealVec* fracToCart, int ix, int iy,
+t_complex AmoebaReferencePmeMultipoleForce::computeInducedDipoleGridValue(const int2& particleGridIndices, const RealVec* cartToFrac, int ix, int iy,
                                                                           const IntVec& gridPoint,
                                                                           const vector<RealVec>& inputInducedDipole,
                                                                           const vector<RealVec>& inputInducedDipolePolar) const 
@@ -5430,16 +5430,12 @@ t_complex AmoebaReferencePmeMultipoleForce::computeInducedDipoleGridValue(const
         if (iz >= _pmeGridDimensions[2]) {
             iz -= _pmeGridDimensions[2];
         }
-        RealVec inducedDipole = RealVec(inputInducedDipole[atomIndex][0]*fracToCart[0][0] + inputInducedDipole[atomIndex][1]*fracToCart[0][1] + inputInducedDipole[atomIndex][2]*fracToCart[0][2],
-                                        inputInducedDipole[atomIndex][0]*fracToCart[1][0] + inputInducedDipole[atomIndex][1]*fracToCart[1][1] + inputInducedDipole[atomIndex][2]*fracToCart[1][2],
-                                        inputInducedDipole[atomIndex][0]*fracToCart[2][0] + inputInducedDipole[atomIndex][1]*fracToCart[2][1] + inputInducedDipole[atomIndex][2]*fracToCart[2][2]);
-        RealVec inducedDipolePolar = RealVec(inputInducedDipolePolar[atomIndex][0]*fracToCart[0][0] + inputInducedDipolePolar[atomIndex][1]*fracToCart[0][1] + inputInducedDipolePolar[atomIndex][2]*fracToCart[0][2],
-                                             inputInducedDipolePolar[atomIndex][0]*fracToCart[1][0] + inputInducedDipolePolar[atomIndex][1]*fracToCart[1][1] + inputInducedDipolePolar[atomIndex][2]*fracToCart[1][2],
-                                             inputInducedDipolePolar[atomIndex][0]*fracToCart[2][0] + inputInducedDipolePolar[atomIndex][1]*fracToCart[2][1] + inputInducedDipolePolar[atomIndex][2]*fracToCart[2][2]);
-
-//        RealVec inducedDipolePolar  = RealVec(scale[0]*inputInducedDipolePolar[atomIndex][0],
-//                                               scale[1]*inputInducedDipolePolar[atomIndex][1],
-//                                               scale[2]*inputInducedDipolePolar[atomIndex][2]);
+        RealVec inducedDipole = RealVec(inputInducedDipole[atomIndex][0]*cartToFrac[0][0] + inputInducedDipole[atomIndex][1]*cartToFrac[0][1] + inputInducedDipole[atomIndex][2]*cartToFrac[0][2],
+                                        inputInducedDipole[atomIndex][0]*cartToFrac[1][0] + inputInducedDipole[atomIndex][1]*cartToFrac[1][1] + inputInducedDipole[atomIndex][2]*cartToFrac[1][2],
+                                        inputInducedDipole[atomIndex][0]*cartToFrac[2][0] + inputInducedDipole[atomIndex][1]*cartToFrac[2][1] + inputInducedDipole[atomIndex][2]*cartToFrac[2][2]);
+        RealVec inducedDipolePolar = RealVec(inputInducedDipolePolar[atomIndex][0]*cartToFrac[0][0] + inputInducedDipolePolar[atomIndex][1]*cartToFrac[0][1] + inputInducedDipolePolar[atomIndex][2]*cartToFrac[0][2],
+                                             inputInducedDipolePolar[atomIndex][0]*cartToFrac[1][0] + inputInducedDipolePolar[atomIndex][1]*cartToFrac[1][1] + inputInducedDipolePolar[atomIndex][2]*cartToFrac[1][2],
+                                             inputInducedDipolePolar[atomIndex][0]*cartToFrac[2][0] + inputInducedDipolePolar[atomIndex][1]*cartToFrac[2][1] + inputInducedDipolePolar[atomIndex][2]*cartToFrac[2][2]);
 
         RealOpenMM4 t = _thetai[0][atomIndex*AMOEBA_PME_ORDER+ix];
         RealOpenMM4 u = _thetai[1][atomIndex*AMOEBA_PME_ORDER+iy];
@@ -5460,10 +5456,10 @@ t_complex AmoebaReferencePmeMultipoleForce::computeInducedDipoleGridValue(const
 void AmoebaReferencePmeMultipoleForce::spreadInducedDipolesOnGrid(const vector<RealVec>& inputInducedDipole,
                                                                   const vector<RealVec>& inputInducedDipolePolar)
 {
-    RealVec fracToCart[3];
+    RealVec cartToFrac[3];
     for (int i = 0; i < 3; i++)
         for (int j = 0; j < 3; j++)
-            fracToCart[i][j] = _pmeGridDimensions[j]*_recipBoxVectors[i][j];
+            cartToFrac[j][i] = _pmeGridDimensions[j]*_recipBoxVectors[i][j];
 
     for (int gridIndex = 0; gridIndex < _totalGridSize; gridIndex++)
     {
@@ -5486,13 +5482,13 @@ void AmoebaReferencePmeMultipoleForce::spreadInducedDipolesOnGrid(const vector<R
                 int2 particleGridIndices;
                 particleGridIndices[0] = x*_pmeGridDimensions[1]*_pmeGridDimensions[2]+y*_pmeGridDimensions[2]+z1;
                 particleGridIndices[1] = x*_pmeGridDimensions[1]*_pmeGridDimensions[2]+y*_pmeGridDimensions[2]+z2;
-                gridValue             += computeInducedDipoleGridValue(particleGridIndices, fracToCart, ix, iy, gridPoint, inputInducedDipole, inputInducedDipolePolar);
+                gridValue             += computeInducedDipoleGridValue(particleGridIndices, cartToFrac, ix, iy, gridPoint, inputInducedDipole, inputInducedDipolePolar);
 
                 if (z1 > gridPoint[2])
                 {
                     particleGridIndices[0]  =  x*_pmeGridDimensions[1]*_pmeGridDimensions[2]+y*_pmeGridDimensions[2];
                     particleGridIndices[1]  =  x*_pmeGridDimensions[1]*_pmeGridDimensions[2]+y*_pmeGridDimensions[2]+gridPoint[2];
-                    gridValue              +=  computeInducedDipoleGridValue(particleGridIndices, fracToCart, ix, iy, gridPoint, inputInducedDipole, inputInducedDipolePolar);
+                    gridValue              +=  computeInducedDipoleGridValue(particleGridIndices, cartToFrac, ix, iy, gridPoint, inputInducedDipole, inputInducedDipolePolar);
                 }
             }
         }
@@ -5714,8 +5710,6 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::computeReciprocalSpaceFixedMultipol
     const int deriv3[] = {3, 8, 9, 6, 14, 16, 12, 19, 17, 18};
     vector<RealOpenMM> cphi(10*_numParticles);
     transformPotentialToCartesianCoordinates(_phi, cphi);
-//    RealVec scale;
-//    getPmeScale(scale);
     RealVec fracToCart[3];
     for (int i = 0; i < 3; i++)
         for (int j = 0; j < 3; j++)
@@ -5739,22 +5733,6 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::computeReciprocalSpaceFixedMultipol
         multipole[8] = particleData[i].quadrupole[QXZ]*2.0;
         multipole[9] = particleData[i].quadrupole[QYZ]*2.0;
 
-//        const RealOpenMM* phi = &_phi[20*i];
-//        torques[i][0] += _electric*(multipole[3]*scale[1]*phi[2] - multipole[2]*scale[2]*phi[3]
-//                      + 2.0*(multipole[6]-multipole[5])*scale[1]*scale[2]*phi[9]
-//                      + multipole[8]*scale[0]*scale[1]*phi[7] + multipole[9]*scale[1]*scale[1]*phi[5]
-//                      - multipole[7]*scale[0]*scale[2]*phi[8] - multipole[9]*scale[2]*scale[2]*phi[6]);
-//
-//        torques[i][1] += _electric*(multipole[1]*scale[2]*phi[3] - multipole[3]*scale[0]*phi[1]
-//                      + 2.0*(multipole[4]-multipole[6])*scale[0]*scale[2]*phi[8]
-//                      + multipole[7]*scale[1]*scale[2]*phi[9] + multipole[8]*scale[2]*scale[2]*phi[6]
-//                      - multipole[8]*scale[0]*scale[0]*phi[4] - multipole[9]*scale[0]*scale[1]*phi[7]);
-//
-//        torques[i][2] += _electric*(multipole[2]*scale[0]*phi[1] - multipole[1]*scale[1]*phi[2]
-//                      + 2.0*(multipole[5]-multipole[4])*scale[0]*scale[1]*phi[7]
-//                      + multipole[7]*scale[0]*scale[0]*phi[4] + multipole[9]*scale[0]*scale[2]*phi[8]
-//                      - multipole[7]*scale[1]*scale[1]*phi[5] - multipole[8]*scale[1]*scale[2]*phi[9]);
-
         const RealOpenMM* phi = &cphi[10*i];
         torques[i][0] += _electric*(multipole[3]*phi[2] - multipole[2]*phi[3]
                       + 2.0*(multipole[6]-multipole[5])*phi[9]
@@ -5782,15 +5760,6 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::computeReciprocalSpaceFixedMultipol
         multipole[7] = _transformed[i].quadrupole[QXY];
         multipole[8] = _transformed[i].quadrupole[QXZ];
         multipole[9] = _transformed[i].quadrupole[QYZ];
-//        multipole[1] *= scale[0];
-//        multipole[2] *= scale[1];
-//        multipole[3] *= scale[2];
-//        multipole[4] *= scale[0]*scale[0];
-//        multipole[5] *= scale[1]*scale[1];
-//        multipole[6] *= scale[2]*scale[2];
-//        multipole[7] *= scale[0]*scale[1];
-//        multipole[8] *= scale[0]*scale[2];
-//        multipole[9] *= scale[1]*scale[2];
 
         RealVec f = RealVec(0.0, 0.0, 0.0);
         for (int k = 0; k < 10; k++) {
@@ -5799,9 +5768,6 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::computeReciprocalSpaceFixedMultipol
             f[1]   += multipole[k]*_phi[20*i+deriv2[k]];
             f[2]   += multipole[k]*_phi[20*i+deriv3[k]];
         }
-//        f[0]           *= scale[0];
-//        f[1]           *= scale[1];
-//        f[2]           *= scale[2];
         f              *= (_electric);
         forces[i]      -= RealVec(f[0]*fracToCart[0][0] + f[1]*fracToCart[0][1] + f[2]*fracToCart[0][2],
                                   f[0]*fracToCart[1][0] + f[1]*fracToCart[1][1] + f[2]*fracToCart[1][2],
@@ -5832,8 +5798,6 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::computeReciprocalSpaceInducedDipole
     for (int i = 0; i < 3; i++)
         for (int j = 0; j < 3; j++)
             cartToFrac[j][i] = fracToCart[i][j] = _pmeGridDimensions[j]*_recipBoxVectors[i][j];
-//    RealVec scale;
-//    getPmeScale(scale);
     RealOpenMM energy = 0.0;
     for (int i = 0; i < _numParticles; i++) {
 
@@ -5881,16 +5845,6 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::computeReciprocalSpaceInducedDipole
         multipole[7] = _transformed[i].quadrupole[QXY];
         multipole[8] = _transformed[i].quadrupole[QXZ];
         multipole[9] = _transformed[i].quadrupole[QYZ];
-//        multipole[1] *= scale[0];
-//        multipole[2] *= scale[1];
-//        multipole[3] *= scale[2];
-//
-//        multipole[4] *= scale[0]*scale[0];
-//        multipole[5] *= scale[1]*scale[1];
-//        multipole[6] *= scale[2]*scale[2];
-//        multipole[7] *= scale[0]*scale[1];
-//        multipole[8] *= scale[0]*scale[2];
-//        multipole[9] *= scale[1]*scale[2];
 
         inducedDipole[0] = _inducedDipole[i][0]*cartToFrac[0][0] + _inducedDipole[i][1]*cartToFrac[0][1] + _inducedDipole[i][2]*cartToFrac[0][2];
         inducedDipole[1] = _inducedDipole[i][0]*cartToFrac[1][0] + _inducedDipole[i][1]*cartToFrac[1][1] + _inducedDipole[i][2]*cartToFrac[1][2];
@@ -5930,9 +5884,6 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::computeReciprocalSpaceInducedDipole
             f[2] += multipole[k]*_phidp[20*i+deriv3[k]];
         }
 
-//        f[0]           *= scale[0];
-//        f[1]           *= scale[1];
-//        f[2]           *= scale[2];
         f              *= (0.5*_electric);
         forces[iIndex] -= RealVec(f[0]*fracToCart[0][0] + f[1]*fracToCart[0][1] + f[2]*fracToCart[0][2],
                                   f[0]*fracToCart[1][0] + f[1]*fracToCart[1][1] + f[2]*fracToCart[1][2],
@@ -5943,9 +5894,6 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::computeReciprocalSpaceInducedDipole
 
 void AmoebaReferencePmeMultipoleForce::recordFixedMultipoleField()
 {
-//    RealVec scale;
-//    getPmeScale(scale);
-//    scale *= -1.0;
     RealVec fracToCart[3];
     for (int i = 0; i < 3; i++)
         for (int j = 0; j < 3; j++)
@@ -5968,9 +5916,6 @@ void AmoebaReferencePmeMultipoleForce::initializeInducedDipoles(vector<UpdateInd
 
 void AmoebaReferencePmeMultipoleForce::recordInducedDipoleField(vector<RealVec>& field, vector<RealVec>& fieldPolar)
 {
-//    RealVec scale;
-//    getPmeScale(scale);
-//    scale *= -1.0;
     RealVec fracToCart[3];
     for (int i = 0; i < 3; i++)
         for (int j = 0; j < 3; j++)
@@ -6152,7 +6097,6 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::calculatePmeSelfEnergy(const vector
     RealOpenMM term      = 2.0*_alphaEwald*_alphaEwald;
     RealOpenMM energy    = (cii + term*(dii/3.0 + 2.0*term*qii/5.0));
                energy   *= -(_electric*_alphaEwald/(_dielectric*SQRT_PI));
-
     return energy;
 }
 
@@ -6898,7 +6842,7 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::calculateElectrostatic(const vector
         }
     }
 
-    calculatePmeSelfTorque(particleData, torques); 
+    calculatePmeSelfTorque(particleData, torques);
     energy += computeReciprocalSpaceInducedDipoleForceAndEnergy(getPolarizationType(), particleData, forces, torques);
     energy += computeReciprocalSpaceFixedMultipoleForceAndEnergy(particleData, forces, torques);
     energy += calculatePmeSelfEnergy(particleData);

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceMultipoleForce.h

-/* Portions copyright (c) 2006-2014 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2015 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -973,7 +973,7 @@ protected:
      * @param particleK         positions and parameters (charge, labFrame dipoles, quadrupoles, ...) for particle K
      * @param scalingFactors    scaling factors for interaction
      * @param forces            vector of particle forces to be updated
-     * @param torques           vector of particle torques to be updated
+     * @param torque            vector of particle torques to be updated
      */
     RealOpenMM calculateElectrostaticPairIxn(const MultipoleParticleData& particleI, const MultipoleParticleData& particleK,
                                              const std::vector<RealOpenMM>& scalingFactors, std::vector<OpenMM::RealVec>& forces, std::vector<RealVec>& torque) const;
@@ -987,6 +987,7 @@ protected:
      * @param particleW               particle3 of lab frame for particleI 
      * @param axisType                axis type (Bisector/Z-then-X, ...)
      * @param torque                  torque on particle I
+     * @param forces                  vector of particle forces to be updated
      */
     void mapTorqueToForceForParticle(const MultipoleParticleData& particleI,
                                      const MultipoleParticleData& particleU,
@@ -1004,8 +1005,6 @@ protected:
      * @param axisType                vector of axis types (Bisector/Z-then-X, ...) for particles
      * @param torques                 output torques
      * @param forces                  output forces 
-     *
-     * @return energy
      */
     void mapTorqueToForce(std::vector<MultipoleParticleData>& particleData, 
                           const std::vector<int>& multipoleAtomXs,
@@ -1388,7 +1387,7 @@ private:
     /**
      * Modify input vector of differences in particle positions for periodic boundary conditions.
      * 
-     * @param delta                   input vector of difference in particle positios; on output adjusted for
+     * @param delta                   input vector of difference in particle positions; on output adjusted for
      *                                periodic boundary conditions
      */
     void getPeriodicDelta(RealVec& deltaR) const;
@@ -1466,12 +1465,9 @@ private:
      * Compute induced dipole grid value.
      *
      * @param particleGridIndices     particle grid indices
-     * @param scale                   integer grid dimension/box size for each dimension
      * @param ix                      x-dimension offset value
      * @param iy                      y-dimension offset value
      * @param gridPoint               grid point for which value is to be computed
-     * @param inputInducedDipole      induced dipole value
-     * @param inputInducedDipolePolar induced dipole value
      */
      RealOpenMM computeFixedMultipolesGridValue(const int2& particleGridIndices, int ix, int iy, const IntVec& gridPoint) const;
 
@@ -1542,7 +1538,7 @@ private:
      * @param jIndex        particle J index
      * @param preFactor1    first factor used in calculating field
      * @param preFactor2    second factor used in calculating field
-     * @param deltaR        delta in particle positions after adjusting for periodic boundary conditions
+     * @param delta         delta in particle positions after adjusting for periodic boundary conditions
      * @param inducedDipole vector of induced dipoles
      * @param field         vector of field at each particle due induced dipole of other particles
      */
@@ -1574,14 +1570,14 @@ private:
      * Compute induced dipole grid value.
      *
      * @param atomIndices             indices of first and last atom contiputing to grid point value
-     * @param scale                   integer grid dimension/box size for each dimension
+     * @param cartToFrac              transformation matrix from Cartesian to fractional coordinates
      * @param ix                      x-dimension offset value
      * @param iy                      y-dimension offset value
      * @param gridPoint               grid point for which value is to be computed
      * @param inputInducedDipole      induced dipole value
      * @param inputInducedDipolePolar induced dipole polar value
      */
-    t_complex computeInducedDipoleGridValue(const int2& atomIndices, const RealVec* fracToCart, int ix, int iy, const IntVec& gridPoint,
+    t_complex computeInducedDipoleGridValue(const int2& atomIndices, const RealVec* cartToFrac, int ix, int iy, const IntVec& gridPoint,
                                             const std::vector<RealVec>& inputInducedDipole,
                                             const std::vector<RealVec>& inputInducedDipolePolar) const;
 

plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaMultipoleForce.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -39,6 +39,7 @@
 #include "openmm/System.h"
 #include "openmm/AmoebaMultipoleForce.h"
 #include "openmm/LangevinIntegrator.h"
+#include "openmm/Vec3.h"
 #include <iostream>
 #include <vector>
 #include <stdlib.h>
@@ -50,6 +51,7 @@
 
 
 using namespace OpenMM;
+using namespace std;
 
 extern "C" OPENMM_EXPORT void registerAmoebaReferenceKernelFactories();
 
@@ -2802,6 +2804,133 @@ static void testMultipoleGridPotential( FILE* log ) {
 
 }
 
+void testTriclinic() {
+    // Create a triclinic box containing eight particles.
+
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(1.8643, 0, 0), Vec3(-0.16248445120445926, 1.8572057756524414, 0), Vec3(0.16248445120445906, -0.14832299817478897, 1.8512735025730875));
+    for (int i = 0; i < 24; i++)
+        system.addParticle(1.0);
+    AmoebaMultipoleForce* force = new AmoebaMultipoleForce();
+    system.addForce(force);
+    force->setNonbondedMethod(AmoebaMultipoleForce::PME);
+    force->setCutoffDistance(0.7);
+    force->setMutualInducedTargetEpsilon(1e-6);
+    vector<int> grid(3);
+    grid[0] = grid[1] = grid[2] = 24;
+    force->setPmeGridDimensions(grid);
+    force->setAEwald(5.4459051633620055);
+    double o_charge = -0.42616, h_charge = 0.21308;
+    vector<double> o_dipole(3), h_dipole(3), o_quadrupole(9, 0.0), h_quadrupole(9, 0.0);
+    o_dipole[0] = 0;
+    o_dipole[1] = 0;
+    o_dipole[2] = 0.0033078867454609203;
+    h_dipole[0] = -0.0053536858428776405;
+    h_dipole[1] = 0;
+    h_dipole[2] = -0.014378273997907321;
+
+    o_quadrupole[0] = 0.00016405937591036892;
+    o_quadrupole[4] = -0.00021618201787005826;
+    o_quadrupole[8] = 5.212264195968935e-05;
+    h_quadrupole[0] = 0.00011465301060008312;
+    h_quadrupole[4] = 8.354184196619263e-05;
+    h_quadrupole[8] = -0.00019819485256627578;
+    h_quadrupole[2] = h_quadrupole[6] = -6.523731100577879e-05;
+
+    for (int i = 0; i < 8; i++) {
+        int atom1 = 3*i, atom2 = 3*i+1, atom3 = 3*i+2;
+        force->addMultipole(o_charge, o_dipole, o_quadrupole, 1, atom2, atom3, -1, 0.39, pow(0.001*0.92, 1.0/6.0), 0.001*0.92);
+        force->addMultipole(h_charge, h_dipole, h_quadrupole, 0, atom1, atom3, -1, 0.39, pow(0.001*0.539, 1.0/6.0), 0.001*0.539);
+        force->addMultipole(h_charge, h_dipole, h_quadrupole, 0, atom1, atom2, -1, 0.39, pow(0.001*0.539, 1.0/6.0), 0.001*0.539);
+        vector<int> coval1_12(2);
+        coval1_12[0] = atom2;
+        coval1_12[1] = atom3;
+        force->setCovalentMap(atom1, AmoebaMultipoleForce::Covalent12, coval1_12);
+        vector<int> coval2_12(1);
+        coval2_12[0] = atom1;
+        force->setCovalentMap(atom2, AmoebaMultipoleForce::Covalent12, coval2_12);
+        force->setCovalentMap(atom3, AmoebaMultipoleForce::Covalent12, coval2_12);
+        vector<int> coval2_13(1);
+        coval2_13[0] = atom3;
+        force->setCovalentMap(atom2, AmoebaMultipoleForce::Covalent13, coval2_13);
+        vector<int> coval3_13(1);
+        coval3_13[0] = atom2;
+        force->setCovalentMap(atom3, AmoebaMultipoleForce::Covalent13, coval3_13);
+        vector<int> polar(3);
+        polar[0] = atom1;
+        polar[1] = atom2;
+        polar[2] = atom3;
+        force->setCovalentMap(atom1, AmoebaMultipoleForce::PolarizationCovalent11, polar);
+        force->setCovalentMap(atom2, AmoebaMultipoleForce::PolarizationCovalent11, polar);
+        force->setCovalentMap(atom3, AmoebaMultipoleForce::PolarizationCovalent11, polar);
+    }
+    vector<Vec3> positions(24);
+    positions[0] = Vec3(0.867966, 0.708769, -0.0696862);
+    positions[1] = Vec3(0.780946, 0.675579, -0.0382259);
+    positions[2] = Vec3(0.872223, 0.681424, -0.161756);
+    positions[3] = Vec3(-0.0117313, 0.824445, 0.683762);
+    positions[4] = Vec3(0.0216892, 0.789544, 0.605003);
+    positions[5] = Vec3(0.0444268, 0.782601, 0.75302);
+    positions[6] = Vec3(0.837906, -0.0092611, 0.681463);
+    positions[7] = Vec3(0.934042, 0.0098069, 0.673406);
+    positions[8] = Vec3(0.793962, 0.0573676, 0.626984);
+    positions[9] = Vec3(0.658995, 0.184432, -0.692317);
+    positions[10] = Vec3(0.588543, 0.240231, -0.671793);
+    positions[11] = Vec3(0.618153, 0.106275, -0.727368);
+    positions[12] = Vec3(0.71466, 0.575358, 0.233152);
+    positions[13] = Vec3(0.636812, 0.612604, 0.286268);
+    positions[14] = Vec3(0.702502, 0.629465, 0.15182);
+    positions[15] = Vec3(-0.242658, -0.850419, -0.250483);
+    positions[16] = Vec3(-0.169206, -0.836825, -0.305829);
+    positions[17] = Vec3(-0.279321, -0.760247, -0.24031);
+    positions[18] = Vec3(-0.803838, -0.360559, 0.230369);
+    positions[19] = Vec3(-0.811375, -0.424813, 0.301849);
+    positions[20] = Vec3(-0.761939, -0.2863, 0.270962);
+    positions[21] = Vec3(-0.148063, 0.824409, -0.827221);
+    positions[22] = Vec3(-0.20902, 0.868798, -0.7677);
+    positions[23] = Vec3(-0.0700878, 0.882333, -0.832221);
+
+    // Compute the forces and energy.
+
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+    Context context(system, integrator, Platform::getPlatformByName("Reference"));
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+
+    // Compare them to values computed by Gromacs.
+
+    double expectedEnergy = 6.8278696;
+    vector<Vec3> expectedForce(24);
+    expectedForce[0] = Vec3(-104.755, 14.0833, 34.359);
+    expectedForce[1] = Vec3(97.324, -1.41419, -142.976);
+    expectedForce[2] = Vec3(29.3968, -6.31784, -3.8702);
+    expectedForce[3] = Vec3(39.1915, 66.852, -28.5767);
+    expectedForce[4] = Vec3(-8.17554, -6.71532, 7.63162);
+    expectedForce[5] = Vec3(-87.3745, -91.8639, 35.4761);
+    expectedForce[6] = Vec3(-19.7568, -40.8693, 5.60238);
+    expectedForce[7] = Vec3(0.840984, 26.878, 10.7822);
+    expectedForce[8] = Vec3(14.3469, 12.0583, -12.075);
+    expectedForce[9] = Vec3(13.757, 16.9954, 46.9403);
+    expectedForce[10] = Vec3(-5.04172, -14.008, -15.3804);
+    expectedForce[11] = Vec3(-2.1715, 3.7405, -37.2209);
+    expectedForce[12] = Vec3(-70.2284, -9.10438, -40.1287);
+    expectedForce[13] = Vec3(4.46014, 3.89949, 4.64842);
+    expectedForce[14] = Vec3(43.045, -4.79905, 151.879);
+    expectedForce[15] = Vec3(20.2129, -0.895376, -27.2086);
+    expectedForce[16] = Vec3(-5.10448, 3.57732, 17.0498);
+    expectedForce[17] = Vec3(-13.7695, -1.03345, 12.3093);
+    expectedForce[18] = Vec3(2.94972, 0.338904, -10.9914);
+    expectedForce[19] = Vec3(0.69036, 1.22591, 4.50198);
+    expectedForce[20] = Vec3(-4.61495, -2.76981, 3.57732);
+    expectedForce[21] = Vec3(73.1489, 16.167, -99.5834);
+    expectedForce[22] = Vec3(-31.8235, 6.11282, -21.125);
+    expectedForce[23] = Vec3(13.167, 7.42242, 103.102);
+    for (int i = 0; i < 24; i++) {
+        ASSERT_EQUAL_VEC(expectedForce[i], state.getForces()[i], 1e-2);
+    }
+    ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-2);
+}
+
 int main( int numberOfArguments, char* argv[] ) {
 
     try {
@@ -2849,6 +2978,10 @@ int main( int numberOfArguments, char* argv[] ) {
         // large box of water
 
         testPMEMutualPolarizationLargeWater( log );
+        
+        // triclinic box of water
+        
+        testTriclinic();
 
     } catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;