source/wham/src/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       character*3 seq,atom,res
  16       character*80 card
  17       double precision sccor(3,20)
  18       integer i,j,iii,ibeg,ishift,ishift1,ity,ires,ires_old
  19       double precision dcj
  20       integer rescode
  21       ibeg=1
  22       ishift1=0
  23       do i=1,10000
  24         read (ipdbin,'(a80)',end=10) card
  25         if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10
  26 C Fish out the ATOM cards.
  27         if (index(card(1:4),'ATOM').gt.0) then
  28           read (card(14:16),'(a3)') atom
  29           if (atom.eq.'CA' .or. atom.eq.'CH3') then
  30 C Calculate the CM of the preceding residue.
  31             if (ibeg.eq.0) call sccenter(ires,iii,sccor)
  32 C Start new residue.
  33             ires_old=ires+ishift-ishift1
  34             read (card(23:26),*) ires
  35 c            print *,"ires_old",ires_old," ires",ires
  36             if (card(27:27).eq."A" .or. card(27:27).eq."B") then
  37 c              ishift1=ishift1+1
  38             endif
  39             read (card(18:20),'(a3)') res
  40             if (ibeg.eq.1) then
  41               ishift=ires-1
  42               if (res.ne.'GLY' .and. res.ne. 'ACE') then
  43                 ishift=ishift-1
  44                 itype(1)=21
  45               endif
  46               ibeg=0
  47             else
  48               ishift=ishift+ires-ires_old-1
  49             endif
  50             ires=ires-ishift+ishift1
  51             if (res.eq.'ACE') then
  52               ity=10
  53             else
  54               itype(ires)=rescode(ires,res,0)
  55             endif
  56             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
  57             write (iout,'(2i3,2x,a,3f8.3)')
  58      &      ires,itype(ires),res,(c(j,ires),j=1,3)
  59             iii=1
  60             do j=1,3
  61               sccor(j,iii)=c(j,ires)
  62             enddo
  63 c            write (*,*) card(23:27),ires,itype(ires)
  64           else if (atom.ne.'O  '.and.atom(1:1).ne.'H' .and.
  65      &             atom.ne.'N  ' .and. atom.ne.'C   ') then
  66             iii=iii+1
  67             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
  68           endif
  69         endif
  70       enddo
  71    10 write (iout,'(a,i5)') ' Nres: ',ires
  72 C Calculate the CM of the last side chain.
  73       call sccenter(ires,iii,sccor)
  74       nres=ires
  75       nsup=nres
  76       nstart_sup=1
  77       if (itype(nres).ne.10) then
  78         nres=nres+1
  79         itype(nres)=21
  80         do j=1,3
  81           dcj=c(j,nres-2)-c(j,nres-3)
  82           c(j,nres)=c(j,nres-1)+dcj
  83           c(j,2*nres)=c(j,nres)
  84         enddo
  85       endif
  86       do i=2,nres-1
  87         do j=1,3
  88           c(j,i+nres)=dc(j,i)
  89         enddo
  90       enddo
  91       do j=1,3
  92         c(j,nres+1)=c(j,1)
  93         c(j,2*nres)=c(j,nres)
  94       enddo
  95       if (itype(1).eq.21) then
  96         nsup=nsup-1
  97         nstart_sup=2
  98         do j=1,3
  99           dcj=c(j,4)-c(j,3)
 100           c(j,1)=c(j,2)-dcj
 101           c(j,nres+1)=c(j,1)
 102         enddo
 103       endif
 104 C Copy the coordinates to reference coordinates
 105       do i=1,2*nres
 106         do j=1,3
 107           cref(j,i)=c(j,i)
 108         enddo
 109       enddo
 110 C Calculate internal coordinates.
 111       do ires=1,nres
 112         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
 113      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
 114      &    (c(j,ires+nres),j=1,3)
 115       enddo
 116       call flush(iout)
 117       call int_from_cart(.true.,.true.)
 118       do i=1,nres
 119         phi_ref(i)=phi(i)
 120         theta_ref(i)=theta(i)
 121         alph_ref(i)=alph(i)
 122         omeg_ref(i)=omeg(i)
 123       enddo
 124       ishift_pdb=ishift
 125       return
 126       end
 127 c---------------------------------------------------------------------------
 128       subroutine int_from_cart(lside,lprn)
 129       implicit none
 130       include 'DIMENSIONS'
 131       include 'DIMENSIONS.ZSCOPT'
 132       include 'COMMON.LOCAL'
 133       include 'COMMON.VAR'
 134       include 'COMMON.CHAIN'
 135       include 'COMMON.INTERACT'
 136       include 'COMMON.IOUNITS'
 137       include 'COMMON.GEO'
 138       include 'COMMON.NAMES'
 139       character*3 seq,atom,res
 140       character*80 card
 141       double precision sccor(3,20)
 142       integer rescode
 143       double precision dist,alpha,beta,di
 144       integer i,j,iti
 145       logical lside,lprn
 146       if (lprn) then
 147         write (iout,'(/a)')
 148      &  'Internal coordinates calculated from crystal structure.'
 149         if (lside) then
 150           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
 151      & '       Phi','    Dsc_id','       Dsc','     Alpha',
 152      & '     Omega'
 153         else
 154           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
 155      & '       Phi'
 156         endif
 157       endif
 158       do i=2,nres
 159         iti=itype(i)
 160         write (iout,*) i,i-1,(c(j,i),j=1,3),(c(j,i-1),j=1,3),dist(i,i-1)
 161         if (itype(i-1).ne.21 .and. itype(i).ne.21 .and.
 162      &    (dist(i,i-1).lt.2.0D0 .or. dist(i,i-1).gt.5.0D0)) then
 163           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 164           stop
 165         endif
 166         theta(i+1)=alpha(i-1,i,i+1)
 167         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
 168       enddo
 169       if (itype(1).eq.21) then
 170         do j=1,3
 171           c(j,1)=c(j,2)+(c(j,3)-c(j,4))
 172         enddo
 173       endif
 174       if (itype(nres).eq.21) then
 175         do j=1,3
 176           c(j,nres)=c(j,nres-1)+(c(j,nres-2)-c(j,nres-3))
 177         enddo
 178       endif
 179       if (lside) then
 180         do i=2,nres-1
 181           do j=1,3
 182             c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
 183           enddo
 184           iti=itype(i)
 185           di=dist(i,nres+i)
 186           if (iti.ne.10) then
 187             alph(i)=alpha(nres+i,i,maxres2)
 188             omeg(i)=beta(nres+i,i,maxres2,i+1)
 189           endif
 190           if (lprn)
 191      &    write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 192      &    rad2deg*theta(i),rad2deg*phi(i),dsc(iti),di,
 193      &    rad2deg*alph(i),rad2deg*omeg(i)
 194         enddo
 195       else if (lprn) then
 196         do i=2,nres
 197           iti=itype(i)
 198           write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 199      &    rad2deg*theta(i),rad2deg*phi(i)
 200         enddo
 201       endif
 202       return
 203       end
 204 c---------------------------------------------------------------------------
 205       subroutine sccenter(ires,nscat,sccor)
 206       implicit none
 207       include 'DIMENSIONS'
 208       include 'COMMON.CHAIN'
 209       integer ires,nscat,i,j
 210       double precision sccor(3,20),sccmj
 211       do j=1,3
 212         sccmj=0.0D0
 213         do i=1,nscat
 214           sccmj=sccmj+sccor(j,i)
 215         enddo
 216         dc(j,ires)=sccmj/nscat
 217       enddo
 218       return
 219       end