source/maxlik/src_FPy/readpdb.F

   1       subroutine readpdb(lprn,iprot,efree_temp,*)
   2 C Read the PDB file and convert the peptide geometry into virtual-chain
   3 C geometry.
   4       implicit none
   5       include 'DIMENSIONS'
   6       include 'DIMENSIONS.ZSCOPT'
   7       include 'COMMON.CONTROL'
   8       include 'COMMON.LOCAL'
   9       include 'COMMON.VAR'
  10       include 'COMMON.CHAIN'
  11       include 'COMMON.INTERACT'
  12       include 'COMMON.IOUNITS'
  13       include 'COMMON.GEO'
  14       include 'COMMON.NAMES'
  15       include 'COMMON.EREF'
  16       integer i,j,iprot,ibeg,ishift1,ires,iires,iii,ires_old,ishift,ity,
  17      &  ishift_pdb
  18       logical lprn,fail
  19       double precision e1(3),e2(3),e3(3)
  20       double precision dcj,efree_temp
  21       character*3 seq,atom,res
  22       character*80 card
  23       double precision sccor(3,20)
  24       integer rescode
  25       efree_temp=0.0d0
  26       ibeg=1
  27       ishift1=0
  28       ishift=0
  29 c      write (iout,*) "READPDB: UNRES_PDB",unres_pdb
  30 c      call flush(iout)
  31       ires=0
  32       ires_old=0
  33       iii=0
  34       read (ipdbin,'(a80)',end=10) card
  35       do while (card(:3).eq.'END' .or. card(:3).eq.'TER')
  36         read (ipdbin,'(a80)',end=10) card
  37       enddo
  38       do i=1,10000
  39 c        write (iout,'(a)') card
  40 c        call flush(iout)
  41 c        if (card(:3).eq.'END') goto 10
  42         if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10
  43 c Read free energy
  44         if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
  45 C Fish out the ATOM cards.
  46         if (index(card(1:4),'ATOM').gt.0) then
  47           read (card(14:16),'(a3)') atom
  48 c          if (atom.eq.'CA' .or. atom.eq.'CH3') then
  49           read (card(23:26),*) ires
  50           read (card(18:20),'(a3)') res
  51 c          write (iout,*) "ires",ires,ires-ishift+ishift1,
  52 c     &      " ires_old",ires_old
  53 c          write (iout,*) "ishift",ishift," ishift1",ishift1
  54 c          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
  55 c          call flush(iout)
  56           iires=ires-ishift+ishift1
  57           if (iires.ne.ires_old) then
  58 C Calculate the CM of the preceding residue.
  59             if (ibeg.eq.0) then
  60 c              write (iout,*) "Calculating sidechain center iii",iii
  61 c              call flush(iout)
  62               if (unres_pdb) then
  63                 do j=1,3
  64                   dc(j,ires)=sccor(j,iii)
  65                 enddo
  66               else
  67                 call sccenter(iires-1,iii,sccor)
  68               endif
  69 c              write (iout,*) "sidechain center position calculated"
  70 c              call flush(iout)
  71               iii=0
  72             endif
  73 C Start new residue.
  74             if (res.eq.'Cl-' .or. res.eq.'Na+') then
  75               ires=ires_old
  76               read (ipdbin,'(a80)',end=10) card
  77               cycle
  78             else if (ibeg.eq.1) then
  79 c              write (iout,*) "BEG ires",ires," iires",iires
  80 c              call flush(iout)
  81               ishift=ires-1
  82               if (res.ne.'GLY' .and. res.ne. 'ACE') then
  83                 ishift=ishift-1
  84                 itype(1)=ntyp1
  85               endif
  86               ires=ires-ishift+ishift1
  87               ires_old=ires
  88 c              write (iout,*) "ishift",ishift," ires",ires,
  89 c     &         " iires",iires," ires_old",ires_old
  90 c              call flush(iout)
  91               ibeg=0
  92             else
  93               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
  94               ires=ires-ishift+ishift1
  95               ires_old=ires
  96 c              write (iout,*) "RESSHIFT",ires,iires,ishift,ishift1
  97 c              call flush(iout)
  98             endif
  99 c            write (iout,*) "ires",ires," iires",iires
 100 c            call flush(iout)
 101             if (res.eq.'ACE' .or. res.eq.'NHE') then
 102               itype(ires)=10
 103             else
 104               itype(ires)=rescode(ires,res,0)
 105             endif
 106           else
 107             ires=ires-ishift+ishift1
 108           endif
 109 c          write (iout,*) "ires_old",ires_old," ires",ires," iires",iires
 110 c          call flush(iout)
 111           if (card(27:27).eq."A" .or. card(27:27).eq."B") then
 112 c            ishift1=ishift1+1
 113           endif
 114 c          write (2,*) "iires",iires," res ",res," ity",ity
 115           if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
 116      &       res.eq.'NHE'.and.atom(:2).eq.'HN') then
 117             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
 118 #ifdef DEBUG
 119             write (iout,'(2i3,2x,a,3f8.3)')
 120      &      iires,itype(iires),res,(c(j,ires),j=1,3)
 121 #endif
 122             iii=iii+1
 123             do j=1,3
 124               sccor(j,iii)=c(j,ires)
 125             enddo
 126 c            write (*,*) card(23:27),ires,itype(ires)
 127           else if (atom.ne.'O  '.and.atom(1:1).ne.'H' .and.
 128      &             atom.ne.'N  ' .and. atom.ne.'C   ' .and.
 129      &             atom.ne.'OXT') then
 130             iii=iii+1
 131             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
 132           endif
 133         endif
 134         read (ipdbin,'(a80)',end=10) card
 135       enddo
 136    10 if (lprn) write (iout,'(a,i5)') ' Number of residues found: ',ires
 137       if (ires.eq.0) return1
 138 C Calculate the CM of the last side chain.
 139       if (iii.gt.0)  then
 140       if (unres_pdb) then
 141         do j=1,3
 142           dc(j,ires)=sccor(j,iii)
 143         enddo
 144       else
 145         call sccenter(iires,iii,sccor)
 146       endif
 147       endif
 148       nres=ires
 149       nsup(iprot)=nres
 150       nstart_sup(iprot)=1
 151       if (itype(nres).ne.10) then
 152         nres=nres+1
 153         itype(nres)=ntyp1
 154         do j=1,3
 155           dcj=c(j,nres-2)-c(j,nres-3)
 156           c(j,nres)=c(j,nres-1)+dcj
 157           c(j,2*nres)=c(j,nres)
 158         enddo
 159       endif
 160       do i=2,nres-1
 161         do j=1,3
 162           c(j,i+nres)=dc(j,i)
 163         enddo
 164       enddo
 165       do j=1,3
 166         c(j,nres+1)=c(j,1)
 167         c(j,2*nres)=c(j,nres)
 168       enddo
 169       if (itype(1).eq.ntyp1) then
 170         nsup(iprot)=nsup(iprot)-1
 171         nstart_sup(iprot)=2
 172         if (unres_pdb) then
 173 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 174           call refsys(2,3,4,e1,e2,e3,fail)
 175           if (fail) then
 176             e2(1)=0.0d0
 177             e2(2)=1.0d0
 178             e2(3)=0.0d0
 179           endif
 180           do j=1,3
 181             c(j,1)=c(j,2)-3.8d0*e2(j)
 182           enddo
 183         else
 184         do j=1,3
 185           dcj=c(j,4)-c(j,3)
 186           c(j,1)=c(j,2)-dcj
 187           c(j,nres+1)=c(j,1)
 188         enddo
 189         endif
 190       endif
 191 C Copy the coordinates to reference coordinates
 192 c      do i=1,2*nres
 193 c        do j=1,3
 194 c          cref(j,i,iprot)=c(j,i)
 195 c        enddo
 196 c      enddo
 197 C Calculate internal coordinates.
 198       if (lprn) then
 199       write (iout,'(/a)')
 200      &  "Cartesian coordinates of the reference structure"
 201       write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
 202      & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
 203       do ires=1,nres
 204         write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
 205      &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
 206      &    (c(j,ires+nres),j=1,3)
 207       enddo
 208       endif
 209       call int_from_cart(.true.,lprn)
 210 c      do i=1,nres
 211 c        phi_ref(i,iprot)=phi(i)
 212 c        theta_ref(i,iprot)=theta(i)
 213 c        alph_ref(i,iprot)=alph(i)
 214 c        omeg_ref(i,iprot)=omeg(i)
 215 c      enddo
 216       ishift_pdb=ishift
 217       return
 218       end
 219 c---------------------------------------------------------------------------
 220       subroutine int_from_cart(lside,lprn)
 221       implicit none
 222       include 'DIMENSIONS'
 223       include 'DIMENSIONS.ZSCOPT'
 224       include 'COMMON.LOCAL'
 225       include 'COMMON.VAR'
 226       include 'COMMON.CHAIN'
 227       include 'COMMON.INTERACT'
 228       include 'COMMON.IOUNITS'
 229       include 'COMMON.GEO'
 230       include 'COMMON.NAMES'
 231       integer i,j,iti
 232       double precision dist,alpha,beta,di
 233       character*3 seq,atom,res
 234       character*80 card
 235       double precision sccor(3,20)
 236       integer rescode
 237       logical lside,lprn
 238       if (lprn) then
 239         write (iout,'(/a)')
 240      &  'Internal coordinates of the reference structure.'
 241         if (lside) then
 242           write (iout,'(8a)') 'Residue','       dvb','     Theta',
 243      & '       Phi','    Dsc_id','       Dsc','     Alpha',
 244      & '     Omega'
 245         else
 246           write (iout,'(4a)') 'Residue','       dvb','     Theta',
 247      & '       Phi'
 248         endif
 249       endif
 250       do i=1,nres-1
 251         iti=itype(i)
 252         if (iti.ne.ntyp1 .and. itype(i+1).ne.ntyp1 .and.
 253      &    (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0)) then
 254           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 255           stop
 256         endif
 257         vbld(i+1)=dist(i,i+1)
 258         vbld_inv(i+1)=1.0d0/vbld(i+1)
 259         if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
 260         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
 261         do j=1,3
 262           dc(j,i)=c(j,i+1)-c(j,i)
 263           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 264         enddo
 265       enddo
 266       if (lside) then
 267         do i=2,nct
 268           do j=1,3
 269             c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
 270           enddo
 271           iti=itype(i)
 272           di=dist(i,nres+i)
 273           vbld(i+nres)=di
 274           if (itype(i).ne.10) then
 275             vbld_inv(i+nres)=1.0d0/di
 276             do j=1,3
 277               dc(j,i+nres)=c(j,i+nres)-c(j,i)
 278               dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 279             enddo
 280           else
 281             vbld_inv(i+nres)=0.0d0
 282             do j=1,3
 283               dc(j,i+nres)=0.0d0
 284               dc_norm(j,i+nres)=0.0d0
 285             enddo
 286           endif
 287           if (iti.ne.10) then
 288             alph(i)=alpha(nres+i,i,maxres2)
 289             omeg(i)=beta(nres+i,i,maxres2,i+1)
 290           endif
 291           if (lprn)
 292      &    write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 293      &    rad2deg*theta(i),rad2deg*phi(i),dsc(iti),di,
 294      &    rad2deg*alph(i),rad2deg*omeg(i)
 295           if (iti.ne.10) then
 296           endif
 297         enddo
 298       else if (lprn) then
 299         do i=2,nres
 300           iti=itype(i)
 301           write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 302      &    rad2deg*theta(i),rad2deg*phi(i)
 303         enddo
 304       endif
 305       return
 306       end
 307 c---------------------------------------------------------------------------
 308       subroutine sccenter(ires,nscat,sccor)
 309       implicit none
 310       include 'DIMENSIONS'
 311       include 'DIMENSIONS.ZSCOPT'
 312       include 'COMMON.CHAIN'
 313       integer ires,nscat,i,j
 314       double precision sccmj
 315       double precision sccor(3,20)
 316       do j=1,3
 317         sccmj=0.0D0
 318         do i=1,nscat
 319           sccmj=sccmj+sccor(j,i)
 320         enddo
 321         dc(j,ires)=sccmj/nscat
 322       enddo
 323       return
 324       end