From 03d5b99bfd6aa4d243c1a8972093f6e0105a3bf3 Mon Sep 17 00:00:00 2001 From: Mathieu Cancade <152990525+mcancade@users.noreply.github.com> Date: Tue, 25 Feb 2025 18:30:40 +0100 Subject: [PATCH 1/2] Write a standard OpenMM script when -x option is selected --- fftool | 237 +++++++++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 189 insertions(+), 48 deletions(-) diff --git a/fftool b/fftool index f849c84..16ce8c3 100755 --- a/fftool +++ b/fftool @@ -1,6 +1,6 @@ #!/usr/bin/env python # fftool - generate force field parameters for molecular system -# Agilio Padua , version 2024/10/17 +# Agilio Padua , version 2024/02/20 # Copyright 2013 Agilio Padua # @@ -82,7 +82,7 @@ def atomic_number(name): def hexiflarge(i, ndig): """ - Convert number to hex starting with A000... if more than n digits. + Convert number to hex starting with A if more than n digits. VMD, OpenMM convention in PDB files with natom > 99999, nmol > 9999 """ @@ -95,22 +95,6 @@ def hexiflarge(i, ndig): i = i - declim + hexlim return f'{i:X}' -def intfromhexiflarge(a): - """ - Convert string containing hex number (large systems) to int. - VMD, OpenMM convention in PDB files with natom > 99999, nmol > 9999 - """ - - if a.strip().isdigit(): - return int(a) - - ndig = len(a) - declim = 10**ndig - hexlim = int('A' + (ndig-1)*'0', 16) - - i = int(a, 16) - return i - hexlim + declim - # -------------------------------------- class vec3(object): @@ -1969,6 +1953,172 @@ 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\n\n') + + f.write('import sys\n') + f.write('import datetime\n') + f.write('import numpy as np\n\n') + + f.write('import openmm\n') + f.write('from openmm import app\n') + f.write('from openmm import unit\n\n') + + f.write('from mdtraj.reporters import HDF5Reporter\n\n') + + f.write("field = 'field.xml'\n") + f.write("config = 'config.pdb'\n") + f.write("#statefile = 'state-eq.xml'\n\n") + + f.write('temperature = 300.0*unit.kelvin\n') + f.write('pressure = 1*unit.bar\n\n') + + f.write("print('#', datetime.datetime.now())\n") + f.write("print()\n\n") + + f.write("### Input data ###\n") + f.write("print('#', field, config)\n") + f.write("forcefield = app.ForceField(field)\n") + f.write("pdb = app.PDBFile(config)\n") + f.write("# If PDBx/mmCIF format is used as config file:\n") + f.write("#pdb = app.PDBxFile(config)\n\n") + + f.write("### Modeller (add extra particles if necessary) ###\n") + f.write("modeller = app.Modeller(pdb.topology, pdb.positions)\n") + f.write("#print('# adding extra particles')\n") + f.write("#modeller.addExtraParticles(forcefield)\n\n") + + f.write("### Print topology (molecules, atoms and bonds) ###\n") + f.write("print('# ', modeller.topology.getNumResidues(), 'molecules',\n") + f.write(" modeller.topology.getNumAtoms(), 'atoms',\n") + f.write(" modeller.topology.getNumBonds(), 'bonds')\n\n") + + f.write("### Print box dimensions ###\n") + f.write("lx = modeller.topology.getUnitCellDimensions().x\n") + f.write("ly = modeller.topology.getUnitCellDimensions().y\n") + f.write("lz = modeller.topology.getUnitCellDimensions().z\n") + f.write("print('# box', lx, ly, lz, 'nm')\n\n") + + f.write("### System creation ###\n") + f.write("system = forcefield.createSystem(modeller.topology, nonbondedMethod=app.PME,\n") + f.write(" nonbondedCutoff=12.0*unit.angstrom, constraints=app.HBonds,\n") + f.write(" ewaldErrorTolerance=1.0e-5)\n\n") + + f.write("### Integrator selection ###\n") + f.write("print('# Nose-Hoover integrator')\n") + f.write("integrator = openmm.NoseHooverIntegrator(temperature, 5/unit.picosecond, 1*unit.femtosecond)\n") + + f.write("#print('# Langevin integrator', temperature)\n") + f.write("#integrator = openmm.LangevinIntegrator(temperature, 5/unit.picosecond, 1*unit.femtosecond)\n\n") + + f.write("### Barostat selection ###\n") + f.write("print('# barostat', pressure)\n") + f.write("barostat = openmm.MonteCarloBarostat(pressure, temperature)\n") + f.write("system.addForce(barostat)\n\n") + + f.write("### Plateform selection ###\n") + f.write("platform = openmm.Platform.getPlatformByName('CUDA')\n") + f.write("#platform = openmm.Platform.getPlatformByName('OpenCL')\n") + f.write("#platform.setPropertyDefaultValue('Precision', 'mixed')\n") + f.write("#properties = {'DeviceIndex': '1', 'Precision': 'mixed'}\n") + f.write("#properties = {'Precision': 'mixed'}\n") + f.write("properties = {'Precision': 'single'}\n\n") + + f.write("### Forces ###\n") + f.write("# Force settings before creating Simulation\n") + f.write("for i, f in enumerate(system.getForces()):\n") + f.write(" f.setForceGroup(i)\n\n") + + f.write("### Simulation creation ###\n") + f.write("sim = app.Simulation(modeller.topology, system, integrator, platform, properties)\n") + f.write("sim.context.setPositions(modeller.positions)\n") + f.write("sim.context.setVelocitiesToTemperature(temperature)\n\n") + + f.write("# Uncomment if loading a statefile:\n") + f.write("#print('# coordinates and velocities from', statefile)\n") + f.write("#sim.loadState(statefile)\n\n") + + f.write("# Uncomment if loading a restart file:\n") + f.write("#print('# coordinates and velocities from restart.chk')\n") + f.write("#sim.loadCheckpoint('restart.chk')\n\n") + + f.write("state = sim.context.getState()\n") + f.write("sim.topology.setPeriodicBoxVectors(state.getPeriodicBoxVectors())\n\n") + + f.write("platform = sim.context.getPlatform()\n") + f.write("print('# platform', platform.getName())\n") + f.write("for prop in platform.getPropertyNames():\n") + f.write(" print('# ', prop, platform.getPropertyValue(sim.context, prop))\n\n") + + f.write("state = sim.context.getState(getEnergy=True)\n") + f.write("print('# PotentielEnergy', state.getPotentialEnergy())\n\n") + + f.write("for i, f in enumerate(system.getForces()):\n") + f.write(" state = sim.context.getState(getEnergy=True, groups={i})\n") + f.write(" print('# ', f.getName(), state.getPotentialEnergy())\n\n") + + f.write("### Reporters ###\n") + f.write("# Introduce h5 reporter format (if necessary)\n") + f.write("h5reporter = (HDF5Reporter('traj.h5', 1000, enforcePeriodicBox=False))\n\n") + + f.write("sim.reporters = []\n") + f.write("sim.reporters.append(app.StateDataReporter(sys.stdout, 1000, step=True,\n") + f.write(" speed=True, temperature=True, separator='\t',\n") + f.write(" totalEnergy=True, potentialEnergy=True, density=True))\n\n") + + f.write("# Remove comments to select format (default: PDB)\n") + f.write("sim.reporters.append(app.PDBReporter('traj.pdb', 1000))\n") + f.write("#sim.reporters.append(h5reporter)\n") + f.write("#sim.reporters.append(app.DCDReporter('traj.dcd', 1000))\n") + f.write("#sim.reporters.append(app.CheckpointReporter('restart.chk', 10000))\n\n") + + f.write("### Run simulation ###\n") + f.write("for i in range(10):\n") + f.write(" sim.step(1000)\n\n") + + f.write("### Close reporters ###\n") + f.write("# Close h5 reporter (if used)\n") + f.write("h5reporter.close()\n\n") + + f.write("### End-of-simulation data ###\n") + f.write("print('# Force groups')\n") + f.write("for i, f in enumerate(system.getForces()):\n") + f.write(" state = sim.context.getState(getEnergy=True, groups={i})\n") + f.write(" print('# ', f.getName(), state.getPotentialEnergy())\n\n") + + f.write("state = sim.context.getState(getPositions=True)\n") + f.write("coords = state.getPositions()\n") + f.write("sim.topology.setPeriodicBoxVectors(state.getPeriodicBoxVectors())\n\n") + + f.write("### Store last configuration (PDB or PDBx format) ###\n") + f.write("app.PDBFile.writeFile(sim.topology, coords, open('last.pdb', 'w'))\n") + f.write("#app.PDBxFile.writeFile(sim.topology, coords, open('last.mmcif', 'w'))\n\n") + + f.write("sim.context.setTime(0)\n") + f.write("sim.context.setStepCount(0)\n") + f.write("sim.saveState('state-run.xml')\n") + f.write("print('# state saved to state-run.xml')\n\n") + + f.write("print()\n") + f.write("print('#', datetime.datetime.now())\n") def writepdb(self): with open('config.pdb', 'w') as f: @@ -2095,9 +2245,8 @@ class system(object): self.checkdupresidues() - nmol = natom = nbond = nangle = ndihed = 0 + natom = nbond = nangle = ndihed = 0 for sp in self.spec: - nmol += sp.nmol natom += sp.nmol * len(sp.atom) nbond += sp.nmol * len(sp.bond) nangle += sp.nmol * len(sp.angle) @@ -2186,10 +2335,9 @@ class system(object): c4 = -4 * dit.par[3] + eps c5 = 0.0 dimpr = ET.SubElement(torsionforce, 'Improper') - # in OpenMM atom 1 (not 3) is the central atom of an improper dihedral - dimpr.set('class1', dit.katp) - dimpr.set('class2', dit.iatp) - dimpr.set('class3', dit.jatp) + dimpr.set('class1', dit.iatp) + dimpr.set('class2', dit.jatp) + dimpr.set('class3', dit.katp) dimpr.set('class4', dit.latp) dimpr.set('c0', f'{c0:.5f}') dimpr.set('c1', f'{c1:.5f}') @@ -2246,9 +2394,8 @@ class system(object): ET.indent(root, space=' ') ftree.write('field.xml') - - if nmol <= intfromhexiflarge('FFFF') and natom <= intfromhexiflarge('FFFFF'): - self.writepdb() + + self.writepdb() self.writepdbxmmcif() @@ -2267,9 +2414,8 @@ class system(object): self.checkdupresidues() - nmol = natom = nbond = nangle = ndihed = 0 + natom = nbond = nangle = ndihed = 0 for sp in self.spec: - nmol += sp.nmol natom += sp.nmol * len(sp.atom) nbond += sp.nmol * len(sp.bond) nangle += sp.nmol * len(sp.angle) @@ -2379,10 +2525,9 @@ class system(object): c4 = -4 * di.par[3] + eps c5 = 0.0 dimpr = ET.SubElement(torsionforce, 'Improper') - # in OpenMM atom 1 (not 3) is the central atom of an improper dihedral - dimpr.set('class1', tp[2]) - dimpr.set('class2', tp[0]) - dimpr.set('class3', tp[1]) + dimpr.set('class1', tp[0]) + dimpr.set('class2', tp[1]) + dimpr.set('class3', tp[2]) dimpr.set('class4', tp[3]) dimpr.set('c0', f'{c0:.5f}') dimpr.set('c1', f'{c1:.5f}') @@ -2450,9 +2595,8 @@ class system(object): ET.indent(root, space=' ') ftree.write('field.xml') - - if nmol <= intfromhexiflarge('FFFF') and natom <= intfromhexiflarge('FFFFF'): - self.writepdb() + + self.writepdb() self.writepdbxmmcif() @@ -2994,7 +3138,7 @@ def main(): help = 'connect bonds across periodic boundaries in '\ 'x, xy, xyz, etc. (default: none)') parser.add_argument('-x', '--xml', action = 'store_true', - help = 'create OpenMM files .xml .pdb or .mmcif '\ + help = 'create OpenMM files .xml .pdb '\ '(needs simbox.xyz created by Packmol)') parser.add_argument('--types', action='store_true', help='unique atom types for non-bonded pairs in xml') @@ -3023,6 +3167,10 @@ def main(): files = args.infiles[1::2] # odd elements are molecule files nmol = sum(int(n) for n in nmols) + if args.charmm or args.gmx: + if nmol > 34575: + raise ValueError('the maximum number of molecules in a PDB file is 34575') + if args.box and args.rho != 0.0: raise RuntimeError('supply density or box dimensions, not both') @@ -3057,7 +3205,7 @@ def main(): 'natom nbond source charge') connect = True spec = [] - i = natom = 0 + i = 0 for zfile in files: spec.append(mol(zfile, connect, box)) spec[i].nmol = int(nmols[i]) @@ -3066,7 +3214,6 @@ def main(): .format(zfile, spec[i].name, spec[i].nmol, spec[i].ffile, len(spec[i].atom), len(spec[i].bond), spec[i].topol, spec[i].charge())) - natom += spec[i].nmol * len(spec[i].atom) i += 1 sim = system(spec, box, args.mix) @@ -3085,22 +3232,17 @@ 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') - if nmol > intfromhexiflarge('FFFF') or natom > intfromhexiflarge('FFFFF'): - print('maximum number of molecules or atoms in a PDB file exceeded, pdb file not written') + sim.writeopenmm() + print('xml files\n field.xml\n config.mmcif\n omm.py') if args.types: sim.writexmltypes(args.mix, args.allpairs) else: sim.writexml(args.mix, args.allpairs) if args.charmm: - if nmol > intfromhexiflarge('FFFF') or natom > intfromhexiflarge('FFFFF'): - raise ValueError('maximum number of molecules or atoms in a PDB file exceeded') sim.readcoords('simbox.xyz') print('charmm files\n field.str\n data.psf\n config.pdb') sim.writecharmm(args.mix, args.allpairs) if args.gmx: - if nmol > intfromhexiflarge('FFFF') or natom > intfromhexiflarge('FFFFF'): - raise ValueError('maximum number of molecules or atoms in a PDB file exceeded') sim.readcoords('simbox.xyz') print('gromacs files\n run.mdp\n field.top\n config.pdb') sim.writegmx(args.mix) @@ -3118,5 +3260,4 @@ def main(): if __name__ == '__main__': - main() - + main() \ No newline at end of file From 8b23fc5aff44c9dcd0ff8cc9b82b9452039e8c88 Mon Sep 17 00:00:00 2001 From: Mathieu Cancade <152990525+mcancade@users.noreply.github.com> Date: Thu, 6 Mar 2025 10:31:45 +0100 Subject: [PATCH 2/2] Up-to-date fftool version --- fftool | 348 +++++++++++++++++++++++++++++++-------------------------- 1 file changed, 188 insertions(+), 160 deletions(-) diff --git a/fftool b/fftool index 16ce8c3..50e44b3 100755 --- a/fftool +++ b/fftool @@ -1,6 +1,6 @@ #!/usr/bin/env python # fftool - generate force field parameters for molecular system -# Agilio Padua , version 2024/02/20 +# Agilio Padua , version 2025/03/06 # Copyright 2013 Agilio Padua # @@ -82,7 +82,7 @@ def atomic_number(name): def hexiflarge(i, ndig): """ - Convert number to hex starting with A if more than n digits. + Convert number to hex starting with A000... if more than n digits. VMD, OpenMM convention in PDB files with natom > 99999, nmol > 9999 """ @@ -95,6 +95,22 @@ def hexiflarge(i, ndig): i = i - declim + hexlim return f'{i:X}' +def intfromhexiflarge(a): + """ + Convert string containing hex number (large systems) to int. + VMD, OpenMM convention in PDB files with natom > 99999, nmol > 9999 + """ + + if a.strip().isdigit(): + return int(a) + + ndig = len(a) + declim = 10**ndig + hexlim = int('A' + (ndig-1)*'0', 16) + + i = int(a, 16) + return i - hexlim + declim + # -------------------------------------- class vec3(object): @@ -1972,153 +1988,155 @@ class system(object): f.write(' {:d} dihedrals /'.format(ndihed)) f.write('/\n\n') - f.write('# #!/usr/bin/env python\n\n') + f.write( +"""# #!/usr/bin/env python - f.write('import sys\n') - f.write('import datetime\n') - f.write('import numpy as np\n\n') +import sys +import datetime +import numpy as np - f.write('import openmm\n') - f.write('from openmm import app\n') - f.write('from openmm import unit\n\n') +import openmm +from openmm import app +from openmm import unit - f.write('from mdtraj.reporters import HDF5Reporter\n\n') +from mdtraj.reporters import HDF5Reporter - f.write("field = 'field.xml'\n") - f.write("config = 'config.pdb'\n") - f.write("#statefile = 'state-eq.xml'\n\n") +field = 'field.xml' +config = 'config.pdb' +#statefile = 'state-eq.xml' - f.write('temperature = 300.0*unit.kelvin\n') - f.write('pressure = 1*unit.bar\n\n') +temperature = 300.0*unit.kelvin +pressure = 1*unit.bar - f.write("print('#', datetime.datetime.now())\n") - f.write("print()\n\n") +print('#', datetime.datetime.now()) +print() - f.write("### Input data ###\n") - f.write("print('#', field, config)\n") - f.write("forcefield = app.ForceField(field)\n") - f.write("pdb = app.PDBFile(config)\n") - f.write("# If PDBx/mmCIF format is used as config file:\n") - f.write("#pdb = app.PDBxFile(config)\n\n") +### 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) - f.write("### Modeller (add extra particles if necessary) ###\n") - f.write("modeller = app.Modeller(pdb.topology, pdb.positions)\n") - f.write("#print('# adding extra particles')\n") - f.write("#modeller.addExtraParticles(forcefield)\n\n") - - f.write("### Print topology (molecules, atoms and bonds) ###\n") - f.write("print('# ', modeller.topology.getNumResidues(), 'molecules',\n") - f.write(" modeller.topology.getNumAtoms(), 'atoms',\n") - f.write(" modeller.topology.getNumBonds(), 'bonds')\n\n") - - f.write("### Print box dimensions ###\n") - f.write("lx = modeller.topology.getUnitCellDimensions().x\n") - f.write("ly = modeller.topology.getUnitCellDimensions().y\n") - f.write("lz = modeller.topology.getUnitCellDimensions().z\n") - f.write("print('# box', lx, ly, lz, 'nm')\n\n") - - f.write("### System creation ###\n") - f.write("system = forcefield.createSystem(modeller.topology, nonbondedMethod=app.PME,\n") - f.write(" nonbondedCutoff=12.0*unit.angstrom, constraints=app.HBonds,\n") - f.write(" ewaldErrorTolerance=1.0e-5)\n\n") - - f.write("### Integrator selection ###\n") - f.write("print('# Nose-Hoover integrator')\n") - f.write("integrator = openmm.NoseHooverIntegrator(temperature, 5/unit.picosecond, 1*unit.femtosecond)\n") +### 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) - f.write("#print('# Langevin integrator', temperature)\n") - f.write("#integrator = openmm.LangevinIntegrator(temperature, 5/unit.picosecond, 1*unit.femtosecond)\n\n") - - f.write("### Barostat selection ###\n") - f.write("print('# barostat', pressure)\n") - f.write("barostat = openmm.MonteCarloBarostat(pressure, temperature)\n") - f.write("system.addForce(barostat)\n\n") - - f.write("### Plateform selection ###\n") - f.write("platform = openmm.Platform.getPlatformByName('CUDA')\n") - f.write("#platform = openmm.Platform.getPlatformByName('OpenCL')\n") - f.write("#platform.setPropertyDefaultValue('Precision', 'mixed')\n") - f.write("#properties = {'DeviceIndex': '1', 'Precision': 'mixed'}\n") - f.write("#properties = {'Precision': 'mixed'}\n") - f.write("properties = {'Precision': 'single'}\n\n") - - f.write("### Forces ###\n") - f.write("# Force settings before creating Simulation\n") - f.write("for i, f in enumerate(system.getForces()):\n") - f.write(" f.setForceGroup(i)\n\n") - - f.write("### Simulation creation ###\n") - f.write("sim = app.Simulation(modeller.topology, system, integrator, platform, properties)\n") - f.write("sim.context.setPositions(modeller.positions)\n") - f.write("sim.context.setVelocitiesToTemperature(temperature)\n\n") - - f.write("# Uncomment if loading a statefile:\n") - f.write("#print('# coordinates and velocities from', statefile)\n") - f.write("#sim.loadState(statefile)\n\n") - - f.write("# Uncomment if loading a restart file:\n") - f.write("#print('# coordinates and velocities from restart.chk')\n") - f.write("#sim.loadCheckpoint('restart.chk')\n\n") - - f.write("state = sim.context.getState()\n") - f.write("sim.topology.setPeriodicBoxVectors(state.getPeriodicBoxVectors())\n\n") - - f.write("platform = sim.context.getPlatform()\n") - f.write("print('# platform', platform.getName())\n") - f.write("for prop in platform.getPropertyNames():\n") - f.write(" print('# ', prop, platform.getPropertyValue(sim.context, prop))\n\n") - - f.write("state = sim.context.getState(getEnergy=True)\n") - f.write("print('# PotentielEnergy', state.getPotentialEnergy())\n\n") - - f.write("for i, f in enumerate(system.getForces()):\n") - f.write(" state = sim.context.getState(getEnergy=True, groups={i})\n") - f.write(" print('# ', f.getName(), state.getPotentialEnergy())\n\n") - - f.write("### Reporters ###\n") - f.write("# Introduce h5 reporter format (if necessary)\n") - f.write("h5reporter = (HDF5Reporter('traj.h5', 1000, enforcePeriodicBox=False))\n\n") - - f.write("sim.reporters = []\n") - f.write("sim.reporters.append(app.StateDataReporter(sys.stdout, 1000, step=True,\n") - f.write(" speed=True, temperature=True, separator='\t',\n") - f.write(" totalEnergy=True, potentialEnergy=True, density=True))\n\n") - - f.write("# Remove comments to select format (default: PDB)\n") - f.write("sim.reporters.append(app.PDBReporter('traj.pdb', 1000))\n") - f.write("#sim.reporters.append(h5reporter)\n") - f.write("#sim.reporters.append(app.DCDReporter('traj.dcd', 1000))\n") - f.write("#sim.reporters.append(app.CheckpointReporter('restart.chk', 10000))\n\n") - - f.write("### Run simulation ###\n") - f.write("for i in range(10):\n") - f.write(" sim.step(1000)\n\n") - - f.write("### Close reporters ###\n") - f.write("# Close h5 reporter (if used)\n") - f.write("h5reporter.close()\n\n") - - f.write("### End-of-simulation data ###\n") - f.write("print('# Force groups')\n") - f.write("for i, f in enumerate(system.getForces()):\n") - f.write(" state = sim.context.getState(getEnergy=True, groups={i})\n") - f.write(" print('# ', f.getName(), state.getPotentialEnergy())\n\n") - - f.write("state = sim.context.getState(getPositions=True)\n") - f.write("coords = state.getPositions()\n") - f.write("sim.topology.setPeriodicBoxVectors(state.getPeriodicBoxVectors())\n\n") - - f.write("### Store last configuration (PDB or PDBx format) ###\n") - f.write("app.PDBFile.writeFile(sim.topology, coords, open('last.pdb', 'w'))\n") - f.write("#app.PDBxFile.writeFile(sim.topology, coords, open('last.mmcif', 'w'))\n\n") - - f.write("sim.context.setTime(0)\n") - f.write("sim.context.setStepCount(0)\n") - f.write("sim.saveState('state-run.xml')\n") - f.write("print('# state saved to state-run.xml')\n\n") - - f.write("print()\n") - f.write("print('#', datetime.datetime.now())\n") +#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: @@ -2245,8 +2263,9 @@ class system(object): self.checkdupresidues() - natom = nbond = nangle = ndihed = 0 + nmol = natom = nbond = nangle = ndihed = 0 for sp in self.spec: + nmol += sp.nmol natom += sp.nmol * len(sp.atom) nbond += sp.nmol * len(sp.bond) nangle += sp.nmol * len(sp.angle) @@ -2335,9 +2354,10 @@ class system(object): c4 = -4 * dit.par[3] + eps c5 = 0.0 dimpr = ET.SubElement(torsionforce, 'Improper') - dimpr.set('class1', dit.iatp) - dimpr.set('class2', dit.jatp) - dimpr.set('class3', dit.katp) + # in OpenMM atom 1 (not 3) is the central atom of an improper dihedral + dimpr.set('class1', dit.katp) + dimpr.set('class2', dit.iatp) + dimpr.set('class3', dit.jatp) dimpr.set('class4', dit.latp) dimpr.set('c0', f'{c0:.5f}') dimpr.set('c1', f'{c1:.5f}') @@ -2394,8 +2414,9 @@ class system(object): ET.indent(root, space=' ') ftree.write('field.xml') - - self.writepdb() + + if nmol <= intfromhexiflarge('FFFF') and natom <= intfromhexiflarge('FFFFF'): + self.writepdb() self.writepdbxmmcif() @@ -2414,8 +2435,9 @@ class system(object): self.checkdupresidues() - natom = nbond = nangle = ndihed = 0 + nmol = natom = nbond = nangle = ndihed = 0 for sp in self.spec: + nmol += sp.nmol natom += sp.nmol * len(sp.atom) nbond += sp.nmol * len(sp.bond) nangle += sp.nmol * len(sp.angle) @@ -2525,9 +2547,10 @@ class system(object): c4 = -4 * di.par[3] + eps c5 = 0.0 dimpr = ET.SubElement(torsionforce, 'Improper') - dimpr.set('class1', tp[0]) - dimpr.set('class2', tp[1]) - dimpr.set('class3', tp[2]) + # in OpenMM atom 1 (not 3) is the central atom of an improper dihedral + dimpr.set('class1', tp[2]) + dimpr.set('class2', tp[0]) + dimpr.set('class3', tp[1]) dimpr.set('class4', tp[3]) dimpr.set('c0', f'{c0:.5f}') dimpr.set('c1', f'{c1:.5f}') @@ -2595,8 +2618,9 @@ class system(object): ET.indent(root, space=' ') ftree.write('field.xml') - - self.writepdb() + + if nmol <= intfromhexiflarge('FFFF') and natom <= intfromhexiflarge('FFFFF'): + self.writepdb() self.writepdbxmmcif() @@ -3138,7 +3162,7 @@ def main(): help = 'connect bonds across periodic boundaries in '\ 'x, xy, xyz, etc. (default: none)') parser.add_argument('-x', '--xml', action = 'store_true', - help = 'create OpenMM files .xml .pdb '\ + help = 'create OpenMM files .xml .pdb or .mmcif '\ '(needs simbox.xyz created by Packmol)') parser.add_argument('--types', action='store_true', help='unique atom types for non-bonded pairs in xml') @@ -3167,10 +3191,6 @@ def main(): files = args.infiles[1::2] # odd elements are molecule files nmol = sum(int(n) for n in nmols) - if args.charmm or args.gmx: - if nmol > 34575: - raise ValueError('the maximum number of molecules in a PDB file is 34575') - if args.box and args.rho != 0.0: raise RuntimeError('supply density or box dimensions, not both') @@ -3205,7 +3225,7 @@ def main(): 'natom nbond source charge') connect = True spec = [] - i = 0 + i = natom = 0 for zfile in files: spec.append(mol(zfile, connect, box)) spec[i].nmol = int(nmols[i]) @@ -3214,6 +3234,7 @@ def main(): .format(zfile, spec[i].name, spec[i].nmol, spec[i].ffile, len(spec[i].atom), len(spec[i].bond), spec[i].topol, spec[i].charge())) + natom += spec[i].nmol * len(spec[i].atom) i += 1 sim = system(spec, box, args.mix) @@ -3233,16 +3254,22 @@ def main(): if args.xml: sim.readcoords('simbox.xyz') sim.writeopenmm() - print('xml files\n field.xml\n config.mmcif\n omm.py') + 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: sim.writexmltypes(args.mix, args.allpairs) else: sim.writexml(args.mix, args.allpairs) if args.charmm: + if nmol > intfromhexiflarge('FFFF') or natom > intfromhexiflarge('FFFFF'): + raise ValueError('maximum number of molecules or atoms in a PDB file exceeded') sim.readcoords('simbox.xyz') print('charmm files\n field.str\n data.psf\n config.pdb') sim.writecharmm(args.mix, args.allpairs) if args.gmx: + if nmol > intfromhexiflarge('FFFF') or natom > intfromhexiflarge('FFFFF'): + raise ValueError('maximum number of molecules or atoms in a PDB file exceeded') sim.readcoords('simbox.xyz') print('gromacs files\n run.mdp\n field.top\n config.pdb') sim.writegmx(args.mix) @@ -3260,4 +3287,5 @@ def main(): if __name__ == '__main__': - main() \ No newline at end of file + main() +