subroutine contact_cp(var,var2,iff,ieval,in_pdb)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.SBRIDGE'
      include 'COMMON.FFIELD'
      include 'COMMON.IOUNITS'
      include 'COMMON.FRAG'
      include 'COMMON.VAR'
      include 'COMMON.CHAIN'
      include 'COMMON.MINIM'

      character*50 linia
      integer nf,ij(4)
      double precision energy(0:n_ene)
      double precision var(maxvar),var2(maxvar)
      double precision time0,time1
      integer iff(maxres),ieval      
      double precision theta1(maxres),phi1(maxres),alph1(maxres),     
     &                 omeg1(maxres)                             
      logical debug
      
      debug=.false.
c      debug=.true.
      if (ieval.eq.-1) debug=.true.


c
c store selected dist. constrains from 1st structure
c
#ifdef OSF
c     Intercept NaNs in the coordinates
c      write(iout,*) (var(i),i=1,nvar)
      x_sum=0.D0
      do i=1,nvar
        x_sum=x_sum+var(i)
      enddo
      if (x_sum.ne.x_sum) then
        write(iout,*)" *** contact_cp : Found NaN in coordinates"
        call flush(iout) 
        print *," *** contact_cp : Found NaN in coordinates"
        return
      endif
#endif
 

       call var_to_geom(nvar,var)
       call chainbuild                                                           
       nhpb0=nhpb
       ind=0                                                                     
       do i=1,nres-3                                                             
         do j=i+3,nres                                                           
          ind=ind+1                                                              
          if ( iff(i).eq.1.and.iff(j).eq.1 ) then                                           
c            d0(ind)=DIST(i,j)                                                     
c            w(ind)=10.0                                                           
            nhpb=nhpb+1                                                           
            ihpb(nhpb)=i                                                          
            jhpb(nhpb)=j                                                          
            forcon(nhpb)=10.0                                                     
            dhpb(nhpb)=DIST(i,j)
          else
c            w(ind)=0.0
          endif                                                                  
         enddo                                                                   
       enddo                                    
       call hpb_partition

       do i=1,nres                                                               
        theta1(i)=theta(i)                                                      
        phi1(i)=phi(i)                                                          
        alph1(i)=alph(i)                                                        
        omeg1(i)=omeg(i)                                                        
       enddo                      

c
c  freeze sec.elements from 2nd structure 
c
       do i=1,nres
         mask_phi(i)=1
         mask_theta(i)=1
         mask_side(i)=1
       enddo

       call var_to_geom(nvar,var2)
       call secondary2(debug)
       do j=1,nbfrag
        do i=bfrag(1,j),bfrag(2,j)
c         mask(i)=0
         mask_phi(i)=0
         mask_theta(i)=0
        enddo
        if (bfrag(3,j).le.bfrag(4,j)) then 
         do i=bfrag(3,j),bfrag(4,j)
c          mask(i)=0
          mask_phi(i)=0
          mask_theta(i)=0
         enddo
        else
         do i=bfrag(4,j),bfrag(3,j)
c          mask(i)=0
          mask_phi(i)=0
          mask_theta(i)=0
         enddo
        endif
       enddo
       do j=1,nhfrag
        do i=hfrag(1,j),hfrag(2,j)
c         mask(i)=0
         mask_phi(i)=0
         mask_theta(i)=0
        enddo
       enddo
       mask_r=.true.

c
c      copy selected res from 1st to 2nd structure
c

       do i=1,nres                                                             
          if ( iff(i).eq.1 ) then                                           
              theta(i)=theta1(i)                                                      
              phi(i)=phi1(i)                                                          
              alph(i)=alph1(i)                                                        
              omeg(i)=omeg1(i)                       
          endif
       enddo

      if(debug) then   
