I have been assigned a research project using Openmm to model Liquid-Liquid interfaces. I was able to successfully simulate bulk water, but now that Ive moved on to hexane and water I cannot get it to work. I keep getting the error: openmm.OpenMMException: Called setPositions() on a Context with the wrong number of positions, which I believe results from the disconnect with the number of molecules in my PDB and my PSF file.
Here is the Openmm code:
Code: Select all
# This script was generated by OpenMM-Setup on 2023-02-28.
import sys
from openmm import *
from openmm.app import *
from simtk.unit import *
import time
import configparser
from params_validate import *
import os
from datetime import datetime
def get_cutoff_distance(dims):
# args are box dimensions in A. Returns half the box distance in nm or 1.4 nm per Singer
cutoff_distance = 1.4
for i in dims:
if i/20 < cutoff_distance:
cutoff_distance = i/20
return cutoff_distance*nanometers
if __name__ == '__main__':
# start the clock
start = time.time()
# read in parameters. files and params are just dictionaries with the parameters that I am consistently changing
# files, params = get_files_params(sys.argv[1], 'openmm')
files = {
'str': ['MDSim/water/toppar_water_ions.str', 'MDSim/Hexane/C6H14.str'],
'rtf': ['MDSim/water/top_interface.rtf', 'MDSim/Hexane/top_all36_cgenff.rtf'],
'prm': ['MDSim/water/par_interface.prm', 'MDSim/Hexane/par_all36_cgenff.prm'],
'dcd': 'MDSim/water/results.dcd',
'results': 'MDSim/water/results.txt',
'openmm_log': 'MDSim/water/results.log',
'psf': 'MDSim/water_hexane_intf/interface_vmd.psf',
'pdb': 'MDSim/water_hexane_intf/interface_packmol_sim.pdb'
}
params = {
# equilibration parameters
'equil_steps': 10000,
# simulation parameters
'sim_steps': 50000,
'v_to_t_steps': 500, # number to steps between setting the velocities to temperature
'gpu': True,
'temp': 300
}
# get correct surface_dims
surface_dims = get_surface_dims(files['str'])
# box_dims is in nanometers
box_dims = [40, 40, 40]
# Input Files
psf = CharmmPsfFile(files['psf'])
pdb = PDBFile(files['pdb'])
ff = CharmmParameterSet(*files['rtf'], *files['prm'], *files['str'], 'MDSim/water/par_interface.prm', 'MDSim/water/top_interface.rtf')
print(f"Number of atoms in PSF file: {len(list(psf.topology.atoms()))}")# debugging
print(f"Number of positions in PDB file: {len(list(pdb.positions))}") # debugging
# System Configuration
nonbondedMethod = PME
constraints = None
rigidWater = False
# Integration Options
dt = 0.001*picoseconds
friction = 1.0/picosecond
pressure = 1.0*atmospheres
barostatInterval = 25
# intdroduce boolean to turn GPU on or off
if params['gpu']:
platform = Platform.getPlatformByName('Reference') #For CUDA, use 'CUDA'
#platformProperties = {'Precision': 'double'} #For CUDA, use 'Precision': 'double' or 'mixed' (for accuracy and speed, respectively.
else:
platform = Platform.getPlatformByName('CPU')
platformProperties = {}
# pdbReporter = PDBReporter(sys.argv[5], steps)
dcd_interval = 100
dcdReporter = DCDReporter(files['dcd'], dcd_interval)
dataReporter = StateDataReporter(files['openmm_log'], 1000, totalSteps=params['equil_steps'],
step=True, time=True, speed=True, progress=True, elapsedTime=True, remainingTime=True, potentialEnergy=True, kineticEnergy=True, totalEnergy=True, temperature=True, volume=True, density=True, separator='\t')
# Prepare the Simulation
print('Building system..')
# set the box dimensions for periodic boundaries
psf.setBox((box_dims[0])*nanometers, (box_dims[1])*nanometers, (box_dims[2])*nanometers)
# create they system object
system = psf.createSystem(ff, nonbondedMethod=nonbondedMethod,
nonbondedCutoff=get_cutoff_distance([40,40,40]),
constraints=constraints, rigidWater=rigidWater)
# use a VerletIntegrator - this is supposed to give a constant Energy Density ensemble. See OpenMM User Guide
integrator = VerletIntegrator(dt)
# create simulation object
simulation = Simulation(psf.topology, system, integrator)
# set the initial positions using PDB
simulation.context.setPositions(pdb.positions)
# Minimize and Equilibrate
print('Performing energy minimization...')
simulation.minimizeEnergy()
print('Equilibrating...')
simulation.reporters.append(dataReporter)
# use the temperature and maxwell boltzmann velocity distribution to set correct initial velocities
for i in range(int(params['equil_steps']/params['v_to_t_steps'])):
simulation.step(params['v_to_t_steps'])
simulation.context.setVelocitiesToTemperature(params['temp'])
# Simulate
print('Simulating...')
simulation.reporters.append(dcdReporter)
# simulation.reporters.append(pdbReporter)
simulation.currentStep = 0
simulation.step(params['sim_steps'])
# end simulation section -- now to write the simulation parameters to the results file
params['num_trajectory_frames'] = params['sim_steps']/dcd_interval
params['run_time'] = params['sim_steps']*dt
params['frame_spacing'] = params['run_time']/params['num_trajectory_frames']
params['equil_time'] = params['equil_steps']*dt
with open(files['results'], 'w') as f:
f.write('Openmm\n')
f.write('All temperatures are in Kelvin.\n')
for key in params:
f.write(key+' = '+str(params[key])+'\n')
dt_string = datetime.now().strftime('%d/%m/%Y %H:%M:%S')
f.write('appended on '+dt_string+'\n')
f.write('script run time = '+time.strftime('%H:%M:%S', time.gmtime(time.time()-start))+'\n\n')
The rest of our files are here: https://drive.google.com/drive/folders/ ... drive_link.
I apologize it is all a bit messy, weve been struggling and troubleshooting quite a bit. If you need us to clean up the files we can try but we dont understand enough about these programs to know what clean files look like.
Me and another student have inherited this project from a student who has now graduated. It looks like he used Openmm for bulk simulations, and solid-liquid interfaces, but never liquid-liquid interfaces. He automated a lot of his process, and I cannot follow his code very well. The instructions he left were not the most specific and were not sure where to begin. I am concerned that our current workflow may only work for the solid - liquid interfaces and liquid liquid interfaces will need a different workflow.
Our workflow is:
Use gaussian to create a PDB for a single molecule
Plug this PDB into Charmm to get Charmm forcefield files
Take this PDB and pack a 40*40*40 cell with one type of molecules using packmol
Repeat steps 1-3 for the other molecule
Use packmol to pack two 40*40*40 cells with molecules, this time stacking the two on top of each other. We now have three pdb files for packed cells.
Take the packmol outputs and our STRs from charmm and plug that into VMD to create a psf file.
VMD code:
Code: Select all
pdb MDSim/water_hexane_intf/interface_packmol_sim.pdb
package require psfgen
topology MDSim/water_hexane_intf/toppar_water_ions.str
topology MDSim/water_hexane_intf/C6H14.str
segment A {pdb MDSim/water/water_packmol_sim.pdb}
segment B {pdb MDSim/Hexane/hexane_packmol_sim.pdb}
coordpdb MDSim/water/water_packmol_sim.pdb A
coordpdb MDSim/Hexane/hexane_packmol_sim.pdb B
guesscoord
writepsf MDSim/water_hexane_intf/interface_vmd.psf
Take the PSF and the dcd resulting from openmm and view it in VMD.
For the bulk water simulation, we skipped steps 4 and 5, and were able to successfully simulate water.
When we tried a bulk hexane simulation everything worked, except the molecules did not seem packed into a cube but more of a bubble (our nonbondedmethod was PME and not noCutoff).
Our advisor asked us to move ahead with the interface and now weve encountered the issue i described in the first paragraph.
Could anybody advise, either on the error with set positions or on our workflow or if we should consider using something other than openmm?
Thanks!