ExTASY wrappers version of lsdmap

.gitignore

@@ -6,5 +6,3 @@ build/
 *.h
 *.cpp
 *.swp
-bin/*
-examples/*

README.md

@@ -4,15 +4,15 @@ LSDMap
 LSDMap package is used to compute Locally Scaled Diffusion Map. 
 Typical usage is to call the "lsdmap" script:
 
-	lsdmap -f <configuration_file> -c <structure_file> <other_options>
+	lsdmap -f <configuration_file> -c <structure_file> -t <topology_file> <other_options>
 
 or using MPI:
 
-	mpiexec -n <number_of_processors> lsdmap -f <configuration_file> -c <structure_file> <other_options>
+	mpiexec -n <number_of_processors> lsdmap -f <configuration_file> -c <structure_file> -t <topology_file> <other_options>
 
 A typical example of configuration file is ./examples/lsdmap/config.ini
 The structure file should contain all the configurations needed to compute
-LSDMap. After execution of the script, a .ev and a .eg file should have
+LSDMap (e.g., a multi-model .gro file). After execution of the script, a .ev and a .eg file should have
 been generated containing the eigenvectors and eigenvalues, respectively. 
 
 See the paper W. Zheng, M. A. Rohrdanz, M. Maggioni and C. Clementi, J. Chem. Phys., 2011, 134, 144109 for more information on how LSDMap works.
@@ -32,6 +32,16 @@ packages installed:
 
 * cython; version 0.20 or later
 
+If not already available, installing LSDMap will automatically install:
+
+* MDTraj
+
+* h5py
+
+* dask
+
+* cloudpickle
+
 Version 2.6.x or 2.7.x of python should be used. 
 
 
@@ -60,7 +70,7 @@ which contains 1000 configurations of alanine dipeptide in vacuum and
 an example of configuration file (.ini) that should be used to compute
 LSDMap. To test the program, simply type in this folder:
 
-	lsdmap -f config.ini -c aladip_1000.gro
+	lsdmap -f config.ini -c aladip_1000.gro -t aladip.gro
 
 After execution, a file ".ev" and a file ".eg" must have been generated.
 They contain the eigenvectors and eigenvalues of the Fokker-Planck operator,
@@ -73,7 +83,7 @@ and so on.
 
 LSDMap can be computed using MPI using a command similar to:
 
-	mpiexec -n <number_of_processors> lsdmap -f <configuration_file> -c <structure_file>
+	mpiexec -n <number_of_processors> lsdmap -f <configuration_file> -c <structure_file> -t <topology_file>
 
 For more information on lsdmap command, simply type:
 
@@ -88,15 +98,8 @@ DM-d-MD (Diffusion-Map directed Molecular Dynamics) is an adaptive sampling
 algorithm based on LSDMap. For an introduction to DM-d-MD, see the paper
 J.Preto and C. Clementi, Phys. Chem. Chem. Phys., 2014, 16, 19181-19191.
 
-Besides LSDMap, DM-d-MD requires GROMACS to be correctly installed.
-
-DM-d-MD is automatically installed when installing LSDMap via the command:
-
-        python setup.py install
-
-A typical usage of DM-d-MD is to call:
-
-	dmdmd -f <configuration_file>
+DM-d-MD jobs are run through python scripts. For examples, see the examples/dmdmd_ala12 and examples/dmdmd directories.
+Other examples, using AMBER as the MD engine as opposed to GROMACS, are available in the examples folder that comes with Extasy wrappers (https://bitbucket.org/claughton/wrappers)
 
 
 Prerequisites
@@ -105,19 +108,24 @@ Prerequisites
 In order to use DM-d-MD, it is required that GROMACS has been correctly
 installed and that "grompp" and "mdrun" commands are working properly
 for a serial utilization. If not, please visit http://www.gromacs.org/.
+The ExTASY wrappers python library is also required (https://bitbucket.org/claughton/wrappers)
 
 
 Testing
 -------
 
-The Folder examples/dmdmd contains an example of DM-d-MD configuration
-file (dmdmd.ini) as well as files required to run GROMACS MD simulations
+The Folder examples/dmdmd contains an example of DM-d-MD python script (rundmdmd.py) as well as files required to run GROMACS MD simulations
 for the photoactive yellow protein (PYP). DM-d-MD can be launched by 
 executing the command:
 
-	dmdmd -f dmdmd.ini
+	python rundmdmd.py
 
 within the specified folder.
+Note: it may be wise to issue the command:
+
+        setenv GMX_MAXBACKUP -1
+
+first.
 
 
 Diffusion-Map Sampling

dmdmd/reweighting.py

@@ -18,11 +18,11 @@ class ReweightingStep(object):
     def initialize(self, args):
 
         # read structure file
-        struct_file = reader.open(args.struct_file, velocities=True)
-        self.struct_filename = struct_file.filename
-        self.coords = struct_file.readlines()
+        with open(args.struct_file[0]) as f:
+
+            self.coords = f.readlines()
 
-        self.npoints = self.coords.shape[0]
+        self.npoints = len(self.coords)
 
         # read number of copies
         ncfile = reader.open(args.ncfile)
@@ -118,9 +118,10 @@ def get_nneighbors(self, nnfile_name, line):
     def save(self, args):
 
         # save new coordinates
-        format_output_file = os.path.splitext(args.output_file[0])[1]
-        struct_file_writer = writer.open(format_output_file, pattern=args.struct_file[0])
-        struct_file_writer.write(self.new_coords, args.output_file[0])
+        struct_file_writer = open(args.output_file[0], 'w')
+        for s in self.new_coords:
+            struct_file_writer.write(s)
+        struct_file_writer.close()
 
         # save wfile
         wfile_writer = writer.open('.w')
@@ -202,6 +203,8 @@ def run(self):
 
         self.new_coords = np.array(new_coords)
         self.new_weights = np.array(new_weights)
+
+        print "number of simulations in next step:", self.new_weights.shape[0]
 
         sum_old_weights=int(round(np.sum(self.weights)))
         sum_new_weights=int(round(np.sum(self.new_weights)))

dmdmd/tools/p_mdrun.py

@@ -9,7 +9,7 @@
 
 class ParallelMDruns(object):
     ##TICA
-    def write_script(self, script, directory, rank, mdpfile, grofile, topfile, output_grofile, start_coord, num_coord, tprfile='topol.tpr', trrfile='traj.trr', edrfile='ener.edr', shebang='/bin/bash', ndxfile='', grompp_options='', mdrun_options=''):
+    def write_script(self, script, directory, rank, mdpfile, grofile, topfile, output_grofile, start_coord, num_coord, tprfile='topol.tpr', trrfile='traj.trr', edrfile='ener.edr', shebang='/bin/bash', ndxfile='', grompp_options='', mdrun_options='', gmx_suffix=''):
         ##TICA
         frame_designation_path = directory + '/' + 'frame_desig.txt'
         with open(frame_designation_path, 'w') as fiile:
@@ -48,8 +48,8 @@ def write_script(self, script, directory, rank, mdpfile, grofile, topfile, outpu
   sed "$start"','"$end"'!d' $startgro > $tmpstartgro
 
   # gromacs preprocessing & MD
-  grompp %(grompp_options)s -f %(mdpfile)s -c $tmpstartgro -p %(topfile)s %(ndxfile_option)s -o %(tprfile)s 1>/dev/null 2>/dev/null
-  mdrun -nt 1 %(mdrun_options)s -s %(tprfile)s -o %(trrfile)s -e %(edrfile)s 1>/dev/null 2>/dev/null
+  grompp%(gmx_suffix)s %(grompp_options)s -f %(mdpfile)s -c $tmpstartgro -p %(topfile)s %(ndxfile_option)s -o %(tprfile)s 1>/dev/null 2>/dev/null
+  mdrun%(gmx_suffix)s -nt 1 %(mdrun_options)s -s %(tprfile)s -o %(trrfile)s -e %(edrfile)s 1>/dev/null 2>/dev/null
 
   # store data
   cat confout.gro >> $outgro
@@ -128,6 +128,10 @@ def create_arg_parser(self):
             type=str,
             dest="tmpdir")
 
+        parser.add_argument("--gmx_suffix",
+            type=str,
+            default="")
+
         return parser
 
     def copy_additional_files(self, afiles, rundir):
@@ -299,7 +303,7 @@ def run(self):
 
         script = 'run.sh'
         self.write_script(script, rundir, rank, 'grompp.mdp', 'start.gro', 'topol.top', 'out.gro', thread_start_coord, thread_num_coord,\
-                     ndxfile=args.ndxfile, grompp_options=args.grompp_options, mdrun_options=args.mdrun_options) ##TICA
+                     ndxfile=args.ndxfile, grompp_options=args.grompp_options, mdrun_options=args.mdrun_options, gmx_suffix=args.gmx_suffix) ##TICA
 
         comm.Barrier()
 
@@ -326,6 +330,6 @@ def run(self):
                             print >> output_file, line.replace("\n", "")
 
             print "Output data have been saved in %s" %args.output_file
-
+        return
 if __name__ == '__main__':
     ParallelMDruns().run()

examples/dmdmd/README

@@ -0,0 +1,7 @@
+Example DMDMD run on PYP. Two cycles, 50 replicas per cycle.
+
+To run:
+
+setenv GMX_MAXBACKUP -1
+python rundmdmd.py
+

examples/dmdmd/ffG43a1.itp

@@ -1,9 +1,9 @@
 #define _FF_GROMOS96
 #define _FF_GROMOS43A1
 
-[ defaults ]
+;[ defaults ]
 ; nbfunc	comb-rule	gen-pairs	fudgeLJ	fudgeQQ
-  1		1		no		1.0	1.0
+;  1		1		no		1.0	1.0
 
 #include "ffG43a1nb.itp"
 #include "ffG43a1bon.itp"