subroutine sc_move(n_start,n_end,n_maxtry,e_drop,
     +     n_fun,etot)
c     Perform a quick search over side-chain arrangments (over
c     residues n_start to n_end) for a given (frozen) CA trace
c     Only side-chains are minimized (at most n_maxtry times each),
c     not CA positions
c     Stops if energy drops by e_drop, otherwise tries all residues
c     in the given range
c     If there is an energy drop, full minimization may be useful
c     n_start, n_end CAN be modified by this routine, but only if
c     out of bounds (n_start <= 1, n_end >= nres, n_start < n_end)
c     NOTE: this move should never increase the energy
crc      implicit none

c     Includes
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
#ifdef MPI
      include 'mpif.h'
#endif
      include 'COMMON.GEO'
      include 'COMMON.VAR'
      include 'COMMON.HEADER'
      include 'COMMON.IOUNITS'
      include 'COMMON.CHAIN'
      include 'COMMON.FFIELD'

c     External functions
      integer iran_num
      external iran_num

c     Input arguments
      integer n_start,n_end,n_maxtry
      double precision e_drop

c     Output arguments
      integer n_fun
      double precision etot

c     Local variables
      double precision energy(0:n_ene)
      double precision cur_alph(2:nres-1),cur_omeg(2:nres-1)
      double precision orig_e,cur_e
      integer n,n_steps,n_first,n_cur,n_tot,i
      double precision orig_w(n_ene)
      double precision wtime


c     Set non side-chain weights to zero (minimization is faster)
c     NOTE: e(2) does not actually depend on the side-chain, only CA
      orig_w(2)=wscp
      orig_w(3)=welec
      orig_w(4)=wcorr
      orig_w(5)=wcorr5
      orig_w(6)=wcorr6
      orig_w(7)=wel_loc
      orig_w(8)=wturn3
      orig_w(9)=wturn4
      orig_w(10)=wturn6
      orig_w(11)=wang
      orig_w(13)=wtor
      orig_w(14)=wtor_d
      orig_w(15)=wvdwpp

      wscp=0.D0
      welec=0.D0
      wcorr=0.D0
      wcorr5=0.D0
      wcorr6=0.D0
      wel_loc=0.D0
      wturn3=0.D0
      wturn4=0.D0
      wturn6=0.D0
      wang=0.D0
      wtor=0.D0
      wtor_d=0.D0
      wvdwpp=0.D0

c     Make sure n_start, n_end are within proper range
      if (n_start.lt.2) n_start=2
      if (n_end.gt.nres-1) n_end=nres-1
crc      if (n_start.lt.n_end) then
      if (n_start.gt.n_end) then
        n_start=2
        n_end=nres-1
      endif

c     Save the initial values of energy and coordinates
cd      call chainbuild
cd      call etotal(energy)
cd      write (iout,*) 'start sc ene',energy(0)
cd      call enerprint(energy(0))
crc      etot=energy(0)
       n_fun=0
crc      orig_e=etot
crc      cur_e=orig_e
crc      do i=2,nres-1
crc        cur_alph(i)=alph(i)
crc        cur_omeg(i)=omeg(i)
crc      enddo

ct      wtime=MPI_WTIME()
c     Try (one by one) all specified residues, starting from a
c     random position in sequence
c     Stop early if the energy has decreased by at least e_drop
      n_tot=n_end-n_start+1
      n_first=iran_num(0,n_tot-1)
      n_steps=0
      n=0
crc      do while (n.lt.n_tot .and. orig_e-etot.lt.e_drop)
      do while (n.lt.n_tot)
        n_cur=n_start+mod(n_first+n,n_tot)
        call single_sc_move(n_cur,n_maxtry,e_drop,
     +       n_steps,n_fun,etot)
c     If a lower energy was found, update the current structure...
crc        if (etot.lt.cur_e) then
crc          cur_e=etot
crc          do i=2,nres-1
crc            cur_alph(i)=alph(i)
crc            cur_omeg(i)=omeg(i)
crc          enddo
crc        else
c     ...else revert to the previous one
crc          etot=cur_e
crc          do i=2,nres-1
crc            alph(i)=cur_alph(i)
crc            omeg(i)=cur_omeg(i)
crc          enddo
crc        endif
        n=n+1
