Minimization not updating positions

[Pujan Ajmera]

Hi all,

I am attempting to minimize the lipids in my structure (this is part 1 of a 5 step minimization I'm doing). For some reason, the energies aren't really changing, and the coordinates are barely changing - my initial system is definitely not close to a minimum, but this did run for around an hour which makes me thing something was happening. Please let me know if I am doing something wrong in my code. Here are the energies:

Initial Energy: 4.182592537913685e+20 kJ/mol
Minimized Energy: 4.182592537913686e+20 kJ/mol

And here is my script:

Code: Select all

from openmm import *
from openmm.app import *
from openmm.unit import *
import time

# Input Files

psf = CharmmPsfFile('step5_assembly.psf')
psf.setBox(251.166872*angstroms,251.166872*angstroms,192.089708*angstroms)
crd = CharmmCrdFile('step5_assembly.crd')
params = CharmmParameterSet('toppar_all36_lipid_archaeal.str', 'f4s4.prm', 'f4s4.rtf', 'par_all36m_prot.prm', 'toppar_dum_noble_gases.str', 'cam.str', 'toppar_all36_lipid_cardiolipin.str', 'top_all36_prot.rtf', 'toppar_all36_prot_modify_res.str', 'par_interface.prm', 'par_all36_carb.prm', 'toppar_water_ions.str', 'par_all36_na.prm', 'toppar_all36_prot_model.str', 'toppar_ions_won.str', 'toppar_all36_prot_retinol.str', 'par_all36_lipid.prm', 'top_all36_lipid.rtf', 'toppar_all36_prot_c36m_d_aminoacids.str', 'top_interface.rtf','top_all36_carb.rtf', 'top_all36_na.rtf', 'par_all36_cgenff.prm', 'toppar_all36_prot_arg0.str', 'top_all36_cgenff.rtf', 'step5_assembly.str')

# System Configuration

nonbondedMethod = PME
nonbondedCutoff = 1.0*nanometers
ewaldErrorTolerance = 0.0005
constraints = HBonds
rigidWater = True
constraintTolerance = 0.000001
hydrogenMass = 1.5*amu

# Integration Options

dt = 0.002*picoseconds
temperature = 310*kelvin
friction = 1.0/picosecond
pressure = 1.0*atmospheres
barostatInterval = 25

# Simulation Options

platform = Platform.getPlatformByName('CUDA')
platformProperties = {'Precision': 'double'}

print('Building system...')
sys.stdout.flush()
time.sleep(1.0)
topology = psf.topology
positions = crd.positions

system_lipids = psf.createSystem(params, nonbondedMethod=nonbondedMethod, nonbondedCutoff=nonbondedCutoff,
                                 constraints=constraints, rigidWater=rigidWater, ewaldErrorTolerance=ewaldErrorTolerance, hydrogenMass=hydrogenMass)
for atom in topology.atoms():
    if atom.residue.name not in ('POPC', 'POPE', 'TYCL2'):
        system_lipids.setParticleMass(atom.index, 0)
integrator = LangevinMiddleIntegrator(temperature, friction, dt)
integrator.setConstraintTolerance(constraintTolerance)
simulation = Simulation(topology, system_lipids, integrator, platform, platformProperties)
simulation.context.setPositions(positions)
initial_energy = simulation.context.getState(getEnergy=True).getPotentialEnergy()
print(f"Initial Energy: {initial_energy}, minimizing lipids")
sys.stdout.flush()
time.sleep(1.0)
simulation.minimizeEnergy()
minimized_energy = simulation.context.getState(getEnergy=True).getPotentialEnergy()
print(f"Minimized Energy: {minimized_energy}")
sys.stdout.flush()
time.sleep(1.0)
if minimized_energy == initial_energy:
    print("No change in energy, minimization might not have moved the atoms.")
simulation.saveState(f'minimized_state_lipids.xml')
positions = simulation.context.getState(getPositions=True).getPositions()
with open(f'minimized_positions_lipids.pdb', 'w') as f:
    PDBFile.writeFile(simulation.topology, positions, f)

[Peter Eastman]

That means the minimization failed. It tried valiantly for an hour to find a lower energy conformation and then gave up.

That can happen when the starting conformation has really serious problems, which seems to be the case in your system given the insanely high energy. You need to examine your system and figure out where that energy is coming from. It likely means you have some molecules right on top of each other.

[Pujan Ajmera]

Hi Peter,

Thank you for your response. I did an energy decomposition analysis on my initial structure, and also tried letting parts of the system go iteratively to see which was causing the overly large energies:

Initial structure energy decomposition:
Force HarmonicBondForce: 697298.8917322994 kJ/mol
Force HarmonicAngleForce: 430281.2467599694 kJ/mol
Force HarmonicBondForce: 561974.1027477044 kJ/mol
Force PeriodicTorsionForce: 153909.5221186617 kJ/mol
Force CustomTorsionForce: 4388.548054359069 kJ/mol
Force CMAPTorsionForce: 2713.8386978680783 kJ/mol
Force NonbondedForce: 4.1825925379136664e+20 kJ/mol
Force CustomNonbondedForce: 0.0 kJ/mol
Force CMMotionRemover: 0.0 kJ/mol

Energies after each minimization step:
Initial Energy: 4.182592537913685e+20 kJ/mol, minimizing lipids
After minimizing lipids, freezing else: 4.182592537913686e+20 kJ/mol
After minimizing lipids, waters, and ions, freezing protein+substrate+cofactor: 4.182592537913684e+20 kJ/mol
After minimizing everything but substrate+cofactor: 4.182592537841562e+20 kJ/mol
After minimizing everything: -15811651.106276177 kJ/mol

So it seems that the substrates and cofactors are causing some insanely high non-bonded force, which is odd because in my structure they are not overlapping with anything else. I am seeing that my f4s4 clusters (the cofactors) are getting distorted during this last minimization process significantly, but my substrate seems stable - it seems that those must be the cause of high energy then, but again I don't see any overlap of atoms. Do you know if anything else in my parameters can be causing this high nonbonded force?

[Peter Eastman]

Try querying the forces and see which particular atoms have exceptionally high forces. That will tell you where the problem is.