Reader of Tinker files (#4769)

* Add basic version of TinkerFiles

* Refactor TinkerFiles

* Update docstring, type hints, and fix bug when setting box vectors

* Small fixes

* Add unit tests for the TinkerFiles class

* Fixes and updates to TinkerFiles

* Add simuteTinker example

* Update Modeller to work with AMOEBA force fields

* Small fixes

* Relax type hinting

* Fix indices in modeller

* Fix modeller indices

* Fix type hints and usage of Quantity

* Remove numpy protector

* Add reader of .seq files

* Add topology parsing of some protein residues, waters, ions, and generic molecules.

* Miscellaneous improvements

* Update amino acids and nucleotides list

* Various fixes to XML writing, and separate XML writing into a new class

* Comments/warnings

* Add nucleic topological definitions

* Improved handling of peptide residues

* Fix for CYX (disulfide bonds)

* Refactor the topology creation methods

* General improvements, and add support for nucleic-like residues

* No need to handle MP, DP, TP

* Minor improvements

* General refactoring, add automatic determination of topology

* Add TinkerAtomType dataclass, and remove references to biotypes as they are not needed

* Re-add missing parsing of forces and scalars

* Updates to createSystem()

* Add AMOEBA forces

* Add angle-related forces to createSystem

* Add placeholders for missing forces

* Beginning of support for AmoebaMultipoleForce

* Finished support for AmoebaMultipoleForce

* Support for AmoebaVdwForce

* TinkerFiles supports vdw

* Misc updates, and add AmoebaTorsionTorsion, AmoebaWcaDispersion, and AmoebaGeneralizedKirkwood

* Remove XML writer

* Fixes

* Fix wrong indentation in _findBitorsions

* Remove pdb debugging

* Documentation and fixes

* Remove files

* Revert checks in AmoebaVdwForceBuilder and ## @private  markers

* Remove duplicated static methods _getChiralAtomIndex

* Fix GK force

* Fix WcaDispersion force

* Fix WcaDisp

* Fixes and updates

* Cleanup and removing duplicated code

* Bug fixes

* A few more unit conversions

* Minor cleanup

* Misc fixes and updates

* Fix Add AmoebaStretchBendForce

* Simplify force builders

* Update ForceField

* Fix AmoebaPiTorsionForce

* Only add AmoebaWcaDispersionForce if using implicitSolvent

* Simplify amoebaforces

* Stretch torsion and angle torsion

* Misc. fixes

* Improve tests

* Fix cap group identification

* Add/improve tests

* Remove whitespaces from residue names

* Improve tests

* Consistent use of atomClasses list

* Fix match condition in AmoebaOutOfPlaneBendForceBuilder

* Fix AmoebaStretchBendForce

* Final fix for AmoebaStretchBendForce

* Fix AmoebaAngleForce

* Small fixes and improvements

* Update assertion tolerances

* Simplify torsion-torsion force creation

* Small fixes in the tests

* Review comments, type hints, docs for tinkerfiles.py

* Only use standard PDB for AA

* Type hint and docs for amoebaforces

* Reduce tolerances for failing tests

* Fixed error with ZOnly axis type when x particle is not specified

---------
Co-authored-by: peastman <peastman@stanford.edu>

docs-source/usersguide/application/02_running_sims.rst

@@ -369,6 +369,40 @@ on the :class:`CharmmPsfFile`.
 
 Note that both the CHARMM and XPLOR versions of the :file:`psf` file format are supported.
 
+.. _using-tinker-files:
+
+Using Tinker Files
+******************
+
+OpenMM can also load files created by Tinker.  At present, only the AMOEBA force
+field is supported.  Tinker uses an :file:`xyz` file to store the coordinates and
+topology, and one or more files ending in :file:`key` or :file:`prm` to store
+the force field and simulation parameters.  To load them use a :class:`TinkerFiles`
+object.  This is illustrated in the following script.  It can be found in
+OpenMM’s :file:`examples` folder with the name :file:`simulateTinker.py`.
+
+.. samepage::
+    ::
+
+        from openmm.app import *
+        from openmm import *
+        from openmm.unit import *
+        from sys import stdout
+
+        tinker = TinkerFiles('amoeba_solvated_phenol.xyz', ['amoeba_phenol.prm', 'amoebabio18.prm'])
+        system = tinker.createSystem(nonbondedMethod=PME, nonbondedCutoff=0.7*nanometer, vdwCutoff=0.9*nanometer)
+        integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.001*picoseconds)
+        simulation = Simulation(tinker.topology, system, integrator)
+        simulation.context.setPositions(tinker.positions)
+        simulation.minimizeEnergy()
+        simulation.reporters.append(DCDReporter('output.dcd', 1000))
+        simulation.reporters.append(StateDataReporter(stdout, 1000, step=True, potentialEnergy=True, temperature=True))
+        simulation.step(10000)
+
+    .. caption::
+
+        :autonumber:`Example,Tinker example`
+
 .. _the-script-builder-application:
 
 The OpenMM-Setup Application