c
c     prepare description in linia variable
c
        iwsk=0
        nf=0
        if (iff(1).eq.1) then
          iwsk=1
          nf=nf+1
          ij(nf)=1
        endif
        do i=2,nres
           if ( iwsk.eq.0.and.iff(i-1).eq.0.and.iff(i).eq.1 ) then
             iwsk=1
             nf=nf+1
             ij(nf)=i
           endif
           if ( iwsk.eq.1.and.iff(i-1).eq.1.and.iff(i).eq.0 ) then
             iwsk=0
             nf=nf+1
             ij(nf)=i-1
           endif
        enddo
        if (iff(nres).eq.1) then
          nf=nf+1
          ij(nf)=nres
        endif

        write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)')
     &                     "SELECT",ij(1)-1,"-",ij(2)-1,
     &                         ",",ij(3)-1,"-",ij(4)-1

      endif
c
c     run optimization
c
      call contact_cp_min(var,ieval,in_pdb,linia,debug)

      return
      end

      subroutine contact_cp_min(var,ieval,in_pdb,linia,debug)
c
c    input : theta,phi,alph,omeg,in_pdb,linia,debug
c    output : var,ieval
c
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
#ifdef MPI
      include 'mpif.h'
#endif
      include 'COMMON.SBRIDGE'
      include 'COMMON.FFIELD'
      include 'COMMON.IOUNITS'
      include 'COMMON.FRAG'
      include 'COMMON.VAR'
      include 'COMMON.CHAIN'
      include 'COMMON.MINIM'

      character*50 linia
      integer nf,ij(4)
      double precision energy(0:n_ene)
      double precision var(maxvar)
      double precision time0,time1
      integer ieval,info(3)      
      logical debug,fail,check_var,reduce,change

       write(iout,'(a20,i6,a20)')
     &             '------------------',in_pdb,'-------------------'

       if (debug) then
        call chainbuild
        call write_pdb(1000+in_pdb,'combined structure',0d0)
#ifdef MPI
        time0=MPI_WTIME()
#else
        time0=tcpu()
#endif
       endif
       
c
c     run optimization of distances
c     
c     uses d0(),w() and mask() for frozen 2D
c
ctest---------------------------------------------
ctest       NX=NRES-3                                                                 
ctest       NY=((NRES-4)*(NRES-5))/2 
ctest       call distfit(debug,5000)

       do i=1,nres
         mask_side(i)=0
       enddo
 
       ipot01=ipot
       maxmin01=maxmin
       maxfun01=maxfun
c       wstrain01=wstrain
       wsc01=wsc
       wscp01=wscp
       welec01=welec
       wvdwpp01=wvdwpp
c      wang01=wang
       wscloc01=wscloc
       wtor01=wtor
       wtor_d01=wtor_d

       ipot=6
       maxmin=2000
       maxfun=4000
c       wstrain=1.0
       wsc=0.0
       wscp=0.0
       welec=0.0
       wvdwpp=0.0
c      wang=0.0
       wscloc=0.0
       wtor=0.0
       wtor_d=0.0

       call geom_to_var(nvar,var)
cde       change=reduce(var)
       if (check_var(var,info)) then
          write(iout,*) 'cp_min error in input'
          print *,'cp_min error in input'
          return
       endif

cd       call etotal(energy(0))
cd       call enerprint(energy(0))
cd       call check_eint

#ifdef MPI
       time0=MPI_WTIME()
#else
       time0=tcpu()
#endif
cdtest       call minimize(etot,var,iretcode,nfun)                               
cdtest       write(iout,*)'SUMSL return code is',iretcode,' eval SDIST',nfun   
#ifdef MPI
       time1=MPI_WTIME()
#else
       time1=tcpu()
#endif

cd       call etotal(energy(0))
cd       call enerprint(energy(0))
cd       call check_eint 

       do i=1,nres
         mask_side(i)=1
       enddo
 
       ipot=ipot01
       maxmin=maxmin01
       maxfun=maxfun01
c       wstrain=wstrain01
       wsc=wsc01
       wscp=wscp01
       welec=welec01
       wvdwpp=wvdwpp01
c      wang=wang01
       wscloc=wscloc01
       wtor=wtor01
       wtor_d=wtor_d01
ctest--------------------------------------------------
        
       if(debug) then
#ifdef MPI
        time1=MPI_WTIME()
#else
        time1=tcpu()
