From ec8ccdccaada091babaf0499371d4343e8874b9e Mon Sep 17 00:00:00 2001 From: Mathieu Cancade <152990525+mcancade@users.noreply.github.com> Date: Tue, 25 Feb 2025 18:35:45 +0100 Subject: [PATCH 1/2] Write an OpenMM polarisable script with -s option --- polxml | 128 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 127 insertions(+), 1 deletion(-) diff --git a/polxml b/polxml index cc111f9..65aac3e 100755 --- a/polxml +++ b/polxml @@ -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''' @@ -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') @@ -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) @@ -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() From b20d9f57036f677ec5dc8c2ea29ec515ec9c2486 Mon Sep 17 00:00:00 2001 From: Mathieu Cancade <152990525+mcancade@users.noreply.github.com> Date: Tue, 25 Feb 2025 18:37:55 +0100 Subject: [PATCH 2/2] Write Tang-Toennis parameters directly in the omm-p.py script --- coulttxml | 101 ++++++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 80 insertions(+), 21 deletions(-) diff --git a/coulttxml b/coulttxml index e18fbfb..d26431d 100755 --- a/coulttxml +++ b/coulttxml @@ -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))" @@ -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) @@ -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) @@ -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()