I might be confused but once I add Water into my system it seems that it has reached a stable state. I keep adding heat over some time but the molecules just move in their positions.
I am setting the lipids in a spherical formation with two layers (inner and outer). It seems that they break apart and are stuck in the four corners of the system.
I am attaching the code for examination or replication if you are curious. I am also attaching some screenshots to provide some visual guidance on what I mean. (one is constructing the molecule, the second after minimization and adding solvent, and the last after running 1m simulation steps)
Do you have any suggestions on what am I doing wrong?
Code: Select all
from openmm import *
from openmm.app import *
from openmm.unit import *
import MDAnalysis as mda
# Input Files
pdb = PDBFile('lipids/HMDB0000067_Cholesterol_3D.pdb')
forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')
modeller = Modeller(pdb.topology, pdb.positions)
modeller.addSolvent(forcefield=forcefield, model='tip3p', padding=3.0*nanometers, ionicStrength=0.1*molar)
# 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.004*picoseconds
temperature = 300*kelvin
friction = 1.0/picosecond
pressure = 1.0*atmospheres
barostatInterval = 25
# Simulation Options
steps = 1_000_000
equilibrationSteps = 1000
platform = Platform.getPlatformByName('CUDA')
platformProperties = {'Precision': 'single'}
dcdReporter = DCDReporter('trajectory.dcd', 10000)
dataReporter = StateDataReporter('log.txt', 1000, totalSteps=steps,
step=True, speed=True, progress=True, potentialEnergy=True, temperature=True, separator='\t')
energyReporter = StateDataReporter(
'energy.csv',
100,
step=True,
potentialEnergy=True,
kineticEnergy=True,
totalEnergy=True,
temperature=True
)
checkpointReporter = CheckpointReporter('checkpoint.chk', 10000)
# Prepare the Simulation
print('adding solvent')
print('Building system...')
system = forcefield.createSystem(modeller.topology, nonbondedMethod=nonbondedMethod, nonbondedCutoff=nonbondedCutoff,
constraints=constraints, rigidWater=rigidWater, ewaldErrorTolerance=ewaldErrorTolerance, hydrogenMass=hydrogenMass)
system.addForce(MonteCarloBarostat(pressure, temperature, barostatInterval))
integrator = LangevinMiddleIntegrator(temperature, friction, dt)
integrator.setConstraintTolerance(constraintTolerance)
simulation = Simulation(modeller.topology, system, integrator, platform, platformProperties)
simulation.context.setPositions(modeller.positions)
# Minimize and Equilibrate
print('Performing energy minimization...')
simulation.minimizeEnergy()
print('Equilibrating...')
simulation.context.setVelocitiesToTemperature(temperature)
simulation.step(equilibrationSteps)
# Simulate
print('Simulating...')
simulation.reporters.append(dcdReporter)
simulation.reporters.append(dataReporter)
simulation.reporters.append(checkpointReporter)
simulation.reporters.append(energyReporter)
simulation.currentStep = 0
simulation.step(steps)
state = simulation.context.getState(getPositions=True, enforcePeriodicBox=system.usesPeriodicBoundaryConditions())
with open("final_state.pdb", mode="w") as file:
PDBFile.writeFile(simulation.topology, state.getPositions(), file)
def plot_energy():
import pandas as pd
import matplotlib.pyplot as plt
# Load the energy data
data = pd.read_csv('energy.csv')
# Plotting
plt.figure(figsize=(10, 6))
# Plot potential energy
plt.plot(data['#"Step"'], data['Potential Energy (kJ/mole)'], label='Potential Energy')
# Plot kinetic energy
plt.plot(data['#"Step"'], data['Kinetic Energy (kJ/mole)'], label='Kinetic Energy')
# Plot total energy
plt.plot(data['#"Step"'], data['Total Energy (kJ/mole)'], label='Total Energy')
plt.xlabel('Step')
plt.ylabel('Energy (kJ/mole)')
plt.title('Energy vs. Simulation Step')
plt.legend()
plt.show()
def write_trajectory_vid():
u = mda.Universe('final_state.pdb', 'trajectory.dcd')
with mda.Writer("trajectory_trajectory.pdb", multiframe=True, bonds=None, n_atoms=u.atoms.n_atoms) as W:
for ts in u.trajectory:
W.write(u.atoms)
write_trajectory_vid()
plot_energy()