Skip to content

Na Cl Association with AMBER 16

atbogetti edited this page Sep 3, 2018 · 27 revisions

Table of Contents

AMBER Tutorial: Molecular Dynamics of Na+/Cl- Association

by Alex DeGrave, Karl Debiec, Daniel Hatfield

Overview

Requirements: ~36 hr wallclock time on a 4-core, 3.0 GHz Intel Xeon node (one walker per core); ~18 GB disk space

In this introductory tutorial we will use the weighted ensemble strategy to simulate Na+/Cl- association in TIP3P explicit solvent. The system consists of single Na+ and Cl- ions modeled with the Joung and Cheatham ion parameters[1], using the distance between the two ions as the progress coordinate. Amber will be used to run the molecular dynamics, and familiarity with it is a prerequisite (see tutorials). Basic knowledge of python and bash scripting is also necessary.

Choosing Weighted Ensemble Parameters

Running a successful weighted ensemble (WE) simulation involves the choice of a number of parameters. In addition to usual simulation concerns, such as what force field to use, WE simulations depend on one or more progress coordinates, the resampling interval, the binning scheme, and the number of trajectories per bin.

  • Progress coordinate: A coordinate that tracks the progress of the system toward one or more target states. Examples could include the distance between two binding partners, an RMSD from a reference structure, the radius of gyration, or a backbone dihedral angle. More than one progress coordinate may be combined, if desired.
  • Resampling interval: Commonly denoted τ, the resampling interval indicates the time between combination/replication events.
  • Binning scheme: The binning scheme defines a partition of the progress coordinate(s). The WE algorithm performs combination/replication events to achieve some target number of trajectories per bin, and trajectores that are in the same bin may be combined. Commonly, bins are spaced more finely on steep energetic barriers, in order to give the system more chances to climb over the barrier.
  • Number of trajectories per bin: The WE algorithm performs replication/combination events to maintain this number of trajectories.
Thus, a common question when setting up a WE simulation is: what parameters should I choose? While automated procedures are under development, good parameters can often be chosen by using one's chemical intuition. For example, since we wish to simulate the association of Na+ and Cl-, a natural choice of progress coordinate is the distance between the ions. One might imagine that the energetic barriers to association become steepest as the ions displace water molecules of the first and second hydration shells, so we will use a binning scheme that becomes finer as the ions approach each other. Given the size of the ions, we also have an idea of the distance (in progress-coordinate space) the ions can move per unit time, giving us a starting place for the choice of the resampling interval (e.g., something on the order of ps); ideally, most trajectories will not move through more than one or two bins per iteration.

Since we are only concerned with the association of Na+ and Cl-, we will perform a steady-state simulation, in which trajectories are recycled after reaching the bound state. That is, trajectories that reach the bound state will be restarted from our initial (unbound) state, enabling us to focus computational power on the direction of interest. To recycle simulations at the bound state, we need to define the bound state before starting the WE simulation.

Defining the bound state

For some systems, the definition of a target state may be straightforward. For example, if a crystal structure is available for the bound state of a protein complex, we might say that any conformations within reasonable deviation of that crystal structure are bound. In our case, the position of the bound free energy basin may be sensitive to the choice of ion parameters, so we will determine the position of the basin using a method that takes the specific ion parameters into account. After the weighted ensemble simulation, the definition of the target state may be updated, as long as the new definition is a superset of the target state used for recycling (i.e., the new definition must encompass all of the definition used for recycling, and may be larger). Therefore, it is a good rule of thumb to go with the most stringent definition of the target state. In our case, we will estimate the position of the target state by minimizing the energy of Na+ and Cl- ions placed next to each other; following the minimization, the attraction of the electrostatic interaction should balance out repulsion from the van der Waals term.