#endif
        write (iout,'(a,f6.2,a)')' Time for distfit ',time1-time0,' sec'
        call write_pdb(2000+in_pdb,'distfit structure',0d0)
       endif


       ipot0=ipot
       maxmin0=maxmin
       maxfun0=maxfun
       wstrain0=wstrain
c
c      run soft pot. optimization 
c         with constrains:
c             nhpb,ihpb(),jhpb(),forcon(),dhpb() and hpb_partition 
c         and frozen 2D:
c             mask_phi(),mask_theta(),mask_side(),mask_r
c
       ipot=6
       maxmin=2000
       maxfun=4000

cde       change=reduce(var)
cde       if (check_var(var,info)) write(iout,*) 'error before soft'
#ifdef MPI
       time0=MPI_WTIME()
#else
       time0=tcpu()
#endif
       call minimize(etot,var,iretcode,nfun)                               

       write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun   
#ifdef MPI
       time1=MPI_WTIME()
#else
       time1=tcpu()
#endif
       write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,
     &         nfun/(time1-time0),' SOFT eval/s'
       if (debug) then
         call var_to_geom(nvar,var)
         call chainbuild
         call write_pdb(3000+in_pdb,'soft structure',etot)
       endif
c
c      run full UNRES optimization with constrains and frozen 2D
c      the same variables as soft pot. optimizatio
c
       ipot=ipot0
       maxmin=maxmin0
       maxfun=maxfun0
c
c check overlaps before calling full UNRES minim
c
       call var_to_geom(nvar,var)
       call chainbuild
       call etotal(energy(0))
#ifdef OSF
       write(iout,*) 'N7 ',energy(0)
       if (energy(0).ne.energy(0)) then
        write(iout,*) 'N7 error - gives NaN',energy(0)
       endif
#endif
       ieval=1
       if (energy(1).eq.1.0d20) then
         write (iout,'(a,1pe14.5)')'#N7_OVERLAP evdw=1d20',energy(1)
         call overlap_sc(fail)
         if(.not.fail) then
           call etotal(energy(0))
           ieval=ieval+1
           write (iout,'(a,1pe14.5)')'#N7_OVERLAP evdw after',energy(1)
         else
           mask_r=.false.
           nhpb= nhpb0
           link_start=1
           link_end=nhpb
           wstrain=wstrain0
           return
         endif
       endif
       call flush(iout)
c
cdte       time0=MPI_WTIME()
cde       change=reduce(var)
cde       if (check_var(var,info)) then 
cde         write(iout,*) 'error before mask dist'
cde         call var_to_geom(nvar,var)
cde         call chainbuild
cde         call write_pdb(10000+in_pdb,'before mask dist',etot)
cde       endif
cdte       call minimize(etot,var,iretcode,nfun)
cdte       write(iout,*)'SUMSL MASK DIST return code is',iretcode,
cdte     &                          ' eval ',nfun
cdte       ieval=ieval+nfun
cdte
cdte       time1=MPI_WTIME()
cdte       write (iout,'(a,f6.2,f8.2,a)') 
cdte     &        ' Time for mask dist min.',time1-time0,
cdte     &         nfun/(time1-time0),'  eval/s'
cdte       call flush(iout)
       if (debug) then
         call var_to_geom(nvar,var)
         call chainbuild
         call write_pdb(4000+in_pdb,'mask dist',etot)
       endif
c
c      switch off freezing of 2D and 
c      run full UNRES optimization with constrains 
c
       mask_r=.false.
#ifdef MPI
       time0=MPI_WTIME()
#else
       time0=tcpu()
#endif
cde       change=reduce(var)
cde       if (check_var(var,info)) then 
cde         write(iout,*) 'error before dist'
cde         call var_to_geom(nvar,var)
cde         call chainbuild
cde         call write_pdb(11000+in_pdb,'before dist',etot)
cde       endif

       call minimize(etot,var,iretcode,nfun)

cde        change=reduce(var)
cde        if (check_var(var,info)) then 
cde          write(iout,*) 'error after dist',ico
cde          call var_to_geom(nvar,var)
cde          call chainbuild
cde          call write_pdb(12000+in_pdb+ico*1000,'after dist',etot)
cde        endif
       write(iout,*)'SUMSL DIST return code is',iretcode,' eval ',nfun
       ieval=ieval+nfun

