source/cluster/wham/src-M/readpdb.f

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