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