Merge branch 'master' of github.com:leeping/openmm

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

@@ -1563,6 +1563,28 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateForceAndEnergy( const std::ve
     return energy;
 }
 
+void AmoebaReferenceMultipoleForce::calculateInducedDipoles(const std::vector<RealVec>& particlePositions,
+                                                            const std::vector<RealOpenMM>& charges,
+                                                            const std::vector<RealOpenMM>& dipoles,
+                                                            const std::vector<RealOpenMM>& quadrupoles,
+                                                            const std::vector<RealOpenMM>& tholes,
+                                                            const std::vector<RealOpenMM>& dampingFactors,
+                                                            const std::vector<RealOpenMM>& polarity,
+                                                            const std::vector<int>& axisTypes,
+                                                            const std::vector<int>& multipoleAtomZs,
+                                                            const std::vector<int>& multipoleAtomXs,
+                                                            const std::vector<int>& multipoleAtomYs,
+                                                            const std::vector< std::vector< std::vector<int> > >& multipoleAtomCovalentInfo,
+                                                            std::vector<RealVec>& outputInducedDipoles) {
+    // setup, including calculating induced dipoles
+
+    std::vector<MultipoleParticleData> particleData;
+    setup( particlePositions, charges, dipoles, quadrupoles, tholes,
+           dampingFactors, polarity, axisTypes, multipoleAtomZs, multipoleAtomXs, multipoleAtomYs,
+           multipoleAtomCovalentInfo, particleData );
+    outputInducedDipoles = _inducedDipole;
+}
+
 void AmoebaReferenceMultipoleForce::calculateAmoebaSystemMultipoleMoments( const std::vector<RealOpenMM>& masses,
                                                                            const std::vector<RealVec>& particlePositions,
                                                                            const std::vector<RealOpenMM>& charges,

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

@@ -471,6 +471,38 @@ public:
                                         const std::vector< std::vector< std::vector<int> > >& multipoleAtomCovalentInfo,
                                         std::vector<OpenMM::RealVec>& forces );
 
