ctest for thermostats: Langevin, Nose-Hoover and Berendsen

with plotting of temperature distribution using python matplotlib

diff --git a/ctest/1L2Y_B.inp b/ctest/1L2Y_B.inp
new file mode 100644 (file)
index 0000000..f0ce67a
--- /dev/null
+++ b/ctest/1L2Y_B.inp
@@ -0,0 +1,15 @@
+1L2Y with Berendsen thermostat in ff_1l2y MD simulation
+SEED=-3059743 PDBREF MD EXTCONF RESCALE_MODE=2 RESPA
+nstep=200000  ntwe=10    ntwx=10000    dt=0.20 damax=10.0 lang=0 tbf           &
+tau_bath=1.0 t_bath=300 reset_vel=0     respa ntime_split=1 maxtime_split=512
+WLONG=1.00000 WSCP=1.23315 WELEC=0.84476 WBOND=1.00000 WANG=0.62954            &
+WSCLOC=0.10554 WTOR=1.84316 WTORD=1.26571 WCORRH=0.19212 WCORR5=0.00000        &
+WCORR6=0.00000 WEL_LOC=0.37357 WTURN3=1.40323 WTURN4=0.64673 WTURN6=0.00000    &
+WVDWPP=0.23173 WHPB=1.00000 WSCCOR=0.0                                         &
+CUTOFF=7.00000 WCORR4=0.00000
+1L2Y.pdb
+22
+ D   ASN LEU TYR ILE GLN TRP LEU LYS ASP GLY GLY PRO SER SER GLY ARG PRO PRO PRO
+ SER D  
+ 0
+ 0

diff --git a/ctest/1L2Y_L.inp b/ctest/1L2Y_L.inp
new file mode 100644 (file)
index 0000000..93f634f
--- /dev/null
+++ b/ctest/1L2Y_L.inp
@@ -0,0 +1,15 @@
+1L2Y with Langevin thermostat in ff_1l2y MD simulation
+SEED=-3059743 PDBREF MD EXTCONF RESCALE_MODE=2 RESPA
+nstep=200000 ntwe=10  ntwx=10000 dt=0.20 damax=10.0 lang=1                     &
+t_bath=300 reset_vel=0 respa ntime_split=1 maxtime_split=512
+WLONG=1.00000 WSCP=1.23315 WELEC=0.84476 WBOND=1.00000 WANG=0.62954            &
+WSCLOC=0.10554 WTOR=1.84316 WTORD=1.26571 WCORRH=0.19212 WCORR5=0.00000        &
+WCORR6=0.00000 WEL_LOC=0.37357 WTURN3=1.40323 WTURN4=0.64673 WTURN6=0.00000    &
+WVDWPP=0.23173 WHPB=1.00000 WSCCOR=0.0                                         &
+CUTOFF=7.00000 WCORR4=0.00000
+1L2Y.pdb
+22
+ D   ASN LEU TYR ILE GLN TRP LEU LYS ASP GLY GLY PRO SER SER GLY ARG PRO PRO PRO
+ SER D  
+ 0
+ 0

diff --git a/ctest/1L2Y_NH.inp b/ctest/1L2Y_NH.inp
new file mode 100644 (file)
index 0000000..3a1d388
--- /dev/null
+++ b/ctest/1L2Y_NH.inp
@@ -0,0 +1,15 @@
+1L2Y with Nose-Hoover thermostat in ff_1l2y
+SEED=-3059743 PDBREF MD EXTCONF RESCALE_MODE=2 RESPA
+nstep=200000 ntwe=10 ntwx=10000 dt=0.10 damax=10.0 lang=0                      &
+t_bath=300 reset_vel=0 NOSEHOOVER96   Q_NP=1.0
+WLONG=1.00000 WSCP=1.23315 WELEC=0.84476 WBOND=1.00000 WANG=0.62954            &
+WSCLOC=0.10554 WTOR=1.84316 WTORD=1.26571 WCORRH=0.19212 WCORR5=0.00000        &
+WCORR6=0.00000 WEL_LOC=0.37357 WTURN3=1.40323 WTURN4=0.64673 WTURN6=0.00000    &
+WVDWPP=0.23173 WHPB=1.00000 WSCCOR=0.0                                         &
+CUTOFF=7.00000 WCORR4=0.00000
+1L2Y.pdb
+22
+ D   ASN LEU TYR ILE GLN TRP LEU LYS ASP GLY GLY PRO SER SER GLY ARG PRO PRO PRO
+ SER D  
+ 0
+ 0

