NAMD LAB 2: BUILDING AND SIMULATING A PROTEIN-BILAYER SYSTEM

by Mitch Gleed

BETA

Objectives:

The first portion of this lab using CHARMM-GUI is something you can do at home, so if you would like more lab time for coding and simulation, we would recommend you complete section 1 of the lab beforehand.

In this lab we will be working with the same protein that we simulated in CHARMM labs 7 and 8, the Influenza A M2 protein. We construct the protein-bilayer system using CHARMM-GUI, a CHARMM web interface developed at Lehigh University. We will then use NAMD to equilibrate, CHARMM to insert a drug of choice, VMD to create comparison sets for NAMD restraints, and then perform production dynamics in NAMD.

This lab is a preparation lab for the Umbrella Sampling lab, which will use trajectories you simulate in this lab and rely heavily on code you create as well. To this end, there are actually a few points of flexibility in the lab, so depending on if you are doing this lab with other students in a class setting, you might each pick slightly different, but related, options so that you can compare results in the end (like a small research study!). Take a moment with your class and decide what each class member will do differently and what variables should remain the same. Here are just a few of many examples of components where you will have flexibility:

1. CHARMM-GUI

Navigate to http://www.charmm-gui.org. Select Input Generator in the navigation bar, then select Membrane Builder > Bilayer Builder in the updated navigation bar. Scroll down to Protein/Membrane System and enter the M2 PDB protein code of your choice. Average structures have been built for all of the following structures–for the sake of ease you may choose between these, otherwise you will need to modify the restraint protocol used in this lab involving averaged PDB structures.

Once you have entered the PDBID of your choice, ensure Download Source is set to “OPM”, and then scroll down to and click Next Step: Select Model/Chain.

Select Model/Chain

Manipulate PDB

All of the boxes here should be unchecked except for the following:

The protonation feature is limited, instead we prefer to use the mutation feature to set His protonation states. You must set His protonation for all occurences in the peptide. For example, the 2KQT structure has one His residue on each M2 segment, for a total of four; 2L0J, on the other hand, has two on each segment, for a total of eight. We have used the HSE epsilon protonation state for His37 in our lab publications, so we suggest you use that here.

For 2KQT, it should look like the following:

SEGID RESID AA Mut
PROA 37 HIS HSE
PROB 37 HIS HSE
PROC 37 HIS HSE
PROD 37 HIS HSE

Amino acid mutations such as S31N would be performed in a similar fashion.

Generate PDB and Orient Molecule

Calculate Cross-Sectional Area

In this portion, you will decide which type of lipids to use and what type and size of system you want.

As far as the membrane lipid composition goes, you will generally want to choose a lipid according to the lipid used in the study connected to the PDB. If you want to use this as a comparison in the umbrella sampling study, you could use homogenous membranes with varying lipids or varying mixtures of lipids, cholesterol, etc. by changing the lipid ratios in the upper and lower membrane leaflets.

Determine the System Size

If you have time and want to publish using this CHARMM-GUI structure, consider following the extra instructions under Pore Water Options. Otherwise move on to the next step.

This step typically takes the longest for CHARMM-GUI to process, and you may even get a blank white screen for a while. You might consider copying the link that shows just below the yellow CHARMM box in case of a time out.

Build Components

Click to continue

Assemble Components (1 of 2)

Click to continue

Assemble Components (2 of 2)

You can now view the ultimate system size. If you are hoping to publish using this structure, you might want to document some of these statistics.

Generate Equilibration and Dynamics Inputs

Click the red box labeled download .tgz, located near the top right corner of the Membrane Builder interface. This will download a file named charmm-gui.tgz to your system. Upload this file to your personal NAMD lab 2 directory on the supercomputer.

You may close CHARMM-GUI now.

2. Running generated equilibration input files in NAMD

CHARMM-GUI provides output in multiple MD suite formats to automatically begin equilibrating the systems, and for this lab we will use the generated NAMD input files. The simulation parameters and restraints vary slightly between MD packages, but the simulations follow essentially the same pattern.