#ifdef MPI
       time1=MPI_WTIME()
#else
       time1=tcpu()
#endif
       write (iout,'(a,f6.2,f8.2,a)')' Time for dist min.',time1-time0,
     &         nfun/(time1-time0),'  eval/s'
cde       call etotal(energy(0))
cde       write(iout,*) 'N7 after dist',energy(0)
c       call flush(iout)

       if (debug) then
        call var_to_geom(nvar,var)
        call chainbuild
        call write_pdb(in_pdb,linia,etot)
       endif
c
c      reset constrains
c
       nhpb= nhpb0                                                                 
       link_start=1                                                            
       link_end=nhpb     
       wstrain=wstrain0

       return
       end
c--------------------------------------------------------
      subroutine softreg
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
#ifdef MPI
      include 'mpif.h'
#endif
      include 'COMMON.GEO'
      include 'COMMON.CHAIN'
      include 'COMMON.IOUNITS'
      include 'COMMON.VAR'
      include 'COMMON.CONTROL'
      include 'COMMON.SBRIDGE'
      include 'COMMON.FFIELD'
      include 'COMMON.MINIM'
      include 'COMMON.INTERACT'
c
      include 'COMMON.FRAG'       
      integer iff(maxres)
      double precision time0,time1
      double precision energy(0:n_ene),ee
      double precision var(maxvar)
      integer ieval
c
      logical debug,ltest,fail
      character*50 linia
c
      linia='test'
      debug=.true.
      in_pdb=0



c------------------------
c
c  freeze sec.elements 
c
       do i=1,nres
         mask_phi(i)=1
         mask_theta(i)=1
         mask_side(i)=1
         iff(i)=0
       enddo

       do j=1,nbfrag
        do i=bfrag(1,j),bfrag(2,j)
         mask_phi(i)=0
         mask_theta(i)=0
         iff(i)=1
        enddo
        if (bfrag(3,j).le.bfrag(4,j)) then 
         do i=bfrag(3,j),bfrag(4,j)
          mask_phi(i)=0
          mask_theta(i)=0
          iff(i)=1
         enddo
        else
         do i=bfrag(4,j),bfrag(3,j)
          mask_phi(i)=0
          mask_theta(i)=0
          iff(i)=1
         enddo
        endif
       enddo
       do j=1,nhfrag
        do i=hfrag(1,j),hfrag(2,j)
         mask_phi(i)=0
         mask_theta(i)=0
         iff(i)=1
        enddo
       enddo
       mask_r=.true.



       nhpb0=nhpb
c
c store dist. constrains
c
       do i=1,nres-3                                                             
         do j=i+3,nres                                                           
           if ( iff(i).eq.1.and.iff(j).eq.1 ) then
            nhpb=nhpb+1                                                           
            ihpb(nhpb)=i                                                          
            jhpb(nhpb)=j                                                          
            forcon(nhpb)=0.1                                                     
            dhpb(nhpb)=DIST(i,j)
           endif
         enddo                                                                   
       enddo                                    
       call hpb_partition

       if (debug) then
        call chainbuild
        call write_pdb(100+in_pdb,'input reg. structure',0d0)
       endif
       

       ipot0=ipot
       maxmin0=maxmin
       maxfun0=maxfun
       wstrain0=wstrain
       wang0=wang
c
c      run soft pot. optimization 
c
       ipot=6
       wang=3.0
       maxmin=2000
       maxfun=4000
       call geom_to_var(nvar,var)
#ifdef MPI
       time0=MPI_WTIME()
#else
       time0=tcpu()
#endif
       call minimize(etot,var,iretcode,nfun)                               

       write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun   
#ifdef MPI
       time1=MPI_WTIME()
#else
       time1=tcpu()
#endif
       write (iout,'(a,f6.2,f8.2,a)')'  Time for soft min.',time1-time0,
     &         nfun/(time1-time0),' SOFT eval/s'
       if (debug) then
         call var_to_geom(nvar,var)
         call chainbuild
         call write_pdb(300+in_pdb,'soft structure',etot)
       endif