In your WESTPA installation directory, referred to from hereon as $WEST_ROOT, find the directory nacl_amb in lib/examples. This directory contains all the files needed for this tutorial. In order to keep the installation clean, make a copy of nacl_amb somewhere outside of your WESTPA installation. Inside nacl_amb, you will find a subdirectory named bound_state_definition. This directory contains the scripts and files for our minimization:

  • 1_leap: Here, we prepare initial coordinates for the system. As a guess of the location of the bound state, we start with the Na+ and Cl- ions positioned at bottom of the well of the Lennard-Jones potential between the ions, which can be deduced from the ion parameters specified by Joung and Cheatham.[1] If you haven't already, change directories to 1_leap, and after making sure that Amber is in your $PATH, run the script run.sh to build a topology file and solvate the system.
  • 2_min: Now we minimize the system. Change directories to 2_min and run run.sh to perform an energy minimization. To find the distance between the two ions, call analyze.sh. The script uses cpptraj to calculate the distance between the centers of the ions (in Å), and prints it to the terminal. Take note of this value, as you will need it when preparing the weighted ensemble configuration.

Preparing the WESTPA Configuration

Now we are ready to decide upon the weighted ensemble parameters, which are specified in the configuration file west.cfg.

Begin by examining the file west.cfg, which is found in the nacl_amb directory.

This file defines the WESTPA configuration, including:

  • the number of dimensions of the progress coordinate
  • the number of data points per trajectory segment
  • the binning scheme and number of trajectories per bin
  • the number of weighted ensemble iterations to run
  • the locations of various scripts that WESTPA needs to execute the simulation
For this system, the progress coordinate will be the distance between the Na+ and Cl- ions, so the progress coordinate is one-dimensional. For the binning scheme, we will place bin boundaries between the position of bound (target) state, at which trajectories will be recycled, and the unbound state. Here, we define the unbound state as >10 Å separation, which is (somewhat arbitrarily) taken to be the point at which the ions no longer interact to an appreciable extent. We space bin boundaries more finely near the bound state to help the system climb the energetic barriers imposed by the hydration shells. We will start by using 5 trajectories per bin; this gives the system a few chances to progress to a new bin during each iteration, while not requiring an undue amount of computational power. While we do not need to set the resampling interval τ--it is determined by the input files passed to the dynamics engine--it is helpful to consider now what this interval will be. For this tutorial, we choose 5 ps, which is long enough for some but usually not all trajectories to progress to a new bin. During each iteration we will collect 51 datapoints (50 datapoints plus the initial point, which should be the same as the final timepoint from the previous iteration), giving our data 100 fs time resolution.

The file west.cfg already contains an example set of bin boundaries. Update the boundaries to reflect the value you recorded after minimizing the system in the previous section of the tutorial. Your value should be similar to the value already in west.cfg.

The binning scheme used for the weighted ensemble simulation. Na+ is depicted in red, and Cl- is depicted in blue.

Note that our bins are based loosely on those defined by Zwier, Kaus, and Chong[2]; for more information on choosing weighted ensemble parameters, you may find it helpful to read this paper.

Preparing the system for simulation

Before starting the weighted ensemble simulation, we will prepare a system that includes Na+, Cl-, and TIP3P waters, with the ions in the "unbound" state. Following minimization and equilibration, this system will become a "basis state" for the WESTPA simulation, i.e., the system will be used to start new trajectories. If you wish to skip the minimization and equilibration, you may proceed directly to the next section; completed output files are provided in the example directory.

Scripts and input files for each step of preparing the system are provided in nacl_amb/prep:

  • 1_leap Begin by examining the files in 1_leap. The file nacl_no_solvent.pdb contains the initial atomic coordinates for the Na+ and Cl-, which are separated by 12 Å; such a file can be generated by hand or using a variety of modeling utilities. Call run.sh from the command line, which will use tleap to solvate the system and make a topology file.
  • 2_min Next we perform an energy minimization to prepare for equilibration. Change directories to 2_min and call run.sh to perform the minimization. Note that fairly stiff restraints are applied to hold the Na+ and Cl- ions near the desired distance.
  • 3_eq1 Now that we have a minimized system, we heat the system in a 20 ps NVT equilibration. Change directories to 3_eq1 and call run.sh to perform the equilibration. We apply a Langevin thermostat, which is the same as we will apply for the main weighted ensemble simulation.
  • 4_eq2 Next we perform a 1 ns NPT equilibration to bring the system to the desired pressure and relax the water molecules. Change directories to 4_eq2 and call run.sh to perform the equilibration. Given that the ions have no internal degrees of freedom, we will not perform an additional, unrestrained equilibration.
