2D replica exchange with homology constraints
[unres.git] / source / unres / src_MD / readrtns.F
index 88b7242..18e2c2e 100644 (file)
@@ -96,6 +96,7 @@ c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
 C Set up the time limit (caution! The time must be input in minutes!)
       read_cart=index(controlcard,'READ_CART').gt.0
       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
+      call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
@@ -749,6 +750,10 @@ C 12/1/95 Added weight for the multi-body term WCORR
        call reada(weightcard,'WTORD',wtor_d,1.0D0)
        call reada(weightcard,'WANG',wang,1.0D0)
        call reada(weightcard,'WSCLOC',wscloc,1.0D0)
+       call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
+       call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
+       call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
+       call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
        call reada(weightcard,'SCAL14',scal14,0.4D0)
        call reada(weightcard,'SCALSCP',scalscp,1.0d0)
        call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
@@ -778,11 +783,16 @@ C 12/1/95 Added weight for the multi-body term WCORR
        weights(18)=scal14
        weights(21)=wsccor
       endif
+       weights(25)=wdfa_dist
+       weights(26)=wdfa_tor
+       weights(27)=wdfa_nei
+       weights(28)=wdfa_beta
 
       if(me.eq.king.or..not.out1file)
      & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
-     &  wturn4,wturn6
+     &  wturn4,wturn6,
+     &  wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
    10 format (/'Energy-term weights (unscaled):'//
      & 'WSCC=   ',f10.6,' (SC-SC)'/
      & 'WSCP=   ',f10.6,' (SC-p)'/
@@ -801,7 +811,11 @@ C 12/1/95 Added weight for the multi-body term WCORR
      & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
-     & 'WTURN6= ',f10.6,' (turns, 6th order)')
+     & 'WTURN6= ',f10.6,' (turns, 6th order)'/
+     & 'WDFA_D= ',f10.6,' (DFA, distance)'   /
+     & 'WDFA_T= ',f10.6,' (DFA, torsional)'   /
+     & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)'   /
+     & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
       if(me.eq.king.or..not.out1file)then
        if (wcorr4.gt.0.0d0) then
         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
@@ -829,7 +843,8 @@ C 12/1/95 Added weight for the multi-body term WCORR
       if(me.eq.king.or..not.out1file)
      & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
-     &  wturn4,wturn6
+     &  wturn4,wturn6,
+     &  wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
    22 format (/'Energy-term weights (scaled):'//
      & 'WSCC=   ',f10.6,' (SC-SC)'/
      & 'WSCP=   ',f10.6,' (SC-p)'/
@@ -848,7 +863,11 @@ C 12/1/95 Added weight for the multi-body term WCORR
      & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
-     & 'WTURN6= ',f10.6,' (turns, 6th order)')
+     & 'WTURN6= ',f10.6,' (turns, 6th order)'/
+     & 'WDFA_D= ',f10.6,' (DFA, distance)'   /
+     & 'WDFA_T= ',f10.6,' (DFA, torsional)'   /
+     & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)'   /
+     & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
       if(me.eq.king.or..not.out1file)
      & write (iout,*) "Reference temperature for weights calculation:",
      &  temp0
@@ -940,8 +959,8 @@ C 10/03/12 Adam: Recalculate coordinates with new side chain positions
          call chainbuild
         endif  
 C Following 2 lines for diagnostics; comment out if not needed
-        write (iout,*) "After sideadd"
-        call intout
+c        write (iout,*) "After sideadd"
+c        call intout
       endif
       if (indpdb.eq.0) then
 C Read sequence if not taken from the pdb file.
@@ -1039,6 +1058,21 @@ C 8/13/98 Set limits to generating the dihedral angles
 cd      print *,'NNT=',NNT,' NCT=',NCT
       if (itype(1).eq.21) nnt=2
       if (itype(nres).eq.21) nct=nct-1
+
+C     Bartek:READ init_vars
+C     Initialize variables!
+C     Juyong:READ read_info
+C     READ fragment information!!
+C     both routines should be in dfa.F file!!
+
+      if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
+     &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
+       call init_dfa_vars
+       print*, 'init_dfa_vars finished!'
+       call read_dfa_info
+       print*, 'read_dfa_info finished!'
+      endif
+C
       if (pdbref) then
         if(me.eq.king.or..not.out1file)
      &   write (iout,'(a,i3)') 'nsup=',nsup
