Merge branch 'master' of https://github.com/peastman/openmm

.travis.yml

-language: c
+language: cpp
+compiler:
+  - clang
 
-install:
-  - source tools/ci/install.sh
-  - export PYTHONUNBUFFERED=true
+before_install:
+  - sudo apt-get update -qq
+  - sudo apt-get install -qq libpcre3 libpcre3-dev gromacs
+  - sudo apt-get install -qq swig doxygen llvm-3.3
+  - sudo apt-get install -qq python-numpy python-scipy python-pip
+  - sudo pip install nose
+  - export ASAN_SYMBOLIZER_PATH=/usr/bin/llvm-symbolizer-3.3
 
 script:
-  - export CC="clang++"
-  - source deactivate
-  - conda install --yes conda-build
-  - # Build the conda package, testing build before packaging.
-  - conda build tools/conda-recipe
-  - # Install the conda package locally.
-  - source activate $python
-  - conda install $HOME/miniconda/conda-bld/linux-64/openmm-dev-*
-  - conda list -e
-  - # Run the Python tests.
-  - pushd .
-  - cd wrappers/python/tests
+  - cmake -DCMAKE_INSTALL_PREFIX=$HOME/OpenMM .
+  - make -j2
+  - make -j2 install
+  - sudo make PythonInstall
+  - # run all of the tests
+  - ctest -j2 -V
+  - # get a list of all of the failed tests into this stupid ctest format
+  - python -c 'fn = "Testing/Temporary/LastTestsFailed.log"; import os; os.path.exists(fn) or exit(0); l = [line.split(":")[0] for line in open(fn)]; triplets = zip(l, l, [","]*len(l)); print "".join(",".join(t) for t in triplets)' > FailedTests.log
+  - # rerun all of the failed tests
+  - if [ -s FailedTests.log ]; then ctest -V -I FailedTests.log; fi;
+  - # run the python tests too
+  - cd python/tests
   - nosetests -vv --processes=-1 --process-timeout=200
-  - popd
-
-env:
-  global:
-    # encrypted BINSTAR_TOKEN for push of dev package to binstar
-    - secure: Qz3pEYXXFnNQ/WK+15ad4cdbLJvzgCIZRwKD9fLiS3CDO2ldAQWxzaz8RQOwqbFtZUWu7lQpr+GukNJz5p0w18QEto+BxLYG9aW5mjoc+F2vCjyWFjkwnJ/Z/3uBKTcr5x9Y7HKaPGivaJ4BNACifjt7cCpeVJzV6u2+bBgSoHc=
-
-  matrix:
-    - python=2.7  CONDA_PY=27
-    #- python=3.3  CONDA_PY=33
-
-after_success:
-  - echo "after_success"
-  - source tools/ci/after_success.sh
-
-

openmmapi/include/openmm/CustomNonbondedForce.h

@@ -157,6 +157,7 @@ public:
      *                  of r, the distance between them, as well as any global and per-particle parameters
      */
     explicit CustomNonbondedForce(const std::string& energy);
+    CustomNonbondedForce(const CustomNonbondedForce& rhs); // copy constructor
     ~CustomNonbondedForce();
     /**
      * Get the number of particles for which force field parameters have been defined.
@@ -466,6 +467,7 @@ public:
 protected:
     ForceImpl* createImpl() const;
 private:
+    // REMEMBER TO UPDATE THE COPY CONSTRUCTOR IF YOU ADD ANY NEW FIELDS !!
     class ParticleInfo;
     class PerParticleParameterInfo;
     class GlobalParameterInfo;

openmmapi/include/openmm/TabulatedFunction.h

@@ -59,6 +59,7 @@ class OPENMM_EXPORT TabulatedFunction {
 public:
     virtual ~TabulatedFunction() {
     }
+    virtual TabulatedFunction* Copy() const = 0;
 };
 
 /**
@@ -96,6 +97,10 @@ public:
      * @param max            the value of x corresponding to the last element of values
      */
     void setFunctionParameters(const std::vector<double>& values, double min, double max);
+    /**
+     * Create a deep copy of the tabulated function.
+     */
+    Continuous1DFunction* Copy() const;
 private:
     std::vector<double> values;
     double min, max;
@@ -151,6 +156,10 @@ public:
      * @param ymax       the value of y corresponding to the last element of values
      */
     void setFunctionParameters(int xsize, int ysize, const std::vector<double>& values, double xmin, double xmax, double ymin, double ymax);
+    /**
+     * Create a deep copy of the tabulated function
+     */
+    Continuous2DFunction* Copy() const;
 private:
     std::vector<double> values;
     int xsize, ysize;
@@ -222,6 +231,10 @@ public:
      * @param zmax       the value of z corresponding to the last element of values
      */
     void setFunctionParameters(int xsize, int ysize, int zsize, const std::vector<double>& values, double xmin, double xmax, double ymin, double ymax, double zmin, double zmax);
