source/cluster/wham/src-M-SAXS/readpdb.f.safe

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