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