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