At this point, the system is equilibrated and ready for simulation with WESTPA. The necessary input files will be copied over from the prep directory when you call init.sh, which will be detailed below.

Preparing the Amber input file

Next, we will prepare an input file for the Amber dynamics program (pmemd or sander), which WESTPA will use to propagate the dynamics during the weighted ensemble simulation. The configuration file, md.in, is found in amber_config.

Several settings in this configuration file are important to WESTPA. First is the segment duration, τ, which is equal to the timestep times the number of steps. The appropriate τ depends on the system and progress coordinate, and for this system we have chosen 5.0 ps. For information on choosing τ, see the Guide to Choosing WE Parameters and the publication on which this tutorial is based.[2] The second setting of particular importance to WESTPA is the use of a stochastic thermostat; here we use the Langevin thermostat. WESTPA requires a stochastic element such that trajectories branched from a common point diverge. In order to ensure that this occurs, it is also necessary to set the random seed to a different value for each segment. Setting the random seed based on the time, as is common for brute-force simulations, may cause problems as simulations started simultaneously from the same parent will not diverge. For an explanation of other settings, please see the annotated version of the file, md.in.annotated.

Preparing for dynamics propagation

The job of WESTPA is to orchestrate the weighted ensemble simulation: WESTPA decides which trajectories need to be propagated, keeps track of the weight of each segment, keeps track of the relationships between segments, and stores auxiliary data. Calculations such as the propagation of dynamics and the calculation of progress coordinates is handled by the user. In this tutorial, we use the executable simulation manager, which enables users to perform these calculations through user-defined scripts. Working example scripts are provided in the example directory. Before delving into the specifics of the scripts, it is helpful to have an overview of how WESTPA uses the scripts. Following this overview, it is suggested to read through the individual scripts; additional comments are provided in the scripts.

During each weighted ensemble iteration, WESTPA produces a number of trajectory segments, each of which is a short simulation. Each trajectory segment can be either a new trajectory, or a continuation of a prior trajectory segment. For example, at the start of a WESTPA simulation, all trajectory segments are new. New trajectory segments may also arise after recycling events, in which a trajectory reaching the target state is terminated and a new trajectory is started from the initial state. Usually, the majority of trajectory segments after the start of a WESTPA simulation will be continuations of prior segments.

To run a trajectory segment of either type, WESTPA launches a user-specified script. Via environment variables, it supplies to the script information such as the location of a restart file representing the initial coordinates for that segment. The user-specified script must process this information, run the dynamics for the simulation segment, and return data to WESTPA (via environment variables). In this tutorial, the script that performs these functions is runseg.sh, found in westpa_scripts.

In the case of a new trajectory segments, three scripts are involved. To generate a new segment, WESTPA starts with a basis state. Before using this configuration as the starting point for a new trajectory segment, the user has the option to modify the configuration, forming an initial state from the basis state. For example, in a simulation of protein-protein association, one might consider starting with coordinates for the two binding partners (the basis state). An initial state could then be generated by randomly orienting the binding partners, solvating the system, and equilibrating; this would be performed by a script such as gen_istate.sh. In this tutorial--in order to make things simpler--we generate initial states from the basis state by simply symlinking the restart file. Given an initial state, WESTPA next runs a script to calculate the progress coordinate for that initial state. This task is performed by get_pcoord.sh. Note that this script is only called for newly-generated iniitial states; progress coordinates for trajectory segments should be calculated in runseg.sh. Finally, after calculating the progress coordinate for the initial state, runseg.sh is called to propagate the trajectory.

In the case that trajectory segment is continued from a prior trajectory segment, only runseg.sh is called. For both continuations of prior segments and new segments, it is very important to include the system's starting configuration in the progress coodinate calculation. In other words, the first progress coordinate value of a segment should be the same as the final progress coordinate value of the parent segment. For Amber, this means we load the restart file in addition to the coordinate file (.nc) when performing the calculation.

In addition to gen_istate.sh, get_pcoord.sh, and runseg.sh,, a few helper scripts are used:

aux_functions.py

