update

[unres.git] / source / cluster / wham / src / readrtns.F

diff --git a/source/cluster/wham/src/readrtns.F b/source/cluster/wham/src/readrtns.F
index 8e63ff8..4a6d6e7 100644 (file)
--- a/source/cluster/wham/src/readrtns.F
+++ b/source/cluster/wham/src/readrtns.F
@@ -30,12 +30,14 @@ C
       refstr=(index(controlcard,'REFSTR').gt.0)
       write (iout,*) "REFSTR",refstr
       pdbref=(index(controlcard,'PDBREF').gt.0)
+      dyn_ss=(index(controlcard,'DYN_SS').gt.0)
       iscode=index(controlcard,'ONE_LETTER')
       tree=(index(controlcard,'MAKE_TREE').gt.0)
       min_var=(index(controlcard,'MINVAR').gt.0)
       plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
       punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
-      call readi(controlcard,'NCUT',ncut,1)
+      call readi(controlcard,'NCUT',ncut,0)
+      call readi(controlcard,'NCLUST',nclust,5)
       call readi(controlcard,'NSTART',nstart,0)
       call readi(controlcard,'NEND',nend,0)
       call reada(controlcard,'ECUT',ecut,10.0d0)
@@ -44,7 +46,8 @@ C
       lgrp=(index(controlcard,'LGRP').gt.0)
       caonly=(index(controlcard,'CA_ONLY').gt.0)
       print_dist=(index(controlcard,'PRINT_DIST').gt.0)
-      call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
+      if (ncut.gt.0) 
+     & call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
       call readi(controlcard,'IOPT',iopt,2) 
       lside = index(controlcard,"SIDE").gt.0
       efree = index(controlcard,"EFREE").gt.0
@@ -126,9 +129,43 @@ C Read weights of the subsequent energy terms.
       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
+      call reada(weightcard,"D0CM",d0cm,3.78d0)
+      call reada(weightcard,"AKCM",akcm,15.1d0)
+      call reada(weightcard,"AKTH",akth,11.0d0)
+      call reada(weightcard,"AKCT",akct,12.0d0)
+      call reada(weightcard,"V1SS",v1ss,-1.08d0)
+      call reada(weightcard,"V2SS",v2ss,7.61d0)
+      call reada(weightcard,"V3SS",v3ss,13.7d0)
+      call reada(weightcard,"EBR",ebr,-5.50D0)
       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)
+      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
+      write (iout,'(/a)') "Disulfide bridge parameters:"
+      write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
+        write (iout,'(a,f10.2)') 'S-S depth: ',ss_depth
+      write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
+      write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
+      write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
+     & ' v3ss:',v3ss
+      write (iout,'(2(a,f10.2))') 'ht:',ht,' eps:', eps(1,1)
       if (wcorr4.gt.0.0d0) wcorr=wcorr4
       weights(1)=wsc
       weights(2)=wscp
@@ -205,6 +242,19 @@ C Convert sequence to numeric code
       do i=1,nres
         itype(i)=rescode(i,sequence(i),iscode)
       enddo
+      if (itype(2).eq.10.and.itype(1).eq.ntyp1) then
+        write (iout,*)
+     &   "Glycine is the first full residue, initial dummy deleted"
+        do i=1,nres
+          itype(i)=itype(i+1)
+        enddo
+        nres=nres-1
+      endif
+      if (itype(nres-1).eq.10.and.itype(nres).eq.ntyp1) then
+        write (iout,*)
+     &   "Glycine is the last full residue, terminal dummy deleted"
+        nres=nres-1
+      endif
       print *,nres
       print '(20i4)',(itype(i),i=1,nres)
 
@@ -327,7 +377,7 @@ c      endif
             enddo
           enddo
         endif
-        call contact(.true.,ncont_ref,icont_ref)
+        call contact(print_contact_map,ncont_ref,icont_ref)
       endif
 c Read distance restraints
       if (constr_dist.gt.0) then
@@ -431,6 +481,22 @@ c          forcon(i)=fbr
         enddo
       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.
+c          write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU"
+          enddo
+      endif
+      print *, "Leaving brigde read"
       return
       end
 c----------------------------------------------------------------------------