source/unres/src-HCD-5D/readpdb-mult.F

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