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