X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Freadrtns.F;h=af416e8bf70e6e45c788e47968138a32b4eccec7;hb=c1b8e24ad8ccf66d0617b0055089ea0c1193840d;hp=5b74f874d9c470a0f146b83d5f54bff5ea61f5ef;hpb=4bb2cdc2f60d154b9170e189d6a36992128e81f9;p=unres.git

diff --git a/source/unres/src_MD/readrtns.F b/source/unres/src_MD/readrtns.F
index 5b74f87..af416e8 100644
--- a/source/unres/src_MD/readrtns.F
+++ b/source/unres/src_MD/readrtns.F
@@ -860,12 +860,40 @@ 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.
+      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
        print *,'indpdb=',indpdb,' pdbref=',pdbref
       endif
       if (indpdb.gt.0 .or. pdbref) then
@@ -916,8 +944,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.
@@ -948,9 +976,9 @@ c      print *,nres
 c      print '(20i4)',(itype(i),i=1,nres)
       do i=1,nres
 #ifdef PROCOR
-        if (itype(i).eq.21 .or. itype(i+1).eq.21) then
+        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
 #else
-        if (itype(i).eq.21) then
+        if (itype(i).eq.ntyp1) then
 #endif
           itel(i)=0
 #ifdef PROCOR
@@ -1013,8 +1041,8 @@ C 8/13/98 Set limits to generating the dihedral angles
 #endif
       nct=nres
 cd      print *,'NNT=',NNT,' NCT=',NCT
-      if (itype(1).eq.21) nnt=2
-      if (itype(nres).eq.21) nct=nct-1
+      if (itype(1).eq.ntyp1) nnt=2
+      if (itype(nres).eq.ntyp1) nct=nct-1
       if (pdbref) then
         if(me.eq.king.or..not.out1file)
      &   write (iout,'(a,i3)') 'nsup=',nsup
@@ -1225,18 +1253,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)
@@ -1300,10 +1345,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)
@@ -1314,7 +1361,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
@@ -1908,7 +1956,12 @@ c----------------------------------------------------------------------------
       do i=1,dim
         tablica(i)=default
       enddo
+#ifdef G77
+      aux=lancuch(:ilen(lancuch))//"="	
+      iread=index(rekord,aux)
+#else
       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
+#endif
       if (iread.eq.0) return
       iread=iread+ilen(lancuch)+1
       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
@@ -1926,7 +1979,12 @@ c----------------------------------------------------------------------------
       do i=1,dim
         tablica(i)=default
       enddo
+#ifdef G77
+      aux=lancuch(:ilen(lancuch))//"="
+      iread=index(rekord,aux)
+#else
       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
+#endif
       if (iread.eq.0) return
       iread=iread+ilen(lancuch)+1
       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
@@ -2442,7 +2500,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
@@ -2567,12 +2625,14 @@ c------------------------------------------------------------------------------
       include "DIMENSIONS"
       include "COMMON.IOUNITS"
       character*(*) source
+      character* 256 tmpfile
       integer ilen
       external ilen
+      tmpfile=source(:ilen(source))
       write (*,*) "Moving ",source(:ilen(source)),
      & " from temporary directory to working directory"
-      write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
-      call system("/bin/mv "//source(:ilen(source))//" "//curdir)
+      write (*,*) "/bin/mv "//tmpfile//" "//curdir
+      call system("/bin/mv "//tmpfile//" "//curdir)
       return
       end
 c------------------------------------------------------------------------------