When progress coordinates and other data is sent to WESTPA, we pipe the information to a temporary file, which WESTPA will load. By default, WESTPA calls numpy.loadtxt on these files, but other methods can be defined by the user. In aux_functions.py, one such method is defined; it is used to load the atomic coordinates of Na+ and Cl-.

These atomic coordinates are not needed by WESTPA, and are therefore an example of auxiliary data. WESTPA includes the ability to calculate and store such auxiliary data as the simulation is run, which is often easier than looping over iterations and segments afterwards. The choice of what auxiliary data to store is up to the user and will depend on the goals of the project and the constraints of the system. For example, while storing atomic coordinates in the WESTPA data file west.h5 is convenient, it may not be feasible for large systems.

The location (module and function name) of the function for loading atomic coordinates is specified to WESTPA in west.cfg.

env.sh

This script sets environment variables that may be used across the simulation. These include the root directory for WESTPA, the root directory for the simulation, the name for the simulation, and the python executable to use. It also sets the executables for Amber; using an environment variable for this purpose makes it easier to transition code to different hardware or test different builds or flags of an MD code without editing multiple files.

Make sure that your $AMBERHOME environment variable is correctly set and that $AMBERHOME/lib is in $LD_LIBRARY_PATH. This should be the case if you have sourced amber.sh from your Amber installation directory.

init.sh

This script initializes the WESTPA system. It copies the restart file from the end of the equilibration to bstates, and it copies the topology file to amber_config. It removes files from previous runs and uses gen_istates.sh and get_pcoord.sh to generate initial states. This is also where the basis states and target states are defined. For this system we define the bound target state as 1.6 Å separation; by specifying --tstate bound,1.0" we indicate that the target state (named "bound") comprises the bin that contains progress coordinate value 1.0. Once walkers reach the bin containing this value (i.e. our first bin), they are recycled. The file init.sh is also where we specify the number of initial states to initialize per bin, as --segs-per-state.

post_iter.sh

This script cleans up after each iteration. WESTPA simulations can generate large numbers of files, potentially conflicting with filesystem restrictions. After each iteration, post_iter.sh packages the segment logs as a tar file.

run.sh

This script is used to run WESTPA.

tar_segs.sh

This script is used to tar segments after the WESTPA simulation has been run, in order to reduce the number of files produced. In order to allow extension of the simulation, the last segment is not tarred. Typically, it is advisable not to tar segments after each iteration (i.e. in post_iter.sh); while the main WESTPA process is tarring, other cores are idle, potentially wasting CPU time.

A summary of the scripts is as follows:

Input File Description
env.sh set environment variables
gen_istate.sh generate initial states from basis states
get_pcoord.sh calculate progress coordinate for initial states
aux_functions.py functions for loading auxilliary data
runseg.sh segment implementation
post_iter.sh post-segment cleanup
init.sh initialize WESTPA
run.sh run WESTPA
tar_segs.sh tar segments

The above files are listed roughly in the order in which it is appropriate to configure them. gen_istate.sh, get_pcoord.sh, runseg.sh, and post_iter.sh are located in the westpa_scripts subfolder.

Running the simulation

Make sure that the environment variable $WEST_ROOT is set and points your WESTPA installation. To set this variable, enter export WEST_ROOT=/path/to/your/westpa/installation at the terminal, modifying the command with the appropriate path. You can check if $WEST_ROOT is set by typing echo $WEST_ROOT at the terminal.

From the simulation root directory ($WEST_SIM_ROOT), initialize the simulation using the command:

./init.sh

The main job of this script is to call w_init, which creates the main WESTPA data file west.h5. In addition, it moves around a few files, creates directories, and clears out data from previous runs.

After init.sh has completed, call ./run.sh to start the WESTPA simulation. This will run the simulation using as many cores as are available on the host machine. If you wish to run a WESTPA simulation on a cluster, a script similar to run.sh can be executed from a batch script; see the WESTPA wiki for examples.

Note that this simulation will take some time to run! Currently, the simulation is set to run for 10 iterations, which will take a couple hours. Later in the tutorial, we will extend the simulation to 100 iterations, which could take a few days on typical hardware. If you do not wish to run the whole simulation, a complete WESTPA data file (west.h5) can be downloaded here; this file is sufficient to complete the rest of the tutorial.

