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