introduction of different ftors for each site
[unres.git] / source / unres / src_MD-M / readrtns_CSA.F
index 265b705..1cda5ea 100644 (file)
@@ -295,11 +295,11 @@ cd       endif
           do i=1,nrep
            iremd_m_total=iremd_m_total+remd_m(i)
           enddo
-          write (iout,*) 'Total number of replicas ',iremd_m_total
+           write (iout,*) 'Total number of replicas ',iremd_m_total
+          endif
          endif
-      endif
       if(me.eq.king.or..not.out1file) 
-     & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
+     &   write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
       return
       end
 c--------------------------------------------------------------------------
@@ -345,6 +345,7 @@ C
       rest = index(controlcard,"REST").gt.0
       tbf = index(controlcard,"TBF").gt.0
       usampl = index(controlcard,"USAMPL").gt.0
+
       mdpdb = index(controlcard,"MDPDB").gt.0
       call reada(controlcard,"T_BATH",t_bath,300.0d0)
       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1) 
@@ -538,7 +539,7 @@ C
       integer rescode
       double precision x(maxvar)
       character*256 pdbfile
-      character*320 weightcard
+      character*400 weightcard
       character*80 weightcard_t,ucase
       dimension itype_pdb(maxres)
       common /pizda/ itype_pdb
@@ -550,54 +551,54 @@ C
 C Body
 C
 C Read weights of the subsequent energy terms.
-      call card_concat(weightcard)
-      call reada(weightcard,'WLONG',wlong,1.0D0)
-      call reada(weightcard,'WSC',wsc,wlong)
-      call reada(weightcard,'WSCP',wscp,wlong)
-      call reada(weightcard,'WELEC',welec,1.0D0)
-      call reada(weightcard,'WVDWPP',wvdwpp,welec)
-      call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
-      call reada(weightcard,'WCORR4',wcorr4,0.0D0)
-      call reada(weightcard,'WCORR5',wcorr5,0.0D0)
-      call reada(weightcard,'WCORR6',wcorr6,0.0D0)
-      call reada(weightcard,'WTURN3',wturn3,1.0D0)
-      call reada(weightcard,'WTURN4',wturn4,1.0D0)
-      call reada(weightcard,'WTURN6',wturn6,1.0D0)
-      call reada(weightcard,'WSCCOR',wsccor,1.0D0)
-      call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
-      call reada(weightcard,'WBOND',wbond,1.0D0)
-      call reada(weightcard,'WTOR',wtor,1.0D0)
-      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,'SCAL14',scal14,0.4D0)
-      call reada(weightcard,'SCALSCP',scalscp,1.0d0)
-      call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
-      call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
-      call reada(weightcard,'TEMP0',temp0,300.0d0)
-      if (index(weightcard,'SOFT').gt.0) ipot=6
+       call card_concat(weightcard)
+       call reada(weightcard,'WLONG',wlong,1.0D0)
+       call reada(weightcard,'WSC',wsc,wlong)
+       call reada(weightcard,'WSCP',wscp,wlong)
+       call reada(weightcard,'WELEC',welec,1.0D0)
+       call reada(weightcard,'WVDWPP',wvdwpp,welec)
+       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
+       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
+       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
+       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
+       call reada(weightcard,'WTURN3',wturn3,1.0D0)
+       call reada(weightcard,'WTURN4',wturn4,1.0D0)
+       call reada(weightcard,'WTURN6',wturn6,1.0D0)
+       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
+       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
+       call reada(weightcard,'WBOND',wbond,1.0D0)
+       call reada(weightcard,'WTOR',wtor,1.0D0)
+       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,'SCAL14',scal14,0.4D0)
+       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
+       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
+       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
+       call reada(weightcard,'TEMP0',temp0,300.0d0)
+       if (index(weightcard,'SOFT').gt.0) ipot=6
 C 12/1/95 Added weight for the multi-body term WCORR
