From: Cezary Czaplewski Date: Tue, 8 Mar 2016 01:28:55 +0000 (+0100) Subject: ctest for thermostats: Langevin, Nose-Hoover and Berendsen X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=40bf66aba1eb110e45e3c01a2001428f6864f300 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 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 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 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 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 --- 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 "`pwd`/${file_stat}.png" + + 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 "`pwd`/${file_stat}.png" + + 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 --- 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)