source/wham/src-HCD-5D/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       integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity
  18       logical lprn /.false./,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
  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)=ntyp1
  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 continue
 148 #ifdef DEBUG
 149       write (iout,'(a,i5)') ' Number of residues found: ',ires
 150 #endif
 151       if (ires.eq.0) return
 152 C Calculate the CM of the last side chain.
 153       if (iii.gt.0)  then
 154       if (unres_pdb) then
 155         do j=1,3
 156           dc(j,ires)=sccor(j,iii)
 157         enddo
 158       else
 159         call sccenter(ires,iii,sccor)
 160       endif
 161       endif
 162       nres=ires
 163       nsup=nres
 164       nstart_sup=1
 165       if (itype(nres).ne.10) then
 166         nres=nres+1
 167         itype(nres)=ntyp1
 168         if (unres_pdb) then
 169 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
 170           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
 171           if (fail) then
 172             e2(1)=0.0d0
 173             e2(2)=1.0d0
 174             e2(3)=0.0d0
 175           endif
 176           do j=1,3
 177             c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
 178           enddo
 179         else
 180         do j=1,3
 181           dcj=c(j,nres-2)-c(j,nres-3)
 182           c(j,nres)=c(j,nres-1)+dcj
 183           c(j,2*nres)=c(j,nres)
 184         enddo
 185         endif
 186       endif
 187       do i=2,nres-1
 188         do j=1,3
 189           c(j,i+nres)=dc(j,i)
 190         enddo
 191       enddo
 192       do j=1,3
 193         c(j,nres+1)=c(j,1)
 194         c(j,2*nres)=c(j,nres)
 195       enddo
 196       if (itype(1).eq.ntyp1) then
 197         nsup=nsup-1
 198         nstart_sup=2
 199         if (unres_pdb) then
 200 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 201           call refsys(2,3,4,e1,e2,e3,fail)
 202           if (fail) then
 203             e2(1)=0.0d0
 204             e2(2)=1.0d0
 205             e2(3)=0.0d0
 206           endif
 207           do j=1,3
 208             c(j,1)=c(j,2)-3.8d0*e2(j)
 209           enddo
 210         else
 211         do j=1,3
 212           dcj=c(j,4)-c(j,3)
 213           c(j,1)=c(j,2)-dcj
 214           c(j,nres+1)=c(j,1)
 215         enddo
 216         endif
 217       endif
 218 C Copy the coordinates to reference coordinates
 219 c      do i=1,2*nres
 220 c        do j=1,3
 221 c          cref(j,i)=c(j,i)
 222 c        enddo
 223 c      enddo
 224 C Calculate internal coordinates.
 225       if (lprn) then
 226       write (iout,'(/a)')
 227      &  "Cartesian coordinates of the reference structure"
 228       write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
 229      & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
 230       do ires=1,nres
 231         write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
 232      &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
 233      &    (c(j,ires+nres),j=1,3)
 234       enddo
 235       endif
 236 C Calculate internal coordinates.
 237        write (iout,'(a)')
 238      &   "Backbone and SC coordinates as read from the PDB"
 239        do ires=1,nres
 240         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
 241      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
 242      &    (c(j,nres+ires),j=1,3)
 243        enddo
 244       call int_from_cart(.true.,.false.)
 245       call sc_loc_geom(.false.)
 246       do i=1,nres
 247         thetaref(i)=theta(i)
 248         phiref(i)=phi(i)
 249       enddo
 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 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       include 'COMMON.LOCAL'
 293       include 'COMMON.VAR'
 294       include 'COMMON.CHAIN'
 295       include 'COMMON.INTERACT'
 296       include 'COMMON.IOUNITS'
 297       include 'COMMON.GEO'
 298       include 'COMMON.NAMES'
 299       include 'COMMON.CONTROL'
 300       character*3 seq,res
 301 c      character*5 atom
 302       character*80 card
 303       dimension sccor(3,20)
 304       integer rescode
 305       logical lside,lprn
 306        if (lprn) then
 307         write (iout,'(/a)')
 308      &  'Internal coordinates calculated from crystal structure.'
 309         if (lside) then
 310           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
 311      & '     Gamma','    Dsc_id','       Dsc','     Alpha',
 312      & '     Beta '
 313         else
 314           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
 315      & '     Gamma'
 316         endif
 317        endif
 318       do i=1,nres-1
 319         iti=itype(i)
 320         if (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0) then
 321           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 322 ctest          stop
 323         endif
 324         vbld(i+1)=dist(i,i+1)
 325         vbld_inv(i+1)=1.0d0/vbld(i+1)
 326         if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
 327         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
 328       enddo
 329 c      if (unres_pdb) then
 330 c        if (itype(1).eq.ntyp1) then
 331 c          theta(3)=90.0d0*deg2rad
 332 c          phi(4)=180.0d0*deg2rad
 333 c          vbld(2)=3.8d0
 334 c          vbld_inv(2)=1.0d0/vbld(2)
 335 c        endif
 336 c        if (itype(nres).eq.ntyp1) then
 337 c          theta(nres)=90.0d0*deg2rad
 338 c          phi(nres)=180.0d0*deg2rad
 339 c          vbld(nres)=3.8d0
 340 c          vbld_inv(nres)=1.0d0/vbld(2)
 341 c        endif
 342 c      endif
 343       if (lside) then
 344         do i=2,nres-1
 345           do j=1,3
 346             c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
 347      &     +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
 348           enddo
 349           iti=itype(i)
 350           di=dist(i,nres+i)
 351 C 10/03/12 Adam: Correction for zero SC-SC bond length
 352           if (itype(i).ne.10 .and. itype(i).ne.ntyp1. and. di.eq.0.0d0)
 353      &     di=dsc(itype(i))
 354           vbld(i+nres)=di
 355           if (itype(i).ne.10) then
 356             vbld_inv(i+nres)=1.0d0/di
 357           else
 358             vbld_inv(i+nres)=0.0d0
 359           endif
 360           if (iti.ne.10) then
 361             alph(i)=alpha(nres+i,i,maxres2)
 362             omeg(i)=beta(nres+i,i,maxres2,i+1)
 363           endif
 364            if (lprn)
 365      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
 366      &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
 367      &     rad2deg*alph(i),rad2deg*omeg(i)
 368         enddo
 369       else if (lprn) then
 370         do i=2,nres
 371           iti=itype(i)
 372           write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 373      &     rad2deg*theta(i),rad2deg*phi(i)
 374         enddo
 375       endif
 376       return
 377       end
 378 c-------------------------------------------------------------------------------
 379       subroutine sc_loc_geom(lprn)
 380       implicit real*8 (a-h,o-z)
 381       include 'DIMENSIONS'
 382       include 'DIMENSIONS.ZSCOPT'
 383       include 'COMMON.LOCAL'
 384       include 'COMMON.VAR'
 385       include 'COMMON.CHAIN'
 386       include 'COMMON.INTERACT'
 387       include 'COMMON.IOUNITS'
 388       include 'COMMON.GEO'
 389       include 'COMMON.NAMES'
 390       include 'COMMON.CONTROL'
 391       double precision x_prime(3),y_prime(3),z_prime(3)
 392       logical lprn
 393       do i=1,nres-1
 394         do j=1,3
 395           dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
 396         enddo
 397       enddo
 398       do i=2,nres-1
 399         if (itype(i).ne.10) then
 400           do j=1,3
 401             dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
 402           enddo
 403         else
 404           do j=1,3
 405             dc_norm(j,i+nres)=0.0d0
 406           enddo
 407         endif
 408       enddo
 409       do i=2,nres-1
 410         costtab(i+1) =dcos(theta(i+1))
 411         sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
 412         cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
 413         sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
 414         cosfac2=0.5d0/(1.0d0+costtab(i+1))
 415         cosfac=dsqrt(cosfac2)
 416         sinfac2=0.5d0/(1.0d0-costtab(i+1))
 417         sinfac=dsqrt(sinfac2)
 418         it=itype(i)
 419         if (it.ne.10) then
 420 c
 421 C  Compute the axes of tghe local cartesian coordinates system; store in
 422 c   x_prime, y_prime and z_prime
 423 c
 424         do j=1,3
 425           x_prime(j) = 0.00
 426           y_prime(j) = 0.00
 427           z_prime(j) = 0.00
 428         enddo
 429         do j = 1,3
 430           x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
 431           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
 432         enddo
 433         call vecpr(x_prime,y_prime,z_prime)
 434 c
 435 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
 436 C to local coordinate system. Store in xx, yy, zz.
 437 c
 438         xx=0.0d0
 439         yy=0.0d0
 440         zz=0.0d0
 441         do j = 1,3
 442           xx = xx + x_prime(j)*dc_norm(j,i+nres)
 443           yy = yy + y_prime(j)*dc_norm(j,i+nres)
 444           zz = zz + z_prime(j)*dc_norm(j,i+nres)
 445         enddo
 446
 447         xxref(i)=xx
 448         yyref(i)=yy
 449         zzref(i)=zz
 450         else
 451         xxref(i)=0.0d0
 452         yyref(i)=0.0d0
 453         zzref(i)=0.0d0
 454         endif
 455       enddo
 456       if (lprn) then
 457         do i=2,nres
 458           iti=itype(i)
 459           write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),
 460      &      yyref(i),zzref(i)
 461         enddo
 462       endif
 463       return
 464       end
 465 c---------------------------------------------------------------------------
 466       subroutine sccenter(ires,nscat,sccor)
 467       implicit real*8 (a-h,o-z)
 468       include 'DIMENSIONS'
 469       include 'COMMON.CHAIN'
 470       dimension sccor(3,20)
 471       do j=1,3
 472         sccmj=0.0D0
 473         do i=1,nscat
 474           sccmj=sccmj+sccor(j,i)
 475         enddo
 476         dc(j,ires)=sccmj/nscat
 477       enddo
 478       return
 479       end
 480 c---------------------------------------------------------------------------
 481       subroutine bond_regular
 482       implicit real*8 (a-h,o-z)
 483       include 'DIMENSIONS'
 484       include 'COMMON.VAR'
 485       include 'COMMON.LOCAL'
 486       include 'COMMON.CALC'
 487       include 'COMMON.INTERACT'
 488       include 'COMMON.CHAIN'
 489       do i=1,nres-1
 490        vbld(i+1)=vbl
 491        vbld_inv(i+1)=1.0d0/vbld(i+1)
 492        vbld(i+1+nres)=dsc(itype(i+1))
 493        vbld_inv(i+1+nres)=dsc_inv(itype(i+1))
 494 c       print *,vbld(i+1),vbld(i+1+nres)
 495       enddo
 496       return
 497       end
 498       subroutine readpdb_template(k)
 499 C Read the PDB file with gaps for read_constr_homology with read2sigma
 500 C and convert the peptide geometry into virtual-chain geometry.
 501       implicit real*8 (a-h,o-z)
 502       include 'DIMENSIONS'
 503       include 'DIMENSIONS.ZSCOPT'
 504       include 'DIMENSIONS.FREE'
 505       include 'COMMON.LOCAL'
 506       include 'COMMON.VAR'
 507       include 'COMMON.CHAIN'
 508       include 'COMMON.INTERACT'
 509       include 'COMMON.IOUNITS'
 510       include 'COMMON.GEO'
 511       include 'COMMON.NAMES'
 512       include 'COMMON.CONTROL'
 513       integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity
 514       logical lprn /.false./,fail
 515       double precision e1(3),e2(3),e3(3)
 516       double precision dcj,efree_temp
 517       character*3 seq,res
 518       character*5 atom
 519       character*80 card
 520       double precision sccor(3,20)
 521       integer rescode
 522       efree_temp=0.0d0
 523       ibeg=1
 524       ishift1=0
 525       ishift=0
 526 c      write (2,*) "UNRES_PDB",unres_pdb
 527       ires=0
 528       ires_old=0
 529       iii=0
 530       lsecondary=.false.
 531       nhfrag=0
 532       nbfrag=0
 533       do
 534         read (ipdbin,'(a80)',end=10) card
 535 c        write (iout,'(a)') card
 536         if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10
 537 C Fish out the ATOM cards.
 538         if (index(card(1:4),'ATOM').gt.0) then
 539           read (card(12:16),*) atom
 540 c          write (iout,*) "! ",atom," !",ires
 541 c          if (atom.eq.'CA' .or. atom.eq.'CH3') then
 542           read (card(23:26),*) ires
 543           read (card(18:20),'(a3)') res
 544 c          write (iout,*) "ires",ires,ires-ishift+ishift1,
 545 c     &      " ires_old",ires_old
 546 c          write (iout,*) "ishift",ishift," ishift1",ishift1
 547 c          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
 548           if (ires-ishift+ishift1.ne.ires_old) then
 549 C Calculate the CM of the preceding residue.
 550 c            if (ibeg.eq.0) call sccenter(ires,iii,sccor)
 551             if (ibeg.eq.0) then
 552 c              write (iout,*) "Calculating sidechain center iii",iii
 553               if (unres_pdb) then
 554                 do j=1,3
 555                   dc(j,ires)=sccor(j,iii)
 556                 enddo
 557               else
 558                 call sccenter(ires_old,iii,sccor)
 559               endif
 560               iii=0
 561             endif
 562 C Start new residue.
 563             if (res.eq.'Cl-' .or. res.eq.'Na+') then
 564               ires=ires_old
 565               cycle
 566             else if (ibeg.eq.1) then
 567 c              write (iout,*) "BEG ires",ires
 568               ishift=ires-1
 569               if (res.ne.'GLY' .and. res.ne. 'ACE') then
 570                 ishift=ishift-1
 571                 itype(1)=ntyp1
 572               endif
 573               ires=ires-ishift+ishift1
 574               ires_old=ires
 575 c              write (iout,*) "ishift",ishift," ires",ires,
 576 c     &         " ires_old",ires_old
 577               ibeg=0
 578             else
 579               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
 580               ires=ires-ishift+ishift1
 581               ires_old=ires
 582             endif
 583             if (res.eq.'ACE' .or. res.eq.'NHE') then
 584               itype(ires)=10
 585             else
 586               itype(ires)=rescode(ires,res,0)
 587             endif
 588           else
 589             ires=ires-ishift+ishift1
 590           endif
 591 c          write (iout,*) "ires_old",ires_old," ires",ires
 592           if (card(27:27).eq."A" .or. card(27:27).eq."B") then
 593 c            ishift1=ishift1+1
 594           endif
 595 c          write (2,*) "ires",ires," res ",res," ity",ity
 596           if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
 597      &       res.eq.'NHE'.and.atom(:2).eq.'HN') then
 598             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
 599 c            write (iout,*) "backbone ",atom
 600 #ifdef DEBUG
 601             write (iout,'(2i3,2x,a,3f8.3)')
 602      &      ires,itype(ires),res,(c(j,ires),j=1,3)
 603 #endif
 604             iii=iii+1
 605             do j=1,3
 606               sccor(j,iii)=c(j,ires)
 607             enddo
 608             if (ishift.ne.0) then
 609               ires_ca=ires+ishift-ishift1
 610             else
 611               ires_ca=ires
 612             endif
 613 c            write (*,*) card(23:27),ires,itype(ires)
 614           else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.
 615      &             atom.ne.'N' .and. atom.ne.'C' .and.
 616      &             atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.
 617      &             atom.ne.'OXT' .and. atom(:2).ne.'3H') then
 618 c            write (iout,*) "sidechain ",atom
 619             iii=iii+1
 620             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
 621           endif
 622         endif
 623       enddo
 624    10 continue
 625 #ifdef DEBUG
 626       write (iout,'(a,i5)') ' Number of residues found: ',ires
 627 #endif
 628       if (ires.eq.0) return
 629 C Calculate the CM of the last side chain.
 630       if (iii.gt.0)  then
 631       if (unres_pdb) then
 632         do j=1,3
 633           dc(j,ires)=sccor(j,iii)
 634         enddo
 635       else
 636         call sccenter(ires,iii,sccor)
 637       endif
 638       endif
 639       nres=ires
 640       nsup=nres
 641       nstart_sup=1
 642       if (itype(nres).ne.10) then
 643         nres=nres+1
 644         itype(nres)=ntyp1
 645         do j=1,3
 646           dcj=c(j,nres-2)-c(j,nres-3)
 647           c(j,nres)=c(j,nres-1)+dcj
 648           c(j,2*nres)=c(j,nres)
 649         enddo
 650       endif
 651       do i=2,nres-1
 652         do j=1,3
 653           c(j,i+nres)=dc(j,i)
 654         enddo
 655       enddo
 656       do j=1,3
 657         c(j,nres+1)=c(j,1)
 658         c(j,2*nres)=c(j,nres)
 659       enddo
 660       if (itype(1).eq.ntyp1) then
 661         nsup=nsup-1
 662         nstart_sup=2
 663         if (unres_pdb) then
 664 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 665           call refsys(2,3,4,e1,e2,e3,fail)
 666           if (fail) then
 667             e2(1)=0.0d0
 668             e2(2)=1.0d0
 669             e2(3)=0.0d0
 670           endif
 671           do j=1,3
 672             c(j,1)=c(j,2)-3.8d0*e2(j)
 673           enddo
 674         else
 675         do j=1,3
 676           dcj=c(j,4)-c(j,3)
 677           c(j,1)=c(j,2)-dcj
 678           c(j,nres+1)=c(j,1)
 679         enddo
 680         endif
 681       endif
 682 C Calculate internal coordinates.
 683       if (lprn) then
 684       write (iout,'(/a)')
 685      &  "Cartesian coordinates of the reference structure"
 686       write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
 687      & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
 688       do ires=1,nres
 689         write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
 690      &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
 691      &    (c(j,ires+nres),j=1,3)
 692       enddo
 693       endif
 694 C Calculate internal coordinates.
 695        write (iout,'(a)')
 696      &   "Backbone and SC coordinates as read from the PDB"
 697        do ires=1,nres
 698         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
 699      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
 700      &    (c(j,nres+ires),j=1,3)
 701        enddo
 702       call int_from_cart(.true.,.false.)
 703       call sc_loc_geom(.false.)
 704       do i=1,nres
 705         thetaref(i)=theta(i)
 706         phiref(i)=phi(i)
 707       enddo
 708       do i=1,nres-1
 709         do j=1,3
 710           dc(j,i)=c(j,i+1)-c(j,i)
 711           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 712         enddo
 713       enddo
 714       do i=2,nres-1
 715         do j=1,3
 716           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 717           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 718         enddo
 719 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
 720 c     &   vbld_inv(i+nres)
 721       enddo
 722 c      call chainbuild
 723 C Copy the coordinates to reference coordinates
 724       do i=1,2*nres
 725         do j=1,3
 726           cref(j,i)=c(j,i)
 727           chomo(j,i,k)=c(j,i)
 728         enddo
 729       enddo
 730
 731
 732       ishift_pdb=ishift
 733       return
 734       end