adjusted uninfluential energy tail corrections

chignolin/Analyze_folded.py

@@ -1,20 +1,20 @@
 #! /usr/bin/env python3
 
-# Used for Fig.3
+# Default arguments used for Fig.3
 # Calculates the folded fraction as a function of temperature and pressure
+# CAUTION: run Prepare_analysis.sh  before this
 
 import sys
 import numpy as np
 import pandas as pd
 import subprocess
 import argparse
 
-# CAUTION: run label_and_combine_traj.sh  before this, to combine the trajectories
 #set columns
 ene_col=1
 vol_col=2
 bias_col=4
-basin_col=6
+basin_col=5
 
 #parser
 parser = argparse.ArgumentParser(description='calculate fraction folded and DeltaG over a range of temperatures and pressures')
@@ -25,7 +25,7 @@
 parser.add_argument('--minpres',dest='minpres',type=float,default=1,required=False,help='the minimum preserature')
 parser.add_argument('--maxpres',dest='maxpres',type=float,default=4000,required=False,help='the maximum preserature')
 parser.add_argument('--nbins',dest='nbins',type=int,default=50,required=False,help='number of bins')
-parser.add_argument('--tran',dest='tran',type=int,default=0,required=False,help='transient to be skipped')
+parser.add_argument('--tran',dest='tran',type=int,default=400000,required=False,help='transient to be skipped')
 parser.add_argument('--bck',dest='bck',type=str,default='',required=False,help='backup prefix, e.g. \"bck.0.\"')
 parser.add_argument('-f',dest='filename',type=str,default='all_Colvar.data',required=False,help='input file name')
 parser.add_argument('-o',dest='outfilename',type=str,default='fraction_folded.data',required=False,help='output file name')

chignolin/Analyze_histo-2D.py

@@ -1,20 +1,19 @@
 #! /usr/bin/env python3
 
-# Used for Fig.2a
+# Default arguments used for Fig.2a
 # Calculates the histogram of the sampled target distribution as a function of energy and volume
+# CAUTION: run Prepare_analysis.sh  before this
 
 import sys
 import numpy as np
 import pandas as pd
 import subprocess
 import argparse
 
-# CAUTION: run label_and_combine_traj.sh  before this, to combine the trajectories
-
 #parser
 parser = argparse.ArgumentParser(description='calculate Histogram of energies and volumes')
 parser.add_argument('--nbins',dest='nbins',type=int,default=300,required=False,help='number of bins')
-parser.add_argument('--tran',dest='tran',type=int,default=0,required=False,help='transient to be skipped')
+parser.add_argument('--tran',dest='tran',type=int,default=400000,required=False,help='transient to be skipped')
 parser.add_argument('--bck',dest='bck',type=str,default='',required=False,help='backup prefix, e.g. \"bck.0.\"')
 parser.add_argument('-f',dest='filename',type=str,default='all_Colvar.data',required=False,help='input file name')
 parser.add_argument('-o',dest='outfilename',type=str,default='Histo-2D.data',required=False,help='output file name')

chignolin/Analyze_neff-2D.py

@@ -1,16 +1,15 @@
 #! /usr/bin/env python3
 
-# Used for Fig.2b
+# Default arguments used for Fig.2b
 # Calculates the effective sample size as a function of temperature and pressure
+# CAUTION: run Prepare_analysis.sh  before this
 
 import sys
 import numpy as np
 import pandas as pd
 import subprocess
 import argparse
 
-# CAUTION: run label_and_combine_traj.sh  before this, to combine the trajectories
-
 #parser
 parser = argparse.ArgumentParser(description='calculate Neff over a range of temperatures and pressures')
 parser.add_argument('--temp',dest='temp',type=float,default=500,required=False,help='the simulation temperature')
@@ -20,7 +19,7 @@
 parser.add_argument('--minpres',dest='minpres',type=float,default=1,required=False,help='the minimum preserature')
 parser.add_argument('--maxpres',dest='maxpres',type=float,default=4000,required=False,help='the maximum preserature')
 parser.add_argument('--nbins',dest='nbins',type=int,default=50,required=False,help='number of bins')
-parser.add_argument('--tran',dest='tran',type=int,default=0,required=False,help='transient to be skipped')
+parser.add_argument('--tran',dest='tran',type=int,default=400000,required=False,help='transient to be skipped')
 parser.add_argument('--bck',dest='bck',type=str,default='',required=False,help='backup prefix, e.g. \"bck.0.\"')
 parser.add_argument('-f',dest='filename',type=str,default='all_Colvar.data',required=False,help='input file name')
 parser.add_argument('-o',dest='outfilename',type=str,default='Neff-2D.data',required=False,help='output file name')

chignolin/Analyze_folded-blocks.py → chignolin/Analyze_temp_folding.py

@@ -1,32 +1,32 @@
 #! /usr/bin/env python3
 