@@ -1127,8 +1161,17 @@ czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
 c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
       if (constr_dist.gt.0) then
         call read_dist_constr
-        call hpb_partition
       endif
+
+
+      if (constr_homology.gt.0) then
+        call read_constr_homology
+      else
+        homol_nset=0
+      endif
+
+
+      if (nhpb.gt.0) call hpb_partition
 c      write (iout,*) "After read_dist_constr nhpb",nhpb
 c      call flush(iout)
       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
@@ -1251,18 +1294,6 @@ C Generate distance constraints, if the PDB structure is to be regularized.
         write (iout,'(20i4)') (iss(i),i=1,ns)
        if (dyn_ss) then
           write(iout,*)"Running with dynamic disulfide-bond formation"
-          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
        else
         write (iout,'(/a/)') 'Pre-formed links are:' 
        do i=1,nss
@@ -1270,14 +1301,27 @@ C Generate distance constraints, if the PDB structure is to be regularized.
          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))
@@ -1340,10 +1384,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)
@@ -1354,7 +1400,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
@@ -2367,6 +2414,7 @@ c-------------------------------------------------------------------------------
       include 'COMMON.CHAIN'
       include 'COMMON.IOUNITS'
       include 'COMMON.MD'
+      include 'COMMON.CONTROL'
       open(irest2,file=rest2name,status='unknown')
       read(irest2,*) totT,EK,potE,totE,t_bath
       do i=1,2*nres
@@ -2375,7 +2423,7 @@ c-------------------------------------------------------------------------------
       do i=1,2*nres
          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
       enddo
-      if(usampl) then
+      if(usampl.or.homol_nset.gt.1) then
              read (irest2,*) iset
       endif
       close(irest2)
@@ -2482,7 +2530,7 @@ c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
         if (wfrag_(i).gt.0.0d0) then
         do j=ifrag_(1,i),ifrag_(2,i)-1
           do k=j+1,ifrag_(2,i)
-            write (iout,*) "j",j," k",k
+c            write (iout,*) "j",j," k",k
             ddjk=dist(j,k)
             if (constr_dist.eq.1) then
             nhpb=nhpb+1
@@ -2572,6 +2620,135 @@ c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
       return
       end
 c-------------------------------------------------------------------------------
