source/cluster/wham/src-old/readpdb.F.org

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