source/maxlik/src_MD_T_maxlik-NEWCORR-PMF/readpdb.F.old

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