REMD with AMBER - Revision history

import>Csw34: /* Analysing the results */

2008-10-24T12:48:19Z

Analysing the results

@@ Line 96: / Line 96: @@
 </pre>
-Note that if you have large files, this may take a while - but it is less boring than opening each in vi and adding a line at the top!
+Note that if you have large files, this may take a while - but it is less boring than opening each in vi and adding a line at the top! Once you have done that, you can use '''ptraj''' to your hearts content. For example, here is a simple script that takes the trajectory files generated during the REX run, and produces input for a WHAM analysis - where the order parameters are the C-alpha RMSD from the crystal structure (1LE1.inpcrd) and the radius of gyration. Note that the trajectory is read in from exchange 11+1(offset)=12 to allow each replica to have equilibrated.
 ===Clustering structures===

import>Csw34: /* Clustering structures */

2008-10-24T12:40:59Z

Clustering structures

← Older revision		Revision as of 12:40, 24 October 2008
Line 100:		Line 100:
	===Clustering structures===		===Clustering structures===

	If you also have the MMTSB Toolset installed - you can cluster the output of your REX runs quite easily. Clustering provides information on the most populated conformational groups, and makes it easy to pick out common structures from the run without manually trawling through every structure one at a time.		If you also have the MMTSB Toolset installed - you can cluster the output of your REX runs quite easily. Clustering provides information on the most populated conformational groups, and makes it easy to pick out common structures from the run without manually trawling through every structure one at a time. Here is an example script that you can use to perform such a clustering separately for each replica in a run. You will need to modify it before you use it due to differing file names and windows over which you are sampling as the MMTSB algorithm can only cope with <50000 structures.

			<pre>
			#!/bin/bash
			# Script to extract PDBs from an AMBER9 REX run and perform
			# a simple clustering analysis for each replica
			# You need to edit the numbers and file names
			# when using more replicas or different systems.
			# The clustering uses the MMTSB script kclust and as a result,
			# MMTSB must be installed and working.

			# NOTE: this script is designed to work in a directory one level
			# below that containing the mdcrd files from your REX run e.g.
			# if the REX was done in /home/csw34/REX then this script should
			# reside in /home/csw34/REX/clustering for example

			# WARNING! The output files from AMBER REX are NOT ok to use directly
			# They lack the blank line at the top which you see in normal AMBER
			# trajectories and so cannot be safely read into ptraj, vmd or well
			# anything. I suggest you add a line using for example this perl:
			# perl -pi -e 'print "\n" if $. == 1' mdcrd.rep1

			# Echo the awk script used later to evaluate the centroids. Credit
			# goes to Ross Walker for this one.
			echo 'BEGIN{b0=2;}' > analyse_centroid.awk
			echo '' >> analyse_centroid.awk
			echo '{centind=index($1,"#Centroid");' >> analyse_centroid.awk
			echo ' c=$2;' >> analyse_centroid.awk
			echo ' getline;centind=index($0,"#Centroid");' >> analyse_centroid.awk
			echo ' FIL0 = sprintf("centroid%2.2d.member.dat",c)' >> analyse_centroid.awk
			echo ' while(centind != 1){' >> analyse_centroid.awk
			echo ' print $1,$3 > FIL0 ;' >> analyse_centroid.awk
			echo ' getline;centind=index($0,"#Centroid");' >> analyse_centroid.awk
			echo ' }' >> analyse_centroid.awk
			echo '' >> analyse_centroid.awk
			echo '' >> analyse_centroid.awk
			echo ' numcent=0; print $2,$1,NR-b0; ' >> analyse_centroid.awk
			echo ' c=$2;' >> analyse_centroid.awk
			echo ' getline; endrec = index("End",$0);' >> analyse_centroid.awk
			echo ' while( endrec != 1 ){' >> analyse_centroid.awk
			echo ' FIL = sprintf("centroid%2.2d.pdb",c);' >> analyse_centroid.awk
			echo ' print > FIL;' >> analyse_centroid.awk
			echo ' getline; endrec = index($0,"#End");' >> analyse_centroid.awk
			echo ' }' >> analyse_centroid.awk
			echo ' b0=NR+2;' >> analyse_centroid.awk
			echo ' } ' >> analyse_centroid.awk

			# Loop over each replica
			# EDIT TO CHANGE REPLICA NUMBERS
			for i in 1 2 3 4 5 6 7 8
			do
			mkdir rep$i
			mkdir rep$i/pdbs
			# Echo the ptraj input file for the replica. Notice that I'm
			# reading in the mdcrd file from frame 24540. This is to allow
			# the replica to equilibrate before we start looking at structures
			# and also - kclust can only deal with up to 49999 structures :(
			# EDIT TO CHANGE FILE NAMES OR STRUCTURE SAMPLING WINDOW
			echo 'trajin ../mdcrd.rep'$i' 24540 74530' > extractrep$i.trajin
			echo 'trajout rep'$i'/pdbs/rep'$i'.pdb pdb'>> extractrep$i.trajin
			echo 'reference ../1LE1.inpcrd' >> extractrep$i.trajin
			echo 'rms reference mass :1-12@CA' >> extractrep$i.trajin
			# Run ptraj to produce the output PDBs
			echo 'Extracting PDBs for replica '$i
			# EDIT TO CHANGE TOPOLOGY FILE NAME OR LOCATION
			ptraj ../1LE1_ff99SB.prmtop < extractrep$i.trajin
			# Produce a file containing all the PDB names for clustering
			ls -1 rep$i/pdbs > rep$i/clustfils
			# Perform the clustering analysis using the MMTSB toolset kclust script
			# CHECK THE OPTIONS OF KCLUST SO YOU KNOW WHAT IT IS DOING
			cd rep$i/pdbs
			echo 'Clustering structures using kclust...'
			kclust -mode rmsd -centroid -cdist -heavy -lsqfit -radius 6 -maxerr 1 -iterate ../clustfils > ../centroid_rep$i
			# Extract info on each cluster centroid using awk script
			cd ..
			awk -f ../analyse_centroid.awk centroid_rep$i \| tee centroid_rep$i.stats
			# How many clusters are there?
			max=`wc centroid_rep$i.stats \| awk '{print $1}'`
			# Now we want to extract the structure closest to the centroid of each
			echo 'Extracting PDB closest to each centroid'
			for ((n=1; n <= $max ; n++))
			do
			# Loop to get around a file naming convention thing
			if [[ "$n" -lt "10" ]]
			then closest=`sort -n -k2 centroid0$n.member.dat \| head -1 \| awk '{print $1}'`
			else
			closest=`sort -n -k2 centroid$n.member.dat \| head -1 \| awk '{print $1}'`
			fi
			# As we didn't start at the first structure, this must be incremented
			# by the appropriate offset. In this case 24539 (24540-1)
			# EDIT IN ACCORDANCE WITH THE START OF THE SAMPLING WINDOW
			closest=$(($closest + 24539))
			# Now, extract the PDB
			cp pdbs/rep$i.pdb.$closest ./bestmember_centroid$n.pdb
			done
			# Remove the thousands of PDBs that are taking up disk space
			rm -r pdbs
			# Return to the correct directory to process the next replica
			cd ..
			echo 'Done with replica '$i
			echo ''
			done
			#Cleaning up
			rm .trajin .awk
			</pre>