cd
cd      call chainbuild
cd      call etotal(energy)
cd      print *,'running',n,energy(0)
      enddo

cd      call chainbuild
cd      call etotal(energy)
cd      write (iout,*) 'end   sc ene',energy(0)

c     Put the original weights back to calculate the full energy
      wscp=orig_w(2)
      welec=orig_w(3)
      wcorr=orig_w(4)
      wcorr5=orig_w(5)
      wcorr6=orig_w(6)
      wel_loc=orig_w(7)
      wturn3=orig_w(8)
      wturn4=orig_w(9)
      wturn6=orig_w(10)
      wang=orig_w(11)
      wtor=orig_w(13)
      wtor_d=orig_w(14)
      wvdwpp=orig_w(15)

crc      n_fun=n_fun+1
ct      write (iout,*) 'sc_local time= ',MPI_WTIME()-wtime
      return
      end

c-------------------------------------------------------------

      subroutine single_sc_move(res_pick,n_maxtry,e_drop,
     +     n_steps,n_fun,e_sc)
c     Perturb one side-chain (res_pick) and minimize the
c     neighbouring region, keeping all CA's and non-neighbouring
c     side-chains fixed
c     Try until e_drop energy improvement is achieved, or n_maxtry
c     attempts have been made
c     At the start, e_sc should contain the side-chain-only energy(0)
c     nsteps and nfun for this move are ADDED to n_steps and n_fun
crc      implicit none

c     Includes
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.VAR'
      include 'COMMON.INTERACT'
      include 'COMMON.CHAIN'
      include 'COMMON.MINIM'
      include 'COMMON.FFIELD'
      include 'COMMON.IOUNITS'

c     External functions
      double precision dist
      external dist

c     Input arguments
      integer res_pick,n_maxtry
      double precision e_drop

c     Input/Output arguments
      integer n_steps,n_fun
      double precision e_sc

c     Local variables
      logical fail
      integer i,j
      integer nres_moved
      integer iretcode,loc_nfun,orig_maxfun,n_try
      double precision sc_dist,sc_dist_cutoff
      double precision energy(0:n_ene),orig_e,cur_e
      double precision evdw,escloc
      double precision cur_alph(2:nres-1),cur_omeg(2:nres-1)
      double precision var(maxvar)

      double precision orig_theta(1:nres),orig_phi(1:nres),
     +     orig_alph(1:nres),orig_omeg(1:nres)


c     Define what is meant by "neighbouring side-chain"
      sc_dist_cutoff=5.0D0

c     Don't do glycine or ends
      i=itype(res_pick)
      if (i.eq.10 .or. i.eq.21) return

c     Freeze everything (later will relax only selected side-chains)
      mask_r=.true.
      do i=1,nres
        mask_phi(i)=0
        mask_theta(i)=0
        mask_side(i)=0
      enddo

c     Find the neighbours of the side-chain to move
c     and save initial variables
crc      orig_e=e_sc
crc      cur_e=orig_e
      nres_moved=0
      do i=2,nres-1
c     Don't do glycine (itype(j)==10)
        if (itype(i).ne.10) then
          sc_dist=dist(nres+i,nres+res_pick)
        else
          sc_dist=sc_dist_cutoff
        endif
        if (sc_dist.lt.sc_dist_cutoff) then
          nres_moved=nres_moved+1
          mask_side(i)=1
          cur_alph(i)=alph(i)
          cur_omeg(i)=omeg(i)
        endif
      enddo

      call chainbuild
      call egb1(evdw)
      call esc(escloc)
      e_sc=wsc*evdw+wscloc*escloc
cd      call etotal(energy)
cd      print *,'new       ',(energy(k),k=0,n_ene)
      orig_e=e_sc
      cur_e=orig_e

      n_try=0
      do while (n_try.lt.n_maxtry .and. orig_e-cur_e.lt.e_drop)
