Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Write a standard polarisable OpenMM script (with polxml) + add directly the Tang-Toennis parameters inside (with coulttxml) #42

Open
wants to merge 2 commits into
base: main
Choose a base branch
from
+207 −22
Open

Write a standard polarisable OpenMM script (with polxml) + add directly the Tang-Toennis parameters inside (with coulttxml) #42

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ec8ccdc
Write an OpenMM polarisable script with -s option
mcancade Feb 25, 2025
b20d9f5
Write Tang-Toennis parameters directly in the omm-p.py script
mcancade Feb 25, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
101 changes: 80 additions & 21 deletions coulttxml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -88,8 +88,6 @@ def count_molecules_from_pdbx(pdbx_filename):
# -------------------------------------

before = '''
# OpenMM code to add charge-dipole damping in polarizable simulations

# CoulTT damping function
condition = "(don1*acc2 + don2*acc1) *"
coulttexp = condition + "q1*q2/(fpe0*r)*(- c*exp(-b*r) * (1 + b*r + (b*r)^2/2 + (b*r)^3/6 + (b*r)^4/24))"
Expand Down Expand Up @@ -129,6 +127,8 @@ def main():
help = 'PDB file with configuration')
parser.add_argument('--core', action = 'store_true',
help = 'Use core charge instead of induced dipole charge in TT damping')
parser.add_argument('--s', metavar = 'omm-p.py',
help = 'OpenMM input script')
args = parser.parse_args()

print("force field from", args.xml)
Expand All @@ -138,6 +138,11 @@ def main():
residues_xml = root.find('Residues').findall('Residue')
residue_names_xml = [residue.attrib['name'] for residue in residues_xml]

atomtypes_xml = root.find('AtomTypes').findall('Type')
atomtypes_names_xml = [type.attrib['name'] for type in atomtypes_xml]
atomtypes_masses_xml = [type.attrib['mass'] for type in atomtypes_xml]

light_atoms_mass_cutoff = 2.00

if 'pdb' in args.pdb:
num_molecules = count_molecules_from_pdb(args.pdb)
Expand All @@ -149,25 +154,79 @@ def main():
print("Error: The number of different residues in the XML does not match the count from the PDB or PDBx file.")
sys.exit(1)

with open('addCoulTT.py', 'w') as f:
f.write(before)
# Loop over each residue and its corresponding count from the PDB
f.write("# Put flag 1 in 1st or 2nd field in atoms to be screened\n")
f.write("# donor (1st field) is the charge, acceptor (2nd field) the induced dipole\n")
for residue_xml, count in zip(residue_names_xml, num_molecules):
f.write(f"for i in range({count}):\n")
atoms = next((res for res in residues_xml if res.attrib['name'] == residue_xml), None).findall('Atom')
for i, atom in enumerate(atoms):
name = atom.attrib['name']
charge = float(atom.attrib['charge'])
if not args.core and i < len(atoms) - 1:
nxtname = atoms[i+1].attrib['name']
if nxtname.startswith('D'):
charge = - float(atoms[i+1].attrib['charge'])
f.write(f" CoulTT.addParticle([0, 0, {charge:8.5f}]) # {name:4s} {atom.attrib['type']:12s}\n")
f.write(after)

print("OpenMM code for CoulTT written to addCoulTT.py")
if args.s is not None:
with open("omm-p-sc.py", 'w') as outputfile:
with open(args.s, 'r') as file:
data = file.readlines()
line = 0
while line < len(data):
count = 0
if data[line] == "### Forces ###\n":
del data[line:line+1]
outputfile.write("### Forces ###\n")
outputfile.write("# Force settings before creating Simulation\n")

outputfile.write(before)
outputfile.write("# Tang-Toennis screening\n")
outputfile.write("# Put flag 1 in 1st or 2nd field in atoms to be screened\n")
outputfile.write("# donor (1st field) is the charge, acceptor (2nd field) the induced dipole\n")
for residue_xml, count in zip(residue_names_xml, num_molecules):
outputfile.write(f"for i in range({count}):\n")
atoms = next((res for res in residues_xml if res.attrib['name'] == residue_xml), None).findall('Atom')
for i, atom in enumerate(atoms):
name = atom.attrib['name']
charge = float(atom.attrib['charge'])
if not args.core and i < len(atoms) - 1:
nxtname = atoms[i+1].attrib['name']
if nxtname.startswith('D'):
charge = - float(atoms[i+1].attrib['charge'])
for k in range(0, len(atomtypes_xml)):
if atom.attrib['type'] == atomtypes_names_xml[k]:
if atomtypes_xml[k].attrib["class"] == "DRUD":
outputfile.write(f" CoulTT.addParticle([0, 1, {charge:8.5f}]) # {name:4s} {atom.attrib['type']:12s}\n")
elif float(atomtypes_xml[k].attrib["mass"]) < light_atoms_mass_cutoff:
outputfile.write(f" CoulTT.addParticle([0, 0, {charge:8.5f}]) # {name:4s} {atom.attrib['type']:12s}\n")
elif float(atomtypes_xml[k].attrib["mass"]) > light_atoms_mass_cutoff:
outputfile.write(f" CoulTT.addParticle([0, 1, {charge:8.5f}]) # {name:4s} {atom.attrib['type']:12s}\n")

outputfile.write(after)
count = count+1

elif count == 0:
outputfile.write(data[line])

line = line +1

print("Tang-Toennis parameters written to omm-p-sc.py")

elif args.s is None:
with open('addCoulTT.py', 'w') as f:
f.write("# OpenMM code to add charge-dipole damping in polarizable simulations\n")
f.write(before)
# Loop over each residue and its corresponding count from the PDB
f.write("# Put flag 1 in 1st or 2nd field in atoms to be screened\n")
f.write("# donor (1st field) is the charge, acceptor (2nd field) the induced dipole\n")
for residue_xml, count in zip(residue_names_xml, num_molecules):
f.write(f"for i in range({count}):\n")
atoms = next((res for res in residues_xml if res.attrib['name'] == residue_xml), None).findall('Atom')
for i, atom in enumerate(atoms):
name = atom.attrib['name']
charge = float(atom.attrib['charge'])
if not args.core and i < len(atoms) - 1:
nxtname = atoms[i+1].attrib['name']
if nxtname.startswith('D'):
charge = - float(atoms[i+1].attrib['charge'])
for k in range(0, len(atomtypes_xml)):
if atom.attrib['type'] == atomtypes_names_xml[k]:
if atomtypes_xml[k].attrib["class"] == "DRUD":
f.write(f" CoulTT.addParticle([0, 1, {charge:8.5f}]) # {name:4s} {atom.attrib['type']:12s}\n")
elif float(atomtypes_xml[k].attrib["mass"]) < light_atoms_mass_cutoff:
f.write(f" CoulTT.addParticle([0, 0, {charge:8.5f}]) # {name:4s} {atom.attrib['type']:12s}\n")
elif float(atomtypes_xml[k].attrib["mass"]) > light_atoms_mass_cutoff:
f.write(f" CoulTT.addParticle([0, 1, {charge:8.5f}]) # {name:4s} {atom.attrib['type']:12s}\n")
f.write(after)

print("OpenMM code for CoulTT written to addCoulTT.py")

if __name__ == '__main__':
main()
128 changes: 127 additions & 1 deletion polxml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -366,6 +366,126 @@ class Topology(object):
f'{label_atom_id_j:5s} {label_atom_id_i:5s} 0 0 covale\n')
f.write('data_\n')

def rewrite_openmm(self, omm_script, config_input):
with open("omm-p.py", 'w') as outputfile:
with open(omm_script, 'r') as file:
data = file.readlines()
line = 0
while line < len(data):
count = 0
if data[line] == "# created by fftool\n":
outputfile.write("# created by fftool and polxml\n")
outputfile.write('# {:d} atoms /'.format(self.natom))
outputfile.write(' {:d} Drude particles //\n\n'.format(self.ndrude))
del data[line:line+3]

if data[line] == "field = 'field.xml'\n":
outputfile.write("field = 'field-p.xml'\n")
count = count+1

if data[line] == "config = 'config.pdb'\n":
if '.mmcif' in config_input:
outputfile.write("config = 'config-p.mmcif'\n")
count = count+1
else :
outputfile.write("config = 'config-p.pdb'\n")
count = count+1

if data[line] == "config = 'config.mmcif'\n":
outputfile.write("config = 'config-p.mmcif'\n")
count = count+1

if data[line] == "pdb = app.PDBFile(config)\n":
if '.mmcif' in config_input:
outputfile.write("#pdb = app.PDBFile(config)\n")
outputfile.write("# If PDBx/mmCIF format is used as config file:\n")
outputfile.write("pdb = app.PDBxFile(config)\n\n")
del data[line:line+3]
count = count+1

if data[line] == "### Integrator selection ###\n":
del data[line:line+5]
outputfile.write("### Integrator selection ###\n")
outputfile.write("#print('# Nose-Hoover integrator')\n")
outputfile.write("#integrator = openmm.NoseHooverIntegrator(temperature, 5/unit.picosecond, 1*unit.femtosecond)\n")
outputfile.write("print('Drude Nose-Hoover integrator', temperature)\n")
outputfile.write("integrator = openmm.DrudeNoseHooverIntegrator(temperature, 5/unit.picosecond,\n")
outputfile.write(" 1*unit.kelvin, 1/unit.picosecond, 1*unit.femtosecond)\n")
outputfile.write("#print('# Drude Langevin integrator')\n")
outputfile.write("#integrator = openmm.DrudeLangevinIntegrator(temperature, 5/unit.picosecond,\n")
outputfile.write("# 1*unit.kelvin, 20/unit.picosecond, 1*unit.femtosecond)\n")
outputfile.write("integrator.setMaxDrudeDistance(0.2*unit.angstrom)\n")
outputfile.write("print('# max Drude distance', integrator.getMaxDrudeDistance())\n\n")
count = count+1

if data[line] == "### Run simulation ###\n":
del data[line:line+4]
outputfile.write("### Constants ###\n")
outputfile.write("kB = unit.BOLTZMANN_CONSTANT_kB/(unit.joule/unit.kelvin)\n")
outputfile.write("NA = unit.AVOGADRO_CONSTANT_NA*unit.mole\n\n")

outputfile.write("### Drude section ###\n")
outputfile.write("iat = [ i for i, atom in enumerate(modeller.topology.atoms()) if atom.name[0] != 'D' ]\n")
outputfile.write("idr = [ i for i, atom in enumerate(modeller.topology.atoms()) if atom.name[0] == 'D' ]\n")
outputfile.write("nat = len(iat)\n")
outputfile.write("ndr = len(idr)\n\n")

outputfile.write("nall = modeller.topology.getNumAtoms()\n")
outputfile.write("mall = np.array([ system.getParticleMass(i)/unit.dalton for i in range(nall) ])\n\n")

outputfile.write("ncons = system.getNumConstraints()\n\n")

outputfile.write("# Reduced mass of DC-DP pairs\n")
outputfile.write("mu = np.zeros(nall)\n")
outputfile.write("for i in idr:\n")
outputfile.write(" mu[i] = 1.0/(1.0/mall[i-1] + 1.0/mall[i])\n")
outputfile.write("mu = mu.reshape((nall, 1))\n")
outputfile.write("vdr = np.zeros((nall, 3))\n\n")

outputfile.write("# Add Drude masses back to cores\n")
outputfile.write("mat = np.copy(mall)\n")
outputfile.write("for i in idr:\n")
outputfile.write(" mat[i-1] += 0.4\n")
outputfile.write("mat = mat.take(iat)\n")
outputfile.write("mat = mat.reshape((nat, 1))\n\n")

outputfile.write("mall = mall.reshape((nall, 1))\n\n")

outputfile.write("print('#', nat, 'atoms', ndr, 'DP', ncons, 'constraints')\n")
outputfile.write("print('# running...')\n\n")

outputfile.write("dof_all = 3*nall - ncons\n")
outputfile.write("dof_at = 3*nat - ncons\n")
outputfile.write("dof_dr = 3*ndr\n\n")

outputfile.write("### Run simulation ###\n")
outputfile.write("for i in range(10):\n")
outputfile.write(" sim.step(1000)\n")
outputfile.write(" state = sim.context.getState(getVelocities=True)\n")
outputfile.write(" vel = state.getVelocities(asNumpy=True)/(unit.nanometer/unit.picosecond)\n")
outputfile.write(" Tall = np.sum(mall*vel**2)/(dof_all*kB)*(1e3/NA)*unit.kelvin\n")
outputfile.write(" vat = vel.take(iat, axis=0)\n")
outputfile.write(" Tat = np.sum(mat*vat**2)/(dof_at*kB)*(1e3/NA)*unit.kelvin\n\n")

outputfile.write(" for i in idr:\n")
outputfile.write(" vdr[i] = vel[i] - vel[i-1]\n")
outputfile.write(" Tdr = np.sum(mu*vdr**2)/(dof_dr*kB)*(1e3/NA)*unit.kelvin\n")
outputfile.write(" print('# Tall', Tall, 'Tatoms', Tat, 'Tdrude', Tdr)\n\n")
count = count+1

if data[line] == "### Store last configuration (PDB or PDBx format) ###\n":
if '.mmcif' in config_input:
outputfile.write("### Store last configuration (PDB or PDBx format) ###\n")
outputfile.write("#app.PDBFile.writeFile(sim.topology, coords, open('last.pdb', 'w'))\n")
outputfile.write("app.PDBxFile.writeFile(sim.topology, coords, open('last.mmcif', 'w'))\n\n")
del data[line:line+3]
count = count+1

elif count == 0:
outputfile.write(data[line])

line = line +1

def settypes(self, ff):
'''assign types from ff to atoms in topology'''

Expand Down Expand Up @@ -644,7 +764,7 @@ class Forcefield(object):
break

ljforce = self.root.find('LennardJonesForce')
if ljforce:
if ljforce is not None:
dp = ET.SubElement(ljforce, 'Atom')
dp.set('class', 'DRUD')
dp.set('sigma', '1.0')
Expand All @@ -668,6 +788,8 @@ def main():
help = 'PDB or PDBx/mmcif file with configuration (default: config.pdb)')
parser.add_argument('-op', '--outpdb', default = 'config-p.pdb',
help = 'PDB or PDBx/mmcif file with configuration (default: config-p.pdb)')
parser.add_argument('-s', '--inscript',
help = 'input OpenMM script')
args = parser.parse_args()

config = Topology(args.inpdb)
Expand All @@ -690,6 +812,10 @@ def main():
config.writepdb(args.outpdb)

print('configuration written to', args.outpdb)

if args.inscript is not None:
config.rewrite_openmm(args.inscript, args.inpdb)
print('OpenMM python script modified to omm-p.py')

if __name__ == '__main__':
main()