-      call reada(weightcard,'WCORRH',wcorr,1.0D0)
-      if (wcorr4.gt.0.0d0) wcorr=wcorr4
-      weights(1)=wsc
-      weights(2)=wscp
-      weights(3)=welec
-      weights(4)=wcorr
-      weights(5)=wcorr5
-      weights(6)=wcorr6
-      weights(7)=wel_loc
-      weights(8)=wturn3
-      weights(9)=wturn4
-      weights(10)=wturn6
-      weights(11)=wang
-      weights(12)=wscloc
-      weights(13)=wtor
-      weights(14)=wtor_d
-      weights(15)=wstrain
-      weights(16)=wvdwpp
-      weights(17)=wbond
-      weights(18)=scal14
-      weights(21)=wsccor
+       call reada(weightcard,'WCORRH',wcorr,1.0D0)
+       if (wcorr4.gt.0.0d0) wcorr=wcorr4
+       weights(1)=wsc
+       weights(2)=wscp
+       weights(3)=welec
+       weights(4)=wcorr
+       weights(5)=wcorr5
+       weights(6)=wcorr6
+       weights(7)=wel_loc
+       weights(8)=wturn3
+       weights(9)=wturn4
+       weights(10)=wturn6
+       weights(11)=wang
+       weights(12)=wscloc
+       weights(13)=wtor
+       weights(14)=wtor_d
+       weights(15)=wstrain
+       weights(16)=wvdwpp
+       weights(17)=wbond
+       weights(18)=scal14
+       weights(21)=wsccor
       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,
@@ -679,13 +680,49 @@ 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,"ATRISS",atriss,0.301D0)
+      call reada(weightcard,"BTRISS",btriss,0.021D0)
+      call reada(weightcard,"CTRISS",ctriss,1.001D0)
+      call reada(weightcard,"DTRISS",dtriss,1.001D0)
+      write (iout,*) "ATRISS=", atriss
+      write (iout,*) "BTRISS=", btriss
+      write (iout,*) "CTRISS=", ctriss
+      write (iout,*) "DTRISS=", dtriss
+      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
+        if (wstrain.ne.0.0) then
+         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
+        else
+          ss_depth=0.0
+        endif
+      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
-c       print *,'indpdb=',indpdb,' pdbref=',pdbref
+       write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
+       write (iout,*)" HT",Ht
+       print *,'indpdb=',indpdb,' pdbref=',pdbref
       endif
       if (indpdb.gt.0 .or. pdbref) then
         read(inp,'(a)') pdbfile
@@ -697,9 +734,11 @@ c       print *,'indpdb=',indpdb,' pdbref=',pdbref
   33    write (iout,'(a)') 'Error opening PDB file.'
         stop
   34    continue
-c        print *,'Begin reading pdb data'
+c        write (iout,*) 'Begin reading pdb data'
+c        call flush(iout)
         call readpdb
-c        print *,'Finished reading pdb data'
+c        write (iout,*) 'Finished reading pdb data'
+c        call flush(iout)
         if(me.eq.king.or..not.out1file)
      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
      &   ' nstart_sup=',nstart_sup
@@ -789,21 +828,23 @@ C 8/13/98 Set limits to generating the dihedral angles
       enddo
       read (inp,*) ndih_constr
       if (ndih_constr.gt.0) then
-        read (inp,*) ftors
-        read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
+C        read (inp,*) ftors
+        read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
+     &  i=1,ndih_constr)
         if(me.eq.king.or..not.out1file)then
          write (iout,*) 
      &   'There are',ndih_constr,' constraints on phi angles.'
          do i=1,ndih_constr
-          write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
+          write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
+     &    ftors(i)
          enddo
         endif
         do i=1,ndih_constr
           phi0(i)=deg2rad*phi0(i)
           drange(i)=deg2rad*drange(i)
         enddo