The equilibration protocol involves progressive relaxation of restraints. Each step lasts 50 - 100 ps. The relaxation is performed as follows (per the CHARMM-GUI site):

  Step 1 Step 2 Step 3 Step 4 Step 5 Step 6
BB 10.0 5.0 2.5 1.0 0.5 0.1
SC 5.0 2.5 1.0 0.5 0.1 0.0
tforce 2.5 1.0 0.5 0.0 0.0 0.0
mforce 2.5 0.0 0.0 0.0 0.0 0.0
ion 10.0 0.0 0.0 0.0 0.0 0.0

Other restraints implemented include a barrier restraint to keep water molecules away from the hydrophobic lipid regions, dihedral restraints on lipids, etc.

Extracting CHARMM-GUI output

Now that charmm-gui.tgz is in your NAMD lab 2 directory, execute the following command: tar -xvf charmm-gui.tgz. This will extract the contents of the compressed file (analagous to unzipping a .zip file in Windows).

This will create a directory called “charmm-gui” that contains all of the files used by CHARMM-GUI and the automatically generated input files for NAMD.

Submitting NAMD GPU jobs

NAMD performs very well on standard Linux clusters, but is particularly fast when used on nodes with GPU’s. A sample submission script titled bash.pt1_equilibrate.sh is provided for you, let’s go ahead and start digging into it. Open the file and examine its contents.

You will see some familiar lines near the start that tell the scheduler what resources to ask for:

#!/bin/bash

#SBATCH --time=18:00:00   # walltime
#SBATCH --ntasks=24   # number of processor cores (i.e. tasks)
#SBATCH --nodes=1   # number of nodes
#SBATCH --gres=gpu:4
#SBATCH --mem=64G   # memory per CPU core

This script is tailored to the architecture of the m8g cluster, which is comprised of nodes each with 24 processing cores, 4 GPU’s, and 64 GB of RAM. The option --gres=gpu:4 requests 4 GPU’s.

We will use this script to modify the restraint protocol provided by CHARMM-GUI. To save time for the meat of the lab, we won’t delve into the syntax of the section labeled # Add a channel backbone tilt restraint..., but it is appending an extra restraint to the NAMD equilibration protocol given by CHARMM-GUI to keep the channel axis lightly restrained to the Z axis using the tilt collective variable.

To allow NAMD to use the GPU’s, we need to give Linux an environment ready for GPU simulations.

Append the following code beneath the line # NAMD variables:

# NAMD variables 
module purge
module load compiler_gnu/6.4
module load cuda/8.0
export dir=/fslhome/mgleed/software/namd/exec/NAMD_Git-2017-11-04_Linux-x86_64-multicore-CUDA

This will reset the environment with module purge, and give the libraries needed for GPU simulations with compiler_gnu/6.4 and cuda/8.0. Finally, the directory containing the NAMD GPU-ready executable is given by the variable dir.

Now we need to launch NAMD and load the input files. Append the following code beneath the line # Run NAMD:

# Run NAMD 
cd $SLURM_SUBMIT_DIR/charmm-gui/namd
$dir/namd2 +p${SLURM_CPUS_ON_NODE} +idlepoll +devices "0,1,2,3" step6.1_equilibration.inp > step6.1_equilibration.inp.log

The cd command moves to the folder containing the NAMD input files. Using cd in a submission script will help prevent creating output files in a location you don’t want. The next line is where NAMD is invoked.

There are simpler ways to run GPU jobs, but this is the fastest we’ve discovered for the m8g cluster based on a variety of benchmark results.

This submission script is nearly complete, but while we have the GPU node scheduled, let’s go ahead and run all of the steps of equilibration, not just 6.1. Create new lines for step 6.2 through 6.6, based on the last line, and append them to the end of the script.

If you don’t want to use GPU’s for your simulation (such as when the GPU nodes have 100% utilization, check the site), remove the line containing --gres, reduce the number of processors, RAM, etc., and use the following pattern to launch jobs:

