source/wham/src-HCD-5D/read_dist_constr.F

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