Trying to complete simulation.context.setPositions(molecule.positions)
(molecule from sdf file)
The code below runs nicely (so can be used if you are beginner like me to "get started"). However, for a molecule several molecular file formats (sdf and pdb) seem to be needed which appears to duplicate data)
Quite a basic question but has taken a few days already so I am sure I am missing something obvious:
Using the code from link
import os
print (os.environ['CONDA_DEFAULT_ENV'])
# check I am using openmm environment
from openmm.app import *
from openmm import *
from openmm.unit import *
from sys import stdout
from openff.toolkit import ForceField, Molecule
from openff.units import unit
from openff.interchange import Interchange
molecule = Molecule.from_file("PEE1.sdf")
# Convert this molecule to a topology
Topology = molecule.to_topology()
pdb = PDBFile('PEE1.pdb')
positions = pdb.getPositions()
# Define periodicity via box vectors. This seems essential but I was #hoping this could be read from pdb file or sdf file
Topology.box_vectors = unit.Quantity([4, 4, 4], unit.nanometer)
# Load OpenFF 2.0.0 "Sage"
sage = ForceField("openff-2.0.0.offxml")
interchange = Interchange.from_smirnoff(topology=[molecule], force_field=sage)
openmm_system = interchange.to_openmm()
integrator = LangevinIntegrator(300 * kelvin, 1 / picosecond, 2 * femtoseconds)
simulation = Simulation(Topology, openmm_system, integrator)
simulation.context.setPositions(positions)
print(openmm_system.getNumParticles())
# Write the trajectory to a file called "output.pdb"
#simulation.reporters.append(PDBReporter('output.pdb', 1000))
simulation.saveCheckpoint("saved.txt")
# Report infomation to the screen as the simulation runs
#simulation.reporters.append(StateDataReporter(stdout, 100, step=True,potentialEnergy=True, temperature=True))
simulation.minimizeEnergy()
# Reporting parameters
#simulation.reporters.append(app.DCDReporter('trajectory.dcd', 1000))
simulation.step(10000)
with open('system.xml', 'w') as output:
output.write(XmlSerializer.serialize(openmm_system))
# Visualize the molecule
molecule.visualize(show_all_hydrogens=True)
The steps that I followed are:
- I have my molecule (for example phenol) as sdf. I can create the molecule.topology from the sdf file (and pdb file)
Issue : I can create the initial starting positions" only from the pdb version. Are positions only made available from a pdb file? (If I could write out the molecule.positions as a file I might be able to work out the format: it is some sort of binary file if I just "print" it out)
- Currently, I'm trying to save the sdf file as a pdb file (via openBabel, but when opened by openmm, openmm reports error as I have many identical atoms in the pdb file.
Issue: I can correct the pdb file manually by changing HETEROATOM 1 C
carbons to HETEROATOM C1, C2
etc, but if I had 1000 phenols (say) this would be tedious. ( I have written a script to change names of atom to atom ++ number (just add index to atom name))
I would have thought "positions" (molecule.positions) could be derived from a(ny)file type such as .xyz,.sdf,.mol etc.
Have I missed something really obvious or an easy way to generate molecule.postions? (I think that this is my last hurdle to running an md). Following tutorial does work but not when I work with eg non-aminoacid structures.
This is a heavily edited re-write of my question. I can get this example to "run" (need to check scientific validity) by manual edit of pdb file (rename ATOM to ATOM++index)
sdf file:
test.xyz
OpenBabel05202321543D
18 17 0 0 0 0 0 0 0 0999 V2000
1.1816 2.2735 -0.7735 C 0 0 0 0 0 0 0 0 0 0 0 0
0.5464 1.3284 0.2322 C 0 0 0 0 0 0 0 0 0 0 0 0
0.8836 2.0043 -1.7916 H 0 0 0 0 0 0 0 0 0 0 0 0
2.2741 2.2300 -0.7168 H 0 0 0 0 0 0 0 0 0 0 0 0
0.8696 3.3078 -0.5908 H 0 0 0 0 0 0 0 0 0 0 0 0
0.9348 1.6888 1.6606 C 0 0 0 0 0 0 0 0 0 0 0 0
-0.5442 1.3519 0.1293 H 0 0 0 0 0 0 0 0 0 0 0 0
0.8491 0.2961 0.0154 H 0 0 0 0 0 0 0 0 0 0 0 0
2.0228 1.6530 1.7769 H 0 0 0 0 0 0 0 0 0 0 0 0
0.3324 0.7544 2.5567 O 0 0 0 0 0 0 0 0 0 0 0 0
0.5843 2.7007 1.8989 H 0 0 0 0 0 0 0 0 0 0 0 0
0.6282 1.0771 3.9148 C 0 0 0 0 0 0 0 0 0 0 0 0
0.0002 0.0365 4.8207 C 0 0 0 0 0 0 0 0 0 0 0 0
0.2233 2.0663 4.1524 H 0 0 0 0 0 0 0 0 0 0 0 0
1.7137 1.0876 4.0632 H 0 0 0 0 0 0 0 0 0 0 0 0
0.3587 -0.9664 4.5689 H 0 0 0 0 0 0 0 0 0 0 0 0
0.2366 0.2433 5.8686 H 0 0 0 0 0 0 0 0 0 0 0 0
-1.0885 0.0273 4.7018 H 0 0 0 0 0 0 0 0 0 0 0 0
1 2 1 0 0 0 0
1 5 1 0 0 0 0
2 7 1 0 0 0 0
2 8 1 0 0 0 0
3 1 1 0 0 0 0
4 1 1 0 0 0 0
6 10 1 0 0 0 0
6 2 1 0 0 0 0
9 6 1 0 0 0 0
10 12 1 0 0 0 0
11 6 1 0 0 0 0
12 15 1 0 0 0 0
13 12 1 0 0 0 0
13 17 1 0 0 0 0
13 18 1 0 0 0 0
14 12 1 0 0 0 0
16 13 1 0 0 0 0
M END
$$$$
code block
. Also add a link to the tutorial which you mentioned. $\endgroup$