Adding support for PLUMED: based on xiangda's efforts

.gitignore

@@ -22,4 +22,10 @@ CMakeFiles/
 .ycm_extra_conf.py
 compile_commands.json
 source.sh
+src/Plumed.cmake
+src/Plumed.h
+src/Plumed.inc
+example/*/inputs
+example/*/outputs
+example/*/results
 

README.md

@@ -156,7 +156,7 @@ the trajectory will be discontinuous due to replica swapping.
 ### Advanced usage
 
 Explore the `upside/example/` folder for examples on running with restraints, membrane
-simulations, pulling simulations, HDX prediction, and more.
+simulations, pulling simulations, HDX prediction, PLUMED usaged and more.
 
 ## Simulation analysis and visualization
 

example/14.Plumed/0.run.py

@@ -0,0 +1,137 @@
+#!/usr/bin/env python
+import sys, os, shutil
+import subprocess as sp
+import numpy as np
+
+upside_path = os.environ['UPSIDE_HOME']
+upside_utils_dir = os.path.expanduser(upside_path+"/py")
+sys.path.insert(0, upside_utils_dir)
+import run_upside as ru
+
+#----------------------------------------------------------------------
+## General Settings and Path
+#----------------------------------------------------------------------
+
+pdb_id         = 'chig' # switch to 1dfn for multi-chain showing
+pdb_dir        = './pdb'
+plumed_dir     = 'plumed_input'
+plumed_file    = 'rmsd_metad.in'  # change this to include different plumed input files
+PLUMED_JUST_PRINT = False
+sim_id         = 'metad'
+is_native      = True
+ff             = 'ff_2.1'
+T              = 0.8
+duration       = 10000
+frame_interval = 100
+base_dir       = './'
+
+#----------------------------------------------------------------------
+## Initialization
+#----------------------------------------------------------------------
+
+input_dir  = "{}/inputs".format(base_dir)
+output_dir = "{}/outputs".format(base_dir)
+run_dir    = "{}/{}".format(output_dir, sim_id)
+
+make_dirs = [input_dir, output_dir, run_dir]
+for direc in make_dirs:
+    if not os.path.exists(direc):
+        os.makedirs(direc)
+
+#----------------------------------------------------------------------
+## Generate Upside readable initial structure (and fasta) from PDB 
+#----------------------------------------------------------------------
+
+print ("Initial structure gen...")
+cmd = (
+       "python {0}/PDB_to_initial_structure.py "
+       "{1}/{2}.pdb "
+       "{3}/{2} "
+       "--record-chain-breaks "
+      ).format(upside_utils_dir, pdb_dir, pdb_id, input_dir )
+print (cmd)
+sp.check_output(cmd.split())
+
+
+#----------------------------------------------------------------------
+## Configure
+#----------------------------------------------------------------------
+
+# parameters
+param_dir_base = os.path.expanduser(upside_path+"/parameters/")
+param_dir_common = param_dir_base + "common/"
+param_dir_ff = param_dir_base + '{}/'.format(ff)
+
+# options
+print ("Configuring...")
+fasta = "{}/{}.fasta".format(input_dir, pdb_id)
+kwargs = dict(
+               rama_library              = param_dir_common + "rama.dat",
+               rama_sheet_mix_energy     = param_dir_ff + "sheet",
+               reference_state_rama      = param_dir_common + "rama_reference.pkl",
+               hbond_energy              = param_dir_ff + "hbond.h5",
+               rotamer_placement         = param_dir_ff + "sidechain.h5",
+               dynamic_rotamer_1body     = True,
+               rotamer_interaction       = param_dir_ff + "sidechain.h5",
+               environment_potential     = param_dir_ff + "environment.h5",
+               bb_environment_potential  = param_dir_ff + "bb_env.dat",
+               chain_break_from_file     = "{}/{}.chain_breaks".format(input_dir, pdb_id),
+             )
+
+if is_native:
+    kwargs['initial_structure'] =  "{}/{}.initial.npy".format(input_dir, pdb_id)
+
+config_base = "{}/{}.up".format( input_dir, pdb_id)
+config_stdout = ru.upside_config(fasta, config_base, **kwargs)
+
+print ("Config commandline options:")
+print (config_stdout)
+
+
+#----------------------------------------------------------------------
+## Advanced Configure
+#----------------------------------------------------------------------
+
+# Here, we use the run_upside.advanced_config function to add more advanced configuration
+# on the config file obtained in the previous step. advanced_config() allows us to add many
+# features, including but not limited to:
+#     pulling energy
+#     restraints (spring energy)
+#     wall potential
+#     contact energy
+# Example 8 will show more details.
+
+kwargs = dict(
+               #restraint_groups = ['0-9'],
+               #restraint_spring_constant = 2.0,
+                plumed = f'{plumed_dir}/{plumed_file}',
+                plumed_just_print = PLUMED_JUST_PRINT,  # True: no bias, just print
+             )
+
+config_stdout = ru.advanced_config(config_base, **kwargs)
+print ("Advanced Config commandline options:")
+print (config_stdout)
+
+#----------------------------------------------------------------------
+## Run Settings
+#----------------------------------------------------------------------
+
+
+randomseed = np.random.randint(0,100000) # Might want to change the fixed seed for the random number
+randomseed = 1
+
+upside_opts = (
+                 "--duration {} "
+                 "--frame-interval {} "
+                 "--temperature {} "
+                 "--seed {} "
+              )
+upside_opts = upside_opts.format(duration, frame_interval, T, randomseed)
+
+h5_file  = "{}/{}.run.up".format(run_dir, pdb_id)
+log_file = "{}/{}.run.log".format(run_dir, pdb_id)
+shutil.copyfile(config_base, h5_file)
+
+print ("Running...")
+cmd = "{}/obj/upside {} {} | tee {}".format(upside_path, upside_opts, h5_file, log_file)
+sp.check_call(cmd, shell=True)

example/14.Plumed/compare_plumed_upside.py

@@ -0,0 +1,103 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+
+
+
+import numpy as np
+import mdtraj_upside as mu
+import mdtraj as md
+import os
+import matplotlib.pyplot as plt
+
+POS = False
+RG = False
+RMSD = True
+
+# ## compare coordinates
+
+
+
+
+#traj = mu.load_upside_traj("outputs/simple_test/chig.run.up")
+traj = mu.load_upside_traj("outputs/metad/chig.run.up")
+
+
+
+
+
+def quick_xyz_parser(file, natoms):
+    coords_all = []
+    i=0
+    with open(file, 'r') as f:
+        for line in f.readlines():
+            if line.rstrip("\n") == str(natoms):
+                i+=1
+                try:
+                    coords
+                except NameError:
+                    coords = []
+                else:
+                    coords_all.append(coords)
+                    coords = []
+            elif line[0] != 'X':
+                pass
+            else:
+                X, Y, Z = line.lstrip("X ").rstrip("\n").split(" ")
+
+                coords.append([float(X), float(Y), float(Z)])
+        coords_all.append(coords)
+    return coords_all
+
+
+if not os.path.exists("results"):
+    os.mkdir("results")
+
+
+
+if POS:
+    coords_all = np.array(quick_xyz_parser('outputs/simple_test/chig.xyz.plumed', natoms=30))
+    core_ids = traj.topology.select("name N or name CA or name C")
+
+    plt.plot(traj.time, np.sum((coords_all - traj.xyz[:, core_ids])**2, axis=(1,2)))
+    plt.xlabel("Time (in Upside unit)")
+    plt.ylabel("RMS difference in coordinates")
+    plt.savefig("results/compare_pos_upside_plumed.png", dpi=200, bbox_inches='tight')
+    plt.show()
+
+
+# ## compare Rg
+
+def calc_rg(xyz):
+    return np.sqrt(np.sum((xyz - xyz.mean(axis=(0,1)))**2)/ len(xyz[0, :]))
+
+if RG:
+    rgs_upside = [calc_rg(t.xyz[:, core_ids]) for t in traj]
+    rgs_plumed = np.loadtxt("outputs/simple_test/chig.rg.plumed")[:, 1]
+    plt.plot(traj.time, rgs_upside, '-', label='mdtraj calculation')
+    plt.plot(traj.time, rgs_plumed, '--', label='plumed calculation')
+    plt.legend()
+    plt.xlabel("Time (in Upside unit)")
+    plt.ylabel("Rg (nm)")
+    plt.savefig("results/compare_Rg_upside_plumed.png", dpi=200, bbox_inches='tight')
+    plt.show()
+
+
+# ## Compare RMSD 
+if RMSD:
+    sele = traj.top.select("name CA")
+    rmsd_upside = md.rmsd(traj, traj, atom_indices=sele)
+    rmsd_plumed = np.loadtxt("outputs/metad/chig.rmsd.plumed")[:,1]
+    plt.plot(rmsd_upside, '-', label='mdtraj calculation')
+    plt.plot(rmsd_plumed, '--', label='plumed calculation')
+    plt.legend()
+    plt.xlabel("Time (in Upside unit)")
+    plt.ylabel("Rmsd (nm)")
+    plt.savefig("results/compare_Rmsd_upside_plumed.png", dpi=200, bbox_inches='tight')
+    plt.show()
+
+
+
+
+
+