c     Move the selected residue (don't worry if it fails)
        call gen_side(itype(res_pick),theta(res_pick+1),
     +       alph(res_pick),omeg(res_pick),fail)

c     Minimize the side-chains starting from the new arrangement
        call geom_to_var(nvar,var)
        orig_maxfun=maxfun
        maxfun=7

crc        do i=1,nres
crc          orig_theta(i)=theta(i)
crc          orig_phi(i)=phi(i)
crc          orig_alph(i)=alph(i)
crc          orig_omeg(i)=omeg(i)
crc        enddo

        call minimize_sc1(e_sc,var,iretcode,loc_nfun)
        
cv        write(*,'(2i3,2f12.5,2i3)') 
cv     &       res_pick,nres_moved,orig_e,e_sc-cur_e,
cv     &       iretcode,loc_nfun

c$$$        if (iretcode.eq.8) then
c$$$          write(iout,*)'Coordinates just after code 8'
c$$$          call chainbuild
c$$$          call all_varout
c$$$          call flush(iout)
c$$$          do i=1,nres
c$$$            theta(i)=orig_theta(i)
c$$$            phi(i)=orig_phi(i)
c$$$            alph(i)=orig_alph(i)
c$$$            omeg(i)=orig_omeg(i)
c$$$          enddo
c$$$          write(iout,*)'Coordinates just before code 8'
c$$$          call chainbuild
c$$$          call all_varout
c$$$          call flush(iout)
c$$$        endif

        n_fun=n_fun+loc_nfun
        maxfun=orig_maxfun
        call var_to_geom(nvar,var)

c     If a lower energy was found, update the current structure...
        if (e_sc.lt.cur_e) then
cv              call chainbuild
cv              call etotal(energy)
cd              call egb1(evdw)
cd              call esc(escloc)
cd              e_sc1=wsc*evdw+wscloc*escloc
cd              print *,'     new',e_sc1,energy(0)
cv              print *,'new       ',energy(0)
cd              call enerprint(energy(0))
          cur_e=e_sc
          do i=2,nres-1
            if (mask_side(i).eq.1) then
              cur_alph(i)=alph(i)
              cur_omeg(i)=omeg(i)
            endif
          enddo
        else
c     ...else revert to the previous one
          e_sc=cur_e
          do i=2,nres-1
            if (mask_side(i).eq.1) then
              alph(i)=cur_alph(i)
              omeg(i)=cur_omeg(i)
            endif
          enddo
        endif
        n_try=n_try+1

      enddo
      n_steps=n_steps+n_try

c     Reset the minimization mask_r to false
      mask_r=.false.

      return
      end

c-------------------------------------------------------------

      subroutine sc_minimize(etot,iretcode,nfun)
c     Minimizes side-chains only, leaving backbone frozen
crc      implicit none

c     Includes
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.VAR'
      include 'COMMON.CHAIN'
      include 'COMMON.FFIELD'

c     Output arguments
      double precision etot
      integer iretcode,nfun

c     Local variables
      integer i
      double precision orig_w(n_ene),energy(0:n_ene)
      double precision var(maxvar)


c     Set non side-chain weights to zero (minimization is faster)
c     NOTE: e(2) does not actually depend on the side-chain, only CA
      orig_w(2)=wscp
      orig_w(3)=welec
      orig_w(4)=wcorr
      orig_w(5)=wcorr5
      orig_w(6)=wcorr6
      orig_w(7)=wel_loc
      orig_w(8)=wturn3
      orig_w(9)=wturn4
      orig_w(10)=wturn6
      orig_w(11)=wang
      orig_w(13)=wtor
      orig_w(14)=wtor_d

      wscp=0.D0
      welec=0.D0
      wcorr=0.D0
      wcorr5=0.D0
      wcorr6=0.D0
      wel_loc=0.D0
      wturn3=0.D0
      wturn4=0.D0
      wturn6=0.D0
      wang=0.D0
      wtor=0.D0
      wtor_d=0.D0

c     Prepare to freeze backbone
      do i=1,nres
        mask_phi(i)=0
        mask_theta(i)=0
        mask_side(i)=1
      enddo

c     Minimize the side-chains
      mask_r=.true.
      call geom_to_var(nvar,var)
      call minimize(etot,var,iretcode,nfun)
      call var_to_geom(nvar,var)
      mask_r=.false.

c     Put the original weights back and calculate the full energy
      wscp=orig_w(2)
      welec=orig_w(3)
      wcorr=orig_w(4)
      wcorr5=orig_w(5)
      wcorr6=orig_w(6)
      wel_loc=orig_w(7)
      wturn3=orig_w(8)
      wturn4=orig_w(9)
      wturn6=orig_w(10)
      wang=orig_w(11)
      wtor=orig_w(13)
      wtor_d=orig_w(14)

      call chainbuild
      call etotal(energy)
      etot=energy(0)

      return
      end

c-------------------------------------------------------------
      subroutine minimize_sc1(etot,x,iretcode,nfun)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2)) 
      include 'COMMON.IOUNITS'
      include 'COMMON.VAR'
      include 'COMMON.GEO'
      include 'COMMON.MINIM'
      common /srutu/ icall
      dimension iv(liv)                                               
      double precision minval,x(maxvar),d(maxvar),v(1:lv),xx(maxvar)
      double precision energia(0:n_ene)
      external func,gradient,fdum
      external func_restr1,grad_restr1
      logical not_done,change,reduce 
      common /przechowalnia/ v

      call deflt(2,iv,liv,lv,v)                                         