Analyzing the data

Output

Output File Remarks
traj_segs output from each iteration and segment
seg_logs log files from each iteration and segment
west.h5 WESTPA output in hdf5 database
west.log WESTPA log file

traj_segs

This folder stores the results of the WESTPA simulation, organized by iteration and segment. This includes all files generated by runseg.sh, including those generated by Amber. For this system, the only files saved are seg.rst, seg.nc, and seg.log corresponding to the restart file (coordinates and velocities of the final timepoint), the coordinates for the full trajectory, and the log file. After the simulation has been run, tar_segs.sh may be used to reduce each iteration to a single tar file.

seg_logs

This folder stores logs from each iteration and segment. post_iter.sh has been used to combine each segment into a single tar file.

west.h5

This file stores the simulation output in an hdf5 database. This includes the relationships between successive walkers, bin weights, progress coordinates, and auxiliary data.

west.log

This file contains a brief log of simulation progress. As WESTPA runs, it outputs information such as the current iteration number, the number of populated bins, and the time needed for each iteration in this log. This is also where errors are output.

Since only 10 iterations have been run, we do not yet have enough data to analyze. Edit west.cfg and change max_total_iterations to 100. Extend using the command:

./run.sh

Computing the association rate

WESTPA includes a variety of tools for analysis, all of which are parallellized and may be run on a single thread, multiple processors, or across multiple nodes on a computing cluster. WESTPA offers two modalities for analysis: (i) a set of command line tools, each designed for a specific task, and (ii) an automated, all-in-one analysis tool supporting multiple analysis schemes and interactive analyses via an "ipython"-like interface. In this tutorial, we will use the latter option, called w_ipa.

Analysis schemes are defined in the WESTPA configuration file west.cfg. At the bottom of west.cfg, you should see a section similar to the following:

  analysis:
     directory: ANALYSIS                
     kinetics:                          
       step_iter: 1
       evolution: cumulative
       extra: [ 'disable-correl' ]
     analysis_schemes:                  
       TEST:
         enabled: True
         bins:
           - type: RectilinearBinMapper
             boundaries: 
               - [0.0,2.6,10.0,'inf']
         states:
           - label: unbound
             coords: 
               - [10.1] 
           - label: bound
             coords: 
               - [0]

The w_ipa tool will analyze the simulation by figuring out which bin each trajectory segment exists in, at each time point. Then the tool will calculate fluxes and estimate rate constants between the states. Options regarding this analysis should be specified in the WESTPA configuration file west.cfg, as exemplified above.

Option Description
directory Specifies the directory where analysis files are stored
analysis_schemes Under this subheading, one or more "analysis schemes" may be specified; each analysis scheme comprises a set of analysis parameters. Using multiple analysis schemes is a convenient way to test out the effect of different parameters on your results.
bins Here we specify binning parameters for analysis of the simulation, which may be the same or different from binning parameters used during the course of the weighted ensemble simulation. Both the type of bin mapper and bin boundaries should be specified.
states Here we specify states, which are defined as collections of one or more bins. During the analysis, rate constants and fluxes will be calculated from each state to every other state. Often times, you may want the states to correpond to a physically meaningful partition of the state space, such as bound and unbound states. Each state should have a label. Boundaries of the state are specified using a list of coordinates; each of these coordinates is assigned to a bin, as defined in the bins section of the yaml file, under analysis_schemes. Each bin to which a coordinate is assigned is then a member of the state being defined. For example, - [0,1] would specify that all bins including progress coordinates 0 or 1 would be included in the corresponding state.
kinetics Options specific to analyzing the kinetics
step_iter When analyzing and averaging kinetic data, skip ahead in blocks of the specified number of WE iterations.
evolution Specify how to average kinetic data. Here, we cumulatively average.

Parameters may be specified globally or for each analysis scheme individually. In the example at hand, the options under the keyword kinetics are applied to every analysis scheme (though here we only have one analysis scheme). Alternatively, the same options could be specified under each analysis scheme.

