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