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)=ntyp1
  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)=ntyp1
  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.ntyp1) 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)=ntyp1
 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.ntyp1) 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 c      do i=1,2*nres
 148 c        do j=1,3
 149 c          cref(j,i)=c(j,i)
 150 c        enddo
 151 c      enddo
 152
 153       kkk=1
 154       lll=0
 155       cou=1
 156       do i=1,nres
 157       lll=lll+1
 158 cc      write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
 159       if (i.gt.1) then
 160       if ((itype(i-1).eq.ntyp1).and.(i.gt.2)) then
 161       chain_length=lll-1
 162       kkk=kkk+1
 163 c       write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
 164       lll=1
 165       endif
 166       endif
 167         do j=1,3
 168           cref(j,i,cou)=c(j,i)
 169           cref(j,i+nres,cou)=c(j,i+nres)
 170           if (i.le.nres) then
 171           chain_rep(j,lll,kkk)=c(j,i)
 172           chain_rep(j,lll+nres,kkk)=c(j,i+nres)
 173           endif
 174          enddo
 175       enddo
 176       do j=1,3
 177       chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
 178       chain_rep(j,chain_length+nres,symetr)
 179      &=chain_rep(j,chain_length+nres,1)
 180       enddo
 181
 182       if (symetr.gt.1) then
 183        call permut(symetr)
 184        nperm=1
 185        do i=1,symetr
 186        nperm=nperm*i
 187        enddo
 188 c       do i=1,nperm
 189 c       write(iout,*) "tabperm", (tabperm(i,kkk),kkk=1,4)
 190 c       enddo
 191        do i=1,nperm
 192         cou=0
 193         do kkk=1,symetr
 194          icha=tabperm(i,kkk)
 195 c         write (iout,*) i,icha
 196          do lll=1,chain_length
 197           cou=cou+1
 198            if (cou.le.nres) then
 199            do j=1,3
 200             kupa=mod(lll,chain_length)
 201             iprzes=(kkk-1)*chain_length+lll
 202             if (kupa.eq.0) kupa=chain_length
 203 c            write (iout,*) "kupa", kupa
 204             cref(j,iprzes,i)=chain_rep(j,kupa,icha)
 205             cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
 206           enddo
 207           endif
 208          enddo
 209         enddo
 210        enddo
 211        endif
 212
 213 C-koniec robienia kopidm
 214       do kkk=1,nperm
 215       write (iout,*) "nowa struktura", nperm
 216       do i=1,nres
 217         write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),
 218      &cref(2,i,kkk),
 219      &cref(3,i,kkk),cref(1,nres+i,kkk),
 220      &cref(2,nres+i,kkk),cref(3,nres+i,kkk)
 221       enddo
 222   100 format (//'              alpha-carbon coordinates       ',
 223      &          '     centroid coordinates'/
 224      1          '       ', 6X,'X',11X,'Y',11X,'Z',
 225      &                          10X,'X',11X,'Y',11X,'Z')
 226   110 format (a,'(',i3,')',6f12.5)
 227        enddo
 228
 229
 230       ishift_pdb=ishift
 231       return
 232       end
 233 c---------------------------------------------------------------------------
 234       subroutine int_from_cart(lside,lprn)
 235       implicit real*8 (a-h,o-z)
 236       include 'DIMENSIONS'
 237       include 'COMMON.LOCAL'
 238       include 'COMMON.VAR'
 239       include 'COMMON.CHAIN'
 240       include 'COMMON.INTERACT'
 241       include 'COMMON.IOUNITS'
 242       include 'COMMON.GEO'
 243       include 'COMMON.NAMES'
 244       character*3 seq,atom,res
 245       character*80 card
 246       dimension sccor(3,20)
 247       integer rescode
 248       logical lside,lprn
 249       if (lprn) then
 250         write (iout,'(/a)')
 251      &  'Internal coordinates calculated from crystal structure.'
 252         if (lside) then
 253           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
 254      & '       Phi','    Dsc_id','       Dsc','     Alpha',
 255      & '     Omega'
 256         else
 257           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
 258      & '       Phi'
 259         endif
 260       endif
 261       call flush(iout)
 262       do i=nnt+1,nct
 263         iti=itype(i)
 264 c        write (iout,*) i,dist(i,i-1)
 265         if (dist(i,i-1).lt.2.0D0 .or. dist(i,i-1).gt.5.0D0) then
 266           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 267           stop
 268         endif
 269         theta(i+1)=alpha(i-1,i,i+1)
 270         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
 271       enddo
 272       if (lside) then
 273         do i=2,nres-1
 274           do j=1,3
 275             c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
 276           enddo
 277           iti=itype(i)
 278           di=dist(i,nres+i)
 279           if (iti.ne.10) then
 280             alph(i)=alpha(nres+i,i,maxres2)
 281             omeg(i)=beta(nres+i,i,maxres2,i+1)
 282           endif
 283           if (lprn)
 284      &    write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 285      &    rad2deg*theta(i),rad2deg*phi(i),dsc(iti),di,rad2deg*alph(i),
 286      &    rad2deg*omeg(i)
 287         enddo
 288       else if (lprn) then
 289         do i=2,nres
 290           iti=itype(i)
 291           write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 292      &    rad2deg*theta(i),rad2deg*phi(i)
 293         enddo
 294       endif
 295       return
 296       end
 297 c---------------------------------------------------------------------------
 298       subroutine sccenter(ires,nscat,sccor)
 299       implicit real*8 (a-h,o-z)
 300       include 'DIMENSIONS'
 301       include 'COMMON.CHAIN'
 302       dimension sccor(3,20)
 303       do j=1,3
 304         sccmj=0.0D0
 305         do i=1,nscat
 306           sccmj=sccmj+sccor(j,i)
 307         enddo
 308         dc(j,ires)=sccmj/nscat
 309       enddo
 310       return
 311       end