source/maxlik/src_MD_T_maxlik-NEWCORR-PMF-PDB/readpdb.F.safe

   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,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       write (2,*) "UNRES_PDB",unres_pdb
  29       ires=0
  30       do i=1,10000
  31         read (ipdbin,'(a80)',end=10) card
  32         if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10
  33 c Read free energy
  34         if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
  35 C Fish out the ATOM cards.
  36         if (index(card(1:4),'ATOM').gt.0) then
  37           read (card(14:16),'(a3)') atom
  38           if (atom.eq.'CA' .or. atom.eq.'CH3') then
  39 C Calculate the CM of the preceding residue.
  40 c            if (ibeg.eq.0) call sccenter(ires,iii,sccor)
  41             if (ibeg.eq.0) then
  42               if (unres_pdb) then
  43                 do j=1,3
  44                   dc(j,ires)=sccor(j,iii)
  45                 enddo
  46               else
  47                 call sccenter(ires,iii,sccor)
  48               endif
  49             endif
  50 C Start new residue.
  51             ires_old=ires+ishift-ishift1
  52             read (card(23:26),*) ires
  53 c            print *,"ires_old",ires_old," ires",ires
  54             if (card(27:27).eq."A" .or. card(27:27).eq."B") then
  55 c              ishift1=ishift1+1
  56             endif
  57             read (card(18:20),'(a3)') res
  58             if (ibeg.eq.1) then
  59               ishift=ires-1
  60               if (res.ne.'GLY' .and. res.ne. 'ACE') then
  61                 ishift=ishift-1
  62                 itype(1)=21
  63               endif
  64               ibeg=0
  65             else
  66               ishift=ishift+ires-ires_old-1
  67             endif
  68             ires=ires-ishift+ishift1
  69             if (res.eq.'ACE' .or. res.eq.'NHE') then
  70               itype(ires)=10
  71             else
  72               itype(ires)=rescode(ires,res,0)
  73             endif
  74 c            write (2,*) "ires",ires," res ",res," ity",ity
  75             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
  76 #ifdef DEBUG
  77             write (iout,'(2i3,2x,a,3f8.3)')
  78      &      ires,itype(ires),res,(c(j,ires),j=1,3)
  79 #endif
  80             iii=1
  81             do j=1,3
  82               sccor(j,iii)=c(j,ires)
  83             enddo
  84 c            write (*,*) card(23:27),ires,itype(ires)
  85           else if (atom.ne.'O  '.and.atom(1:1).ne.'H' .and.
  86      &             atom.ne.'N  ' .and. atom.ne.'C   ') then
  87             iii=iii+1
  88             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
  89           endif
  90         endif
  91       enddo
  92    10 write (iout,'(a,i5)') ' Number of residues found: ',ires
  93       if (ires.eq.0) return1
  94 C Calculate the CM of the last side chain.
  95       if (unres_pdb) then
  96         do j=1,3
  97           dc(j,ires)=sccor(j,iii)
  98         enddo
  99       else
 100         call sccenter(ires,iii,sccor)
 101       endif
 102       nres=ires
 103       nsup(iprot)=nres
 104       nstart_sup(iprot)=1
 105       if (itype(nres).ne.10) then
 106         nres=nres+1
 107         itype(nres)=21
 108         do j=1,3
 109           dcj=c(j,nres-2)-c(j,nres-3)
 110           c(j,nres)=c(j,nres-1)+dcj
 111           c(j,2*nres)=c(j,nres)
 112         enddo
 113       endif
 114       do i=2,nres-1
 115         do j=1,3
 116           c(j,i+nres)=dc(j,i)
 117         enddo
 118       enddo
 119       do j=1,3
 120         c(j,nres+1)=c(j,1)
 121         c(j,2*nres)=c(j,nres)
 122       enddo
 123       if (itype(1).eq.21) then
 124         nsup(iprot)=nsup(iprot)-1
 125         nstart_sup(iprot)=2
 126         if (unres_pdb) then
 127 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 128           call refsys(2,3,4,e1,e2,e3,fail)
 129           if (fail) then
 130             e2(1)=0.0d0
 131             e2(2)=1.0d0
 132             e2(3)=0.0d0
 133           endif
 134           do j=1,3
 135             c(j,1)=c(j,2)-3.8d0*e2(j)
 136           enddo
 137         else
 138         do j=1,3
 139           dcj=c(j,4)-c(j,3)
 140           c(j,1)=c(j,2)-dcj
 141           c(j,nres+1)=c(j,1)
 142         enddo
 143         endif
 144       endif
 145 C Copy the coordinates to reference coordinates
 146 c      do i=1,2*nres
 147 c        do j=1,3
 148 c          cref(j,i,iprot)=c(j,i)
 149 c        enddo
 150 c      enddo
 151 C Calculate internal coordinates.
 152       if (lprn) then
 153       write (iout,'(/a)')
 154      &  "Cartesian coordinates of the reference structure"
 155       write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
 156      & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
 157       do ires=1,nres
 158         write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
 159      &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
 160      &    (c(j,ires+nres),j=1,3)
 161       enddo
 162       endif
 163       call int_from_cart(.true.,lprn)
 164 c      do i=1,nres
 165 c        phi_ref(i,iprot)=phi(i)
 166 c        theta_ref(i,iprot)=theta(i)
 167 c        alph_ref(i,iprot)=alph(i)
 168 c        omeg_ref(i,iprot)=omeg(i)
 169 c      enddo
 170       ishift_pdb=ishift
 171       return
 172       end
 173 c---------------------------------------------------------------------------
 174       subroutine int_from_cart(lside,lprn)
 175       implicit none
 176       include 'DIMENSIONS'
 177       include 'DIMENSIONS.ZSCOPT'
 178       include 'COMMON.LOCAL'
 179       include 'COMMON.VAR'
 180       include 'COMMON.CHAIN'
 181       include 'COMMON.INTERACT'
 182       include 'COMMON.IOUNITS'
 183       include 'COMMON.GEO'
 184       include 'COMMON.NAMES'
 185       integer i,j,iti
 186       double precision dist,alpha,beta,di
 187       character*3 seq,atom,res
 188       character*80 card
 189       double precision sccor(3,20)
 190       integer rescode
 191       logical lside,lprn
 192       if (lprn) then
 193         write (iout,'(/a)')
 194      &  'Internal coordinates of the reference structure.'
 195         if (lside) then
 196           write (iout,'(8a)') 'Residue','       dvb','     Theta',
 197      & '       Phi','    Dsc_id','       Dsc','     Alpha',
 198      & '     Omega'
 199         else
 200           write (iout,'(4a)') 'Residue','       dvb','     Theta',
 201      & '       Phi'
 202         endif
 203       endif
 204       do i=1,nres-1
 205         iti=itype(i)
 206         if (iti.ne.21 .and. itype(i+1).ne.21 .and.
 207      &    (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0)) then
 208           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 209           stop
 210         endif
 211         vbld(i+1)=dist(i,i+1)
 212         vbld_inv(i+1)=1.0d0/vbld(i+1)
 213         if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
 214         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
 215         do j=1,3
 216           dc(j,i)=c(j,i+1)-c(j,i)
 217           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 218         enddo
 219       enddo
 220       if (lside) then
 221         do i=2,nct
 222           do j=1,3
 223             c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
 224           enddo
 225           iti=itype(i)
 226           di=dist(i,nres+i)
 227           vbld(i+nres)=di
 228           if (itype(i).ne.10) then
 229             vbld_inv(i+nres)=1.0d0/di
 230             do j=1,3
 231               dc(j,i+nres)=c(j,i+nres)-c(j,i)
 232               dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 233             enddo
 234           else
 235             vbld_inv(i+nres)=0.0d0
 236             do j=1,3
 237               dc(j,i+nres)=0.0d0
 238               dc_norm(j,i+nres)=0.0d0
 239             enddo
 240           endif
 241           if (iti.ne.10) then
 242             alph(i)=alpha(nres+i,i,maxres2)
 243             omeg(i)=beta(nres+i,i,maxres2,i+1)
 244           endif
 245           if (lprn)
 246      &    write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 247      &    rad2deg*theta(i),rad2deg*phi(i),dsc(iti),di,
 248      &    rad2deg*alph(i),rad2deg*omeg(i)
 249           if (iti.ne.10) then
 250           endif
 251         enddo
 252       else if (lprn) then
 253         do i=2,nres
 254           iti=itype(i)
 255           write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 256      &    rad2deg*theta(i),rad2deg*phi(i)
 257         enddo
 258       endif
 259       return
 260       end
 261 c---------------------------------------------------------------------------
 262       subroutine sccenter(ires,nscat,sccor)
 263       implicit none
 264       include 'DIMENSIONS'
 265       include 'DIMENSIONS.ZSCOPT'
 266       include 'COMMON.CHAIN'
 267       integer ires,nscat,i,j
 268       double precision sccmj
 269       double precision sccor(3,20)
 270       do j=1,3
 271         sccmj=0.0D0
 272         do i=1,nscat
 273           sccmj=sccmj+sccor(j,i)
 274         enddo
 275         dc(j,ires)=sccmj/nscat
 276       enddo
 277       return
 278       end