examples/python-examples/simulateTinker.py

+from openmm.app import *
+from openmm import *
+from openmm.unit import *
+from sys import stdout
+
+tinker = TinkerFiles('amoeba_solvated_phenol.xyz', ['amoeba_phenol.prm', 'amoebabio18.prm'])
+system = tinker.createSystem(nonbondedMethod=PME, nonbondedCutoff=0.7*nanometer, vdwCutoff=0.9*nanometer)
+integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.001*picoseconds)
+simulation = Simulation(tinker.topology, system, integrator)
+simulation.context.setPositions(tinker.positions)
+simulation.minimizeEnergy()
+simulation.reporters.append(DCDReporter('output.dcd', 1000))
+simulation.reporters.append(StateDataReporter(stdout, 1000, step=True, potentialEnergy=True, temperature=True))
+simulation.step(10000)

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

@@ -1480,7 +1480,7 @@ double AmoebaReferenceMultipoleForce::calculateElectrostaticPairIxn(const Multip
 
 void AmoebaReferenceMultipoleForce::mapTorqueToForceForParticle(const MultipoleParticleData& particleI,
                                                                 const MultipoleParticleData& particleU,
-                                                                const MultipoleParticleData& particleV,
+                                                                      MultipoleParticleData* particleV,
                                                                       MultipoleParticleData* particleW,
                                                                       int axisType, const Vec3& torque,
                                                                       vector<Vec3>& forces) const
@@ -1521,7 +1521,15 @@ void AmoebaReferenceMultipoleForce::mapTorqueToForceForParticle(const MultipoleP
     Vec3 vectorU = particleU.position - particleI.position;
     norms[U] = normalizeVec3(vectorU);
 
-    Vec3 vectorV = particleV.position - particleI.position;
+    Vec3 vectorV;
+    if (axisType != AmoebaMultipoleForce::ZOnly)
+        vectorV = particleV->position - particleI.position;
+    else {
+        if (fabs(vectorU[0]/norms[U]) < 0.866)
+            vectorV = Vec3(1, 0, 0);
+        else
+            vectorV = Vec3(0, 1, 0);
+    }
     norms[V] = normalizeVec3(vectorV);
 
     Vec3 vectorW;
@@ -1585,7 +1593,7 @@ void AmoebaReferenceMultipoleForce::mapTorqueToForceForParticle(const MultipoleP
             forces[particleU.particleIndex][ii]                 -=  forceU;
 
             double forceV                                        =  vectorUV[ii]*factor3 + factor4*vectorVW[ii];
-            forces[particleV.particleIndex][ii]                 -=  forceV;
+            forces[particleV->particleIndex][ii]                 -=  forceV;
 
             forces[particleI.particleIndex][ii]                 +=  (forceU + forceV);
         }
@@ -1645,7 +1653,7 @@ void AmoebaReferenceMultipoleForce::mapTorqueToForceForParticle(const MultipoleP
         forces[particleU.particleIndex]        -= forceU;
 
         Vec3 forceV               = (vectorS*angles[VS][1] - t1*angles[VS][0])*factor3;
-        forces[particleV.particleIndex]        -= forceV;
+        forces[particleV->particleIndex]        -= forceV;
 
         Vec3 forceW               = (vectorS*angles[WS][1] - t2*angles[WS][0])*factor4;
         forces[particleW->particleIndex]       -= forceW;
@@ -1679,7 +1687,7 @@ void AmoebaReferenceMultipoleForce::mapTorqueToForceForParticle(const MultipoleP
             dw /= 3.0;
 
             forces[particleU.particleIndex][ii] -= du;
-            forces[particleV.particleIndex][ii] -= dv;
+            forces[particleV->particleIndex][ii] -= dv;
             if (particleW)
                 forces[particleW->particleIndex][ii] -= dw;
             forces[particleI.particleIndex][ii] += (du + dv + dw);
@@ -1712,7 +1720,8 @@ void AmoebaReferenceMultipoleForce::mapTorqueToForce(vector<MultipoleParticleDat
     for (unsigned int ii = 0; ii < particleData.size(); ii++) {
         if (axisTypes[ii] != AmoebaMultipoleForce::NoAxisType) {
              mapTorqueToForceForParticle(particleData[ii],
-                                         particleData[multipoleAtomZs[ii]], particleData[multipoleAtomXs[ii]],
+                                         particleData[multipoleAtomZs[ii]],
+                                         multipoleAtomXs[ii] > -1 ? &particleData[multipoleAtomXs[ii]] : NULL,
                                          multipoleAtomYs[ii] > -1 ? &particleData[multipoleAtomYs[ii]] : NULL,
                                          axisTypes[ii], torques[ii], forces);
         }

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

-/* Portions copyright (c) 2006-2022 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2025 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -1104,7 +1104,7 @@ protected:
      */
     void mapTorqueToForceForParticle(const MultipoleParticleData& particleI,
                                      const MultipoleParticleData& particleU,
-                                     const MultipoleParticleData& particleV,
+                                           MultipoleParticleData* particleV,
                                            MultipoleParticleData* particleW,
                                      int axisType, const Vec3& torque, std::vector<OpenMM::Vec3>& forces) const;
 

wrappers/python/CMakeLists.txt

@@ -78,6 +78,7 @@ foreach(SUBDIR ${SUBDIRS})
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*psf"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/charmm22.*"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/par*.inp"
+        "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.xyz"
     )
     foreach(file ${STAGING_INPUT_FILES1})
         set(STAGING_INPUT_FILES ${STAGING_INPUT_FILES} "${file}")

wrappers/python/openmm/app/__init__.py

@@ -22,6 +22,7 @@ from .pdbreporter import PDBReporter, PDBxReporter
 from .xtcreporter import XTCReporter
 from .amberprmtopfile import AmberPrmtopFile, HCT, OBC1, OBC2, GBn, GBn2
 from .amberinpcrdfile import AmberInpcrdFile
+from .tinkerfiles import TinkerFiles
 from .dcdfile import DCDFile
 from .gromacsgrofile import GromacsGroFile
 from .gromacstopfile import GromacsTopFile

wrappers/python/openmm/app/modeller.py

@@ -38,7 +38,7 @@ from openmm.app import Topology, PDBFile, ForceField
 from openmm.app.forcefield import AllBonds, CutoffNonPeriodic, CutoffPeriodic, DrudeGenerator, _getDataDirectories
 from openmm.app.internal import compiled
 from openmm.vec3 import Vec3
-from openmm import System, Context, NonbondedForce, CustomNonbondedForce, HarmonicBondForce, HarmonicAngleForce, VerletIntegrator, LangevinIntegrator, LocalEnergyMinimizer
+from openmm import System, Context, NonbondedForce, AmoebaVdwForce, AmoebaMultipoleForce, CustomNonbondedForce, HarmonicBondForce, HarmonicAngleForce, VerletIntegrator, LangevinIntegrator, LocalEnergyMinimizer
 from openmm.unit import nanometer, molar, elementary_charge, degree, acos, is_quantity, dot, norm, kilojoules_per_mole
 import openmm.unit as unit
 from . import element as elem
@@ -310,16 +310,20 @@ class Modeller(object):
         # Determine the total charge of the system
         system = forcefield.createSystem(self.topology, residueTemplates=residueTemplates)
         for i in range(system.getNumForces()):
-            if isinstance(system.getForce(i), NonbondedForce):
+            if isinstance(system.getForce(i), (NonbondedForce, AmoebaMultipoleForce)):
                 nonbonded = system.getForce(i)
                 break
         else:
             raise ValueError('The ForceField does not specify a NonbondedForce')
 
         totalCharge = 0.0
-        for i in range(nonbonded.getNumParticles()):
-            nb_i = nonbonded.getParticleParameters(i)
-            totalCharge += nb_i[0].value_in_unit(elementary_charge)
+        if isinstance(nonbonded, AmoebaMultipoleForce):
+            for i in range(nonbonded.getNumMultipoles()):
+                totalCharge += nonbonded.getMultipoleParameters(i)[0].value_in_unit(elementary_charge)
+        else:
+            for i in range(nonbonded.getNumParticles()):
+                totalCharge += nonbonded.getParticleParameters(i)[0].value_in_unit(elementary_charge)
+
         # Round to nearest integer
         totalCharge = int(floor(0.5 + totalCharge))
 
@@ -519,15 +523,18 @@ class Modeller(object):
         system = forcefield.createSystem(self.topology, residueTemplates=residueTemplates)
         nonbonded = None
         for i in range(system.getNumForces()):
-            if isinstance(system.getForce(i), NonbondedForce):
+            if isinstance(system.getForce(i), (NonbondedForce, AmoebaVdwForce)):
                 nonbonded = system.getForce(i)
+      
         if nonbonded is None:
             raise ValueError('The ForceField does not specify a NonbondedForce')
         cutoff = [waterRadius]*system.getNumParticles()
+
         for i in range(system.getNumParticles()):
             params = nonbonded.getParticleParameters(i)
             if params[2] != 0*kilojoules_per_mole:
                 cutoff[i] += params[1].value_in_unit(nanometer)*vdwRadiusPerSigma
+
         waterCutoff = waterRadius
         if len(cutoff) == 0:
             maxCutoff = waterCutoff

wrappers/python/tests/systems/alanine-dipeptide-implicit.xyz

+    22
+     1  C      2.000000    2.090000    0.000000   221     2     4     5     6
+     2  C      3.427000    2.641000   -0.000000   223     1     3     7
+     3  O      4.391000    1.877000   -0.000000   224     2
+     4  H      2.000000    1.000000   -0.000000   222     1
+     5  H      1.486000    2.454000    0.890000   222     1
+     6  H      1.486000    2.454000   -0.890000   222     1
+     7  N      3.555000    3.970000   -0.000000     7     2     8    11
+     8  CA     4.853000    4.614000   -0.000000     8     7     9    12    13
+     9  C      4.713000    6.129000    0.000000     9     8    10    17
+    10  O      3.601000    6.653000    0.000000    11     9
+    11  HN     2.733000    4.556000   -0.000000    10     7
+    12  H      5.408000    4.316000    0.890000    12     8
+    13  C      5.661000    4.221000   -1.232000    13     8    14    15    16
+    14  H      5.123000    4.521000   -2.131000    14    13
+    15  H      6.630000    4.719000   -1.206000    14    13
+    16  H      5.809000    3.141000   -1.241000    14    13
+    17  N      5.846000    6.835000    0.000000   227     9    18    19
+    18  C      5.846000    8.284000    0.000000   229    17    20    21    22
+    19  HN     6.737000    6.359000   -0.000000   228    17
+    20  H      4.819000    8.648000    0.000000   230    18
+    21  H      6.360000    8.648000    0.890000   230    18
+    22  H      6.360000    8.648000   -0.890000   230    18