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