Getting Started with MELD
Running your first REMD run
In essence MELD provides a customizable way of running REMD. We can decide to change temperature, hamiltonian or both along the replica. In this tutorial you will learn the basic elements to run your first Replica exchange molecular Dynamics. We will create a setup.py script with all the python commands to setup a simulation.
Loading Some Helpful Modules
import numpy as np import meld from meld import unit as u
Setting up the system
You can setup your protein either from a pdb file or from sequence.
# Option 1: from a pdb file protein = meld.AmberSubSystemFromPdbFile("example.pdb") # Option 2: from sequence in a file protein = meld.AmberSubSystemFromSequence("NALA ALA CALA") # Option 3: from a sequence sequence = meld.get_sequence_from_AA1(filename='sequence.dat') protein = meld.AmberSubSystemFromSequence(sequence)
Once we have the protein system we have to specify a force field. Current
ff14sb. We can also chose the
implicit solvent model, in this case
We also need to setup various options about how the simulation will be run, including
the timestep and any nonbonded cutoffs. By setting
use_big_testep, we will be using
a 3.5 fs timestep with hydrogen mass repartitioning.
build_options = meld.AmberOptions( forcefield="ff14sbside", implicit_solvent_model="gbNeck2", use_big_timestep=True, cutoff=1.8*u.nanometer ) builder = system.AmberSystemBuilder(build_options)
Now we generate the system:
system = builder.build_system([protein]).finalize()
At this point we can start defining different ways to setup the replica ladders.
In MELD we have a parameter called alpha, with values between [0,1] which map on
to the replica ladder. The lowest replica will always have
the highest replica will have
alpha=1.0. The intermediate replcias will
initial linearly interpolate between these values and will be modified
adaptively during the simualtion to attain approximately equal average
acceptance probabilities between adjacent replicas.
The various properties of the system, e.g. temperatures or force constants, are
set based on
alpha through the use of scalers.
We will use a geometric scaling of the temperatures. As an initial example let us setup a temparature replica ladder that will cover the whole replica space going from 300K to 450K.
system.temperature_scaler = meld.GeometricTemperatureScaler(0, 1.0, 300., 450.)
Next, we will specify some options for the system. The
controls the number of molecular dynamics steps between exchanges. Since we set
use_big_timestep the value of
14_286 corresponds to 50 ps between exchanges.
The system will be minimized for
minimize_steps before the simulation starts.
options = meld.RunOptions( timesteps = 14_286 minimize_steps = 20_000, )
Setting up replica exchange
Next, we need to setup replica exchange. The main parameters to set are: the
total number of rounds of simulation to run
n_steps, the number of replicas
n_replicas, the number of trials per replica exchange step
n_replicas**2). There are also various options for setting
alpha, but the defaults should generally work fine.
In this case, we’re doing 5000 steps of replica exchange, each with 50 ps of molecular dynamics, for a total of 250 ns.
remd = meld.setup_replica_exchange(system, n_replicas=4, n_steps=5000)
Storing simulation output
Finally, we need to store all of this setup information to disk in preparation for running the calculation
meld.setup_data_store(system, options, remd)
Full setup script
The full setup script should be saved as something like
import numpy as np import meld from meld import unit as u protein = meld.AmberSubSystemFromSequence("NALA ALA CALA") build_options = meld.AmberOptions( forcefield="ff14sbside", implicit_solvent_model="gbNeck2", use_big_timestep=True, cutoff=1.8 * u.nanometer ) builder = system.AmberSystemBuilder(build_options) system = builder.build_system([protein]).finalize() system.temperature_scaler = meld.GeometricTemperatureScaler(0, 1.0, 300. * u.kelvin, 450. * u.kelvin) options = meld.RunOptions( timesteps = 14_286 minimize_steps = 20_000, ) remd = meld.setup_replica_exchange(system, n_replicas=4, n_steps=5000) meld.setup_data_store(system, options, remd)
Running the system
python setup_simulation.py you should get a Data directory
with all the files needed to run MELD. Use your queing system to submit an MPI
job with the number of replicas you have indicated. Currently, we need one GPU
for each replica.
aprun -n 2 -N 1 launch_remd --debug
For debugging purposes, it is also possible to run locally on a laptop or workstation, although this can be very slow.
launch_remd_multiplex --debug --platform CUDA
You can replace
Reference if you do not have a GPU, athough this
will be quite slow.