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