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