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         do j=1,3
 169           dcj=c(j,nres-2)-c(j,nres-3)
 170           c(j,nres)=c(j,nres-1)+dcj
 171           c(j,2*nres)=c(j,nres)
 172         enddo
 173       endif
 174       do i=2,nres-1
 175         do j=1,3
 176           c(j,i+nres)=dc(j,i)
 177         enddo
 178       enddo
 179       do j=1,3
 180         c(j,nres+1)=c(j,1)
 181         c(j,2*nres)=c(j,nres)
 182       enddo
 183       if (itype(1).eq.21) then
 184         nsup=nsup-1
 185         nstart_sup=2
 186         if (unres_pdb) then
 187 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 188           call refsys(2,3,4,e1,e2,e3,fail)
 189           if (fail) then
 190             e2(1)=0.0d0
 191             e2(2)=1.0d0
 192             e2(3)=0.0d0
 193           endif
 194           do j=1,3
 195             c(j,1)=c(j,2)-3.8d0*e2(j)
 196           enddo
 197         else
 198         do j=1,3
 199           dcj=c(j,4)-c(j,3)
 200           c(j,1)=c(j,2)-dcj
 201           c(j,nres+1)=c(j,1)
 202         enddo
 203         endif
 204       endif
 205 C Copy the coordinates to reference coordinates
 206 c      do i=1,2*nres
 207 c        do j=1,3
 208 c          cref(j,i)=c(j,i)
 209 c        enddo
 210 c      enddo
 211 C Calculate internal coordinates.
 212       if (lprn) then
 213       write (iout,'(/a)')
 214      &  "Cartesian coordinates of the reference structure"
 215       write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
 216      & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
 217       do ires=1,nres
 218         write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
 219      &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
 220      &    (c(j,ires+nres),j=1,3)
 221       enddo
 222       endif
 223 C Calculate internal coordinates.
 224       if(me.eq.king.or..not.out1file)then
 225        write (iout,'(a)')
 226      &   "Backbone and SC coordinates as read from the PDB"
 227        do ires=1,nres
 228         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
 229      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
 230      &    (c(j,nres+ires),j=1,3)
 231        enddo
 232       endif
 233       call int_from_cart(.true.,.false.)
 234       call sc_loc_geom(.false.)
 235       do i=1,nres
 236         thetaref(i)=theta(i)
 237         phiref(i)=phi(i)
 238       enddo
 239       do i=1,nres-1
 240         do j=1,3
 241           dc(j,i)=c(j,i+1)-c(j,i)
 242           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 243         enddo
 244       enddo
 245       do i=2,nres-1
 246         do j=1,3
 247           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 248           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 249         enddo
 250 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
 251 c     &   vbld_inv(i+nres)
 252       enddo
 253 c      call chainbuild
 254 C Copy the coordinates to reference coordinates
 255       do i=1,2*nres
 256         do j=1,3
 257           cref(j,i)=c(j,i)
 258         enddo
 259       enddo
 260
 261
 262       do j=1,nbfrag
 263         do i=1,4
 264          bfrag(i,j)=bfrag(i,j)-ishift
 265         enddo
 266       enddo
 267
 268       do j=1,nhfrag
 269         do i=1,2
 270          hfrag(i,j)=hfrag(i,j)-ishift
 271         enddo
 272       enddo
 273       ishift_pdb=ishift
 274       return
 275       end
 276 c---------------------------------------------------------------------------
 277       subroutine int_from_cart(lside,lprn)
 278       implicit real*8 (a-h,o-z)
 279       include 'DIMENSIONS'
 280 #ifdef MPI
 281       include "mpif.h"
 282 #endif
 283       include 'COMMON.LOCAL'
 284       include 'COMMON.VAR'
 285       include 'COMMON.CHAIN'
 286       include 'COMMON.INTERACT'
 287       include 'COMMON.IOUNITS'
 288       include 'COMMON.GEO'
 289       include 'COMMON.NAMES'
 290       include 'COMMON.CONTROL'
 291       include 'COMMON.SETUP'
 292       character*3 seq,res
 293 c      character*5 atom
 294       character*80 card
 295       dimension sccor(3,20)
 296       integer rescode
 297       logical lside,lprn
 298       if(me.eq.king.or..not.out1file)then
 299        if (lprn) then
 300         write (iout,'(/a)')
 301      &  'Internal coordinates calculated from crystal structure.'
 302         if (lside) then
 303           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
 304      & '     Gamma','    Dsc_id','       Dsc','     Alpha',
 305      & '     Beta '
 306         else
 307           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
 308      & '     Gamma'
 309         endif
 310        endif
 311       endif
 312       do i=1,nres-1
 313         iti=itype(i)
 314         if (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0) then
 315           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 316 ctest          stop
 317         endif
 318         vbld(i+1)=dist(i,i+1)
 319         vbld_inv(i+1)=1.0d0/vbld(i+1)
 320         if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
 321         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
 322       enddo
 323 c      if (unres_pdb) then
 324 c        if (itype(1).eq.21) then
 325 c          theta(3)=90.0d0*deg2rad
 326 c          phi(4)=180.0d0*deg2rad
 327 c          vbld(2)=3.8d0
 328 c          vbld_inv(2)=1.0d0/vbld(2)
 329 c        endif
 330 c        if (itype(nres).eq.21) then
 331 c          theta(nres)=90.0d0*deg2rad
 332 c          phi(nres)=180.0d0*deg2rad
 333 c          vbld(nres)=3.8d0
 334 c          vbld_inv(nres)=1.0d0/vbld(2)
 335 c        endif
 336 c      endif
 337       if (lside) then
 338         do i=2,nres-1
 339           do j=1,3
 340             c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
 341      &     +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
 342           enddo
 343           iti=itype(i)
 344           di=dist(i,nres+i)
 345 C 10/03/12 Adam: Correction for zero SC-SC bond length
 346           if (itype(i).ne.10 .and. itype(i).ne.21. and. di.eq.0.0d0)
 347      &     di=dsc(itype(i))
 348           vbld(i+nres)=di
 349           if (itype(i).ne.10) then
 350             vbld_inv(i+nres)=1.0d0/di
 351           else
 352             vbld_inv(i+nres)=0.0d0
 353           endif
 354           if (iti.ne.10) then
 355             alph(i)=alpha(nres+i,i,maxres2)
 356             omeg(i)=beta(nres+i,i,maxres2,i+1)
 357           endif
 358           if(me.eq.king.or..not.out1file)then
 359            if (lprn)
 360      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
 361      &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
 362      &     rad2deg*alph(i),rad2deg*omeg(i)
 363           endif
 364         enddo
 365       else if (lprn) then
 366         do i=2,nres
 367           iti=itype(i)
 368           if(me.eq.king.or..not.out1file)
 369      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 370      &     rad2deg*theta(i),rad2deg*phi(i)
 371         enddo
 372       endif
 373       return
 374       end
 375 c-------------------------------------------------------------------------------
 376       subroutine sc_loc_geom(lprn)
 377       implicit real*8 (a-h,o-z)
 378       include 'DIMENSIONS'
 379 #ifdef MPI
 380       include "mpif.h"
 381 #endif
 382       include 'COMMON.LOCAL'
 383       include 'COMMON.VAR'
 384       include 'COMMON.CHAIN'
 385       include 'COMMON.INTERACT'
 386       include 'COMMON.IOUNITS'
 387       include 'COMMON.GEO'
 388       include 'COMMON.NAMES'
 389       include 'COMMON.CONTROL'
 390       include 'COMMON.SETUP'
 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           if(me.eq.king.or..not.out1file)
 460      &     write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),
 461      &      yyref(i),zzref(i)
 462         enddo
 463       endif
 464       return
 465       end
 466 c---------------------------------------------------------------------------
 467       subroutine sccenter(ires,nscat,sccor)
 468       implicit real*8 (a-h,o-z)
 469       include 'DIMENSIONS'
 470       include 'COMMON.CHAIN'
 471       dimension sccor(3,20)
 472       do j=1,3
 473         sccmj=0.0D0
 474         do i=1,nscat
 475           sccmj=sccmj+sccor(j,i)
 476         enddo
 477         dc(j,ires)=sccmj/nscat
 478       enddo
 479       return
 480       end
 481 c---------------------------------------------------------------------------
 482       subroutine bond_regular
 483       implicit real*8 (a-h,o-z)
 484       include 'DIMENSIONS'
 485       include 'COMMON.VAR'
 486       include 'COMMON.LOCAL'
 487       include 'COMMON.CALC'
 488       include 'COMMON.INTERACT'
 489       include 'COMMON.CHAIN'
 490       do i=1,nres-1
 491        vbld(i+1)=vbl
 492        vbld_inv(i+1)=1.0d0/vbld(i+1)
 493        vbld(i+1+nres)=dsc(itype(i+1))
 494        vbld_inv(i+1+nres)=dsc_inv(itype(i+1))
 495 c       print *,vbld(i+1),vbld(i+1+nres)
 496       enddo
 497       return
 498       end
 499       subroutine readpdb_template(k)
 500 C Read the PDB file with gaps for read_constr_homology with read2sigma
 501 C and convert the peptide geometry into virtual-chain geometry.
 502       implicit real*8 (a-h,o-z)
 503       include 'DIMENSIONS'
 504       include 'COMMON.LOCAL'
 505       include 'COMMON.VAR'
 506       include 'COMMON.CHAIN'
 507       include 'COMMON.INTERACT'
 508       include 'COMMON.IOUNITS'
 509       include 'COMMON.GEO'
 510       include 'COMMON.NAMES'
 511       include 'COMMON.CONTROL'
 512       include 'COMMON.DISTFIT'
 513       include 'COMMON.SETUP'
 514       include 'COMMON.MD'
 515       integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity,
 516      &  ishift_pdb
 517       logical lprn /.false./,fail
 518       double precision e1(3),e2(3),e3(3)
 519       double precision dcj,efree_temp
 520       character*3 seq,res
 521       character*5 atom
 522       character*80 card
 523       double precision sccor(3,20)
 524       integer rescode
 525       efree_temp=0.0d0
 526       ibeg=1
 527       ishift1=0
 528       ishift=0
 529 c      write (2,*) "UNRES_PDB",unres_pdb
 530       ires=0
 531       ires_old=0
 532       iii=0
 533       lsecondary=.false.
 534       nhfrag=0
 535       nbfrag=0
 536       do i=1,10000
 537         read (ipdbin,'(a80)',end=10) card
 538 c        write (iout,'(a)') card
 539         if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10
 540 C Fish out the ATOM cards.
 541         if (index(card(1:4),'ATOM').gt.0) then
 542           read (card(12:16),*) atom
 543 c          write (iout,*) "! ",atom," !",ires
 544 c          if (atom.eq.'CA' .or. atom.eq.'CH3') then
 545           read (card(23:26),*) ires
 546           read (card(18:20),'(a3)') res
 547 c          write (iout,*) "ires",ires,ires-ishift+ishift1,
 548 c     &      " ires_old",ires_old
 549 c          write (iout,*) "ishift",ishift," ishift1",ishift1
 550 c          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
 551           if (ires-ishift+ishift1.ne.ires_old) then
 552 C Calculate the CM of the preceding residue.
 553 c            if (ibeg.eq.0) call sccenter(ires,iii,sccor)
 554             if (ibeg.eq.0) then
 555 c              write (iout,*) "Calculating sidechain center iii",iii
 556               if (unres_pdb) then
 557                 do j=1,3
 558                   dc(j,ires)=sccor(j,iii)
 559                 enddo
 560               else
 561                 call sccenter(ires_old,iii,sccor)
 562               endif
 563               iii=0
 564             endif
 565 C Start new residue.
 566             if (res.eq.'Cl-' .or. res.eq.'Na+') then
 567               ires=ires_old
 568               cycle
 569             else if (ibeg.eq.1) then
 570 c              write (iout,*) "BEG ires",ires
 571               ishift=ires-1
 572               if (res.ne.'GLY' .and. res.ne. 'ACE') then
 573                 ishift=ishift-1
 574                 itype(1)=21
 575               endif
 576               ires=ires-ishift+ishift1
 577               ires_old=ires
 578 c              write (iout,*) "ishift",ishift," ires",ires,
 579 c     &         " ires_old",ires_old
 580               ibeg=0
 581             else
 582               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
 583               ires=ires-ishift+ishift1
 584               ires_old=ires
 585             endif
 586             if (res.eq.'ACE' .or. res.eq.'NHE') then
 587               itype(ires)=10
 588             else
 589               itype(ires)=rescode(ires,res,0)
 590             endif
 591           else
 592             ires=ires-ishift+ishift1
 593           endif
 594 c          write (iout,*) "ires_old",ires_old," ires",ires
 595           if (card(27:27).eq."A" .or. card(27:27).eq."B") then
 596 c            ishift1=ishift1+1
 597           endif
 598 c          write (2,*) "ires",ires," res ",res," ity",ity
 599           if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
 600      &       res.eq.'NHE'.and.atom(:2).eq.'HN') then
 601             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
 602 c            write (iout,*) "backbone ",atom
 603 #ifdef DEBUG
 604             write (iout,'(2i3,2x,a,3f8.3)')
 605      &      ires,itype(ires),res,(c(j,ires),j=1,3)
 606 #endif
 607             iii=iii+1
 608             do j=1,3
 609               sccor(j,iii)=c(j,ires)
 610             enddo
 611             if (ishift.ne.0) then
 612               ires_ca=ires+ishift-ishift1
 613             else
 614               ires_ca=ires
 615             endif
 616 c            write (*,*) card(23:27),ires,itype(ires)
 617           else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.
 618      &             atom.ne.'N' .and. atom.ne.'C' .and.
 619      &             atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.
 620      &             atom.ne.'OXT' .and. atom(:2).ne.'3H') then
 621 c            write (iout,*) "sidechain ",atom
 622             iii=iii+1
 623             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
 624           endif
 625         endif
 626       enddo
 627    10 continue
 628 #ifdef DEBUG
 629       write (iout,'(a,i5)') ' Number of residues found: ',ires
 630 #endif
 631       if (ires.eq.0) return
 632 C Calculate the CM of the last side chain.
 633       if (iii.gt.0)  then
 634       if (unres_pdb) then
 635         do j=1,3
 636           dc(j,ires)=sccor(j,iii)
 637         enddo
 638       else
 639         call sccenter(ires,iii,sccor)
 640       endif
 641       endif
 642       nres=ires
 643       nsup=nres
 644       nstart_sup=1
 645       if (itype(nres).ne.10) then
 646         nres=nres+1
 647         itype(nres)=21
 648         do j=1,3
 649           dcj=c(j,nres-2)-c(j,nres-3)
 650           c(j,nres)=c(j,nres-1)+dcj
 651           c(j,2*nres)=c(j,nres)
 652         enddo
 653       endif
 654       do i=2,nres-1
 655         do j=1,3
 656           c(j,i+nres)=dc(j,i)
 657         enddo
 658       enddo
 659       do j=1,3
 660         c(j,nres+1)=c(j,1)
 661         c(j,2*nres)=c(j,nres)
 662       enddo
 663       if (itype(1).eq.21) then
 664         nsup=nsup-1
 665         nstart_sup=2
 666         if (unres_pdb) then
 667 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 668           call refsys(2,3,4,e1,e2,e3,fail)
 669           if (fail) then
 670             e2(1)=0.0d0
 671             e2(2)=1.0d0
 672             e2(3)=0.0d0
 673           endif
 674           do j=1,3
 675             c(j,1)=c(j,2)-3.8d0*e2(j)
 676           enddo
 677         else
 678         do j=1,3
 679           dcj=c(j,4)-c(j,3)
 680           c(j,1)=c(j,2)-dcj
 681           c(j,nres+1)=c(j,1)
 682         enddo
 683         endif
 684       endif
 685 C Calculate internal coordinates.
 686       if (lprn) then
 687       write (iout,'(/a)')
 688      &  "Cartesian coordinates of the reference structure"
 689       write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
 690      & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
 691       do ires=1,nres
 692         write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
 693      &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
 694      &    (c(j,ires+nres),j=1,3)
 695       enddo
 696       endif
 697 C Calculate internal coordinates.
 698       if(me.eq.king.or..not.out1file)then
 699        write (iout,'(a)')
 700      &   "Backbone and SC coordinates as read from the PDB"
 701        do ires=1,nres
 702         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
 703      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
 704      &    (c(j,nres+ires),j=1,3)
 705        enddo
 706       endif
 707       call int_from_cart(.true.,.false.)
 708       call sc_loc_geom(.false.)
 709       do i=1,nres
 710         thetaref(i)=theta(i)
 711         phiref(i)=phi(i)
 712       enddo
 713       do i=1,nres-1
 714         do j=1,3
 715           dc(j,i)=c(j,i+1)-c(j,i)
 716           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 717         enddo
 718       enddo
 719       do i=2,nres-1
 720         do j=1,3
 721           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 722           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 723         enddo
 724 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
 725 c     &   vbld_inv(i+nres)
 726       enddo
 727 c      call chainbuild
 728 C Copy the coordinates to reference coordinates
 729       do i=1,2*nres
 730         do j=1,3
 731           cref(j,i)=c(j,i)
 732         enddo
 733       enddo
 734
 735
 736       ishift_pdb=ishift
 737       return
 738       end