From aa3325db240484cc3d5142ab70985dde83d66ef3 Mon Sep 17 00:00:00 2001 From: Cezary Czaplewski Date: Thu, 10 Mar 2016 17:29:30 +0100 Subject: [PATCH] ctest unres src_MIN and changes for the same energy as src_MD --- source/unres/src_MIN/CMakeLists.txt | 151 +++++++++++++++++++++++++++ source/unres/src_MIN/energy_p_new_barrier.F | 72 +++++++++---- 2 files changed, 205 insertions(+), 18 deletions(-) diff --git a/source/unres/src_MIN/CMakeLists.txt b/source/unres/src_MIN/CMakeLists.txt index 6138673..fadf127 100644 --- a/source/unres/src_MIN/CMakeLists.txt +++ b/source/unres/src_MIN/CMakeLists.txt @@ -222,3 +222,154 @@ set_property(TARGET UNRES_MIN_BIN PROPERTY RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINA #========================================= install(TARGETS UNRES_MIN_BIN DESTINATION ${CMAKE_INSTALL_PREFIX}/unres/MINIM) +#========================================= +# TESTS +#========================================= + +#========================================= +# Generate data test files +#========================================= +# test_single_ala.sh +#========================================= + +FILE(WRITE ${CMAKE_CURRENT_BINARY_DIR}/scripts/test_single_ala.sh +"#!/bin/sh +export POT=GB +export PREFIX=ala10 +#----------------------------------------------------------------------------- +UNRES_BIN=${CMAKE_BINARY_DIR}/bin/${UNRES_BIN} +#----------------------------------------------------------------------------- +DD=${CMAKE_SOURCE_DIR}/PARAM +export BONDPAR=$DD/bond.parm +export THETPAR=$DD/thetaml.5parm +export ROTPAR=$DD/scgauss.parm +export TORPAR=$DD/torsion_631Gdp.parm +export TORDPAR=$DD/torsion_double_631Gdp.parm +export ELEPAR=$DD/electr_631Gdp.parm +export SIDEPAR=$DD/sc_GB_opt.1gab_3S_qclass5no310-shan2-sc-16-10-8k +export FOURIER=$DD/fourier_opt.parm.1igd_hc_iter3_3 +export SCPPAR=$DD/scp.parm +export SCCORPAR=$DD/sccor_pdb_shelly.dat +export PATTERN=$DD/patterns.cart +#----------------------------------------------------------------------------- +$UNRES_BIN +") + +# +# File permissions workaround +# +FILE( COPY ${CMAKE_CURRENT_BINARY_DIR}/scripts/test_single_ala.sh + DESTINATION ${CMAKE_CURRENT_BINARY_DIR} + FILE_PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE +) + + + +#========================================= +# ala10.inp +#========================================= + +file(WRITE ${CMAKE_CURRENT_BINARY_DIR}/ala10.inp +"ala10 unblocked +SEED=-1111333 MD ONE_LETTER rescale_mode=2 +nstep=15000 ntwe=100 ntwx=1000 dt=0.1 lang=0 tbf t_bath=300 damax=1.0 & +reset_moment=1000 reset_vel=1000 +WLONG=1.35279 WSCP=1.59304 WELEC=0.71534 WBOND=1.00000 WANG=1.13873 & +WSCLOC=0.16258 WTOR=1.98599 WTORD=1.57069 WCORRH=0.42887 WCORR5=0.00000 & +WCORR6=0.00000 WEL_LOC=0.16036 WTURN3=1.68722 WTURN4=0.66230 WTURN6=0.00000 & +WVDWPP=0.11371 WHPB=1.00000 & +CUTOFF=7.00000 WCORR4=0.00000 +12 +XAAAAAAAAAAX + 0 + 0 + 90.0000 90.0000 90.0000 90.000 90.000 90.000 90.000 90.000 + 90.0000 90.0000 + 180.0000 180.0000 180.0000 180.000 180.000 180.000 180.000 180.000 + 180.0000 + 110.0000 110.0000 110.0000 100.000 110.000 100.000 110.000 110.000 + 110.0000 110.0000 + -120.0000 -120.0000 -120.000 -120.000 -120.000 -120.000 -120.000 -120.000 + -120.0000 -120.0000 +") + +if(UNRES_MD_FF STREQUAL "E0LL2Y") + +FILE(COPY ${CMAKE_SOURCE_DIR}/ctest/prota.pdb + DESTINATION ${CMAKE_CURRENT_BINARY_DIR} ) + +FILE(COPY ${CMAKE_SOURCE_DIR}/ctest/prota_ENE.inp + DESTINATION ${CMAKE_CURRENT_BINARY_DIR} ) + +FILE(COPY ${CMAKE_SOURCE_DIR}/ctest/prota_MIN_CART.inp + DESTINATION ${CMAKE_CURRENT_BINARY_DIR} ) + +FILE(COPY ${CMAKE_SOURCE_DIR}/ctest/1l2y_MIN_INT.inp + DESTINATION ${CMAKE_CURRENT_BINARY_DIR} ) + +FILE(COPY ${CMAKE_SOURCE_DIR}/ctest/1L2Y.pdb + DESTINATION ${CMAKE_CURRENT_BINARY_DIR} ) + +FILE(COPY ${CMAKE_SOURCE_DIR}/ctest/prota_unres_energy_check.sh + 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 +#========================================= + +FILE(WRITE ${CMAKE_CURRENT_BINARY_DIR}/scripts/test_prota_E0LL2Y.sh +"#!/bin/sh +export POT=GB +export FGPROCS=$2 +export PREFIX=$1 +#----------------------------------------------------------------------------- +UNRES_BIN=${CMAKE_BINARY_DIR}/bin/${UNRES_BIN} +#----------------------------------------------------------------------------- +DD=${CMAKE_SOURCE_DIR}/PARAM +export BONDPAR=$DD/bond_AM1.parm +export THETPAR=$DD/theta_abinitio.parm +export ROTPAR=$DD/rotamers_AM1_aura.10022007.parm +export TORPAR=$DD/torsion_631Gdp.parm +export TORDPAR=$DD/torsion_double_631Gdp.parm +export ELEPAR=$DD/electr_631Gdp.parm +export SIDEPAR=$DD/scinter_$POT.parm +export FOURIER=$DD/fourier_opt.parm.1igd_hc_iter3_3 +export SCPPAR=$DD/scp.parm +export SCCORPAR=$DD/sccor_am1_pawel.dat +export THETPARPDB=$DD/thetaml.5parm +export ROTPARPDB=$DD/scgauss.parm +export PATTERN=$DD/patterns.cart +#----------------------------------------------------------------------------- +echo CTEST_FULL_OUTPUT +$UNRES_BIN +./prota_unres_energy_check.sh $1 ${UNRES_BIN} +") + +# +# File permissions workaround +# +FILE( COPY ${CMAKE_CURRENT_BINARY_DIR}/scripts/test_prota_E0LL2Y.sh + DESTINATION ${CMAKE_CURRENT_BINARY_DIR} + FILE_PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE +) + +endif(UNRES_MD_FF STREQUAL "E0LL2Y") + + + +# Add tests + + if(UNRES_MD_FF STREQUAL "GAB") + add_test(NAME UNRES-MIN_Ala10 COMMAND sh ${CMAKE_CURRENT_BINARY_DIR}/test_single_ala.sh ) + endif(UNRES_MD_FF STREQUAL "GAB") + + if(UNRES_MD_FF STREQUAL "E0LL2Y") + add_test(NAME UNRES-MIN_ENE_prota COMMAND sh ${CMAKE_CURRENT_BINARY_DIR}/test_prota_E0LL2Y.sh prota_ENE 1 ) + add_test(NAME UNRES-MIN_MIN_prota COMMAND sh ${CMAKE_CURRENT_BINARY_DIR}/test_prota_E0LL2Y.sh prota_MIN_CART 1 ) + add_test(NAME UNRES-MIN_MIN_INT COMMAND sh ${CMAKE_CURRENT_BINARY_DIR}/test_prota_E0LL2Y.sh 1l2y_MIN_INT 1 ) + endif(UNRES_MD_FF STREQUAL "E0LL2Y") + + diff --git a/source/unres/src_MIN/energy_p_new_barrier.F b/source/unres/src_MIN/energy_p_new_barrier.F index 3f19648..7095298 100644 --- a/source/unres/src_MIN/energy_p_new_barrier.F +++ b/source/unres/src_MIN/energy_p_new_barrier.F @@ -4314,6 +4314,8 @@ c do i=ibondp_start,ibondp_end diff = vbld(i)-vbldp0 c write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff + if (energy_dec) write (iout,'(a7,i5,4f7.3)') + & "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff estr=estr+diff*diff do j=1,3 gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i) @@ -4332,6 +4334,12 @@ c diff=vbld(i+nres)-vbldsc0(1,iti) c write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff, c & AKSC(1,iti),AKSC(1,iti)*diff*diff + if (energy_dec) then + write (iout,*) + & "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff, + & AKSC(1,iti),AKSC(1,iti)*diff*diff + call flush(iout) + endif estr=estr+0.5d0*AKSC(1,iti)*diff*diff do j=1,3 gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres) @@ -4610,6 +4618,8 @@ C logical lprn /.false./, lprn1 /.false./ etheta=0.0D0 do i=ithet_start,ithet_end + if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or. + &(itype(i).eq.ntyp1)) cycle dethetai=0.0d0 dephii=0.0d0 dephii1=0.0d0 @@ -4619,7 +4629,8 @@ C coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) enddo - if (i.gt.3) then +C if (i.gt.3) then + if (i.gt.3 .and. itype(imax0(i-3,1)).ne.ntyp1) then #ifdef OSF phii=phi(i) if (phii.ne.phii) phii=150.0 @@ -4633,13 +4644,13 @@ C enddo else phii=0.0d0 - ityp1=nthetyp+1 + ityp1=ithetyp(itype(i-2)) do k=1,nsingle cosph1(k)=0.0d0 sinph1(k)=0.0d0 enddo endif - if (i.lt.nres) then + if ((i.lt.nres).and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) if (phii1.ne.phii1) phii1=150.0 @@ -4654,7 +4665,7 @@ C enddo else phii1=0.0d0 - ityp3=nthetyp+1 + ityp3=ithetyp(itype(i)) do k=1,nsingle cosph2(k)=0.0d0 sinph2(k)=0.0d0 @@ -4764,6 +4775,8 @@ C & i,theta(i)*rad2deg,phii*rad2deg, & phii1*rad2deg,ethetai etheta=etheta+ethetai + if (energy_dec) write (iout,'(a6,i5,0pf7.3)') + & 'ebend',i,ethetai if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1 gloc(nphi+i-2,icg)=wang*dethetai @@ -5704,12 +5717,14 @@ C 6/23/01 Compute double torsional energy include 'COMMON.IOUNITS' include 'COMMON.FFIELD' include 'COMMON.TORCNSTR' + include 'COMMON.CONTROL' logical lprn C Set lprn=.true. for debugging lprn=.false. c lprn=.true. etors_d=0.0D0 do i=iphid_start,iphid_end + etors_d_ii=0.0D0 itori=itortyp(itype(i-2)) itori1=itortyp(itype(i-1)) itori2=itortyp(itype(i)) @@ -5729,6 +5744,8 @@ C Regular cosine and sine terms sinphi2=dsin(j*phii1) etors_d=etors_d+v1cij*cosphi1+v1sij*sinphi1+ & v2cij*cosphi2+v2sij*sinphi2 + if (energy_dec) etors_d_ii=etors_d_ii+ + & v1cij*cosphi1+v1sij*sinphi1+v2cij*cosphi2+v2sij*sinphi2 gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1) gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2) enddo @@ -5744,12 +5761,17 @@ C Regular cosine and sine terms sinphi1m2=dsin(l*phii-(k-l)*phii1) etors_d=etors_d+v1cdij*cosphi1p2+v2cdij*cosphi1m2+ & v1sdij*sinphi1p2+v2sdij*sinphi1m2 + if (energy_dec) etors_d_ii=etors_d_ii+ + & v1cdij*cosphi1p2+v2cdij*cosphi1m2+ + & v1sdij*sinphi1p2+v2sdij*sinphi1m2 gloci1=gloci1+l*(v1sdij*cosphi1p2+v2sdij*cosphi1m2 & -v1cdij*sinphi1p2-v2cdij*sinphi1m2) gloci2=gloci2+(k-l)*(v1sdij*cosphi1p2-v2sdij*cosphi1m2 & -v1cdij*sinphi1p2+v2cdij*sinphi1m2) enddo enddo + if (energy_dec) write (iout,'(a6,i5,0pf7.3)') + & 'etor_d',i,etors_d_ii gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*gloci1 gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*gloci2 enddo @@ -5787,11 +5809,21 @@ c write (iout,*) "EBACK_SC_COR",itau_start,itau_end esccor=0.0D0 do i=itau_start,itau_end esccor_ii=0.0D0 + if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle isccori=isccortyp(itype(i-2)) isccori1=isccortyp(itype(i-1)) c write (iout,*) "EBACK_SC_COR",i,nterm_sccor(isccori,isccori1) phii=phi(i) - do intertyp=1,3 !intertyp + +cccc Added 9 May 2012 +cc Tauangle is torsional engle depending on the value of first digit +c(see comment below) +cc Omicron is flat angle depending on the value of first digit +c(see comment below) +C print *,i,tauangle(1,i) + +c do intertyp=1,3 !intertyp + do intertyp=2,2 !intertyp cc Added 09 May 2012 (Adasko) cc Intertyp means interaction type of backbone mainchain correlation: c 1 = SC...Ca...Ca...Ca @@ -5799,18 +5831,16 @@ c 2 = Ca...Ca...Ca...SC c 3 = SC...Ca...Ca...SCi gloci=0.0D0 if (((intertyp.eq.3).and.((itype(i-2).eq.10).or. - & (itype(i-1).eq.10).or.(itype(i-2).eq.ntyp1).or. - & (itype(i-1).eq.ntyp1))) + & (itype(i-1).eq.10).or.(itype(i-2).eq.21).or. + & (itype(i-1).eq.21))) & .or. ((intertyp.eq.1).and.((itype(i-2).eq.10) - & .or.(itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1) - & .or.(itype(i).eq.ntyp1))) + & .or.(itype(i-2).eq.21))) & .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or. - & (itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or. - & (itype(i-3).eq.ntyp1)))) cycle - if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.ntyp1)) cycle - if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.ntyp1)) + & (itype(i-1).eq.21)))) cycle + if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.21)) cycle + if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.21)) & cycle - do j=1,nterm_sccor(isccori,isccori1) + do j=1,nterm_sccor(isccori,isccori1) v1ij=v1sccor(j,intertyp,isccori,isccori1) v2ij=v2sccor(j,intertyp,isccori,isccori1) cosphi=dcos(j*tauangle(intertyp,i)) @@ -5820,14 +5850,20 @@ c 3 = SC...Ca...Ca...SCi enddo c write (iout,*) "EBACK_SC_COR",i,esccor,intertyp gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci +c write (iout,*) "WTF",intertyp,i,itype(i),v1ij*cosphi+v2ij*sinphi +c &gloc_sc(intertyp,i-3,icg) if (lprn) & write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') - & restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,isccori,isccori1, - & (v1sccor(j,intertyp,isccori,isccori1),j=1,6) - & ,(v2sccor(j,intertyp,isccori,isccori1),j=1,6) -C gsccor_loc(i-3)=gsccor_loc(i-3)+gloci + & restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1, + & (v1sccor(j,intertyp,itori,itori1),j=1,6) + & ,(v2sccor(j,intertyp,itori,itori1),j=1,6) + gsccor_loc(i-3)=gsccor_loc(i-3)+gloci enddo !intertyp enddo +c do i=1,nres +c write (iout,*) "W@T@F", gloc_sc(1,i,icg),gloc_sc(2,i,icg), +c & gloc_sc(3,i,icg) +c enddo return end c---------------------------------------------------------------------------- -- 1.7.9.5