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