diff --git a/ctest/matplotlib_fit_hist.py b/ctest/matplotlib_fit_hist.py
new file mode 100755 (executable)
index 0000000..8601325
--- /dev/null
+++ b/ctest/matplotlib_fit_hist.py
@@ -0,0 +1,47 @@
+#! /usr/bin/env python
+
+import matplotlib
+#matplotlib.use('GTK')
+matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.optimize import curve_fit
+from math import sqrt
+import sys
+
+def prob_T(x,a):
+    g=np.float128(105.)
+    Tr=np.float128(300.)
+    return a * (  x**((g-2)/2) * np.exp( -g*x/(2*Tr) ) )   
+
+#x,y= np.loadtxt('1L2Y_L_GB000.stat',usecols=(0,10),skiprows=30,unpack=True)
+#x,y= np.loadtxt('1L2Y_NH_GB000.stat',usecols=(0,11),skiprows=10000,unpack=True)
+#x,y= np.loadtxt('1L2Y_B_GB000.stat',usecols=(0,10),skiprows=30,unpack=True)
+if (sys.argv[1] == '1L2Y_NH_GB000.stat'):
+ x,y= np.loadtxt(sys.argv[1],usecols=(0,11),skiprows=10000,unpack=True)
+else: 
+ x,y= np.loadtxt(sys.argv[1],usecols=(0,10),skiprows=10000,unpack=True)
+
+h,bin=np.histogram(y,bins=50,density=True)
+
+plt.bar(bin[:-1], h, width = bin[2]-bin[1])
+plt.xlim(min(bin), max(bin))
+plt.ylim(0,max(h))
+plt.ylabel('probality')
+plt.xlabel('temperature')
+
+center = (bin[:-1] + bin[1:]) / 2
+#print bin
+#print center
+popt, pcov = curve_fit(prob_T, center, h, p0=10E-107)
+xfine = np.linspace(min(bin), max(bin), 100)
+#print popt
+
+chi_squared = np.sum((prob_T(center, *popt)-h)**2)
+print '%15.10f' % chi_squared
+
+#print xfine
+#print prob_T(xfine,popt[0])
+plt.plot(xfine,prob_T(xfine,popt[0]),'-',c='red')
+plt.savefig(sys.argv[1]+'.png')
+#plt.show()   

diff --git a/ctest/prota_unres_energy_check.sh b/ctest/prota_unres_energy_check.sh
index 45759b9..3aa859d 100755 (executable)
--- a/ctest/prota_unres_energy_check.sh
+++ b/ctest/prota_unres_energy_check.sh
@@ -70,6 +70,29 @@ elif [ "$1" == "1l2y_micro" ]; then
  else
   exit 0
  fi
+elif [ "$1" == "1L2Y_L" ] || [ "$1" == "1L2Y_NH" ]; then
+ chi2=`./matplotlib_fit_hist.py $file_stat`
+ echo 'Chi2 for fitting theoretical temperature distribution ' ${chi2}
+ echo  "<DartMeasurementFile name=\"Temperature distribution $1\" type=\"image/png\">`pwd`/${file_stat}.png</DartMeasurementFile>"
+ 
+ if [ `echo "a=${chi2};a>0.001"|bc -l` != "0" ]; then
+    echo 'chi2 greater than 0.001'
+    exit 1
+ else
+    exit 0
+ fi
+elif [ "$1" == "1L2Y_B" ]; then
+ chi2=`./matplotlib_fit_hist.py $file_stat`
+ echo 'Chi2 for fitting theoretical temperature distribution ' ${chi2}
+ echo  "<DartMeasurementFile name=\"Temperature distribution $1\" type=\"image/png\">`pwd`/${file_stat}.png</DartMeasurementFile>"
+ 
+ if [ `echo "a=${chi2};a>0.01"|bc -l` != "0" ]; then
+    echo 'chi2 greater than 0.01'
+    exit 1
+ else
+    exit 0
+ fi
+
 else
  exit 1
 fi

