source/cluster/wham/src-HCD-5D/readpdb.f.safe

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