+    /**
+     * Create a deep copy of the tabulated function
+     */
+    Continuous3DFunction* Copy() const;
 private:
     std::vector<double> values;
     int xsize, ysize, zsize;
@@ -253,6 +266,10 @@ public:
      * @param values         the tabulated values of the function f(x)
      */
     void setFunctionParameters(const std::vector<double>& values);
+    /**
+     * Create a deep copy of the tabulated function
+     */
+    Discrete1DFunction* Copy() const;
 private:
     std::vector<double> values;
 };
@@ -291,6 +308,10 @@ public:
      *                  values[i+xsize*j] = f(i,j).  This must be of length xsize*ysize.
      */
     void setFunctionParameters(int xsize, int ysize, const std::vector<double>& values);
+    /**
+     * Create a deep copy of the tabulated function
+     */
+    Discrete2DFunction* Copy() const;
 private:
     int xsize, ysize;
     std::vector<double> values;
@@ -333,6 +354,10 @@ public:
      *                  values[i+xsize*j+xsize*ysize*k] = f(i,j,k).  This must be of length xsize*ysize*zsize.
      */
     void setFunctionParameters(int xsize, int ysize, int zsize, const std::vector<double>& values);
+    /**
+     * Create a deep copy of the tabulated function
+     */
+    Discrete3DFunction* Copy() const;
 private:
     int xsize, ysize, zsize;
     std::vector<double> values;

openmmapi/src/CustomNonbondedForce.cpp

@@ -51,6 +51,23 @@ CustomNonbondedForce::CustomNonbondedForce(const string& energy) : energyExpress
     switchingDistance(-1.0), useSwitchingFunction(false), useLongRangeCorrection(false) {
 }
 
