Openmm setup MD with intitial positions

Question

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:

1. 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)

2. 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
$$$$

Answer 1

I tried running the code that you provided: and I did observe the following issues with it:

CRYST1    8.664    8.295    8.300  90.00  90.00  90.00 P 1           1
ATOM      1  C1  PHEEX   1     -10.331   2.354   0.160  0.00  0.00      PHEE C
ATOM      2  C2  PHEEX   1      -8.818   2.260   0.255  0.00  0.00      PHEE C
ATOM      3  H3  PHEEX   1     -10.734   1.553  -0.467  0.00  0.00      PHEE H
ATOM      4  H4  PHEEX   1     -10.638   3.314  -0.268  0.00  0.00      PHEE H
ATOM      5  H5  PHEEX   1     -10.779   2.266   1.154  0.00  0.00      PHEE H
ATOM      6  C6  PHEEX   1      -8.159   2.359  -1.118  0.00  0.00      PHEE C
ATOM      7  H7  PHEEX   1      -8.549   1.314   0.738  0.00  0.00      PHEE H
ATOM      8  H8  PHEEX   1      -8.433   3.054   0.906  0.00  0.00      PHEE H
ATOM      9  H9  PHEEX   1      -8.466   1.535  -1.772  0.00  0.00      PHEE H
ATOM     10  O10 PHEEX   1      -6.735   2.409  -1.025  0.00  0.00      PHEE O
ATOM     11  H11 PHEEX   1      -8.472   3.290  -1.603  0.00  0.00      PHEE H
ATOM     12  C12 PHEEX   1      -6.140   1.152  -0.713  0.00  0.00      PHEE C
ATOM     13  C13 PHEEX   1      -4.635   1.297  -0.828  0.00  0.00      PHEE C
ATOM     14  H14 PHEEX   1      -6.491   0.384  -1.410  0.00  0.00      PHEE H
ATOM     15  H15 PHEEX   1      -6.398   0.859   0.309  0.00  0.00      PHEE H
ATOM     16  H16 PHEEX   1      -4.352   1.582  -1.847  0.00  0.00      PHEE H
ATOM     17  H17 PHEEX   1      -4.132   0.360  -0.573  0.00  0.00      PHEE H
ATOM     18  H18 PHEEX   1      -4.270   2.087  -0.163  0.00  0.00      PHEE H
END

This is the PDB file that was generated from the updated SDF file

The first issue was regarding the SDF file which was provided as a text block. I have since cleaned it up and pasted it as a code block in the question.

from openmm             import *
from openmm.app         import *
from simtk.unit         import *

from openff.toolkit     import ForceField, Molecule
from openff.units       import unit
from openff.interchange import Interchange

molecule    =   Molecule.from_file("cleaned.sdf") # use the sdf file that I have placed
Topology    =   molecule.to_topology()

coordinates =   PDBFile("phenolethylether.pdb")
positions   =   coordinates.getPositions()

Topology.box_vectors = Quantity([15, 15, 15], nanometer)

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)

This is the code that i used to check if your code was actually creating a system of 15x15x15 molecules. However, the code does not do that and instead sets the PBC conditions in x-,y- and z-axis to 15 nanometers. This is probably not what you want.

So to investigate if i can create a system consisting of 10x10x10 molecules, i used the following code to generate a PDB file containing 1000 molecules.

set pdb phenolethylether.pdb

set a 9
set b 9
set c 9

set xdist 10
set ydist 10
set zdist 10

set xend [expr $a / 2]
set xstart [expr -1 * $xend]

set yend [expr $b / 2]
set ystart [expr -1 * $yend]

set zend [expr $c / 2]
set zstart [expr -1 * $zend]

set counter 1

for { set xcounter $xstart } { $xcounter <= $xend } { incr xcounter } {
    set xtrans [expr $xdist * $xcounter]
    for { set ycounter $ystart } { $ycounter <= $yend } { incr ycounter } {
        set ytrans [expr $ydist * $ycounter]
        for { set zcounter $zstart } { $zcounter <= $zend } { incr zcounter } {
            set ztrans [expr $zdist * $zcounter]
            mol new $pdb type pdb
            set sel [atomselect top all]
            set com [measure center $sel weight mass]
            $sel moveby [vecscale -1 $com]
            $sel moveby [list $xtrans $ytrans $ztrans]
            $sel set segname M$counter
            $sel set resid $counter
            $sel writepdb test.$xcounter.$ycounter.$zcounter.pdb
            mol delete all
            incr counter
        }
    }
}

package require topotools 1.6
set molList {}

for { set xcounter $xstart } { $xcounter <= $xend } { incr xcounter } {
for {set ycounter $ystart } { $ycounter <= $yend } { incr ycounter } {
    for {set zcounter $zstart } { $zcounter <= $zend } { incr zcounter } {
        set mol [mol new test.$xcounter.$ycounter.$zcounter.pdb type pdb]
        lappend molList $mol
        }
    }
}

set mol [::TopoTools::mergemols $molList]
animate write pdb final.pdb


mol delete all
mol new final.pdb

package require pbctools 2.8

set sel [atomselect top all]
set minmaxs [measure minmax $sel]
set pbcvec [vecsub [lindex $minmaxs 1] [lindex $minmaxs 0]]

set pbcx [expr {round[lindex $pbcvec 0]]}]
set pbcy [expr {round[lindex $pbcvec 1]]}]
set pbcz [expr {round[lindex $pbcvec 2]]}]

pbc set $pbcvec
$sel writepdb final.pbc.pdb

exit

However, When i tried loading this file to provide the positions, I got an error saying openmm.OpenMMException: Called setPositions() on a Context with the wrong number of positions which is expected since the SDF file only contains 18 atoms and the PDB file now contains 18,000 atoms.

To troubleshoot this, I tried creating an SDF file from the generated PDB file. This was recognized by the code, and the positions were set. However, we now get a new error /home/ssm_grp/anaconda3/envs/openmm/lib/python3.10/site-packages/openff/interchange/smirnoff/_nonbonded.py:385: UserWarning: Warning! Partial charge method 'am1bcc' is not designed for use on large (i.e. > 150 atoms) molecules and may crash or take hours to run on this molecule (found 13122 atoms). For more, see https://docs.openforcefield.org/projects/toolkit/en/stable/faq.html#parameterizing-my-system-which-contains-a-large-molecule-is-taking-forever-whats-wrong molecule.assign_partial_charges(method) and this crashed immediately after this.

Note The reason why the code ran for you is because you were simulating only one molecule placed in a 4x4x4 $nm^3$ box, and the assignment of partial charges will proceed without any issues.

Suggestions

1. If you can, try generating files for AmberFF or CHARMM FF and used the corresponding files to setup the calculations. This would take away the issues related to partial charge calculations and topologies.
2. If not, try reducing the system size to something that is less than 150 (but this defeats the purpose of doing an MD simulation)