+
+      subroutine read_constr_homology
+
+      include 'DIMENSIONS'
+#ifdef MPI
+      include 'mpif.h'
+#endif
+      include 'COMMON.SETUP'
+      include 'COMMON.CONTROL'
+      include 'COMMON.CHAIN'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.MD'
+      include 'COMMON.GEO'
+      include 'COMMON.INTERACT'
+      double precision odl_temp,sigma_odl_temp
+      common /przechowalnia/ odl_temp(maxres,maxres,max_template),
+     &    sigma_odl_temp(maxres,maxres,max_template)
+      character*2 kic2
+      character*24 model_ki_dist, model_ki_angle
+      character*500 controlcard
+      character*3200 controlcard1
+      integer ki, i, j, k, l
+      logical lprn /.true./
+
+      call card_concat(controlcard)
+      call readi(controlcard,"HOMOL_NSET",homol_nset,1)       
+      if (homol_nset.gt.1)then
+         call card_concat(controlcard)
+         read(controlcard,*) (waga_dist(i),i=1,homol_nset) 
+         call card_concat(controlcard)
+         read(controlcard,*) (waga_angle(i),i=1,homol_nset) 
+         write(iout,*) "iset distance_weight angle_weight"
+         do i=1,homol_nset
+           write(iout,*) i,waga_dist(i),waga_angle(i)
+         enddo
+      else
+       iset=1
+       call reada(controlcard,"HOMOL_DIST",waga_dist(1),1.0d0)
+       call reada(controlcard,"HOMOL_ANGLE",waga_angle(1),1.0d0)
+      endif
+
+      write (iout,*) "nnt",nnt," nct",nct
+      call flush(iout)
+      lim_odl=0
+      lim_dih=0
+      do i=1,nres
+        do j=i+2,nres
+          do ki=1,constr_homology
+            sigma_odl_temp(i,j,ki)=0.0d0
+            odl_temp(i,j,ki)=0.0d0
+          enddo
+        enddo
+      enddo
+      do i=1,nres-3
+        do ki=1,constr_homology
+          dih(ki,i)=0.0d0
+          sigma_dih(ki,i)=0.0d0
+        enddo
+      enddo
+      do ki=1,constr_homology
+          write(kic2,'(i2)') ki
+          if (ki.le.9) kic2="0"//kic2(2:2)
+
+          model_ki_dist="model"//kic2//".dist"
+          model_ki_angle="model"//kic2//".angle"
+        open (ientin,file=model_ki_dist,status='old')
+        do irec=1,maxdim !petla do czytania wiezow na odleglosc
+          read (ientin,*,end=1401) i, j, odl_temp(i+nnt-1,j+nnt-1,ki),
+     &       sigma_odl_temp(i+nnt-1,j+nnt-1,ki)
+          odl_temp(j+nnt-1,i+nnt-1,ki)=odl_temp(i+nnt-1,j+nnt-1,ki)
+          sigma_odl_temp(j+nnt-1,i+nnt-1,ki)=
+     &     sigma_odl_temp(i+nnt-1,j+nnt-1,ki)
+        enddo
+ 1401 continue
+        close (ientin)
+        open (ientin,file=model_ki_angle,status='old')
+        do irec=1,maxres-3 !petla do czytania wiezow na katy torsyjne
+          read (ientin,*,end=1402) i, j, k,l,dih(ki,i+nnt-1),
+     &      sigma_dih(ki,i+nnt-1)
+          if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1
+          sigma_dih(ki,i+nnt-1)=1.0d0/sigma_dih(ki,i+nnt-1)**2
+        enddo
+ 1402 continue
+        close (ientin)
+      enddo
+      ii=0
+      write (iout,*) "nnt",nnt," nct",nct
+      do i=nnt,nct-2
+        do j=i+2,nct
+          ki=1
+c          write (iout,*) "i",i," j",j," constr_homology",constr_homology
+          do while (ki.le.constr_homology .and.
+     &        sigma_odl_temp(i,j,ki).le.0.0d0)
+c            write (iout,*) "ki",ki," sigma_odl",sigma_odl_temp(i,j,ki)
+            ki=ki+1
+          enddo
+c          write (iout,*) "ki",ki
+          if (ki.gt.constr_homology) cycle
+          ii=ii+1
+          ires_homo(ii)=i
+          jres_homo(ii)=j
+          do ki=1,constr_homology
+            odl(ki,ii)=odl_temp(i,j,ki)
+            sigma_odl(ki,ii)=1.0d0/sigma_odl_temp(i,j,ki)**2
+          enddo      
+        enddo
+      enddo
+      lim_odl=ii
+      if (constr_homology.gt.0) call homology_partition
+c Print restraints
+      if (.not.lprn) return
+      write (iout,*) "Distance restraints from templates"
+      do ii=1,lim_odl
+        write(iout,'(3i5,10(2f8.2,4x))') ii,ires_homo(ii),jres_homo(ii),
+     &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
+      enddo
+      write (iout,*) "Dihedral angle restraints from templates"
+      do i=nnt,lim_dih
+        write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
+     &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
+      enddo
+c      write(iout,*) "TEST CZYTANIA1",odl(1,2,1),odl(1,3,1),odl(1,4,1)
+c      write(iout,*) "TEST CZYTANIA2",dih(1,1),dih(2,1),dih(3,1)
+
+
+      return
+      end
+c----------------------------------------------------------------------
+
 #ifdef WINIFL
       subroutine flush(iu)
       return
@@ -2583,6 +2760,7 @@ c-------------------------------------------------------------------------------
       return
       end
 #endif
+
 c------------------------------------------------------------------------------
       subroutine copy_to_tmp(source)
       include "DIMENSIONS"