Update RAMD to version 5.0.4 54/4554/2
authorDavid <dhardy@ks.uiuc.edu>
Fri, 7 Sep 2018 19:29:48 +0000 (14:29 -0500)
committerDavid <dhardy@ks.uiuc.edu>
Fri, 7 Sep 2018 20:42:02 +0000 (15:42 -0500)
Updated scripts for RAMD to version 5.0.4 that has been tested
successfully with NAMD versions 2.10 and 2.12.  Retained earlier scripts
in directory ramd-4.1.

From the website: t-random acceleration molecular dynamics (tRAMD) is a
protocol for the ranking of drug candidates by their residence time and
obtaining insights into ligand-target dissociation mechanism.
https://www.h-its.org/downloads/ramd/

Contributed by Stefan Richter and Rebecca Wade.

Change-Id: If9c0decf899570d0eb04dd04d03c088701f965fc

14 files changed:
lib/ramd-4.1/README [moved from lib/ramd/README with 100% similarity]
lib/ramd-4.1/examples/README.examples [moved from lib/ramd/examples/README.examples with 100% similarity]
lib/ramd-4.1/examples/example_ramd-md_output/example_ramd-md.namdout [moved from lib/ramd/examples/example_ramd-md_output/example_ramd-md.namdout with 100% similarity]
lib/ramd-4.1/examples/example_ramd_output/example_ramd.namdout [moved from lib/ramd/examples/example_ramd_output/example_ramd.namdout with 100% similarity]
lib/ramd-4.1/examples/run_ramd-md_example.sh [moved from lib/ramd/examples/run_ramd-md_example.sh with 100% similarity]
lib/ramd-4.1/examples/run_ramd_example.sh [moved from lib/ramd/examples/run_ramd_example.sh with 100% similarity]
lib/ramd-4.1/scripts/ramd-4.1.tcl [moved from lib/ramd/scripts/ramd-4.1.tcl with 100% similarity]
lib/ramd-4.1/scripts/ramd-4.1_script.tcl [moved from lib/ramd/scripts/ramd-4.1_script.tcl with 100% similarity]
lib/ramd-4.1/scripts/vectors.tcl [moved from lib/ramd/scripts/vectors.tcl with 100% similarity]
lib/ramd-5.0.4/README [new file with mode: 0644]
lib/ramd-5.0.4/make_distribution.sh [new file with mode: 0755]
lib/ramd-5.0.4/tcl/README [new file with mode: 0644]
lib/ramd-5.0.4/tcl/ramd-5.tcl [new file with mode: 0644]
lib/ramd-5.0.4/tcl/ramd-5_script.tcl [new file with mode: 0644]

