source/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 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)=21
  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)=21
  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.21) 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)=21
 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.21) 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.21) 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       do j=1,3
 183       chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
 184       chain_rep(j,chain_length+nres,symetr)
 185      &=chain_rep(j,chain_length+nres,1)
 186       enddo
 187 c diagnostic
 188
 189 c diagnostic
 190 c       write (iout,*) "spraw lancuchy",chain_length,symetr
 191 c       do i=1,24
 192 c         do kkk=1,chain_length
 193 c           write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
 194 c         enddo
 195 c        enddo
 196 c enddiagnostic
 197 C makes copy of chains
 198 c      write (iout,*) "symetr", symetr
 199
 200       if (symetr.gt.1) then
 201        call permut(symetr)
 202        nperm=1
 203        do i=1,symetr
 204        nperm=nperm*i
 205        enddo
 206 c       do i=1,nperm
 207 c       write(iout,*) "tabperm", (tabperm(i,kkk),kkk=1,4)
 208 c       enddo
 209        do i=1,nperm
 210         cou=0
 211         do kkk=1,symetr
 212          icha=tabperm(i,kkk)
 213 c         write (iout,*) i,icha
 214          do lll=1,chain_length
 215           cou=cou+1
 216            if (cou.le.nres) then
 217            do j=1,3
 218             kupa=mod(lll,chain_length)
 219             iprzes=(kkk-1)*chain_length+lll
 220             if (kupa.eq.0) kupa=chain_length
 221 c            write (iout,*) "kupa", kupa
 222             cref(j,iprzes,i)=chain_rep(j,kupa,icha)
 223             cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
 224           enddo
 225           endif
 226          enddo
 227         enddo
 228        enddo
 229        endif
 230 C-koniec robienia kopidm
 231       do kkk=1,nperm
 232       write (iout,*) "nowa struktura", nperm
 233       do i=1,nres
 234         write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),
 235      &cref(2,i,kkk),
 236      &cref(3,i,kkk),cref(1,nres+i,kkk),
 237      &cref(2,nres+i,kkk),cref(3,nres+i,kkk)
 238       enddo
 239   100 format (//'              alpha-carbon coordinates       ',
 240      &          '     centroid coordinates'/
 241      1          '       ', 6X,'X',11X,'Y',11X,'Z',
 242      &                          10X,'X',11X,'Y',11X,'Z')
 243   110 format (a,'(',i3,')',6f12.5)
 244
 245       enddo
 246
 247       ishift_pdb=ishift
 248       return
 249       end
 250 c---------------------------------------------------------------------------
 251       subroutine int_from_cart(lside,lprn)
 252       implicit none
 253       include 'DIMENSIONS'
 254       include 'DIMENSIONS.ZSCOPT'
 255       include 'COMMON.LOCAL'
 256       include 'COMMON.VAR'
 257       include 'COMMON.CHAIN'
 258       include 'COMMON.INTERACT'
 259       include 'COMMON.IOUNITS'
 260       include 'COMMON.GEO'
 261       include 'COMMON.NAMES'
 262       character*3 seq,atom,res
 263       character*80 card
 264       double precision sccor(3,20)
 265       integer rescode
 266       double precision dist,alpha,beta,di
 267       integer i,j,iti
 268       logical lside,lprn
 269       if (lprn) then
 270         write (iout,'(/a)')
 271      &  'Internal coordinates calculated from crystal structure.'
 272         if (lside) then
 273           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
 274      & '       Phi','    Dsc_id','       Dsc','     Alpha',
 275      & '     Omega'
 276         else
 277           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
 278      & '       Phi'
 279         endif
 280       endif
 281       do i=2,nres
 282         iti=itype(i)
 283         write (iout,*) i,i-1,(c(j,i),j=1,3),(c(j,i-1),j=1,3),dist(i,i-1)
 284         if (itype(i-1).ne.21 .and. itype(i).ne.21 .and.
 285      &    (dist(i,i-1).lt.2.0D0 .or. dist(i,i-1).gt.5.0D0)) then
 286           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 287           stop
 288         endif
 289         theta(i+1)=alpha(i-1,i,i+1)
 290         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
 291       enddo
 292       if (itype(1).eq.21) then
 293         do j=1,3
 294           c(j,1)=c(j,2)+(c(j,3)-c(j,4))
 295         enddo
 296       endif
 297       if (itype(nres).eq.21) then
 298         do j=1,3
 299           c(j,nres)=c(j,nres-1)+(c(j,nres-2)-c(j,nres-3))
 300         enddo
 301       endif
 302       if (lside) then
 303         do i=2,nres-1
 304           do j=1,3
 305             c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
 306           enddo
 307           iti=itype(i)
 308           di=dist(i,nres+i)
 309           if (iti.ne.10) then
 310             alph(i)=alpha(nres+i,i,maxres2)
 311             omeg(i)=beta(nres+i,i,maxres2,i+1)
 312           endif
 313           if (lprn)
 314      &    write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 315      &    rad2deg*theta(i),rad2deg*phi(i),dsc(iti),di,
 316      &    rad2deg*alph(i),rad2deg*omeg(i)
 317         enddo
 318       else if (lprn) then
 319         do i=2,nres
 320           iti=itype(i)
 321           write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 322      &    rad2deg*theta(i),rad2deg*phi(i)
 323         enddo
 324       endif
 325       return
 326       end
 327 c---------------------------------------------------------------------------
 328       subroutine sccenter(ires,nscat,sccor)
 329       implicit none
 330       include 'DIMENSIONS'
 331       include 'COMMON.CHAIN'
 332       integer ires,nscat,i,j
 333       double precision sccor(3,20),sccmj
 334       do j=1,3
 335         sccmj=0.0D0
 336         do i=1,nscat
 337           sccmj=sccmj+sccor(j,i)
 338         enddo
 339         dc(j,ires)=sccmj/nscat
 340       enddo
 341       return
 342       end