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