To run the analysis, simply run w_ipa from the simulation directory. w_ipa will automatically perform analyses that need to be performed or load up the analysis files if they have already been completed. After running the analyses, w_ipa will drop you to an ipython-like interactive python interpreter. From here, you have access to WESTPA datasets and a variety of analysis routines through the w object.

We can plot the cumulatively averaged flux into the bound state using the built-in terminal-based plotter of w_ipa. First, we will check what state index should be passed to the plotting function. To check the state labels and their corresponding indices, type the following command:

In [1]: w.state_labels

Out [1]:
State labels and definitions!
0: bound
1: unbound
2: Unknown

We can see that the bound state corresponds to index 0. Therefore, to plot the flux into the bound state, type:

w.current.direct.target_flux_evolution.plot(0)

The x-axis represents the iteration number, and the y-axis represents the flux into the bound state in units of τ-1.

Alternatively, we can plot the flux into the bound state using pyplot:

import matplotlib.pyplot as pyplot
pyplot.plot(w.direct['target_flux_evolution']['expected'][:,0], color='black')
pyplot.plot(w.direct['target_flux_evolution']['ci_lbound'][:,0], color='gray')
pyplot.plot(w.direct['target_flux_evolution']['ci_ubound'][:,0], color='gray')
pyplot.xlabel('Iteration')
pyplot.ylabel(u'Mean Flux (1/\u03c4)')
pyplot.show()

We can see that the flux has plateaued, indicating that the simulation has reached steady-state conditions. When calculating the rate, we discard the portion of data prior to reaching a steady-state. In the above example, we might say that a steady-state has been achieved by 20 iterations. We will therefore calculate the association using only the last 80 iterations (the final 80% of the simulation). To do so, exit w_ipa and add another line to west.cfg, under the kinetics section of analysis:

     kinetics:                          
       step_iter: 1
       first_iter: 20 ##<<< add this to start the analysis at iteration 20. 
       evolution: cumulative
       extra: [ 'disable-correl' ]

Then run w_ipa again. While we could again examine target_flux_evolution, this time we will look at a different dataset. Type the following from the w_ipa interactive interpreter:

In [1]: w.current.direct.rate_evolution

Out [1]:
rate_evolution data:                                
bound   -> unbound: mean=4.986332950426400e-02 CI=(1.880728819138737e-02, 8.701677126259222e-02) * tau^-1                                                                                                       
unbound -> bound  : mean=6.865117023515006e-04 CI=(3.586815972209720e-04, 1.060736136259851e-03) * tau^-1                                                                                                       
To access data, index via the following names:      
['ci_lbound', 'ci_ubound', 'corr_len', 'error', 'expected', 'iter_start', 'iter_stop', 'name', 'plot', 'raw', 'sterr']

Unlike the target_flux_evolution, the rate_evolution takes into account the amount of probability in the bound and unbound states. Since we have performed a steady state simulation in which all trajectories begin in the unbound state, the unbound → bound rate should be the same as the flux into the bound state. In other words, for this simulation, it doesn't matter if we look at the target_flux_evolution or rate_evolution dataset. Still, in other cases the rate_evolution dataset is useful. In particular, this is the case for equilibrium WE simulations.

The association rate in your simulation should be within an order of magnitude of the results shown above. Since τ for our simulation was 5.0 ps, in order to determine the association rate in units of μs-1, the flux should be multiplied by 2×105, giving an association rate of 137 μs-1, with a 95% confidence interval of 72 μs-1 to 212 μs-1. We can then convert the association rate to units of M-1s-1, taking into account the box's volume (in this example, it was 40716 Å3):

(137 × 106 s-1) × (40716 Å3) × (1 liter/1027 Å3) x (6.02 x 1023 ions/ 1 mol) = 3.4 × 109 M-1s-1

Finally, we can repeat the calculation using a different averaging method. The above examples all employ a cumulative average, in which the rate estimate for a given iteration employs all of the data starting from the beginning of the simulation (or the iteration specified by first-iter), up to and including the given iteration. Now we will try using a block averaging method, in which an estimate is obtained for consecutive windows of a given size. For example, we can obtain an estimate for the period between 0 and 10 iterations, then the period between 11 and 20 iterations, then the period between 21 and 30 iterations, and so on. Open west.cfg and add an additional analysis scheme, as follows:

       BLOCKED:
         enabled: True
         kinetics:                  
           step_iter: 10
           first_iter: 1
           evolution: blocked
           extra: [ 'disable-correl' ]
         bins:
           - type: RectilinearBinMapper
             boundaries:
               - [0.0,2.6,10.0,'inf']
         states:
           - label: bound
             coords:
               - [0]
           - label: unbound
             coords:
               - [10.1]