c
c      run full UNRES optimization with constrains and frozen 2D
c      the same variables as soft pot. optimizatio
c
       ipot=ipot0
       wang=wang0
       maxmin=maxmin0
       maxfun=maxfun0
#ifdef MPI
       time0=MPI_WTIME()
#else
       time0=tcpu()
#endif
       call minimize(etot,var,iretcode,nfun)
       write(iout,*)'SUMSL MASK DIST return code is',iretcode,
     &                          ' eval ',nfun
       ieval=nfun

#ifdef MPI
       time1=MPI_WTIME()
#else
       time1=tcpu()
#endif
       write (iout,'(a,f6.2,f8.2,a)') 
     &        '  Time for mask dist min.',time1-time0,
     &         nfun/(time1-time0),'  eval/s'
       if (debug) then
         call var_to_geom(nvar,var)
         call chainbuild
         call write_pdb(400+in_pdb,'mask & dist',etot)
       endif
c
c      switch off constrains and 
c      run full UNRES optimization with frozen 2D 
c

c
c      reset constrains
c
       nhpb_c=nhpb
       nhpb=nhpb0                                                                  
       link_start=1                                                            
       link_end=nhpb     
       wstrain=wstrain0


#ifdef MPI
       time0=MPI_WTIME()
#else
       time0=tcpu()
#endif
       call minimize(etot,var,iretcode,nfun)
       write(iout,*)'SUMSL MASK return code is',iretcode,' eval ',nfun
       ieval=ieval+nfun

#ifdef MPI
       time1=MPI_WTIME()
#else
       time1=tcpu()
#endif
       write (iout,'(a,f6.2,f8.2,a)')'  Time for mask min.',time1-time0,
     &         nfun/(time1-time0),'  eval/s'


       if (debug) then
        call var_to_geom(nvar,var)
        call chainbuild
        call write_pdb(500+in_pdb,'mask 2d frozen',etot)
       endif

       mask_r=.false.


c
c      run full UNRES optimization with constrains and NO frozen 2D
c

       nhpb=nhpb_c                                                                  
       link_start=1                                                            
       link_end=nhpb     
       maxfun=maxfun0/5

       do ico=1,5

       wstrain=wstrain0/ico

#ifdef MPI
       time0=MPI_WTIME()
#else
       time0=tcpu()
#endif
       call minimize(etot,var,iretcode,nfun)
       write(iout,'(a10,f6.3,a14,i3,a6,i5)')
     &   ' SUMSL DIST',wstrain,' return code is',iretcode,
     &                          ' eval ',nfun
       ieval=nfun

#ifdef MPI
       time1=MPI_WTIME()
#else
       time0=tcpu()
#endif
       write (iout,'(a,f6.2,f8.2,a)') 
     &        '  Time for dist min.',time1-time0,
     &         nfun/(time1-time0),'  eval/s'
       if (debug) then
         call var_to_geom(nvar,var)
         call chainbuild
         call write_pdb(600+in_pdb+ico,'dist cons',etot)
       endif

       enddo
c
       nhpb=nhpb0                                                                  
       link_start=1                                                            
       link_end=nhpb     
       wstrain=wstrain0
       maxfun=maxfun0


c
      if (minim) then
#ifdef MPI
       time0=MPI_WTIME()
#else
       time0=tcpu()
#endif
       call minimize(etot,var,iretcode,nfun)
       write(iout,*)'------------------------------------------------'
       write(iout,*)'SUMSL return code is',iretcode,' eval ',nfun,
     &  '+ DIST eval',ieval
      
#ifdef MPI
       time1=MPI_WTIME()
#else
       time1=tcpu()
#endif
       write (iout,'(a,f6.2,f8.2,a)')'  Time for full min.',time1-time0,
     &         nfun/(time1-time0),' eval/s'


       call var_to_geom(nvar,var)
       call chainbuild        
       call write_pdb(999,'full min',etot)
      endif
       
      return
      end


      subroutine beta_slide(i1,i2,i3,i4,i5,ieval,ij)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
