diff --git a/fftool b/fftool index f849c84..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/10/17 +# Agilio Padua , 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: