simtk.org

Code: Select all

#always include at the top of simulation
import openmm.app as app
import openmm as mm
import openmm.unit as unit
import os
os.environ['JAX_ENABLE_X64'] = '1' #set JAX to 64bits by default
import openmmtools as tools
from openmmtools import integrators
from openmmtools import forces #still under construction, check for updates occasionally
from openmmtools import alchemy
from sys import stdout
from time import perf_counter as pf
import MDAnalysis as mda
import numpy as np

start = pf()

pdb = app.PDBFile('ZGY_1p5nmPadWaterCube.pdb')
forcefield = app.ForceField('./amber99sb-ZGY-hack.xml','/home/leo/.conda/envs/openmmEnv/lib/python3.10/site-packages/openmm/app/data/tip3p.xml')
#forcefield = app.ForceField('/home/hall2401/ForceFields/amber99sb-ZGY-hack.xml','/home/hall2401/ForceFields/tip3p.xml')

# System Configuration
nonbondedMethod = app.PME
nonbondedCutoff = 1.2*unit.nanometers
ewaldErrorTolerance = 0.0005
constraints = None
lambda_value = 0.8

# Integration Options
dt = 0.0005*unit.picoseconds # 0.5 femtoseconds
temperature = 298*unit.kelvin 
friction = 1.0/unit.picosecond
pressure = 1.0*unit.atmospheres
barostatInterval = 25

# Simulation Options
steps = 20000#00 #1ns with the 0.5fs timestep
equilibrationSteps = 10000#000 #0.5ns equilibration
CPUplatform = mm.Platform.getPlatformByName('CPU')
GPUplatform = mm.Platform.getPlatformByName('CUDA')
GPUplatformProperties = {'Precision': 'single'}
pdbReporter = app.PDBReporter('trajectory-solvation_CONCGly_CONCNaCl_TEMP_LAMBDA_RUN0p1.pdb', 1000)
dataReporter = app.StateDataReporter('trajectory-solvation_CONCGly_CONCNaCl_TEMP_LAMBDA_RUN0p1.csv', 1000, totalSteps=steps,
    step=True, speed=True, progress=True, potentialEnergy=True, temperature=True, separator=',')
checkpointReporter = app.CheckpointReporter('checkpoint-solvation_CONCGly_CONCNaCl_TEMP_LAMBDA_RUN0p1.chk', 10000)

# Create the system
print('Building system...')
topology = pdb.topology
positions = pdb.positions
system = forcefield.createSystem(topology=topology, nonbondedMethod=nonbondedMethod, nonbondedCutoff=nonbondedCutoff,
    constraints=constraints, ewaldErrorTolerance=ewaldErrorTolerance)

#Get residue of interest
target_residue_index = 0  # Index of the target glycine residue
target_atoms = [atom.index for atom in pdb.topology.atoms() if atom.residue.index == target_residue_index]

#create system where target atoms are alchemically modifiable
alchemical_region = alchemy.AlchemicalRegion(alchemical_atoms=target_atoms)
factory = alchemy.AbsoluteAlchemicalFactory()
alchemical_system = factory.create_alchemical_system(system,alchemical_region)

integrator = integrators.LangevinIntegrator(temperature, friction, dt)
#integrator = mm.LangevinMiddleIntegrator(temperature, friction, dt)
simulation = app.Simulation(topology, alchemical_system, integrator, GPUplatform, GPUplatformProperties)
simulation.context.setPositions(positions)


print("performing energy minimisation...")
simulation.minimizeEnergy()
simulation.context.setVelocitiesToTemperature(temperature)

print("equilibrating system...")
simulation.reporters.append(dataReporter)
simulation.step(equilibrationSteps)

print("applying alchemical transformation")
#apply the Lambda parametre to the system
alchemical_state = alchemy.AlchemicalState.from_system(alchemical_system)
alchemical_state.lambda_electrostatics = 0.9
alchemical_state.lambda_sterics = 0.9
#alchemical_state.set_alchemical_parameters(lambda_value)
alchemical_state.apply_to_context(simulation.context)

print("equilibrating system again...")
#simulation.reporters.append(dataReporter)
simulation.currentStep = 0
simulation.step(equilibrationSteps)

print("applying alchemical transformation")
#apply the Lambda parametre to the system
alchemical_state = alchemy.AlchemicalState.from_system(alchemical_system)
alchemical_state.lambda_electrostatics = lambda_value
alchemical_state.lambda_sterics = lambda_value
#alchemical_state.set_alchemical_parameters(lambda_value)
alchemical_state.apply_to_context(simulation.context)

# Simulate
print('Simulating...')
simulation.reporters.append(pdbReporter)
#simulation.reporters.append(dataReporter) #already added during equilibration
simulation.reporters.append(checkpointReporter)
simulation.currentStep = 0
simulation.step(steps)

end = pf()
print('runtime (mn) = ',(end-start)/60)

Code: Select all

#always include at the top of simulation
import openmm.app as app
import openmm as mm
import openmm.unit as unit
import os
os.environ['JAX_ENABLE_X64'] = '1' #set JAX to 64bits by default
import openmmtools as tools
from openmmtools import integrators
from openmmtools import forces #still under construction, check for updates occasionally
from openmmtools import alchemy
from sys import stdout
from time import perf_counter as pf
import MDAnalysis as mda
import numpy as np

start = pf()
LAMBDA = 10 #LAMBDA
TEMP = 298#TEMP
pdb = app.PDBFile('ZGY_1p5nmPadWaterCube.pdb')
forcefield = app.ForceField('./amber99sb-ZGY-hack.xml','/home/leo/.conda/envs/openmmEnv/lib/python3.10/site-packages/openmm/app/data/tip3p.xml')
#forcefield = app.ForceField('/home/hall2401/ForceFields/amber99sb-ZGY-hack.xml','/home/hall2401/ForceFields/tip3p.xml')

# System Configuration
nonbondedMethod = app.PME
nonbondedCutoff = 1.2*unit.nanometers
ewaldErrorTolerance = 0.0005
constraints = None
lambda_value = LAMBDA/100 #LAMBDA will be given as a percentage for purpose of filenames

# Integration Options
dt = 0.0005*unit.picoseconds # 0.5 femtoseconds
temperature = TEMP*unit.kelvin 
friction = 1.0/unit.picosecond
pressure = 1.0*unit.atmospheres
barostatInterval = 25

# Simulation Options
steps = 20000#00 #1ns with the 0.5fs timestep
equilibrationSteps = 10000#000 #0.5ns equilibration
CPUplatform = mm.Platform.getPlatformByName('CPU')
GPUplatform = mm.Platform.getPlatformByName('CUDA')
GPUplatformProperties = {'Precision': 'single'}
pdbReporter = app.PDBReporter('trajectory-solvation_CONCGly_CONCNaCl_TEMP_LAMBDA_RUN0p1.pdb', 1000)
dataReporter = app.StateDataReporter('trajectory-solvation_CONCGly_CONCNaCl_TEMP_LAMBDA_RUN0p1.csv', 1000, totalSteps=steps,
    step=True, speed=True, progress=True, potentialEnergy=True, temperature=True, separator=',')
checkpointReporter = app.CheckpointReporter('checkpoint-solvation_CONCGly_CONCNaCl_TEMP_LAMBDA_RUN0p1.chk', 10000)

# Create the system
print('Building system...')
topology = pdb.topology
positions = pdb.positions
system = forcefield.createSystem(topology=topology, nonbondedMethod=nonbondedMethod, nonbondedCutoff=nonbondedCutoff,
    constraints=constraints, ewaldErrorTolerance=ewaldErrorTolerance)

#Get residue of interest
target_residue_index = 0  # Index of the target glycine residue
target_atoms = [atom.index for atom in pdb.topology.atoms() if atom.residue.index == target_residue_index]

integrator = mm.LangevinMiddleIntegrator(temperature, friction, dt)
simulation = app.Simulation(topology, system, integrator, GPUplatform, GPUplatformProperties)

simulation.context.setPositions(positions)


# Get the NonbondedForce object
#we get the first element of a list that presumably contains just one element as there should be only one nonbonded force object
nonbonded_force = [force for force in system.getForces() if isinstance(force, mm.NonbondedForce)][0] 
    
# Lambda parametre already defined
#dielectric constant is not considered as we're doing an explicit solvent simulation
# Create a CustomNonbondedForce that scales both Lennard-Jones and Coulombic interactions
custom_force = mm.CustomNonbondedForce(
    'lambda*(4*epsilon*((sigma/r)^12 - (sigma/r)^6) + 138.935456*charge1*charge2/(1/r)); ' 
    'sigma=0.5*(sigma1+sigma2); epsilon=sqrt(epsilon1*epsilon2)')