* 12 means fresh start, dont call deflt                                 
      iv(1)=12                                                          
* max num of fun calls                                                  
      if (maxfun.eq.0) maxfun=500
      iv(17)=maxfun
* max num of iterations                                                 
      if (maxmin.eq.0) maxmin=1000
      iv(18)=maxmin
* controls output                                                       
      iv(19)=2                                                          
* selects output unit                                                   
c     iv(21)=iout                                                       
      iv(21)=0
* 1 means to print out result                                           
      iv(22)=0                                                          
* 1 means to print out summary stats                                    
      iv(23)=0                                                          
* 1 means to print initial x and d                                      
      iv(24)=0                                                          
* min val for v(radfac) default is 0.1                                  
      v(24)=0.1D0                                                       
* max val for v(radfac) default is 4.0                                  
      v(25)=2.0D0                                                       
c     v(25)=4.0D0                                                       
* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)    
* the sumsl default is 0.1                                              
      v(26)=0.1D0
* false conv if (act fnctn decrease) .lt. v(34)                         
* the sumsl default is 100*machep                                       
      v(34)=v(34)/100.0D0                                               
* absolute convergence                                                  
      if (tolf.eq.0.0D0) tolf=1.0D-4
      v(31)=tolf
* relative convergence                                                  
      if (rtolf.eq.0.0D0) rtolf=1.0D-4
      v(32)=rtolf
* controls initial step size                                            
       v(35)=1.0D-1                                                    
* large vals of d correspond to small components of step                
      do i=1,nphi
        d(i)=1.0D-1
      enddo
      do i=nphi+1,nvar
        d(i)=1.0D-1
      enddo
      IF (mask_r) THEN
       call x2xx(x,xx,nvar_restr)
       call sumsl(nvar_restr,d,xx,func_restr1,grad_restr1,
     &                    iv,liv,lv,v,idum,rdum,fdum)      
       call xx2x(x,xx)
      ELSE
       call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum)      
      ENDIF
      etot=v(10)                                                      
      iretcode=iv(1)
      nfun=iv(6)

      return  
      end  
************************************************************************
      subroutine func_restr1(n,x,nf,f,uiparm,urparm,ufparm)  
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.DERIV'
      include 'COMMON.IOUNITS'
      include 'COMMON.GEO'
      include 'COMMON.FFIELD'
      include 'COMMON.INTERACT'
      include 'COMMON.TIME1'
      common /chuju/ jjj
      double precision energia(0:n_ene),evdw,escloc
      integer jjj
      double precision ufparm,e1,e2
      external ufparm                                                   
      integer uiparm(1)                                                 
      real*8 urparm(1)                                                    
      dimension x(maxvar)
      nfl=nf
      icg=mod(nf,2)+1

#ifdef OSF
c     Intercept NaNs in the coordinates, before calling etotal
      x_sum=0.D0
      do i=1,n
        x_sum=x_sum+x(i)
      enddo
      FOUND_NAN=.false.
      if (x_sum.ne.x_sum) then
        write(iout,*)"   *** func_restr1 : Found NaN in coordinates"
        f=1.0D+73
        FOUND_NAN=.true.
        return
      endif