import>Csw34 at 12:29, 24 October 2008

2008-10-24T12:29:33Z

← Older revision		Revision as of 12:29, 24 October 2008
Line 1:		Line 1:
			===Setting up AMBER REX===

	To run REMD with AMBER you need '''sander.MPI'''. You can either copy it		To run REMD with AMBER you need '''sander.MPI'''. You can either copy it
	(clust:~bs360/bin/sander.MPI) or create your own. To do so you have to copy the directory		(clust:~bs360/bin/sander.MPI) or create your own. To do so you have to copy the directory

import>Csw34 at 12:29, 24 October 2008

2008-10-24T12:29:00Z

← Older revision		Revision as of 12:29, 24 October 2008
Line 78:		Line 78:
	You can find an example for REMD with AMBER in		You can find an example for REMD with AMBER in
	clust:/sharedscratch/bs360/amberrex_test/		clust:/sharedscratch/bs360/amberrex_test/

			===Analysing the results===

			[[AMBER]] provides a multitude of tools to analyse MD trajectories, all of which are well documented in the AMBER manual under '''ptraj'''. However, before you can use ptraj to analyse any of the output mdcrd.rep* files, you need to make a slight modification. [[AMBER]] differentiates between replica exchange and CTMD output by NOT putting in a blank line at the start of REX trajectory files. This means that if you try to load them into '''ptraj''' for analysis, you will in fact get complete garbage out. Thankfully, this is easily solved by the addition of a blank line using, for example, the following script:

			<pre>
			#!/bin/bash
			# Script to add a blank line at the top of each file. This is important as when doing REX
			# AMBER uses a different output format (without the blank line) for trajectories.
			# Sadly this means you can't use ptraj, vmd, or anything to read them in until you add it.
			for i in 1 2 3 4 5 6 7 8
			do
			perl -pi -e 'print "\n" if $. == 1' mdcrd.rep$i
			done
			</pre>

			Note that if you have large files, this may take a while - but it is less boring than opening each in vi and adding a line at the top!

			===Clustering structures===

			If you also have the MMTSB Toolset installed - you can cluster the output of your REX runs quite easily. Clustering provides information on the most populated conformational groups, and makes it easy to pick out common structures from the run without manually trawling through every structure one at a time.

