source/unres/src-HCD-5D/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       include 'COMMON.DISTFIT'
  15       include 'COMMON.SETUP'
  16       include 'COMMON.FRAG'
  17       integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity,
  18      &  ishift_pdb
  19       logical lprn /.false./,fail
  20       double precision e1(3),e2(3),e3(3)
  21       double precision dcj,efree_temp
  22       character*3 seq,res
  23       character*5 atom
  24       character*80 card
  25       double precision sccor(3,50)
  26       integer rescode
  27       efree_temp=0.0d0
  28       ibeg=1
  29       ishift1=0
  30       ishift=0
  31 c      write (2,*) "UNRES_PDB",unres_pdb
  32       ires=0
  33       ires_old=0
  34       iii=0
  35       lsecondary=.false.
  36       nhfrag=0
  37       nbfrag=0
  38       do
  39         read (ipdbin,'(a80)',end=10) card
  40 c        write (iout,'(a)') card
  41         if (card(:5).eq.'HELIX') then
  42          nhfrag=nhfrag+1
  43          lsecondary=.true.
  44          read(card(22:25),*) hfrag(1,nhfrag)
  45          read(card(34:37),*) hfrag(2,nhfrag)
  46         endif
  47         if (card(:5).eq.'SHEET') then
  48          nbfrag=nbfrag+1
  49          lsecondary=.true.
  50          read(card(24:26),*) bfrag(1,nbfrag)
  51          read(card(35:37),*) bfrag(2,nbfrag)
  52 crc----------------------------------------
  53 crc  to be corrected !!!
  54          bfrag(3,nbfrag)=bfrag(1,nbfrag)
  55          bfrag(4,nbfrag)=bfrag(2,nbfrag)
  56 crc----------------------------------------
  57         endif
  58         if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10
  59 c Read free energy
  60         if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
  61 C Fish out the ATOM cards.
  62         if (index(card(1:4),'ATOM').gt.0) then
  63           read (card(12:16),*) atom
  64 c          write (iout,*) "! ",atom," !",ires
  65 c          if (atom.eq.'CA' .or. atom.eq.'CH3') then
  66           read (card(23:26),*) ires
  67           read (card(18:20),'(a3)') res
  68 c          write (iout,*) "ires",ires,ires-ishift+ishift1,
  69 c     &      " ires_old",ires_old
  70 c          write (iout,*) "ishift",ishift," ishift1",ishift1
  71 c          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
  72           if (ires-ishift+ishift1.ne.ires_old) then
  73 C Calculate the CM of the preceding residue.
  74 c            if (ibeg.eq.0) call sccenter(ires,iii,sccor)
  75             if (ibeg.eq.0) then
  76 c              write (iout,*) "Calculating sidechain center iii",iii
  77               if (unres_pdb) then
  78                 do j=1,3
  79                   dc(j,ires)=sccor(j,iii)
  80                 enddo
  81               else
  82                 call sccenter(ires_old,iii,sccor)
  83               endif
  84               iii=0
  85             endif
  86 C Start new residue.
  87             if (res.eq.'Cl-' .or. res.eq.'Na+') then
  88               ires=ires_old
  89               cycle
  90             else if (ibeg.eq.1) then
  91 c              write (iout,*) "BEG ires",ires
  92               ishift=ires-1
  93               if (res.ne.'GLY' .and. res.ne. 'ACE') then
  94                 ishift=ishift-1
  95                 itype(1)=ntyp1
  96               endif
  97               ires=ires-ishift+ishift1
  98               ires_old=ires
  99 c              write (iout,*) "ishift",ishift," ires",ires,
 100 c     &         " ires_old",ires_old
 101               ibeg=0
 102             else
 103               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
 104               ires=ires-ishift+ishift1
 105               ires_old=ires
 106             endif
 107             if (res.eq.'ACE' .or. res.eq.'NHE') then
 108               itype(ires)=10
 109             else
 110               itype(ires)=rescode(ires,res,0)
 111             endif
 112           else
 113             ires=ires-ishift+ishift1
 114           endif
 115 c          write (iout,*) "ires_old",ires_old," ires",ires
 116           if (card(27:27).eq."A" .or. card(27:27).eq."B") then
 117 c            ishift1=ishift1+1
 118           endif
 119 c          write (2,*) "ires",ires," res ",res," ity",ity
 120           if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
 121      &       res.eq.'NHE'.and.atom(:2).eq.'HN') then
 122             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
 123 c            write (iout,*) "backbone ",atom
 124 #ifdef DEBUG
 125             write (iout,'(2i3,2x,a,3f8.3)')
 126      &      ires,itype(ires),res,(c(j,ires),j=1,3)
 127 #endif
 128             iii=iii+1
 129             do j=1,3
 130               sccor(j,iii)=c(j,ires)
 131             enddo
 132             if (ishift.ne.0) then
 133               ires_ca=ires+ishift-ishift1
 134             else
 135               ires_ca=ires
 136             endif
 137 c            write (*,*) card(23:27),ires,itype(ires)
 138           else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.
 139      &             atom.ne.'N' .and. atom.ne.'C' .and.
 140      &             atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.
 141      &             atom.ne.'OXT' .and. atom(:2).ne.'3H') then
 142 c            write (iout,*) "sidechain ",atom
 143             iii=iii+1
 144             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
 145           endif
 146         endif
 147       enddo
 148    10 continue
 149 #ifdef DEBUG
 150       write (iout,'(a,i5)') ' Number of residues found: ',ires
 151 #endif
 152       if (ires.eq.0) return
 153 C Calculate the CM of the last side chain.
 154       if (iii.gt.0)  then
 155       if (unres_pdb) then
 156         do j=1,3
 157           dc(j,ires)=sccor(j,iii)
 158         enddo
 159       else
 160         call sccenter(ires,iii,sccor)
 161       endif
 162       endif
 163       nres=ires
 164       nsup=nres
 165       nstart_sup=1
 166       if (itype(nres).ne.10) then
 167         nres=nres+1
 168         itype(nres)=ntyp1
 169         if (unres_pdb) then
 170 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
 171           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
 172           if (fail) then
 173             e2(1)=0.0d0
 174             e2(2)=1.0d0
 175             e2(3)=0.0d0
 176           endif
 177           do j=1,3
 178             c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
 179           enddo
 180         else
 181         do j=1,3
 182           dcj=c(j,nres-2)-c(j,nres-3)
 183           c(j,nres)=c(j,nres-1)+dcj
 184           c(j,2*nres)=c(j,nres)
 185         enddo
 186         endif
 187       endif
 188       do i=2,nres-1
 189         do j=1,3
 190           c(j,i+nres)=dc(j,i)
 191         enddo
 192       enddo
 193       do j=1,3
 194         c(j,nres+1)=c(j,1)
 195         c(j,2*nres)=c(j,nres)
 196       enddo
 197       if (itype(1).eq.ntyp1) then
 198         nsup=nsup-1
 199         nstart_sup=2
 200         if (unres_pdb) then
 201 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 202           call refsys(2,3,4,e1,e2,e3,fail)
 203           if (fail) then
 204             e2(1)=0.0d0
 205             e2(2)=1.0d0
 206             e2(3)=0.0d0
 207           endif
 208           do j=1,3
 209             c(j,1)=c(j,2)-3.8d0*e2(j)
 210           enddo
 211         else
 212         do j=1,3
 213           dcj=c(j,4)-c(j,3)
 214           c(j,1)=c(j,2)-dcj
 215           c(j,nres+1)=c(j,1)
 216         enddo
 217         endif
 218       endif
 219 C Copy the coordinates to reference coordinates
 220 c      do i=1,2*nres
 221 c        do j=1,3
 222 c          cref(j,i)=c(j,i)
 223 c        enddo
 224 c      enddo
 225 C Calculate internal coordinates.
 226       if (lprn) then
 227       write (iout,'(/a)')
 228      &  "Cartesian coordinates of the reference structure"
 229       write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
 230      & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
 231       do ires=1,nres
 232         write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
 233      &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
 234      &    (c(j,ires+nres),j=1,3)
 235       enddo
 236       endif
 237 C Calculate internal coordinates.
 238       if(me.eq.king.or..not.out1file)then
 239        write (iout,'(a)')
 240      &   "Backbone and SC coordinates as read from the PDB"
 241        do ires=1,nres
 242         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
 243      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
 244      &    (c(j,nres+ires),j=1,3)
 245        enddo
 246       endif
 247       call int_from_cart(.true.,.false.)
 248       call sc_loc_geom(.false.)
 249       do i=1,nres
 250         thetaref(i)=theta(i)
 251         phiref(i)=phi(i)
 252       enddo
 253       do i=1,nres-1
 254         do j=1,3
 255           dc(j,i)=c(j,i+1)-c(j,i)
 256           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 257         enddo
 258       enddo
 259       do i=2,nres-1
 260         do j=1,3
 261           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 262           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 263         enddo
 264 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
 265 c     &   vbld_inv(i+nres)
 266       enddo
 267 c      call chainbuild
 268 C Copy the coordinates to reference coordinates
 269       do i=1,2*nres
 270         do j=1,3
 271           cref(j,i)=c(j,i)
 272         enddo
 273       enddo
 274
 275
 276       do j=1,nbfrag
 277         do i=1,4
 278          bfrag(i,j)=bfrag(i,j)-ishift
 279         enddo
 280       enddo
 281
 282       do j=1,nhfrag
 283         do i=1,2
 284          hfrag(i,j)=hfrag(i,j)-ishift
 285         enddo
 286       enddo
 287       ishift_pdb=ishift
 288       return
 289       end
 290 c-----------------------------------------------------------------------
 291       subroutine readpdb_template(k)
 292 C Read the PDB file with gaps for read_constr_homology with read2sigma
 293 C and convert the peptide geometry into virtual-chain geometry.
 294       implicit real*8 (a-h,o-z)
 295       include 'DIMENSIONS'
 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       include 'COMMON.CONTROL'
 304       include 'COMMON.DISTFIT'
 305       include 'COMMON.SETUP'
 306       include 'COMMON.MD'
 307       integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity,
 308      &  ishift_pdb
 309       logical lprn /.false./,fail
 310       double precision e1(3),e2(3),e3(3)
 311       double precision dcj,efree_temp
 312       character*3 seq,res
 313       character*5 atom
 314       character*80 card
 315       double precision sccor(3,50)
 316       integer rescode
 317       efree_temp=0.0d0
 318       ibeg=1
 319       ishift1=0
 320       ishift=0
 321 c      write (2,*) "UNRES_PDB",unres_pdb
 322       ires=0
 323       ires_old=0
 324       iii=0
 325       lsecondary=.false.
 326       nhfrag=0
 327       nbfrag=0
 328       do
 329         read (ipdbin,'(a80)',end=10) card
 330 c        write (iout,'(a)') card
 331         if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10
 332 C Fish out the ATOM cards.
 333         if (index(card(1:4),'ATOM').gt.0) then
 334           read (card(12:16),*) atom
 335 c          write (iout,*) "! ",atom," !",ires
 336 c          if (atom.eq.'CA' .or. atom.eq.'CH3') then
 337           read (card(23:26),*) ires
 338           read (card(18:20),'(a3)') res
 339 c          write (iout,*) "ires",ires,ires-ishift+ishift1,
 340 c     &      " ires_old",ires_old
 341 c          write (iout,*) "ishift",ishift," ishift1",ishift1
 342 c          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
 343           if (ires-ishift+ishift1.ne.ires_old) then
 344 C Calculate the CM of the preceding residue.
 345 c            if (ibeg.eq.0) call sccenter(ires,iii,sccor)
 346             if (ibeg.eq.0) then
 347 c              write (iout,*) "Calculating sidechain center iii",iii
 348               if (unres_pdb) then
 349                 do j=1,3
 350                   dc(j,ires)=sccor(j,iii)
 351                 enddo
 352               else
 353                 call sccenter(ires_old,iii,sccor)
 354               endif
 355               iii=0
 356             endif
 357 C Start new residue.
 358             if (res.eq.'Cl-' .or. res.eq.'Na+') then
 359               ires=ires_old
 360               cycle
 361             else if (ibeg.eq.1) then
 362 c              write (iout,*) "BEG ires",ires
 363               ishift=ires-1
 364               if (res.ne.'GLY' .and. res.ne. 'ACE') then
 365                 ishift=ishift-1
 366                 itype(1)=ntyp1
 367               endif
 368               ires=ires-ishift+ishift1
 369               ires_old=ires
 370 c              write (iout,*) "ishift",ishift," ires",ires,
 371 c     &         " ires_old",ires_old
 372               ibeg=0
 373             else
 374               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
 375               ires=ires-ishift+ishift1
 376               ires_old=ires
 377             endif
 378             if (res.eq.'ACE' .or. res.eq.'NHE') then
 379               itype(ires)=10
 380             else
 381               itype(ires)=rescode(ires,res,0)
 382             endif
 383           else
 384             ires=ires-ishift+ishift1
 385           endif
 386 c          write (iout,*) "ires_old",ires_old," ires",ires
 387           if (card(27:27).eq."A" .or. card(27:27).eq."B") then
 388 c            ishift1=ishift1+1
 389           endif
 390 c          write (2,*) "ires",ires," res ",res," ity",ity
 391           if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
 392      &       res.eq.'NHE'.and.atom(:2).eq.'HN') then
 393             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
 394 c            write (iout,*) "backbone ",atom
 395 #ifdef DEBUG
 396             write (iout,'(2i3,2x,a,3f8.3)')
 397      &      ires,itype(ires),res,(c(j,ires),j=1,3)
 398 #endif
 399             iii=iii+1
 400             do j=1,3
 401               sccor(j,iii)=c(j,ires)
 402             enddo
 403             if (ishift.ne.0) then
 404               ires_ca=ires+ishift-ishift1
 405             else
 406               ires_ca=ires
 407             endif
 408 c            write (*,*) card(23:27),ires,itype(ires)
 409           else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.
 410      &             atom.ne.'N' .and. atom.ne.'C' .and.
 411      &             atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.
 412      &             atom.ne.'OXT' .and. atom(:2).ne.'3H') then
 413 c            write (iout,*) "sidechain ",atom
 414             iii=iii+1
 415             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
 416           endif
 417         endif
 418       enddo
 419    10 continue
 420 #ifdef DEBUG
 421       write (iout,'(a,i5)') ' Number of residues found: ',ires
 422 #endif
 423       if (ires.eq.0) return
 424 C Calculate the CM of the last side chain.
 425       if (iii.gt.0)  then
 426       if (unres_pdb) then
 427         do j=1,3
 428           dc(j,ires)=sccor(j,iii)
 429         enddo
 430       else
 431         call sccenter(ires,iii,sccor)
 432       endif
 433       endif
 434       nres=ires
 435       nsup=nres
 436       nstart_sup=1
 437       if (itype(nres).ne.10) then
 438         nres=nres+1
 439         itype(nres)=ntyp1
 440         do j=1,3
 441           dcj=c(j,nres-2)-c(j,nres-3)
 442           c(j,nres)=c(j,nres-1)+dcj
 443           c(j,2*nres)=c(j,nres)
 444         enddo
 445       endif
 446       do i=2,nres-1
 447         do j=1,3
 448           c(j,i+nres)=dc(j,i)
 449         enddo
 450       enddo
 451       do j=1,3
 452         c(j,nres+1)=c(j,1)
 453         c(j,2*nres)=c(j,nres)
 454       enddo
 455       if (itype(1).eq.ntyp1) then
 456         nsup=nsup-1
 457         nstart_sup=2
 458         if (unres_pdb) then
 459 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 460           call refsys(2,3,4,e1,e2,e3,fail)
 461           if (fail) then
 462             e2(1)=0.0d0
 463             e2(2)=1.0d0
 464             e2(3)=0.0d0
 465           endif
 466           do j=1,3
 467             c(j,1)=c(j,2)-3.8d0*e2(j)
 468           enddo
 469         else
 470         do j=1,3
 471           dcj=c(j,4)-c(j,3)
 472           c(j,1)=c(j,2)-dcj
 473           c(j,nres+1)=c(j,1)
 474         enddo
 475         endif
 476       endif
 477 C Calculate internal coordinates.
 478       if (lprn) then
 479       write (iout,'(/a)')
 480      &  "Cartesian coordinates of the reference structure"
 481       write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
 482      & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
 483       do ires=1,nres
 484         write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
 485      &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
 486      &    (c(j,ires+nres),j=1,3)
 487       enddo
 488       endif
 489 C Calculate internal coordinates.
 490       if(me.eq.king.or..not.out1file)then
 491        write (iout,'(a)')
 492      &   "Backbone and SC coordinates as read from the PDB"
 493        do ires=1,nres
 494         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
 495      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
 496      &    (c(j,nres+ires),j=1,3)
 497        enddo
 498       endif
 499       call int_from_cart(.true.,.false.)
 500       call sc_loc_geom(.false.)
 501       do i=1,nres
 502         thetaref(i)=theta(i)
 503         phiref(i)=phi(i)
 504       enddo
 505       do i=1,nres-1
 506         do j=1,3
 507           dc(j,i)=c(j,i+1)-c(j,i)
 508           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 509         enddo
 510       enddo
 511       do i=2,nres-1
 512         do j=1,3
 513           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 514           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 515         enddo
 516 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
 517 c     &   vbld_inv(i+nres)
 518       enddo
 519 c      call chainbuild
 520 C Copy the coordinates to reference coordinates
 521       do i=1,2*nres
 522         do j=1,3
 523           cref(j,i)=c(j,i)
 524           chomo(j,i,k)=c(j,i)
 525         enddo
 526       enddo
 527
 528
 529       ishift_pdb=ishift
 530       return
 531       end