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