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

   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
  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 i=1,10000
  39          write (iout,'(a)') card
  40 c        call flush(iout)
  41 c        if (card(:3).eq.'END') goto 10
  42         if (card(:3).eq.'END' .or.
  43      &    card(:3).eq.'TER'.and.prevcard(:5).eq.'REMARK') goto 10
  44         if (card(:3).eq.'TER') then
  45 C End current chain
  46           write (iout,*) "Im here!"
  47           write (iout,*) "ires_old",ires_old
  48           ires_old=ires+2
  49           itype(ires_old-1)=ntyp1
  50           iterter(ires_old-1)=1
  51           itype(ires_old)=ntyp1
  52           iterter(ires_old)=1
  53           ibeg=2
  54           write (iout,*) "ires_old",ires_old
  55           write (iout,*) "Chain ended",ires,ishift,ires_old
  56           if (unres_pdb) then
  57             do j=1,3
  58               dc(j,ires)=sccor(j,iii)
  59             enddo
  60           else
  61             call sccenter(ires,iii,sccor)
  62           endif
  63         endif
  64         prevcard=card
  65 c Read free energy
  66         if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
  67 C Fish out the ATOM cards.
  68         if (index(card(1:4),'ATOM').gt.0) then
  69           read (card(14:16),'(a3)') atom
  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 c          call flush(iout)
  78           iires=ires-ishift+ishift1
  79           if (iires.ne.ires_old) then
  80 C Calculate the CM of the preceding residue.
  81             if (ibeg.eq.0) then
  82 c              write (iout,*) "Calculating sidechain center iii",iii
  83 c              call flush(iout)
  84               if (unres_pdb) then
  85                 do j=1,3
  86                   dc(j,ires)=sccor(j,iii)
  87                 enddo
  88               else
  89                 call sccenter(iires-1,iii,sccor)
  90               endif
  91 c              write (iout,*) "sidechain center position calculated"
  92 c              call flush(iout)
  93               iii=0
  94             endif
  95 C Start new residue.
  96             if (res.eq.'Cl-' .or. res.eq.'Na+') then
  97               ires=ires_old
  98               read (ipdbin,'(a80)',end=10) card
  99               cycle
 100             else if (ibeg.eq.1) then
 101 c              write (iout,*) "BEG ires",ires," iires",iires
 102 c              call flush(iout)
 103               ishift=ires-1
 104               if (res.ne.'GLY' .and. res.ne. 'ACE') then
 105                 ishift=ishift-1
 106                 itype(1)=ntyp1
 107               endif
 108               ires=ires-ishift+ishift1
 109 c              ires_old=ires
 110 c              write (iout,*) "ishift",ishift," ires",ires,
 111 c     &         " iires",iires," ires_old",ires_old
 112 c              call flush(iout)
 113               ibeg=0
 114             else
 115               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
 116               ires=ires-ishift+ishift1
 117 c              ires_old=ires
 118 c              write (iout,*) "RESSHIFT",ires,iires,ishift,ishift1
 119 c              call flush(iout)
 120             endif
 121 c            write (iout,*) "ires",ires," iires",iires
 122 c            call flush(iout)
 123             if (res.eq.'ACE' .or. res.eq.'NHE') then
 124               itype(ires)=10
 125             else
 126               itype(ires)=rescode(ires,res,0)
 127             endif
 128           else
 129             ires=ires-ishift+ishift1
 130           endif
 131 c          write (iout,*) "ires_old",ires_old," ires",ires," iires",iires
 132 c          call flush(iout)
 133           if (card(27:27).eq."A" .or. card(27:27).eq."B") then
 134 c            ishift1=ishift1+1
 135           endif
 136 c          write (2,*) "iires",iires," res ",res," ity",ity
 137           if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
 138      &       res.eq.'NHE'.and.atom(:2).eq.'HN') then
 139             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
 140 #ifdef DEBUG
 141             write (iout,'(2i3,2x,a,3f8.3)')
 142      &      iires,itype(iires),res,(c(j,ires),j=1,3)
 143 #endif
 144             iii=iii+1
 145             do j=1,3
 146               sccor(j,iii)=c(j,ires)
 147             enddo
 148 c            write (*,*) card(23:27),ires,itype(ires)
 149           else if (atom.ne.'O  '.and.atom(1:1).ne.'H' .and.
 150      &             atom.ne.'N  ' .and. atom.ne.'C   ' .and.
 151      &             atom.ne.'OXT') then
 152             iii=iii+1
 153             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
 154           endif
 155         endif
 156         read (ipdbin,'(a80)',end=10) card
 157       enddo
 158    10 write (iout,'(a,i5)') ' Number of residues found: ',ires
 159       prevcard=card
 160       if (ires.eq.0) return1
 161 C Calculate the CM of the last side chain.
 162       if (iii.gt.0)  then
 163       if (unres_pdb) then
 164         do j=1,3
 165           dc(j,ires)=sccor(j,iii)
 166         enddo
 167       else
 168         call sccenter(iires,iii,sccor)
 169       endif
 170       endif
 171       nres=ires
 172       nsup(iprot)=nres
 173       nstart_sup(iprot)=1
 174       if (itype(nres).ne.10) then
 175         nres=nres+1
 176         itype(nres)=ntyp1
 177         do j=1,3
 178           dcj=c(j,nres-2)-c(j,nres-3)
 179           c(j,nres)=c(j,nres-1)+dcj
 180           c(j,2*nres)=c(j,nres)
 181         enddo
 182       endif
 183       do i=2,nres-1
 184         do j=1,3
 185           c(j,i+nres)=dc(j,i)
 186         enddo
 187       enddo
 188       do j=1,3
 189         c(j,nres+1)=c(j,1)
 190         c(j,2*nres)=c(j,nres)
 191       enddo
 192       if (itype(1).eq.ntyp1) then
 193         nsup(iprot)=nsup(iprot)-1
 194         nstart_sup(iprot)=2
 195         if (unres_pdb) then
 196 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 197           call refsys(2,3,4,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,1)=c(j,2)-3.8d0*e2(j)
 205           enddo
 206         else
 207         do j=1,3
 208           dcj=c(j,4)-c(j,3)
 209           c(j,1)=c(j,2)-dcj
 210           c(j,nres+1)=c(j,1)
 211         enddo
 212         endif
 213       endif
 214 C Copy the coordinates to reference coordinates
 215 c      do i=1,2*nres
 216 c        do j=1,3
 217 c          cref(j,i,iprot)=c(j,i)
 218 c        enddo
 219 c      enddo
 220 C Calculate internal coordinates.
 221       if (lprn) then
 222       write (iout,'(/a)')
 223      &  "Cartesian coordinates of the reference structure"
 224       write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
 225      & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
 226       do ires=1,nres
 227         write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
 228      &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
 229      &    (c(j,ires+nres),j=1,3)
 230       enddo
 231       endif
 232       call int_from_cart(.true.,lprn)
 233       call sc_loc_geom(lprn)
 234 c      do i=1,nres
 235 c        phi_ref(i,iprot)=phi(i)
 236 c        theta_ref(i,iprot)=theta(i)
 237 c        alph_ref(i,iprot)=alph(i)
 238 c        omeg_ref(i,iprot)=omeg(i)
 239 c      enddo
 240       ishift_pdb=ishift
 241       stop
 242       return
 243       end
 244 c---------------------------------------------------------------------------
 245       subroutine int_from_cart(lside,lprn)
 246       implicit none
 247       include 'DIMENSIONS'
 248       include 'DIMENSIONS.ZSCOPT'
 249       include 'COMMON.LOCAL'
 250       include 'COMMON.VAR'
 251       include 'COMMON.CHAIN'
 252       include 'COMMON.INTERACT'
 253       include 'COMMON.IOUNITS'
 254       include 'COMMON.GEO'
 255       include 'COMMON.NAMES'
 256       integer i,j,iti
 257       double precision dist,alpha,beta,di
 258       character*3 seq,atom,res
 259       character*80 card
 260       double precision sccor(3,20)
 261       integer rescode
 262       logical lside,lprn
 263       if (lprn) then
 264         write (iout,'(/a)')
 265      &  'Internal coordinates of the reference structure.'
 266         if (lside) then
 267           write (iout,'(8a)') 'Residue','       dvb','     Theta',
 268      & '       Phi','    Dsc_id','       Dsc','     Alpha',
 269      & '     Omega'
 270         else
 271           write (iout,'(4a)') 'Residue','       dvb','     Theta',
 272      & '       Phi'
 273         endif
 274       endif
 275       do i=1,nres-1
 276         iti=itype(i)
 277         if (iti.ne.ntyp1 .and. itype(i+1).ne.ntyp1 .and.
 278      &    (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.6.5D0)) then
 279           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 280           stop
 281         endif
 282         vbld(i+1)=dist(i,i+1)
 283         vbld_inv(i+1)=1.0d0/vbld(i+1)
 284         if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
 285         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
 286         do j=1,3
 287           dc(j,i)=c(j,i+1)-c(j,i)
 288           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 289         enddo
 290       enddo
 291       if (lside) then
 292         do i=2,nct
 293           do j=1,3
 294             c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
 295           enddo
 296           iti=itype(i)
 297           di=dist(i,nres+i)
 298           vbld(i+nres)=di
 299           if (itype(i).ne.10) then
 300             vbld_inv(i+nres)=1.0d0/di
 301             do j=1,3
 302               dc(j,i+nres)=c(j,i+nres)-c(j,i)
 303               dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 304             enddo
 305           else
 306             vbld_inv(i+nres)=0.0d0
 307             do j=1,3
 308               dc(j,i+nres)=0.0d0
 309               dc_norm(j,i+nres)=0.0d0
 310             enddo
 311           endif
 312           if (iti.ne.10) then
 313             alph(i)=alpha(nres+i,i,maxres2)
 314             omeg(i)=beta(nres+i,i,maxres2,i+1)
 315           endif
 316           if (lprn)
 317      &    write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 318      &    rad2deg*theta(i),rad2deg*phi(i),dsc(iti),di,
 319      &    rad2deg*alph(i),rad2deg*omeg(i)
 320           if (iti.ne.10) then
 321           endif
 322         enddo
 323       else if (lprn) then
 324         do i=2,nres
 325           iti=itype(i)
 326           write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 327      &    rad2deg*theta(i),rad2deg*phi(i)
 328         enddo
 329       endif
 330       return
 331       end
 332 c-------------------------------------------------------------------------------
 333       subroutine sc_loc_geom(lprn)
 334       implicit real*8 (a-h,o-z)
 335       include 'DIMENSIONS'
 336       include 'DIMENSIONS.ZSCOPT'
 337       include 'COMMON.LOCAL'
 338       include 'COMMON.VAR'
 339       include 'COMMON.CHAIN'
 340       include 'COMMON.INTERACT'
 341       include 'COMMON.IOUNITS'
 342       include 'COMMON.GEO'
 343       include 'COMMON.NAMES'
 344       include 'COMMON.CONTROL'
 345       double precision x_prime(3),y_prime(3),z_prime(3)
 346       logical lprn
 347       do i=1,nres-1
 348         do j=1,3
 349           dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
 350         enddo
 351       enddo
 352       do i=2,nres-1
 353         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 354           do j=1,3
 355             dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
 356           enddo
 357         else
 358           do j=1,3
 359             dc_norm(j,i+nres)=0.0d0
 360           enddo
 361         endif
 362       enddo
 363       do i=2,nres-1
 364         costtab(i+1) =dcos(theta(i+1))
 365         sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
 366         cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
 367         sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
 368         cosfac2=0.5d0/(1.0d0+costtab(i+1))
 369         cosfac=dsqrt(cosfac2)
 370         sinfac2=0.5d0/(1.0d0-costtab(i+1))
 371         sinfac=dsqrt(sinfac2)
 372         it=itype(i)
 373         if (it.ne.10 .and. itype(i).ne.ntyp1) then
 374 c
 375 C  Compute the axes of tghe local cartesian coordinates system; store in
 376 c   x_prime, y_prime and z_prime
 377 c
 378         do j=1,3
 379           x_prime(j) = 0.00
 380           y_prime(j) = 0.00
 381           z_prime(j) = 0.00
 382         enddo
 383         do j = 1,3
 384           x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
 385           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
 386         enddo
 387         call vecpr(x_prime,y_prime,z_prime)
 388 c
 389 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
 390 C to local coordinate system. Store in xx, yy, zz.
 391 c
 392         xx=0.0d0
 393         yy=0.0d0
 394         zz=0.0d0
 395         do j = 1,3
 396           xx = xx + x_prime(j)*dc_norm(j,i+nres)
 397           yy = yy + y_prime(j)*dc_norm(j,i+nres)
 398           zz = zz + z_prime(j)*dc_norm(j,i+nres)
 399         enddo
 400
 401         xxtab(i)=xx
 402         yytab(i)=yy
 403         zztab(i)=zz
 404         else
 405         xxtab(i)=0.0d0
 406         yytab(i)=0.0d0
 407         zztab(i)=0.0d0
 408         endif
 409       enddo
 410       if (lprn) then
 411         do i=2,nres
 412           iti=itype(i)
 413           write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxtab(i),yytab(i),
 414      &     zztab(i)
 415         enddo
 416       endif
 417       return
 418       end
 419 c---------------------------------------------------------------------------
 420       subroutine sccenter(ires,nscat,sccor)
 421       implicit none
 422       include 'DIMENSIONS'
 423       include 'DIMENSIONS.ZSCOPT'
 424       include 'COMMON.CHAIN'
 425       integer ires,nscat,i,j
 426       double precision sccmj
 427       double precision sccor(3,20)
 428       do j=1,3
 429         sccmj=0.0D0
 430         do i=1,nscat
 431           sccmj=sccmj+sccor(j,i)
 432         enddo
 433         dc(j,ires)=sccmj/nscat
 434       enddo
 435       return
 436       end