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