Merge remote-tracking branch 'upstream/master'

plugins/rpmd/openmmapi/src/RPMDMonteCarloBarostat.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2010-2015 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/RPMDMonteCarloBarostat.h"
+#include "openmm/internal/RPMDMonteCarloBarostatImpl.h"
+
+using namespace OpenMM;
+
+RPMDMonteCarloBarostat::RPMDMonteCarloBarostat(double defaultPressure, int frequency) :
+        defaultPressure(defaultPressure), frequency(frequency) {
+    setRandomNumberSeed(0);
+}
+
+ForceImpl* RPMDMonteCarloBarostat::createImpl() const {
+    return new RPMDMonteCarloBarostatImpl(*this);
+}

plugins/rpmd/openmmapi/src/RPMDMonteCarloBarostatImpl.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2010-2015 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/internal/RPMDMonteCarloBarostatImpl.h"
+#include "openmm/internal/ContextImpl.h"
+#include "openmm/internal/OSRngSeed.h"
+#include "openmm/Context.h"
+#include "openmm/kernels.h"
+#include "openmm/OpenMMException.h"
+#include "openmm/RPMDIntegrator.h"
+#include <cmath>
+#include <vector>
+#include <algorithm>
+
+using namespace OpenMM;
+using namespace OpenMM_SFMT;
+using std::vector;
+
+const float BOLTZMANN = 1.380658e-23f; // (J/K)
+const float AVOGADRO = 6.0221367e23f;
+const float RGAS = BOLTZMANN*AVOGADRO; // (J/(mol K))
+const float BOLTZ = RGAS/1000;         // (kJ/(mol K))
+
+RPMDMonteCarloBarostatImpl::RPMDMonteCarloBarostatImpl(const RPMDMonteCarloBarostat& owner) : owner(owner), step(0) {
+}
+
+void RPMDMonteCarloBarostatImpl::initialize(ContextImpl& context) {
+    RPMDIntegrator* integrator = dynamic_cast<RPMDIntegrator*>(&context.getIntegrator());
+    if (integrator == NULL)
+        throw OpenMMException("RPMDMonteCarloBarostat must be used with an RPMDIntegrator");;
+    if (!integrator->getApplyThermostat())
+        throw OpenMMException("RPMDMonteCarloBarostat requires the integrator's thermostat to be enabled");;
+    kernel = context.getPlatform().createKernel(ApplyMonteCarloBarostatKernel::Name(), context);
+    kernel.getAs<ApplyMonteCarloBarostatKernel>().initialize(context.getSystem(), owner);
+    savedPositions.resize(integrator->getNumCopies());
+    Vec3 box[3];
+    context.getPeriodicBoxVectors(box[0], box[1], box[2]);
+    double volume = box[0][0]*box[1][1]*box[2][2];
+    volumeScale = 0.01*volume;
+    numAttempted = 0;
+    numAccepted = 0;
+    int randSeed = owner.getRandomNumberSeed();
+    // A random seed of 0 means use a unique one
+    if (randSeed == 0) randSeed = osrngseed();
+    init_gen_rand(randSeed, random);
+}
+
+void RPMDMonteCarloBarostatImpl::updateRPMDState(ContextImpl& context) {
+    if (++step < owner.getFrequency() || owner.getFrequency() == 0)
+        return;
+    step = 0;
+
+    // Compute the current potential energy.
+
+    RPMDIntegrator& integrator = dynamic_cast<RPMDIntegrator&>(context.getIntegrator());
+
+    // Record the initial positions and energy
+
+    double initialEnergy = 0;
+    int numCopies = integrator.getNumCopies();
+    for (int i = 0; i < numCopies; i++) {
+        State state = integrator.getState(i, State::Positions | State::Energy);
+        savedPositions[i] = state.getPositions();
+        initialEnergy += state.getPotentialEnergy();
+    }
+
+    // Compute the centroid.
+
+    int numParticles = context.getSystem().getNumParticles();
+    vector<Vec3> centroid(numParticles, Vec3());
+    for (int i = 0; i < numParticles; i++) {
+        for (int j = 0; j < numCopies; j++)
+            centroid[i] += savedPositions[j][i];
+        centroid[i] *= 1.0/numCopies;
+    }
+
+    // Modify the periodic box size and scale the coordinates of the centroid.
+
+    Vec3 box[3];
+    context.getPeriodicBoxVectors(box[0], box[1], box[2]);
+    double volume = box[0][0]*box[1][1]*box[2][2];
+    double deltaVolume = volumeScale*2*(genrand_real2(random)-0.5);
+    double newVolume = volume+deltaVolume;
+    double lengthScale = std::pow(newVolume/volume, 1.0/3.0);
+    context.setPositions(centroid);
+    kernel.getAs<ApplyMonteCarloBarostatKernel>().scaleCoordinates(context, lengthScale, lengthScale, lengthScale);
+    context.getOwner().setPeriodicBoxVectors(box[0]*lengthScale, box[1]*lengthScale, box[2]*lengthScale);
+    State scaledState = context.getOwner().getState(State::Positions);
+
+    // Now apply the same offset to all the copies.
+
+    vector<Vec3> delta(numParticles);
+    for (int i = 0; i < numParticles; i++)
+        delta[i] = scaledState.getPositions()[i]-centroid[i];
+    double finalEnergy = 0;
+    vector<Vec3> positions(numParticles);
+    for (int copy = 0; copy < numCopies; copy++) {
+        for (int i = 0; i < numParticles; i++)
+            positions[i] = savedPositions[copy][i]+delta[i];
+        integrator.setPositions(copy, positions);
+        finalEnergy += integrator.getState(copy, State::Energy).getPotentialEnergy();
+    }
+
+    // Compute the energy of the modified system.
+
+    double pressure = context.getParameter(RPMDMonteCarloBarostat::Pressure())*(AVOGADRO*1e-25);
+    double kT = BOLTZ*integrator.getTemperature();
+    double w = (finalEnergy-initialEnergy)/numCopies + pressure*deltaVolume - context.getMolecules().size()*kT*std::log(newVolume/volume);
+    if (w > 0 && genrand_real2(random) > std::exp(-w/kT)) {
+        // Reject the step.
+
+        for (int copy = 0; copy < numCopies; copy++)
+            integrator.setPositions(copy, savedPositions[copy]);
+        context.getOwner().setPeriodicBoxVectors(box[0], box[1], box[2]);
+        volume = newVolume;
+    }
+    else
+        numAccepted++;
+    numAttempted++;
+    if (numAttempted >= 10) {
+        if (numAccepted < 0.25*numAttempted) {
+            volumeScale /= 1.1;
+            numAttempted = 0;
+            numAccepted = 0;
+        }
+        else if (numAccepted > 0.75*numAttempted) {
+            volumeScale = std::min(volumeScale*1.1, volume*0.3);
+            numAttempted = 0;
+            numAccepted = 0;
+        }
+    }
+}
+
+std::map<std::string, double> RPMDMonteCarloBarostatImpl::getDefaultParameters() {
+    std::map<std::string, double> parameters;
+    parameters[RPMDMonteCarloBarostat::Pressure()] = getOwner().getDefaultPressure();
+    return parameters;
+}
+
+std::vector<std::string> RPMDMonteCarloBarostatImpl::getKernelNames() {
+    std::vector<std::string> names;
+    names.push_back(ApplyMonteCarloBarostatKernel::Name());
+    return names;
+}