#endif

      call var_to_geom_restr(n,x)
      call zerograd
      call chainbuild
cd    write (iout,*) 'ETOTAL called from FUNC'
      call egb1(evdw)
      call esc(escloc)
      f=wsc*evdw+wscloc*escloc
cd      call etotal(energia(0))
cd      f=wsc*energia(1)+wscloc*energia(12)
cd      print *,f,evdw,escloc,energia(0)
C
C Sum up the components of the Cartesian gradient.
C
      do i=1,nct
        do j=1,3
          gradx(j,i,icg)=wsc*gvdwx(j,i)
        enddo
      enddo

      return                                                            
      end                                                               
c-------------------------------------------------------
      subroutine grad_restr1(n,x,nf,g,uiparm,urparm,ufparm)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.CHAIN'
      include 'COMMON.DERIV'
      include 'COMMON.VAR'
      include 'COMMON.INTERACT'
      include 'COMMON.FFIELD'
      include 'COMMON.IOUNITS'
      external ufparm
      integer uiparm(1)
      double precision urparm(1)
      dimension x(maxvar),g(maxvar)

      icg=mod(nf,2)+1
      if (nf-nfl+1) 20,30,40
   20 call func_restr1(n,x,nf,f,uiparm,urparm,ufparm)
c     write (iout,*) 'grad 20'
      if (nf.eq.0) return
      goto 40
   30 call var_to_geom_restr(n,x)
      call chainbuild 
C
C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
C
   40 call cartder
C
C Convert the Cartesian gradient into internal-coordinate gradient.
C

      ig=0
      ind=nres-2                                                                    
      do i=2,nres-2                
       IF (mask_phi(i+2).eq.1) THEN                                             
        gphii=0.0D0                                                             
        do j=i+1,nres-1                                                         
          ind=ind+1                                 
          do k=1,3                                                              
            gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)                            
            gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)                           
          enddo                                                                 
        enddo                                                                   
        ig=ig+1
        g(ig)=gphii
       ELSE
        ind=ind+nres-1-i
       ENDIF
      enddo                                        


      ind=0
      do i=1,nres-2
       IF (mask_theta(i+2).eq.1) THEN
        ig=ig+1
	gthetai=0.0D0
	do j=i+1,nres-1
          ind=ind+1
	  do k=1,3
            gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
            gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
          enddo
        enddo
        g(ig)=gthetai
       ELSE
        ind=ind+nres-1-i
       ENDIF
      enddo

      do i=2,nres-1
	if (itype(i).ne.10) then
         IF (mask_side(i).eq.1) THEN
          ig=ig+1
          galphai=0.0D0
	  do k=1,3
	    galphai=galphai+dxds(k,i)*gradx(k,i,icg)
          enddo
          g(ig)=galphai
         ENDIF
        endif
      enddo

      
      do i=2,nres-1
        if (itype(i).ne.10) then
         IF (mask_side(i).eq.1) THEN
          ig=ig+1
	  gomegai=0.0D0
	  do k=1,3
	    gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
          enddo
	  g(ig)=gomegai
         ENDIF
        endif
      enddo

C
C Add the components corresponding to local energy terms.
C

      ig=0
      igall=0
      do i=4,nres
        igall=igall+1
        if (mask_phi(i).eq.1) then
          ig=ig+1
          g(ig)=g(ig)+gloc(igall,icg)
        endif
      enddo

      do i=3,nres
        igall=igall+1
        if (mask_theta(i).eq.1) then
          ig=ig+1
          g(ig)=g(ig)+gloc(igall,icg)
        endif
      enddo
     
      do ij=1,2
      do i=2,nres-1
        if (itype(i).ne.10) then
          igall=igall+1
          if (mask_side(i).eq.1) then
            ig=ig+1
            g(ig)=g(ig)+gloc(igall,icg)
          endif
        endif
      enddo
      enddo

cd      do i=1,ig
cd        write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
cd      enddo
      return
      end
C-----------------------------------------------------------------------------
      subroutine egb1(evdw)