#138.9... is the value of the Coulomb constant in Openmm in kJ/nm
custom_force.addGlobalParameter('lambda', lambda_value)
#charge, sigma and epsilon are named parametres whose value will be input for each relevant particle
custom_force.addPerParticleParameter('charge')
custom_force.addPerParticleParameter('sigma')
custom_force.addPerParticleParameter('epsilon')
#give it the same parametres as other nonbonded forces
custom_force.setCutoffDistance(nonbondedCutoff)
custom_force.setNonbondedMethod(mm.CustomNonbondedForce.CutoffPeriodic)

# Add particles to the custom force using parameters from the existing force field
#the whole system must be included in the custom force because the residue can interact with the whole system
for i in range(system.getNumParticles()):
    charge, sigma, epsilon = nonbonded_force.getParticleParameters(i)
    custom_force.addParticle([charge, sigma, epsilon])
    
    # Enable lambda scaling for interactions involving the target glycine residue
    if i in target_atoms:
        #setParticleParametre means interactions will be calculated for this particle
        custom_force.setParticleParameters(i, [charge, sigma, epsilon])
        nonbonded_force.setParticleParameters(i, 0, 0, 0) #remove target atoms from the nonbonded force

# Add exclusions to avoid double-counting interactions that are already included in bonded interactions
#the residue presumably only interacts with itself through bonded interactions
# Extract bonds from the topology and convert to the required format
bonds = [(bond[0].index, bond[1].index) for bond in pdb.topology.bonds()]

# Add exclusions for 1-2 interactions (atoms directly bonded to each other - this is what the second argument is for)
custom_force.createExclusionsFromBonds(bonds, 1)

# Add the custom force to the system
system.addForce(custom_force)

# Disable interactions in the default NonbondedForce for target atoms
#for atom_index in target_atoms:
 #   for other_index in range(system.getNumParticles()):
        #if other_index not in target_atoms: #interactions between atoms from the target should also be excluded as they're taken care of by the custom force
  #      nonbonded_force.addException(atom_index, other_index, 0.0, 0.0, 0.0,replace=False) #there might be exceptions already for bonded atoms




print("performing energy minimisation...")
simulation.minimizeEnergy()
simulation.context.setVelocitiesToTemperature(temperature)

print("equilibrating system...")
simulation.reporters.append(dataReporter)
simulation.step(equilibrationSteps)


# Simulate
print('Simulating...')
simulation.reporters.append(pdbReporter)
#simulation.reporters.append(dataReporter) #already added during equilibration
simulation.reporters.append(checkpointReporter)
simulation.currentStep = 0
simulation.step(steps)

end = pf()
print('runtime (mn) = ',(end-start)/60)

Code: Select all

for i in range(nonbonded_force.getNumExceptions()):
    params = nonbonded_force.getExceptionParameters(i)
    custom_force.addExclusion(params[0], params[1])

Code: Select all

nonbonded_force.setParticleParameters(i, 0, 0, 0) #remove target atoms from the nonbonded force

Code: Select all

#always include at the top of simulation
import openmm.app as app
import openmm as mm
import openmm.unit as unit
import os
os.environ['JAX_ENABLE_X64'] = '1' #set JAX to 64bits by default
from time import perf_counter as pf
import MDAnalysis as mda
import numpy as np

start = pf()
LAMBDA = 10 #LAMBDA
TEMP = 298#TEMP
pdb = app.PDBFile('ZGY_1p5nmPadWaterCube.pdb')
forcefield = app.ForceField('./amber99sb-ZGY-hack.xml','/home/leo/.conda/envs/openmmEnv/lib/python3.10/site-packages/openmm/app/data/tip3p.xml')
#forcefield = app.ForceField('/home/hall2401/ForceFields/amber99sb-ZGY-hack.xml','/home/hall2401/ForceFields/tip3p.xml')

# System Configuration
nonbondedMethod = app.PME
nonbondedCutoff = 1.2*unit.nanometers
ewaldErrorTolerance = 0.0005
constraints = None
lambda_value = LAMBDA/100 #LAMBDA will be given as a percentage for purpose of filenames

# Integration Options
dt = 0.0005*unit.picoseconds # 0.5 femtoseconds
temperature = TEMP*unit.kelvin 
friction = 1.0/unit.picosecond
pressure = 1.0*unit.atmospheres
barostatInterval = 25

