source/cluster/wham/src/readpdb.F

   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       logical lsecondary
  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 c
 247       enddo
 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(lside1,lprn)
 289       implicit real*8 (a-h,o-z)
 290       include 'DIMENSIONS'
 291 #ifdef MPI
 292       include "mpif.h"
 293 #endif
 294       include 'COMMON.LOCAL'
 295       include 'COMMON.VAR'
 296       include 'COMMON.CHAIN'
 297       include 'COMMON.INTERACT'
 298       include 'COMMON.IOUNITS'
 299       include 'COMMON.GEO'
 300       include 'COMMON.NAMES'
 301       include 'COMMON.CONTROL'
 302       include 'COMMON.SETUP'
 303       character*3 seq,atom,res
 304 c      character*5 atom
 305       character*80 card
 306       double precision sccor(3,20)
 307 c     dimension sccor(3,20)
 308       integer rescode
 309       logical lside1,lprn
 310       double precision dist,alpha,beta,di
 311       if(me.eq.king.or..not.out1file)then
 312        if (lprn) then
 313         write (iout,'(/a)')
 314      &  'Internal coordinates calculated from crystal structure.'
 315         if (lside) then
 316           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
 317      & '     Gamma','    Dsc_id','       Dsc','     Alpha',
 318      & '     Beta '
 319         else
 320           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
 321      & '     Gamma'
 322         endif
 323        endif
 324       endif
 325       do i=1,nres-1
 326         iti=itype(i)
 327         if (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0) then
 328           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 329 ctest          stop
 330         endif
 331         vbld(i+1)=dist(i,i+1)
 332         vbld_inv(i+1)=1.0d0/vbld(i+1)
 333         if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
 334         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
 335       enddo
 336 c      if (unres_pdb) then
 337 c        if (itype(1).eq.21) then
 338 c          theta(3)=90.0d0*deg2rad
 339 c          phi(4)=180.0d0*deg2rad
 340 c          vbld(2)=3.8d0
 341 c          vbld_inv(2)=1.0d0/vbld(2)
 342 c        endif
 343 c        if (itype(nres).eq.21) then
 344 c          theta(nres)=90.0d0*deg2rad
 345 c          phi(nres)=180.0d0*deg2rad
 346 c          vbld(nres)=3.8d0
 347 c          vbld_inv(nres)=1.0d0/vbld(2)
 348 c        endif
 349 c      endif
 350       if (lside) then
 351         do i=2,nres-1
 352           do j=1,3
 353             c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
 354      &     +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
 355           enddo
 356           iti=itype(i)
 357           di=dist(i,nres+i)
 358 C 10/03/12 Adam: Correction for zero SC-SC bond length
 359           if (itype(i).ne.10 .and. itype(i).ne.21. and. di.eq.0.0d0)
 360      &     di=dsc(itype(i))
 361           vbld(i+nres)=di
 362           if (itype(i).ne.10) then
 363             vbld_inv(i+nres)=1.0d0/di
 364           else
 365             vbld_inv(i+nres)=0.0d0
 366           endif
 367           if (iti.ne.10) then
 368             alph(i)=alpha(nres+i,i,maxres2)
 369             omeg(i)=beta(nres+i,i,maxres2,i+1)
 370           endif
 371           if(me.eq.king.or..not.out1file)then
 372            if (lprn)
 373      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
 374      &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
 375      &     rad2deg*alph(i),rad2deg*omeg(i)
 376           endif
 377         enddo
 378       else if (lprn) then
 379         do i=2,nres
 380           iti=itype(i)
 381           if(me.eq.king.or..not.out1file)
 382      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 383      &     rad2deg*theta(i),rad2deg*phi(i)
 384         enddo
 385       endif
 386       return
 387       end
 388 c-------------------------------------------------------------------------------
 389       subroutine sc_loc_geom(lprn)
 390       implicit real*8 (a-h,o-z)
 391       include 'DIMENSIONS'
 392 #ifdef MPI
 393       include "mpif.h"
 394 #endif
 395       include 'COMMON.LOCAL'
 396       include 'COMMON.VAR'
 397       include 'COMMON.CHAIN'
 398       include 'COMMON.INTERACT'
 399       include 'COMMON.IOUNITS'
 400       include 'COMMON.GEO'
 401       include 'COMMON.NAMES'
 402       include 'COMMON.CONTROL'
 403       include 'COMMON.SETUP'
 404       double precision x_prime(3),y_prime(3),z_prime(3)
 405       logical lprn
 406       do i=1,nres-1
 407         do j=1,3
 408           dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
 409         enddo
 410       enddo
 411       do i=2,nres-1
 412         if (itype(i).ne.10) then
 413           do j=1,3
 414             dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
 415           enddo
 416         else
 417           do j=1,3
 418             dc_norm(j,i+nres)=0.0d0
 419           enddo
 420         endif
 421       enddo
 422       do i=2,nres-1
 423         costtab(i+1) =dcos(theta(i+1))
 424         sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
 425         cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
 426         sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
 427         cosfac2=0.5d0/(1.0d0+costtab(i+1))
 428         cosfac=dsqrt(cosfac2)
 429         sinfac2=0.5d0/(1.0d0-costtab(i+1))
 430         sinfac=dsqrt(sinfac2)
 431         it=itype(i)
 432         if (it.ne.10) then
 433 c
 434 C  Compute the axes of tghe local cartesian coordinates system; store in
 435 c   x_prime, y_prime and z_prime
 436 c
 437         do j=1,3
 438           x_prime(j) = 0.00
 439           y_prime(j) = 0.00
 440           z_prime(j) = 0.00
 441         enddo
 442         do j = 1,3
 443           x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
 444           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
 445         enddo
 446         call vecpr(x_prime,y_prime,z_prime)
 447 c
 448 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
 449 C to local coordinate system. Store in xx, yy, zz.
 450 c
 451         xx=0.0d0
 452         yy=0.0d0
 453         zz=0.0d0
 454         do j = 1,3
 455           xx = xx + x_prime(j)*dc_norm(j,i+nres)
 456           yy = yy + y_prime(j)*dc_norm(j,i+nres)
 457           zz = zz + z_prime(j)*dc_norm(j,i+nres)
 458         enddo
 459
 460         xxref(i)=xx
 461         yyref(i)=yy
 462         zzref(i)=zz
 463         else
 464         xxref(i)=0.0d0
 465         yyref(i)=0.0d0
 466         zzref(i)=0.0d0
 467         endif
 468       enddo
 469       if (lprn) then
 470         do i=2,nres
 471           iti=itype(i)
 472           if(me.eq.king.or..not.out1file)
 473      &     write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),
 474      &      yyref(i),zzref(i)
 475         enddo
 476       endif
 477       return
 478       end
 479 c---------------------------------------------------------------------------
 480       subroutine sccenter(ires,nscat,sccor)
 481       implicit real*8 (a-h,o-z)
 482       include 'DIMENSIONS'
 483       include 'COMMON.CHAIN'
 484       dimension sccor(3,20)
 485       do j=1,3
 486         sccmj=0.0D0
 487         do i=1,nscat
 488           sccmj=sccmj+sccor(j,i)
 489         enddo
 490         dc(j,ires)=sccmj/nscat
 491       enddo
 492       return
 493       end
 494 c---------------------------------------------------------------------------
 495       subroutine bond_regular
 496       implicit real*8 (a-h,o-z)
 497       include 'DIMENSIONS'
 498       include 'COMMON.VAR'
 499       include 'COMMON.LOCAL'
 500       include 'COMMON.CALC'
 501       include 'COMMON.INTERACT'
 502       include 'COMMON.CHAIN'
 503       do i=1,nres-1
 504        vbld(i+1)=vbl
 505        vbld_inv(i+1)=1.0d0/vbld(i+1)
 506        vbld(i+1+nres)=dsc(itype(i+1))
 507        vbld_inv(i+1+nres)=dsc_inv(itype(i+1))
 508 c       print *,vbld(i+1),vbld(i+1+nres)
 509       enddo
 510       return
 511       end