+    /**
+     * Calculate particle induced dipoles.
+     *
+     * @param masses                    particle masses
+     * @param particlePositions         Cartesian coordinates of particles
+     * @param charges                   scalar charges for each particle
+     * @param dipoles                   molecular frame dipoles for each particle
+     * @param quadrupoles               molecular frame quadrupoles for each particle
+     * @param tholes                    Thole factors for each particle
+     * @param dampingFactors            dampling factors for each particle
+     * @param polarity                  polarity for each particle
+     * @param axisTypes                 axis type (Z-then-X, ... ) for each particle
+     * @param multipoleAtomZs           indicies of particle specifying the molecular frame z-axis for each particle
+     * @param multipoleAtomXs           indicies of particle specifying the molecular frame x-axis for each particle
+     * @param multipoleAtomYs           indicies of particle specifying the molecular frame y-axis for each particle
+     * @param multipoleAtomCovalentInfo covalent info needed to set scaling factors
+     * @param outputMultipoleMoments    output multipole moments
+     */
+    void calculateInducedDipoles(const std::vector<OpenMM::RealVec>& particlePositions,
+                                 const std::vector<RealOpenMM>& charges,
+                                 const std::vector<RealOpenMM>& dipoles,
+                                 const std::vector<RealOpenMM>& quadrupoles,
+                                 const std::vector<RealOpenMM>& tholes,
+                                 const std::vector<RealOpenMM>& dampingFactors,
+                                 const std::vector<RealOpenMM>& polarity,
+                                 const std::vector<int>& axisTypes,
+                                 const std::vector<int>& multipoleAtomZs,
+                                 const std::vector<int>& multipoleAtomXs,
+                                 const std::vector<int>& multipoleAtomYs,
+                                 const std::vector< std::vector< std::vector<int> > >& multipoleAtomCovalentInfo,
+                                 std::vector<RealVec>& outputInducedDipoles);
+
     /**
      * Calculate system multipole moments.
      *

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

@@ -2605,6 +2605,40 @@ static void testPMEMutualPolarizationLargeWater( FILE* log ) {
 
 }
 
+// test querying particle induced dipoles
+
+static void testParticleInducedDipoles() {
+    int numberOfParticles     = 8;
+    int inputPmeGridDimension = 0;
+    double cutoff             = 9000000.0;
+    std::vector<Vec3> forces;
+    double energy;
+
+    System system;
+    AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();;
+    setupMultipoleAmmonia(system, amoebaMultipoleForce, AmoebaMultipoleForce::NoCutoff, AmoebaMultipoleForce::Mutual, 
+                                             cutoff, inputPmeGridDimension);
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+    Context context(system, integrator, Platform::getPlatformByName("Reference"));
+    getForcesEnergyMultipoleAmmonia(context, forces, energy);
+    std::vector<Vec3> dipole;
+    amoebaMultipoleForce->getInducedDipoles(context, dipole);
+    
+    // Compare to values calculated by TINKER.
+    
+    std::vector<Vec3> expectedDipole(numberOfParticles);
+    expectedDipole[0] = Vec3(0.0031710288, 9.3687453e-7, -0.0006919963);
+    expectedDipole[1] = Vec3(8.0279737504e-5, -0.000279376, 4.778060103e-5);
+    expectedDipole[2] = Vec3(0.000079322, 0.0002789804, 4.8696656126e-5);
+    expectedDipole[3] = Vec3(-0.0001407394, 1.540638116e-6, -0.0007077775);
+    expectedDipole[4] = Vec3(0.0019564439, -1.0409717e-7, 0.0007332188);
+    expectedDipole[5] = Vec3(0.0008213891, -0.0007749618, -0.0003883865);
+    expectedDipole[6] = Vec3(0.0046133992, -7.2868019e-7, 0.0002500622);
+    expectedDipole[7] = Vec3(0.0008204731, 0.0007772727, -0.0003856176);
+    for (int i = 0; i < numberOfParticles; i++)
+        ASSERT_EQUAL_VEC(expectedDipole[i], dipole[i], 1e-4);
+}
+
 // test computation of system multipole moments
 
 static void testSystemMultipoleMoments( FILE* log ) {
@@ -2778,6 +2812,10 @@ int main( int numberOfArguments, char* argv[] ) {
 
         testMultipoleAmmoniaDirectPolarization( log );
 
+        // test querying induced dipoles
+        
+        testParticleInducedDipoles();
+
         // test mutual polarization, no cutoff
 
         testMultipoleAmmoniaMutualPolarization( log );

plugins/cpupme/src/CpuPmeKernels.cpp

@@ -35,9 +35,9 @@
 #include "CpuPmeKernels.h"
 #include "SimTKOpenMMRealType.h"
 #include "openmm/internal/hardware.h"
+#include "openmm/internal/vectorize.h"
 #include <cmath>
 #include <cstring>
-#include <smmintrin.h>
 
 using namespace OpenMM;
 using namespace std;
@@ -47,52 +47,50 @@ static const int PME_ORDER = 5;
 bool CpuCalcPmeReciprocalForceKernel::hasInitializedThreads = false;
 int CpuCalcPmeReciprocalForceKernel::numThreads = 0;
 
-#define EXTRACT_FLOAT(v, element) _mm_cvtss_f32(_mm_shuffle_ps(v, v, _MM_SHUFFLE(0, 0, 0, element)))
-
 static void spreadCharge(int start, int end, float* posq, float* grid, int gridx, int gridy, int gridz, int numParticles, Vec3 periodicBoxSize) {
     float temp[4];
-    __m128 boxSize = _mm_set_ps(0, (float) periodicBoxSize[2], (float) periodicBoxSize[1], (float) periodicBoxSize[0]);
-    __m128 invBoxSize = _mm_set_ps(0, (float) (1/periodicBoxSize[2]), (float) (1/periodicBoxSize[1]), (float) (1/periodicBoxSize[0]));
-    __m128 gridSize = _mm_set_ps(0, gridz, gridy, gridx);
-    __m128i gridSizeInt = _mm_set_epi32(0, gridz, gridy, gridx);
-    __m128 one  = _mm_set1_ps(1);
-    __m128 scale = _mm_set1_ps(1.0f/(PME_ORDER-1));
+    fvec4 boxSize((float) periodicBoxSize[0], (float) periodicBoxSize[1], (float) periodicBoxSize[2], 0);
+    fvec4 invBoxSize((float) (1/periodicBoxSize[0]), (float) (1/periodicBoxSize[1]), (float) (1/periodicBoxSize[2]), 0);
+    fvec4 gridSize(gridx, gridy, gridz, 0);
+    ivec4 gridSizeInt(gridx, gridy, gridz, 0);
+    fvec4 one(1);
+    fvec4 scale(1.0f/(PME_ORDER-1));
     const float epsilonFactor = sqrt(ONE_4PI_EPS0);
     memset(grid, 0, sizeof(float)*gridx*gridy*gridz);
     for (int i = start; i < end; i++) {
         // Find the position relative to the nearest grid point.
         
-        __m128 pos = _mm_loadu_ps(&posq[4*i]);
-        __m128 posFloor = _mm_floor_ps(_mm_mul_ps(pos, invBoxSize));
-        __m128 posInBox = _mm_sub_ps(pos, _mm_mul_ps(boxSize, posFloor));
-        __m128 t = _mm_mul_ps(_mm_mul_ps(posInBox, invBoxSize), gridSize);
-        __m128i ti = _mm_cvttps_epi32(t);
-        __m128 dr = _mm_sub_ps(t, _mm_cvtepi32_ps(ti));
-        __m128i gridIndex = _mm_sub_epi32(ti, _mm_and_si128(gridSizeInt, _mm_cmpeq_epi32(ti, gridSizeInt)));
+        fvec4 pos(&posq[4*i]);
+        fvec4 posFloor = floor(pos*invBoxSize);
+        fvec4 posInBox = pos-boxSize*posFloor;
+        fvec4 t = posInBox*invBoxSize*gridSize;
+        ivec4 ti = t;
+        fvec4 dr = t-ti;
+        ivec4 gridIndex = ti-(gridSizeInt&ti==gridSizeInt);
         
         // Compute the B-spline coefficients.
         
-        __m128 data[PME_ORDER];
-        data[PME_ORDER-1] = _mm_setzero_ps();
+        fvec4 data[PME_ORDER];
+        data[PME_ORDER-1] = 0.0f;
         data[1] = dr;
-        data[0] = _mm_sub_ps(one, dr);
+        data[0] = one-dr;
         for (int j = 3; j < PME_ORDER; j++) {
-            __m128 div = _mm_set1_ps(1.0f/(j-1));
-            data[j-1] = _mm_mul_ps(_mm_mul_ps(div, dr), data[j-2]);
+            fvec4 div(1.0f/(j-1));
+            data[j-1] = div*dr*data[j-2];
             for (int k = 1; k < j-1; k++)
-                data[j-k-1] = _mm_mul_ps(div, _mm_add_ps(_mm_mul_ps(_mm_add_ps(dr, _mm_set1_ps(k)), data[j-k-2]), _mm_mul_ps(_mm_sub_ps(_mm_set1_ps(j-k), dr), data[j-k-1])));
-            data[0] = _mm_mul_ps(_mm_mul_ps(div, _mm_sub_ps(one, dr)), data[0]);
+                data[j-k-1] = div*((dr+k)*data[j-k-2]+(fvec4(j-k)-dr)*data[j-k-1]);
+            data[0] = div*(one-dr)*data[0];
         }
-        data[PME_ORDER-1] = _mm_mul_ps(_mm_mul_ps(scale, dr), data[PME_ORDER-2]);
+        data[PME_ORDER-1] = scale*dr*data[PME_ORDER-2];
         for (int j = 1; j < (PME_ORDER-1); j++)
-            data[PME_ORDER-j-1] = _mm_mul_ps(scale, _mm_add_ps(_mm_mul_ps(_mm_add_ps(dr, _mm_set1_ps(j)), data[PME_ORDER-j-2]), _mm_mul_ps(_mm_sub_ps(_mm_set1_ps(PME_ORDER-j), dr), data[PME_ORDER-j-1])));
-        data[0] = _mm_mul_ps(_mm_mul_ps(scale, _mm_sub_ps(one, dr)), data[0]);
+            data[PME_ORDER-j-1] = scale*((dr+j)*data[PME_ORDER-j-2]+(fvec4(PME_ORDER-j)-dr)*data[PME_ORDER-j-1]);
+        data[0] = scale*(one-dr)*data[0];
         
         // Spread the charges.
         
-        int gridIndexX = _mm_extract_epi32(gridIndex, 0);
-        int gridIndexY = _mm_extract_epi32(gridIndex, 1);
-        int gridIndexZ = _mm_extract_epi32(gridIndex, 2);
+        int gridIndexX = gridIndex[0];
+        int gridIndexY = gridIndex[1];
+        int gridIndexZ = gridIndex[2];
         if (gridIndexX < 0)
             return; // This happens when a simulation blows up and coordinates become NaN.
         int zindex[PME_ORDER];
@@ -101,21 +99,21 @@ static void spreadCharge(int start, int end, float* posq, float* grid, int gridx
             zindex[j] -= (zindex[j] >= gridz ? gridz : 0);
         }
         float charge = epsilonFactor*posq[4*i+3];
-        __m128 zdata0to3 = _mm_set_ps(EXTRACT_FLOAT(data[3], 2), EXTRACT_FLOAT(data[2], 2), EXTRACT_FLOAT(data[1], 2), EXTRACT_FLOAT(data[0], 2));
-        float zdata4 = EXTRACT_FLOAT(data[4], 2);
+        fvec4 zdata0to3(data[0][2], data[1][2], data[2][2], data[3][2]);
+        float zdata4 = data[4][2];
         if (gridIndexZ+4 < gridz) {
             for (int ix = 0; ix < PME_ORDER; ix++) {
                 int xbase = gridIndexX+ix;
                 xbase -= (xbase >= gridx ? gridx : 0);
                 xbase = xbase*gridy*gridz;
-                float xdata = charge*EXTRACT_FLOAT(data[ix], 0);
+                float xdata = charge*data[ix][0];
                 for (int iy = 0; iy < PME_ORDER; iy++) {
                     int ybase = gridIndexY+iy;
                     ybase -= (ybase >= gridy ? gridy : 0);
                     ybase = xbase + ybase*gridz;
-                    float multiplier = xdata*EXTRACT_FLOAT(data[iy], 1);
-                    __m128 add0to3 = _mm_mul_ps(zdata0to3, _mm_set1_ps(multiplier));
-                    _mm_storeu_ps(&grid[ybase+gridIndexZ], _mm_add_ps(_mm_loadu_ps(&grid[ybase+gridIndexZ]), add0to3));
+                    float multiplier = xdata*data[iy][1];
+                    fvec4 add0to3 = zdata0to3*multiplier;
+                    (fvec4(&grid[ybase+gridIndexZ])+add0to3).store(&grid[ybase+gridIndexZ]);
                     grid[ybase+zindex[4]] += multiplier*zdata4;
                 }
             }
@@ -125,14 +123,14 @@ static void spreadCharge(int start, int end, float* posq, float* grid, int gridx
                 int xbase = gridIndexX+ix;
                 xbase -= (xbase >= gridx ? gridx : 0);
                 xbase = xbase*gridy*gridz;
-                float xdata = charge*EXTRACT_FLOAT(data[ix], 0);
+                float xdata = charge*data[ix][0];
                 for (int iy = 0; iy < PME_ORDER; iy++) {
                     int ybase = gridIndexY+iy;
                     ybase -= (ybase >= gridy ? gridy : 0);
                     ybase = xbase + ybase*gridz;
-                    float multiplier = xdata*EXTRACT_FLOAT(data[iy], 1);
-                    __m128 add0to3 = _mm_mul_ps(zdata0to3, _mm_set1_ps(multiplier));
-                    _mm_store_ps(temp, add0to3);
+                    float multiplier = xdata*data[iy][1];
+                    fvec4 add0to3 = zdata0to3*multiplier;
+                    add0to3.store(temp);
                     grid[ybase+zindex[0]] += temp[0];
                     grid[ybase+zindex[1]] += temp[1];
                     grid[ybase+zindex[2]] += temp[2];
@@ -245,51 +243,51 @@ static void reciprocalConvolution(int start, int end, fftwf_complex* grid, int g
 }
 
 static void interpolateForces(int start, int end, float* posq, float* force, float* grid, int gridx, int gridy, int gridz, int numParticles, Vec3 periodicBoxSize) {
-    __m128 boxSize = _mm_set_ps(0, (float) periodicBoxSize[2], (float) periodicBoxSize[1], (float) periodicBoxSize[0]);
-    __m128 invBoxSize = _mm_set_ps(0, (float) (1/periodicBoxSize[2]), (float) (1/periodicBoxSize[1]), (float) (1/periodicBoxSize[0]));
-    __m128 gridSize = _mm_set_ps(0, gridz, gridy, gridx);
-    __m128i gridSizeInt = _mm_set_epi32(0, gridz, gridy, gridx);
-    __m128 one  = _mm_set1_ps(1);
-    __m128 scale = _mm_set1_ps(1.0f/(PME_ORDER-1));
+    fvec4 boxSize((float) periodicBoxSize[0], (float) periodicBoxSize[1], (float) periodicBoxSize[2], 0);
+    fvec4 invBoxSize((float) (1/periodicBoxSize[0]), (float) (1/periodicBoxSize[1]), (float) (1/periodicBoxSize[2]), 0);
+    fvec4 gridSize(gridx, gridy, gridz, 0);
+    ivec4 gridSizeInt(gridx, gridy, gridz, 0);
+    fvec4 one(1);
+    fvec4 scale(1.0f/(PME_ORDER-1));
     const float epsilonFactor = sqrt(ONE_4PI_EPS0);
     for (int i = start; i < end; i++) {
         // Find the position relative to the nearest grid point.
         
-        __m128 pos = _mm_loadu_ps(&posq[4*i]);
-        __m128 posFloor = _mm_floor_ps(_mm_mul_ps(pos, invBoxSize));
-        __m128 posInBox = _mm_sub_ps(pos, _mm_mul_ps(boxSize, posFloor));
-        __m128 t = _mm_mul_ps(_mm_mul_ps(posInBox, invBoxSize), gridSize);
-        __m128i ti = _mm_cvttps_epi32(t);
-        __m128 dr = _mm_sub_ps(t, _mm_cvtepi32_ps(ti));
-        __m128i gridIndex = _mm_sub_epi32(ti, _mm_and_si128(gridSizeInt, _mm_cmpeq_epi32(ti, gridSizeInt)));
+        fvec4 pos(&posq[4*i]);
+        fvec4 posFloor = floor(pos*invBoxSize);
+        fvec4 posInBox = pos-boxSize*posFloor;
+        fvec4 t = posInBox*invBoxSize*gridSize;
+        ivec4 ti = t;
+        fvec4 dr = t-ti;
+        ivec4 gridIndex = ti-(gridSizeInt&ti==gridSizeInt);
         
         // Compute the B-spline coefficients.
         
-        __m128 data[PME_ORDER];
-        __m128 ddata[PME_ORDER];
-        data[PME_ORDER-1] = _mm_setzero_ps();
+        fvec4 data[PME_ORDER];
+        fvec4 ddata[PME_ORDER];
+        data[PME_ORDER-1] = 0.0f;
         data[1] = dr;
-        data[0] = _mm_sub_ps(one, dr);
+        data[0] = one-dr;
         for (int j = 3; j < PME_ORDER; j++) {
-            __m128 div = _mm_set1_ps(1.0f/(j-1));
-            data[j-1] = _mm_mul_ps(_mm_mul_ps(div, dr), data[j-2]);
+            fvec4 div(1.0f/(j-1));
+            data[j-1] = div*dr*data[j-2];
             for (int k = 1; k < j-1; k++)
-                data[j-k-1] = _mm_mul_ps(div, _mm_add_ps(_mm_mul_ps(_mm_add_ps(dr, _mm_set1_ps(k)), data[j-k-2]), _mm_mul_ps(_mm_sub_ps(_mm_set1_ps(j-k), dr), data[j-k-1])));
-            data[0] = _mm_mul_ps(_mm_mul_ps(div, _mm_sub_ps(one, dr)), data[0]);
+                data[j-k-1] = div*((dr+k)*data[j-k-2]+(fvec4(j-k)-dr)*data[j-k-1]);
+            data[0] = div*(one-dr)*data[0];
         }
-        ddata[0] = _mm_sub_ps(_mm_set1_ps(0), data[0]);
+        ddata[0] = -data[0];
         for (int j = 1; j < PME_ORDER; j++)
-            ddata[j] = _mm_sub_ps(data[j-1], data[j]);
-        data[PME_ORDER-1] = _mm_mul_ps(_mm_mul_ps(scale, dr), data[PME_ORDER-2]);
+            ddata[j] = data[j-1]-data[j];
+        data[PME_ORDER-1] = scale*dr*data[PME_ORDER-2];
         for (int j = 1; j < (PME_ORDER-1); j++)
-            data[PME_ORDER-j-1] = _mm_mul_ps(scale, _mm_add_ps(_mm_mul_ps(_mm_add_ps(dr, _mm_set1_ps(j)), data[PME_ORDER-j-2]), _mm_mul_ps(_mm_sub_ps(_mm_set1_ps(PME_ORDER-j), dr), data[PME_ORDER-j-1])));
-        data[0] = _mm_mul_ps(_mm_mul_ps(scale, _mm_sub_ps(one, dr)), data[0]);
+            data[PME_ORDER-j-1] = scale*((dr+j)*data[PME_ORDER-j-2]+(fvec4(PME_ORDER-j)-dr)*data[PME_ORDER-j-1]);
+        data[0] = scale*(one-dr)*data[0];
                 
         // Compute the force on this atom.
         
-        int gridIndexX = _mm_extract_epi32(gridIndex, 0);
-        int gridIndexY = _mm_extract_epi32(gridIndex, 1);
-        int gridIndexZ = _mm_extract_epi32(gridIndex, 2);
+        int gridIndexX = gridIndex[0];
+        int gridIndexY = gridIndex[1];
+        int gridIndexZ = gridIndex[2];
         if (gridIndexX < 0)
             return; // This happens when a simulation blows up and coordinates become NaN.
         int zindex[PME_ORDER];
@@ -297,34 +295,34 @@ static void interpolateForces(int start, int end, float* posq, float* force, flo
             zindex[j] = gridIndexZ+j;
             zindex[j] -= (zindex[j] >= gridz ? gridz : 0);
         }
-        __m128 zdata[PME_ORDER];
+        fvec4 zdata[PME_ORDER];
         for (int j = 0; j < PME_ORDER; j++)
-            zdata[j] = _mm_set_ps(0, EXTRACT_FLOAT(ddata[j], 2), EXTRACT_FLOAT(data[j], 2), EXTRACT_FLOAT(data[j], 2));
-        __m128 f = _mm_set1_ps(0);
+            zdata[j] = fvec4(data[j][2], data[j][2], ddata[j][2], 0);
+        fvec4 f = 0.0f;
         for (int ix = 0; ix < PME_ORDER; ix++) {
             int xbase = gridIndexX+ix;
             xbase -= (xbase >= gridx ? gridx : 0);
             xbase = xbase*gridy*gridz;
-            float dx = EXTRACT_FLOAT(data[ix], 0);
-            float ddx = EXTRACT_FLOAT(ddata[ix], 0);
-            __m128 xdata = _mm_set_ps(0, dx, dx, ddx);
+            float dx = data[ix][0];
+            float ddx = ddata[ix][0];
+            fvec4 xdata(ddx, dx, dx, 0);
 
             for (int iy = 0; iy < PME_ORDER; iy++) {
                 int ybase = gridIndexY+iy;
                 ybase -= (ybase >= gridy ? gridy : 0);
                 ybase = xbase + ybase*gridz;
-                float dy = EXTRACT_FLOAT(data[iy], 1);
-                float ddy = EXTRACT_FLOAT(ddata[iy], 1);
-                __m128 xydata = _mm_mul_ps(xdata, _mm_set_ps(0, dy, ddy, dy));
+                float dy = data[iy][1];
+                float ddy = ddata[iy][1];
+                fvec4 xydata = xdata*fvec4(dy, ddy, dy, 0);
 
                 for (int iz = 0; iz < PME_ORDER; iz++) {
-                    __m128 gridValue = _mm_set1_ps(grid[ybase+zindex[iz]]);
-                    f = _mm_add_ps(f, _mm_mul_ps(xydata, _mm_mul_ps(zdata[iz], gridValue)));
+                    fvec4 gridValue(grid[ybase+zindex[iz]]);
+                    f = f+xydata*zdata[iz]*gridValue;
                 }
             }
         }
-        f = _mm_mul_ps(invBoxSize, _mm_mul_ps(gridSize, _mm_mul_ps(f, _mm_set1_ps(-epsilonFactor*posq[4*i+3]))));
-        _mm_store_ps(&force[4*i], f);        
+        f = invBoxSize*gridSize*f*(-epsilonFactor*posq[4*i+3]);
+        f.store(&force[4*i]);
     }
 }
 
@@ -509,10 +507,10 @@ void CpuCalcPmeReciprocalForceKernel::runThread(int index) {
             threadWait();
             int numGrids = threadData.size();
             for (int i = gridStart; i < gridEnd; i += 4) {
-                __m128 sum = _mm_load_ps(&realGrid[i]);
+                fvec4 sum(&realGrid[i]);
                 for (int j = 1; j < numGrids; j++)
-                    sum = _mm_add_ps(sum, _mm_load_ps(&threadData[j]->tempGrid[i]));
-                _mm_store_ps(&realGrid[i], sum);
+                    sum += fvec4(&threadData[j]->tempGrid[i]);
+                sum.store(&realGrid[i]);
             }
             threadWait();
             if (lastBoxSize != periodicBoxSize) {

wrappers/python/simtk/openmm/app/amberprmtopfile.py

@@ -61,6 +61,11 @@ class GBn(object):
         return 'GBn'
 GBn = GBn()
 
+class GBn2(object):
+    def __repr__(self):
+        return 'GBn2'
+GBn2 = GBn2()
+
 class AmberPrmtopFile(object):
     """AmberPrmtopFile parses an AMBER prmtop file and constructs a Topology and (optionally) an OpenMM System from it."""
 
@@ -69,6 +74,7 @@ class AmberPrmtopFile(object):
         top = Topology()
         ## The Topology read from the prmtop file
         self.topology = top
+        self.elements = []
 
         # Load the prmtop file
 
@@ -96,7 +102,14 @@ class AmberPrmtopFile(object):
             if atomName in atomReplacements:
                 atomName = atomReplacements[atomName]
 
-            # Try to guess the element.
+            # Get the element from the prmtop file if available
+            if prmtop.has_atomic_number:
+                try:
+                    element = elem.Element.getByAtomicNumber(int(prmtop._raw_data['ATOMIC_NUMBER'][index]))
+                except KeyError:
+                    element = None
+            else:
+                # Try to guess the element from the atom name.
 
                 upper = atomName.upper()
                 if upper.startswith('CL'):
@@ -110,7 +123,9 @@ class AmberPrmtopFile(object):
                         element = elem.get_by_symbol(atomName[0])
                     except KeyError:
                         element = None
+
             top.addAtom(atomName, element, r)
+            self.elements.append(element)
 
         # Add bonds to the topology
 
@@ -137,7 +152,7 @@ class AmberPrmtopFile(object):
          - constraints (object=None) Specifies which bonds angles should be implemented with constraints.
            Allowed values are None, HBonds, AllBonds, or HAngles.
          - rigidWater (boolean=True) If true, water molecules will be fully rigid regardless of the value passed for the constraints argument
-         - implicitSolvent (object=None) If not None, the implicit solvent model to use.  Allowed values are HCT, OBC1, OBC2, or GBn.
+         - implicitSolvent (object=None) If not None, the implicit solvent model to use.  Allowed values are HCT, OBC1, OBC2, GBn, or GBn2.
          - soluteDielectric (float=1.0) The solute dielectric constant to use in the implicit solvent model.
          - solventDielectric (float=78.5) The solvent dielectric constant to use in the implicit solvent model.
          - removeCMMotion (boolean=True) If true, a CMMotionRemover will be added to the System
@@ -167,19 +182,22 @@ class AmberPrmtopFile(object):
             raise ValueError('Illegal value for constraints')
         if implicitSolvent is None:
             implicitString = None
-        elif implicitSolvent == HCT:
+        elif implicitSolvent is HCT:
             implicitString = 'HCT'
-        elif implicitSolvent == OBC1:
+        elif implicitSolvent is OBC1:
             implicitString = 'OBC1'
-        elif implicitSolvent == OBC2:
+        elif implicitSolvent is OBC2:
             implicitString = 'OBC2'
-        elif implicitSolvent == GBn:
+        elif implicitSolvent is GBn:
             implicitString = 'GBn'
+        elif implicitSolvent is GBn2:
+            implicitString = 'GBn2'
         else:
             raise ValueError('Illegal value for implicit solvent model')
         sys = amber_file_parser.readAmberSystem(prmtop_loader=self._prmtop, shake=constraintString, nonbondedCutoff=nonbondedCutoff,
                                                  nonbondedMethod=methodMap[nonbondedMethod], flexibleConstraints=False, gbmodel=implicitString,
-                                                 soluteDielectric=soluteDielectric, solventDielectric=solventDielectric, rigidWater=rigidWater)
+                                                 soluteDielectric=soluteDielectric, solventDielectric=solventDielectric, rigidWater=rigidWater,
+                                                 elements=self.elements)
         if hydrogenMass is not None:
             for atom1, atom2 in self.topology.bonds():
                 if atom1.element == elem.hydrogen:

wrappers/python/simtk/openmm/app/element.py

@@ -203,3 +203,8 @@ ununtrium =      Element(113, "ununtrium",      "Uut", 284*daltons)
 ununquadium =    Element(114, "ununquadium",    "Uuq", 289*daltons)
 ununpentium =    Element(115, "ununpentium",    "Uup", 288*daltons)
 ununhexium =     Element(116, "ununhexium",     "Uuh", 292*daltons)
+
+# Aliases to recognize common alternative spellings. Both the '==' and 'is'
+# relational operators will work with any chosen name
+sulphur = sulfur
+aluminium = aluminum

wrappers/python/simtk/openmm/app/internal/amber_file_parser.py

@@ -52,6 +52,7 @@ except:
 
 import simtk.unit as units
 import simtk.openmm
+from simtk.openmm.app import element as elem
 import customgbforces as customgb
 
 #=============================================================================================
@@ -472,34 +473,74 @@ class PrmtopLoader(object):
             self._excludedAtoms.append(atomList)
         return self._excludedAtoms
 
-    def getGBParms(self, symbls=None):
+    def getGBParms(self, gbmodel, elements):
         """Return list giving GB params, Radius and screening factor"""
-        try:
-            return self._gb_List
-        except AttributeError:
-            pass
-        self._gb_List=[]
-        radii=self._raw_data["RADII"]
-        screen=self._raw_data["SCREEN"]
+        gb_List=[]
+        radii=[float(r) for r in self._raw_data["RADII"]]
+        screen=[float(s) for s in self._raw_data["SCREEN"]]
+        if gbmodel == 'GBn2':
+            alpha = [0 for i in self._raw_data['RADII']]
+            beta = [0 for i in self._raw_data['RADII']]
+            gamma = [0 for i in self._raw_data['RADII']]
         # Update screening parameters for GBn if specified
-        if symbls:
-            for (i, symbl) in enumerate(symbls):
-                if symbl[0] in ('c', 'C'):
+        if gbmodel == 'GBn':
+            if elements is None:
+                raise Exception('GBn model requires element information')
+            for i, element in enumerate(elements):
+                if element is elem.carbon:
                     screen[i] = 0.48435382330
-                elif symbl[0] in ('h', 'H'):
+                elif element is elem.hydrogen:
                     screen[i] = 1.09085413633
-                elif symbl[0] in ('n', 'N'):
+                elif element is elem.nitrogen:
                     screen[i] = 0.700147318409
-                elif symbl[0] in ('o', 'O'):
+                elif element is elem.oxygen:
                     screen[i] = 1.06557401132
-                elif symbl[0] in ('s', 'S'):
+                elif element is elem.sulfur:
                     screen[i] = 0.602256336067
                 else:
                     screen[i] = 0.5
+        if gbmodel == 'GBn2':
+            if elements is None:
+                raise Exception('GBn2 model requires element information')
+            for i, element in enumerate(elements):
+                if element is elem.carbon:
+                    screen[i] = 1.058554
+                    alpha[i] = 0.733756
+                    beta[i] = 0.506378
+                    gamma[i] = 0.205844
+                elif element is elem.hydrogen:
+                    screen[i] = 1.425952
+                    alpha[i] = 0.788440
+                    beta[i] = 0.798699
+                    gamma[i] = 0.437334
+                elif element is elem.nitrogen:
+                    screen[i] = 0.733599
+                    alpha[i] = 0.503364
+                    beta[i] = 0.316828
+                    gamma[i] = 0.192915
+                elif element is elem.oxygen:
+                    screen[i] = 1.061039
+                    alpha[i] = 0.867814
+                    beta[i] = 0.876635
+                    gamma[i] = 0.387882
+                elif element is elem.sulfur:
+                    screen[i] = -0.703469
+                    alpha[i] = 0.867814
+                    beta[i] = 0.876635
+                    gamma[i] = 0.387882
+                else: # not optimized
+                    screen[i] = 0.5
+                    alpha[i] = 1.0
+                    beta[i] = 0.8
+                    gamma[i] = 4.85
         lengthConversionFactor = units.angstrom.conversion_factor_to(units.nanometer)
-        for iAtom in range(len(radii)):
-            self._gb_List.append((float(radii[iAtom])*lengthConversionFactor, float(screen[iAtom])))
-        return self._gb_List
+        if gbmodel == 'GBn2':
+            for rad, scr, alp, bet, gam in zip(radii, screen, alpha, beta, gamma):
+                gb_List.append((rad*lengthConversionFactor, scr, alp, bet, gam))
+        else:
+            for rad, scr in zip(radii, screen):
+                gb_List.append((rad*lengthConversionFactor, scr))
+        return gb_List
 
     def getBoxBetaAndDimensions(self):
         """Return periodic boundary box beta angle and dimensions"""
@@ -516,11 +557,17 @@ class PrmtopLoader(object):
     def has_scee_scnb(self):
         return ("SCEE_SCALE_FACTOR" in self._raw_data and "SCNB_SCALE_FACTOR" in self._raw_data)
 
+    @property
+    def has_atomic_number(self):
+        return 'ATOMIC_NUMBER' in self._raw_data
+
 #=============================================================================================
 # AMBER System builder (based on, but not identical to, systemManager from 'zander')
 #=============================================================================================
 
-def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmodel=None, soluteDielectric=1.0, solventDielectric=78.5, nonbondedCutoff=None, nonbondedMethod='NoCutoff', scee=None, scnb=None, mm=None, verbose=False, EwaldErrorTolerance=None, flexibleConstraints=True, rigidWater=True):
+def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmodel=None, soluteDielectric=1.0, solventDielectric=78.5,
+                    nonbondedCutoff=None, nonbondedMethod='NoCutoff', scee=None, scnb=None, mm=None, verbose=False,
+                    EwaldErrorTolerance=None, flexibleConstraints=True, rigidWater=True, elements=None):
     """
     Create an OpenMM System from an Amber prmtop file.
 
@@ -813,15 +860,12 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
     if gbmodel is not None:
         if verbose: print "Adding GB parameters..."
         charges = prmtop.getCharges()
-        symbls = None
         cutoff = None
         if nonbondedMethod != 'NoCutoff':
             cutoff = nonbondedCutoff
             if units.is_quantity(cutoff):
                 cutoff = cutoff.value_in_unit(units.nanometers)
-        if gbmodel == 'GBn':
-            symbls = prmtop.getAtomTypes()
-        gb_parms = prmtop.getGBParms(symbls)
+        gb_parms = prmtop.getGBParms(gbmodel, elements)
         if gbmodel == 'HCT':
             gb = customgb.GBSAHCTForce(solventDielectric, soluteDielectric, 'ACE', cutoff)
         elif gbmodel == 'OBC1':
@@ -832,11 +876,16 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
             gb.setSolventDielectric(solventDielectric)
         elif gbmodel == 'GBn':
             gb = customgb.GBSAGBnForce(solventDielectric, soluteDielectric, 'ACE', cutoff)
+        elif gbmodel == 'GBn2':
+            gb = customgb.GBSAGBn2Force(solventDielectric, soluteDielectric, 'ACE', cutoff)
         else:
             raise Exception("Illegal value specified for implicit solvent model")
         for iAtom in range(prmtop.getNumAtoms()):
             if gbmodel == 'OBC2':
                 gb.addParticle(charges[iAtom], gb_parms[iAtom][0], gb_parms[iAtom][1])
+            elif gbmodel == 'GBn2':
+                gb.addParticle([charges[iAtom], gb_parms[iAtom][0], gb_parms[iAtom][1],
+                                gb_parms[iAtom][2], gb_parms[iAtom][3], gb_parms[iAtom][4]])
             else:
                 gb.addParticle([charges[iAtom], gb_parms[iAtom][0], gb_parms[iAtom][1]])
         system.addForce(gb)

wrappers/python/simtk/openmm/app/internal/customgbforces.py

@@ -8,7 +8,7 @@ Medical Research, grant U54 GM072970. See https://simtk.org.
 
 Portions copyright (c) 2012 University of Virginia and the Authors.
 Authors: Christoph Klein, Michael R. Shirts
-Contributors:
+Contributors: Jason M. Swails
 
 Permission is hereby granted, free of charge, to any person obtaining a
 copy of this software and associated documentation files (the "Software"),
@@ -29,8 +29,9 @@ OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
 USE OR OTHER DEALINGS IN THE SOFTWARE.
 """
 
+from __future__ import division
+
 from simtk.openmm import CustomGBForce
-import sys, pdb, pickle
 
 d0=[2.26685,2.32548,2.38397,2.44235,2.50057,2.55867,2.61663,2.67444,
     2.73212,2.78965,2.84705,2.9043,2.96141,3.0184,3.07524,3.13196,
@@ -212,11 +213,11 @@ def GBSAHCTForce(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=Non
 
     custom = CustomGBForce()
 
-    custom.addPerParticleParameter("q");
-    custom.addPerParticleParameter("radius");
-    custom.addPerParticleParameter("scale");
-    custom.addGlobalParameter("solventDielectric", solventDielectric);
-    custom.addGlobalParameter("soluteDielectric", soluteDielectric);
+    custom.addPerParticleParameter("q")
+    custom.addPerParticleParameter("radius")
+    custom.addPerParticleParameter("scale")
+    custom.addGlobalParameter("solventDielectric", solventDielectric)
+    custom.addGlobalParameter("soluteDielectric", soluteDielectric)
     custom.addGlobalParameter("offset", 0.009)
     custom.addComputedValue("I", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
                                   "U=r+sr2;"
@@ -237,11 +238,11 @@ def GBSAOBC1Force(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=No
 
     custom = CustomGBForce()
 
-    custom.addPerParticleParameter("q");
-    custom.addPerParticleParameter("radius");
-    custom.addPerParticleParameter("scale");
-    custom.addGlobalParameter("solventDielectric", solventDielectric);
-    custom.addGlobalParameter("soluteDielectric", soluteDielectric);
+    custom.addPerParticleParameter("q")
+    custom.addPerParticleParameter("radius")
+    custom.addPerParticleParameter("scale")
+    custom.addGlobalParameter("solventDielectric", solventDielectric)
+    custom.addGlobalParameter("soluteDielectric", soluteDielectric)
     custom.addGlobalParameter("offset", 0.009)
     custom.addComputedValue("I",  "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
                                   "U=r+sr2;"
@@ -262,11 +263,11 @@ def GBSAOBC2Force(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=No
 
     custom = CustomGBForce()
 
-    custom.addPerParticleParameter("q");
-    custom.addPerParticleParameter("radius");
-    custom.addPerParticleParameter("scale");
-    custom.addGlobalParameter("solventDielectric", solventDielectric);
-    custom.addGlobalParameter("soluteDielectric", soluteDielectric);
+    custom.addPerParticleParameter("q")
+    custom.addPerParticleParameter("radius")
+    custom.addPerParticleParameter("scale")
+    custom.addGlobalParameter("solventDielectric", solventDielectric)
+    custom.addGlobalParameter("soluteDielectric", soluteDielectric)
     custom.addGlobalParameter("offset", 0.009)
     custom.addComputedValue("I",  "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
                                   "U=r+sr2;"
@@ -296,12 +297,12 @@ def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=Non
 
     custom = CustomGBForce()
 
-    custom.addPerParticleParameter("q");
-    custom.addPerParticleParameter("radius");
-    custom.addPerParticleParameter("scale");
+    custom.addPerParticleParameter("q")
+    custom.addPerParticleParameter("radius")
+    custom.addPerParticleParameter("scale")
 
-    custom.addGlobalParameter("solventDielectric", solventDielectric);
-    custom.addGlobalParameter("soluteDielectric", soluteDielectric);
+    custom.addGlobalParameter("solventDielectric", solventDielectric)
+    custom.addGlobalParameter("soluteDielectric", soluteDielectric)
     custom.addGlobalParameter("offset", 0.009)
     custom.addGlobalParameter("neckScale", 0.361825)
     custom.addGlobalParameter("neckCut", 0.68)
@@ -323,3 +324,50 @@ def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=Non
                               "psi=I*or; or=radius-offset", CustomGBForce.SingleParticle)
     _createEnergyTerms(custom, SA, cutoff)
     return custom
+
+"""
+Amber Equivalents: igb = 8
+"""
+def GBSAGBn2Force(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=None):
+
+
+    """
+    Indexing for tables:
+        input: radius1, radius2
+        index = (radius2*200-20)*21 + (radius1*200-20)
+        output: index of desired value in row-by-row, 1D version of Tables 3 & 4
+    """
+
+
+    custom = CustomGBForce()
+
+    custom.addPerParticleParameter("q")
+    custom.addPerParticleParameter("radius")
+    custom.addPerParticleParameter("scale")
+    custom.addPerParticleParameter("alpha")
+    custom.addPerParticleParameter("beta")
+    custom.addPerParticleParameter("gamma")
+
+    custom.addGlobalParameter("solventDielectric", solventDielectric)
+    custom.addGlobalParameter("soluteDielectric", soluteDielectric)
+    custom.addGlobalParameter("offset", 0.0195141)
+    custom.addGlobalParameter("neckScale", 0.826836)
+    custom.addGlobalParameter("neckCut", 0.68)
+
+    custom.addFunction("getd0", d0, 0, 440)
+    custom.addFunction("getm0", m0, 0, 440)
+
+    custom.addComputedValue("I",  "Ivdw+neckScale*Ineck;"
+                                  "Ineck=step(radius1+radius2+neckCut-r)*getm0(index)/(1+100*(r-getd0(index))^2+0.3*1000000*(r-getd0(index))^6);"
+                                  "index = (radius2*200-20)*21 + (radius1*200-20);"
+                                  "Ivdw=step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
+                                  "U=r+sr2;"
+                                  "L=max(or1, D);"
+                                  "D=abs(r-sr2);"
+                                  "sr2 = scale2*or2;"
+                                  "or1 = radius1-offset; or2 = radius2-offset", CustomGBForce.ParticlePairNoExclusions)
+
+    custom.addComputedValue("B", "1/(1/or-tanh(alpha*psi-beta*psi^2+gamma*psi^3)/radius);"
+                              "psi=I*or; or=radius-offset", CustomGBForce.SingleParticle)
+    _createEnergyTerms(custom, SA, cutoff)
+    return custom

wrappers/python/simtk/openmm/app/internal/pdbstructure.py

@@ -436,7 +436,7 @@ class Chain(object):
         self._finalize()
 
     def get_residue(self, residue_number, insertion_code=' '):
-        return residues_by_num_icode[str(residue_number) + insertion_code]
+        return self.residues_by_num_icode[str(residue_number) + insertion_code]
 
     def __contains__(self, residue_number):
         return self.residues_by_number.__contains__(residue_number)
@@ -916,7 +916,6 @@ if __name__=='__main__':
     import doctest, sys
     doctest.testmod(sys.modules[__name__])
 
-    import sys
     import os
     import gzip
     import re

wrappers/python/simtk/openmm/app/statedatareporter.py

@@ -31,9 +31,10 @@ USE OR OTHER DEALINGS IN THE SOFTWARE.
 __author__ = "Peter Eastman"
 __version__ = "1.0"
 
+import bz2
+import gzip
 import simtk.openmm as mm
 import simtk.unit as unit
-from simtk.openmm.app import PDBFile
 import math
 
 class StateDataReporter(object):
@@ -67,6 +68,13 @@ class StateDataReporter(object):
         self._reportInterval = reportInterval
         self._openedFile = isinstance(file, str)
         if self._openedFile:
+            # Detect the desired compression scheme from the filename extension
+            # and open all files unbuffered
+            if file.endswith('.gz'):
+                self._out = gzip.GzipFile(fileobj=open(file, 'wb', 0))
+            elif file.endswith('.bz2'):
+                self._out = bz2.BZ2File(file, 'w', 0)
+            else:
                 self._out = open(file, 'w', 0)
         else:
             self._out = file
@@ -109,7 +117,10 @@ class StateDataReporter(object):
             self._initializeConstants(simulation)
             headers = self._constructHeaders()
             print >>self._out, '#"%s"' % ('"'+self._separator+'"').join(headers)
+            try:
                 self._out.flush()
+            except AttributeError:
+                pass
             self._hasInitialized = True
 
         # Check for errors.
@@ -120,7 +131,10 @@ class StateDataReporter(object):
 
         # Write the values.
         print >>self._out, self._separator.join(str(v) for v in values)
+        try:
             self._out.flush()
+        except AttributeError:
+            pass
 
     def _constructReportValues(self, simulation, state):
         """Query the simulation for the current state of our observables of interest.

wrappers/python/src/swig_doxygen/swigInputConfig.py

@@ -145,6 +145,8 @@ SKIP_METHODS = [('State',),
                 ('IntegrateDrudeSCFStepKernel',),
                 ('XmlSerializer',  'serialize'),
                 ('XmlSerializer',  'deserialize'),
+                ('fvec4',),
+                ('ivec4',),
 ]
 
 # The build script assumes method args that are non-const references are