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