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