diff --git a/source/unres/src_MD/CMakeLists.txt b/source/unres/src_MD/CMakeLists.txt
index ce3b264..c90d228 100644 (file)
--- a/source/unres/src_MD/CMakeLists.txt
+++ b/source/unres/src_MD/CMakeLists.txt
@@ -400,6 +400,15 @@ FILE(COPY ${CMAKE_SOURCE_DIR}/ctest/1l2y_MIN_INT.inp
 FILE(COPY ${CMAKE_SOURCE_DIR}/ctest/1l2y_MIN_REGULAR_INT.inp
         DESTINATION ${CMAKE_CURRENT_BINARY_DIR} )
 
+FILE(COPY ${CMAKE_SOURCE_DIR}/ctest/1L2Y_B.inp
+        DESTINATION ${CMAKE_CURRENT_BINARY_DIR} )
+
+FILE(COPY ${CMAKE_SOURCE_DIR}/ctest/1L2Y_L.inp
+        DESTINATION ${CMAKE_CURRENT_BINARY_DIR} )
+
+FILE(COPY ${CMAKE_SOURCE_DIR}/ctest/1L2Y_NH.inp
+        DESTINATION ${CMAKE_CURRENT_BINARY_DIR} )
+
 
 FILE(COPY ${CMAKE_SOURCE_DIR}/ctest/1L2Y.pdb
         DESTINATION ${CMAKE_CURRENT_BINARY_DIR} )
@@ -410,6 +419,12 @@ FILE(COPY ${CMAKE_SOURCE_DIR}/ctest/prota_unres_energy_check.sh
         FILE_PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE
 )
 
+FILE(COPY ${CMAKE_SOURCE_DIR}/ctest/matplotlib_fit_hist.py
+        DESTINATION ${CMAKE_CURRENT_BINARY_DIR} 
+        FILE_PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE
+)
+
+
 #=========================================
 #  test_prota_E0LL2Y.sh
 #=========================================
@@ -467,6 +482,9 @@ if(NOT UNRES_WITH_MPI)
     add_test(NAME UNRES_MIN_INT COMMAND sh ${CMAKE_CURRENT_BINARY_DIR}/test_prota_E0LL2Y.sh 1l2y_MIN_INT 1 )
     add_test(NAME UNRES_REGULAR COMMAND sh ${CMAKE_CURRENT_BINARY_DIR}/test_prota_E0LL2Y.sh 1l2y_MIN_REGULAR_INT 1 )
     add_test(NAME UNRES_MD_microcanonical COMMAND sh ${CMAKE_CURRENT_BINARY_DIR}/test_prota_E0LL2Y.sh 1l2y_micro 1 )
+    add_test(NAME UNRES_Langevin COMMAND sh ${CMAKE_CURRENT_BINARY_DIR}/test_prota_E0LL2Y.sh 1L2Y_L 1 )
+    add_test(NAME UNRES_NoseHoover COMMAND sh ${CMAKE_CURRENT_BINARY_DIR}/test_prota_E0LL2Y.sh 1L2Y_NH 1 )
+    add_test(NAME UNRES_Berendsen COMMAND sh ${CMAKE_CURRENT_BINARY_DIR}/test_prota_E0LL2Y.sh 1L2Y_B 1 )
   endif(UNRES_MD_FF STREQUAL "E0LL2Y")
 
 else(NOT UNRES_WITH_MPI)
@@ -489,6 +507,9 @@ else(NOT UNRES_WITH_MPI)
     add_test(NAME UNRES_MIN_INT COMMAND mpiexec ${boot_lam} -np 1 ${CMAKE_CURRENT_BINARY_DIR}/test_prota_E0LL2Y.sh 1l2y_MIN_INT 1 )
     add_test(NAME UNRES_REGULAR COMMAND mpiexec ${boot_lam} -np 1 ${CMAKE_CURRENT_BINARY_DIR}/test_prota_E0LL2Y.sh 1l2y_MIN_REGULAR_INT 1 )
     add_test(NAME UNRES_MD_microcanonical COMMAND mpiexec ${boot_lam} -np 2 ${CMAKE_CURRENT_BINARY_DIR}/test_prota_E0LL2Y.sh 1l2y_micro 2 )
+    add_test(NAME UNRES_Langevin COMMAND mpiexec ${boot_lam} -np 2 ${CMAKE_CURRENT_BINARY_DIR}/test_prota_E0LL2Y.sh 1L2Y_L 2 )
+    add_test(NAME UNRES_NoseHoover COMMAND mpiexec ${boot_lam} -np 2 ${CMAKE_CURRENT_BINARY_DIR}/test_prota_E0LL2Y.sh 1L2Y_NH 2 )
+    add_test(NAME UNRES_Berendsen COMMAND mpiexec ${boot_lam} -np 2 ${CMAKE_CURRENT_BINARY_DIR}/test_prota_E0LL2Y.sh 1L2Y_B 2 )
   endif(UNRES_MD_FF STREQUAL "E0LL2Y")
 
 endif(NOT UNRES_WITH_MPI)