-        if(me.eq.king.or..not.out1file)
-     &   write (iout,*) 'FTORS',ftors
+C        if(me.eq.king.or..not.out1file)
+C     &   write (iout,*) 'FTORS',ftors
         do i=1,ndih_constr
           ii = idih_constr(i)
           phibound(1,ii) = phi0(i)-drange(i)
@@ -896,7 +937,9 @@ czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
           enddo
           call contact(.true.,ncont_ref,icont_ref,co)
         endif
-c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
+        endif
+        print *, "A TU"
+        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
         call flush(iout)
         if (constr_dist.gt.0) call read_dist_constr
         write (iout,*) "After read_dist_constr nhpb",nhpb
@@ -916,7 +959,7 @@ c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
         enddo
         endif
-      endif
+C      endif
       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and. 
      &    modecalc.ne.10) then
@@ -1055,18 +1098,35 @@ 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)
@@ -1082,6 +1142,7 @@ cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)') 
      &  'Processor',myrank,': end reading molecular data.'
 #endif
+      print *,"A TU?"
       return
       end
 c--------------------------------------------------------------------------
@@ -1123,17 +1184,19 @@ C Read information about disulfide bridges.
       include 'COMMON.SETUP'
 C Read bridging residues.
       read (inp,*) ns,(iss(i),i=1,ns)
-c      print *,'ns=',ns
+      print *,'ns=',ns
       if(me.eq.king.or..not.out1file)
      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
 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 +1207,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
@@ -2202,7 +2266,8 @@ c-------------------------------------------------------------------------------
       integer ifrag_(2,100),ipair_(2,100)
       double precision wfrag_(100),wpair_(100)
       character*500 controlcard
-c      write (iout,*) "Calling read_dist_constr"
+      print *, "WCHODZE"
+      write (iout,*) "Calling read_dist_constr"
 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 c      call flush(iout)
       call card_concat(controlcard)
@@ -2236,11 +2301,11 @@ c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 c            write (iout,*) "j",j," k",k
             ddjk=dist(j,k)
             if (constr_dist.eq.1) then
-              nhpb=nhpb+1
-              ihpb(nhpb)=j
-              jhpb(nhpb)=k
+            nhpb=nhpb+1
+            ihpb(nhpb)=j
+            jhpb(nhpb)=k
               dhpb(nhpb)=ddjk
-              forcon(nhpb)=wfrag_(i) 
+            forcon(nhpb)=wfrag_(i) 
             else if (constr_dist.eq.2) then
               if (ddjk.le.dist_cut) then
                 nhpb=nhpb+1
@@ -2296,11 +2361,30 @@ c            write (iout,*) "j",j," k",k
         enddo
         endif
       enddo 
+      print *,ndist_
       do i=1,ndist_
-        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
+        if (constr_dist.eq.11) then
+        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
+     &     ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
+        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
+        else
+C        print *,"in else"
+        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
+     &     ibecarb(i),forcon(nhpb+1)
+        endif
         if (forcon(nhpb+1).gt.0.0d0) then
           nhpb=nhpb+1
-          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
+          if (ibecarb(i).gt.0) then
+            ihpb(i)=ihpb(i)+nres
+            jhpb(i)=jhpb(i)+nres
+          endif
+          if (dhpb(nhpb).eq.0.0d0)
+     &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
+        endif
+C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
+C        if (forcon(nhpb+1).gt.0.0d0) then
+C          nhpb=nhpb+1
+C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 #ifdef MPI
           if (.not.out1file .or. me.eq.king)
      &    write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
@@ -2309,7 +2393,7 @@ c            write (iout,*) "j",j," k",k
           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 #endif
-        endif
+
       enddo
       call flush(iout)
       return
@@ -2418,7 +2502,7 @@ C
 #endif
       if (OKRandom) then
 c        r1 = prng_next(me)
-         r1=ran_number(0.0D0,1.0D0)
+        r1=ran_number(0.0D0,1.0D0)
         if(me.eq.king)
      &   write (iout,*) 'ran_num',r1
         if (r1.lt.0.0d0) OKRandom=.false.