C
C This subroutine calculates the interaction energy of nonbonded side chains
C assuming the Gay-Berne potential of interaction.
C
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.GEO'
      include 'COMMON.VAR'
      include 'COMMON.LOCAL'
      include 'COMMON.CHAIN'
      include 'COMMON.DERIV'
      include 'COMMON.NAMES'
      include 'COMMON.INTERACT'
      include 'COMMON.IOUNITS'
      include 'COMMON.CALC'
      include 'COMMON.CONTROL'
      logical lprn
      evdw=0.0D0
c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
      evdw=0.0D0
      lprn=.false.
c     if (icall.eq.0) lprn=.true.
      ind=0
      do i=iatsc_s,iatsc_e


        itypi=itype(i)
        itypi1=itype(i+1)
        xi=c(1,nres+i)
        yi=c(2,nres+i)
        zi=c(3,nres+i)
        dxi=dc_norm(1,nres+i)
        dyi=dc_norm(2,nres+i)
        dzi=dc_norm(3,nres+i)
        dsci_inv=dsc_inv(itypi)
C
C Calculate SC interaction energy.
C
        do iint=1,nint_gr(i)
          do j=istart(i,iint),iend(i,iint)
          IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
            ind=ind+1
            itypj=itype(j)
            dscj_inv=dsc_inv(itypj)
            sig0ij=sigma(itypi,itypj)
            chi1=chi(itypi,itypj)
            chi2=chi(itypj,itypi)
            chi12=chi1*chi2
            chip1=chip(itypi)
            chip2=chip(itypj)
            chip12=chip1*chip2
            alf1=alp(itypi)
            alf2=alp(itypj)
            alf12=0.5D0*(alf1+alf2)
C For diagnostics only!!!
c           chi1=0.0D0
c           chi2=0.0D0
c           chi12=0.0D0
c           chip1=0.0D0
c           chip2=0.0D0
c           chip12=0.0D0
c           alf1=0.0D0
c           alf2=0.0D0
c           alf12=0.0D0
            xj=c(1,nres+j)-xi
            yj=c(2,nres+j)-yi
            zj=c(3,nres+j)-zi
            dxj=dc_norm(1,nres+j)
            dyj=dc_norm(2,nres+j)
            dzj=dc_norm(3,nres+j)
            rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
            rij=dsqrt(rrij)
C Calculate angle-dependent terms of energy and contributions to their
C derivatives.
            call sc_angular
            sigsq=1.0D0/sigsq
            sig=sig0ij*dsqrt(sigsq)
            rij_shift=1.0D0/rij-sig+sig0ij
C I hate to put IF's in the loops, but here don't have another choice!!!!
            if (rij_shift.le.0.0D0) then
              evdw=1.0D20
cd              write (iout,'(2(a3,i3,2x),17(0pf7.3))')
cd     &        restyp(itypi),i,restyp(itypj),j,
cd     &        rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 
              return
            endif
            sigder=-sig*sigsq
c---------------------------------------------------------------
            rij_shift=1.0D0/rij_shift 
            fac=rij_shift**expon
            e1=fac*fac*aa(itypi,itypj)
            e2=fac*bb(itypi,itypj)
            evdwij=eps1*eps2rt*eps3rt*(e1+e2)
            eps2der=evdwij*eps3rt
            eps3der=evdwij*eps2rt
            evdwij=evdwij*eps2rt*eps3rt
            evdw=evdw+evdwij
            if (lprn) then
            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
cd            write (iout,'(2(a3,i3,2x),17(0pf7.3))')
cd     &        restyp(itypi),i,restyp(itypj),j,
cd     &        epsi,sigm,chi1,chi2,chip1,chip2,
cd     &        eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
cd     &        om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
cd     &        evdwij
            endif

            if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') 
     &                        'evdw',i,j,evdwij

C Calculate gradient components.
            e1=e1*eps1*eps2rt**2*eps3rt**2
            fac=-expon*(e1+evdwij)*rij_shift
            sigder=fac*sigder
            fac=rij*fac
C Calculate the radial part of the gradient
            gg(1)=xj*fac
            gg(2)=yj*fac
            gg(3)=zj*fac
C Calculate angular part of the gradient.
            call sc_grad
          ENDIF
          enddo      ! j
        enddo        ! iint
      enddo          ! i
      end
C-----------------------------------------------------------------------------