# Simulation Options
steps = 20000#00 #1ns with the 0.5fs timestep
equilibrationSteps = 10000#00 #0.5ns equilibration
CPUplatform = mm.Platform.getPlatformByName('CPU')
GPUplatform = mm.Platform.getPlatformByName('CUDA')
GPUplatformProperties = {'Precision': 'single'}
pdbReporter = app.PDBReporter('trajectory-solvation_water_TEMP_LAMBDA_RUN.pdb', 1000)
dataReporter = app.StateDataReporter('trajectory-solvation_water_TEMP_LAMBDA_RUN.csv', 1000, totalSteps=steps,
    step=True, speed=True, progress=True, potentialEnergy=True, temperature=True, separator=',')
checkpointReporter = app.CheckpointReporter('checkpoint-solvation_water_TEMP_LAMBDA_RUN.chk', 10000)

# Create the system
print('Building system...')
topology = pdb.topology
positions = pdb.positions
system = forcefield.createSystem(topology=topology, nonbondedMethod=nonbondedMethod, nonbondedCutoff=nonbondedCutoff,
    constraints=constraints, ewaldErrorTolerance=ewaldErrorTolerance)

#Get residue of interest
target_residue_index = 0  # Index of the target glycine residue
target_atoms = [atom.index for atom in pdb.topology.atoms() if atom.residue.index == target_residue_index]

integrator = mm.LangevinMiddleIntegrator(temperature, friction, dt)
simulation = app.Simulation(topology, system, integrator, GPUplatform, GPUplatformProperties)

simulation.context.setPositions(positions)


# Get the NonbondedForce object
#we get the first element of a list that presumably contains just one element as there should be only one nonbonded force object
nonbonded_force = [force for force in system.getForces() if isinstance(force, mm.NonbondedForce)][0] 
    
# Lambda parametre already defined
#dielectric constant is not considered as we're doing an explicit solvent simulation
# Create a CustomNonbondedForce that scales both Lennard-Jones and Coulombic interactions
custom_force = mm.CustomNonbondedForce(
    'lambda*(4*epsilon*((sigma/r)^12 - (sigma/r)^6) + 138.935456*charge1*charge2/(1/r)); ' 
    'sigma=0.5*(sigma1+sigma2); epsilon=sqrt(epsilon1*epsilon2)')
#138.9... is the value of the Coulomb constant in Openmm in kJ/nm
custom_force.addGlobalParameter('lambda', lambda_value)
#charge, sigma and epsilon are named parametres whose value will be input for each relevant particle
custom_force.addPerParticleParameter('charge')
custom_force.addPerParticleParameter('sigma')
custom_force.addPerParticleParameter('epsilon')
#give it the same parametres as other nonbonded forces
custom_force.setCutoffDistance(nonbondedCutoff)
custom_force.setNonbondedMethod(mm.CustomNonbondedForce.CutoffPeriodic)

# Add particles to the custom force using parameters from the existing force field
#the whole system must be included in the custom force because the residue can interact with the whole system
for i in range(system.getNumParticles()):
    charge, sigma, epsilon = nonbonded_force.getParticleParameters(i)
    custom_force.addParticle([charge, sigma, epsilon]) #when called for the i'th time, sets the parametres for the i'th particle
    
    # Enable lambda scaling for interactions involving the target glycine residue
    if i in target_atoms:
        #setParticleParametre means interactions will be calculated for this particle
        custom_force.setParticleParameters(i, [charge, sigma, epsilon])
        nonbonded_force.setParticleParameters(i, 0, 0, 0) #remove target atoms from the nonbonded force

# Add exclusions to avoid double-counting interactions that are already included in bonded interactions
#the residue presumably only interacts with itself through bonded interactions
# Extract bonds from the topology and convert to the required format
bonds = [(bond[0].index, bond[1].index) for bond in pdb.topology.bonds()]

# Add exclusions for 1-2 interactions (atoms directly bonded to each other - this is what the second argument is for)
custom_force.createExclusionsFromBonds(bonds, 1)

# Add the custom force to the system
system.addForce(custom_force)


print("performing energy minimisation...")
simulation.minimizeEnergy()
simulation.context.setVelocitiesToTemperature(temperature)

print("equilibrating system...")
simulation.reporters.append(dataReporter)
simulation.step(equilibrationSteps)


# Simulate
print('Simulating...')
simulation.reporters.append(pdbReporter)
#simulation.reporters.append(dataReporter) #already added during equilibration
simulation.reporters.append(checkpointReporter)
simulation.currentStep = 0
simulation.step(steps)