# NAMD variables 
module purge 
module load namd/2.12_openmpi-1.8.5_gnu-5.2.0

# Run NAMD 
mpirun $(which namd2) step6.1_equilibration.inp > step6.1_equilibration.inp.log

To learn more, see also MaryLou compute resources to get an idea of what resources you can request and the Job Script Generator to get the scheduler commands you need.

Now, submit this script to the scheduler using the sbatch command. Watch the progress of the job using watch squeue -u [username] and type Control+C to exit watching.

If you have followed the default lab suggestions, this portion will take roughly 45 minutes to complete after the GPU node is scheduled. However, you can continue through the lab while waiting for the job to finish!

3. Adding a drug to the protein-bilayer system using CHARMM

We will now add a drug to the protein-bilayer system using CHARMM. There are several important steps that must be accomplished to be successful:

This is a surprisingly daunting task to perform properly. Rather than spend our time and energy making you create this file on our own, we will provide the file for you. In this way we will continue our NAMD focus in this lab.

Insert drug

Let’s hit the highlights of our the next script, charmm.adddrug.str.

! load topology
	stream @topparloc/load_toppar.str
	read rtf card append name alm/45drugs.rtf
	read para card flex append name alm/45drugs.prm

This code block demonstrates how to append topology and parameters to already loaded topology and parameters. Using append after read rtf causes the alm/45drugs.rtf topology to appended; flex append performs the same function for the parameter file. The 45drugs files contain topology and parameters for 45 adamantyl-derived drugs studied in the Busath lab, some of which vary only by protonation state or chirality.

! loads coordinates from charmm-gui NAMD restart file
	read namd file "charmm-gui/namd/step6.6_equilibration.coor"

