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