+CustomNonbondedForce::CustomNonbondedForce(const CustomNonbondedForce& rhs) {
+    // Copy everything and deep copy the tabulated functions
+    energyExpression = rhs.energyExpression;
+    nonbondedMethod = rhs.nonbondedMethod;
+    cutoffDistance = rhs.cutoffDistance;
+    switchingDistance = rhs.switchingDistance;
+    useSwitchingFunction = rhs.useSwitchingFunction;
+    useLongRangeCorrection = rhs.useLongRangeCorrection;
+    parameters = rhs.parameters;
+    globalParameters = rhs.globalParameters;
+    particles = rhs.particles;
+    exclusions = rhs.exclusions;
+    interactionGroups = rhs.interactionGroups;
+    for (vector<FunctionInfo>::const_iterator it = rhs.functions.begin(); it != rhs.functions.end(); it++)
+        functions.push_back(FunctionInfo(it->name, it->function->Copy()));
+}
+
 CustomNonbondedForce::~CustomNonbondedForce() {
     for (int i = 0; i < (int) functions.size(); i++)
         delete functions[i].function;

openmmapi/src/TabulatedFunction.cpp

@@ -61,6 +61,13 @@ void Continuous1DFunction::setFunctionParameters(const vector<double>& values, d
     this->max = max;
 }
 
+Continuous1DFunction* Continuous1DFunction::Copy() const {
+    vector<double> new_vec(values.size());
+    for (size_t i = 0; i < values.size(); i++)
+        new_vec[i] = values[i];
+    return new Continuous1DFunction(new_vec, min, max);
+}
+
 Continuous2DFunction::Continuous2DFunction(int xsize, int ysize, const vector<double>& values, double xmin, double xmax, double ymin, double ymax) {
     if (xsize < 2 || ysize < 2)
         throw OpenMMException("Continuous2DFunction: must have at least two points along each axis");
@@ -107,6 +114,13 @@ void Continuous2DFunction::setFunctionParameters(int xsize, int ysize, const vec
     this->ymax = ymax;
 }
 
+Continuous2DFunction* Continuous2DFunction::Copy() const {
+    vector<double> new_vec(values.size());
+    for (size_t i = 0; i < values.size(); i++)
+        new_vec[i] = values[i];
+    return new Continuous2DFunction(xsize, ysize, new_vec, xmin, xmax, ymin, ymax);
+}
+
 Continuous3DFunction::Continuous3DFunction(int xsize, int ysize, int zsize, const vector<double>& values, double xmin, double xmax, double ymin, double ymax, double zmin, double zmax) {
     if (xsize < 2 || ysize < 2 || zsize < 2)
         throw OpenMMException("Continuous3DFunction: must have at least two points along each axis");
@@ -166,6 +180,14 @@ void Continuous3DFunction::setFunctionParameters(int xsize, int ysize, int zsize
     this->zmax = zmax;
 }
 
+Continuous3DFunction* Continuous3DFunction::Copy() const {
+    vector<double> new_vec(values.size());
+    for (size_t i = 0; i < values.size(); i++)
+        new_vec[i] = values[i];
+    return new Continuous3DFunction(xsize, ysize, zsize, new_vec, xmin, xmax, ymin, ymax, zmin, zmax);
+}
+
+
 Discrete1DFunction::Discrete1DFunction(const vector<double>& values) {
     this->values = values;
 }
@@ -178,6 +200,13 @@ void Discrete1DFunction::setFunctionParameters(const vector<double>& values) {
     this->values = values;
 }
 
+Discrete1DFunction* Discrete1DFunction::Copy() const {
+    vector<double> new_vec(values.size());
+    for (size_t i = 0; i < values.size(); i++)
+        new_vec[i] = values[i];
+    return new Discrete1DFunction(new_vec);
+}
+
 Discrete2DFunction::Discrete2DFunction(int xsize, int ysize, const vector<double>& values) {
     if (values.size() != xsize*ysize)
         throw OpenMMException("Discrete2DFunction: incorrect number of values");
@@ -200,6 +229,13 @@ void Discrete2DFunction::setFunctionParameters(int xsize, int ysize, const vecto
     this->values = values;
 }
 
+Discrete2DFunction* Discrete2DFunction::Copy() const {
+    vector<double> new_vec(values.size());
+    for (size_t i = 0; i < values.size(); i++)
+        new_vec[i] = values[i];
+    return new Discrete2DFunction(xsize, ysize, new_vec);
+}
+
 Discrete3DFunction::Discrete3DFunction(int xsize, int ysize, int zsize, const vector<double>& values) {
     if (values.size() != xsize*ysize*zsize)
         throw OpenMMException("Discrete3DFunction: incorrect number of values");
@@ -224,3 +260,10 @@ void Discrete3DFunction::setFunctionParameters(int xsize, int ysize, int zsize,
     this->zsize = zsize;
     this->values = values;
 }
+
+Discrete3DFunction* Discrete3DFunction::Copy() const {
+    vector<double> new_vec(values.size());
+    for (size_t i = 0; i < values.size(); i++)
+        new_vec[i] = values[i];
+    return new Discrete3DFunction(xsize, ysize, zsize, new_vec);
+}

wrappers/python/CMakeLists.txt

@@ -55,6 +55,7 @@ foreach(SUBDIR ${SUBDIRS})
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.xml"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.pdb"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.prmtop"
+        "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.inpcrd"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.parm7"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.rst7"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.ncrst"

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

@@ -100,7 +100,7 @@ class ForceField(object):
         """Load one or more XML files and create a ForceField object based on them.
 
         Parameters:
-         - files A list of XML files defining the force field.  Each entry may
+         - files (list) A list of XML files defining the force field.  Each entry may
            be an absolute file path, a path relative to the current working
            directory, a path relative to this module's data subdirectory
            (for built in force fields), or an open file-like object with a
@@ -113,47 +113,59 @@ class ForceField(object):
         self._forces = []
         self._scripts = []
         for file in files:
-            try:
-                # this handles either filenames or open file-like objects
-                tree = etree.parse(file)
-            except IOError:
-                tree = etree.parse(os.path.join(os.path.dirname(__file__), 'data', file))
-            root = tree.getroot()
-
-            # Load the atom types.
-
-            if tree.getroot().find('AtomTypes') is not None:
-                for type in tree.getroot().find('AtomTypes').findall('Type'):
-                    self.registerAtomType(type.attrib)
-
-            # Load the residue templates.
-
-            if tree.getroot().find('Residues') is not None:
-                for residue in root.find('Residues').findall('Residue'):
-                    resName = residue.attrib['name']
-                    template = ForceField._TemplateData(resName)
-                    for atom in residue.findall('Atom'):
-                        template.atoms.append(ForceField._TemplateAtomData(atom.attrib['name'], atom.attrib['type'], self._atomTypes[atom.attrib['type']][2]))
-                    for site in residue.findall('VirtualSite'):
-                        template.virtualSites.append(ForceField._VirtualSiteData(site))
-                    for bond in residue.findall('Bond'):
-                        template.addBond(int(bond.attrib['from']), int(bond.attrib['to']))
-                    for bond in residue.findall('ExternalBond'):
-                        b = int(bond.attrib['from'])
-                        template.externalBonds.append(b)
-                        template.atoms[b].externalBonds += 1
-                    self.registerResidueTemplate(template)
-
-            # Load force definitions
-
-            for child in root:
-                if child.tag in parsers:
-                    parsers[child.tag](child, self)
-
-            # Load scripts
-
-            for node in tree.getroot().findall('Script'):
-                self.registerScript(node.text)
+            self.loadFile(file)
+    
+    def loadFile(self, file):
+        """Load an XML file and add the definitions from it to this FieldField.
+        
+        Parameters:
+         - file (string or file) An XML file containing force field definitions.  It may
+           be either an absolute file path, a path relative to the current working
+           directory, a path relative to this module's data subdirectory
+           (for built in force fields), or an open file-like object with a
+           read() method from which the forcefield XML data can be loaded.
+        """
+        try:
+            # this handles either filenames or open file-like objects
+            tree = etree.parse(file)
+        except IOError:
+            tree = etree.parse(os.path.join(os.path.dirname(__file__), 'data', file))
+        root = tree.getroot()
+
+        # Load the atom types.
+
+        if tree.getroot().find('AtomTypes') is not None:
+            for type in tree.getroot().find('AtomTypes').findall('Type'):
+                self.registerAtomType(type.attrib)
+
+        # Load the residue templates.
+
+        if tree.getroot().find('Residues') is not None:
+            for residue in root.find('Residues').findall('Residue'):
+                resName = residue.attrib['name']
+                template = ForceField._TemplateData(resName)
+                for atom in residue.findall('Atom'):
+                    template.atoms.append(ForceField._TemplateAtomData(atom.attrib['name'], atom.attrib['type'], self._atomTypes[atom.attrib['type']][2]))
+                for site in residue.findall('VirtualSite'):
+                    template.virtualSites.append(ForceField._VirtualSiteData(site))
+                for bond in residue.findall('Bond'):
+                    template.addBond(int(bond.attrib['from']), int(bond.attrib['to']))
+                for bond in residue.findall('ExternalBond'):
+                    b = int(bond.attrib['from'])
+                    template.externalBonds.append(b)
+                    template.atoms[b].externalBonds += 1
+                self.registerResidueTemplate(template)
+
+        # Load force definitions
+
+        for child in root:
+            if child.tag in parsers:
+                parsers[child.tag](child, self)
+
+        # Load scripts
+
+        for node in tree.getroot().findall('Script'):
+            self.registerScript(node.text)
 
     def getGenerators(self):
         """Get the list of all registered generators."""

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
+from simtk.openmm.vec3 import Vec3
 import customgbforces as customgb
 
 #=============================================================================================
@@ -116,6 +117,7 @@ class PrmtopLoader(object):
         self._flags=[]
         self._raw_format={}
         self._raw_data={}
+        self._has_nbfix_terms = False
 
         fIn=open(inFilename)
         for line in fIn:
@@ -218,28 +220,16 @@ class PrmtopLoader(object):
         try:
             return self._massList
         except AttributeError:
-            pass
-
-        self._massList=[]
-        raw_masses=self._raw_data['MASS']
-        for ii in range(self.getNumAtoms()):
-            self._massList.append(float(raw_masses[ii]))
-        self._massList = self._massList
-        return self._massList
+            self._massList = [float(x) for x in self._raw_data['MASS']]
+            return self._massList
 
     def getCharges(self):
         """Return a list of atomic charges in the system"""
         try:
             return self._chargeList
         except AttributeError:
-            pass
-
-        self._chargeList=[]
-        raw_charges=self._raw_data['CHARGE']
-        for ii in range(self.getNumAtoms()):
-            self._chargeList.append(float(raw_charges[ii])/18.2223)
-        self._chargeList = self._chargeList
-        return self._chargeList
+            self._chargeList = [float(x)/18.2223 for x in self._raw_data['CHARGE']]
+            return self._chargeList
 
     def getAtomName(self, iAtom):
         """Return the atom name for iAtom"""
@@ -254,11 +244,8 @@ class PrmtopLoader(object):
         try:
             return self._atomTypeIndexes
         except AttributeError:
-            pass
-        self._atomTypeIndexes=[]
-        for atomTypeIndex in  self._raw_data['ATOM_TYPE_INDEX']:
-            self._atomTypeIndexes.append(int(atomTypeIndex))
-        return self._atomTypeIndexes
+            self._atomTypeIndexes = [int(x) for x in self._raw_data['ATOM_TYPE_INDEX']]
+            return self._atomTypeIndexes
 
     def getAtomType(self, iAtom):
         """Return the AMBER atom type for iAtom"""
@@ -275,11 +262,11 @@ class PrmtopLoader(object):
 
     def getResidueLabel(self, iAtom=None, iRes=None):
         """Return residue label for iAtom OR iRes"""
-        if iRes==None and iAtom==None:
+        if iRes is None and iAtom is None:
             raise Exception("only specify iRes or iAtom, not both")
-        if iRes!=None and iAtom!=None:
+        if iRes is not None and iAtom is not None:
             raise Exception("iRes or iAtom must be set")
-        if iRes!=None:
+        if iRes is not None:
             return self._raw_data['RESIDUE_LABEL'][iRes]
         else:
             return self.getResidueLabel(iRes=self._getResiduePointer(iAtom))
@@ -306,6 +293,9 @@ class PrmtopLoader(object):
         elements of the Lennard-Jones A and B coefficient matrices are found,
         NbfixPresent exception is raised
         """
+        if self._has_nbfix_terms:
+            raise NbfixPresent('Off-diagonal Lennard-Jones elements found. '
+                        'Cannot determine LJ parameters for individual atoms.')
         try:
             return self._nonbondTerms
         except AttributeError:
@@ -344,13 +334,13 @@ class PrmtopLoader(object):
                 b = float(self._raw_data['LENNARD_JONES_BCOEF'][index])
                 if a == 0 or b == 0:
                     if a != 0 or b != 0 or (wdij != 0 and rij != 0):
-                        del self._nonbondTerms
+                        self._has_nbfix_terms = True
                         raise NbfixPresent('Off-diagonal Lennard-Jones elements'
                                            ' found. Cannot determine LJ '
                                            'parameters for individual atoms.')
                 elif (abs((a - (wdij * rij ** 12)) / a) > 1e-6 or
                       abs((b - (2 * wdij * rij**6)) / b) > 1e-6):
-                    del self._nonbondTerms
+                    self._has_nbfix_terms = True
                     raise NbfixPresent('Off-diagonal Lennard-Jones elements '
                                        'found. Cannot determine LJ parameters '
                                        'for individual atoms.')
@@ -712,6 +702,10 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
         warnings.warn("1-4 scaling parameters in topology file are being ignored. "
             "This is not recommended unless you know what you are doing.")
 
+    has_1264 = 'LENNARD_JONES_CCOEF' in prmtop._raw_data.keys()
+    if has_1264:
+        parm_ccoef = [float(x) for x in prmtop._raw_data['LENNARD_JONES_CCOEF']]
+
     # Use pyopenmm implementation of OpenMM by default.
     if mm is None:
         mm = simtk.openmm
@@ -875,32 +869,52 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
                 idx = nbidx[numTypes*i+j] - 1
                 acoef[i+numTypes*j] = sqrt(parm_acoef[idx]) * afac
                 bcoef[i+numTypes*j] = parm_bcoef[idx] * bfac
-        cforce = mm.CustomNonbondedForce('(a/r6)^2-b/r6; r6=r^6;'
-                                         'a=acoef(type1, type2);'
-                                         'b=bcoef(type1, type2);')
+        if has_1264:
+            cfac = ene_conv * length_conv**4
+            ccoef = [0 for i in range(numTypes*numTypes)]
+            for i in range(numTypes):
+                for j in range(numTypes):
+                    idx = nbidx[numTypes*i+j] - 1
+                    ccoef[i+numTypes*j] = parm_ccoef[idx] * cfac
+            cforce = mm.CustomNonbondedForce('(a/r6)^2-b/r6-c/r^4; r6=r^6;'
+                                             'a=acoef(type1, type2);'
+                                             'b=bcoef(type1, type2);'
+                                             'c=ccoef(type1, type2);')
+        else:
+            cforce = mm.CustomNonbondedForce('(a/r6)^2-b/r6; r6=r^6;'
+                                             'a=acoef(type1, type2);'
+                                             'b=bcoef(type1, type2);')
         cforce.addTabulatedFunction('acoef',
                     mm.Discrete2DFunction(numTypes, numTypes, acoef))
         cforce.addTabulatedFunction('bcoef',
                     mm.Discrete2DFunction(numTypes, numTypes, bcoef))
+        if has_1264:
+            cforce.addTabulatedFunction('ccoef',
+                        mm.Discrete2DFunction(numTypes, numTypes, ccoef))
         cforce.addPerParticleParameter('type')
         for atom in prmtop._getAtomTypeIndexes():
             cforce.addParticle((atom-1,))
-        # Now set the various properties based on the NonbondedForce object
-        if nonbondedMethod in ('PME', 'Ewald', 'CutoffPeriodic'):
-            cforce.setNonbondedMethod(cforce.CutoffPeriodic)
-            cforce.setCutoffDistance(nonbondedCutoff)
-            cforce.setUseLongRangeCorrection(True)
-        elif nonbondedMethod == 'CutoffNonPeriodic':
-            cforce.setNonbondedMethod(cforce.CutoffNonPeriodic)
-            cforce.setCutoffDistance(nonbondedCutoff)
-        elif nonbondedMethod == 'NoCutoff':
-            cforce.setNonbondedMethod(cforce.NoCutoff)
-        else:
-            raise ValueError('Unrecognized cutoff option %s' % nonbondedMethod)
     else:
         for (charge, (rVdw, epsilon)) in zip(prmtop.getCharges(), nonbondTerms):
             sigma = rVdw * sigmaScale
             force.addParticle(charge, sigma, epsilon)
+        if has_1264:
+            numTypes = prmtop.getNumTypes()
+            nbidx = [int(x) for x in prmtop._raw_data['NONBONDED_PARM_INDEX']]
+            ccoef = [0 for i in range(numTypes*numTypes)]
+            ene_conv = units.kilocalories_per_mole.conversion_factor_to(units.kilojoules_per_mole)
+            length_conv = units.angstroms.conversion_factor_to(units.nanometers)
+            cfac = ene_conv * length_conv**4
+            for i in range(numTypes):
+                for j in range(numTypes):
+                    idx = nbidx[numTypes*i+j] - 1
+                    ccoef[i+numTypes*j] = parm_ccoef[idx] * cfac
+            cforce = mm.CustomNonbondedForce('-c/r^4; c=ccoef(type1, type2)')
+            cforce.addTabulatedFunction('ccoef',
+                        mm.Discrete2DFunction(numTypes, numTypes, ccoef))
+            cforce.addPerParticleParameter('type')
+            for atom in prmtop._getAtomTypeIndexes():
+                cforce.addParticle((atom-1,))
 
 
     # Add 1-4 Interactions
@@ -926,10 +940,23 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
 
     # Copy the exceptions as exclusions to the CustomNonbondedForce if we have
     # NBFIX terms
-    if nbfix:
+    if nbfix or has_1264:
         for i in range(force.getNumExceptions()):
             ii, jj, chg, sig, eps = force.getExceptionParameters(i)
             cforce.addExclusion(ii, jj)
+        # Now set the various properties based on the NonbondedForce object
+        if nonbondedMethod in ('PME', 'Ewald', 'CutoffPeriodic'):
+            cforce.setNonbondedMethod(cforce.CutoffPeriodic)
+            cforce.setCutoffDistance(nonbondedCutoff)
+            cforce.setUseLongRangeCorrection(True)
+        elif nonbondedMethod == 'CutoffNonPeriodic':
+            cforce.setNonbondedMethod(cforce.CutoffNonPeriodic)
+            cforce.setCutoffDistance(nonbondedCutoff)
+        elif nonbondedMethod == 'NoCutoff':
+            cforce.setNonbondedMethod(cforce.NoCutoff)
+        else:
+            raise ValueError('Unrecognized cutoff option %s' % nonbondedMethod)
+        # Add this force to the system
         system.addForce(cforce)
     system.addForce(force)
 
@@ -1191,9 +1218,9 @@ class AmberAsciiRestart(object):
             if hasbox:
                 boxVectors = np.zeros((3, 3), np.float32)
         else:
-            coordinates = [[0.0, 0.0, 0.0] for i in range(self.natom)]
+            coordinates = [Vec3(0.0, 0.0, 0.0) for i in range(self.natom)]
             if hasvels:
-                velocities = [[0.0, 0.0, 0.0] for i in range(self.natom)]
+                velocities = [Vec3(0.0, 0.0, 0.0) for i in range(self.natom)]
             if hasbox:
                 boxVectors = [[0.0, 0.0, 0.0] for i in range(3)]
 
@@ -1203,14 +1230,16 @@ class AmberAsciiRestart(object):
         idx = 0
         for i in range(startline, endline):
             line = lines[i]
-            coordinates[idx][0] = float(line[ 0:12])
-            coordinates[idx][1] = float(line[12:24])
-            coordinates[idx][2] = float(line[24:36])
+            x = float(line[ 0:12])
+            y = float(line[12:24])
+            z = float(line[24:36])
+            coordinates[idx] = Vec3(x, y, z)
             idx += 1
             if idx < self.natom:
-                coordinates[idx][0] = float(line[36:48])
-                coordinates[idx][1] = float(line[48:60])
-                coordinates[idx][2] = float(line[60:72])
+                x = float(line[36:48])
+                y = float(line[48:60])
+                z = float(line[60:72])
+                coordinates[idx] = Vec3(x, y, z)
                 idx += 1
         self.coordinates = units.Quantity(coordinates, units.angstroms)
         startline = endline
@@ -1220,14 +1249,16 @@ class AmberAsciiRestart(object):
             idx = 0
             for i in range(startline, endline):
                 line = lines[i]
-                velocities[idx][0] = float(line[ 0:12]) * VELSCALE
-                velocities[idx][1] = float(line[12:24]) * VELSCALE
-                velocities[idx][2] = float(line[24:36]) * VELSCALE
+                x = float(line[ 0:12]) * VELSCALE
+                y = float(line[12:24]) * VELSCALE
+                z = float(line[24:36]) * VELSCALE
+                velocities[idx] = Vec3(x, y, z)
                 idx += 1
                 if idx < self.natom:
-                    velocities[idx][0] = float(line[36:48]) * VELSCALE
-                    velocities[idx][1] = float(line[48:60]) * VELSCALE
-                    velocities[idx][2] = float(line[60:72]) * VELSCALE
+                    x = float(line[36:48]) * VELSCALE
+                    y = float(line[48:60]) * VELSCALE
+                    z = float(line[60:72]) * VELSCALE
+                    velocities[idx] = Vec3(x, y, z)
                     idx += 1
             startline = endline
             self.velocities = units.Quantity(velocities,
@@ -1242,7 +1273,7 @@ class AmberAsciiRestart(object):
             lengths = tmp[:3]
             angles = tmp[3:]
             _box_vectors_from_lengths_angles(lengths, angles, boxVectors)
-            self.boxVectors = units.Quantity(boxVectors, units.angstroms)
+            self.boxVectors = [units.Quantity(Vec3(*x), units.angstrom) for x in boxVectors]
 
 class AmberNetcdfRestart(object):
     """
@@ -1325,11 +1356,11 @@ class AmberNetcdfRestart(object):
 
         # They are already numpy -- convert to list if we don't want numpy
         if not asNumpy:
-            self.coordinates = self.coordinates.tolist()
+            self.coordinates = [Vec3(*x) for x in self.coordinates]
             if self.velocities is not None:
-                self.velocities = self.velocities.tolist()
+                self.velocities = [Vec3(*x) for x in self.velocities]
             if self.boxVectors is not None:
-                self.boxVectors = self.boxVectors.tolist()
+                self.boxVectors = [Vec3(*x) for x in self.boxVectors]
 
         # Now add the units
         self.coordinates = units.Quantity(self.coordinates, units.angstroms)
@@ -1337,7 +1368,7 @@ class AmberNetcdfRestart(object):
             self.velocities = units.Quantity(self.velocities,
                                              units.angstroms/units.picoseconds)
         if self.boxVectors is not None:
-            self.boxVectors = units.Quantity(self.boxVectors, units.angstroms)
+            self.boxVectors = [units.Quantity(x, units.angstroms) for x in self.boxVectors]
         self.time = units.Quantity(self.time, units.picosecond)
 
 def _box_vectors_from_lengths_angles(lengths, angles, boxVectors):
@@ -1427,6 +1458,7 @@ def readAmberCoordinates(filename, asNumpy=False):
         try:
             crdfile = AmberAsciiRestart(filename)
         except TypeError:
+            raise
             raise TypeError('Problem parsing %s as an ASCII Amber restart file'
                             % filename)
         except (IndexError, ValueError):

wrappers/python/tests/TestAmberPrmtopFile.py

@@ -8,7 +8,9 @@ import simtk.openmm.app.element as elem
 prmtop1 = AmberPrmtopFile('systems/alanine-dipeptide-explicit.prmtop')
 prmtop2 = AmberPrmtopFile('systems/alanine-dipeptide-implicit.prmtop')
 prmtop3 = AmberPrmtopFile('systems/ff14ipq.parm7')
+prmtop4 = AmberPrmtopFile('systems/Mg_water.prmtop')
 inpcrd3 = AmberInpcrdFile('systems/ff14ipq.rst7')
+inpcrd4 = AmberInpcrdFile('systems/Mg_water.inpcrd')
 
 class TestAmberPrmtopFile(unittest.TestCase):
 
@@ -145,35 +147,35 @@ class TestAmberPrmtopFile(unittest.TestCase):
         totalMass2 = sum([system2.getParticleMass(i) for i in range(system2.getNumParticles())]).value_in_unit(amu)
         self.assertAlmostEqual(totalMass1, totalMass2)
 
-#   def test_NBFIX_LongRange(self):
-#       """Test prmtop files with NBFIX LJ modifications w/ long-range correction"""
-#       system = prmtop3.createSystem(nonbondedMethod=PME,
-#                                     nonbondedCutoff=8*angstroms)
-#       # Check the forces
-#       has_nonbond_force = has_custom_nonbond_force = False
-#       nonbond_exceptions = custom_nonbond_exclusions = 0
-#       for force in system.getForces():
-#           if isinstance(force, NonbondedForce):
-#               has_nonbond_force = True
-#               nonbond_exceptions = force.getNumExceptions()
-#           elif isinstance(force, CustomNonbondedForce):
-#               has_custom_nonbond_force = True
-#               custom_nonbond_exceptions = force.getNumExclusions()
-#       self.assertTrue(has_nonbond_force)
-#       self.assertTrue(has_custom_nonbond_force)
-#       self.assertEqual(nonbond_exceptions, custom_nonbond_exceptions)
-#       integrator = VerletIntegrator(1.0*femtoseconds)
-#       # Use reference platform, since it should always be present and
-#       # 'working', and the system is plenty small so this won't be too slow
-#       sim = Simulation(prmtop3.topology, system, integrator, Platform.getPlatformByName('Reference'))
-#       # Check that the energy is about what we expect it to be
-#       sim.context.setPeriodicBoxVectors(*inpcrd3.boxVectors)
-#       sim.context.setPositions(inpcrd3.positions)
-#       ene = sim.context.getState(getEnergy=True, enforcePeriodicBox=True).getPotentialEnergy()
-#       ene = ene.value_in_unit(kilocalories_per_mole)
-#       # Make sure the energy is relatively close to the value we get with
-#       # Amber using this force field.
-#       self.assertAlmostEqual(-7099.44989739/ene, 1, places=3)
+    def test_NBFIX_LongRange(self):
+        """Test prmtop files with NBFIX LJ modifications w/ long-range correction"""
+        system = prmtop3.createSystem(nonbondedMethod=PME,
+                                      nonbondedCutoff=8*angstroms)
+        # Check the forces
+        has_nonbond_force = has_custom_nonbond_force = False
+        nonbond_exceptions = custom_nonbond_exclusions = 0
+        for force in system.getForces():
+            if isinstance(force, NonbondedForce):
+                has_nonbond_force = True
+                nonbond_exceptions = force.getNumExceptions()
+            elif isinstance(force, CustomNonbondedForce):
+                has_custom_nonbond_force = True
+                custom_nonbond_exceptions = force.getNumExclusions()
+        self.assertTrue(has_nonbond_force)
+        self.assertTrue(has_custom_nonbond_force)
+        self.assertEqual(nonbond_exceptions, custom_nonbond_exceptions)
+        integrator = VerletIntegrator(1.0*femtoseconds)
+        # Use reference platform, since it should always be present and
+        # 'working', and the system is plenty small so this won't be too slow
+        sim = Simulation(prmtop3.topology, system, integrator, Platform.getPlatformByName('Reference'))
+        # Check that the energy is about what we expect it to be
+        sim.context.setPeriodicBoxVectors(*inpcrd3.boxVectors)
+        sim.context.setPositions(inpcrd3.positions)
+        ene = sim.context.getState(getEnergy=True, enforcePeriodicBox=True).getPotentialEnergy()
+        ene = ene.value_in_unit(kilocalories_per_mole)
+        # Make sure the energy is relatively close to the value we get with
+        # Amber using this force field.
+        self.assertAlmostEqual(-7099.44989739/ene, 1, places=3)
 
     def test_NBFIX_noLongRange(self):
         """Test prmtop files with NBFIX LJ modifications w/out long-range correction"""
@@ -198,13 +200,46 @@ class TestAmberPrmtopFile(unittest.TestCase):
         # 'working', and the system is plenty small so this won't be too slow
         sim = Simulation(prmtop3.topology, system, integrator, Platform.getPlatformByName('Reference'))
         # Check that the energy is about what we expect it to be
-        sim.context.setPeriodicBoxVectors(*inpcrd3.boxVectors)
-        sim.context.setPositions(inpcrd3.positions)
+        sim.context.setPeriodicBoxVectors(*inpcrd3.getBoxVectors())
+        sim.context.setPositions(inpcrd3.getPositions())
         ene = sim.context.getState(getEnergy=True, enforcePeriodicBox=True).getPotentialEnergy()
         ene = ene.value_in_unit(kilocalories_per_mole)
         # Make sure the energy is relatively close to the value we get with
         # Amber using this force field.
         self.assertAlmostEqual(-7042.3903307/ene, 1, places=3)
 
+    def test_LJ1264(self):
+        """Test prmtop with 12-6-4 vdW potential implemented"""
+        system = prmtop4.createSystem(nonbondedMethod=PME,
+                                      nonbondedCutoff=8*angstroms)
+        # Check the forces
+        has_nonbond_force = has_custom_nonbond_force = False
+        nonbond_exceptions = custom_nonbond_exclusions = 0
+        for force in system.getForces():
+            if isinstance(force, NonbondedForce):
+                has_nonbond_force = True
+                nonbond_exceptions = force.getNumExceptions()
+                force.setUseDispersionCorrection(False)
+            elif isinstance(force, CustomNonbondedForce):
+                self.assertTrue(force.getUseLongRangeCorrection())
+                has_custom_nonbond_force = True
+                custom_nonbond_exceptions = force.getNumExclusions()
+                force.setUseLongRangeCorrection(False)
+        self.assertTrue(has_nonbond_force)
+        self.assertTrue(has_custom_nonbond_force)
+        self.assertEqual(nonbond_exceptions, custom_nonbond_exceptions)
+        integrator = VerletIntegrator(1.0*femtoseconds)
+        # Use reference platform, since it should always be present and
+        # 'working', and the system is plenty small so this won't be too slow
+        sim = Simulation(prmtop4.topology, system, integrator, Platform.getPlatformByName('Reference'))
+        # Check that the energy is about what we expect it to be
+        sim.context.setPeriodicBoxVectors(*inpcrd4.boxVectors)
+        sim.context.setPositions(inpcrd4.positions)
+        ene = sim.context.getState(getEnergy=True, enforcePeriodicBox=True).getPotentialEnergy()
+        ene = ene.value_in_unit(kilocalories_per_mole)
+        # Make sure the energy is relatively close to the value we get with
+        # Amber using this force field.
+        self.assertAlmostEqual(-7307.2735621/ene, 1, places=3)
+
 if __name__ == '__main__':
     unittest.main()