#ifdef MPI
      include 'mpif.h'
#endif
      include 'COMMON.GEO'
      include 'COMMON.VAR'
      include 'COMMON.INTERACT'
      include 'COMMON.IOUNITS'
      include 'COMMON.FRAG'
      include 'COMMON.SBRIDGE'
      include 'COMMON.CONTROL'
      include 'COMMON.FFIELD'
      include 'COMMON.MINIM'
      include 'COMMON.CHAIN'
      double precision time0,time1
      double precision energy(0:n_ene),ee
      double precision var(maxvar)
      integer jdata(5),isec(maxres)
c
      jdata(1)=i1
      jdata(2)=i2
      jdata(3)=i3
      jdata(4)=i4
      jdata(5)=i5

      call secondary2(.false.)

      do i=1,nres
          isec(i)=0
      enddo
      do j=1,nbfrag
       do i=bfrag(1,j),bfrag(2,j)
          isec(i)=1
       enddo
       do i=bfrag(4,j),bfrag(3,j),sign(1,bfrag(3,j)-bfrag(4,j))
          isec(i)=1
       enddo
      enddo
      do j=1,nhfrag
       do i=hfrag(1,j),hfrag(2,j)
          isec(i)=2
       enddo
      enddo

c
c cut strands at the ends
c
      if (jdata(2)-jdata(1).gt.3) then
       jdata(1)=jdata(1)+1
       jdata(2)=jdata(2)-1
       if (jdata(3).lt.jdata(4)) then
          jdata(3)=jdata(3)+1
          jdata(4)=jdata(4)-1
       else
          jdata(3)=jdata(3)-1
          jdata(4)=jdata(4)+1    
       endif
      endif

cv      call chainbuild
cv      call etotal(energy(0))
cv      etot=energy(0)
cv      write(iout,*) nnt,nct,etot
cv      call write_pdb(ij*100,'first structure',etot)
cv      write(iout,*) 'N16 test',(jdata(i),i=1,5)

c------------------------
c      generate constrains 
c
       ishift=jdata(5)-2
       if(ishift.eq.0) ishift=-2
       nhpb0=nhpb
       call chainbuild                                                           
       do i=jdata(1),jdata(2)                                                             
        isec(i)=-1
        if(jdata(4).gt.jdata(3))then
         do j=jdata(3)+i-jdata(1)-2,jdata(3)+i-jdata(1)+2
            isec(j)=-1
cd            print *,i,j,j+ishift
            nhpb=nhpb+1                                                           
            ihpb(nhpb)=i                                                          
            jhpb(nhpb)=j                                                          
            forcon(nhpb)=1000.0                                                     
            dhpb(nhpb)=DIST(i,j+ishift)
         enddo               
        else
         do j=jdata(3)-i+jdata(1)+2,jdata(3)-i+jdata(1)-2,-1
            isec(j)=-1
cd            print *,i,j,j+ishift
            nhpb=nhpb+1                                                           
            ihpb(nhpb)=i                                                          
            jhpb(nhpb)=j                                                          
            forcon(nhpb)=1000.0                                                     
            dhpb(nhpb)=DIST(i,j+ishift)
         enddo
        endif                                                    
       enddo      

       do i=nnt,nct-2
         do j=i+2,nct
           if(isec(i).gt.0.or.isec(j).gt.0) then
cd            print *,i,j
            nhpb=nhpb+1
            ihpb(nhpb)=i
            jhpb(nhpb)=j
            forcon(nhpb)=0.1
            dhpb(nhpb)=DIST(i,j)
           endif
         enddo
       enddo
                              
       call hpb_partition

       call geom_to_var(nvar,var)       
       maxfun0=maxfun
       wstrain0=wstrain
       maxfun=4000/5

       do ico=1,5

        wstrain=wstrain0/ico

cv        time0=MPI_WTIME()
        call minimize(etot,var,iretcode,nfun)
        write(iout,'(a10,f6.3,a14,i3,a6,i5)')
     &   ' SUMSL DIST',wstrain,' return code is',iretcode,
     &                          ' eval ',nfun
        ieval=ieval+nfun
