Write a standard OpenMM script when -x option is selected

fftool

@@ -1,6 +1,6 @@
 #!/usr/bin/env python
 # fftool - generate force field parameters for molecular system
-# Agilio Padua <[email protected]>, version 2024/10/17
+# Agilio Padua <[email protected]>, version 2025/03/06
 
 # Copyright 2013 Agilio Padua
 #
@@ -1969,6 +1969,174 @@ class system(object):
 
             # f.write('\n')
 
+    def writeopenmm(self):
+        natom = nbond = nangle = ndihed = 0
+        for sp in self.spec:
+            natom += sp.nmol * len(sp.atom)
+            nbond += sp.nmol * len(sp.bond)
+            nangle += sp.nmol * len(sp.angle)
+            ndihed += sp.nmol * (len(sp.dihed) + len(sp.dimpr))
+
+        with open('omm.py', 'w') as f:
+            f.write('# created by fftool\n')
+            f.write('# {:d} atoms /'.format(natom))
+            if nbond is not None:
+                f.write(' {:d} bonds /'.format(nbond))
+            if nangle is not None:
+                f.write(' {:d} angles /'.format(nangle))
+            if ndihed is not None:
+                f.write(' {:d} dihedrals /'.format(ndihed))
+            f.write('/\n\n')
+
+            f.write(
+"""# #!/usr/bin/env python
+
+import sys
+import datetime
+import numpy as np
+
+import openmm
+from openmm import app
+from openmm import unit
+
+from mdtraj.reporters import HDF5Reporter
+
+field = 'field.xml'
+config = 'config.pdb'
+#statefile = 'state-eq.xml'
+
+temperature = 300.0*unit.kelvin
+pressure = 1*unit.bar
+
+print('#', datetime.datetime.now())
+print()
+
+### Input data ###
+print('#', field, config)
+forcefield = app.ForceField(field)
+pdb = app.PDBFile(config)
+# If PDBx/mmCIF format is used as config file:
+#pdb = app.PDBxFile(config)
+
+### Modeller (add extra particles if necessary) ###
+modeller = app.Modeller(pdb.topology, pdb.positions)
+#print('#   adding extra particles')
+#modeller.addExtraParticles(forcefield)
+
+### Print topology (molecules, atoms and bonds) ###
+print('#  ', modeller.topology.getNumResidues(), 'molecules',
+    modeller.topology.getNumAtoms(), 'atoms',
+    modeller.topology.getNumBonds(), 'bonds')
+
+### Print box dimensions ###
+lx = modeller.topology.getUnitCellDimensions().x
+ly = modeller.topology.getUnitCellDimensions().y
+lz = modeller.topology.getUnitCellDimensions().z
+print('#   box', lx, ly, lz, 'nm')
+
+### System creation ###
+system = forcefield.createSystem(modeller.topology, nonbondedMethod=app.PME,
+    nonbondedCutoff=12.0*unit.angstrom, constraints=app.HBonds,
+    ewaldErrorTolerance=1.0e-5)
+
+### Integrator selection ###
+print('# Nose-Hoover integrator')
+integrator = openmm.NoseHooverIntegrator(temperature, 5/unit.picosecond, 1*unit.femtosecond)
+
+#print('# Langevin integrator', temperature)
+#integrator = openmm.LangevinIntegrator(temperature, 5/unit.picosecond, 1*unit.femtosecond)
+
+### Barostat selection ###
+print('#   barostat', pressure)
+barostat = openmm.MonteCarloBarostat(pressure, temperature)
+system.addForce(barostat)
+
+### Plateform selection ###
+platform = openmm.Platform.getPlatformByName('CUDA')
+#platform = openmm.Platform.getPlatformByName('OpenCL')
+#platform.setPropertyDefaultValue('Precision', 'mixed')
+#properties = {'DeviceIndex': '1', 'Precision': 'mixed'}
+#properties = {'Precision': 'mixed'}
+properties = {'Precision': 'single'}
+
+### Forces ###
+# Force settings before creating Simulation
+for i, f in enumerate(system.getForces()):
+    f.setForceGroup(i)
+
+### Simulation creation ###
+sim = app.Simulation(modeller.topology, system, integrator, platform, properties)
+sim.context.setPositions(modeller.positions)
+sim.context.setVelocitiesToTemperature(temperature)
+
+# Uncomment if loading a statefile:
+#print('# coordinates and velocities from', statefile)
+#sim.loadState(statefile)
+
+# Uncomment if loading a restart file:
+#print('# coordinates and velocities from restart.chk')
+#sim.loadCheckpoint('restart.chk')
+
+state = sim.context.getState()
+sim.topology.setPeriodicBoxVectors(state.getPeriodicBoxVectors())
+
+platform = sim.context.getPlatform()
+print('# platform', platform.getName())
+for prop in platform.getPropertyNames():
+    print('#  ', prop, platform.getPropertyValue(sim.context, prop))
+
+state = sim.context.getState(getEnergy=True)
+print('# PotentielEnergy', state.getPotentialEnergy())
+
+for i, f in enumerate(system.getForces()):
+    state = sim.context.getState(getEnergy=True, groups={i})
+    print('#  ', f.getName(), state.getPotentialEnergy())
+
+### Reporters ###
+# Introduce h5 reporter format (if necessary)
+h5reporter = (HDF5Reporter('traj.h5', 1000, enforcePeriodicBox=False))
+
+sim.reporters = []
+sim.reporters.append(app.StateDataReporter(sys.stdout, 1000, step=True,
+    speed=True, temperature=True, separator='\t',
+    totalEnergy=True, potentialEnergy=True, density=True))
+
+# Remove comments to select format (default: PDB)
+sim.reporters.append(app.PDBReporter('traj.pdb', 1000))
+#sim.reporters.append(h5reporter)
+#sim.reporters.append(app.DCDReporter('traj.dcd', 1000))
+#sim.reporters.append(app.CheckpointReporter('restart.chk', 10000))
+
+### Run simulation ###
+for i in range(10):
+    sim.step(1000)
+
+### Close reporters ###
+# Close h5 reporter (if used)
+h5reporter.close()
+
+### End-of-simulation data ###
+print('# Force groups')
+for i, f in enumerate(system.getForces()):
+    state = sim.context.getState(getEnergy=True, groups={i})
+    print('#  ', f.getName(), state.getPotentialEnergy())
+
+state = sim.context.getState(getPositions=True)
+coords = state.getPositions()
+sim.topology.setPeriodicBoxVectors(state.getPeriodicBoxVectors())
+
+### Store last configuration (PDB or PDBx format) ###
+app.PDBFile.writeFile(sim.topology, coords, open('last.pdb', 'w'))
+#app.PDBxFile.writeFile(sim.topology, coords, open('last.mmcif', 'w'))
+
+sim.context.setTime(0)
+sim.context.setStepCount(0)
+sim.saveState('state-run.xml')
+print('# state saved to state-run.xml')
+
+print()
+print('#', datetime.datetime.now())"""
+            )
 
     def writepdb(self):
         with open('config.pdb', 'w') as f:
@@ -3085,7 +3253,8 @@ def main():
         sim.writelmp(args.mix, args.allpairs, args.units)
     if args.xml:
         sim.readcoords('simbox.xyz')
-        print('xml files\n  field.xml\n  config.pdb\n  config.mmcif')
+        sim.writeopenmm()
+        print('xml files\n  field.xml\n  config.pdb\n  config.mmcif\n  omm.py')
         if nmol > intfromhexiflarge('FFFF') or natom > intfromhexiflarge('FFFFF'):
             print('maximum number of molecules or atoms in a PDB file exceeded, pdb file not written')
         if args.types: