source/maxlik/src_MD_T_maxlik-NEWCORR-PMF-5/readpdb.F

   1       subroutine readpdb(lprn,iprot,efree_temp,*)
   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       include 'COMMON.EREF'
  16       integer i,j,iprot,ibeg,ishift1,ires,iires,iii,ires_old,ishift,ity,
  17      &  ishift_pdb,ires_ca
  18       logical lprn,fail
  19       double precision e1(3),e2(3),e3(3)
  20       double precision dcj,efree_temp
  21       character*3 seq,atom,res
  22       character*80 card,prevcard
  23       double precision sccor(3,20)
  24       integer rescode,iterter(maxres)
  25       efree_temp=0.0d0
  26       ibeg=1
  27       ishift1=0
  28       ishift=0
  29 c      write (iout,*) "READPDB: UNRES_PDB",unres_pdb
  30 c      call flush(iout)
  31       ires=0
  32       ires_old=0
  33       iii=0
  34       read (ipdbin,'(a80)',end=10) card
  35       do while (card(:3).eq.'END' .or. card(:3).eq.'TER')
  36         read (ipdbin,'(a80)',end=10) card
  37       enddo
  38       do
  39         if (card(:3).eq.'END' .or.
  40      &    card(:3).eq.'TER'.and.prevcard(:5).eq.'REMARK') goto 10
  41         if (card(:3).eq.'TER') then
  42 C End current chain
  43           ires_old=ires+2
  44           itype(ires_old-1)=ntyp1
  45           iterter(ires_old-1)=1
  46           itype(ires_old)=ntyp1
  47           iterter(ires_old)=1
  48           ibeg=2
  49 c          write (iout,*) "Chain ended",ires,ishift,ires_old
  50           if (unres_pdb) then
  51             do j=1,3
  52               dc(j,ires)=sccor(j,iii)
  53             enddo
  54           else
  55             call sccenter(ires,iii,sccor)
  56           endif
  57         endif
  58         prevcard=card
  59 C Fish out the ATOM cards.
  60         if (index(card(1:4),'ATOM').gt.0) then
  61           read (card(12:16),*) atom
  62 c          write (iout,*) "! ",atom," !",ires
  63 c          if (atom.eq.'CA' .or. atom.eq.'CH3') then
  64           read (card(23:26),*) ires
  65           read (card(18:20),'(a3)') res
  66 c          write (iout,*) "ires",ires,ires-ishift+ishift1,
  67 c     &      " ires_old",ires_old
  68 c          write (iout,*) "ishift",ishift," ishift1",ishift1
  69 c          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
  70           if (ires-ishift+ishift1.ne.ires_old) then
  71 C Calculate the CM of the preceding residue.
  72             if (ibeg.eq.0) then
  73               if (unres_pdb) then
  74                 do j=1,3
  75                   dc(j,ires)=sccor(j,iii)
  76                 enddo
  77               else
  78                 call sccenter(ires_old,iii,sccor)
  79               endif
  80               iii=0
  81             endif
  82 C Start new residue.
  83             if (res.eq.'Cl-' .or. res.eq.'Na+') then
  84               ires=ires_old
  85               cycle
  86             else if (ibeg.eq.1) then
  87 c              write (iout,*) "BEG ires",ires
  88               ishift=ires-1
  89               if (res.ne.'GLY' .and. res.ne. 'ACE') then
  90                 ishift=ishift-1
  91                 itype(1)=ntyp1
  92               endif
  93               ires=ires-ishift+ishift1
  94               ires_old=ires
  95 c              write (iout,*) "ishift",ishift," ires",ires,
  96 c     &         " ires_old",ires_old
  97 c              write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
  98               ibeg=0
  99             else if (ibeg.eq.2) then
 100 c Start a new chain
 101               ishift=-ires_old+ires-1
 102               ires=ires_old+1
 103 c              write (iout,*) "New chain started",ires,ishift
 104               ibeg=0
 105             else
 106               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
 107               ires=ires-ishift+ishift1
 108               ires_old=ires
 109             endif
 110             if (res.eq.'ACE' .or. res.eq.'NHE') then
 111               itype(ires)=10
 112             else
 113               itype(ires)=rescode(ires,res,0,.true.)
 114               if (itype(ires).eq.0) stop
 115             endif
 116           else
 117             ires=ires-ishift+ishift1
 118           endif
 119 c          write (iout,*) "ires_old",ires_old," ires",ires
 120 c          if (card(27:27).eq."A" .or. card(27:27).eq."B") then
 121 c            ishift1=ishift1+1
 122 c          endif
 123 c          write (2,*) "ires",ires," res ",res," ity",ity
 124           if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
 125      &       res.eq.'NHE'.and.atom(:2).eq.'HN') then
 126             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
 127 c            write (iout,*) "backbone ",atom
 128 #ifdef DEBUG
 129             write (iout,'(2i3,2x,a,3f8.3)')
 130      &      ires,itype(ires),res,(c(j,ires),j=1,3)
 131 #endif
 132             iii=iii+1
 133             do j=1,3
 134               sccor(j,iii)=c(j,ires)
 135             enddo
 136             if (ishift.ne.0) then
 137               ires_ca=ires+ishift-ishift1
 138             else
 139               ires_ca=ires
 140             endif
 141 c            write (*,*) card(23:27),ires,itype(ires)
 142           else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.
 143      &             atom.ne.'N' .and. atom.ne.'C' .and.
 144      &             atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.
 145      &             atom.ne.'OXT' .and. atom(:2).ne.'3H') then
 146 c            write (iout,*) "sidechain ",atom
 147             iii=iii+1
 148             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
 149           endif
 150         endif
 151         read (ipdbin,'(a80)',end=10) card
 152       enddo
 153    10 write (iout,'(a,i5)') ' Nres: ',ires
 154       if (ires.eq.0) return1
 155 C Calculate dummy residue coordinates inside the "chain" of a multichain
 156 C system
 157       nres=ires
 158       do i=2,nres-1
 159 c        write (iout,*) i,itype(i),itype(i+1)
 160         if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
 161          if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
 162 C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
 163 C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
 164 C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
 165            if (unres_pdb) then
 166 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
 167             call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
 168             if (fail) then
 169               e2(1)=0.0d0
 170               e2(2)=1.0d0
 171               e2(3)=0.0d0
 172             endif !fail
 173             do j=1,3
 174              c(j,i)=c(j,i-1)-1.9d0*e2(j)
 175             enddo
 176            else   !unres_pdb
 177            do j=1,3
 178              dcj=(c(j,i-2)-c(j,i-3))/2.0
 179             if (dcj.eq.0) dcj=1.23591524223
 180              c(j,i)=c(j,i-1)+dcj
 181              c(j,nres+i)=c(j,i)
 182            enddo
 183           endif   !unres_pdb
 184          else     !itype(i+1).eq.ntyp1
 185           if (unres_pdb) then
 186 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 187             call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
 188             if (fail) then
 189               e2(1)=0.0d0
 190               e2(2)=1.0d0
 191               e2(3)=0.0d0
 192             endif
 193             do j=1,3
 194               c(j,i)=c(j,i+1)-1.9d0*e2(j)
 195             enddo
 196           else !unres_pdb
 197            do j=1,3
 198             dcj=(c(j,i+3)-c(j,i+2))/2.0
 199             if (dcj.eq.0) dcj=1.23591524223
 200             c(j,i)=c(j,i+1)-dcj
 201             c(j,nres+i)=c(j,i)
 202            enddo
 203           endif !unres_pdb
 204          endif !itype(i+1).eq.ntyp1
 205         endif  !itype.eq.ntyp1
 206       enddo
 207 C Calculate the CM of the last side chain.
 208       if (unres_pdb) then
 209         do j=1,3
 210           dc(j,ires)=sccor(j,iii)
 211         enddo
 212       else
 213         call sccenter(ires,iii,sccor)
 214       endif
 215       nsup=nres
 216       nstart_sup=1
 217       if (itype(nres).ne.10) then
 218         nres=nres+1
 219         itype(nres)=ntyp1
 220         if (unres_pdb) then
 221 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
 222           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
 223           if (fail) then
 224             e2(1)=0.0d0
 225             e2(2)=1.0d0
 226             e2(3)=0.0d0
 227           endif
 228           do j=1,3
 229             c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
 230           enddo
 231         else
 232         do j=1,3
 233           dcj=(c(j,nres-2)-c(j,nres-3))/2.0
 234             if (dcj.eq.0) dcj=1.23591524223
 235           c(j,nres)=c(j,nres-1)+dcj
 236           c(j,2*nres)=c(j,nres)
 237         enddo
 238       endif
 239       endif
 240       do i=2,nres-1
 241         do j=1,3
 242           c(j,i+nres)=dc(j,i)
 243         enddo
 244       enddo
 245       do j=1,3
 246         c(j,nres+1)=c(j,1)
 247         c(j,2*nres)=c(j,nres)
 248       enddo
 249       if (itype(1).eq.ntyp1) then
 250         nsup=nsup-1
 251         nstart_sup=2
 252         if (unres_pdb) then
 253 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 254           call refsys(2,3,4,e1,e2,e3,fail)
 255           if (fail) then
 256             e2(1)=0.0d0
 257             e2(2)=1.0d0
 258             e2(3)=0.0d0
 259           endif
 260           do j=1,3
 261             c(j,1)=c(j,2)-1.9d0*e2(j)
 262           enddo
 263         else
 264         do j=1,3
 265           dcj=(c(j,4)-c(j,3))/2.0
 266           c(j,1)=c(j,2)-dcj
 267           c(j,nres+1)=c(j,1)
 268         enddo
 269         endif
 270       endif
 271 C Copy the coordinates to reference coordinates
 272 c      do i=1,2*nres
 273 c        do j=1,3
 274 c          cref(j,i)=c(j,i)
 275 c        enddo
 276 c      enddo
 277 C Calculate internal coordinates.
 278       if (lprn) then
 279       write (iout,'(/a)')
 280      &  "Cartesian coordinates of the reference structure"
 281       write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
 282      & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
 283       do ires=1,nres
 284         write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
 285      &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
 286      &    (c(j,ires+nres),j=1,3)
 287       enddo
 288       endif
 289 C Calculate internal coordinates.
 290       write (iout,'(a)')
 291      &   "Backbone and SC coordinates as read from the PDB"
 292       do ires=1,nres
 293         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
 294      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
 295      &    (c(j,nres+ires),j=1,3)
 296       enddo
 297       call int_from_cart(.true.,.false.)
 298       call sc_loc_geom(.false.)
 299       do i=1,nres-1
 300         do j=1,3
 301           dc(j,i)=c(j,i+1)-c(j,i)
 302           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 303         enddo
 304       enddo
 305       do i=2,nres-1
 306         do j=1,3
 307           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 308           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 309         enddo
 310 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
 311 c     &   vbld_inv(i+nres)
 312       enddo
 313
 314 C Calculate internal coordinates.
 315       do ires=1,nres
 316         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
 317      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
 318      &    (c(j,nres+ires),j=1,3)
 319       enddo
 320 C      print *,"before int_from_cart"
 321       call int_from_cart(.true.,.false.)
 322       call sc_loc_geom(.true.)
 323 c      call chainbuild
 324 C Copy the coordinates to reference coordinates
 325 C Splits to single chain if occurs
 326
 327       if (lprn) then
 328       write (iout,'(/a)')
 329      &  "Cartesian coordinates of the reference structure"
 330       write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
 331      & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
 332       do ires=1,nres
 333         write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
 334      &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
 335      &    (c(j,ires+nres),j=1,3)
 336       enddo
 337       endif
 338       call int_from_cart(.true.,lprn)
 339       call sc_loc_geom(lprn)
 340       ishift_pdb=ishift
 341       return
 342       end
 343 c---------------------------------------------------------------------------
 344       subroutine int_from_cart(lside,lprn)
 345       implicit none
 346       include 'DIMENSIONS'
 347       include 'DIMENSIONS.ZSCOPT'
 348       include 'COMMON.LOCAL'
 349       include 'COMMON.VAR'
 350       include 'COMMON.CHAIN'
 351       include 'COMMON.INTERACT'
 352       include 'COMMON.IOUNITS'
 353       include 'COMMON.GEO'
 354       include 'COMMON.NAMES'
 355       integer i,j,iti
 356       double precision dist,alpha,beta,di
 357       character*3 seq,atom,res
 358       character*80 card
 359       double precision sccor(3,20)
 360       integer rescode
 361       logical lside,lprn
 362       if (lprn) then
 363         write (iout,'(/a)')
 364      &  'Internal coordinates of the reference structure.'
 365         if (lside) then
 366           write (iout,'(8a)') 'Residue','       dvb','     Theta',
 367      & '       Phi','    Dsc_id','       Dsc','     Alpha',
 368      & '     Omega'
 369         else
 370           write (iout,'(4a)') 'Residue','       dvb','     Theta',
 371      & '       Phi'
 372         endif
 373       endif
 374       do i=1,nres-1
 375         iti=itype(i)
 376         if (iti.ne.ntyp1 .and. itype(i+1).ne.ntyp1 .and.
 377      &    (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.6.5D0)) then
 378           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 379           stop
 380         endif
 381         vbld(i+1)=dist(i,i+1)
 382         vbld_inv(i+1)=1.0d0/vbld(i+1)
 383         if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
 384         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
 385         do j=1,3
 386           dc(j,i)=c(j,i+1)-c(j,i)
 387           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 388         enddo
 389       enddo
 390       if (lside) then
 391         do i=2,nct
 392           do j=1,3
 393             c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
 394           enddo
 395           iti=itype(i)
 396           di=dist(i,nres+i)
 397           vbld(i+nres)=di
 398           if (itype(i).ne.10) then
 399             vbld_inv(i+nres)=1.0d0/di
 400             do j=1,3
 401               dc(j,i+nres)=c(j,i+nres)-c(j,i)
 402               dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 403             enddo
 404           else
 405             vbld_inv(i+nres)=0.0d0
 406             do j=1,3
 407               dc(j,i+nres)=0.0d0
 408               dc_norm(j,i+nres)=0.0d0
 409             enddo
 410           endif
 411           if (iti.ne.10) then
 412             alph(i)=alpha(nres+i,i,maxres2)
 413             omeg(i)=beta(nres+i,i,maxres2,i+1)
 414           endif
 415           if (lprn)
 416      &    write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 417      &    rad2deg*theta(i),rad2deg*phi(i),dsc(iti),di,
 418      &    rad2deg*alph(i),rad2deg*omeg(i)
 419           if (iti.ne.10) then
 420           endif
 421         enddo
 422       else if (lprn) then
 423         do i=2,nres
 424           iti=itype(i)
 425           write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 426      &    rad2deg*theta(i),rad2deg*phi(i)
 427         enddo
 428       endif
 429       return
 430       end
 431 c-------------------------------------------------------------------------------
 432       subroutine sc_loc_geom(lprn)
 433       implicit real*8 (a-h,o-z)
 434       include 'DIMENSIONS'
 435       include 'DIMENSIONS.ZSCOPT'
 436       include 'COMMON.LOCAL'
 437       include 'COMMON.VAR'
 438       include 'COMMON.CHAIN'
 439       include 'COMMON.INTERACT'
 440       include 'COMMON.IOUNITS'
 441       include 'COMMON.GEO'
 442       include 'COMMON.NAMES'
 443       include 'COMMON.CONTROL'
 444       double precision x_prime(3),y_prime(3),z_prime(3)
 445       logical lprn
 446       do i=1,nres-1
 447         do j=1,3
 448           dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
 449         enddo
 450       enddo
 451       do i=2,nres-1
 452         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 453           do j=1,3
 454             dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
 455           enddo
 456         else
 457           do j=1,3
 458             dc_norm(j,i+nres)=0.0d0
 459           enddo
 460         endif
 461       enddo
 462       do i=2,nres-1
 463         costtab(i+1) =dcos(theta(i+1))
 464         sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
 465         cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
 466         sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
 467         cosfac2=0.5d0/(1.0d0+costtab(i+1))
 468         cosfac=dsqrt(cosfac2)
 469         sinfac2=0.5d0/(1.0d0-costtab(i+1))
 470         sinfac=dsqrt(sinfac2)
 471         it=itype(i)
 472         if (it.ne.10 .and. itype(i).ne.ntyp1) then
 473 c
 474 C  Compute the axes of tghe local cartesian coordinates system; store in
 475 c   x_prime, y_prime and z_prime
 476 c
 477         do j=1,3
 478           x_prime(j) = 0.00
 479           y_prime(j) = 0.00
 480           z_prime(j) = 0.00
 481         enddo
 482         do j = 1,3
 483           x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
 484           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
 485         enddo
 486         call vecpr(x_prime,y_prime,z_prime)
 487 c
 488 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
 489 C to local coordinate system. Store in xx, yy, zz.
 490 c
 491         xx=0.0d0
 492         yy=0.0d0
 493         zz=0.0d0
 494         do j = 1,3
 495           xx = xx + x_prime(j)*dc_norm(j,i+nres)
 496           yy = yy + y_prime(j)*dc_norm(j,i+nres)
 497           zz = zz + z_prime(j)*dc_norm(j,i+nres)
 498         enddo
 499
 500         xxtab(i)=xx
 501         yytab(i)=yy
 502         zztab(i)=zz
 503         else
 504         xxtab(i)=0.0d0
 505         yytab(i)=0.0d0
 506         zztab(i)=0.0d0
 507         endif
 508       enddo
 509       if (lprn) then
 510         do i=2,nres
 511           iti=itype(i)
 512           write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxtab(i),yytab(i),
 513      &     zztab(i)
 514         enddo
 515       endif
 516       return
 517       end
 518 c---------------------------------------------------------------------------
 519       subroutine sccenter(ires,nscat,sccor)
 520       implicit none
 521       include 'DIMENSIONS'
 522       include 'DIMENSIONS.ZSCOPT'
 523       include 'COMMON.CHAIN'
 524       integer ires,nscat,i,j
 525       double precision sccmj
 526       double precision sccor(3,20)
 527       do j=1,3
 528         sccmj=0.0D0
 529         do i=1,nscat
 530           sccmj=sccmj+sccor(j,i)
 531         enddo
 532         dc(j,ires)=sccmj/nscat
 533       enddo
 534       return
 535       end