cv        time1=MPI_WTIME()
cv       write (iout,'(a,f6.2,f8.2,a)') 
cv     &        '  Time for dist min.',time1-time0,
cv     &         nfun/(time1-time0),'  eval/s'
cv         call var_to_geom(nvar,var)
cv         call chainbuild
cv         call write_pdb(ij*100+ico,'dist cons',etot)

       enddo
c
       nhpb=nhpb0                                                                  
       call hpb_partition
       wstrain=wstrain0
       maxfun=maxfun0
c
cd      print *,etot
      wscloc0=wscloc
      wscloc=10.0
      call sc_move(nnt,nct,100,100d0,nft_sc,etot)
      wscloc=wscloc0
cv      call chainbuild
cv      call etotal(energy(0))
cv      etot=energy(0)
cv      call write_pdb(ij*100+10,'sc_move',etot)
cd      call intout
cd      print *,nft_sc,etot

      return
      end

      subroutine beta_zip(i1,i2,ieval,ij)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
#ifdef MPI
      include 'mpif.h'
#endif
      include 'COMMON.GEO'
      include 'COMMON.VAR'
      include 'COMMON.INTERACT'
      include 'COMMON.IOUNITS'
      include 'COMMON.FRAG'
      include 'COMMON.SBRIDGE'
      include 'COMMON.CONTROL'
      include 'COMMON.FFIELD'
      include 'COMMON.MINIM'
      include 'COMMON.CHAIN'
      double precision time0,time1
      double precision energy(0:n_ene),ee
      double precision var(maxvar)
      character*10 test

cv      call chainbuild
cv      call etotal(energy(0))
cv      etot=energy(0)
cv      write(test,'(2i5)') i1,i2
cv      call write_pdb(ij*100,test,etot)
cv      write(iout,*) 'N17 test',i1,i2,etot,ij

c
c      generate constrains 
c
       nhpb0=nhpb
       nhpb=nhpb+1                                                           
       ihpb(nhpb)=i1                                                          
       jhpb(nhpb)=i2                                                          
       forcon(nhpb)=1000.0                                                     
       dhpb(nhpb)=4.0
                              
       call hpb_partition

       call geom_to_var(nvar,var)       
       maxfun0=maxfun
       wstrain0=wstrain
       maxfun=1000/5

       do ico=1,5
        wstrain=wstrain0/ico
cv        time0=MPI_WTIME()
        call minimize(etot,var,iretcode,nfun)
        write(iout,'(a10,f6.3,a14,i3,a6,i5)')
     &   ' SUMSL DIST',wstrain,' return code is',iretcode,
     &                          ' eval ',nfun
        ieval=ieval+nfun
cv        time1=MPI_WTIME()
cv       write (iout,'(a,f6.2,f8.2,a)') 
cv     &        '  Time for dist min.',time1-time0,
cv     &         nfun/(time1-time0),'  eval/s'
c do not comment the next line
         call var_to_geom(nvar,var)
cv         call chainbuild
cv         call write_pdb(ij*100+ico,'dist cons',etot)
       enddo

       nhpb=nhpb0                                                                  
       call hpb_partition
       wstrain=wstrain0
       maxfun=maxfun0

cv      call etotal(energy(0))
cv      etot=energy(0)
cv      write(iout,*) 'N17 test end',i1,i2,etot,ij


      return
      end
c----------------------------------------------------------------------------

      subroutine write_pdb(npdb,titelloc,ee)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.IOUNITS'
      character*50 titelloc1                                                     
      character*(*) titelloc
      character*3 zahl   
      character*5 liczba5
      double precision ee
      integer npdb,ilen
      external ilen

      titelloc1=titelloc
      lenpre=ilen(prefix)
      if (npdb.lt.1000) then
       call numstr(npdb,zahl)
       open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
      else
        if (npdb.lt.10000) then                              
         write(liczba5,'(i1,i4)') 0,npdb
        else   
         write(liczba5,'(i5)') npdb
        endif
        open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
      endif
      call pdbout(ee,titelloc1,ipdb)
      close(ipdb)
      return
      end