import>Bs360 at 11:44, 9 October 2008

2008-10-09T11:44:15Z

← Older revision		Revision as of 11:44, 9 October 2008
Line 75:		Line 75:
	</pre>		</pre>
	Comparing restrt to the original inpcrd file, you'll notice that it is almost twice as long as before! This is because the velocities for each atom are now appended to the bottom of the file. Simply copy them into your inpcrd file and you're good to start the simulation.		Comparing restrt to the original inpcrd file, you'll notice that it is almost twice as long as before! This is because the velocities for each atom are now appended to the bottom of the file. Simply copy them into your inpcrd file and you're good to start the simulation.

	Comparing the restart file produced with

	You can find an example for REMD with AMBER in		You can find an example for REMD with AMBER in
	clust:/sharedscratch/bs360/amberrex_test/		clust:/sharedscratch/bs360/amberrex_test/

	'''THIS EXAMPLE NEEDS UPDATING AS IT DOES NOT HAVE irest=1 AND ntx=5!'''

import>Csw34: /* Determining the replica temperatures */

2008-10-07T14:44:38Z

Determining the replica temperatures

← Older revision		Revision as of 14:44, 7 October 2008
Line 47:		Line 47:

	===Determining the replica temperatures===		===Determining the replica temperatures===
	As we're not using the MMTSB toolset, you will need to calculate the replica temperatures yourself. This is done as follows with <math>n_{tot}=</math>(TOTAL number of replicas ~~- 1)~~, <math>n_{rep}=</math>(current replica number - 1), <math>T_{min}=</math>the lowest temperature used and <math>T_{max}=</math>the highest temperature used.		As we're not using the MMTSB toolset, you will need to calculate the replica temperatures yourself. This is done as follows with <math>n_{tot}=</math>TOTAL number of replicas, <math>n_{rep}=</math>(current replica number - 1), <math>T_{min}=</math>the lowest temperature used and <math>T_{max}=</math>the highest temperature used.

	<math>A</math> is an constant defined as:		<math>A</math> is an constant defined as:

import>Csw34: /* Determining the replica temperatures */

2008-10-07T14:43:32Z

Determining the replica temperatures

← Older revision		Revision as of 14:43, 7 October 2008
Line 47:		Line 47:

	===Determining the replica temperatures===		===Determining the replica temperatures===
	As we're not using the MMTSB toolset, you will need to calculate the replica temperatures yourself. This is done as follows with <math>n_{tot}=</math>TOTAL number of replicas - 1), <math>n_{rep}=</math>(current replica number - 1), <math>T_{min}=</math>the lowest temperature used and <math>T_{max}=</math>the highest temperature used.		As we're not using the MMTSB toolset, you will need to calculate the replica temperatures yourself. This is done as follows with <math>n_{tot}=</math>(TOTAL number of replicas - 1), <math>n_{rep}=</math>(current replica number - 1), <math>T_{min}=</math>the lowest temperature used and <math>T_{max}=</math>the highest temperature used.

	<math>A</math> is an constant defined as:		<math>A</math> is an constant defined as:

import>Csw34: /* Determining the temperatures */

2008-10-07T14:43:02Z

Determining the temperatures

