Bias-exchange metadynamics implementation

colvartools/namd/Common/ammonia_water.pdb

@@ -0,0 +1,9 @@
+CRYST1    0.000    0.000    0.000  90.00  90.00  90.00 P 1           1
+ATOM      1  N1  AMM1A   1       0.597   0.171   0.242  0.00  0.00      AMM1  
+ATOM      2  H11 AMM1A   1       0.413   0.867   0.956  0.00  0.00      AMM1  
+ATOM      3  H12 AMM1A   1      -0.060   0.366  -0.506  0.00  0.00      AMM1  
+ATOM      4  H13 AMM1A   1       0.311  -0.715   0.642  0.00  0.00      AMM1  
+ATOM      5  OH2 TIP3W   1       3.324   0.855  -0.028  0.00  0.00      WAT   
+ATOM      6  H1  TIP3W   1       2.396   0.647   0.073  0.00  0.00      WAT   
+ATOM      7  H2  TIP3W   1       3.640   0.237  -0.687  0.00  0.00      WAT   
+END

colvartools/namd/Common/ammonia_water.psf

@@ -0,0 +1,45 @@
+PSF EXT
+
+         5 !NTITLE
+ REMARKS original generated structure x-plor psf file
+ REMARKS topology charmm-ff/top_all36_cgenff.rtf 
+ REMARKS topology charmm-ff/toppar_water_ions.str 
+ REMARKS segment AMM1 { first NONE; last NONE; auto angles dihedrals }
+ REMARKS segment WAT { first NONE; last NONE; auto none  }
+
+         7 !NATOM
+         1 AMM1     1        AMM1     N1       NG331   -1.125000       14.0070           0
+         2 AMM1     1        AMM1     H11      HGPAM3   0.375000        1.0080           0
+         3 AMM1     1        AMM1     H12      HGPAM3   0.375000        1.0080           0
+         4 AMM1     1        AMM1     H13      HGPAM3   0.375000        1.0080           0
+         5 WAT      1        TIP3     OH2      OT      -0.834000       15.9994           0
+         6 WAT      1        TIP3     H1       HT       0.417000        1.0080           0
+         7 WAT      1        TIP3     H2       HT       0.417000        1.0080           0
+
+         6 !NBOND: bonds
+         1         2         1         3         1         4         5         6
+         5         7         6         7
+
+         4 !NTHETA: angles
+         2         1         4         2         1         3         3         1         4
+         6         5         7
+
+         0 !NPHI: dihedrals
+
+
+         0 !NIMPHI: impropers
+
+
+         0 !NDON: donors
+
+
+         0 !NACC: acceptors
+
+
+         0 !NNB
+
+         0         0         0         0         0         0         0
+
+         1         0 !NGRP
+         0         0         0
+

colvartools/namd/Common/common.namd

@@ -0,0 +1,91 @@
+# -*- tcl -*-
+
+if { [info exists mol_name] == 0 } {
+    set mol_name "ammonia_water"
+}
+
+if { [info exists coor_pdb_file] == 0 } {
+    set coor_pdb_file "../Common/${mol_name}.pdb"
+}
+
+if { ${mol_name} == "ammonia_water" } { 
+    set par_files [list "../Common/par_all36_ammonia_water.prm"]
+}
+
+# Uncomment below to use full C36 force field (needs NAMD-compatible water_ions)
+# set ff_folder "../Common/charmm-ff"
+# set par_files [list \
+#                    "${ff_folder}/par_all36_prot.prm" \
+#                    "${ff_folder}/par_all36_cgenff.prm" \
+#                    "${ff_folder}/par_all36_carb.prm" \
+#                    "${ff_folder}/par_all36_lipid.prm" \
+#                    "${ff_folder}/par_all36_na.prm" \
+#                    "${ff_folder}/par_water_ions.prm"]
+
+if { [info exists job] == 0 } {
+    abort "Undefined variable job"
+}
+
+paraTypeCharmm          on
+structure               ../Common/${mol_name}.psf
+foreach par_file ${par_files} {
+    parameters          ${par_file}
+}
+1-4scaling              1.0
+exclude                 scaled1-4
+rigidBonds              all
+useSettle               on
+
+switchdist              10.0
+cutoff                  12.0
+pairlistdist            14.0
+vdwForceSwitching       on
+fullElectFrequency      1
+nonBondedFreq           1
+
+seed                    87654321
+
+langevin                on
+langevinTemp            300.0
+langevinDamping         10.0
+
+coordinates             ${coor_pdb_file}
+if { [info exists coor_bin_file] > 0 } {
+    binCoordinates      ${coor_bin_file}
+}
+if { [info exists vel_bin_file] > 0 } {
+    binVelocities       ${vel_bin_file}
+} else {
+    if { [info exists vel_pdb_file] > 0 } {
+        velocities      ${vel_pdb_file}
+    } else {
+        temperature     300.0
+    }
+}
+if { [info exists xsc_file] > 0 } {
+    extendedSystem      ${xsc_file}
+}
+
+COMmotion               no
+zeroMomentum            yes
+
+timestep                1.0
+
+outputName              ${job}
+restartName             ${job}.tmp
+restartFreq             ${restart_freq}
+if { [info exists pdb_restart] > 0 } {
+    binaryRestart       no
+    binaryOutput        no
+} else {
+    binaryRestart       yes
+    binaryOutput        yes
+}
+DCDFile                 ${job}.coor.dcd
+DCDFreq                 2500
+DCDUnitCell             yes
+
+outputEnergies          500
+outputMomenta           500
+outputPressure          500
+outputTiming            500