similarity index 100%
rename from lib/ramd/README
rename to lib/ramd-4.1/README
diff --git a/lib/ramd-5.0.4/README b/lib/ramd-5.0.4/README
new file mode 100644 (file)
index 0000000..717bb37
--- /dev/null
@@ -0,0 +1,37 @@
+******************************************************************************************\r
+# Random Acceleration Molecular Dynamics (RAMD)                                          *\r
+# Implementation tested on NAMD v2.10 and 2.12                                           *\r
+# December 2017                                                                          *\r
+#                                                                                        *\r
+# Copyright (c) 2017, HITS gGmbH (Heidelberg Institute for Theoretical Studies), Germany * \r
+# Authors:  Vlad Cojocaru, Stefan Richter, Daria Kokh, Rebecca Wade                      *\r
+# Email: mcmsoft@h-its.org                                                               *\r
+#*****************************************************************************************\r
+\r
+\r
+\r
+---------------Random Acceleration Molecular Dynamics------------------------------------\r
+\r
+RAMD 5.0 implementation tested in NAMD 2.6, 2.7b2, 2.10 and 2.12\r
+\r
+Note that to be able to estimate dissociation rate constants (koff rates) using the tauRAMD method as described in the following publication :\r
+\r
+   Kokh et al., J. Chem. Theory Comput. (2018), DOI: 10.1021/acs.jctc.8b00230\r
+   https://pubs.acs.org/doi/full/10.1021/acs.jctc.8b00230\r
+\r
+   one needs to download additional scripts for processing the RAMD trajectories from the HITS website:\r
+   https://www.h-its.org/downloads/ramd/\r
+\r
+\r
+\r
+\r
+Distribution contains:\r
+\r
+1. tcl                 \r
+                   the latest version of the tcl scripts for running RAMD simulations\r
+2. example/1WDHI        \r
+                   an example for execution RAMD simulations; contains input files and bash script "RAMD-force.sh"\r
+3. TRJ-Analysis-R (if downloaded from https://www.h-its.org/downloads/ramd/)   \r
+                   R script for computing relative residence times from RAMD simulations (tauTAMD method); application  example for HSP90 \r
+\r
+\r
diff --git a/lib/ramd-5.0.4/make_distribution.sh b/lib/ramd-5.0.4/make_distribution.sh
new file mode 100755 (executable)
index 0000000..183e535
--- /dev/null
@@ -0,0 +1,10 @@
+#!/bin/bash
+#To extract a release, do the following:
+export SUBREL=5.0.4
+rm -rf ../ramd-${SUBREL}*
+git archive --format=tar --prefix=ramd-${SUBREL}/ master | gzip -9 > ../ramd-${SUBREL}.tgz
+cd .. &&  rm -rf tmp && mkdir tmp && cd tmp
+tar xf ../ramd-${SUBREL}.tgz
+rm -r ramd-${SUBREL}/TRJ-Analysis-R
+zip -rm ../ramd-only-${SUBREL}.zip ramd-${SUBREL}
+cd .. && rmdir tmp
diff --git a/lib/ramd-5.0.4/tcl/README b/lib/ramd-5.0.4/tcl/README
new file mode 100644 (file)
index 0000000..06c763f
--- /dev/null
@@ -0,0 +1,154 @@
+******************************************************************************************\r
+# Random Acceleration Molecular Dynamics (RAMD)                                          *\r
+# Implementation tested on NAMD v2.10 and 2.12                                           *\r
+# December 2017                                                                          *\r
+#                                                                                        *\r
+# Copyright (c) 2017, HITS gGmbH (Heidelberg Institute for Theoretical Studies), Germany * \r
+# Authors:  Vlad Cojocaru, Stefan Richter, Daria Kokh, Rebecca Wade                      *\r
+# Email: mcmsoft@h-its.org                                                               *\r
+#*****************************************************************************************\r
+\r
+Random Acceleration Molecular Dynamics\r
+RAMD 5.0 implementation tested in NAMD 2.6, 2.7b2, 2.10 and 2.12\r
+\r
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\r
+NOTE: It requires NAMD to be compiled against tcl 8.5\r
+NOTE: The parameter forceRAMD replaces the previous acceleration parameter\r
+NOTE: Option of combined RAMD-MD simulations is no longer supported\r
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\r
+\r
+Scripts (in the './tcl' directory):\r
+   ramd-5.tcl        - 'wrapper' script that is sourced in the \r
+                       NAMD configuration file\r
+   ramd-5_script.tcl - script that performs the calculations\r
+\r
+Note that to be able to estimate dissociation rate constants (koff rates) using the tauRAMD method as described in the following publication :\r
+\r
+   Kokh et al., J. Chem. Theory Comput. (2018), DOI: 10.1021/acs.jctc.8b00230\r
+   https://pubs.acs.org/doi/full/10.1021/acs.jctc.8b00230\r
+\r
+   one needs to download additional scripts for processing the RAMD trajectories from the HITS website:\r
+   https://www.h-its.org/downloads/ramd/\r
+\r
+   \r
+Simulation methods implemented:\r
+   (1) Pure RAMD Simulations (RAMD)\r
+   \r
+Reference describing the original method (pure RAMD simulations):\r
+* Luedemann, S.K., Lounnas, V. and R. C. Wade. How do Substrates Enter and Products Exit the Buried Active Site of Cytochrome P450cam ? 1. Random Expulsion Molecular Dynamics Investigation of Ligand Access Channels and Mechanisms. J Mol Biol, 303:797-811 (2000) \r
+(additional references may be found at 'https://www.h-its.org/en/research/mcm/software/ramd/'\r
+\r
+Notes on the current implementation:\r
+(see examples of NAMD configuration files in the ../examples directory for description of input parameters)\r
+\r
+*** version 5.0.3 and 5.0.4 : Cleanup of code and documentation. Added example case.\r
+\r
+*** version 5.0.2 allows to set the filename of the ramd log file. This is \r
+    useful for cases where several simulations start from one directory.\r
+    Corrected the description of the ramd force. Version 4 ramd required \r
+    acceleration magnitude to be defined in the input. In this version, a force magnitude is required as input instead,\r
+    i.e. the input parameter "forceRAMD" defines an external force given in kcal/mol*A\r
+\r
+*** version 5.0.1 adds support for NAMD 2.12. Since this requires the call to an\r
+    additional function (enalbetotalforces), one needs to set  the version\r
+    of NAMD used, in order to decide if this function should be called or not.\r
+    PLEASE SET THE VALUE OF namdVersion as shown below (default is 2.12):\r
+\r
+    ramd namdVersion 2.12\r
+    or \r
+    ramd namdVersion 2.10\r
+\r
+\r
+\r
+*** version 5.0 replaces the vector.tcl script by the internal NAMD vector implementation and other\r
+                performance enhancements.\r
+\r
+*** version 4.1 fixes the following bugs (the implementation and usage is the same as described for version 4.0):\r
+       - the scripts are now also compatible with tcl 8.3 (version 4.0 was only compatible with tcl 8.4)\r
+       - the exit of NAMD is error-free if a ligand exit event occurred\r
+\r
+*** version 4.0 is the first released version\r
+\r
+*** acceleration is applied on the ligand center of mass ( kcal/mol*A*amu );\r
+    ensures compatibility between parameters used for simulations in AMBER 8 and NAMD\r
+\r
+*** in pure RAMD simulations the acceleration 'a' is applied throughout the simulation;\r
+\r
+    1 RAMD block = 'ramdSteps' (input parameter) steps of RAMD simulation\r
+    threshold for the distance traveled by the ligand in 1 RAMD block is 'rRamdMin' (input parameter)\r
+\r
+    the distance traveled by the ligand in 1 RAMD block is 'dr';\r
+    the vector direction of 'a' is changed between 2 RAMD blocks if 'dr' < 'rRamdMin'\r
+    the vector direction of 'a' is changed if 'dr' < 'rRamdMin'\r
+\r
+\r
+*** in combined RAMD-MD simulations, RAMD blocks alternate with standard MD blocks \r
+    ('ramdSteps' and 'mdSteps' input parameters);\r
+    \r
+    1 RAMD block = 'ramdSteps' (input parameter) steps of RAMD simulation\r
+    1 MD block = 'mdSteps' (input parameter) steps of standard MD simulation\r
+    threshold for the distance traveled by the ligand in 1 RAMD block is 'rRamdMin' (input parameter)\r
+    threshold for the distance traveled by the ligand in 1 MD block is 'rMdMin' (input parameter)\r
+\r
+    the distance between the protein and ligand centers of mass is 'd'\r
+    the distance traveled by the ligand in 1 RAMD block is 'dr';\r
+    the distance traveled by the ligand in 1 MD block is 'dm'\r
+\r
+    switch from RAMD to MD\r
+       if 'dr' > 'rRamdMin'\r
+       no switch if 'dr' < 'rRamdMin'\r
+    switch from MD to RAMD\r
+       if 'dm' < 'rMdMin' and 'd' > 0 (acceleration direction kept from previous RAMD block) \r
+       if 'dm' < 'rMdMin' and 'd' < 0 (acceleration direction changed) \r
+       if 'dm' > 'rMdMin' and 'd' < 0 (acceleration direction changed)\r
+       no switch id 'dm' > 'rMdMin' and 'd' > 0 \r
+\r
+*** combined RAMD-MD simulation may be started either with RAMD or with standard MD\r
+    (see the 'mdStart' input parameter )\r
+\r
+*** forces acting on the ligand (total and external) are recorded and printed\r
+    (see input parameter 'forceOutFreq')\r
+    \r
+    the format of the printing is vectorial {Fx Fy Fz} with the length of vector in brackets \r
+    the corresponding position (x, y, z) of the ligand COM is printed \r
+    \r
+    input parameters 'mdSteps' and 'ramdSteps' have to be multiples of 'forceOutFreq' \r
+    (namd exits with error otherwise)\r
+    \r
+    the forces acting on the ligand due to the simulation system can be obtain by substracting\r
+    the external forces from the total forces (this feature will be implemented in a future release; \r
+    in the meantime scripts can be written to process the detailed output of this version and \r
+    get exactly this information)\r
+     \r
+    in combined RAMD-MD simulation\r
+    NAMD exits with an error if external forces are found during the standard MD steps\r
+    \r
+    the forces are written in the default NAMD configuration file; the output is labeled with 'FORCE'\r
+    \r
+*** the output of RAMD and RAMD-MD simulations is dumped in the default NAMD log file\r
+    and is labeled with  "RAMD" or "MD" depending whether it is written in a RAMD or MD step\r
+    \r
+*** debugging option implemented\r
+    (see input parameter 'debugLevel')\r
+    \r
+    this parameter is strictly only recommended for short test runs. \r
+    The output is very rich in information and for a long run is confusing. \r
+    \r
+    The following information is printed:\r
+       comparison of the added external force with the external force printed by 'loadforces' for each atom\r
+          (they should have equal values)   \r
+       counting of the timestep, md steps and ramd steps per block (to make sure it does it correclty)\r
+    \r
+    the output is labelled with 'RAMD DEBUG' or 'MD DEBUG'\r
+    \r
+*** Simulations stop when a maximum distance between the ligand and protein COMs\r
+    is reached ('maxDist' input parameter)    \r
+\r
+Example runs for testing are provided in the './examples' directory\r
+   File './examples/example_ramd.namdin' is an example of NAMD configuration file for pure RAMD simulation \r
+   File './examples/example_ramd-md.namdin' is an example of NAMD configuration file for combined RAMD-MD simulation \r
+\r
+IMPORTANT NOTE\r
+There is a significant difference between ligand exit times depending on whether Langevin dynamics or\r
+the weak coupling algorithm of Berendsen are used for temperature control. Users are encouraged to run\r
+test runs using both methods and decide what is appropriate for their specific simulation system and setup.\r
diff --git a/lib/ramd-5.0.4/tcl/ramd-5.tcl b/lib/ramd-5.0.4/tcl/ramd-5.tcl
new file mode 100644 (file)
index 0000000..8d4988c
--- /dev/null
@@ -0,0 +1,109 @@
+#***********************************************************************************\r
+#                                                                                  *\r
+# Random Acceleration Molecular Dynamics (RAMD)                                    *\r
+# Implementation for NAMD v2.12                                                    *\r
+# Version 5.0.4 July 2018                                                     *\r
+#                                                                                  *\r
+# Copyright (c) 2017, HITS gGmbH, Heidelberg, Germany                              * \r
+# Authors: Vlad Cojocaru, Stefan Richter, Daria Kokh, Rebecca Wade                 *\r
+# Email: mcmsoft@h-its.org                                                         *\r
+#                                                                                  *\r
+# The first Tcl script to run RAMD in NAMD (v2.5+) was written by Harish Vashisth  *\r
+# Ref: Vashisth H et al, Biophys J. 2008 Nov 1;95(9):4193-204. Epub 2008 Aug 1     *\r
+#                                                                                  *\r
+# The Vashisth script inspired some lines of this script, but mainly this script   *\r
+#    is based on the fortran code written for AMBER 8                              *\r
+#                                                                                  *\r
+# The structure of this script is inspired by the Adaptive Biasing Force module    *\r
+#    distributed with NAMD 2.6                                                     *\r
+#                                                                                  *\r
+# The original RAMD method is described in:                                        *\r
+#    Ref1: Luedemann,S.K.,Lounnas,V.and R.C.Wade.,                                 *\r
+#          J Mol Biol, 303:797-811 (2000)                                          *\r
+#    Ref2: Schleinkofer,K.,Sudarko,Winn,P.,Luedemann,S.K.and R.C.Wade,             *\r
+#          EMBO Reports, 6, 584-589 (2005)                                         *\r
+#                                                                                  *\r
+#                                                                                  *\r
+# Disclaimer:  This script is for research purposes only. EML Research does not    *\r
+#              assume any responsibility for the software or its use.              * \r
+#                                                                                  *\r
+#   The script along with usage examples is available at                           *\r
+#   https://www.h-its.org/en/research/mcm/software/                                *\r
+#***********************************************************************************\r
+\r
+#*******************************************************\r
+# Startup                                              *\r
+#*******************************************************\r
+package provide ramd 5.0\r
+\r
+#*******************************************************\r
+# Parameter definitions\r
+#*******************************************************\r
+\r
+namespace eval ::RAMD {\r
+ set version "5.0"\r
+ if {! [info exists RAMDdir]} { set RAMDdir [file dirname [info script]] }\r
+ # If it fails, try the local directory\r
+ if { $RAMDdir == "" } { set RAMDdir "." }\r
\r
+ TclForces             on\r
+ TclForcesScript       $RAMDdir/ramd-5_script.tcl\r
+ array set defaults {\r
+  ramdSteps               50\r
+  namdVersion             2.12\r
+  forceRAMD                16.0\r
+  rMinRamd                 0.01\r
+  forceOutFreq             10\r
+  firstProtAtom            1\r
+  ramdSeed             14253\r
+  mdSteps                  0\r
+  mdStart                 no\r
+  maxDist                 50\r
+  debugLevel               0\r
+  ramdfilename         "ramd.log"\r
+ }\r
+\r
+ set mandatory "firstRamdAtom lastRamdAtom lastProtAtom"\r
+ set silent "rMinMd"\r
+\r
+ array set capitals {}\r
+ foreach param [concat $mandatory $silent [array names defaults]] {\r
+  set capitals([string tolower $param]) $param\r
+  # not set yet\r
+  set alreadySet($param) 0\r
+ }\r
+} ;# namespace\r
+\r
+\r
+proc ramd { keyword value } {\r
+ set ::RAMD::keyword $keyword\r
+ set ::RAMD::value $value\r
+\r
+ namespace eval ::RAMD {\r
+  # Build list of all allowed parameter names\r
+  set list [array names capitals]\r
+  set lowercase [string tolower $keyword]\r
+  # Process parameters\r
+  if {[lsearch $list $lowercase] != -1} {\r
+   set keyword $capitals($lowercase)\r
+   if { $alreadySet($keyword) } {\r
+    print "RAMD> WARNING - multiple definitions of parameter $keyword" \r
+   }\r
+   set $keyword $value\r
+   set alreadySet($keyword) 1\r
+   return\r
+  } else {\r
+   error [format "Unknown RAMD keyword: %s" $keyword]\r
+  }\r
+ } ;# namespace\r
+\r
+} ;# proc ramd\r
+\r
+# define upper-case synonyms to proc abf \r
+proc RAMD { keyword value } {\r
+    ramd $keyword $value\r
+}\r
+proc Ramd { keyword value } {\r
+    ramd $keyword $value\r
+}\r
+\r
diff --git a/lib/ramd-5.0.4/tcl/ramd-5_script.tcl b/lib/ramd-5.0.4/tcl/ramd-5_script.tcl
new file mode 100644 (file)
index 0000000..99a4d9b
--- /dev/null
@@ -0,0 +1,484 @@
+#***********************************************************************************\r
+#                                                                                  *\r
+# Random Acceleration Molecular Dynamics (RAMD)                                    *\r
+# Implementation for NAMD v2.12                                                    *\r
+# Version 5.03, September 2017                                                     *\r
+#                                                                                  *\r
+# Copyright (c) 2017, HITS gGmbH, Heidelberg, Germany                              * \r
+# Authors: Vlad Cojocaru, Stefan Richter, Daria Kokh, Rebecca Wade                 *\r
+# Email: mcmsoft@h-its.org                                                         *\r
+#                                                                                  *\r
+# The first Tcl script to run RAMD in NAMD (v2.5+) was written by Harish Vashisth  *\r
+# Ref: Vashisth H et al, Biophys J. 2008 Nov 1;95(9):4193-204. Epub 2008 Aug 1     *\r
+#                                                                                  *\r
+# The Vashisth script inspired some lines of this script, but mainly this script   *\r
+#    is based on the fortran code written for AMBER 8                              *\r
+#                                                                                  *\r
+# The structure of this script is inspired by the Adaptive Biasing Force module    *\r
+#    distributed with NAMD 2.6                                                     *\r
+#                                                                                  *\r
+# The original RAMD method is described in:                                        *\r
+#    Ref1: Luedemann,S.K.,Lounnas,V.and R.C.Wade.,                                 *\r
+#          J Mol Biol, 303:797-811 (2000)                                          *\r
+#    Ref2: Schleinkofer,K.,Sudarko,Winn,P.,Luedemann,S.K.and R.C.Wade,             *\r
+#          EMBO Reports, 6, 584-589 (2005)                                         *\r
+#                                                                                  *\r
+# The present script is a slightly modified version of ramd-5_script.tcl           *\r
+#  adjusted to use the random force magnitude as  input for the simulations        *\r
+#  instead of the acceleration                                                     *\r
+#                                                                                  *\r
+# Disclaimer:  This script is for research purposes only. HITS gGmbH does not      *\r
+#              assume any responsibility for the software or its use.              * \r
+#                                                                                  *\r
+# !!! Quantitative reproducibility of the results obtained with AMBER 8            *\r
+#     is not possible due to a numerical errors and such like                      * \r
+#                                                                                  *\r
+# !!! This script is under development.                                            *\r
+#                                                                                  *\r
+#   The script along with usage examples is available at                           *\r
+#   https://www.h-its.org/en/research/mcm/software/                                *\r
+#***********************************************************************************\r
+namespace eval ::RAMD {\r
+# set ramdfilename "ramd.log"\r
+ set ramdfileid [open $ramdfilename "w"]\r
+ fconfigure $ramdfileid -blocking 0\r
+ if {$namdVersion > 2.10} {\r
+   enabletotalforces\r
+ }\r
+ puts $ramdfileid "RAMD:"\r
+ puts $ramdfileid "RAMD:   -------------------------------------------------------------------  "\r
+ puts $ramdfileid "RAMD:   Random Acceleration Molecular Dynamics Simulation version $version"\r
+ puts $ramdfileid "RAMD:   -------------------------------------------------------------------  "\r
+ puts $ramdfileid "RAMD:"\r
+ #***** Assign default values for the parameters not specified in the configuration file\r
+ foreach option [array names defaults] {\r
+  if {! [info exists $option]} {\r
+   set $option $defaults($option)\r
+   puts $ramdfileid [format "RAMD: %25s = %s" $option [expr $$option]]\r
+  } elseif { [info exists $option] } {\r
+   puts $ramdfileid [format "RAMD: %25s = %s" $option [expr $$option]]\r
+  }\r
+ }\r
+ #***** Check if mandatory parameters are specified in the configuration file\r
+ foreach var $mandatory {\r
+  if {! [info exists $var]} {\r
+   error "RAMD: Mandatory parameter '$var' is not set -- cannot start RAMD"\r
+  } else {\r
+   puts $ramdfileid [format "RAMD: %25s = %s" $var [expr $$var]]\r
+  }\r
+ }\r
+ #***** Check if 'forceOutFreq' is equal to 1; exit with error if that's the case\r
+ if { $forceOutFreq == 1 } { error "RAMD: ERROR: 'forceOutFreq' parameter may not be 1" } \r
+\r
+ #***** Check if 'mdSteps' is specified in the configuration file\r
+\r
+ #***** Performed pure RAMD if 'mdSteps' = 0\r
+ if { $mdSteps == 0 } { \r
+  \r
+  #***** Check that the number of ramd steps is a multiple of 'forceOutFreq'; exit with error if not\r
+  set r [expr "$ramdSteps % $forceOutFreq"]\r
+  if { $r != 0 } { error "RAMD: ERROR: The number of RAMD steps is not a multiple of 'forceOutFreq'" } \r
\r
+  puts $ramdfileid "RAMD: Pure RAMD simulation is performed" \r
+  \r
+  #***** If 'mdSteps' is 0 and "mdStart" is yes, give a warning\r
+  if { $mdStart == "yes" } { \r
+   puts $ramdfileid "RAMD: WARNING: 'mdStart' has no meaning for pure RAMD simulation; it will be ignored" \r
+  }\r
+\r
+  #***** If 'mdSteps' is 0 and "rMinMd" is set, give a warning\r
+  if { [info exists rMinMd] } {\r
+   puts $ramdfileid "RAMD: WARNING: 'rMinMd' specified while 'mdSteps' is 0"\r
+   puts $ramdfileid "RAMD: WARNING: For combined RAMD-MD simulation 'mdSteps' must be greater than 0"\r
+   puts $ramdfileid "RAMD: WARNING: Ignore 'rMinMd' and perform pure RAMD simulation"\r
+  }\r
+  \r
+ }\r
+\r
+ #***** Perform combined RAMD with MD simulation if 'mdSteps' is not 0 and 'rMinMd' is specified\r
+ if { $mdSteps != 0  } { \r
+  error "To run ramd setting mdSteps is not supported"\r
+\r
+ }\r
+     \r
+ #***** Make a list of all the atoms on which the force will be applied\r
+ set ramdAtoms {}\r
+ for { set i $firstRamdAtom } { $i <= $lastRamdAtom } { incr i } { lappend ramdAtoms $i }\r
+ puts $ramdfileid "RAMD: Atoms subject to the random force are: $ramdAtoms"\r
+ foreach ramdAtom $ramdAtoms { addatom $ramdAtom }\r
+ #***** Define a group of the ligand atoms; the force will be applied on the center of mass of this group\r
+ set ramdGroup [ addgroup $ramdAtoms ]\r
+\r
+ #***** Define a group containing all protein atoms\r
+ set protAtoms {}\r
+ for { set i $firstProtAtom } { $i <= $lastProtAtom } { incr i } { lappend protAtoms $i }\r
+ foreach protAtom $protAtoms { addatom $protAtom }\r
+ set protGroup [ addgroup $protAtoms ]\r
+\r
+ #***** Some variable initialization \r
+ set timeStep 0; set exitFlag 0; \r
+ set prevLigCOM "0.0 0.0 0.0"; set prevProtCOM "0.0 0.0 0.0"; set prevDist 0;\r
+\r
+ #***** Initialization of simulation flags\r
+ if { $mdSteps == 0 } {\r
+  set ramdFlag 1; set mdFlag 0; set ramdStep 0; set mdStep 0;\r
+ }\r
\r
+} ;# namespace\r
\r
+#***** In root namespace (::) for all procedures we have to add the following procedure definition\r
+#proc veclength {v} {\r
+# return [expr {sqrt([veclength2 $v])}]\r
+#}\r
+#***** Source the vectors and matrices procedures from VMD\r
+#source $RAMD::RAMDdir/vectors.tcl\r
+\r
+#***********************************************************\r
+# PROCEDURE TO GENERATE RANDOMLY ORIENTED ACCELERATION \r
+#***********************************************************\r
+proc genRandAccel { timeStep } {\r
+namespace eval ::RAMD {\r
+\r
+ set pi 3.141592653589793238462643383\r
+# set pi [expr "2.0*asin(1.0)"]\r
+\r
+ #***** Generate new random orientation of the ramd force\r
+# set randTheta [expr "rand()"]\r
+# set randPsi [expr "rand()"]\r
+ set theta [expr "2*$pi*rand()"]\r
+ set psi [expr "$pi*rand()"]\r
+ set sinpsi [expr "sin($psi)"]\r
+ set rx [expr "cos($theta)*$sinpsi"]\r
+ set ry [expr "sin($theta)*$sinpsi"]\r
+ set rz [expr "cos($psi)"]\r
+ set r "$rx $ry $rz"\r
+ set lenr [veclength $r]\r
+\r
+ # Force given in kcal/mol*A  in the NAMD configuration file \r
+ set vecAccel [vecscale [expr "$forceRAMD"] $r ]\r
+  \r
+ return \r
\r
+} ;# namespace\r
+} ;# proc genRandAccel {timestep}\r
+\r
+\r
+#*****************************************************************************\r
+# PROCEDURE TO EVALUATE THE DISTANCE TRAVELLED BY THE LIGAND IN N RAMD STEPS\r
+#*****************************************************************************\r
+proc evalWalkDist { timeStep prevLigCOM prevProtCOM currLigCOM currProtCOM } {\r
+namespace eval ::RAMD {\r
\r
+ #***** Compute the relative position of the ligand com with regard to the protein com\r
+ set prevRelLigCOM [ vecsub $prevLigCOM $prevProtCOM ]\r
+ set currRelLigCOM [ vecsub $currLigCOM $currProtCOM ]\r
\r
+ #***** Compute the distance travelled by the ligand com during a ramd or md stint\r
+ set vecWalkDist [vecsub $currRelLigCOM $prevRelLigCOM]\r
+ set walkDist [veclength $vecWalkDist]\r
+\r
+# set vecWalkDistX [lindex $vecWalkDist 0]\r
+# set vecWalkDistY [lindex $vecWalkDist 1]\r
+# set vecWalkDistZ [lindex $vecWalkDist 2]\r
\r
+ return  \r
+\r
+} ;# namespace\r
+} ;# proc evalWalkDist\r
+  \r
+\r
+#**************************************************************\r
+# PROCEDURE TO APPLY THE FORCE WHICH IS CALLED EVERY TIME STEP \r
+#**************************************************************\r
+proc calcforces {} {\r
+namespace eval ::RAMD {\r
+ #***** Terminate NAMD if the ligand has exited from the protein\r
+ if { $exitFlag == 1 } {\r
+  puts $ramdfileid "EXIT: $timeStep  > MAX DISTANCE LIGAND COM - PROTEIN COM REACHED"\r
+  puts $ramdfileid "EXIT: $timeStep  > LIGAND EXIT EVENT DETECTED: STOP SIMULATION"\r
+  puts $ramdfileid "EXIT: $timeStep  > EXIT NAMD"\r
+  fconfigure $ramdfileid -blocking 1\r
+  close $ramdfileid\r
+  puts "EXIT: $timeStep  > MAX DISTANCE LIGAND COM - PROTEIN COM REACHED"\r
+  puts "EXIT: $timeStep  > LIGAND EXIT EVENT DETECTED: STOP SIMULATION"\r
+  puts "EXIT: $timeStep  > EXIT NAMD"\r
+  set process [pid]\r
+  puts "Killing $process"\r
+  exec kill -15 $process\r
+#  exec kill -15 $$\r
+  exit 0\r
+#  exec /bin/bash qdel_jobid.sh\r
+ } \r
+\r
+ if { [ array exists coords ] } { array unset coords }\r
+# if { [ array exists masses ] } { array unset masses }\r
+# if { [ array exists extForces ] } { array unset extForces }\r
+ if { [ array exists totForces ] } { array unset totForces }\r
+\r
+ #***** Load coordinates for all the atoms and groups defined\r
+ loadcoords coords\r
+ #***** Load masses for all the atoms and groups defined\r
+ if  {[array exists masses] == 0} { loadmasses masses }\r
+ #***** Calculate the mass of the ligand\r
+ if {[info exists ligMass] == 0} {\r
+    set ligMass 0\r
+    foreach ramdAtom $ramdAtoms {\r
+      set ligMass [expr $ligMass + $masses($ramdAtom)]\r
+    }\r
+    puts $ramdfileid "RAMD: calculated ligand mass $ligMass"\r
+    foreach ramdAtom $ramdAtoms {\r
+      set relAtomMass($ramdAtom) [expr "$masses($ramdAtom)/$ligMass"]\r
+    }\r
+ }\r
+\r
+ #***** Load external forces from previous time step for all the atoms and groups defined\r
+# loadforces extForces\r
+ #***** Load total forces from previous time step for all the atoms and groups defined\r
+ loadtotalforces totForces \r
\r
+\r
+ #***** Calculate the position of protein and ligand COM\r
+ set protCOM "$coords($protGroup)"\r
+ set ligCOM "$coords($ramdGroup)"\r
+  \r
+ #***** Initialize ramd simulation or combined ramd-md simulation that begins with ramd\r
+ if { $timeStep == 0 } {\r
+  \r
+  expr "srand($ramdSeed)"\r
+  \r
+  set vMin [ expr "($rMinRamd)/($ramdSteps)" ]\r
+  if { $mdSteps == 0 } { \r
+   puts $ramdfileid "RAMD: $timeStep  ***** INITIALIZE RAMD SIMULATION *****" \r
+  } else { \r
+   puts $ramdfileid "RAMD: $timeStep  ***** INITIALIZE COMBINED RAMD-MD SIMULATION *****" \r
+  }\r
+  puts $ramdfileid "RAMD: $timeStep     >>> minimum travelled distance (A): $rMinRamd"\r
+  puts $ramdfileid "RAMD: $timeStep     >>> minimum velocity (A/fs): $vMin"\r
+\r
+  #***** Initial com positions\r
+  set currLigCOM $ligCOM; set currProtCOM $protCOM  \r
+  puts $ramdfileid "RAMD: $timeStep     >>> LIGAND COM IS: $currLigCOM"\r
+  puts $ramdfileid "RAMD: $timeStep     >>> PROTEIN COM IS: $currProtCOM"\r
+\r
+  #***** Evaluate initial distance between ligand com and protein com\r
+  set currDist [veclength [vecsub $currLigCOM $currProtCOM]]\r
+  puts $ramdfileid "RAMD: $timeStep     >>> DISTANCE LIGAND COM - PPROTEIN COM IS: DIST = $currDist"\r
+\r
+  #***** Generate an initial orientation for the external force to be applied when ramd is switched on\r
+  genRandAccel $timeStep\r
+  puts $ramdfileid "RAMD: $timeStep     >>> INITIAL RANDOM DIRECTION: $r :: ||r|| = $lenr"\r
+\r
+  puts $ramdfileid "RAMD: $timeStep  ***** START WITH $ramdSteps STEPS OF RAMD SIMULATION *****"\r
\r
+  #***** Reset the positions of the ligand and protein COMs and the distance ligand com - protein com\r
+  set prevLigCOM "$currLigCOM"; set prevProtCOM "$currProtCOM"; set prevDist "$currDist"\r
+\r
+  incr timeStep\r
+  return  \r
+\r
+ } \r
+\r
\r
+ #***** Perform ramd simulation \r
+ if { $timeStep != 0 && $exitFlag == 0  } {\r
+\r
+  #***** Count ramd steps\r
+  incr ramdStep\r
+#  if { $debugLevel != 0 } {\r
+#   puts $ramdfileid "RAMD DEBUG: TIMESTEP IS: $timeStep; RAMD STEP IS: $ramdStep; MD STEP IS: $mdStep"\r
+#  }\r
+    \r
+  #***** Define and apply the force to each atom of the ligand\r
+#  foreach ramdAtom $ramdAtoms {\r
+#   set atomMass "$masses($ramdAtom)"\r
+#   set atomForce [ vecscale $atomMass $vecAccel ]\r
+#   set atomForce [ vecscale "$relAtomMass($ramdAtom)" $vecAccel ]\r
+#   addforce $ramdAtom $atomForce\r
+#   if { $debugLevel != 0 } { \r
+#      set atomForceValue [ veclength "$atomForce" ]\r
+#      puts "RAMD DEBUG: ATOM $ramdAtom: MASS $atomMass: ADD FORCE $atomForceValue" \r
+#      unset atomForceValue\r
+#   }\r
+#   unset atomForce\r
+#   unset atomMass\r
+#  }\r
+\r
+  #***** Define the force vector that is applied to the center of mass of the ligand\r
+#  set force [vecscale $ligMass $vecAccel]\r
+#  set force $vecAccel\r
+#  set fx [lindex $force 0]\r
+#  set fy [lindex $force 1]\r
+#  set fz [lindex $force 2]\r
+#  set newaccel [ vecscale [expr "1.0/$ligMass"] $vecAccel ]\r
\r
+  #***** Check the magnitude of the force\r
+#  set f [expr "sqrt((($fx)*($fx)+($fy)*($fy)+($fz)*($fz)))"]\r
+   set f [veclength $vecAccel]\r
+  \r
+  #***** Set flag for writing force output\r
+  if { $forceOutFreq != 0 } { set outputFlag [expr "$ramdStep % $forceOutFreq"] }\r
+  \r
+  #***** Write force output every 'forceOutFreq' steps \r
+  if {  $outputFlag == 0 } {\r
+      \r
+   #***** Write the force vector and direction\r
+   puts $ramdfileid "RAMD FORCE: $timeStep  > LIGAND COM is: $ligCOM\nRAMD FORCE: $timeStep  > PROTEIN COM IS $protCOM\nEXTERNAL FORCE VECTOR (F): $vecAccel; ||F|| = $f\nRAMD FORCE: $timeStep  > EXTERNAL FORCE DIRECTION (r): $r; ||r|| = $lenr"\r
+ #  puts $ramdfileid "RAMD FORCE: $timeStep  > EXTERNAL ACCELERATION VECTOR (F): ||Accel|| = $newaccel"\r
+   }\r
+   #***** Calculate external and total forces acting on the ligand\r
+#   set totLigForceX 0; set totLigForceY 0; set totLigForceZ 0; set totLigForceV 0;\r
+#   set extLigForceX 0; set extLigForceY 0; set extLigForceZ 0; set extLigForceV 0;\r
+   set totLigForce "0 0 0"\r
+#   set extLigForce "0 0 0"\r
+   foreach ramdAtom $ramdAtoms {\r
+     set atomForce [ vecscale "$relAtomMass($ramdAtom)" $vecAccel ]\r
+     addforce $ramdAtom $atomForce\r
+#    set atomMass "$masses($ramdAtom)"\r
+    set totAtomForce "$totForces($ramdAtom)"\r
+#    set totAtomForceX [ lindex $totAtomForce 0 ]\r
+#    set totAtomForceY [ lindex $totAtomForce 1 ]\r
+#    set totAtomForceZ [ lindex $totAtomForce 2 ]\r
+#    set totAtomForceV [ veclength "$totAtomForce" ]\r
+    set totLigForce [vecadd $totLigForce $totAtomForce]\r
+#    set totLigForceX [ expr "$totLigForceX + $totAtomForceX" ]\r
+#    set totLigForceY [ expr "$totLigForceY + $totAtomForceY" ]\r
+#    set totLigForceZ [ expr "$totLigForceZ + $totAtomForceZ" ]\r
+#    set totLigForceV [ expr "$totLigForceV + $totAtomForceV" ]\r
+#    set extAtomForce "$extForces($ramdAtom)"\r
+#    set extAtomForceX [ lindex "$extAtomForce" 0 ]\r
+#    set extAtomForceY [ lindex "$extAtomForce" 1 ]\r
+#    set extAtomForceZ [ lindex "$extAtomForce" 2 ]\r
+#    set extAtomForceV [ veclength "$extAtomForce" ]\r
+#    set extLigForceX [ expr "$extLigForceX + $extAtomForceX" ]\r
+#    set extLigForceY [ expr "$extLigForceY + $extAtomForceY" ]\r
+#    set extLigForceZ [ expr "$extLigForceZ + $extAtomForceZ" ]\r
+#    set extLigForce [vecadd $extLigForce $extAtomForce]\r
+#    set extLigForceV [ expr "$extLigForceV + $extAtomForceV" ]\r
+#    if { $debugLevel != 0 } { \r
+#     puts $ramdfileid "RAMD DEBUG: ATOM $ramdAtom: MASS $atomMass: EXT FORCE $extAtomForceV"\r
+#     puts $ramdfileid "RAMD DEBUG: ATOM $ramdAtom: MASS $atomMass: TOT FORCE $totAtomForceV"\r
+#    }    \r
+#    unset atomMass\r
+#    unset totAtomForceX; unset totAtomForceY; unset totAtomForceZ; unset totAtomForceV; unset totAtomForce\r
+#    unset extAtomForceX; unset extAtomForceY; unset extAtomForceZ; unset extAtomForceV; unset extAtomForce\r
+   }\r
+#   set totLigForce "$totLigForceX $totLigForceY $totLigForceZ"\r
+#   set extLigForce "$extLigForceX $extLigForceY $extLigForceZ"\r
\r
+   #***** Write external forces acting on the ligand com for debugging purposes\r
+#   if { $debugLevel !=0 } {\r
+#    set extLigForceV [veclength $extLigForce]\r
+#    puts $ramdfileid "RAMD DEBUG: $timeStep  > EXTERNAL FORCE ON THE LIGAND COM IS: $extLigForce ($extLigForceV)"\r
+#    unset extLigForceV; \r
+#   }\r
+\r
+   set totLigForceV [veclength $totLigForce]\r
+   #***** Write total forces acting on the ligand com\r
+  if {  $outputFlag == 0 } {\r
+\r
+   puts $ramdfileid "RAMD FORCE: $timeStep  > TOTAL FORCE ON THE LIGAND COM IS: $totLigForce ($totLigForceV)"\r
+  }   \r
+   unset totLigForce; #unset extLigForce; \r
+   unset totLigForceV;\r
+\r
\r
+#  elseif { !  [ array exists totForces ] && $outputFlag == 0 } {\r
+  \r
+#   error "RAMD: $timeStep  > ERROR: EXTERNAL FORCES NOT PRESENT DURING RAMD STEP: EXIT NAMD"\r
+  \r
+#  }\r
+\r
+  #***** Set flag for evaluating ramd simulation\r
+  set evalRamdFlag [expr "$ramdStep % $ramdSteps"]\r
+\r
+  #***** Evaluate ramd stint\r
+  if { $evalRamdFlag == 0 } {\r
+  \r
+   puts $ramdfileid "RAMD: $timeStep  ***** EVALUATE $ramdSteps RAMD STEPS AT TIMESTEP $timeStep *****"\r
+\r
+   #***** com positions\r
+   set currLigCOM $ligCOM; set currProtCOM $protCOM\r
+#   if { $debugLevel !=0 } {\r
+#    puts $ramdfileid "RAMD DEBUG: $timeStep  > PREVIOUS LIGAND COM IS: $prevLigCOM"\r
+#    puts $ramdfileid "RAMD DEBUG: $timeStep  > PREVIOUS PROTEIN COM IS: $prevProtCOM"\r
+#    puts $ramdfileid "RAMD DEBUG: $timeStep  > CURRENT LIGAND COM IS: $currLigCOM"\r
+#    puts $ramdfileid "RAMD DEBUG: $timeStep  > CURRENT PROTEIN COM IS: $currProtCOM"\r
+#   }\r
+  \r
+   #***** Evaluate distance between ligand com and protein com\r
+   set currDist [veclength [vecsub $currLigCOM $currProtCOM]]   \r
+   #***** Evaluate the change in the distance between the protein and the ligand com during the ramd stint \r
+   set diffDist [expr "${currDist}-${prevDist}"]\r
+   puts $ramdfileid "RAMD: $timeStep     >>> DISTANCE LIGAND COM - PPROTEIN COM IS: $currDist; IT CHANGED BY $diffDist"\r
+\r
+   #***** Check if the ligand has exited the protein\r
+   if { $currDist >= $maxDist } { set exitFlag 1; return }\r
+   \r
+   #***** Compute the distance travelled by the ligand com during the ramd stint\r
+   evalWalkDist $timeStep $prevLigCOM $prevProtCOM $currLigCOM $currProtCOM\r
\r
+   #***** Evaluate whether a new force direction will be generated\r
+   if { $walkDist <= $rMinRamd } {    \r
+\r
+    genRandAccel $timeStep\r
+\r
+    puts $ramdfileid "RAMD: $timeStep     >>> THE DISTANCE TRAVELLED BY THE LIGAND IS: $walkDist (< $rMinRamd)\nRAMD: $timeStep     >>> CONTINUE WITH $ramdSteps STEPS OF RAMD SIMULATION\nRAMD: $timeStep     >>> CHANGE ACCELERATION DIRECTION TO: $r; ||r|| = $lenr"\r
+    \r
+    #***** Reset the ramd step count\r
+    #***** Reset the positions of the ligand and protein COMs\r
+    set ramdStep 0; set prevLigCOM "$currLigCOM"; set prevProtCOM "$currProtCOM"; set prevDist "$currDist"\r
\r
+    #***** Increment time step and go to the next time step right now\r
+    incr timeStep\r
+    return\r
+  \r
+   } elseif { $walkDist > $rMinRamd } {\r
+    puts $ramdfileid "RAMD: $timeStep     >>> THE DISTANCE TRAVELLED BY THE LIGAND IS: $walkDist (> $rMinRamd)\nRAMD: $timeStep     >>> CONTINUE WITH $ramdSteps STEPS OF RAMD SIMULATION\nRAMD: $timeStep     >>> KEEP PREVIOUS ACCELERATION DIRECTION: $r; ||r|| = $lenr"\r
+   \r
+    #***** Reset the positions of the ligand and protein COMs\r
+    set prevLigCOM "$currLigCOM"; set prevProtCOM "$currProtCOM"; set prevDist "$currDist"\r
\r
+    #***** Increment time step and go to the next time step right now\r
+    incr timeStep\r
+    return\r
+  \r
+   }  \r
+   #***** Ensure that the positions of the ligand and protein COMs are reset after the evaluation\r
+   set prevLigCOM "$currLigCOM"; set prevProtCOM "$currProtCOM"; set prevDist "$currDist"\r
+\r
+   #***** Increment time step and go to the next time step right now\r
+   incr timeStep\r
+   return\r
+  \r
+  }\r
+\r
+  #***** Unset the force\r
+#  if { [info exists force] } { unset force }\r
+#  if { [info exists fx] } { unset fx }\r
+#  if { [info exists fy] } { unset fy }\r
+#  if { [info exists fz] } { unset fz }\r
+#  if { [info exists f] } { unset f }\r
+#  if { [info exists totLigForce] } { unset totLigForce }\r
+#  if { [info exists extLigForce] } { unset extLigForce }\r
+#  if { [info exists totLigForceV] } { unset totLigForceV }\r
+#  if { [info exists extLigForceV] } { unset extLigForceV }\r
+  \r
+  #***** Increment time step and go to the next time step right now\r
+  incr timeStep\r
+  return\r
+   \r
+ }\r
+\r
+\r
\r
+\r
+ #***** Increment time step and go to the next time step right now\r
+ incr timeStep\r
+ return\r
+\r
+} ;# namespace\r
+} ;# proc calcforces {}\r
\r
+#*****************\r
+# END\r
+#*****************\r