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