plugins/rpmd/platforms/cuda/src/CudaRpmdKernels.cpp

@@ -257,13 +257,8 @@ void CudaIntegrateRPMDStepKernel::computeForces(ContextImpl& context) {
         context.updateContextState();
         Vec3 finalBox[3];
         context.getPeriodicBoxVectors(finalBox[0], finalBox[1], finalBox[2]);
-        if (initialBox[0] != finalBox[0] || initialBox[1] != finalBox[1] || initialBox[2] != finalBox[2]) {
-            // A barostat was applied during updateContextState().  Adjust the particle positions in all the
-            // other copies to match this one.
-            
-            void* args[] = {&positions->getDevicePointer(), &cu.getPosq().getDevicePointer(), &cu.getAtomIndexArray().getDevicePointer(), &i};
-            cu.executeKernel(translateKernel, args, cu.getNumAtoms());
-        }
+        if (initialBox[0] != finalBox[0] || initialBox[1] != finalBox[1] || initialBox[2] != finalBox[2])
+            throw OpenMMException("Standard barostats cannot be used with RPMDIntegrator.  Use RPMDMonteCarloBarostat instead.");
         context.calcForcesAndEnergy(true, false, groupsNotContracted);
         void* copyFromContextArgs[] = {&cu.getForce().getDevicePointer(), &forces->getDevicePointer(), &cu.getVelm().getDevicePointer(),
                 &velocities->getDevicePointer(), &cu.getPosq().getDevicePointer(), &positions->getDevicePointer(), &cu.getAtomIndexArray().getDevicePointer(), &i};

plugins/rpmd/platforms/cuda/tests/TestCudaRpmd.cpp

@@ -38,11 +38,11 @@
 #include "openmm/Context.h"
 #include "openmm/CustomNonbondedForce.h"
 #include "openmm/HarmonicBondForce.h"
-#include "openmm/MonteCarloBarostat.h"
 #include "openmm/NonbondedForce.h"
 #include "openmm/Platform.h"
 #include "openmm/System.h"
 #include "openmm/RPMDIntegrator.h"
+#include "openmm/RPMDMonteCarloBarostat.h"
 #include "openmm/VirtualSite.h"
 #include "SimTKOpenMMUtilities.h"
 #include "sfmt/SFMT.h"
@@ -371,7 +371,7 @@ void testContractions() {
     system.addForce(bonds);
     NonbondedForce* nonbonded = new NonbondedForce();
     nonbonded->setCutoffDistance(cutoff);
-    nonbonded->setNonbondedMethod(NonbondedForce::PME);
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
     nonbonded->setForceGroup(1);
     nonbonded->setReciprocalSpaceForceGroup(2);
     system.addForce(nonbonded);
@@ -491,11 +491,11 @@ void testWithBarostat() {
     system.addForce(bonds);
     NonbondedForce* nonbonded = new NonbondedForce();
     nonbonded->setCutoffDistance(cutoff);
-    nonbonded->setNonbondedMethod(NonbondedForce::PME);
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
     nonbonded->setForceGroup(1);
     nonbonded->setReciprocalSpaceForceGroup(2);
     system.addForce(nonbonded);
-    system.addForce(new MonteCarloBarostat(0.5, temperature));
+    system.addForce(new RPMDMonteCarloBarostat(0.5, 10));
 
     // Create a cloud of molecules.
 

plugins/rpmd/platforms/opencl/src/OpenCLRpmdKernels.cpp

@@ -273,13 +273,8 @@ void OpenCLIntegrateRPMDStepKernel::computeForces(ContextImpl& context) {
         context.updateContextState();
         Vec3 finalBox[3];
         context.getPeriodicBoxVectors(finalBox[0], finalBox[1], finalBox[2]);
-        if (initialBox[0] != finalBox[0] || initialBox[1] != finalBox[1] || initialBox[2] != finalBox[2]) {
-            // A barostat was applied during updateContextState().  Adjust the particle positions in all the
-            // other copies to match this one.
-            
-            translateKernel.setArg<cl_int>(3, i);
-            cl.executeKernel(translateKernel, cl.getNumAtoms());
-        }
+        if (initialBox[0] != finalBox[0] || initialBox[1] != finalBox[1] || initialBox[2] != finalBox[2])
+            throw OpenMMException("Standard barostats cannot be used with RPMDIntegrator.  Use RPMDMonteCarloBarostat instead.");
         context.calcForcesAndEnergy(true, false, groupsNotContracted);
         copyFromContextKernel.setArg<cl_int>(7, i);
         cl.executeKernel(copyFromContextKernel, cl.getNumAtoms());

plugins/rpmd/platforms/opencl/tests/TestOpenCLRpmd.cpp

@@ -38,11 +38,11 @@
 #include "openmm/Context.h"
 #include "openmm/CustomNonbondedForce.h"
 #include "openmm/HarmonicBondForce.h"
-#include "openmm/MonteCarloBarostat.h"
 #include "openmm/NonbondedForce.h"
 #include "openmm/Platform.h"
 #include "openmm/System.h"
 #include "openmm/RPMDIntegrator.h"
+#include "openmm/RPMDMonteCarloBarostat.h"
 #include "openmm/VirtualSite.h"
 #include "SimTKOpenMMUtilities.h"
 #include "sfmt/SFMT.h"
@@ -372,7 +372,7 @@ void testContractions() {
     system.addForce(bonds);
     NonbondedForce* nonbonded = new NonbondedForce();
     nonbonded->setCutoffDistance(cutoff);
-    nonbonded->setNonbondedMethod(NonbondedForce::PME);
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
     nonbonded->setForceGroup(1);
     nonbonded->setReciprocalSpaceForceGroup(2);
     system.addForce(nonbonded);
@@ -492,11 +492,11 @@ void testWithBarostat() {
     system.addForce(bonds);
     NonbondedForce* nonbonded = new NonbondedForce();
     nonbonded->setCutoffDistance(cutoff);
-    nonbonded->setNonbondedMethod(NonbondedForce::PME);
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
     nonbonded->setForceGroup(1);
     nonbonded->setReciprocalSpaceForceGroup(2);
     system.addForce(nonbonded);
-    system.addForce(new MonteCarloBarostat(0.5, temperature));
+    system.addForce(new RPMDMonteCarloBarostat(0.5, 10));
 
     // Create a cloud of molecules.
 

plugins/rpmd/platforms/reference/src/ReferenceRpmdKernels.cpp

@@ -280,17 +280,8 @@ void ReferenceIntegrateRPMDStepKernel::computeForces(ContextImpl& context, const
         context.updateContextState();
         Vec3 finalBox[3];
         context.getPeriodicBoxVectors(finalBox[0], finalBox[1], finalBox[2]);
-        if (initialBox[0] != finalBox[0] || initialBox[1] != finalBox[1] || initialBox[2] != finalBox[2]) {
-            // A barostat was applied during updateContextState().  Adjust the particle positions in all the
-            // other copies to match this one.
-            
-            for (int j = 0; j < numParticles; j++) {
-                Vec3 delta = pos[j]-positions[i][j];
-                for (int k = 0; k < totalCopies; k++)
-                    if (k != i)
-                        positions[k][j] += delta;
-            }
-        }
+        if (initialBox[0] != finalBox[0] || initialBox[1] != finalBox[1] || initialBox[2] != finalBox[2])
+            throw OpenMMException("Standard barostats cannot be used with RPMDIntegrator.  Use RPMDMonteCarloBarostat instead.");
         positions[i] = pos;
         velocities[i] = vel;
         context.calcForcesAndEnergy(true, false, groupsNotContracted);

plugins/rpmd/platforms/reference/tests/TestReferenceRpmd.cpp

@@ -37,11 +37,11 @@
 #include "openmm/CMMotionRemover.h"
 #include "openmm/Context.h"
 #include "openmm/HarmonicBondForce.h"
-#include "openmm/MonteCarloBarostat.h"
 #include "openmm/NonbondedForce.h"
 #include "openmm/Platform.h"
 #include "openmm/System.h"
 #include "openmm/RPMDIntegrator.h"
+#include "openmm/RPMDMonteCarloBarostat.h"
 #include "openmm/VirtualSite.h"
 #include "SimTKOpenMMUtilities.h"
 #include "sfmt/SFMT.h"
@@ -255,7 +255,7 @@ void testContractions() {
     system.addForce(bonds);
     NonbondedForce* nonbonded = new NonbondedForce();
     nonbonded->setCutoffDistance(cutoff);
-    nonbonded->setNonbondedMethod(NonbondedForce::PME);
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
     nonbonded->setForceGroup(1);
     nonbonded->setReciprocalSpaceForceGroup(2);
     system.addForce(nonbonded);
@@ -375,11 +375,11 @@ void testWithBarostat() {
     system.addForce(bonds);
     NonbondedForce* nonbonded = new NonbondedForce();
     nonbonded->setCutoffDistance(cutoff);
-    nonbonded->setNonbondedMethod(NonbondedForce::PME);
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
     nonbonded->setForceGroup(1);
     nonbonded->setReciprocalSpaceForceGroup(2);
     system.addForce(nonbonded);
-    system.addForce(new MonteCarloBarostat(0.5, temperature));
+    system.addForce(new RPMDMonteCarloBarostat(0.5, 10));
 
     // Create a cloud of molecules.
 

serialization/src/StateProxy.cpp

@@ -50,7 +50,7 @@ void StateProxy::serialize(const void* object, SerializationNode& node) const {
     boxVectorsNode.createChildNode("A").setDoubleProperty("x", a[0]).setDoubleProperty("y", a[1]).setDoubleProperty("z", a[2]);
     boxVectorsNode.createChildNode("B").setDoubleProperty("x", b[0]).setDoubleProperty("y", b[1]).setDoubleProperty("z", b[2]);
     boxVectorsNode.createChildNode("C").setDoubleProperty("x", c[0]).setDoubleProperty("y", c[1]).setDoubleProperty("z", c[2]);
-    try {
+    if ((s.getDataTypes()&State::Parameters) != 0) {
         s.getParameters();
         SerializationNode& parametersNode = node.createChildNode("Parameters");
         map<string, double> stateParams = s.getParameters();
@@ -58,46 +58,36 @@ void StateProxy::serialize(const void* object, SerializationNode& node) const {
         for (it = stateParams.begin(); it!=stateParams.end();it++) {
             parametersNode.setDoubleProperty(it->first, it->second);
         }
-    } catch (const OpenMMException &) {
-        // do nothing
     }
-    try {
+    if ((s.getDataTypes()&State::Energy) != 0) {
         s.getPotentialEnergy();
         SerializationNode& energiesNode = node.createChildNode("Energies");
         energiesNode.setDoubleProperty("PotentialEnergy", s.getPotentialEnergy());
         energiesNode.setDoubleProperty("KineticEnergy", s.getKineticEnergy());
-    } catch (const OpenMMException &) {
-        // do nothing
     }
-    try {
+    if ((s.getDataTypes()&State::Positions) != 0) {
         s.getPositions();
         SerializationNode& positionsNode = node.createChildNode("Positions");
         vector<Vec3> statePositions = s.getPositions();
         for (int i=0; i<statePositions.size();i++) {
            positionsNode.createChildNode("Position").setDoubleProperty("x", statePositions[i][0]).setDoubleProperty("y", statePositions[i][1]).setDoubleProperty("z", statePositions[i][2]);
         }
-    } catch (const OpenMMException &) {
-        // do nothing
     }
-    try {
+    if ((s.getDataTypes()&State::Velocities) != 0) {
         s.getVelocities();
         SerializationNode& velocitiesNode = node.createChildNode("Velocities");
         vector<Vec3> stateVelocities = s.getVelocities();
         for (int i=0; i<stateVelocities.size();i++) {
            velocitiesNode.createChildNode("Velocity").setDoubleProperty("x", stateVelocities[i][0]).setDoubleProperty("y", stateVelocities[i][1]).setDoubleProperty("z", stateVelocities[i][2]);
         }
-    } catch (const OpenMMException &) {
-        // do nothing
     }
-    try {
+    if ((s.getDataTypes()&State::Forces) != 0) {
         s.getForces();
         SerializationNode& forcesNode = node.createChildNode("Forces");
         vector<Vec3> stateForces = s.getForces();
         for (int i=0; i<stateForces.size();i++) {
             forcesNode.createChildNode("Force").setDoubleProperty("x", stateForces[i][0]).setDoubleProperty("y", stateForces[i][1]).setDoubleProperty("z", stateForces[i][2]);
         }
-    } catch (const OpenMMException &) {
-        // do nothing
     }
 }
 
@@ -113,83 +103,54 @@ void* StateProxy::deserialize(const SerializationNode& node) const {
     const SerializationNode& CVec = boxVectorsNode.getChildNode("C");
     Vec3 outCVec(CVec.getDoubleProperty("x"),CVec.getDoubleProperty("y"),CVec.getDoubleProperty("z"));
     int types = 0;
-    map<string, double> outStateParams;
-    try {
-        const SerializationNode& parametersNode = node.getChildNode("Parameters");
-        // inStateParams is really a <string,double> pair, where string is the name and double is the value
-        // but we want to avoid casting a string to a double and instead use the built in routines,
-        map<string, string> inStateParams = parametersNode.getProperties();
-        for (map<string, string>::const_iterator pit = inStateParams.begin(); pit != inStateParams.end(); pit++) {
-            outStateParams[pit->first] = parametersNode.getDoubleProperty(pit->first);
+    vector<int> arraySizes;
+    State::StateBuilder builder(outTime);
+    const vector<SerializationNode>& children = node.getChildren();
+    for (int j = 0; j < (int) children.size(); j++) {
+        const SerializationNode& child = children[j];
+        if (child.getName() == "Parameters") {
+            map<string, double> outStateParams;
+            // inStateParams is really a <string,double> pair, where string is the name and double is the value
+            // but we want to avoid casting a string to a double and instead use the built in routines,
+            map<string, string> inStateParams = child.getProperties();
+            for (map<string, string>::const_iterator pit = inStateParams.begin(); pit != inStateParams.end(); pit++) {
+                outStateParams[pit->first] = child.getDoubleProperty(pit->first);
+            }
+            builder.setParameters(outStateParams);
         }
-        types = types | State::Parameters;
-    } catch (const OpenMMException &) {
-        // do nothing
-    }
-    double potentialEnergy;
-    double kineticEnergy;
-    try {
-        const SerializationNode& energiesNode = node.getChildNode("Energies");
-        potentialEnergy = energiesNode.getDoubleProperty("PotentialEnergy");
-        kineticEnergy = energiesNode.getDoubleProperty("KineticEnergy");
-        types = types | State::Energy;
-    } catch (const OpenMMException &) {
-        // do nothing    
-    }
-    vector<Vec3> outPositions;
-    vector<Vec3> outVelocities;
-    vector<Vec3> outForces;
-    try {
-        const SerializationNode& positionsNode = node.getChildNode("Positions");
-        for (int i = 0; i < (int) positionsNode.getChildren().size(); i++) {
-            const SerializationNode& particle = positionsNode.getChildren()[i];
-            outPositions.push_back(Vec3(particle.getDoubleProperty("x"),particle.getDoubleProperty("y"),particle.getDoubleProperty("z")));
+        else if (child.getName() == "Energies") {
+            double potentialEnergy = child.getDoubleProperty("PotentialEnergy");
+            double kineticEnergy = child.getDoubleProperty("KineticEnergy");
+            builder.setEnergy(kineticEnergy, potentialEnergy);
         }
-        types = types | State::Positions;
-    } catch (const OpenMMException &) {
-        // do nothing    
-    }
-    try {
-        const SerializationNode& velocitiesNode = node.getChildNode("Velocities");
-        for (int i = 0; i < (int) velocitiesNode.getChildren().size(); i++) {
-            const SerializationNode& particle = velocitiesNode.getChildren()[i];
-            outVelocities.push_back(Vec3(particle.getDoubleProperty("x"),particle.getDoubleProperty("y"),particle.getDoubleProperty("z")));
+        else if (child.getName() == "Positions") {
+            vector<Vec3> outPositions;
+            for (int i = 0; i < (int) child.getChildren().size(); i++) {
+                const SerializationNode& particle = child.getChildren()[i];
+                outPositions.push_back(Vec3(particle.getDoubleProperty("x"),particle.getDoubleProperty("y"),particle.getDoubleProperty("z")));
+            }
+            builder.setPositions(outPositions);
+            arraySizes.push_back(outPositions.size());
         }
-        types = types | State::Velocities;
-    } catch (const OpenMMException &) {
-        // do nothing    
-    }
-    try {
-        const SerializationNode& forcesNode = node.getChildNode("Forces");
-        for (int i = 0; i < (int) forcesNode.getChildren().size(); i++) {
-            const SerializationNode& particle = forcesNode.getChildren()[i];
-            outForces.push_back(Vec3(particle.getDoubleProperty("x"),particle.getDoubleProperty("y"),particle.getDoubleProperty("z")));
+        else if (child.getName() == "Velocities") {
+            vector<Vec3> outVelocities;
+            for (int i = 0; i < (int) child.getChildren().size(); i++) {
+                const SerializationNode& particle = child.getChildren()[i];
+                outVelocities.push_back(Vec3(particle.getDoubleProperty("x"),particle.getDoubleProperty("y"),particle.getDoubleProperty("z")));
+            }
+            builder.setVelocities(outVelocities);
+            arraySizes.push_back(outVelocities.size());
+        }
+        else if (child.getName() == "Forces") {
+            vector<Vec3> outForces;
+            for (int i = 0; i < (int) child.getChildren().size(); i++) {
+                const SerializationNode& particle = child.getChildren()[i];
+                outForces.push_back(Vec3(particle.getDoubleProperty("x"),particle.getDoubleProperty("y"),particle.getDoubleProperty("z")));
+            }
+            builder.setForces(outForces);
+            arraySizes.push_back(outForces.size());
         }
-        types = types | State::Forces;
-    } catch (const OpenMMException &) {
-        // do nothing
-    }
-    vector<int> arraySizes;
-    State::StateBuilder builder(outTime);
-    if (types & State::Positions) {
-        builder.setPositions(outPositions);
-        arraySizes.push_back(outPositions.size());
-    }  
-    if (types & State::Velocities) {
-        builder.setVelocities(outVelocities);
-        arraySizes.push_back(outVelocities.size());
-    }
-    if (types & State::Forces) {
-        builder.setForces(outForces);
-        arraySizes.push_back(outForces.size());
-    }
-    if (types & State::Energy) {
-        builder.setEnergy(kineticEnergy, potentialEnergy);
-    }
-    if (types & State::Parameters) {
-        builder.setParameters(outStateParams);
     }
-
     for (int i = 1; i < arraySizes.size(); i++) {
         if (arraySizes[i] != arraySizes[i-1]) {
             throw(OpenMMException("State Deserialization Particle Size Mismatch, check number of particles in Forces, Velocities, Positions!"));

serialization/tests/TestSerializeState.cpp

@@ -66,8 +66,8 @@ void testSerialization() {
     system.addForce(nonbonded);
 	system.addForce(new AndersenThermostat(393.3, 19.3));
 	system.addForce(new MonteCarloBarostat(25, 393.3, 25));
-	Integrator *intg = new LangevinIntegrator(300,79,0.002);
-	Context *ctxt = new Context(system, *intg);
+	LangevinIntegrator intg(300,79,0.002);
+	Context context(system, intg);
 	
 	// Set positions, velocities, forces
 	vector<Vec3> positions;
@@ -79,11 +79,11 @@ void testSerialization() {
 		velocities.push_back(Vec3( ((float) rand()/(float) RAND_MAX)*6.2, ((float) rand()/(float) RAND_MAX)*6.2, ((float) rand()/(float) RAND_MAX)*6.2));
 	}
 
-	ctxt->setPositions(positions);
-	ctxt->setVelocities(velocities);
+	context.setPositions(positions);
+	context.setVelocities(velocities);
 
 	// Serialize and then deserialize it.
-	State s1 = ctxt->getState(State::Positions | State::Velocities | State::Forces | State::Energy | State::Parameters);
+	State s1 = context.getState(State::Positions | State::Velocities | State::Forces | State::Energy | State::Parameters);
 
 	stringstream buffer;
     XmlSerializer::serialize<State>(&s1, "State", buffer);
@@ -132,6 +132,55 @@ void testSerialization() {
 		assert((it1->first).compare(it2->first) == 0);
 		ASSERT_EQUAL(it1->second, it2->second);
 	}
+    delete copy;
+
+    // Now create a series of States that include only one type of information.  Verify
+    // that serialization works correctly for them.
+
+    for (int types = 1; types <= 16; types *= 2) {
+        State s3 = context.getState(types);
+        stringstream buffer2;
+        XmlSerializer::serialize<State>(&s3, "State", buffer2);
+        copy = XmlSerializer::deserialize<State>(buffer2);
+        int foundTypes = 0;
+        try {
+            copy->getPositions();
+            foundTypes += State::Positions;
+        }
+        catch (...) {
+            // Ignore
+        }
+        try {
+            copy->getVelocities();
+            foundTypes += State::Velocities;
+        }
+        catch (...) {
+            // Ignore
+        }
+        try {
+            copy->getForces();
+            foundTypes += State::Forces;
+        }
+        catch (...) {
+            // Ignore
+        }
+        try {
+            copy->getPotentialEnergy();
+            foundTypes += State::Energy;
+        }
+        catch (...) {
+            // Ignore
+        }
+        try {
+            copy->getParameters();
+            foundTypes += State::Parameters;
+        }
+        catch (...) {
+            // Ignore
+        }
+        delete copy;
+        ASSERT_EQUAL(types, foundTypes);
+    }
 }
 
 int main() {

tests/TestSystem.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2015 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/HarmonicAngleForce.h"
+#include "openmm/HarmonicBondForce.h"
+#include "openmm/System.h"
+#include "openmm/VirtualSite.h"
+#include <iostream>
+
+using namespace OpenMM;
+using namespace std;
+
+/**
+ * Test the methods for manipulating System objects.
+ */
+void testCreateSystem() {
+    int numParticles = 10;
+    System system;
+    
+    // Test adding and modifying particles.
+    
+    for (int i = 0; i < numParticles; i++)
+        system.addParticle(1.0+0.1*i);
+    system.setParticleMass(5, 100.0);
+    ASSERT_EQUAL(numParticles, system.getNumParticles());
+    for (int i = 0; i < numParticles; i++) {
+        double mass = (i == 5 ? 100.0 : 1.0+0.1*i);
+        ASSERT_EQUAL(mass, system.getParticleMass(i));
+    }
+    
+    // Test adding, removing, and modifying constraints.
+    
+    for (int i = 0; i < numParticles-1; i++)
+        system.addConstraint(i, i+1, 0.2*i);
+    system.removeConstraint(5);
+    system.setConstraintParameters(3, 0, 5, 99.0);
+    ASSERT_EQUAL(numParticles-2, system.getNumConstraints());
+    for (int i = 0; i < numParticles-2; i++) {
+        int p1, p2;
+        double dist;
+        system.getConstraintParameters(i, p1, p2, dist);
+        if (i == 3) {
+            ASSERT_EQUAL(0, p1);
+            ASSERT_EQUAL(5, p2);
+            ASSERT_EQUAL(99.0, dist);
+        }
+        else {
+            int j = (i < 5 ? i : i+1);
+            ASSERT_EQUAL(j, p1);
+            ASSERT_EQUAL(j+1, p2);
+            ASSERT_EQUAL(0.2*j, dist);
+        }
+    }
+    
+    // Test adding and removing forces.
+    
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    system.addForce(bonds);
+    HarmonicAngleForce* angles = new HarmonicAngleForce();
+    system.addForce(angles);
+    ASSERT_EQUAL(2, system.getNumForces());
+    ASSERT_EQUAL(bonds, &system.getForce(0));
+    ASSERT_EQUAL(angles, &system.getForce(1));
+    system.removeForce(0);
+    ASSERT_EQUAL(1, system.getNumForces());
+    ASSERT_EQUAL(angles, &system.getForce(0));
+    
+    // Test adding and removing virtual sites.
+    
+    for (int i = 0; i < numParticles; i++) {
+        ASSERT_EQUAL(false, system.isVirtualSite(i));
+    }
+    TwoParticleAverageSite* site = new TwoParticleAverageSite(2, 3, 0.4, 0.6);
+    system.setVirtualSite(4, site);
+    for (int i = 0; i < numParticles; i++) {
+        ASSERT_EQUAL(i == 4, system.isVirtualSite(i));
+    }
+    ASSERT_EQUAL(site, &system.getVirtualSite(4));
+    system.setVirtualSite(4, NULL);
+    for (int i = 0; i < numParticles; i++) {
+        ASSERT_EQUAL(false, system.isVirtualSite(i));
+    }
+}
+
+int main() {
+    try {
+        testCreateSystem();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}

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

@@ -777,11 +777,12 @@ class CharmmPsfFile(object):
         # Add each chain (separate 'system's) and residue
         for atom in self.atom_list:
             if atom.system != last_chain:
-                chain = topology.addChain()
+                chain = topology.addChain(atom.system)
                 last_chain = atom.system
             if atom.residue.idx != last_residue:
                 last_residue = atom.residue.idx
-                residue = topology.addResidue(atom.residue.resname, chain)
+                residue = topology.addResidue(atom.residue.resname, chain,
+                                              str(atom.residue.idx))
             if atom.type is not None:
                 # This is the most reliable way of determining the element
                 atomic_num = atom.type.atomic_number

wrappers/python/simtk/openmm/app/data/amoeba2013.xml

 <ForceField>
  <Info>
-  <Source>amoebapro13_5-2014.prm</Source>
-  <DateGenerated>2014-05-08</DateGenerated>
+  <Source>amoebapro13.prm</Source>
+  <DateGenerated>2015-02-19</DateGenerated>
   <Reference>Yue Shi, Zhen Xia, Jiajing Zhang, Robert Best, Chuanjie Wu, Jay W. Ponder, and Pengyu Ren. Polarizable Atomic Multipole-Based AMOEBA Force Field for Proteins. Journal of Chemical Theory and Computation, 9(9):4046–4063, 2013.</Reference>
  </Info>
  <AtomTypes>
@@ -170,7 +170,7 @@
   <Type name="163" class="32" element="C" mass="12.011"/>
   <Type name="164" class="33" element="O" mass="15.999"/>
   <Type name="165" class="34" element="O" mass="15.999"/>
-  <Type name="166" class="35" element="H" mass="15.999"/>
+  <Type name="166" class="35" element="H" mass="1.008"/>
   <Type name="167" class="8" element="C" mass="12.011"/>
   <Type name="168" class="9" element="H" mass="1.008"/>
   <Type name="169" class="8" element="C" mass="12.011"/>

wrappers/python/simtk/openmm/app/data/amoeba2013_gk.xml

 <ForceField>
  <Info>
-  <Source>amoebapro13_5-2014.prm</Source>
-  <DateGenerated>2014-05-08</DateGenerated>
+  <Source>amoebapro13.prm</Source>
+  <DateGenerated>2015-02-19</DateGenerated>
   <Reference>Yue Shi, Zhen Xia, Jiajing Zhang, Robert Best, Chuanjie Wu, Jay W. Ponder, and Pengyu Ren. Polarizable Atomic Multipole-Based AMOEBA Force Field for Proteins. Journal of Chemical Theory and Computation, 9(9):4046–4063, 2013.</Reference>
  </Info>
  <AmoebaGeneralizedKirkwoodForce solventDielectric="78.3" soluteDielectric="1.0" includeCavityTerm="1" probeRadius="0.14" surfaceAreaFactor="-170.351730663">

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

@@ -34,6 +34,7 @@ __version__ = "1.0"
 import os
 import itertools
 import xml.etree.ElementTree as etree
+import math
 from math import sqrt, cos
 import simtk.openmm as mm
 import simtk.unit as unit
@@ -2101,7 +2102,7 @@ class AmoebaAngleGenerator:
                     hit = 1
                     if isConstrained and self.k[i] != 0.0:
                         angleDict['idealAngle'] = self.angle[i][0]
-                        addAngleConstraint(angle, self.angle[i][0], data, sys)
+                        addAngleConstraint(angle, self.angle[i][0]*math.pi/180.0, data, sys)
                     elif self.k[i] != 0:
                         lenAngle = len(self.angle[i])
                         if (lenAngle > 1):
@@ -2179,7 +2180,7 @@ class AmoebaAngleGenerator:
                     hit = 1
                     angleDict['idealAngle'] = self.angle[i][0]
                     if (isConstrained and self.k[i] != 0.0):
-                        addAngleConstraint(angle, self.angle[i][0], data, sys)
+                        addAngleConstraint(angle, self.angle[i][0]*math.pi/180.0, data, sys)
                     else:
                         force.addAngle(angle[0], angle[1], angle[2], angle[3], self.angle[i][0], self.k[i])
                     break

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

@@ -135,6 +135,8 @@ class PrmtopLoader(object):
                     format = format[index0+1:index1]
                     m = FORMAT_RE_PATTERN.search(format)
                     self._raw_format[self._flags[-1]] = (format, m.group(1), m.group(2), m.group(3), m.group(4))
+                elif line.startswith('%COMMENT'):
+                    continue
                 elif self._flags \
                      and 'TITLE'==self._flags[-1] \
                      and not self._raw_data['TITLE']:

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

@@ -150,6 +150,8 @@ class PdbStructure(object):
         self._reset_residue_numbers()
         # Read one line at a time
         for pdb_line in input_stream:
+            if not isinstance(pdb_line, str):
+                pdb_line = pdb_line.decode('utf-8')
             # Look for atoms
             if (pdb_line.find("ATOM  ") == 0) or (pdb_line.find("HETATM") == 0):
                 self._add_atom(Atom(pdb_line, self))

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

@@ -79,15 +79,33 @@ def computePeriodicBoxVectors(a_length, b_length, c_length, alpha, beta, gamma):
     b = b - a*round(b[0]/a[0])
     return (a, b, c)*nanometers
 
+def reducePeriodicBoxVectors(periodicBoxVectors):
+    """ Reduces the representation of the PBC. periodicBoxVectors is expected to
+    be an unpackable iterable of length-3 iterables
+    """
+    if is_quantity(periodicBoxVectors):
+        a, b, c = periodicBoxVectors.value_in_unit(nanometers)
+    else:
+        a, b, c = periodicBoxVectors
+    a = Vec3(*a)
+    b = Vec3(*b)
+    c = Vec3(*c)
+
+    c = c - b*round(c[1]/b[1])
+    c = c - a*round(c[0]/a[0])
+    b = b - a*round(b[0]/a[0])
+
+    return (a, b, c) * nanometers
+
 def computeLengthsAndAngles(periodicBoxVectors):
     """Convert periodic box vectors to lengths and angles.
 
     Lengths are returned in nanometers and angles in radians.
     """
     if is_quantity(periodicBoxVectors):
-        (a, b, c) = vectors.value_in_unit(nanometers)
+        (a, b, c) = periodicBoxVectors.value_in_unit(nanometers)
     else:
-        a, b, c = vectors
+        a, b, c = periodicBoxVectors
     a_length = norm(a)
     b_length = norm(b)
     c_length = norm(c)

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

@@ -924,7 +924,7 @@ class Modeller(object):
                     # extra points.
 
                     template = None
-                    residueNoEP = Residue(residue.name, residue.index, residue.chain)
+                    residueNoEP = Residue(residue.name, residue.index, residue.chain, residue.id)
                     residueNoEP._atoms = [atom for atom in residue.atoms() if atom.element is not None]
                     if signature in forcefield._templateSignatures:
                         for t in forcefield._templateSignatures[signature]:

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

@@ -229,16 +229,19 @@ class PDBFile(object):
                 map[atom.attrib[id]] = name
 
     @staticmethod
-    def writeFile(topology, positions, file=sys.stdout, modelIndex=None):
+    def writeFile(topology, positions, file=sys.stdout, keepIds=False):
         """Write a PDB file containing a single model.
 
         Parameters:
          - topology (Topology) The Topology defining the model to write
          - positions (list) The list of atomic positions to write
          - file (file=stdout) A file to write to
+         - keepIds (bool=False) If True, keep the residue and chain IDs specified in the Topology rather than generating
+           new ones.  Warning: It is up to the caller to make sure these are valid IDs that satisfy the requirements of
+           the PDB format.  Otherwise, the output file will be invalid.
         """
         PDBFile.writeHeader(topology, file)
-        PDBFile.writeModel(topology, positions, file)
+        PDBFile.writeModel(topology, positions, file, keepIds=keepIds)
         PDBFile.writeFooter(topology, file)
 
     @staticmethod
@@ -262,7 +265,7 @@ class PDBFile(object):
             print >>file, "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f P 1           1 " % (a_length, b_length, c_length, alpha, beta, gamma)
 
     @staticmethod
-    def writeModel(topology, positions, file=sys.stdout, modelIndex=None):
+    def writeModel(topology, positions, file=sys.stdout, modelIndex=None, keepIds=False):
         """Write out a model to a PDB file.
 
         Parameters:
@@ -270,6 +273,9 @@ class PDBFile(object):
          - positions (list) The list of atomic positions to write
          - file (file=stdout) A file to write the model to
          - modelIndex (int=None) If not None, the model will be surrounded by MODEL/ENDMDL records with this index
+         - keepIds (bool=False) If True, keep the residue and chain IDs specified in the Topology rather than generating
+           new ones.  Warning: It is up to the caller to make sure these are valid IDs that satisfy the requirements of
+           the PDB format.  Otherwise, the output file will be invalid.
         """
         if len(list(topology.atoms())) != len(positions):
             raise ValueError('The number of positions must match the number of atoms')
@@ -284,13 +290,20 @@ class PDBFile(object):
         if modelIndex is not None:
             print >>file, "MODEL     %4d" % modelIndex
         for (chainIndex, chain) in enumerate(topology.chains()):
-            chainName = chr(ord('A')+chainIndex%26)
+            if keepIds:
+                chainName = chain.id
+            else:
+                chainName = chr(ord('A')+chainIndex%26)
             residues = list(chain.residues())
             for (resIndex, res) in enumerate(residues):
                 if len(res.name) > 3:
                     resName = res.name[:3]
                 else:
                     resName = res.name
+                if keepIds:
+                    resId = res.id
+                else:
+                    resId = "%4d" % ((resIndex+1)%10000)
                 for atom in res.atoms():
                     if len(atom.name) < 4 and atom.name[:1].isalpha() and (atom.element is None or len(atom.element.symbol) < 2):
                         atomName = ' '+atom.name
@@ -303,16 +316,15 @@ class PDBFile(object):
                         symbol = atom.element.symbol
                     else:
                         symbol = ' '
-                    line = "ATOM  %5d %-4s %3s %s%4d    %s%s%s  1.00  0.00          %2s  " % (
-                        atomIndex%100000, atomName, resName, chainName,
-                        (resIndex+1)%10000, _format_83(coords[0]),
+                    line = "ATOM  %5d %-4s %3s %s%4s    %s%s%s  1.00  0.00          %2s  " % (
+                        atomIndex%100000, atomName, resName, chainName, resId, _format_83(coords[0]),
                         _format_83(coords[1]), _format_83(coords[2]), symbol)
                     assert len(line) == 80, 'Fixed width overflow detected'
                     print >>file, line
                     posIndex += 1
                     atomIndex += 1
                 if resIndex == len(residues)-1:
-                    print >>file, "TER   %5d      %3s %s%4d" % (atomIndex, resName, chainName, resIndex+1)
+                    print >>file, "TER   %5d      %3s %s%4s" % (atomIndex, resName, chainName, resId)
                     atomIndex += 1
         if modelIndex is not None:
             print >>file, "ENDMDL"