Minimization not updating positions
Posted: Tue Sep 03, 2024 9:46 pm
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:
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)