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