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