source/wham/src-M-homology/read_dist_constr.F

   1       subroutine read_dist_constr
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4       include 'DIMENSIONS.FREE'
   5       include 'COMMON.CONTROL'
   6       include 'COMMON.CHAIN'
   7       include 'COMMON.IOUNITS'
   8       include 'COMMON.SBRIDGE'
   9       integer ifrag_(2,100),ipair_(2,100)
  10       double precision wfrag_(100),wpair_(100)
  11       character*500 controlcard
  12       logical normalize,next
  13       integer restr_type
  14       double precision xlink(4,0:4) /
  15 c           a          b       c     sigma
  16      &   0.0d0,0.0d0,0.0d0,0.0d0,                             ! default, no xlink potential
  17      &   0.00305218d0,9.46638d0,4.68901d0,4.74347d0,          ! ZL
  18      &   0.00214928d0,12.7517d0,0.00375009d0,6.13477d0,       ! ADH
  19      &   0.00184547d0,11.2678d0,0.00140292d0,7.00868d0,       ! PDH
  20      &   0.000161786d0,6.29273d0,4.40993d0,7.13956d0    /     ! DSS
  21       write (iout,*) "Calling read_dist_constr"
  22 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
  23 c      call flush(iout)
  24       next=.true.
  25
  26       DO WHILE (next)
  27
  28       call card_concat(controlcard,.true.)
  29       next = index(controlcard,"NEXT").gt.0
  30       call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
  31       write (iout,*) "restr_type",restr_type
  32       call readi(controlcard,"NFRAG",nfrag_,0)
  33       call readi(controlcard,"NPAIR",npair_,0)
  34       call readi(controlcard,"NDIST",ndist_,0)
  35       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
  36       if (restr_type.eq.10)
  37      &  call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
  38       if (restr_type.eq.12)
  39      &  call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
  40       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
  41       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
  42       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
  43       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
  44       normalize = index(controlcard,"NORMALIZE").gt.0
  45       write (iout,*) "WBOLTZD",wboltzd
  46       write (iout,*) "SCAL_PEAK",scal_peak
  47 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
  48 c      write (iout,*) "IFRAG"
  49 c      do i=1,nfrag_
  50 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
  51 c      enddo
  52 c      write (iout,*) "IPAIR"
  53 c      do i=1,npair_
  54 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
  55 c      enddo
  56       if (nfrag_.gt.0) then
  57         nres0=nres
  58         read(inp,'(a)') pdbfile
  59         write (iout,*)
  60      & "Distance restraints will be constructed from structure ",pdbfile
  61         open(ipdbin,file=pdbfile,status='old',err=11)
  62         call readpdb(.true.)
  63         nres=nres0
  64         close(ipdbin)
  65       endif
  66       do i=1,nfrag_
  67         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
  68         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
  69      &    ifrag_(2,i)=nstart_sup+nsup-1
  70 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
  71 c        call flush(iout)
  72         if (wfrag_(i).eq.0.0d0) cycle
  73         do j=ifrag_(1,i),ifrag_(2,i)-1
  74           do k=j+1,ifrag_(2,i)
  75 c            write (iout,*) "j",j," k",k
  76             ddjk=dist(j,k)
  77             if (restr_type.eq.1) then
  78               nhpb=nhpb+1
  79               irestr_type(nhpb)=1
  80               ihpb(nhpb)=j
  81               jhpb(nhpb)=k
  82               dhpb(nhpb)=ddjk
  83               forcon(nhpb)=wfrag_(i)
  84             else if (constr_dist.eq.2) then
  85               if (ddjk.le.dist_cut) then
  86                 nhpb=nhpb+1
  87                 irestr_type(nhpb)=1
  88                 ihpb(nhpb)=j
  89                 jhpb(nhpb)=k
  90                 dhpb(nhpb)=ddjk
  91                 forcon(nhpb)=wfrag_(i)
  92               endif
  93             else if (restr_type.eq.3) then
  94               nhpb=nhpb+1
  95               irestr_type(nhpb)=1
  96               ihpb(nhpb)=j
  97               jhpb(nhpb)=k
  98               dhpb(nhpb)=ddjk
  99               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 100             endif
 101             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
 102      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 103           enddo
 104         enddo
 105       enddo
 106       do i=1,npair_
 107         if (wpair_(i).eq.0.0d0) cycle
 108         ii = ipair_(1,i)
 109         jj = ipair_(2,i)
 110         if (ii.gt.jj) then
 111           itemp=ii
 112           ii=jj
 113           jj=itemp
 114         endif
 115         do j=ifrag_(1,ii),ifrag_(2,ii)
 116           do k=ifrag_(1,jj),ifrag_(2,jj)
 117             if (restr_type.eq.1) then
 118               nhpb=nhpb+1
 119               irestr_type(nhpb)=1
 120               ihpb(nhpb)=j
 121               jhpb(nhpb)=k
 122               dhpb(nhpb)=ddjk
 123               forcon(nhpb)=wfrag_(i)
 124             else if (constr_dist.eq.2) then
 125               if (ddjk.le.dist_cut) then
 126                 nhpb=nhpb+1
 127                 irestr_type(nhpb)=1
 128                 ihpb(nhpb)=j
 129                 jhpb(nhpb)=k
 130                 dhpb(nhpb)=ddjk
 131                 forcon(nhpb)=wfrag_(i)
 132               endif
 133             else if (restr_type.eq.3) then
 134               nhpb=nhpb+1
 135               irestr_type(nhpb)=1
 136               ihpb(nhpb)=j
 137               jhpb(nhpb)=k
 138               dhpb(nhpb)=ddjk
 139               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 140             endif
 141             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
 142      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 143           enddo
 144         enddo
 145       enddo
 146
 147 c      print *,ndist_
 148       write (iout,*) "Distance restraints as read from input"
 149       do i=1,ndist_
 150         if (restr_type.eq.12) then
 151           read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
 152      &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
 153      &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
 154      &    fordepth_peak(nhpb_peak+1),npeak
 155 c          write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
 156 c     &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
 157 c     &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
 158 c     &    fordepth_peak(nhpb_peak+1),npeak
 159           if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
 160      &      fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
 161           nhpb_peak=nhpb_peak+1
 162           irestr_type_peak(nhpb_peak)=12
 163           if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
 164           ipeak(2,npeak)=i
 165           write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
 166      &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
 167      &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
 168      &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
 169      &     fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
 170           if (ibecarb_peak(nhpb_peak).gt.0) then
 171             ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
 172             jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
 173           endif
 174         else if (restr_type.eq.11) then
 175           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
 176      &     dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
 177 c        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
 178           if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
 179           nhpb=nhpb+1
 180           irestr_type(nhpb)=11
 181           write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
 182      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
 183      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
 184           if (ibecarb(nhpb).gt.0) then
 185             ihpb(nhpb)=ihpb(nhpb)+nres
 186             jhpb(nhpb)=jhpb(nhpb)+nres
 187           endif
 188         else if (constr_dist.eq.10) then
 189 c Cross-lonk Markov-like potential
 190           call card_concat(controlcard,.true.)
 191           call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
 192           call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
 193           ibecarb(nhpb+1)=0
 194           if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
 195           if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
 196           if (index(controlcard,"ZL").gt.0) then
 197             link_type=1
 198           else if (index(controlcard,"ADH").gt.0) then
 199             link_type=2
 200           else if (index(controlcard,"PDH").gt.0) then
 201             link_type=3
 202           else if (index(controlcard,"DSS").gt.0) then
 203             link_type=4
 204           else
 205             link_type=0
 206           endif
 207           call reada(controlcard,"AXLINK",dhpb(nhpb+1),
 208      &       xlink(1,link_type))
 209           call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
 210      &       xlink(2,link_type))
 211           call reada(controlcard,"CXLINK",fordepth(nhpb+1),
 212      &       xlink(3,link_type))
 213           call reada(controlcard,"SIGMA",forcon(nhpb+1),
 214      &       xlink(4,link_type))
 215           call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
 216 c          read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
 217 c     &      dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
 218           if (forcon(nhpb+1).le.0.0d0 .or.
 219      &       (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
 220           nhpb=nhpb+1
 221           irestr_type(nhpb)=10
 222           if (ibecarb(nhpb).gt.0) then
 223             ihpb(nhpb)=ihpb(nhpb)+nres
 224             jhpb(nhpb)=jhpb(nhpb)+nres
 225           endif
 226           write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
 227      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
 228      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
 229      &     irestr_type(nhpb)
 230         else
 231 C        print *,"in else"
 232           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
 233      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
 234           if (forcon(nhpb+1).gt.0.0d0) then
 235           nhpb=nhpb+1
 236           if (dhpb1(nhpb).eq.0.0d0) then
 237             irestr_type(nhpb)=1
 238           else
 239             irestr_type(nhpb)=2
 240           endif
 241           if (ibecarb(nhpb).gt.0) then
 242             ihpb(nhpb)=ihpb(nhpb)+nres
 243             jhpb(nhpb)=jhpb(nhpb)+nres
 244           endif
 245           if (dhpb(nhpb).eq.0.0d0)
 246      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 247         endif
 248         write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
 249      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
 250         endif
 251 C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
 252 C        if (forcon(nhpb+1).gt.0.0d0) then
 253 C          nhpb=nhpb+1
 254 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 255       enddo
 256
 257       ENDDO ! next
 258
 259       fordepthmax=0.0d0
 260       if (normalize) then
 261         do i=nss+1,nhpb
 262           if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
 263      &      fordepthmax=fordepth(i)
 264         enddo
 265         do i=nss+1,nhpb
 266           if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
 267         enddo
 268       endif
 269       if (nhpb.gt.nss)  then
 270         write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
 271      &  "The following",nhpb-nss,
 272      &  " distance restraints have been imposed:",
 273      &  "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V",
 274      &  "  score"," type"
 275         do i=nss+1,nhpb
 276           write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
 277      &  ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
 278      &  irestr_type(i)
 279         enddo
 280       endif
 281       write (iout,*)
 282       call hpb_partition
 283       call flush(iout)
 284       return
 285    11 write (iout,*)"read_dist_restr: error reading reference structure"
 286       stop
 287       end