-# Used for the inset of Fig.3
+# Default arguments used for the inset of Fig.3
 # Calculates the folded fraction at fixed pressure for different temperatures, estimating the uncertainties with block average
+# CAUTION: run Prepare_analysis.sh  before this
 
 import sys
 import numpy as np
 import pandas as pd
 import subprocess
 import argparse
 
-# CAUTION: run label_and_combine_traj.sh  before this, to combine the trajectories
 
 #set columns
 ene_col=1
 vol_col=2
 bias_col=4
-basin_col=6
+basin_col=5
 
 #parser
 parser = argparse.ArgumentParser(description='reweight as a function of a CV for a given temperature and pressure')
-parser.add_argument('--blocks',dest='num_blocks',type=int,required=True,help='number of blocks')
+parser.add_argument('--blocks',dest='num_blocks',type=int,default=4,required=False,help='number of blocks')
 parser.add_argument('--temp',dest='temp',type=float,default=500,required=False,help='the simulation temperature')
 parser.add_argument('--mintemp',dest='mintemp',type=float,default=280,required=False,help='the minimum temperature')
 parser.add_argument('--maxtemp',dest='maxtemp',type=float,default=370,required=False,help='the maximum temperature')
 parser.add_argument('--pres',dest='pres',type=float,default=2000,required=False,help='the simulation pressure (bar)')
 parser.add_argument('--rewpres',dest='rewpres',type=float,default=1,required=False,help='the reweighting pressure (bar)')
 parser.add_argument('--nbins',dest='nbins',type=int,default=50,required=False,help='number of bins')
-parser.add_argument('--tran',dest='tran',type=int,default=0,required=False,help='transient to be skipped')
+parser.add_argument('--tran',dest='tran',type=int,default=400000,required=False,help='transient to be skipped')
 parser.add_argument('--bck',dest='bck',type=str,default='',required=False,help='backup prefix, e.g. \"bck.0.\"')
 parser.add_argument('-f',dest='filename',type=str,default='all_Colvar.data',required=False,help='input file name')
 parser.add_argument('-o',dest='outfilename',type=str,default='temp_folded.data',required=False,help='output file name')

chignolin/Prepare_analysis.sh

@@ -0,0 +1,60 @@
+#! /bin/bash
+
+# This script prepares the file all_Colvar.data, to be used for further analysis:
+# - the full energy is taken from GROMACS edr files (this might be fixed in the future, see https://github.com/plumed/plumed2/issues/567 )
+#   NB: GROMACS energy is very similar to the PLUMED one, and it actually does not make a difference for our observables
+# - each point of the trajectory is assigned to a basin by looking at the C_alpha-RMSD using the same criteria as https://dx.plos.org/10.1371/journal.pone.0032131
+# - all walkers trajectories are combined and sorted according time
+
+end=300000
+if [ $# -eq 1 ]
+then
+  end=$1
+fi
+
+#same criteria as https://dx.plos.org/10.1371/journal.pone.0032131
+low_lim=0.1
+up_lim=0.4
+
+echo "  end = $end"
+bck.meup.sh -i all_Colvar.data
+echo "#! FIELDS time full_ene vol pdb_rmsd opes.bias basin" > all_Colvar.data
+for i in `seq 0 39`
+do 
+#get full energy, with tail corrections
+  ene=energy.$i.xvg
+  gmx_mpi energy -f chignolin.$i.edr -o $ene <<< "12 23 0" #the volume is actually the same...
+  awk -v e=$end 'NR>26{if($1==e+1) exit;print $0}' $ene > tmp.$ene #caution, header length might change
+#add basin labels
+  col=Colvar.${i}.data
+  echo -en "  $f\r"
+  if [ $i -eq 6 ] || [ $i -eq 8 ] || [ $i -eq 9 ] #specific of my initial conditions
+  then 
+    st=0
+  else
+    st=1
+  fi
+  awk -v a=$st -v low=$low_lim -v up=$up_lim -v e=$end '{
+    if($1==e+1) exit;
+    if($1=="#!") skip=1;
+    else{
+      if(skip==1) skip=0;
+      else{
+        if(a==0 && $4>up) a=1;
+        if(a==1 && $4<low) a=0;
+        print $4,$5,a; #rmsd bias basin
+      }
+    }
+  }' $col > tmp.$col
+#put all together
+  if [ `wc -l < tmp.$ene` != `wc -l < tmp.$col` ]
+  then
+    echo "something is wrong, energy file and colvar file have different lengths"
+    exit
+  fi
+  paste tmp.$ene tmp.$col |awk '{print $0}' >> all_Colvar.data
+  rm tmp.$ene tmp.$col
+done
+echo -en "  Sorting...        \r"
+sort -gsk1 all_Colvar.data > tmp.all_Colvar.data
+mv tmp.all_Colvar.data all_Colvar.data