← Older revision		Revision as of 14:43, 7 October 2008
Line 46:		Line 46:
	exponentially distributed between <math>T_{min}</math> and <math>T_{max}</math> as in the MMTSB protocol. Note that with '''irest=1''', each MD cycle will restart the trajectory from the relevant restart file. This is why you can have '''tempi=0''', as AMBER reads the temperature from that file and ignores the MD input.		exponentially distributed between <math>T_{min}</math> and <math>T_{max}</math> as in the MMTSB protocol. Note that with '''irest=1''', each MD cycle will restart the trajectory from the relevant restart file. This is why you can have '''tempi=0''', as AMBER reads the temperature from that file and ignores the MD input.

	===Determining the temperatures===		===Determining the replica temperatures===
	As we're not using the MMTSB toolset, you will need to calculate the replica temperatures yourself. This is done as follows with <math>n_{tot}=</math>TOTAL number of replicas - 1), <math>n_{rep}=</math>(current replica number - 1), <math>T_{min}=</math>the lowest temperature used and <math>T_{max}=</math>the highest temperature used.		As we're not using the MMTSB toolset, you will need to calculate the replica temperatures yourself. This is done as follows with <math>n_{tot}=</math>TOTAL number of replicas - 1), <math>n_{rep}=</math>(current replica number - 1), <math>T_{min}=</math>the lowest temperature used and <math>T_{max}=</math>the highest temperature used.

Line 54:		Line 54:
	The temperature of replica <math>n+1</math> is then defined as:		The temperature of replica <math>n+1</math> is then defined as:

	<math>T_{n+1}=T_{min}\exp{An}</math>		<math>T_{n+1}=T_{min}\exp{(An)}</math>

	===The input.crd files===		===The input.crd files===

import>Csw34 at 14:42, 7 October 2008

2008-10-07T14:42:19Z

← Older revision		Revision as of 14:42, 7 October 2008
Line 44:		Line 44:
	'''ntwr''' determines after now many MD steps the restart files mdinfo.rep* is updated.		'''ntwr''' determines after now many MD steps the restart files mdinfo.rep* is updated.
	The '''temperatures''' for each replica you have to determine yourself, i.e, they are not automatically		The '''temperatures''' for each replica you have to determine yourself, i.e, they are not automatically
	exponentially distributed between ~~T_min~~ and ~~T_max~~ as in the MMTSB protocol. Note that with '''irest=1''', each MD cycle will restart the trajectory from the relevant restart file. This is why you can have '''tempi=0''', as AMBER reads the temperature from that file and ignores the MD input.		exponentially distributed between <math>T_{min}</math> and <math>T_{max}</math> as in the MMTSB protocol. Note that with '''irest=1''', each MD cycle will restart the trajectory from the relevant restart file. This is why you can have '''tempi=0''', as AMBER reads the temperature from that file and ignores the MD input.

			===Determining the temperatures===
			As we're not using the MMTSB toolset, you will need to calculate the replica temperatures yourself. This is done as follows with <math>n_{tot}=</math>TOTAL number of replicas - 1), <math>n_{rep}=</math>(current replica number - 1), <math>T_{min}=</math>the lowest temperature used and <math>T_{max}=</math>the highest temperature used.

			<math>A</math> is an constant defined as:
			<math>A=\frac{ln(T_{max}/T_{min})}{n_{tot}-1}</math>

			The temperature of replica <math>n+1</math> is then defined as:

			<math>T_{n+1}=T_{min}\exp{An}</math>

	===The input.crd files===		===The input.crd files===

import>Csw34 at 13:31, 3 October 2008

2008-10-03T13:31:40Z

← Older revision		Revision as of 13:31, 3 October 2008
Line 6:		Line 6:
	the path names accordingly).		the path names accordingly).
	After producing the config.h file you can compile the MPI version by '''make parallel'''.		After producing the config.h file you can compile the MPI version by '''make parallel'''.
	You should have loaded the pgi64 and mpi/lam/64/pgi/7.1.1 modules.		You should have loaded the pgi64 and mpi/lam/64/pgi/7.1.1 modules. This version has some modifications made to the Sander source code in the files mdfill.f and multisander.f. They are tagged with 'bs360' to make them easy to find. It is not known if they are all necessary, but it works fine with them!

	To run REMD with AMBER you need following '''input files''':		To run REMD with AMBER you need following '''input files''':