You should now have two analysis schemes: TEST (the original scheme), and BLOCKED, the new scheme. While kinetics options are set globally, in the section of west.cfg above the TEST scheme, the kinetics options for the BLOCKED scheme supercede the global kinetics options. Since the bins and states are the same for both analysis schemes, they could also be specified in the area for global parameters.

When you open w_ipa again, the tool will run the additional analysis. At the command prompt, type w.scheme to see what scheme you are currently using. If the scheme is not BLOCKED, change the scheme by entering w.scheme = 'BLOCKED' at the iteractive prompt. To see the results of the rate analysis, type w.current.rate_evolution.direct at the interactive prompt. Now switch the scheme to TEST by typing w.scheme = 'TEST', and look at the results of the corresponding rate analysis by typing w.current.rate_evolution.direct at the prompt. You should see that the results are different for the two analysis methods. You can also compare the results of the analysis at different points in the simulation using the iteration property of of the w object; type w.iteration to display the current iteration, and change the iteration by typing (for example) w.iteration = 50 to change the iteration to 50.

Visualizing a selected pathway

The tool w_ipa also includes methods to make concatenating segments for one of your completed pathways straightforward. If it is not already open, run w_ipa. From the prompt, type the following:

w.current.successful_trajectories

This will return a list of pairs of state indices indicating the observed types of transitions (e.g., [(1,0)] indicates that a transition was observed from the unbound state 1 to the bound state 0). As a reminder, you can check the state labels by typing w.state_labels from the w_ipa prompt.

It may be the case that no successful trajectories where observed during the current iteration. If that is the case, try searching for an iteration that includes successful trajectories by setting the iteration to a different value (e.g., w.iteration = 99), and again examining w.current.successful_trajectories. Once you find an iteration with one or more successful trajectories, type the following to return a list of the corresponding segment indices:

s = w.current.successful_trajectories
s[(1,0)]

Next, trace the trajectory's history using the method w.trace. For example, if segment 10 was successful, type t = w.trace(10).

Using atomic coordinates stored in west.h5, we can generate a complete trajectory viewable using Visual Molecular Dynamics. To do so, we may run the script amberTraj.sh included in the main simulation directory. This script will take the successful trajectories from the completed simulation and stitch them together using cppTraj. The resulting file trace.nc will be placed inside of a new directory that the script will create trajAnalysis.

Open the trajectory file in vmd (or another molecular graphics tool) and watch the movie. If you are using vmd, you will need to change the representation (e.g., to van der Waals representation) in order to see the ions. Based on how we have imaged the trajectory, the Na+ will appear to stay fixed in the center of the box. The Cl- ion may appear to jump around as it crosses the periodic boundaries. At the end of the trajectory, you should see that the ions have associated!


Useful links

Useful hints

  • Make sure your paths are set correctly in env.sh
  • If the simulation doesn't stop properly with CTRL+C, use CTRL+Z.
  • Another method to stop the simulation relatively cleanly is to rename runseg.sh; WESTPA will shut the simulation down and prevent the hdf5 file from becoming corrupted. Some extra steps may be necessary to ensure that the analysis scripts can be run successfully.

References

  1. a b I. S. Joung and T. E. Cheatham. "Determination of alkali and Halide Monovalent Ion Parameters for Use in Explicitly Solvated Biomolecular Simulations" J.Phys. Chem. 2008. 112, 9020-9041.
  2. a b M. C. Zwier, J. W. Kaus, and L. T. Chong. Efficient Explicit-Solvent Molecular Dynamics Simulations of Molecular Association Kinetics: Methane/Methane, Na+/Cl-, Methane/Benzene, and K+/18-Crown-6 Ether. JCTC, 2011. 7, 1189-1197
Clone this wiki locally