From b9025a366733ac0290debbd52ede3527d29c9ed5 Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Mon, 13 Apr 2015 17:07:52 +0200 Subject: [PATCH] Introduction of SS to newcorr and SSS to src_MD-M --- source/unres/src_MD-M-newcorr/CMakeLists.txt | 56 +- source/unres/src_MD-M-newcorr/COMMON.DERIV | 2 +- source/unres/src_MD-M-newcorr/COMMON.SBRIDGE | 17 +- source/unres/src_MD-M-newcorr/MREMD.F | 20 +- .../unres/src_MD-M-newcorr/energy_p_new_barrier.F | 29 +- source/unres/src_MD-M-newcorr/geomout.F | 30 +- source/unres/src_MD-M-newcorr/initialize_p.F | 1 + source/unres/src_MD-M-newcorr/parmread.F | 4 +- source/unres/src_MD-M-newcorr/readrtns_CSA.F | 56 +- source/unres/src_MD-M-newcorr/ssMD.F | 1951 ++++++++++++++++++++ source/unres/src_MD-M/energy_p_new_barrier.F | 14 + source/unres/src_MD-M/parmread.F | 2 +- source/unres/src_MD-M/readrtns_CSA.F | 4 + source/unres/src_MD-M/ssMD.F | 148 ++ 14 files changed, 2279 insertions(+), 55 deletions(-) create mode 100644 source/unres/src_MD-M-newcorr/ssMD.F diff --git a/source/unres/src_MD-M-newcorr/CMakeLists.txt b/source/unres/src_MD-M-newcorr/CMakeLists.txt index b65d077..29291da 100644 --- a/source/unres/src_MD-M-newcorr/CMakeLists.txt +++ b/source/unres/src_MD-M-newcorr/CMakeLists.txt @@ -5,11 +5,6 @@ enable_language (Fortran) #================================ -# build the xdrf library -#================================ -#add_subdirectory(xdrf) - -#================================ # Set source file lists #================================ set(UNRES_MDM_SRC0 @@ -35,6 +30,7 @@ set(UNRES_MDM_SRC0 distfit.f djacob.f econstr_local.F + eigen.f elecont.f energy_split-sep.F entmcm.F @@ -70,6 +66,7 @@ set(UNRES_MDM_SRC0 permut.F pinorm.f printmat.f + prng_32.F q_measure.F ran.f randgens.f @@ -92,17 +89,9 @@ set(UNRES_MDM_SRC0 timing.F together.F unres.F + ssMD.F ) -if (Fortran_COMPILER_NAME STREQUAL "ifort") - set(UNRES_MDM_SRC0 ${UNRES_MDM_SRC0} prng.f ) -elseif(Fortran_COMPILER_NAME STREQUAL "mpif90") - set(UNRES_MDM_SRC0 ${UNRES_MDM_SRC0} prng.f ) -else() - set(UNRES_MDM_SRC0 ${UNRES_MDM_SRC0} prng_32.F ) -endif (Fortran_COMPILER_NAME STREQUAL "ifort") - -set(UNRES_MDM_SRC2 eigen.f) set(UNRES_MDM_SRC3 energy_p_new_barrier.F energy_p_new-sep_barrier.F gradient_p.F ) set(UNRES_MDM_PP_SRC @@ -140,6 +129,7 @@ set(UNRES_MDM_PP_SRC newconf.f parmread.F permut.F + prng_32.F q_measure1.F q_measure3.F q_measure.F @@ -172,10 +162,10 @@ endif(NOT Fortran_COMPILER_NAME STREQUAL "ifort") # Set comipiler flags for different sourcefiles #================================================ if (Fortran_COMPILER_NAME STREQUAL "ifort") - set(FFLAGS0 "-CB -g -ip -w" ) - set(FFLAGS1 "-w -g " ) + set(FFLAGS0 "-ip -w" ) + set(FFLAGS1 "-w -g -d2 -CA -CB" ) set(FFLAGS2 "-w -g -00 ") - set(FFLAGS3 "-CB -g -w -ipo " ) + set(FFLAGS3 "-w -ipo " ) elseif (Fortran_COMPILER_NAME STREQUAL "gfortran") set(FFLAGS0 "-std=legacy -I. " ) set(FFLAGS1 "-std=legacy -g -I. " ) @@ -186,15 +176,15 @@ endif (Fortran_COMPILER_NAME STREQUAL "ifort") # Add MPI compiler flags if(UNRES_WITH_MPI) - set(FFLAGS0 "${FFLAGS0} -I${MPIF_INCLUDE_DIRECTORIES}") - set(FFLAGS1 "${FFLAGS1} -I${MPIF_INCLUDE_DIRECTORIES}") - set(FFLAGS2 "${FFLAGS2} -I${MPIF_INCLUDE_DIRECTORIES}") - set(FFLAGS3 "${FFLAGS3} -I${MPIF_INCLUDE_DIRECTORIES}") + set(FFLAGS0 "${FFLAGS0} -I${MPI_Fortran_INCLUDE_PATH}") + set(FFLAGS1 "${FFLAGS1} -I${MPI_Fortran_INCLUDE_PATH}") + set(FFLAGS2 "${FFLAGS2} -I${MPI_Fortran_INCLUDE_PATH}") + set(FFLAGS3 "${FFLAGS3} -I${MPI_Fortran_INCLUDE_PATH}") endif(UNRES_WITH_MPI) set_property(SOURCE ${UNRES_MDM_SRC0} APPEND PROPERTY COMPILE_FLAGS ${FFLAGS0} ) -set_property(SOURCE ${UNRES_MD_SRC1} PROPERTY COMPILE_FLAGS ${FFLAGS1} ) -set_property(SOURCE ${UNRES_MD_SRC2} PROPERTY COMPILE_FLAGS ${FFLAGS2} ) +#set_property(SOURCE ${UNRES_MD_SRC1} PROPERTY COMPILE_FLAGS ${FFLAGS1} ) +#set_property(SOURCE ${UNRES_MD_SRC2} PROPERTY COMPILE_FLAGS ${FFLAGS2} ) set_property(SOURCE ${UNRES_MDM_SRC3} PROPERTY COMPILE_FLAGS ${FFLAGS3} ) #========================================= @@ -202,7 +192,7 @@ set_property(SOURCE ${UNRES_MDM_SRC3} PROPERTY COMPILE_FLAGS ${FFLAGS3} ) #========================================= if(UNRES_MD_FF STREQUAL "GAB" ) # set preprocesor flags - set(CPPFLAGS "PROCOR -DUNRES -DISNAN -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -DSCCORPDB -DNEWCORR" ) + set(CPPFLAGS "PROCOR -DUNRES -DISNAN -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -DSCCORPDB -DTIMING -DTIMING_ENE" ) #========================================= # Settings for E0LL2Y force field @@ -248,7 +238,9 @@ if (UNRES_WITH_MPI) endif(UNRES_WITH_MPI) +#========================================= # add 64-bit specific preprocessor flags +#========================================= if (architektura STREQUAL "64") set(CPPFLAGS "${CPPFLAGS} -DAMD64") endif (architektura STREQUAL "64") @@ -267,10 +259,10 @@ set_property(SOURCE proc_proc.c PROPERTY COMPILE_DEFINITIONS "SGI" ) #======================================== if(UNRES_WITH_MPI) # binary with mpi - set(UNRES_BIN "unres_${Fortran_COMPILER_NAME}_MPICH_${UNRES_FF}.exe") + set(UNRES_BIN "unresMD-M_${Fortran_COMPILER_NAME}_MPICH_${UNRES_MD_FF}.exe") else(UNRES_WITH_MPI) # binary without mpi - set(UNRES_BIN "unres_${Fortran_COMPILER_NAME}_single_${UNRES_FF}.exe") + set(UNRES_BIN "unresMD-M_${Fortran_COMPILER_NAME}_single_${UNRES_MD_FF}.exe") endif(UNRES_WITH_MPI) #========================================= @@ -308,15 +300,14 @@ set_property(SOURCE ${CMAKE_CURRENT_BINARY_DIR}/cinfo.f PROPERTY COMPILE_FLAGS " #========================================= # Set full unres MD-M sources #========================================= -set(UNRES_MDM_SRCS ${UNRES_MDM_SRC0} ${UNRES_MDM_SRC2} ${UNRES_MDM_SRC3} ${CMAKE_CURRENT_BINARY_DIR}/cinfo.f proc_proc.c ) +set(UNRES_MDM_SRCS ${UNRES_MDM_SRC0} ${UNRES_MDM_SRC3} ${CMAKE_CURRENT_BINARY_DIR}/cinfo.f proc_proc.c ) #========================================= # Build the binary #========================================= add_executable(UNRES_BIN-MD-M ${UNRES_MDM_SRCS} ) set_target_properties(UNRES_BIN-MD-M PROPERTIES OUTPUT_NAME ${UNRES_BIN}) - -#set_property(TARGET ${UNRES_BIN} PROPERTY RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin/unres/MD ) +set_property(TARGET UNRES_BIN-MD-M PROPERTY RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin ) #add_dependencies (${UNRES_BIN} ${UNRES_XDRFLIB}) #========================================= @@ -324,13 +315,18 @@ set_target_properties(UNRES_BIN-MD-M PROPERTIES OUTPUT_NAME ${UNRES_BIN}) #========================================= # link MPI library (libmpich.a) if(UNRES_WITH_MPI) - target_link_libraries( UNRES_BIN-MD-M ${MPIF_LIBRARIES} ) + target_link_libraries( UNRES_BIN-MD-M ${MPI_Fortran_LIBRARIES} ) endif(UNRES_WITH_MPI) # link libxdrf.a #message("UNRES_XDRFLIB=${UNRES_XDRFLIB}") target_link_libraries( UNRES_BIN-MD-M xdrf ) #========================================= +# Install Path +#========================================= +install(TARGETS UNRES_BIN-MD-M DESTINATION ${CMAKE_INSTALL_PREFIX}) + +#========================================= # TESTS #========================================= diff --git a/source/unres/src_MD-M-newcorr/COMMON.DERIV b/source/unres/src_MD-M-newcorr/COMMON.DERIV index 524d72a..d7c98bd 100644 --- a/source/unres/src_MD-M-newcorr/COMMON.DERIV +++ b/source/unres/src_MD-M-newcorr/COMMON.DERIV @@ -2,7 +2,7 @@ & gvdwpp,gel_loc,gel_loc_long,gvdwc_scpp, & gradx_scp,gvdwc_scp,ghpbx,ghpbc,gloc,gloc_x,dtheta,dphi,dalpha, & domega,gscloc,gsclocx,gradcorr,gradcorr_long,gradcorr5_long, - & gradcorr6_long,gcorr6_turn_long + & gradcorr6_long,gcorr6_turn_long,gvdwx integer nfl,icg common /derivat/ dcdv(6,maxdim),dxdv(6,maxdim),dxds(6,maxres), & gradx(3,maxres,2),gradc(3,maxres,2),gvdwx(3,maxres), diff --git a/source/unres/src_MD-M-newcorr/COMMON.SBRIDGE b/source/unres/src_MD-M-newcorr/COMMON.SBRIDGE index 4cc80c8..91dd2cd 100644 --- a/source/unres/src_MD-M-newcorr/COMMON.SBRIDGE +++ b/source/unres/src_MD-M-newcorr/COMMON.SBRIDGE @@ -1,12 +1,17 @@ - double precision ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss + double precision ss_depth,ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss integer ns,nss,nfree,iss - common /sbridge/ ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss, + common /sbridge/ ss_depth,ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss, & ns,nss,nfree,iss(maxss) - double precision dhpb,forcon - integer ihpb,jhpb,nhpb - common /links/ dhpb(maxdim),forcon(maxdim),ihpb(maxdim), - & jhpb(maxdim),nhpb + double precision dhpb,dhpb1,forcon + integer ihpb,jhpb,nhpb,idssb,jdssb + common /links/ dhpb(maxdim),dhpb1(maxdim),forcon(maxdim), + & ihpb(maxdim),jhpb(maxdim),ibecarb(maxdim),nhpb double precision weidis common /restraints/ weidis integer link_start,link_end common /links_split/ link_start,link_end + double precision Ht,dyn_ssbond_ij + logical dyn_ss,dyn_ss_mask + common /dyn_ssbond/ dyn_ssbond_ij(maxres,maxres), + & idssb(maxdim),jdssb(maxdim), + & Ht,dyn_ss,dyn_ss_mask(maxres) diff --git a/source/unres/src_MD-M-newcorr/MREMD.F b/source/unres/src_MD-M-newcorr/MREMD.F index d55a95b..cac85aa 100644 --- a/source/unres/src_MD-M-newcorr/MREMD.F +++ b/source/unres/src_MD-M-newcorr/MREMD.F @@ -1500,8 +1500,15 @@ c end debugging call xdrffloat_(ixdrf, real(t_restart1(4,il)), iret) call xdrfint_(ixdrf, nss, iret) do j=1,nss - call xdrfint_(ixdrf, ihpb(j), iret) - call xdrfint_(ixdrf, jhpb(j), iret) +C do j=1,nss + if (dyn_ss) then + call xdrfint(ixdrf, idssb(j)+nres, iret) + call xdrfint(ixdrf, jdssb(j)+nres, iret) + else + call xdrfint_(ixdrf, ihpb(j), iret) + call xdrfint_(ixdrf, jhpb(j), iret) + endif + enddo enddo call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret) call xdrfint_(ixdrf, iset_restart1(il), iret) @@ -1538,8 +1545,13 @@ c end debugging call xdrffloat(ixdrf, real(t_restart1(4,il)), iret) call xdrfint(ixdrf, nss, iret) do j=1,nss - call xdrfint(ixdrf, ihpb(j), iret) - call xdrfint(ixdrf, jhpb(j), iret) + if (dyn_ss) then + call xdrfint(ixdrf, idssb(j)+nres, iret) + call xdrfint(ixdrf, jdssb(j)+nres, iret) + else + call xdrfint(ixdrf, ihpb(j), iret) + call xdrfint(ixdrf, jhpb(j), iret) + endif enddo call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret) call xdrfint(ixdrf, iset_restart1(il), iret) diff --git a/source/unres/src_MD-M-newcorr/energy_p_new_barrier.F b/source/unres/src_MD-M-newcorr/energy_p_new_barrier.F index 42005e1..5e90c17 100644 --- a/source/unres/src_MD-M-newcorr/energy_p_new_barrier.F +++ b/source/unres/src_MD-M-newcorr/energy_p_new_barrier.F @@ -121,6 +121,10 @@ C C Calculate electrostatic (H-bonding) energy of the main chain. C 107 continue +cmc +cmc Sep-06: egb takes care of dynamic ss bonds too +cmc +C if (dyn_ss) call dyn_set_nss c print *,"Processor",myrank," computed USCSC" #ifdef TIMING time01=MPI_Wtime() @@ -300,6 +304,7 @@ c Here are the energies showed per procesor if the are more processors c per molecule then we sum it up in sum_energy subroutine c print *," Processor",myrank," calls SUM_ENERGY" call sum_energy(energia,.true.) + if (dyn_ss) call dyn_set_nss c print *," Processor",myrank," left SUM_ENERGY" #ifdef TIMING time_sumene=time_sumene+MPI_Wtime()-time00 @@ -1412,6 +1417,7 @@ C include 'COMMON.IOUNITS' include 'COMMON.CALC' include 'COMMON.CONTROL' + include 'COMMON.SBRIDGE' logical lprn evdw=0.0D0 ccccc energy_dec=.false. @@ -1439,6 +1445,12 @@ C Calculate SC interaction energy. C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) + IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN + call dyn_ssbond_ene(i,j,evdwij) + evdw=evdw+evdwij + if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') + & 'evdw',i,j,evdwij,' ss' + ELSE ind=ind+1 itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle @@ -1533,6 +1545,7 @@ C Calculate the radial part of the gradient gg(3)=zj*fac C Calculate angular part of the gradient. call sc_grad + endif ! dyn_ss enddo ! j enddo ! iint enddo ! i @@ -4280,10 +4293,20 @@ C iii and jjj point to the residues for which the distance is assigned. cd write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj C 24/11/03 AL: SS bridges handled separately because of introducing a specific C distance and angle dependent SS bond potential. - if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. - & iabs(itype(jjj)).eq.1) then - call ssbond_ene(iii,jjj,eij) +cmc if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then +C 18/07/06 MC: Use the convention that the first nss pairs are SS bonds +c if (.not.dyn_ss .and. i.le.nss) then +C 15/02/13 CC dynamic SSbond +C if (.not.dyn_ss.and. +C & ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then + + if (.not.dyn_ss .and. i.le.nss) then +C 15/02/13 CC dynamic SSbond - additional check + if (ii.gt.nres + & .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then + call ssbond_ene(iii,jjj,eij) ehpb=ehpb+2*eij + endif cd write (iout,*) "eij",eij else C Calculate the distance between the two points and its difference from the diff --git a/source/unres/src_MD-M-newcorr/geomout.F b/source/unres/src_MD-M-newcorr/geomout.F index c6a3944..e869b4a 100644 --- a/source/unres/src_MD-M-newcorr/geomout.F +++ b/source/unres/src_MD-M-newcorr/geomout.F @@ -80,9 +80,15 @@ cmodel write (iunit,'(a5,i6)') 'MODEL',1 if (nss.gt.0) then do i=1,nss + if (dyn_ss) then + write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') + & 'SSBOND',i,'CYS',idssb(i)-nnt+1, + & 'CYS',jdssb(i)-nnt+1 + else write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') - & 'SSBOND',i,'CYS',ihpb(i)-1-nres, - & 'CYS',jhpb(i)-1-nres + & 'SSBOND',i,'CYS',ihpb(i)-nnt+1-nres, + & 'CYS',jhpb(i)-nnt+1-nres + endif enddo endif @@ -124,7 +130,14 @@ cmodel write (iunit,'(a5,i6)') 'MODEL',1 write (iunit,30) ica(nct),ica(nct)+1 endif do i=1,nss +C write (iunit,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1 +c trzeba uporzdkowac +C write (iunit,30) ica(ihpb(i))+1,ica(jhpb(i))+1 + if (dyn_ss) then + write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1 + else write (iunit,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1 + endif enddo write (iunit,'(a6)') 'ENDMDL' 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3) @@ -279,8 +292,14 @@ c---------------------------------------------------------------- open(icart,file=cartname,access="append") #endif write (icart,'(e15.8,2e15.5,f12.5,$)') time,potE,uconst,t_bath - write (icart,'(i4,$)') +C write (icart,'(i4,$)') + if (dyn_ss) then + write (icart,'(i4,$)') + & nss,(idssb(j)+nres,jdssb(j)+nres,j=1,nss) + else + write (icart,'(i4,$)') & nss,(ihpb(j),jhpb(j),j=1,nss) + endif write (icart,'(i4,20f7.4)') nfrag+npair+3*nfrag_back, & (qfrag(i),i=1,nfrag),(qpair(i),i=1,npair), & (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back) @@ -346,8 +365,13 @@ c----------------------------------------------------------------- call xdrffloat(ixdrf, real(t_bath), iret) call xdrfint(ixdrf, nss, iret) do j=1,nss + if (dyn_ss) then + call xdrfint(ixdrf, idssb(j)+nres, iret) + call xdrfint(ixdrf, jdssb(j)+nres, iret) + else call xdrfint(ixdrf, ihpb(j), iret) call xdrfint(ixdrf, jhpb(j), iret) + endif enddo call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret) do i=1,nfrag diff --git a/source/unres/src_MD-M-newcorr/initialize_p.F b/source/unres/src_MD-M-newcorr/initialize_p.F index a650241..d2474fc 100644 --- a/source/unres/src_MD-M-newcorr/initialize_p.F +++ b/source/unres/src_MD-M-newcorr/initialize_p.F @@ -377,6 +377,7 @@ cd write (iout,*) 'ns=',ns,' nss=',nss,' ihpb,jhpb', cd & (ihpb(i),jhpb(i),i=1,nss) do i=nnt,nct-1 scheck=.false. + if (dyn_ss) goto 10 do ii=1,nss if (ihpb(ii).eq.i+nres) then scheck=.true. diff --git a/source/unres/src_MD-M-newcorr/parmread.F b/source/unres/src_MD-M-newcorr/parmread.F index 03efdb2..e48d010 100644 --- a/source/unres/src_MD-M-newcorr/parmread.F +++ b/source/unres/src_MD-M-newcorr/parmread.F @@ -912,11 +912,13 @@ C read (ifourier,*,end=115,err=115) read (ifourier,*,end=115,err=115) (b(ii),ii=1,13) #ifdef NEWCORR + write (iout,*) "TUTUTU" read (ifourier,*,end=115,err=115) (bnew1(ii,1,i),ii=1,3) read (ifourier,*,end=115,err=115) (bnew2(ii,1,i),ii=1,3) read (ifourier,*,end=115,err=115) (bnew1(ii,2,i),ii=1,1) read (ifourier,*,end=115,err=115) (bnew2(ii,2,i),ii=1,1) read (ifourier,*,end=115,err=115) (eenew(ii,i),ii=1,1) + #endif if (lprint) then write (iout,*) 'Type',i @@ -1266,7 +1268,7 @@ c lprint=.false. C C Define the constants of the disulfide bridge C - ebr=-5.50D0 + ebr=-12.00D0 c c Old arbitrary potential - commented out. c diff --git a/source/unres/src_MD-M-newcorr/readrtns_CSA.F b/source/unres/src_MD-M-newcorr/readrtns_CSA.F index 265b705..078ff5e 100644 --- a/source/unres/src_MD-M-newcorr/readrtns_CSA.F +++ b/source/unres/src_MD-M-newcorr/readrtns_CSA.F @@ -679,12 +679,36 @@ C 12/1/95 Added weight for the multi-body term WCORR call reada(weightcard,"V2SS",v2ss,7.61d0) call reada(weightcard,"V3SS",v3ss,13.7d0) call reada(weightcard,"EBR",ebr,-5.50D0) + dyn_ss=(index(weightcard,'DYN_SS').gt.0) + do i=1,maxres + dyn_ss_mask(i)=.false. + enddo + do i=1,maxres-1 + do j=i+1,maxres + dyn_ssbond_ij(i,j)=1.0d300 + enddo + enddo + call reada(weightcard,"HT",Ht,0.0D0) + if (dyn_ss) then + ss_depth=ebr/wsc-0.25*eps(1,1) + Ht=Ht/wsc-0.25*eps(1,1) + akcm=akcm*wstrain/wsc + akth=akth*wstrain/wsc + akct=akct*wstrain/wsc + v1ss=v1ss*wstrain/wsc + v2ss=v2ss*wstrain/wsc + v3ss=v3ss*wstrain/wsc + else + ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain + endif + if(me.eq.king.or..not.out1file) then write (iout,*) "Parameters of the SS-bond potential:" write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth, & " AKCT",akct write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss - write (iout,*) "EBR",ebr + write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth + write (iout,*)" HT",Ht c print *,'indpdb=',indpdb,' pdbref=',pdbref endif if (indpdb.gt.0 .or. pdbref) then @@ -1055,19 +1079,36 @@ C Generate distance constraints, if the PDB structure is to be regularized. write (iout,'(/a,i3,a)') & 'The chain contains',ns,' disulfide-bridging cysteines.' write (iout,'(20i4)') (iss(i),i=1,ns) + if (dyn_ss) then + write(iout,*)"Running with dynamic disulfide-bond formation" + else write (iout,'(/a/)') 'Pre-formed links are:' do i=1,nss i1=ihpb(i)-nres i2=jhpb(i)-nres it1=itype(i1) it2=itype(i2) - if (me.eq.king.or..not.out1file) - & write (iout,'(2a,i3,3a,i3,a,3f10.3)') + write (iout,'(2a,i3,3a,i3,a,3f10.3)') & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i), & ebr,forcon(i) enddo write (iout,'(a)') endif + endif + if (ns.gt.0.and.dyn_ss) then + do i=nss+1,nhpb + ihpb(i-nss)=ihpb(i) + jhpb(i-nss)=jhpb(i) + forcon(i-nss)=forcon(i) + dhpb(i-nss)=dhpb(i) + enddo + nhpb=nhpb-nss + nss=0 + call hpb_partition + do i=1,ns + dyn_ss_mask(iss(i))=.true. + enddo + endif if (i2ndstr.gt.0) call secstrp2dihc c call geom_to_var(nvar,x) c call etotal(energia(0)) @@ -1130,10 +1171,12 @@ C Check whether the specified bridging residues are cystines. do i=1,ns if (itype(iss(i)).ne.1) then if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') - & 'Do you REALLY think that the residue ',restyp(iss(i)),i, + & 'Do you REALLY think that the residue ', + & restyp(itype(iss(i))),i, & ' can form a disulfide bridge?!!!' write (*,'(2a,i3,a)') - & 'Do you REALLY think that the residue ',restyp(iss(i)),i, + & 'Do you REALLY think that the residue ', + & restyp(itype(iss(i))),i, & ' can form a disulfide bridge?!!!' #ifdef MPI call MPI_Finalize(MPI_COMM_WORLD,ierror) @@ -1144,7 +1187,8 @@ C Check whether the specified bridging residues are cystines. C Read preformed bridges. if (ns.gt.0) then read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss) - write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss) + if(fg_rank.eq.0) + &write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss) if (nss.gt.0) then nhpb=nss C Check if the residues involved in bridges are in the specified list of diff --git a/source/unres/src_MD-M-newcorr/ssMD.F b/source/unres/src_MD-M-newcorr/ssMD.F new file mode 100644 index 0000000..15800ae --- /dev/null +++ b/source/unres/src_MD-M-newcorr/ssMD.F @@ -0,0 +1,1951 @@ +c---------------------------------------------------------------------------- + subroutine check_energies +c implicit none + +c Includes + include 'DIMENSIONS' + include 'COMMON.CHAIN' + include 'COMMON.VAR' + include 'COMMON.IOUNITS' + include 'COMMON.SBRIDGE' + include 'COMMON.LOCAL' + include 'COMMON.GEO' + +c External functions + double precision ran_number + external ran_number + +c Local variables + integer i,j,k,l,lmax,p,pmax + double precision rmin,rmax + double precision eij + + double precision d + double precision wi,rij,tj,pj + + +c return + + i=5 + j=14 + + d=dsc(1) + rmin=2.0D0 + rmax=12.0D0 + + lmax=10000 + pmax=1 + + do k=1,3 + c(k,i)=0.0D0 + c(k,j)=0.0D0 + c(k,nres+i)=0.0D0 + c(k,nres+j)=0.0D0 + enddo + + do l=1,lmax + +ct wi=ran_number(0.0D0,pi) +c wi=ran_number(0.0D0,pi/6.0D0) +c wi=0.0D0 +ct tj=ran_number(0.0D0,pi) +ct pj=ran_number(0.0D0,pi) +c pj=ran_number(0.0D0,pi/6.0D0) +c pj=0.0D0 + + do p=1,pmax +ct rij=ran_number(rmin,rmax) + + c(1,j)=d*sin(pj)*cos(tj) + c(2,j)=d*sin(pj)*sin(tj) + c(3,j)=d*cos(pj) + + c(3,nres+i)=-rij + + c(1,i)=d*sin(wi) + c(3,i)=-rij-d*cos(wi) + + do k=1,3 + dc(k,nres+i)=c(k,nres+i)-c(k,i) + dc_norm(k,nres+i)=dc(k,nres+i)/d + dc(k,nres+j)=c(k,nres+j)-c(k,j) + dc_norm(k,nres+j)=dc(k,nres+j)/d + enddo + + call dyn_ssbond_ene(i,j,eij) + enddo + enddo + + call exit(1) + + return + end + +C----------------------------------------------------------------------------- + + subroutine dyn_ssbond_ene(resi,resj,eij) +c implicit none + +c Includes + include 'DIMENSIONS' + include 'COMMON.SBRIDGE' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.VAR' + include 'COMMON.IOUNITS' + include 'COMMON.CALC' +#ifndef CLUST +#ifndef WHAM + include 'COMMON.MD' +#endif +#endif + +c External functions + double precision h_base + external h_base + +c Input arguments + integer resi,resj + +c Output arguments + double precision eij + +c Local variables + logical havebond +c integer itypi,itypj,k,l + double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi + double precision sig0ij,ljd,sig,fac,e1,e2 + double precision dcosom1(3),dcosom2(3),ed + double precision pom1,pom2 + double precision ljA,ljB,ljXs + double precision d_ljB(1:3) + double precision ssA,ssB,ssC,ssXs + double precision ssxm,ljxm,ssm,ljm + double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3) + double precision f1,f2,h1,h2,hd1,hd2 + double precision omega,delta_inv,deltasq_inv,fac1,fac2 +c-------FIRST METHOD + double precision xm,d_xm(1:3) +c-------END FIRST METHOD +c-------SECOND METHOD +c$$$ double precision ss,d_ss(0:3),ljf,d_ljf(0:3) +c-------END SECOND METHOD + +c-------TESTING CODE + logical checkstop,transgrad + common /sschecks/ checkstop,transgrad + + integer icheck,nicheck,jcheck,njcheck + double precision echeck(-1:1),deps,ssx0,ljx0 +c-------END TESTING CODE + + + i=resi + j=resj + + itypi=itype(i) + dxi=dc_norm(1,nres+i) + dyi=dc_norm(2,nres+i) + dzi=dc_norm(3,nres+i) + dsci_inv=vbld_inv(i+nres) + + itypj=itype(j) + xj=c(1,nres+j)-c(1,nres+i) + yj=c(2,nres+j)-c(2,nres+i) + zj=c(3,nres+j)-c(3,nres+i) + dxj=dc_norm(1,nres+j) + dyj=dc_norm(2,nres+j) + dzj=dc_norm(3,nres+j) + dscj_inv=vbld_inv(j+nres) + + chi1=chi(itypi,itypj) + chi2=chi(itypj,itypi) + chi12=chi1*chi2 + chip1=chip(itypi) + chip2=chip(itypj) + chip12=chip1*chip2 + alf1=alp(itypi) + alf2=alp(itypj) + alf12=0.5D0*(alf1+alf2) + + rrij=1.0D0/(xj*xj+yj*yj+zj*zj) + rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse +c The following are set in sc_angular +c erij(1)=xj*rij +c erij(2)=yj*rij +c erij(3)=zj*rij +c om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3) +c om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3) +c om12=dxi*dxj+dyi*dyj+dzi*dzj + call sc_angular + rij=1.0D0/rij ! Reset this so it makes sense + + sig0ij=sigma(itypi,itypj) + sig=sig0ij*dsqrt(1.0D0/sigsq) + + ljXs=sig-sig0ij + ljA=eps1*eps2rt**2*eps3rt**2 + ljB=ljA*bb(itypi,itypj) + ljA=ljA*aa(itypi,itypj) + ljxm=ljXs+(-2.0D0*aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) + + ssXs=d0cm + deltat1=1.0d0-om1 + deltat2=1.0d0+om2 + deltat12=om2-om1+2.0d0 + cosphi=om12-om1*om2 + ssA=akcm + ssB=akct*deltat12 + ssC=ss_depth + & +akth*(deltat1*deltat1+deltat2*deltat2) + & +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi + ssxm=ssXs-0.5D0*ssB/ssA + +c-------TESTING CODE +c$$$c Some extra output +c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA +c$$$ ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj) +c$$$ ssx0=ssB*ssB-4.0d0*ssA*ssC +c$$$ if (ssx0.gt.0.0d0) then +c$$$ ssx0=ssXs+0.5d0*(-ssB+sqrt(ssx0))/ssA +c$$$ else +c$$$ ssx0=ssxm +c$$$ endif +c$$$ ljx0=ljXs+(-aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) +c$$$ write(iout,'(a,4f8.2,2f15.2,3f6.2)')"SSENERGIES ", +c$$$ & ssxm,ljxm,ssx0,ljx0,ssm,ljm,om1,om2,om12 +c$$$ return +c-------END TESTING CODE + +c-------TESTING CODE +c Stop and plot energy and derivative as a function of distance + if (checkstop) then + ssm=ssC-0.25D0*ssB*ssB/ssA + ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj) + if (ssm.lt.ljm .and. + & dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then + nicheck=1000 + njcheck=1 + deps=0.5d-7 + else + checkstop=.false. + endif + endif + if (.not.checkstop) then + nicheck=0 + njcheck=-1 + endif + + do icheck=0,nicheck + do jcheck=-1,njcheck + if (checkstop) rij=(ssxm-1.0d0)+ + & ((ljxm-ssxm+2.0d0)*icheck)/nicheck+jcheck*deps +c-------END TESTING CODE + + if (rij.gt.ljxm) then + havebond=.false. + ljd=rij-ljXs + fac=(1.0D0/ljd)**expon + e1=fac*fac*aa(itypi,itypj) + e2=fac*bb(itypi,itypj) + eij=eps1*eps2rt*eps3rt*(e1+e2) + eps2der=eij*eps3rt + eps3der=eij*eps2rt + eij=eij*eps2rt*eps3rt + + sigder=-sig/sigsq + e1=e1*eps1*eps2rt**2*eps3rt**2 + ed=-expon*(e1+eij)/ljd + sigder=ed*sigder + eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1 + eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2 + eom12=eij*eps1_om12+eps2der*eps2rt_om12 + & -2.0D0*alf12*eps3der+sigder*sigsq_om12 + else if (rij.lt.ssxm) then + havebond=.true. + ssd=rij-ssXs + eij=ssA*ssd*ssd+ssB*ssd+ssC + + ed=2*akcm*ssd+akct*deltat12 + pom1=akct*ssd + pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi + eom1=-2*akth*deltat1-pom1-om2*pom2 + eom2= 2*akth*deltat2+pom1-om1*pom2 + eom12=pom2 + else + omega=v1ss+2.0d0*v2ss*cosphi+3.0d0*v3ss*cosphi*cosphi + + d_ssxm(1)=0.5D0*akct/ssA + d_ssxm(2)=-d_ssxm(1) + d_ssxm(3)=0.0D0 + + d_ljxm(1)=sig0ij/sqrt(sigsq**3) + d_ljxm(2)=d_ljxm(1)*sigsq_om2 + d_ljxm(3)=d_ljxm(1)*sigsq_om12 + d_ljxm(1)=d_ljxm(1)*sigsq_om1 + +c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE + xm=0.5d0*(ssxm+ljxm) + do k=1,3 + d_xm(k)=0.5d0*(d_ssxm(k)+d_ljxm(k)) + enddo + if (rij.lt.xm) then + havebond=.true. + ssm=ssC-0.25D0*ssB*ssB/ssA + d_ssm(1)=0.5D0*akct*ssB/ssA + d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1) + d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1) + d_ssm(3)=omega + f1=(rij-xm)/(ssxm-xm) + f2=(rij-ssxm)/(xm-ssxm) + h1=h_base(f1,hd1) + h2=h_base(f2,hd2) + eij=ssm*h1+Ht*h2 + delta_inv=1.0d0/(xm-ssxm) + deltasq_inv=delta_inv*delta_inv + fac=ssm*hd1-Ht*hd2 + fac1=deltasq_inv*fac*(xm-rij) + fac2=deltasq_inv*fac*(rij-ssxm) + ed=delta_inv*(Ht*hd2-ssm*hd1) + eom1=fac1*d_ssxm(1)+fac2*d_xm(1)+h1*d_ssm(1) + eom2=fac1*d_ssxm(2)+fac2*d_xm(2)+h1*d_ssm(2) + eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3) + else + havebond=.false. + ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj) + d_ljm(1)=-0.5D0*bb(itypi,itypj)/aa(itypi,itypj)*ljB + d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt) + d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt- + + alf12/eps3rt) + d_ljm(1)=d_ljm(1)*(0.5D0*eps2rt_om1/eps2rt-alf1/eps3rt) + f1=(rij-ljxm)/(xm-ljxm) + f2=(rij-xm)/(ljxm-xm) + h1=h_base(f1,hd1) + h2=h_base(f2,hd2) + eij=Ht*h1+ljm*h2 + delta_inv=1.0d0/(ljxm-xm) + deltasq_inv=delta_inv*delta_inv + fac=Ht*hd1-ljm*hd2 + fac1=deltasq_inv*fac*(ljxm-rij) + fac2=deltasq_inv*fac*(rij-xm) + ed=delta_inv*(ljm*hd2-Ht*hd1) + eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1) + eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2) + eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3) + endif +c-------END FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE + +c-------SECOND METHOD, CONTINUOUS SECOND DERIVATIVE +c$$$ ssd=rij-ssXs +c$$$ ljd=rij-ljXs +c$$$ fac1=rij-ljxm +c$$$ fac2=rij-ssxm +c$$$ +c$$$ d_ljB(1)=ljB*(eps2rt_om1/eps2rt-2.0d0*alf1/eps3rt) +c$$$ d_ljB(2)=ljB*(eps2rt_om2/eps2rt+2.0d0*alf2/eps3rt) +c$$$ d_ljB(3)=ljB*(eps1_om12+eps2rt_om12/eps2rt-2.0d0*alf12/eps3rt) +c$$$ +c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA +c$$$ d_ssm(1)=0.5D0*akct*ssB/ssA +c$$$ d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1) +c$$$ d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1) +c$$$ d_ssm(3)=omega +c$$$ +c$$$ ljm=-0.25D0*bb(itypi,itypj)/aa(itypi,itypj) +c$$$ do k=1,3 +c$$$ d_ljm(k)=ljm*d_ljB(k) +c$$$ enddo +c$$$ ljm=ljm*ljB +c$$$ +c$$$ ss=ssA*ssd*ssd+ssB*ssd+ssC +c$$$ d_ss(0)=2.0d0*ssA*ssd+ssB +c$$$ d_ss(2)=akct*ssd +c$$$ d_ss(1)=-d_ss(2)-2.0d0*akth*deltat1-om2*omega +c$$$ d_ss(2)=d_ss(2)+2.0d0*akth*deltat2-om1*omega +c$$$ d_ss(3)=omega +c$$$ +c$$$ ljf=bb(itypi,itypj)/aa(itypi,itypj) +c$$$ ljf=9.0d0*ljf*(-0.5d0*ljf)**(1.0d0/3.0d0) +c$$$ d_ljf(0)=ljf*2.0d0*ljB*fac1 +c$$$ do k=1,3 +c$$$ d_ljf(k)=d_ljm(k)+ljf*(d_ljB(k)*fac1*fac1- +c$$$ & 2.0d0*ljB*fac1*d_ljxm(k)) +c$$$ enddo +c$$$ ljf=ljm+ljf*ljB*fac1*fac1 +c$$$ +c$$$ f1=(rij-ljxm)/(ssxm-ljxm) +c$$$ f2=(rij-ssxm)/(ljxm-ssxm) +c$$$ h1=h_base(f1,hd1) +c$$$ h2=h_base(f2,hd2) +c$$$ eij=ss*h1+ljf*h2 +c$$$ delta_inv=1.0d0/(ljxm-ssxm) +c$$$ deltasq_inv=delta_inv*delta_inv +c$$$ fac=ljf*hd2-ss*hd1 +c$$$ ed=d_ss(0)*h1+d_ljf(0)*h2+delta_inv*fac +c$$$ eom1=d_ss(1)*h1+d_ljf(1)*h2+deltasq_inv*fac* +c$$$ & (fac1*d_ssxm(1)-fac2*(d_ljxm(1))) +c$$$ eom2=d_ss(2)*h1+d_ljf(2)*h2+deltasq_inv*fac* +c$$$ & (fac1*d_ssxm(2)-fac2*(d_ljxm(2))) +c$$$ eom12=d_ss(3)*h1+d_ljf(3)*h2+deltasq_inv*fac* +c$$$ & (fac1*d_ssxm(3)-fac2*(d_ljxm(3))) +c$$$ +c$$$ havebond=.false. +c$$$ if (ed.gt.0.0d0) havebond=.true. +c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE + + endif + + if (havebond) then +#ifndef CLUST +#ifndef WHAM +c if (dyn_ssbond_ij(i,j).eq.1.0d300) then +c write(iout,'(a15,f12.2,f8.1,2i5)') +c & "SSBOND_E_FORM",totT,t_bath,i,j +c endif +#endif +#endif + dyn_ssbond_ij(i,j)=eij + else if (.not.havebond .and. dyn_ssbond_ij(i,j).lt.1.0d300) then + dyn_ssbond_ij(i,j)=1.0d300 +#ifndef CLUST +#ifndef WHAM +c write(iout,'(a15,f12.2,f8.1,2i5)') +c & "SSBOND_E_BREAK",totT,t_bath,i,j +#endif +#endif + endif + +c-------TESTING CODE + if (checkstop) then + if (jcheck.eq.0) write(iout,'(a,3f15.8,$)') + & "CHECKSTOP",rij,eij,ed + echeck(jcheck)=eij + endif + enddo + if (checkstop) then + write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps + endif + enddo + if (checkstop) then + transgrad=.true. + checkstop=.false. + endif +c-------END TESTING CODE + + do k=1,3 + dcosom1(k)=(dc_norm(k,nres+i)-om1*erij(k))/rij + dcosom2(k)=(dc_norm(k,nres+j)-om2*erij(k))/rij + enddo + do k=1,3 + gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k) + enddo + do k=1,3 + gvdwx(k,i)=gvdwx(k,i)-gg(k) + & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) + & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv + gvdwx(k,j)=gvdwx(k,j)+gg(k) + & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) + & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv + enddo +cgrad do k=i,j-1 +cgrad do l=1,3 +cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l) +cgrad enddo +cgrad enddo + + do l=1,3 + gvdwc(l,i)=gvdwc(l,i)-gg(l) + gvdwc(l,j)=gvdwc(l,j)+gg(l) + enddo + + return + end + +C----------------------------------------------------------------------------- + + double precision function h_base(x,deriv) +c A smooth function going 0->1 in range [0,1] +c It should NOT be called outside range [0,1], it will not work there. + implicit none + +c Input arguments + double precision x + +c Output arguments + double precision deriv + +c Local variables + double precision xsq + + +c Two parabolas put together. First derivative zero at extrema +c$$$ if (x.lt.0.5D0) then +c$$$ h_base=2.0D0*x*x +c$$$ deriv=4.0D0*x +c$$$ else +c$$$ deriv=1.0D0-x +c$$$ h_base=1.0D0-2.0D0*deriv*deriv +c$$$ deriv=4.0D0*deriv +c$$$ endif + +c Third degree polynomial. First derivative zero at extrema + h_base=x*x*(3.0d0-2.0d0*x) + deriv=6.0d0*x*(1.0d0-x) + +c Fifth degree polynomial. First and second derivatives zero at extrema +c$$$ xsq=x*x +c$$$ h_base=x*xsq*(6.0d0*xsq-15.0d0*x+10.0d0) +c$$$ deriv=x-1.0d0 +c$$$ deriv=deriv*deriv +c$$$ deriv=30.0d0*xsq*deriv + + return + end + +c---------------------------------------------------------------------------- + + subroutine dyn_set_nss +c Adjust nss and other relevant variables based on dyn_ssbond_ij +c implicit none + +c Includes + include 'DIMENSIONS' +#ifdef MPI + include "mpif.h" +#endif + include 'COMMON.SBRIDGE' + include 'COMMON.CHAIN' + include 'COMMON.IOUNITS' + include 'COMMON.SETUP' +#ifndef CLUST +#ifndef WHAM + include 'COMMON.MD' +#endif +#endif + +c Local variables + double precision emin + integer i,j,imin + integer diff,allflag(maxdim),allnss, + & allihpb(maxdim),alljhpb(maxdim), + & newnss,newihpb(maxdim),newjhpb(maxdim) + logical found + integer i_newnss(max_fg_procs),displ(0:max_fg_procs) + integer g_newihpb(maxdim),g_newjhpb(maxdim),g_newnss + + allnss=0 + do i=1,nres-1 + do j=i+1,nres + if (dyn_ssbond_ij(i,j).lt.1.0d300) then + allnss=allnss+1 + allflag(allnss)=0 + allihpb(allnss)=i + alljhpb(allnss)=j + endif + enddo + enddo + +cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss) + + 1 emin=1.0d300 + do i=1,allnss + if (allflag(i).eq.0 .and. + & dyn_ssbond_ij(allihpb(i),alljhpb(i)).lt.emin) then + emin=dyn_ssbond_ij(allihpb(i),alljhpb(i)) + imin=i + endif + enddo + if (emin.lt.1.0d300) then + allflag(imin)=1 + do i=1,allnss + if (allflag(i).eq.0 .and. + & (allihpb(i).eq.allihpb(imin) .or. + & alljhpb(i).eq.allihpb(imin) .or. + & allihpb(i).eq.alljhpb(imin) .or. + & alljhpb(i).eq.alljhpb(imin))) then + allflag(i)=-1 + endif + enddo + goto 1 + endif + +cmc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss) + + newnss=0 + do i=1,allnss + if (allflag(i).eq.1) then + newnss=newnss+1 + newihpb(newnss)=allihpb(i) + newjhpb(newnss)=alljhpb(i) + endif + enddo + +#ifdef MPI + if (nfgtasks.gt.1)then + + call MPI_Reduce(newnss,g_newnss,1, + & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR) + call MPI_Gather(newnss,1,MPI_INTEGER, + & i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR) + displ(0)=0 + do i=1,nfgtasks-1,1 + displ(i)=i_newnss(i-1)+displ(i-1) + enddo + call MPI_Gatherv(newihpb,newnss,MPI_INTEGER, + & g_newihpb,i_newnss,displ,MPI_INTEGER, + & king,FG_COMM,IERR) + call MPI_Gatherv(newjhpb,newnss,MPI_INTEGER, + & g_newjhpb,i_newnss,displ,MPI_INTEGER, + & king,FG_COMM,IERR) + if(fg_rank.eq.0) then +c print *,'g_newnss',g_newnss +c print *,'g_newihpb',(g_newihpb(i),i=1,g_newnss) +c print *,'g_newjhpb',(g_newjhpb(i),i=1,g_newnss) + newnss=g_newnss + do i=1,newnss + newihpb(i)=g_newihpb(i) + newjhpb(i)=g_newjhpb(i) + enddo + endif + endif +#endif + + diff=newnss-nss + +cmc write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss) + + do i=1,nss + found=.false. + do j=1,newnss + if (idssb(i).eq.newihpb(j) .and. + & jdssb(i).eq.newjhpb(j)) found=.true. + enddo +#ifndef CLUST +#ifndef WHAM + if (.not.found.and.fg_rank.eq.0) + & write(iout,'(a15,f12.2,f8.1,2i5)') + & "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i) +#endif +#endif + enddo + + do i=1,newnss + found=.false. + do j=1,nss + if (newihpb(i).eq.idssb(j) .and. + & newjhpb(i).eq.jdssb(j)) found=.true. + enddo +#ifndef CLUST +#ifndef WHAM + if (.not.found.and.fg_rank.eq.0) + & write(iout,'(a15,f12.2,f8.1,2i5)') + & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i) +#endif +#endif + enddo + + nss=newnss + do i=1,nss + idssb(i)=newihpb(i) + jdssb(i)=newjhpb(i) + enddo + + return + end + +c---------------------------------------------------------------------------- + +#ifdef WHAM + subroutine read_ssHist + implicit none + +c Includes + include 'DIMENSIONS' + include "DIMENSIONS.FREE" + include 'COMMON.FREE' + +c Local variables + integer i,j + character*80 controlcard + + do i=1,dyn_nssHist + call card_concat(controlcard,.true.) + read(controlcard,*) + & dyn_ssHist(i,0),(dyn_ssHist(i,j),j=1,2*dyn_ssHist(i,0)) + enddo + + return + end +#endif + +c---------------------------------------------------------------------------- + + +C----------------------------------------------------------------------------- +C----------------------------------------------------------------------------- +C----------------------------------------------------------------------------- +C----------------------------------------------------------------------------- +C----------------------------------------------------------------------------- +C----------------------------------------------------------------------------- +C----------------------------------------------------------------------------- + +c$$$c----------------------------------------------------------------------------- +c$$$ +c$$$ subroutine ss_relax(i_in,j_in) +c$$$ implicit none +c$$$ +c$$$c Includes +c$$$ include 'DIMENSIONS' +c$$$ include 'COMMON.VAR' +c$$$ include 'COMMON.CHAIN' +c$$$ include 'COMMON.IOUNITS' +c$$$ include 'COMMON.INTERACT' +c$$$ +c$$$c Input arguments +c$$$ integer i_in,j_in +c$$$ +c$$$c Local variables +c$$$ integer i,iretcode,nfun_sc +c$$$ logical scfail +c$$$ double precision var(maxvar),e_sc,etot +c$$$ +c$$$ +c$$$ mask_r=.true. +c$$$ do i=nnt,nct +c$$$ mask_side(i)=0 +c$$$ enddo +c$$$ mask_side(i_in)=1 +c$$$ mask_side(j_in)=1 +c$$$ +c$$$c Minimize the two selected side-chains +c$$$ call overlap_sc(scfail) ! Better not fail! +c$$$ call minimize_sc(e_sc,var,iretcode,nfun_sc) +c$$$ +c$$$ mask_r=.false. +c$$$ +c$$$ return +c$$$ end +c$$$ +c$$$c------------------------------------------------------------- +c$$$ +c$$$ subroutine minimize_sc(etot_sc,iretcode,nfun) +c$$$c Minimize side-chains only, starting from geom but without modifying +c$$$c bond lengths. +c$$$c If mask_r is already set, only the selected side-chains are minimized, +c$$$c otherwise all side-chains are minimized keeping the backbone frozen. +c$$$ implicit none +c$$$ +c$$$c Includes +c$$$ include 'DIMENSIONS' +c$$$ include 'COMMON.IOUNITS' +c$$$ include 'COMMON.VAR' +c$$$ include 'COMMON.CHAIN' +c$$$ include 'COMMON.GEO' +c$$$ include 'COMMON.MINIM' +c$$$ integer icall +c$$$ common /srutu/ icall +c$$$ +c$$$c Output arguments +c$$$ double precision etot_sc +c$$$ integer iretcode,nfun +c$$$ +c$$$c External functions/subroutines +c$$$ external func_sc,grad_sc,fdum +c$$$ +c$$$c Local variables +c$$$ integer liv,lv +c$$$ parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2)) +c$$$ integer iv(liv) +c$$$ double precision rdum(1) +c$$$ double precision d(maxvar),v(1:lv),x(maxvar),xx(maxvar) +c$$$ integer idum(1) +c$$$ integer i,nvar_restr +c$$$ +c$$$ +c$$$cmc start_minim=.true. +c$$$ call deflt(2,iv,liv,lv,v) +c$$$* 12 means fresh start, dont call deflt +c$$$ iv(1)=12 +c$$$* max num of fun calls +c$$$ if (maxfun.eq.0) maxfun=500 +c$$$ iv(17)=maxfun +c$$$* max num of iterations +c$$$ if (maxmin.eq.0) maxmin=1000 +c$$$ iv(18)=maxmin +c$$$* controls output +c$$$ iv(19)=1 +c$$$* selects output unit +c$$$ iv(21)=0 +c$$$c iv(21)=iout ! DEBUG +c$$$c iv(21)=8 ! DEBUG +c$$$* 1 means to print out result +c$$$ iv(22)=0 +c$$$c iv(22)=1 ! DEBUG +c$$$* 1 means to print out summary stats +c$$$ iv(23)=0 +c$$$c iv(23)=1 ! DEBUG +c$$$* 1 means to print initial x and d +c$$$ iv(24)=0 +c$$$c iv(24)=1 ! DEBUG +c$$$* min val for v(radfac) default is 0.1 +c$$$ v(24)=0.1D0 +c$$$* max val for v(radfac) default is 4.0 +c$$$ v(25)=2.0D0 +c$$$c v(25)=4.0D0 +c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease) +c$$$* the sumsl default is 0.1 +c$$$ v(26)=0.1D0 +c$$$* false conv if (act fnctn decrease) .lt. v(34) +c$$$* the sumsl default is 100*machep +c$$$ v(34)=v(34)/100.0D0 +c$$$* absolute convergence +c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4 +c$$$ v(31)=tolf +c$$$* relative convergence +c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-1 +c$$$ v(32)=rtolf +c$$$* controls initial step size +c$$$ v(35)=1.0D-1 +c$$$* large vals of d correspond to small components of step +c$$$ do i=1,nphi +c$$$ d(i)=1.0D-1 +c$$$ enddo +c$$$ do i=nphi+1,nvar +c$$$ d(i)=1.0D-1 +c$$$ enddo +c$$$ +c$$$ call geom_to_var(nvar,x) +c$$$ IF (mask_r) THEN +c$$$ do i=1,nres ! Just in case... +c$$$ mask_phi(i)=0 +c$$$ mask_theta(i)=0 +c$$$ enddo +c$$$ call x2xx(x,xx,nvar_restr) +c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc, +c$$$ & iv,liv,lv,v,idum,rdum,fdum) +c$$$ call xx2x(x,xx) +c$$$ ELSE +c$$$c When minimizing ALL side-chains, etotal_sc is a little +c$$$c faster if we don't set mask_r +c$$$ do i=1,nres +c$$$ mask_phi(i)=0 +c$$$ mask_theta(i)=0 +c$$$ mask_side(i)=1 +c$$$ enddo +c$$$ call x2xx(x,xx,nvar_restr) +c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc, +c$$$ & iv,liv,lv,v,idum,rdum,fdum) +c$$$ call xx2x(x,xx) +c$$$ ENDIF +c$$$ call var_to_geom(nvar,x) +c$$$ call chainbuild_sc +c$$$ etot_sc=v(10) +c$$$ iretcode=iv(1) +c$$$ nfun=iv(6) +c$$$ return +c$$$ end +c$$$ +c$$$C-------------------------------------------------------------------------- +c$$$ +c$$$ subroutine chainbuild_sc +c$$$ implicit none +c$$$ include 'DIMENSIONS' +c$$$ include 'COMMON.VAR' +c$$$ include 'COMMON.INTERACT' +c$$$ +c$$$c Local variables +c$$$ integer i +c$$$ +c$$$ +c$$$ do i=nnt,nct +c$$$ if (.not.mask_r .or. mask_side(i).eq.1) then +c$$$ call locate_side_chain(i) +c$$$ endif +c$$$ enddo +c$$$ +c$$$ return +c$$$ end +c$$$ +c$$$C-------------------------------------------------------------------------- +c$$$ +c$$$ subroutine func_sc(n,x,nf,f,uiparm,urparm,ufparm) +c$$$ implicit none +c$$$ +c$$$c Includes +c$$$ include 'DIMENSIONS' +c$$$ include 'COMMON.DERIV' +c$$$ include 'COMMON.VAR' +c$$$ include 'COMMON.MINIM' +c$$$ include 'COMMON.IOUNITS' +c$$$ +c$$$c Input arguments +c$$$ integer n +c$$$ double precision x(maxvar) +c$$$ double precision ufparm +c$$$ external ufparm +c$$$ +c$$$c Input/Output arguments +c$$$ integer nf +c$$$ integer uiparm(1) +c$$$ double precision urparm(1) +c$$$ +c$$$c Output arguments +c$$$ double precision f +c$$$ +c$$$c Local variables +c$$$ double precision energia(0:n_ene) +c$$$#ifdef OSF +c$$$c Variables used to intercept NaNs +c$$$ double precision x_sum +c$$$ integer i_NAN +c$$$#endif +c$$$ +c$$$ +c$$$ nfl=nf +c$$$ icg=mod(nf,2)+1 +c$$$ +c$$$#ifdef OSF +c$$$c Intercept NaNs in the coordinates, before calling etotal_sc +c$$$ x_sum=0.D0 +c$$$ do i_NAN=1,n +c$$$ x_sum=x_sum+x(i_NAN) +c$$$ enddo +c$$$c Calculate the energy only if the coordinates are ok +c$$$ if ((.not.(x_sum.lt.0.D0)) .and. (.not.(x_sum.ge.0.D0))) then +c$$$ write(iout,*)" *** func_restr_sc : Found NaN in coordinates" +c$$$ f=1.0D+77 +c$$$ nf=0 +c$$$ else +c$$$#endif +c$$$ +c$$$ call var_to_geom_restr(n,x) +c$$$ call zerograd +c$$$ call chainbuild_sc +c$$$ call etotal_sc(energia(0)) +c$$$ f=energia(0) +c$$$ if (energia(1).eq.1.0D20 .or. energia(0).eq.1.0D99) nf=0 +c$$$ +c$$$#ifdef OSF +c$$$ endif +c$$$#endif +c$$$ +c$$$ return +c$$$ end +c$$$ +c$$$c------------------------------------------------------- +c$$$ +c$$$ subroutine grad_sc(n,x,nf,g,uiparm,urparm,ufparm) +c$$$ implicit none +c$$$ +c$$$c Includes +c$$$ include 'DIMENSIONS' +c$$$ include 'COMMON.CHAIN' +c$$$ include 'COMMON.DERIV' +c$$$ include 'COMMON.VAR' +c$$$ include 'COMMON.INTERACT' +c$$$ include 'COMMON.MINIM' +c$$$ +c$$$c Input arguments +c$$$ integer n +c$$$ double precision x(maxvar) +c$$$ double precision ufparm +c$$$ external ufparm +c$$$ +c$$$c Input/Output arguments +c$$$ integer nf +c$$$ integer uiparm(1) +c$$$ double precision urparm(1) +c$$$ +c$$$c Output arguments +c$$$ double precision g(maxvar) +c$$$ +c$$$c Local variables +c$$$ double precision f,gphii,gthetai,galphai,gomegai +c$$$ integer ig,ind,i,j,k,igall,ij +c$$$ +c$$$ +c$$$ icg=mod(nf,2)+1 +c$$$ if (nf-nfl+1) 20,30,40 +c$$$ 20 call func_sc(n,x,nf,f,uiparm,urparm,ufparm) +c$$$c write (iout,*) 'grad 20' +c$$$ if (nf.eq.0) return +c$$$ goto 40 +c$$$ 30 call var_to_geom_restr(n,x) +c$$$ call chainbuild_sc +c$$$C +c$$$C Evaluate the derivatives of virtual bond lengths and SC vectors in variables. +c$$$C +c$$$ 40 call cartder +c$$$C +c$$$C Convert the Cartesian gradient into internal-coordinate gradient. +c$$$C +c$$$ +c$$$ ig=0 +c$$$ ind=nres-2 +c$$$ do i=2,nres-2 +c$$$ IF (mask_phi(i+2).eq.1) THEN +c$$$ gphii=0.0D0 +c$$$ do j=i+1,nres-1 +c$$$ ind=ind+1 +c$$$ do k=1,3 +c$$$ gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg) +c$$$ gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg) +c$$$ enddo +c$$$ enddo +c$$$ ig=ig+1 +c$$$ g(ig)=gphii +c$$$ ELSE +c$$$ ind=ind+nres-1-i +c$$$ ENDIF +c$$$ enddo +c$$$ +c$$$ +c$$$ ind=0 +c$$$ do i=1,nres-2 +c$$$ IF (mask_theta(i+2).eq.1) THEN +c$$$ ig=ig+1 +c$$$ gthetai=0.0D0 +c$$$ do j=i+1,nres-1 +c$$$ ind=ind+1 +c$$$ do k=1,3 +c$$$ gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg) +c$$$ gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg) +c$$$ enddo +c$$$ enddo +c$$$ g(ig)=gthetai +c$$$ ELSE +c$$$ ind=ind+nres-1-i +c$$$ ENDIF +c$$$ enddo +c$$$ +c$$$ do i=2,nres-1 +c$$$ if (itype(i).ne.10) then +c$$$ IF (mask_side(i).eq.1) THEN +c$$$ ig=ig+1 +c$$$ galphai=0.0D0 +c$$$ do k=1,3 +c$$$ galphai=galphai+dxds(k,i)*gradx(k,i,icg) +c$$$ enddo +c$$$ g(ig)=galphai +c$$$ ENDIF +c$$$ endif +c$$$ enddo +c$$$ +c$$$ +c$$$ do i=2,nres-1 +c$$$ if (itype(i).ne.10) then +c$$$ IF (mask_side(i).eq.1) THEN +c$$$ ig=ig+1 +c$$$ gomegai=0.0D0 +c$$$ do k=1,3 +c$$$ gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg) +c$$$ enddo +c$$$ g(ig)=gomegai +c$$$ ENDIF +c$$$ endif +c$$$ enddo +c$$$ +c$$$C +c$$$C Add the components corresponding to local energy terms. +c$$$C +c$$$ +c$$$ ig=0 +c$$$ igall=0 +c$$$ do i=4,nres +c$$$ igall=igall+1 +c$$$ if (mask_phi(i).eq.1) then +c$$$ ig=ig+1 +c$$$ g(ig)=g(ig)+gloc(igall,icg) +c$$$ endif +c$$$ enddo +c$$$ +c$$$ do i=3,nres +c$$$ igall=igall+1 +c$$$ if (mask_theta(i).eq.1) then +c$$$ ig=ig+1 +c$$$ g(ig)=g(ig)+gloc(igall,icg) +c$$$ endif +c$$$ enddo +c$$$ +c$$$ do ij=1,2 +c$$$ do i=2,nres-1 +c$$$ if (itype(i).ne.10) then +c$$$ igall=igall+1 +c$$$ if (mask_side(i).eq.1) then +c$$$ ig=ig+1 +c$$$ g(ig)=g(ig)+gloc(igall,icg) +c$$$ endif +c$$$ endif +c$$$ enddo +c$$$ enddo +c$$$ +c$$$cd do i=1,ig +c$$$cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i) +c$$$cd enddo +c$$$ +c$$$ return +c$$$ end +c$$$ +c$$$C----------------------------------------------------------------------------- +c$$$ +c$$$ subroutine etotal_sc(energy_sc) +c$$$ implicit none +c$$$ +c$$$c Includes +c$$$ include 'DIMENSIONS' +c$$$ include 'COMMON.VAR' +c$$$ include 'COMMON.INTERACT' +c$$$ include 'COMMON.DERIV' +c$$$ include 'COMMON.FFIELD' +c$$$ +c$$$c Output arguments +c$$$ double precision energy_sc(0:n_ene) +c$$$ +c$$$c Local variables +c$$$ double precision evdw,escloc +c$$$ integer i,j +c$$$ +c$$$ +c$$$ do i=1,n_ene +c$$$ energy_sc(i)=0.0D0 +c$$$ enddo +c$$$ +c$$$ if (mask_r) then +c$$$ call egb_sc(evdw) +c$$$ call esc_sc(escloc) +c$$$ else +c$$$ call egb(evdw) +c$$$ call esc(escloc) +c$$$ endif +c$$$ +c$$$ if (evdw.eq.1.0D20) then +c$$$ energy_sc(0)=evdw +c$$$ else +c$$$ energy_sc(0)=wsc*evdw+wscloc*escloc +c$$$ endif +c$$$ energy_sc(1)=evdw +c$$$ energy_sc(12)=escloc +c$$$ +c$$$C +c$$$C Sum up the components of the Cartesian gradient. +c$$$C +c$$$ do i=1,nct +c$$$ do j=1,3 +c$$$ gradx(j,i,icg)=wsc*gvdwx(j,i) +c$$$ enddo +c$$$ enddo +c$$$ +c$$$ return +c$$$ end +c$$$ +c$$$C----------------------------------------------------------------------------- +c$$$ +c$$$ subroutine egb_sc(evdw) +c$$$C +c$$$C This subroutine calculates the interaction energy of nonbonded side chains +c$$$C assuming the Gay-Berne potential of interaction. +c$$$C +c$$$ implicit real*8 (a-h,o-z) +c$$$ include 'DIMENSIONS' +c$$$ include 'COMMON.GEO' +c$$$ include 'COMMON.VAR' +c$$$ include 'COMMON.LOCAL' +c$$$ include 'COMMON.CHAIN' +c$$$ include 'COMMON.DERIV' +c$$$ include 'COMMON.NAMES' +c$$$ include 'COMMON.INTERACT' +c$$$ include 'COMMON.IOUNITS' +c$$$ include 'COMMON.CALC' +c$$$ include 'COMMON.CONTROL' +c$$$ logical lprn +c$$$ evdw=0.0D0 +c$$$ energy_dec=.false. +c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon +c$$$ evdw=0.0D0 +c$$$ lprn=.false. +c$$$c if (icall.eq.0) lprn=.false. +c$$$ ind=0 +c$$$ do i=iatsc_s,iatsc_e +c$$$ itypi=itype(i) +c$$$ itypi1=itype(i+1) +c$$$ xi=c(1,nres+i) +c$$$ yi=c(2,nres+i) +c$$$ zi=c(3,nres+i) +c$$$ dxi=dc_norm(1,nres+i) +c$$$ dyi=dc_norm(2,nres+i) +c$$$ dzi=dc_norm(3,nres+i) +c$$$c dsci_inv=dsc_inv(itypi) +c$$$ dsci_inv=vbld_inv(i+nres) +c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres) +c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi +c$$$C +c$$$C Calculate SC interaction energy. +c$$$C +c$$$ do iint=1,nint_gr(i) +c$$$ do j=istart(i,iint),iend(i,iint) +c$$$ IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN +c$$$ ind=ind+1 +c$$$ itypj=itype(j) +c$$$c dscj_inv=dsc_inv(itypj) +c$$$ dscj_inv=vbld_inv(j+nres) +c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv, +c$$$c & 1.0d0/vbld(j+nres) +c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) +c$$$ sig0ij=sigma(itypi,itypj) +c$$$ chi1=chi(itypi,itypj) +c$$$ chi2=chi(itypj,itypi) +c$$$ chi12=chi1*chi2 +c$$$ chip1=chip(itypi) +c$$$ chip2=chip(itypj) +c$$$ chip12=chip1*chip2 +c$$$ alf1=alp(itypi) +c$$$ alf2=alp(itypj) +c$$$ alf12=0.5D0*(alf1+alf2) +c$$$C For diagnostics only!!! +c$$$c chi1=0.0D0 +c$$$c chi2=0.0D0 +c$$$c chi12=0.0D0 +c$$$c chip1=0.0D0 +c$$$c chip2=0.0D0 +c$$$c chip12=0.0D0 +c$$$c alf1=0.0D0 +c$$$c alf2=0.0D0 +c$$$c alf12=0.0D0 +c$$$ xj=c(1,nres+j)-xi +c$$$ yj=c(2,nres+j)-yi +c$$$ zj=c(3,nres+j)-zi +c$$$ dxj=dc_norm(1,nres+j) +c$$$ dyj=dc_norm(2,nres+j) +c$$$ dzj=dc_norm(3,nres+j) +c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi +c$$$c write (iout,*) "j",j," dc_norm", +c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j) +c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj) +c$$$ rij=dsqrt(rrij) +c$$$C Calculate angle-dependent terms of energy and contributions to their +c$$$C derivatives. +c$$$ call sc_angular +c$$$ sigsq=1.0D0/sigsq +c$$$ sig=sig0ij*dsqrt(sigsq) +c$$$ rij_shift=1.0D0/rij-sig+sig0ij +c$$$c for diagnostics; uncomment +c$$$c rij_shift=1.2*sig0ij +c$$$C I hate to put IF's in the loops, but here don't have another choice!!!! +c$$$ if (rij_shift.le.0.0D0) then +c$$$ evdw=1.0D20 +c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))') +c$$$cd & restyp(itypi),i,restyp(itypj),j, +c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) +c$$$ return +c$$$ endif +c$$$ sigder=-sig*sigsq +c$$$c--------------------------------------------------------------- +c$$$ rij_shift=1.0D0/rij_shift +c$$$ fac=rij_shift**expon +c$$$ e1=fac*fac*aa(itypi,itypj) +c$$$ e2=fac*bb(itypi,itypj) +c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2) +c$$$ eps2der=evdwij*eps3rt +c$$$ eps3der=evdwij*eps2rt +c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt, +c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 +c$$$ evdwij=evdwij*eps2rt*eps3rt +c$$$ evdw=evdw+evdwij +c$$$ if (lprn) then +c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) +c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj) +c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))') +c$$$ & restyp(itypi),i,restyp(itypj),j, +c$$$ & epsi,sigm,chi1,chi2,chip1,chip2, +c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij, +c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift, +c$$$ & evdwij +c$$$ endif +c$$$ +c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)') +c$$$ & 'evdw',i,j,evdwij +c$$$ +c$$$C Calculate gradient components. +c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2 +c$$$ fac=-expon*(e1+evdwij)*rij_shift +c$$$ sigder=fac*sigder +c$$$ fac=rij*fac +c$$$c fac=0.0d0 +c$$$C Calculate the radial part of the gradient +c$$$ gg(1)=xj*fac +c$$$ gg(2)=yj*fac +c$$$ gg(3)=zj*fac +c$$$C Calculate angular part of the gradient. +c$$$ call sc_grad +c$$$ ENDIF +c$$$ enddo ! j +c$$$ enddo ! iint +c$$$ enddo ! i +c$$$ energy_dec=.false. +c$$$ return +c$$$ end +c$$$ +c$$$c----------------------------------------------------------------------------- +c$$$ +c$$$ subroutine esc_sc(escloc) +c$$$C Calculate the local energy of a side chain and its derivatives in the +c$$$C corresponding virtual-bond valence angles THETA and the spherical angles +c$$$C ALPHA and OMEGA. +c$$$ implicit real*8 (a-h,o-z) +c$$$ include 'DIMENSIONS' +c$$$ include 'COMMON.GEO' +c$$$ include 'COMMON.LOCAL' +c$$$ include 'COMMON.VAR' +c$$$ include 'COMMON.INTERACT' +c$$$ include 'COMMON.DERIV' +c$$$ include 'COMMON.CHAIN' +c$$$ include 'COMMON.IOUNITS' +c$$$ include 'COMMON.NAMES' +c$$$ include 'COMMON.FFIELD' +c$$$ include 'COMMON.CONTROL' +c$$$ double precision x(3),dersc(3),xemp(3),dersc0(3),dersc1(3), +c$$$ & ddersc0(3),ddummy(3),xtemp(3),temp(3) +c$$$ common /sccalc/ time11,time12,time112,theti,it,nlobit +c$$$ delta=0.02d0*pi +c$$$ escloc=0.0D0 +c$$$c write (iout,'(a)') 'ESC' +c$$$ do i=loc_start,loc_end +c$$$ IF (mask_side(i).eq.1) THEN +c$$$ it=itype(i) +c$$$ if (it.eq.10) goto 1 +c$$$ nlobit=nlob(it) +c$$$c print *,'i=',i,' it=',it,' nlobit=',nlobit +c$$$c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad +c$$$ theti=theta(i+1)-pipol +c$$$ x(1)=dtan(theti) +c$$$ x(2)=alph(i) +c$$$ x(3)=omeg(i) +c$$$ +c$$$ if (x(2).gt.pi-delta) then +c$$$ xtemp(1)=x(1) +c$$$ xtemp(2)=pi-delta +c$$$ xtemp(3)=x(3) +c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.) +c$$$ xtemp(2)=pi +c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.) +c$$$ call spline1(x(2),pi-delta,delta,escloci0,escloci1,dersc0(2), +c$$$ & escloci,dersc(2)) +c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1), +c$$$ & ddersc0(1),dersc(1)) +c$$$ call spline2(x(2),pi-delta,delta,dersc0(3),dersc1(3), +c$$$ & ddersc0(3),dersc(3)) +c$$$ xtemp(2)=pi-delta +c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.) +c$$$ xtemp(2)=pi +c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.) +c$$$ call spline1(x(2),pi-delta,delta,esclocbi0,esclocbi1, +c$$$ & dersc0(2),esclocbi,dersc02) +c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1), +c$$$ & dersc12,dersc01) +c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd) +c$$$ dersc0(1)=dersc01 +c$$$ dersc0(2)=dersc02 +c$$$ dersc0(3)=0.0d0 +c$$$ do k=1,3 +c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k) +c$$$ enddo +c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi) +c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci, +c$$$c & esclocbi,ss,ssd +c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi +c$$$c escloci=esclocbi +c$$$c write (iout,*) escloci +c$$$ else if (x(2).lt.delta) then +c$$$ xtemp(1)=x(1) +c$$$ xtemp(2)=delta +c$$$ xtemp(3)=x(3) +c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.) +c$$$ xtemp(2)=0.0d0 +c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.) +c$$$ call spline1(x(2),delta,-delta,escloci0,escloci1,dersc0(2), +c$$$ & escloci,dersc(2)) +c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1), +c$$$ & ddersc0(1),dersc(1)) +c$$$ call spline2(x(2),delta,-delta,dersc0(3),dersc1(3), +c$$$ & ddersc0(3),dersc(3)) +c$$$ xtemp(2)=delta +c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.) +c$$$ xtemp(2)=0.0d0 +c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.) +c$$$ call spline1(x(2),delta,-delta,esclocbi0,esclocbi1, +c$$$ & dersc0(2),esclocbi,dersc02) +c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1), +c$$$ & dersc12,dersc01) +c$$$ dersc0(1)=dersc01 +c$$$ dersc0(2)=dersc02 +c$$$ dersc0(3)=0.0d0 +c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd) +c$$$ do k=1,3 +c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k) +c$$$ enddo +c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi) +c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci, +c$$$c & esclocbi,ss,ssd +c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi +c$$$c write (iout,*) escloci +c$$$ else +c$$$ call enesc(x,escloci,dersc,ddummy,.false.) +c$$$ endif +c$$$ +c$$$ escloc=escloc+escloci +c$$$ if (energy_dec) write (iout,'(a6,i,0pf7.3)') +c$$$ & 'escloc',i,escloci +c$$$c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc +c$$$ +c$$$ gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+ +c$$$ & wscloc*dersc(1) +c$$$ gloc(ialph(i,1),icg)=wscloc*dersc(2) +c$$$ gloc(ialph(i,1)+nside,icg)=wscloc*dersc(3) +c$$$ 1 continue +c$$$ ENDIF +c$$$ enddo +c$$$ return +c$$$ end +c$$$ +c$$$C----------------------------------------------------------------------------- +c$$$ +c$$$ subroutine egb_ij(i_sc,j_sc,evdw) +c$$$C +c$$$C This subroutine calculates the interaction energy of nonbonded side chains +c$$$C assuming the Gay-Berne potential of interaction. +c$$$C +c$$$ implicit real*8 (a-h,o-z) +c$$$ include 'DIMENSIONS' +c$$$ include 'COMMON.GEO' +c$$$ include 'COMMON.VAR' +c$$$ include 'COMMON.LOCAL' +c$$$ include 'COMMON.CHAIN' +c$$$ include 'COMMON.DERIV' +c$$$ include 'COMMON.NAMES' +c$$$ include 'COMMON.INTERACT' +c$$$ include 'COMMON.IOUNITS' +c$$$ include 'COMMON.CALC' +c$$$ include 'COMMON.CONTROL' +c$$$ logical lprn +c$$$ evdw=0.0D0 +c$$$ energy_dec=.false. +c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon +c$$$ evdw=0.0D0 +c$$$ lprn=.false. +c$$$ ind=0 +c$$$c$$$ do i=iatsc_s,iatsc_e +c$$$ i=i_sc +c$$$ itypi=itype(i) +c$$$ itypi1=itype(i+1) +c$$$ xi=c(1,nres+i) +c$$$ yi=c(2,nres+i) +c$$$ zi=c(3,nres+i) +c$$$ dxi=dc_norm(1,nres+i) +c$$$ dyi=dc_norm(2,nres+i) +c$$$ dzi=dc_norm(3,nres+i) +c$$$c dsci_inv=dsc_inv(itypi) +c$$$ dsci_inv=vbld_inv(i+nres) +c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres) +c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi +c$$$C +c$$$C Calculate SC interaction energy. +c$$$C +c$$$c$$$ do iint=1,nint_gr(i) +c$$$c$$$ do j=istart(i,iint),iend(i,iint) +c$$$ j=j_sc +c$$$ ind=ind+1 +c$$$ itypj=itype(j) +c$$$c dscj_inv=dsc_inv(itypj) +c$$$ dscj_inv=vbld_inv(j+nres) +c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv, +c$$$c & 1.0d0/vbld(j+nres) +c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) +c$$$ sig0ij=sigma(itypi,itypj) +c$$$ chi1=chi(itypi,itypj) +c$$$ chi2=chi(itypj,itypi) +c$$$ chi12=chi1*chi2 +c$$$ chip1=chip(itypi) +c$$$ chip2=chip(itypj) +c$$$ chip12=chip1*chip2 +c$$$ alf1=alp(itypi) +c$$$ alf2=alp(itypj) +c$$$ alf12=0.5D0*(alf1+alf2) +c$$$C For diagnostics only!!! +c$$$c chi1=0.0D0 +c$$$c chi2=0.0D0 +c$$$c chi12=0.0D0 +c$$$c chip1=0.0D0 +c$$$c chip2=0.0D0 +c$$$c chip12=0.0D0 +c$$$c alf1=0.0D0 +c$$$c alf2=0.0D0 +c$$$c alf12=0.0D0 +c$$$ xj=c(1,nres+j)-xi +c$$$ yj=c(2,nres+j)-yi +c$$$ zj=c(3,nres+j)-zi +c$$$ dxj=dc_norm(1,nres+j) +c$$$ dyj=dc_norm(2,nres+j) +c$$$ dzj=dc_norm(3,nres+j) +c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi +c$$$c write (iout,*) "j",j," dc_norm", +c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j) +c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj) +c$$$ rij=dsqrt(rrij) +c$$$C Calculate angle-dependent terms of energy and contributions to their +c$$$C derivatives. +c$$$ call sc_angular +c$$$ sigsq=1.0D0/sigsq +c$$$ sig=sig0ij*dsqrt(sigsq) +c$$$ rij_shift=1.0D0/rij-sig+sig0ij +c$$$c for diagnostics; uncomment +c$$$c rij_shift=1.2*sig0ij +c$$$C I hate to put IF's in the loops, but here don't have another choice!!!! +c$$$ if (rij_shift.le.0.0D0) then +c$$$ evdw=1.0D20 +c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))') +c$$$cd & restyp(itypi),i,restyp(itypj),j, +c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) +c$$$ return +c$$$ endif +c$$$ sigder=-sig*sigsq +c$$$c--------------------------------------------------------------- +c$$$ rij_shift=1.0D0/rij_shift +c$$$ fac=rij_shift**expon +c$$$ e1=fac*fac*aa(itypi,itypj) +c$$$ e2=fac*bb(itypi,itypj) +c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2) +c$$$ eps2der=evdwij*eps3rt +c$$$ eps3der=evdwij*eps2rt +c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt, +c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 +c$$$ evdwij=evdwij*eps2rt*eps3rt +c$$$ evdw=evdw+evdwij +c$$$ if (lprn) then +c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) +c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj) +c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))') +c$$$ & restyp(itypi),i,restyp(itypj),j, +c$$$ & epsi,sigm,chi1,chi2,chip1,chip2, +c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij, +c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift, +c$$$ & evdwij +c$$$ endif +c$$$ +c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)') +c$$$ & 'evdw',i,j,evdwij +c$$$ +c$$$C Calculate gradient components. +c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2 +c$$$ fac=-expon*(e1+evdwij)*rij_shift +c$$$ sigder=fac*sigder +c$$$ fac=rij*fac +c$$$c fac=0.0d0 +c$$$C Calculate the radial part of the gradient +c$$$ gg(1)=xj*fac +c$$$ gg(2)=yj*fac +c$$$ gg(3)=zj*fac +c$$$C Calculate angular part of the gradient. +c$$$ call sc_grad +c$$$c$$$ enddo ! j +c$$$c$$$ enddo ! iint +c$$$c$$$ enddo ! i +c$$$ energy_dec=.false. +c$$$ return +c$$$ end +c$$$ +c$$$C----------------------------------------------------------------------------- +c$$$ +c$$$ subroutine perturb_side_chain(i,angle) +c$$$ implicit none +c$$$ +c$$$c Includes +c$$$ include 'DIMENSIONS' +c$$$ include 'COMMON.CHAIN' +c$$$ include 'COMMON.GEO' +c$$$ include 'COMMON.VAR' +c$$$ include 'COMMON.LOCAL' +c$$$ include 'COMMON.IOUNITS' +c$$$ +c$$$c External functions +c$$$ external ran_number +c$$$ double precision ran_number +c$$$ +c$$$c Input arguments +c$$$ integer i +c$$$ double precision angle ! In degrees +c$$$ +c$$$c Local variables +c$$$ integer i_sc +c$$$ double precision rad_ang,rand_v(3),length,cost,sint +c$$$ +c$$$ +c$$$ i_sc=i+nres +c$$$ rad_ang=angle*deg2rad +c$$$ +c$$$ length=0.0 +c$$$ do while (length.lt.0.01) +c$$$ rand_v(1)=ran_number(0.01D0,1.0D0) +c$$$ rand_v(2)=ran_number(0.01D0,1.0D0) +c$$$ rand_v(3)=ran_number(0.01D0,1.0D0) +c$$$ length=rand_v(1)*rand_v(1)+rand_v(2)*rand_v(2)+ +c$$$ + rand_v(3)*rand_v(3) +c$$$ length=sqrt(length) +c$$$ rand_v(1)=rand_v(1)/length +c$$$ rand_v(2)=rand_v(2)/length +c$$$ rand_v(3)=rand_v(3)/length +c$$$ cost=rand_v(1)*dc_norm(1,i_sc)+rand_v(2)*dc_norm(2,i_sc)+ +c$$$ + rand_v(3)*dc_norm(3,i_sc) +c$$$ length=1.0D0-cost*cost +c$$$ if (length.lt.0.0D0) length=0.0D0 +c$$$ length=sqrt(length) +c$$$ rand_v(1)=rand_v(1)-cost*dc_norm(1,i_sc) +c$$$ rand_v(2)=rand_v(2)-cost*dc_norm(2,i_sc) +c$$$ rand_v(3)=rand_v(3)-cost*dc_norm(3,i_sc) +c$$$ enddo +c$$$ rand_v(1)=rand_v(1)/length +c$$$ rand_v(2)=rand_v(2)/length +c$$$ rand_v(3)=rand_v(3)/length +c$$$ +c$$$ cost=dcos(rad_ang) +c$$$ sint=dsin(rad_ang) +c$$$ dc(1,i_sc)=vbld(i_sc)*(dc_norm(1,i_sc)*cost+rand_v(1)*sint) +c$$$ dc(2,i_sc)=vbld(i_sc)*(dc_norm(2,i_sc)*cost+rand_v(2)*sint) +c$$$ dc(3,i_sc)=vbld(i_sc)*(dc_norm(3,i_sc)*cost+rand_v(3)*sint) +c$$$ dc_norm(1,i_sc)=dc(1,i_sc)*vbld_inv(i_sc) +c$$$ dc_norm(2,i_sc)=dc(2,i_sc)*vbld_inv(i_sc) +c$$$ dc_norm(3,i_sc)=dc(3,i_sc)*vbld_inv(i_sc) +c$$$ c(1,i_sc)=c(1,i)+dc(1,i_sc) +c$$$ c(2,i_sc)=c(2,i)+dc(2,i_sc) +c$$$ c(3,i_sc)=c(3,i)+dc(3,i_sc) +c$$$ +c$$$ call chainbuild_cart +c$$$ +c$$$ return +c$$$ end +c$$$ +c$$$c---------------------------------------------------------------------------- +c$$$ +c$$$ subroutine ss_relax3(i_in,j_in) +c$$$ implicit none +c$$$ +c$$$c Includes +c$$$ include 'DIMENSIONS' +c$$$ include 'COMMON.VAR' +c$$$ include 'COMMON.CHAIN' +c$$$ include 'COMMON.IOUNITS' +c$$$ include 'COMMON.INTERACT' +c$$$ +c$$$c External functions +c$$$ external ran_number +c$$$ double precision ran_number +c$$$ +c$$$c Input arguments +c$$$ integer i_in,j_in +c$$$ +c$$$c Local variables +c$$$ double precision energy_sc(0:n_ene),etot +c$$$ double precision org_dc(3),org_dc_norm(3),org_c(3) +c$$$ double precision ang_pert,rand_fact,exp_fact,beta +c$$$ integer n,i_pert,i +c$$$ logical notdone +c$$$ +c$$$ +c$$$ beta=1.0D0 +c$$$ +c$$$ mask_r=.true. +c$$$ do i=nnt,nct +c$$$ mask_side(i)=0 +c$$$ enddo +c$$$ mask_side(i_in)=1 +c$$$ mask_side(j_in)=1 +c$$$ +c$$$ call etotal_sc(energy_sc) +c$$$ etot=energy_sc(0) +c$$$c write(iout,'(a,3d15.5)')" SS_MC_START ",energy_sc(0), +c$$$c + energy_sc(1),energy_sc(12) +c$$$ +c$$$ notdone=.true. +c$$$ n=0 +c$$$ do while (notdone) +c$$$ if (mod(n,2).eq.0) then +c$$$ i_pert=i_in +c$$$ else +c$$$ i_pert=j_in +c$$$ endif +c$$$ n=n+1 +c$$$ +c$$$ do i=1,3 +c$$$ org_dc(i)=dc(i,i_pert+nres) +c$$$ org_dc_norm(i)=dc_norm(i,i_pert+nres) +c$$$ org_c(i)=c(i,i_pert+nres) +c$$$ enddo +c$$$ ang_pert=ran_number(0.0D0,3.0D0) +c$$$ call perturb_side_chain(i_pert,ang_pert) +c$$$ call etotal_sc(energy_sc) +c$$$ exp_fact=exp(beta*(etot-energy_sc(0))) +c$$$ rand_fact=ran_number(0.0D0,1.0D0) +c$$$ if (rand_fact.lt.exp_fact) then +c$$$c write(iout,'(a,3d15.5)')" SS_MC_ACCEPT ",energy_sc(0), +c$$$c + energy_sc(1),energy_sc(12) +c$$$ etot=energy_sc(0) +c$$$ else +c$$$c write(iout,'(a,3d15.5)')" SS_MC_REJECT ",energy_sc(0), +c$$$c + energy_sc(1),energy_sc(12) +c$$$ do i=1,3 +c$$$ dc(i,i_pert+nres)=org_dc(i) +c$$$ dc_norm(i,i_pert+nres)=org_dc_norm(i) +c$$$ c(i,i_pert+nres)=org_c(i) +c$$$ enddo +c$$$ endif +c$$$ +c$$$ if (n.eq.10000.or.etot.lt.30.0D0) notdone=.false. +c$$$ enddo +c$$$ +c$$$ mask_r=.false. +c$$$ +c$$$ return +c$$$ end +c$$$ +c$$$c---------------------------------------------------------------------------- +c$$$ +c$$$ subroutine ss_relax2(etot,iretcode,nfun,i_in,j_in) +c$$$ implicit none +c$$$ include 'DIMENSIONS' +c$$$ integer liv,lv +c$$$ parameter (liv=60,lv=(77+maxres6*(maxres6+17)/2)) +c$$$********************************************************************* +c$$$* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for * +c$$$* the calling subprogram. * +c$$$* when d(i)=1.0, then v(35) is the length of the initial step, * +c$$$* calculated in the usual pythagorean way. * +c$$$* absolute convergence occurs when the function is within v(31) of * +c$$$* zero. unless you know the minimum value in advance, abs convg * +c$$$* is probably not useful. * +c$$$* relative convergence is when the model predicts that the function * +c$$$* will decrease by less than v(32)*abs(fun). * +c$$$********************************************************************* +c$$$ include 'COMMON.IOUNITS' +c$$$ include 'COMMON.VAR' +c$$$ include 'COMMON.GEO' +c$$$ include 'COMMON.MINIM' +c$$$ include 'COMMON.CHAIN' +c$$$ +c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist +c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar), +c$$$ + orig_ss_dist(maxres2,maxres2) +c$$$ +c$$$ double precision etot +c$$$ integer iretcode,nfun,i_in,j_in +c$$$ +c$$$ external dist +c$$$ double precision dist +c$$$ external ss_func,fdum +c$$$ double precision ss_func,fdum +c$$$ +c$$$ integer iv(liv),uiparm(2) +c$$$ double precision v(lv),x(maxres6),d(maxres6),rdum +c$$$ integer i,j,k +c$$$ +c$$$ +c$$$ call deflt(2,iv,liv,lv,v) +c$$$* 12 means fresh start, dont call deflt +c$$$ iv(1)=12 +c$$$* max num of fun calls +c$$$ if (maxfun.eq.0) maxfun=500 +c$$$ iv(17)=maxfun +c$$$* max num of iterations +c$$$ if (maxmin.eq.0) maxmin=1000 +c$$$ iv(18)=maxmin +c$$$* controls output +c$$$ iv(19)=2 +c$$$* selects output unit +c$$$c iv(21)=iout +c$$$ iv(21)=0 +c$$$* 1 means to print out result +c$$$ iv(22)=0 +c$$$* 1 means to print out summary stats +c$$$ iv(23)=0 +c$$$* 1 means to print initial x and d +c$$$ iv(24)=0 +c$$$* min val for v(radfac) default is 0.1 +c$$$ v(24)=0.1D0 +c$$$* max val for v(radfac) default is 4.0 +c$$$ v(25)=2.0D0 +c$$$c v(25)=4.0D0 +c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease) +c$$$* the sumsl default is 0.1 +c$$$ v(26)=0.1D0 +c$$$* false conv if (act fnctn decrease) .lt. v(34) +c$$$* the sumsl default is 100*machep +c$$$ v(34)=v(34)/100.0D0 +c$$$* absolute convergence +c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4 +c$$$ v(31)=tolf +c$$$ v(31)=1.0D-1 +c$$$* relative convergence +c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-4 +c$$$ v(32)=rtolf +c$$$ v(32)=1.0D-1 +c$$$* controls initial step size +c$$$ v(35)=1.0D-1 +c$$$* large vals of d correspond to small components of step +c$$$ do i=1,6*nres +c$$$ d(i)=1.0D0 +c$$$ enddo +c$$$ +c$$$ do i=0,2*nres +c$$$ do j=1,3 +c$$$ orig_ss_dc(j,i)=dc(j,i) +c$$$ enddo +c$$$ enddo +c$$$ call geom_to_var(nvar,orig_ss_var) +c$$$ +c$$$ do i=1,nres +c$$$ do j=i,nres +c$$$ orig_ss_dist(j,i)=dist(j,i) +c$$$ orig_ss_dist(j+nres,i)=dist(j+nres,i) +c$$$ orig_ss_dist(j,i+nres)=dist(j,i+nres) +c$$$ orig_ss_dist(j+nres,i+nres)=dist(j+nres,i+nres) +c$$$ enddo +c$$$ enddo +c$$$ +c$$$ k=0 +c$$$ do i=1,nres-1 +c$$$ do j=1,3 +c$$$ k=k+1 +c$$$ x(k)=dc(j,i) +c$$$ enddo +c$$$ enddo +c$$$ do i=2,nres-1 +c$$$ if (ialph(i,1).gt.0) then +c$$$ do j=1,3 +c$$$ k=k+1 +c$$$ x(k)=dc(j,i+nres) +c$$$ enddo +c$$$ endif +c$$$ enddo +c$$$ +c$$$ uiparm(1)=i_in +c$$$ uiparm(2)=j_in +c$$$ call smsno(k,d,x,ss_func,iv,liv,lv,v,uiparm,rdum,fdum) +c$$$ etot=v(10) +c$$$ iretcode=iv(1) +c$$$ nfun=iv(6)+iv(30) +c$$$ +c$$$ k=0 +c$$$ do i=1,nres-1 +c$$$ do j=1,3 +c$$$ k=k+1 +c$$$ dc(j,i)=x(k) +c$$$ enddo +c$$$ enddo +c$$$ do i=2,nres-1 +c$$$ if (ialph(i,1).gt.0) then +c$$$ do j=1,3 +c$$$ k=k+1 +c$$$ dc(j,i+nres)=x(k) +c$$$ enddo +c$$$ endif +c$$$ enddo +c$$$ call chainbuild_cart +c$$$ +c$$$ return +c$$$ end +c$$$ +c$$$C----------------------------------------------------------------------------- +c$$$ +c$$$ subroutine ss_func(n,x,nf,f,uiparm,urparm,ufparm) +c$$$ implicit none +c$$$ include 'DIMENSIONS' +c$$$ include 'COMMON.DERIV' +c$$$ include 'COMMON.IOUNITS' +c$$$ include 'COMMON.VAR' +c$$$ include 'COMMON.CHAIN' +c$$$ include 'COMMON.INTERACT' +c$$$ include 'COMMON.SBRIDGE' +c$$$ +c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist +c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar), +c$$$ + orig_ss_dist(maxres2,maxres2) +c$$$ +c$$$ integer n +c$$$ double precision x(maxres6) +c$$$ integer nf +c$$$ double precision f +c$$$ integer uiparm(2) +c$$$ real*8 urparm(1) +c$$$ external ufparm +c$$$ double precision ufparm +c$$$ +c$$$ external dist +c$$$ double precision dist +c$$$ +c$$$ integer i,j,k,ss_i,ss_j +c$$$ double precision tempf,var(maxvar) +c$$$ +c$$$ +c$$$ ss_i=uiparm(1) +c$$$ ss_j=uiparm(2) +c$$$ f=0.0D0 +c$$$ +c$$$ k=0 +c$$$ do i=1,nres-1 +c$$$ do j=1,3 +c$$$ k=k+1 +c$$$ dc(j,i)=x(k) +c$$$ enddo +c$$$ enddo +c$$$ do i=2,nres-1 +c$$$ if (ialph(i,1).gt.0) then +c$$$ do j=1,3 +c$$$ k=k+1 +c$$$ dc(j,i+nres)=x(k) +c$$$ enddo +c$$$ endif +c$$$ enddo +c$$$ call chainbuild_cart +c$$$ +c$$$ call geom_to_var(nvar,var) +c$$$ +c$$$c Constraints on all angles +c$$$ do i=1,nvar +c$$$ tempf=var(i)-orig_ss_var(i) +c$$$ f=f+tempf*tempf +c$$$ enddo +c$$$ +c$$$c Constraints on all distances +c$$$ do i=1,nres-1 +c$$$ if (i.gt.1) then +c$$$ tempf=dist(i+nres,i)-orig_ss_dist(i+nres,i) +c$$$ f=f+tempf*tempf +c$$$ endif +c$$$ do j=i+1,nres +c$$$ tempf=dist(j,i)-orig_ss_dist(j,i) +c$$$ if (tempf.lt.0.0D0 .or. j.eq.i+1) f=f+tempf*tempf +c$$$ tempf=dist(j+nres,i)-orig_ss_dist(j+nres,i) +c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf +c$$$ tempf=dist(j,i+nres)-orig_ss_dist(j,i+nres) +c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf +c$$$ tempf=dist(j+nres,i+nres)-orig_ss_dist(j+nres,i+nres) +c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf +c$$$ enddo +c$$$ enddo +c$$$ +c$$$c Constraints for the relevant CYS-CYS +c$$$ tempf=dist(nres+ss_i,nres+ss_j)-8.0D0 +c$$$ f=f+tempf*tempf +c$$$CCCCCCCCCCCCCCCCC ADD SOME ANGULAR STUFF +c$$$ +c$$$c$$$ if (nf.ne.nfl) then +c$$$c$$$ write(iout,'(a,i10,2d15.5)')"IN DIST_FUNC (NF,F,DIST)",nf, +c$$$c$$$ + f,dist(5+nres,14+nres) +c$$$c$$$ endif +c$$$ +c$$$ nfl=nf +c$$$ +c$$$ return +c$$$ end +c$$$ +c$$$C----------------------------------------------------------------------------- diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index e7199c4..b5bfc9e 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -1451,6 +1451,20 @@ C evdw=evdw+evdwij if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') & 'evdw',i,j,evdwij,' ss' +C triple bond artifac removal + do k=j+1,iend(i,iint) +C search over all next residues + if (dyn_ss_mask(k)) then +C check if they are cysteins +C write(iout,*) 'k=',k + call triple_ssbond_ene(i,j,k,evdwij) +C call the energy function that removes the artifical triple disulfide +C bond the soubroutine is located in ssMD.F + evdw=evdw+evdwij + if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') + & 'evdw',i,j,evdwij,'tss' + endif!dyn_ss_mask(k) + enddo! k ELSE ind=ind+1 itypj=iabs(itype(j)) diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index 7853e58..3296e45 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -1256,7 +1256,7 @@ c lprint=.false. C C Define the constants of the disulfide bridge C - ebr=-5.50D0 + ebr=-12.00D0 c c Old arbitrary potential - commented out. c diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index 8dbd3fa..6984aeb 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -680,6 +680,10 @@ C 12/1/95 Added weight for the multi-body term WCORR call reada(weightcard,"V2SS",v2ss,7.61d0) call reada(weightcard,"V3SS",v3ss,13.7d0) call reada(weightcard,"EBR",ebr,-5.50D0) + call reada(weightcard,"DTRISS",dtriss,1.0D0) + call reada(weightcard,"ATRISS",atriss,0.3D0) + call reada(weightcard,"BTRISS",btriss,0.02D0) + call reada(weightcard,"CTRISS",ctriss,1.0D0) dyn_ss=(index(weightcard,'DYN_SS').gt.0) do i=1,maxres dyn_ss_mask(i)=.false. diff --git a/source/unres/src_MD-M/ssMD.F b/source/unres/src_MD-M/ssMD.F index 15800ae..50e1040 100644 --- a/source/unres/src_MD-M/ssMD.F +++ b/source/unres/src_MD-M/ssMD.F @@ -1949,3 +1949,151 @@ c$$$ return c$$$ end c$$$ c$$$C----------------------------------------------------------------------------- +c$$$C----------------------------------------------------------------------------- + subroutine triple_ssbond_ene(resi,resj,resk,eij) + include 'DIMENSIONS' + include 'COMMON.SBRIDGE' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.VAR' + include 'COMMON.IOUNITS' + include 'COMMON.CALC' +#ifndef CLUST +#ifndef WHAM + include 'COMMON.MD' +#endif +#endif + +c External functions + double precision h_base + external h_base + +c Input arguments + integer resi,resj,resk + +c Output arguments + double precision eij,eij1,eij2,eij3 + +c Local variables + logical havebond +c integer itypi,itypj,k,l + double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi + double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij + double precision xik,yik,zik,xjk,yjk,zjk + double precision sig0ij,ljd,sig,fac,e1,e2 + double precision dcosom1(3),dcosom2(3),ed + double precision pom1,pom2 + double precision ljA,ljB,ljXs + double precision d_ljB(1:3) + double precision ssA,ssB,ssC,ssXs + double precision ssxm,ljxm,ssm,ljm + double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3) + + i=resi + j=resj + k=resk +C write(iout,*) resi,resj,resk + itypi=itype(i) + dxi=dc_norm(1,nres+i) + dyi=dc_norm(2,nres+i) + dzi=dc_norm(3,nres+i) + dsci_inv=vbld_inv(i+nres) + xi=c(1,nres+i) + yi=c(2,nres+i) + zi=c(3,nres+i) + + itypj=itype(j) + xj=c(1,nres+j) + yj=c(2,nres+j) + zj=c(3,nres+j) + + dxj=dc_norm(1,nres+j) + dyj=dc_norm(2,nres+j) + dzj=dc_norm(3,nres+j) + dscj_inv=vbld_inv(j+nres) + itypk=itype(k) + xk=c(1,nres+k) + yk=c(2,nres+k) + zk=c(3,nres+k) + + dxk=dc_norm(1,nres+k) + dyk=dc_norm(2,nres+k) + dzk=dc_norm(3,nres+k) + dscj_inv=vbld_inv(k+nres) + xij=xj-xi + xik=xk-xi + xjk=xk-xj + yij=yj-yi + yik=yk-yi + yjk=yk-yj + zij=zj-zi + zik=zk-zi + zjk=zk-zj + rrij=(xij*xij+yij*yij+zij*zij) + rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse + rrik=(xik*xik+yik*yik+zik*zik) + rik=dsqrt(rrik) + rrjk=(xjk*xjk+yjk*yjk+zjk*zjk) + rjk=dsqrt(rrjk) +C there are three combination of distances for each trisulfide bonds +C The first case the ith atom is the center +C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first +C distance y is second distance the a,b,c,d are parameters derived for +C this problem d parameter was set as a penalty currenlty set to 1. + eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**2+ctriss) +C second case jth atom is center + eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**2+ctriss) +C the third case kth atom is the center + eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**2+ctriss) +C eij2=0.0 +C eij3=0.0 +C eij1=0.0 + eij=eij1+eij2+eij3 +C write(iout,*)i,j,k,eij +C The energy penalty calculated now time for the gradient part +C derivative over rij + fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik)) + &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk)) + gg(1)=xij*fac/rij + gg(2)=yij*fac/rij + gg(3)=zij*fac/rij + do m=1,3 + gvdwx(m,i)=gvdwx(m,i)-gg(m) + gvdwx(m,j)=gvdwx(m,j)+gg(m) + enddo + do l=1,3 + gvdwc(l,i)=gvdwc(l,i)-gg(l) + gvdwc(l,j)=gvdwc(l,j)+gg(l) + enddo +C now derivative over rik + fac=-eij1**2/dtriss*(-2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik)) + &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk)) + gg(1)=xik*fac/rik + gg(2)=yik*fac/rik + gg(3)=zik*fac/rik + do m=1,3 + gvdwx(m,i)=gvdwx(m,i)-gg(m) + gvdwx(m,k)=gvdwx(m,k)+gg(m) + enddo + do l=1,3 + gvdwc(l,i)=gvdwc(l,i)-gg(l) + gvdwc(l,k)=gvdwc(l,k)+gg(l) + enddo +C now derivative over rjk + fac=-eij2**2/dtriss*(-2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))- + &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk)) + gg(1)=xjk*fac/rjk + gg(2)=yjk*fac/rjk + gg(3)=zjk*fac/rjk + do m=1,3 + gvdwx(m,j)=gvdwx(m,j)-gg(m) + gvdwx(m,k)=gvdwx(m,k)+gg(m) + enddo + do l=1,3 + gvdwc(l,j)=gvdwc(l,j)-gg(l) + gvdwc(l,k)=gvdwc(l,k)+gg(l) + enddo + return + end -- 1.7.9.5