This code block demonstrates how to load NAMD binary coordinate files into CHARMM. The NAMD import feature is very rudimentary, and formatting is strict. The file name must be encircled by double quotes (") and the coordinates can only be read into the main coordinate set (you can get around this limitation by using coor copy comp and then coor swap).

! reorient entire system to "fit" the simulated protein backbone to the average protein backbone, will be reversed later 
	coor orie rms sele backbone end

This code block causes the main coordinate set, the protein-bilayer system, to be RMS oriented to the backbone atoms of the comparison set, which is the average PDB structure.

The section that starts with ! find alpha carbons of residue S31 and place Amt cage here with N facing inferiorly to roughly match configuration 1 from the 2015 paper is where the drug placement happens. The drug is first moved in the X/Y plane to be centered on the channel, then the drug is placed near the Ser31 residue using a combination of coor stat and calc.

The next block demonstrates how to delete atoms near atoms of interest.

! delete any water overlapping the drug 	
	delete atom sele .byres. ((resn TIP3 .and. type OH2) .and. (fulldrug .around. 2.2)) end

This particular command will delete entire residues (.byres.) containing any water molecule (resn TIP3) oxygen (type OH2) if the oxygen lies within 2.2 angstroms (.around. 2.2) of the drug (fulldrug definined earlier in script).

Submit CHARMM script

Notice the variables pdbid, drugfilename, drugsegid, and topparloc aren’t defined anywhere in the script, but are used multiple times. These are variables that are passed to CHARMM by the submission script. Recall from previous labs the use of environment variables set in Linux and passed to CHARMM or other programs.

Open the file bash.pt2_adddrug.sh and fill in the values for the environment variables that are empty. If you don’t know the proper value for drugsegid, use a text editor or vi to open the PDB file associated with the drug you want to simulate and check the far right column. For example, for protonated Amt, you would use alm034 for drugname and L034 for drugsegid (and alm035 and L035 if deprotonated).

See the alm/ directory for the 45 drugs you can choose from. Amantadine is alm034 (+) and alm035 (neutral) and Rimantadine is alm149 (+) and alm150 (neutral).

Once you are done, check to ensure the CHARMM-GUI job from earlier is complete, and then execute ./bash.pt2_adddrug.sh to run the CHARMM job on the internode (rather than scheduling using sbatch). If the job from earlier isn’t done yet, you can move on to the next section and come back to execute the script.

4. Preparing reference files for colvars using VMD

VMD, like CHARMM, can be used from the command line without any graphic user interface. Scripting is based on the programming language Tcl.

A great way to learn to use VMD scripting is by using VMD on your personal computer. At the start of a session, go to File > Log Tcl Commands to File... and then load molecules, play with representations, atom selections, etc. When your session is complete, load the file you saved at the start of the session in a text editor, and it will give you the Tcl commands used to run your session behind the scenes.

One advantage of VMD is that it supports many plugins and is useful for system building, trajectory analysis, and more. In this lab, we use VMD to construct reference coordinate files for use by the collective variables module in NAMD, as VMD has the ability to easily manipulate alpha and beta columns (to tag atoms of interest) in PDB files.

Create a new file named vmd.compsets.inp. The general idea for this script will be as follows:

Load input files into VMD

Let’s start with the following code:

# load in structure produced by charmm.adddrug.str
	mol new output/2l0j_alm034_oriented.pdb type {pdb} first 0 last -1 step 1 waitfor all

To load a file into VMD, you load it as a new molecule, just as you would in the GUI version you are used to. You start with the keyword mol to utilize the molecule facility, say new to specify you are creating a new molecule, insert a filename to load as the new molecule, use type {pdb} to tell the input reader to expect PDB formatting. The rest are used primarily for trajectories: first 0 means to start with frame 0, last -1 would mean to finish on frame -1, but -1 really means there is no last frame as in this case, step 1 loads every 1 frame of the trajectory, and waitfor all causes VMD to completely load the file before moving on to the next command.

Now that we have loaded our file into VMD, let’s load in average protein structure from the PDB models, which we will use in backbone restraints and to align the simulating protein to the cartesian Z axis as needed. Following the pattern from the previous code block, make a new line that loads the average protein structure as a new molecule. You can find the average structures for several PDB structures in the avg/ directory.

Using set and atomselect to replace protein coordinates

Both molecules are now in VMD, so we can start replacing protein coordinates, etc. Consider the next code block, which you will add to your script:

# replace original backbone coordinates with avg. backbone coordinates
	set repprotein [atomselect 0 "protein"]
	set newprotein [atomselect 1 "protein"]
	$repprotein set {x y z} [$newprotein get {x y z}]

set is used to define a variable. For example, set x 3 would set the variable x to the integer 3. In this case, we use set to define variables as atom selections.

The repprotein variable contains a selection of all the protein atoms in molecule 0, the first structure we loaded. newprotein contains a selection of all the protein atoms in molecule 1, the second structure we loaded.

$repprotein set {x y z} [$newprotein get {x y z}] therefore is interpreted by VMD as atomselect 0 "protein" set {x y z} [atomselect 1 "protein" get {x y z}]. This causes the protein atoms in the first structure to have their coordinates replaced by the coordinates of the atoms in the second structure (the average PDB structure).

Create comparison sets and write PDB’s

Now we can start creating comparison files for restraint references. Append the following code:

# create comparison file sets
	set all [atomselect 0 "all"]
	
	# drug cage
	$all set beta 0
	set slxn [atomselect 0 "segname $::env(drugsegid) and (name $::env(cc1) $::env(cc2) $::env(cc3) $::env(cc4) $::env(cc5) $::env(cc6) $::env(cc7) $::env(cc8) $::env(cc9) $::env(cc10))"]
	$slxn set beta 1
	$all writepdb output/colvar_drugcage.pdb

To use environment variables in VMD, use export variableName = variableValue in your submission script and use $::env(variableName) to use the passed variable in VMD.

The above code will create a restraint file in which the drug cage carbons are marked for restraint. Following this pattern, continue the script and create another restraint reference file to restrain the protein backbone atoms. Set the beta property of all the atoms back to zero, use the atom selector “backbone” to set the beta property of all the backbone atoms to 1, and name your output file output/colvar_bb.pdb.

Run VMD

Create a new file named bash.pt3_compsets.sh. Add the following code at the start:

#!/bin/bash
module purge
module load vmd/1.9.1

This will provide us with the environment required by VMD, and provide a link to the executable binary file, which can be called with $(which vmd).

To launch VMD, add the following code to the submission script:

$(which vmd) < vmd.compsets.inp > output/vmd.compsets.inp.log

Now, before we submit, we need to export the adamantane cage carbon variables for use by the VMD script. The export commands were created by the charmm.adddrug.str script already, we just need to stream the file. Add the following code between the module load ... line and the $(which vmd) ... line:

sed -i "s/ EXPORT CC/export cc/" output/alm034_cc.sh # gotta convert the file to lowercase, charmm output is uppercase by default
source output/alm034_cc.sh # file created by CHARMM containing export commands

Finally, there are a few places in the VMD script where you should be using environment variables, as the code provided here assumes you are using 2L0J and amantadine. Replace these hard-coded values and substitute environment variables–be sure you add the corresponding export lines to the bash.pt3_compsets.sh file as well, as you did in your CHARMM submission script from earlier.

Now, ensure the NAMD simulations and drug insertion in CHARMM from the previous steps worked successfully, and then run the VMD script on the internode with ./bash.pt3_compsets.sh. If the other scripts haven’t finished yet, continue to the next section and revisit these later.

5. Perform dynamics simulations on the protein-bilayer system in NAMD

Now we will perform 1 nanosecond of simulations on the system we have created, utilizing colvars to restrict the protein backbone and track the drug position.

Copy the file step7.1_production.inp found in charmm-gui/namd/ to the lab directory containing the other scripts we’ve created in this lab and name it namd.production.inp. This is the CHARMM-GUI autogenerated script for doing production dynamics on a system after equilibration, and we will use it as a template for our new protein-drug system.

This script needs some important adjustments:

Let’s make these changes now.

Filenames

Replace the values of structure and coordinates at the top of the file to the paths of the protein-drug system PSF and PDB files, respectively. (These are found in output/ and were created in CHARMM earlier in the lab.)

We recommend that you utilize variables, again. For instance, if we are doing Amantadine in 2L0J, you could put set basename 2l0j_alm034 at the top of the file and then substitute portions of the filenames in the structure/coordinates lines with ${basename}. The remainder of the lab will assume Amantadine in 2L0J, but go ahead and replace these patterns with a variable to make your life easier down the line (and in the next lab!)

Next modify outputName, which will give the pattern of all of NAMD’s output.

- outputName         step7.1_production;
+ outputName         output/2l0j_alm034_production;

In this case, the trajectory file will be called 2l0j_alm034_production.dcd in the output/ directory.

Changing from a “restart” script

We are not using atom velocities from the equilibration phase as the PSF has changed. We can, however, use the cell dimensions to get an accurate system volume. Make the following changes to eliminate unecessary loading of restart files:

- set inputname      step6.6_equilibration;
- binCoordinates     $inputname.coor;    # coordinates from last run (binary)
- binVelocities      $inputname.vel;     # velocities from last run (binary)
- extendedSystem     $inputname.xsc;     # cell dimensions from last run (binary)
+ extendedSystem     charmm-gui/namd/step6.6_equilibration.xsc;     # cell dimensions from last run (binary)
+ firsttimestep      0;

dcdfreq represents the number of steps of simulation to perform before writing a frame to the trajectory file. Since the colvars module will be collecting a lot of our important data in real-time during the simulation, we won’t need as much information from the trajectory file. We can safely reduce the trajectory writing frequency by changing the step count:

- dcdfreq           1000;
+ dcdfreq           5000;

Adding drug parameters

NAMD parameter handling can be kind of finicky. In our experience, loading the entire set of parameters, even parameter sets we don’t need, prevents errors without reducing performance. For example, the parameter file that contains water and ions has bond information referencing some of the other redundant parameter files.

We don’t need to read the step5_assembly.namd.str file here, so remove it and then add parameters for amantadine, which is contained in the parameter file for the set of 45 drugs provided in the lab:

- source step5_assembly.namd.str
+ parameters alm/45drugs.prm

Image wrapping and Particle-mesh Ewald

Simplify the following statement, as we are using a rectangular box for the lab:

- if { $boxtype == "hexa" } {
- wrapNearest         on;                # use for non-rectangular cells (wrap to the nearest image)
- } else {
- wrapNearest         off;               # use for non-rectangular cells (wrap to the nearest image)
- }
+ wrapNearest         off;

Make the following changes to let NAMD automatically calculate PME grid dimensions:

- source checkfft.str

- PMEGridSizeX     $fftx;                # should be close to the cell size 
- PMEGridSizeY     $ffty;                # corresponds to the charmm input fftx/y/z
- PMEGridSizeZ     $fftz;
+ PMEGridSpacing      1.0;               # forces NAMD determine the PME dimensions automatically

Dynamics settings and Minimization

Earlier in the lab we equilibrated the system using CHARMM-GUI’s generated scripts. Equilibration was performed using velocity scaling to approximate the canonical ensemble (NVT) in which temperature and volume are held constant.

To use velocity scaling and reassignment to produce an NVT ensemble (for equilibration before dynamics), you use the following lines, as seen in the template:

# Constant Temperature Control
reassignFreq        500;
reassignTemp        $temp;

Since we’ve already equilibrated the system, we will not be using velocity reassignment for temperature control, so delete the above lines from your file.

Now that the system is equilibrated (aside from the drug we added), we will perform dynamics simulations by combining temperature control with a modified Nosé-Hoover Langevin piston to control pressure, resulting in an isobaric-isothermal ensemble (NPT). We will follow the CHARMM-GUI protocol found near the end of the step7.1_production.inp file in charmm-gui/namd/.

The NPT ensemble is used here to prevent large fluctuations in energies due to changes in pressure. One could reason the NVT would be acceptable for the dynamics portion of the simulation as well. Using the NVT ensemble would be a reasonable approach for many systems, so ensure you have explored NVT vs. NPT for anything you wish to publish. NPT ensembles are generally more complex than NVT ensembles and often require more computing power per simulation step.

Remove the remaining lines in the file starting at # Pressure and volume control and insert the following in its place:

# Constant Pressure Control (variable volume)
useGroupPressure       yes;            # use a hydrogen-group based pseudo-molecular viral to calcualte pressure and
                                       # has less fluctuation, is needed for rigid bonds (rigidBonds/SHAKE)
useFlexibleCell        yes;            # yes for anisotropic system like membrane 
useConstantRatio       yes;            # keeps the ratio of the unit cell in the x-y plane constant A=B

langevinPiston          on;            # Nose-Hoover Langevin piston pressure control
langevinPistonTarget  1.01325;         # target pressure in bar 1atm = 1.01325bar 
langevinPistonPeriod  50.0;            # oscillation period in fs. correspond to pgamma T=50fs=0.05ps 
                                       # f=1/T=20.0(pgamma)
langevinPistonDecay   25.0;            # oscillation decay time. smaller value correspons to larger random
                                       # forces and increased coupling to the Langevin temp bath.
                                       # Equall or smaller than piston period
langevinPistonTemp   $temp;            # coupled to heat bath

# Constant Temperature Control
langevin                on;            # langevin dynamics
langevinDamping        1.0;            # damping coefficient of 1/ps (keep low)
langevinTemp         $temp;            # random noise at this level
langevinHydrogen       off;            # don't couple bath to hydrogens

View the documentation on temperature control and pressure control if you would like to learn more.

Minimize and run

Add the following lines to the end of the script:

# run
colvars 		on;
colvarsConfig	colvars.production.inp;
minimize		5000;
run 			500000;

This will turn on the collective variables module, use the configurations set in the colvars.production.inp file, minimize for 5000 steps, and perform 500,000 steps of simulation (that is 1 nanosecond, given 2 fs/step).

Create colvars configuration file

Before we submit the script, we actually need a colvars.production.inp file for the NAMD to read. Here we will set up three collective variables:

Create a file in the lab directory called colvars.production.inp and begin with the following code:

colvarsTrajFrequency    500

This command controls how frequently the collective variable information is reported to the colvars trajectory file (time series) in simulation timesteps. For example, in this case, every 500 steps the protein tilt in cosine units, RMSD in angstroms, and adamantane cage position in angstroms will be written to the text file.

Let’s create our first collective variable using the tilt colvar component. Copy the following code to your file:

colvar { 
  name bb_tilt
  outputvalue off
  tilt { 
	atoms {
		atomsFile          output/colvar_bb.pdb
		atomsCol           B 
		atomsColValue      1.0 
	}
	refPositionsFile      output/colvar_bb.pdb
	refPositionsCol       B
	refPositionsColValue  1.0    
	axis (0.0, 0.0, 1.0) 
  } 
}

We have named the collective variable bb_tilt, turned off output data using outputvalue off (because the protein tilt with respect to the cartesian Z axis will be restrained and is not of particular interest), provide a file containing the atoms of interest atomsFile, provide another file containing the positions of the atoms of interest refPositionsFile, and set the tilt with respect to a vector equal to the cartesian Z axis axis (0.0, 0.0, 1.0).

Next we will create the protein backbone RMSD collective variable. COpy the following code to your file:

colvar {
   name bb_rmsd
   rmsd {
       atoms {
          atomsFile          output/colvar_bb.pdb
          atomsCol           B 
          atomsColValue      1.00
	      centerReference
	      rotateReference		  
          fittingGroup {
             atomsFile      output/colvar_bb.pdb
             atomsCol       B
             atomsColValue  1.00
          }
		  refPositionsFile  output/colvar_bb.pdb   	  
       }
       refPositionsFile      output/colvar_bb.pdb
       refPositionsCol       B
       refPositionsColValue  1.00
   }
}

This collective variable is setup in a similar fashion to the last, with a couple major exceptions. At every step, the collective variable is computed with respect to the fittingGroup atoms. Since we are using the same files to fit, and using centerReference and rotateReference, the RMSD collective variable will track internal distortion of the protein backbone atoms. The reference structure is best-fit aligned at every step, so the RMSD of the simulating structure will always be computed with respect to the average structure “centered” and “rotated” onto itself. Were this relative reference not performed, the RMSD would measure the deviation of the protein in translation and rotation, so drifting even slightly in the membrane would result in very high RMSD values.

Let’s now create the collective variable defining the drug adamantane cage position with respect to the channel axis. This is done using the distanceZ colvars component, which computes the distance of a group of atoms from another group of atoms along a single axis. Copy the following code to your file:

colvar { 
  name drugpos 
  distanceZ { 
    ref { 
      dummyAtom (0, 0, 0)
    } 
    main { 
      atomsFile      output/colvar_drugcage.pdb
      atomsCol       B
      atomsColValue  1.00
	  centerReference
	  rotateReference
      fittingGroup {
         atomsFile      output/colvar_bb.pdb
         atomsCol       B
         atomsColValue  1.00
      }
	  refPositionsFile  output/colvar_bb.pdb   
	} 
	axis (0.0, 0.0, 1.0) 
  } 
}

Since we only care about the drug adamantane cage position with respect to the channel axis, we create a dummyAtom group at the cartesian origin to give us a context for positive and negative values. In this collective variable, we use the file colvar_drugcage.pdb to define the atoms of interest. axis (0.0, 0.0, 1.0) is used to measure only the Z component of the distance between the two groups.

Also, like the last collective variable we discussed, rmsd, we will implement a fittingGroup to keep the drug position data with respect to the channel axis, rather than the cartesian axis. This allows for the distance to be always computed relative to the colvar_bb.pdb atoms–if the channel protein were to translate along the x-y plane, tilt away from the cartesian Z axis, or translate along the Z-axis; the distance from the drug carbon atoms to the cartesian origin would be “centered” and “rotated” to be with respect to the average protein backbone atoms created in VMD.

This concept of a moving frame of reference is more likely to be understood with visual context than in a brief explanation in a lab, so feel free to ask your TA’s or professor if you would like to learn more. You can also learn more in the colvars documentation.

An alternative approach would be to make the ref group Histidine 37 backbone atoms instead of dummyAtom. If the main group were above the His backbone atoms it would result in a positive number, and if it were below the His atoms, a negative number would be recorded. Since you can’t do normal atom selections in colvars, you would need to create a PDB file of the system in which the beta column was flagged exclusively for all His 37 backbone atoms. Careful consideration may need to be given with this method, though, as fluctuations in channel length may distort results.

Finally we implement a harmonic constraint on the ``bb_tilt` collective variable. Copy the following code to your file:

harmonic {
	colvars bb_tilt
	centers 0
	forceConstant 1
}

This will create a harmonic restraint with a force constant of 1 on the collective variable defined by bb_tilt. In this lab we lightly restrain the protein’s tilt to align with the cartesian Z axis.

It could be argued this restraint and collective variable are unnecessary for the simulation, particularly when the other collective variables of interest are all with respect to a best-fit average protein structure, so deviations in protein tilt shouldn’t matter. It was implemented here to prevent excessive tilting while this lab was being designed. Abandoning this restraint was considered, but we concluded it would be a good introduction to the syntax of the important tilt collective variables component.

Save the file and return to the lab directory.

Submit to scheduler

Create a submission script titled bash.pt4_simulation.sh in the lab directory.

Insert the following code into the file and adjust any of the #SBATCH parameters as you deem necessary (currently set to use a full FSL m8 GPU node, 3/2018):

#!/bin/bash

#SBATCH --time=24:00:00   # walltime
#SBATCH --ntasks=24   # number of processor cores (i.e. tasks)
#SBATCH --nodes=1   # number of nodes
#SBATCH --gres=gpu:4
#SBATCH --mem=64G   # total memory request

# NAMD variables
export buildtype=multicorecuda
module purge
module load compiler_gnu/6.4
module load cuda/8.0
export dir=/fslhome/mgleed/software/namd/exec/NAMD_Git-2017-11-04_Linux-x86_64-multicore-CUDA

# Run NAMD 
$dir/namd2 +p${SLURM_CPUS_ON_NODE} +idlepoll +devices $CUDA_VISIBLE_DEVICES namd.production.inp > output/namd.production.inp.log

If you have variables you wish to export for use in NAMD, enter them before the # Run NAMD section. For example, you could put export pdbid=2l0j or export drugname=alm034 if you want to use variable substitution in your NAMD script for amantadine in 2L0J. (As a reminder, use $::env(VARIABLENAME) in your NAMD script to reference an exported variable.) Such an approach will be beneficial to your workflow in later labs.

Run your production script with sbatch. (Following this lab with the m8 GPU architecture, completed 500,000 steps in roughly 40 minutes 3/2018)

Analysis

Once your simulation is complete, open the file named 2l0j_alm034_production.colvars.traj in the output/ directory. Using the data from this file, plot the backbone RMSD and the drug position over time. You might have a few rows where the header labels are duplicated throughout the output file; if you do, just manually delete these lines in your spreadsheet or data analysis software.

Here is a sample graph of Amt in 2L0J following the lab protocol, but with 3 million simulation steps. alt text

Assignment: Summarize the meaning of your results in a couple of sentences. If your classmates varied the lab in any way (different channel, drug, etc.), how do your results compare? Submit your plot(s) to your TA.

NAMD Lab 3

Return to home page