colvartools/namd/Common/par_all36_ammonia_water.prm

@@ -0,0 +1,36 @@
+! Trimmed-down force field file for CGENFF ammonia and TIP3P water
+! Original files at: https://mackerell.umaryland.edu/charmm_ff.shtml#charmm
+
+ATOMS
+!hydrogens
+MASS  -1  HGPAM3     1.00800 ! polar H, NEUTRAL ammonia (#)
+MASS  -1  HT         1.00800 ! TIPS3P WATER HYDROGEN
+MASS  -1  NG331     14.00700 ! neutral ammonia nitrogen
+MASS  -1  OT        15.99940 ! TIPS3P WATER OXYGEN
+
+BONDS
+HT    HT      0.0       1.5139  ! from TIPS3P geometry (for SHAKE w/PARAM)
+HT    OT    450.0       0.9572  ! from TIPS3P geometry
+NG331  HGPAM3  455.50     1.0140 ! AMINE aliphatic amines
+
+ANGLES
+HGPAM3 NG331  HGPAM3   41.50    107.10 ! AMINE aliphatic amines kevo: 29.00 -> 4
+HT   OT   HT     55.0      104.52   ! FROM TIPS3P GEOMETRY
+
+DIHEDRALS
+
+IMPROPERS
+
+NONBONDED nbxmod  5 atom cdiel fshift vatom vdistance vfswitch -
+cutnb 14.0 ctofnb 12.0 ctonnb 10.0 eps 1.0 e14fac 1.0 wmin 1.5
+
+!see mass list above for better description of atom types
+!hydrogens
+HGPAM3   0.0       -0.0120     0.8700 ! aliphatic amines
+NG331    0.0       -0.0700     1.9800 ! aliphatic amines
+HT       0.0       -0.046     0.2245
+OT       0.0       -0.1521    1.7682
+
+NBFIX
+
+END

colvartools/namd/README.md

@@ -0,0 +1,25 @@
+## Tools for replica-exchange computations in NAMD
+
+The Tcl script `colvars_replica_utils.tcl` can be sourced at any time during a NAMD job: inside it are defined Tcl functions that implement bias-exchange protocols; most of these functions rely on the Colvars scripting API.
+
+Unlike in the Tcl scripts in the `lib/replica` folder of a NAMD distribution tree, no `.history` files are written and un-shuffling the DCD trajectory files with `sortreplicas` is not yet supported.  This may change in the future.
+
+In theory, any types of bias can be exchanged, but so far only metadynamics biases have been tested.
+
+### Examples
+
+The folders below contain ready-to-run examples for each specific methodology.  Running them requires a multiple-replica capable NAMD build; for a single workstation/node, this is usually a "netlrts" build (NB: NOT the "multicore" build).
+
+Folder | Description
+------ | -----------
+`eq` | Standard MD run with multiple replicas; use this example to test that the NAMD build being used supports replica-exchange.
+`bemtd` | Example input for bias-exchange metadynamics; requires an updated version of the Colvars module inside it (version >= xxxx-xx-xx).
+
+
+### Post-processing
+
+The Python script `load_trajectories.py` serves the purpose of reading `.colvars.traj` files produced during bias-exchange simulations (i.e. where the variables are continuous but the biases change), and to write out trajectories where the identity of the active bias in each replica is constant.  This allows, for example to monitor the statistical correctness of each replica with the given bias (or lack thereof).
+
+Running this script in either the `eq` or `bemtd` folder should output the desired results.
+
+Internally, this script uses the `plot_colvars_traj.py` script from the parent directory to read Colvars trajectory files.

colvartools/namd/bemtd/README.md

@@ -0,0 +1,27 @@
+## Bias-exchange metadynamics example
+
+The `bemtd` folder contains an example input to run a bias-exchange metadynamics simulation using three collective variables.  These are the distances between the nitrogen atom of an ammonia molecule vs. the oxygen atom and the two hydrogen atoms of a water molecule.
+
+The script `run.sh` will run two consecutive jobs, about an hour long each.  To change this, feel free to adjust the environment variable `numsteps` or range of values for the variable `ijob`.
+
+### Input requirements
+
+The file `ammonia_water-bemtd.colvars.in` contains the Colvars configuration for all biases, including three metadynamics biases that can be exchanged, and one walls restraint that remains fixed.
+
+Besides the required keywords, each metadynamics bias should contain also the following:
+- `outputEnergy yes`, so that the biasing energy can be easily be read during post-processing of the `.colvars.traj` file under the column `E_<bias name>`;
+- `outputFeatures active`, so that the active/inactive status of each bias is reported as 1/0 values in the `<bias name>_active` column; *when the value of this column is 0 for a given bias, its reported energy is out of date and should be ignored even if positive*.
+
+The comments in the NAMD script `bemtd.namd` describe in more detail how to define the biases are subject to being exchanged.  The current implementation requires that *no more than one replica at a time is without biases*.
+
+
+### Analysis
+
+Each replica runs continuously in Cartesian space, with the biasing potential switching from time to time.  This is *different* from what e.g. GROMACS+PLUMED uses, where the coordinates of the replicas are exchanged and represent trajectorioes with a constant bias instead.
+
+To obtain the latter from the former, the `load_trajectories.py` script in the parent folder can be used to extract histograms and time series of selected variables with a constant bias, for example:
+```
+python3 ../load_trajectories.py --first 100000 --skip 10 --variables dist_N_O --biases neutral
+```
+
+Histogram files computed over trajectory segments without bias (e.g. `histogram.neutral-dist_N_O.dat`) can be compared directly with the equivalent files computed on simulations of the same systems without bias (example input in the `eq` folder).

colvartools/namd/bemtd/ammonia_water-bemtd.colvars.in

@@ -0,0 +1,101 @@
+colvar {
+    name dist_N_O
+    outputAppliedForce on
+    upperBoundary 12.0
+    width 0.02
+    distance {
+        group1 {
+            psfSegID AMM1
+            atomNameResidueRange N1 1-1
+        }
+        group2 {
+            psfSegID WAT
+            atomNameResidueRange OH2 1-1
+        }
+    }
+} 
+
+
+metadynamics {
+    name mtd_dist_N_O
+    colvars dist_N_O
+    gaussianSigmas 0.1
+    hillWeight 0.001
+    outputEnergy yes
+    outputFeatures active
+}
+
+
+
+colvar {
+    name dist_N_H1
+    outputAppliedForce on
+    upperBoundary 12.0
+    width 0.02
+    distance {
+        group1 {
+            psfSegID AMM1
+            atomNameResidueRange N1 1-1
+        }
+        group2 {
+            psfSegID WAT
+            atomNameResidueRange H1 1-1
+        }
+    }
+} 
+
+
+metadynamics {
+    name mtd_dist_N_H1
+    colvars dist_N_H1
+    gaussianSigmas 0.1
+    hillWeight 0.001
+    outputEnergy yes
+    outputFeatures active
+}
+
+
+
+colvar {
+    name dist_N_H2
+    outputAppliedForce on
+    upperBoundary 12.0
+    width 0.02
+    distance {
+        group1 {
+            psfSegID AMM1
+            atomNameResidueRange N1 1-1
+        }
+        group2 {
+            psfSegID WAT
+            atomNameResidueRange H2 1-1
+        }
+    }
+} 
+
+
+metadynamics {
+    name mtd_dist_N_H2
+    colvars dist_N_H2
+    gaussianSigmas 0.1
+    hillWeight 0.001
+    outputEnergy yes
+    outputFeatures active
+}
+
+
+
+harmonicWalls {
+    name wall_dist_N_O
+    colvars dist_N_O
+    forceConstant 0.6
+    upperWalls 8.0
+}
+
+
+neutral {
+    name neutral
+    colvars dist_N_O
+    outputEnergy yes
+    outputFeatures active
+}