6.3 Error - Spearate indices parsed into a single index

[George Pantelopulos]

Hi Peter,

Thanks for the explanation! I understand how to use the topology object a bit better well now, too.

That said, there still seems to be some sort of issue going on, despite the fact that the topology object being created has the correct number of chains...

I had originally thought that there was an error due to some particle positions very quickly becoming NAN in simulation. After that I had found the minimized structure to contain these strange "chain-like" structures of water... which led me to think that the system was being treated as one or many big chains that should not have been connected, which lead to my original post.

I've just gone through and done two long minimizations (50,000 steps) with Openmm and with GROMACS on a slightly different CHARMM-GUI-built system and the resulting structures are strikingly different. There is no presence of any "chain-like" water structures in the GROMACS-minimized system, despite these being present in the Openmm-minimized system. I produced PDBs of both of these minimized systems using the PDBWriter, and even so VMD was able to distinguish water molecules as being unique structures in the GROMACS-minimized system, while recognizing water molecules to be chains in the Openmm-Minimized system... It's really bizarre and I can't help but think something else might be up.

I've attached the initial gromacs files and these two pdb files to this message if you'd like to take a look. Thanks again for your help!

sharing.tar.gz

(2.8 MiB) Downloaded 25 times

[George Pantelopulos]

Additionally, it might be useful to try to run a few steps of equilibration to see what's going on before it blows up...

Code: Select all

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

totalsteps = 500
savesteps  = 5
temp = 298.15

gro = GromacsGroFile('gmx_minimization.gro')
top = GromacsTopFile('topol.top', periodicBoxVectors=gro.getPeriodicBoxVectors())

print('Constructing system...')
system = top.createSystem(nonbondedMethod=CutoffNonPeriodic, nonbondedCutoff=1.0*nanometer, constraints=HBonds, removeCMMotion=True)

system.addForce(MonteCarloMembraneBarostat(1.0*bar, 0.0*bar*nanometer, temp*kelvin, MonteCarloMembraneBarostat.XYIsotropic, MonteCarloMembraneBarostat.ZFixed, 25))
integrator = LangevinIntegrator(temp*kelvin, 1.0/picosecond, 2.0*femtoseconds)
platform = Platform.getPlatformByName('CUDA')
properties = {'CudaPrecision': 'mixed'}
simulation = Simulation(top.topology, system, integrator, platform, properties)

simulation.context.setPositions(gro.positions)
simulation.context.setVelocitiesToTemperature(temp*kelvin)

print('Setting up reporters...')
simulation.reporters.append(app.StateDataReporter(stdout, savesteps, step=True, 
    potentialEnergy=True, kineticEnergy=True, totalEnergy=True, temperature=True,
    progress=True, remainingTime=True, speed=True, totalSteps=totalsteps, 
    separator='\t'))
simulation.reporters.append(PDBReporter('equilibration.pdb', savesteps))

print('Running equilibration for %s steps...'%totalsteps)
simulation.step(totalsteps)

If you look at the few frames that come out of this, the problem might be more apparent.

[George Pantelopulos]

Regarding my above post,

If you run a few steps of simulation at 0.2 fs timestep and a very low temperature like ~50K or so the system won't explode... and you should be able to see that the waters are being drawn together into the lipid bilayer. If you save frames frequently enough (like every 50 steps over 5000 steps) it can be seen that waters are being drawn into chain-like conformations. This absolutely should not happen nor should it be expected.

[Peter Eastman]

I see what the problem is. Your TIP3.itp file doesn't include a flexible version of the molecule, and OpenMM relies on that for working out the topology. For example, the following is from the gmx.ff/tip3p.itp file bundled with Gromacs:

Code: Select all

#ifndef FLEXIBLE
[ settles ]
; i	j	funct	length
1	1	0.09572	0.15139

[ exclusions ]
1	2	3
2	1	3
3	1	2
#else
[ bonds ]
; i	j	funct	length	force.c.
1	2	1	0.09572	502416.0 0.09572	502416.0 
1	3	1	0.09572	502416.0 0.09572	502416.0 
	

[ angles ]
; i	j	k	funct	angle	force.c.
2	1	3	1	104.52	628.02	104.52	628.02	
#endif

When OpenMM parses the files, it defines FLEXIBLE and expects to therefore get flexible water molecules. Then when you call createSystem(), if rigidWater is set to True, it replaces any bonds and angles in water molecules with constraints. But in your case, there aren't any. So each water molecule ends up as three separate atoms with no connection between them.

Could you swap in one of the standard tip3p files and see if that works correctly?

By the way, I should probably have asked this before. If you're using CHARMM-GUI to build your system, and if you're then running simulations in OpenMM, why are you going through Gromacs files? CHARMM-GUI can generate OpenMM inputs directly.

Peter

[George Pantelopulos]

Hi Peter,

I added the if else definitions for flexible and rigid waters to the topology and waters appear to be behaving properly now!

I had initially tried using CHARMM-GUI's OpenMM scripts but I was met with a lot of errors with the input files they use... I also tried simply using the CHMARMM files to do this, and was met with some errors regarding the topology files as well. To me, this was the path of least resistance since I have more familiarity with GROMACS.

Thank you again for your patience through this!