source/wham/src-M-SAXS-homology/readpdb.f

   1       subroutine readpdb
   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.SBRIDGE'
  16       character*3 seq,atom,res
  17       character*80 card
  18       double precision sccor(3,50)
  19       integer i,j,iii,ibeg,ishift,ishift1,ity,ires,ires_old
  20       double precision dcj
  21       integer rescode,kkk,lll,icha,cou,kupa,iprzes
  22       ibeg=1
  23       ishift1=0
  24       do
  25         read (ipdbin,'(a80)',end=10) card
  26         if (card(:3).eq.'END') then
  27           goto 10
  28         else if (card(:3).eq.'TER') then
  29 C End current chain
  30 c          ires_old=ires+1
  31           ires_old=ires+2
  32           itype(ires_old-1)=ntyp1
  33           itype(ires_old)=ntyp1
  34           ibeg=2
  35 c          write (iout,*) "Chain ended",ires,ishift,ires_old
  36           call sccenter(ires,iii,sccor)
  37         endif
  38 C Fish out the ATOM cards.
  39         if (index(card(1:4),'ATOM').gt.0) then
  40           read (card(14:16),'(a3)') atom
  41           if (atom.eq.'CA' .or. atom.eq.'CH3') then
  42 C Calculate the CM of the preceding residue.
  43             if (ibeg.eq.0) then
  44               call sccenter(ires,iii,sccor)
  45             endif
  46 C Start new residue.
  47 c            write (iout,'(a80)') card
  48             read (card(23:26),*) ires
  49             read (card(18:20),'(a3)') res
  50             if (ibeg.eq.1) then
  51               ishift=ires-1
  52               if (res.ne.'GLY' .and. res.ne. 'ACE') then
  53                 ishift=ishift-1
  54                 itype(1)=ntyp1
  55               endif
  56 c              write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
  57               ibeg=0
  58             else if (ibeg.eq.2) then
  59 c Start a new chain
  60               ishift=-ires_old+ires-1
  61 c              write (iout,*) "New chain started",ires,ishift
  62               ibeg=0
  63             endif
  64             ires=ires-ishift
  65 c            write (2,*) "ires",ires," ishift",ishift
  66             if (res.eq.'ACE') then
  67               ity=10
  68             else
  69               itype(ires)=rescode(ires,res,0)
  70             endif
  71             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
  72             read(card(61:66),*) bfac(ires)
  73             write (iout,'(2i3,2x,a,3f8.3,5x,f8.3)')
  74      &       ires,itype(ires),res,(c(j,ires),j=1,3),bfac(ires)
  75             iii=1
  76             do j=1,3
  77               sccor(j,iii)=c(j,ires)
  78             enddo
  79           else if (atom.ne.'O  '.and.atom(1:1).ne.'H' .and.
  80      &             atom(1:1).ne.'Q' .and. atom(1:2).ne.'1H' .and.
  81      &             atom(1:2).ne.'2H' .and. atom(1:2).ne.'3H' .and.
  82      &             atom.ne.'N  ' .and. atom.ne.'C   ') then
  83             iii=iii+1
  84             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
  85           endif
  86         endif
  87       enddo
  88    10 write (iout,'(a,i5)') ' Nres: ',ires
  89 C Calculate dummy residue coordinates inside the "chain" of a multichain
  90 C system
  91       nres=ires
  92       do i=2,nres-1
  93 c        write (iout,*) i,itype(i)
  94
  95         if (itype(i).eq.ntyp1) then
  96          if (itype(i+1).eq.ntyp1) then
  97 C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
  98 C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
  99 C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
 100 C           if (unres_pdb) then
 101 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
 102 C            call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
 103 C            if (fail) then
 104 C              e2(1)=0.0d0
 105 C              e2(2)=1.0d0
 106 C              e2(3)=0.0d0
 107 C            endif !fail
 108 C            do j=1,3
 109 C             c(j,i)=c(j,i-1)-1.9d0*e2(j)
 110 C            enddo
 111 C           else   !unres_pdb
 112            do j=1,3
 113              dcj=(c(j,i-2)-c(j,i-3))/2.0
 114              c(j,i)=c(j,i-1)+dcj
 115              c(j,nres+i)=c(j,i)
 116            enddo
 117 C          endif   !unres_pdb
 118          else     !itype(i+1).eq.ntyp1
 119 C          if (unres_pdb) then
 120 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 121 C            call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
 122 C            if (fail) then
 123 C              e2(1)=0.0d0
 124 C              e2(2)=1.0d0
 125 C              e2(3)=0.0d0
 126 C            endif
 127 C            do j=1,3
 128 C              c(j,i)=c(j,i+1)-1.9d0*e2(j)
 129 C            enddo
 130 C          else !unres_pdb
 131            do j=1,3
 132             dcj=(c(j,i+3)-c(j,i+2))/2.0
 133             c(j,i)=c(j,i+1)-dcj
 134             c(j,nres+i)=c(j,i)
 135            enddo
 136 C          endif !unres_pdb
 137          endif !itype(i+1).eq.ntyp1
 138         endif  !itype.eq.ntyp1
 139       enddo
 140 C Calculate the CM of the last side chain.
 141       call sccenter(ires,iii,sccor)
 142       nsup=nres
 143       nstart_sup=1
 144       if (itype(nres).ne.10) then
 145         nres=nres+1
 146         itype(nres)=ntyp1
 147         do j=1,3
 148           dcj=(c(j,nres-2)-c(j,nres-3))/2.0
 149           c(j,nres)=c(j,nres-1)+dcj
 150           c(j,2*nres)=c(j,nres)
 151         enddo
 152       endif
 153       do i=2,nres-1
 154         do j=1,3
 155           c(j,i+nres)=dc(j,i)
 156         enddo
 157       enddo
 158       do j=1,3
 159         c(j,nres+1)=c(j,1)
 160         c(j,2*nres)=c(j,nres)
 161       enddo
 162       if (itype(1).eq.ntyp1) then
 163         nsup=nsup-1
 164         nstart_sup=2
 165         do j=1,3
 166           dcj=(c(j,4)-c(j,3))/2.0
 167           c(j,1)=c(j,2)-dcj
 168           c(j,nres+1)=c(j,1)
 169         enddo
 170       endif
 171 C Calculate internal coordinates.
 172       do ires=1,nres
 173         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
 174      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
 175      &    (c(j,nres+ires),j=1,3)
 176       enddo
 177       call int_from_cart(.true.,.false.)
 178       write (iout,*) "After int_from_cart"
 179       call flush(iout)
 180       do i=1,nres-1
 181         do j=1,3
 182           dc(j,i)=c(j,i+1)-c(j,i)
 183           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 184         enddo
 185       enddo
 186       do i=2,nres-1
 187         do j=1,3
 188           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 189           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 190         enddo
 191 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
 192 c     &   vbld_inv(i+nres)
 193       enddo
 194 c      call chainbuild
 195 C Copy the coordinates to reference coordinates
 196       do i=1,nres
 197         do j=1,3
 198           cref(j,i)=c(j,i)
 199           cref(j,i+nres)=c(j,i+nres)
 200         enddo
 201       enddo
 202   100 format (//'              alpha-carbon coordinates       ',
 203      &          '     centroid coordinates'/
 204      1          '       ', 6X,'X',11X,'Y',11X,'Z',
 205      &                          10X,'X',11X,'Y',11X,'Z')
 206   110 format (a,'(',i3,')',6f12.5)
 207
 208       ishift_pdb=ishift
 209       return
 210       end
 211 c---------------------------------------------------------------------------
 212       subroutine int_from_cart(lside,lprn)
 213       implicit none
 214       include 'DIMENSIONS'
 215       include 'DIMENSIONS.ZSCOPT'
 216       include 'COMMON.LOCAL'
 217       include 'COMMON.VAR'
 218       include 'COMMON.CHAIN'
 219       include 'COMMON.INTERACT'
 220       include 'COMMON.IOUNITS'
 221       include 'COMMON.GEO'
 222       include 'COMMON.NAMES'
 223       character*3 seq,atom,res
 224       character*80 card
 225       double precision sccor(3,50)
 226       integer rescode
 227       double precision dist,alpha,beta,di
 228       integer i,j,iti
 229       logical lside,lprn
 230       if (lprn) then
 231         write (iout,'(/a)')
 232      &  'Internal coordinates calculated from crystal structure.'
 233         if (lside) then
 234           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
 235      & '       Phi','    Dsc_id','       Dsc','     Alpha',
 236      & '     Omega'
 237         else
 238           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
 239      & '       Phi'
 240         endif
 241       endif
 242       do i=2,nres
 243         iti=itype(i)
 244 c        write (iout,*) i,i-1,(c(j,i),j=1,3),(c(j,i-1),j=1,3),dist(i,i-1)
 245         if (itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1 .and.
 246      &    (dist(i,i-1).lt.1.0D0 .or. dist(i,i-1).gt.6.0D0)) then
 247           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 248           stop
 249         endif
 250         theta(i+1)=alpha(i-1,i,i+1)
 251         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
 252       enddo
 253       if (itype(1).eq.ntyp1) then
 254         do j=1,3
 255           c(j,1)=c(j,2)+(c(j,3)-c(j,4))
 256         enddo
 257       endif
 258       if (itype(nres).eq.ntyp1) then
 259         do j=1,3
 260           c(j,nres)=c(j,nres-1)+(c(j,nres-2)-c(j,nres-3))
 261         enddo
 262       endif
 263       if (lside) then
 264         do i=2,nres-1
 265           do j=1,3
 266             c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
 267           enddo
 268           iti=itype(i)
 269           di=dist(i,nres+i)
 270           if (iti.ne.10) then
 271             alph(i)=alpha(nres+i,i,maxres2)
 272             omeg(i)=beta(nres+i,i,maxres2,i+1)
 273           endif
 274           if (lprn)
 275      &    write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 276      &    rad2deg*theta(i),rad2deg*phi(i),dsc(iti),di,
 277      &    rad2deg*alph(i),rad2deg*omeg(i)
 278         enddo
 279       else if (lprn) then
 280         do i=2,nres
 281           iti=itype(i)
 282           write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 283      &    rad2deg*theta(i),rad2deg*phi(i)
 284         enddo
 285       endif
 286       return
 287       end
 288 c---------------------------------------------------------------------------
 289       subroutine sccenter(ires,nscat,sccor)
 290       implicit none
 291       include 'DIMENSIONS'
 292       include 'COMMON.CHAIN'
 293       integer ires,nscat,i,j
 294       double precision sccor(3,50),sccmj
 295       do j=1,3
 296         sccmj=0.0D0
 297         do i=1,nscat
 298           sccmj=sccmj+sccor(j,i)
 299         enddo
 300         dc(j,ires)=sccmj/nscat
 301       enddo
 302       return
 303       end