X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc%2Fmolread_zs.F;h=87a34953ecd4fd1559fb7dffbfe4e260c4f08629;hb=24bbf0f3f06452c69bf4b69b4a9ba620e0260d7b;hp=885c57bfa9d730988c43f6d05d32ccb75ca30486;hpb=fbc0f12bd8022b2c766eb8c350e7bb0451816c1d;p=unres.git

diff --git a/source/wham/src/molread_zs.F b/source/wham/src/molread_zs.F
index 885c57b..87a3495 100644
--- a/source/wham/src/molread_zs.F
+++ b/source/wham/src/molread_zs.F
@@ -94,6 +94,12 @@ C Convert sequence to numeric code
       if (itype(1).eq.21) nnt=2
       if (itype(nres).eq.21) nct=nct-1
       write(iout,*) 'NNT=',NNT,' NCT=',NCT
+c Read distance restraints
+      if (constr_dist.gt.0) then
+        call read_dist_constr
+        call hpb_partition
+      endif
+
       call setup_var
       call init_int_table
       if (ns.gt.0) then
@@ -242,3 +248,129 @@ c      print *,"energy",energ," iscor",iscor
       return
    10 return1
       end
+c-------------------------------------------------------------------------------
+      subroutine read_dist_constr
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.CONTROL'
+      include 'COMMON.CHAIN'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.SBRIDGE'
+      integer ifrag_(2,100),ipair_(2,100)
+      double precision wfrag_(100),wpair_(100)
+      character*500 controlcard
+c      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,.true.)
+      call readi(controlcard,"NFRAG",nfrag_,0)
+      call readi(controlcard,"NPAIR",npair_,0)
+      call readi(controlcard,"NDIST",ndist_,0)
+      call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
+      call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
+      call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
+      call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
+      call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
+      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
+      write (iout,*) "IFRAG"
+      do i=1,nfrag_
+        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
+      enddo
+      write (iout,*) "IPAIR"
+      do i=1,npair_
+        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
+      enddo
+      call flush(iout)
+      if (.not.refstr .and. nfrag_.gt.0) then
+        write (iout,*) 
+     &  "ERROR: no reference structure to compute distance restraints"
+        write (iout,*)
+     &  "Restraints must be specified explicitly (NDIST=number)"
+        stop 
+      endif
+      if (nfrag_.lt.2 .and. npair_.gt.0) then 
+        write (iout,*) "ERROR: Less than 2 fragments specified",
+     &   " but distance restraints between pairs requested"
+        stop 
+      endif 
+      call flush(iout)
+      do i=1,nfrag_
+        if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
+        if (ifrag_(2,i).gt.nstart_sup+nsup-1)
+     &    ifrag_(2,i)=nstart_sup+nsup-1
+c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
+        call flush(iout)
+        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
+            ddjk=dist(j,k)
+            if (constr_dist.eq.1) then
+            nhpb=nhpb+1
+            ihpb(nhpb)=j
+            jhpb(nhpb)=k
+              dhpb(nhpb)=ddjk
+            forcon(nhpb)=wfrag_(i) 
+            else if (constr_dist.eq.2) then
+              if (ddjk.le.dist_cut) then
+                nhpb=nhpb+1
+                ihpb(nhpb)=j
+                jhpb(nhpb)=k
+                dhpb(nhpb)=ddjk
+                forcon(nhpb)=wfrag_(i) 
+              endif
+            else
+              nhpb=nhpb+1
+              ihpb(nhpb)=j
+              jhpb(nhpb)=k
+              dhpb(nhpb)=ddjk
+              forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
+            endif
+            write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
+     &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+          enddo
+        enddo
+        endif
+      enddo
+      do i=1,npair_
+        if (wpair_(i).gt.0.0d0) then
+        ii = ipair_(1,i)
+        jj = ipair_(2,i)
+        if (ii.gt.jj) then
+          itemp=ii
+          ii=jj
+          jj=itemp
+        endif
+        do j=ifrag_(1,ii),ifrag_(2,ii)
+          do k=ifrag_(1,jj),ifrag_(2,jj)
+            nhpb=nhpb+1
+            ihpb(nhpb)=j
+            jhpb(nhpb)=k
+            forcon(nhpb)=wpair_(i)
+            dhpb(nhpb)=dist(j,k)
+            write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
+     &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+          enddo
+        enddo
+        endif
+      enddo 
+      do i=1,ndist_
+        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
+     &     ibecarb(i),forcon(nhpb+1)
+        if (forcon(nhpb+1).gt.0.0d0) then
+          nhpb=nhpb+1
+          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
+      enddo
+      do i=1,nhpb
+          write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
+     &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
+      enddo
+      call flush(iout)
+      return
+      end