subroutine energy(element,form,atom,Natoms,theX,x,L,Eatom) 
  use defs; implicit none
  character(len=20),intent(in)::element,form
  integer(i4b),intent(in)::atom,Natoms
  real(dp),dimension(3),intent(in)::theX
  real(dp),dimension(:,:),intent(in)::x
  real(dp),intent(in)::L
  real(dp),intent(out)::Eatom
  integer(i4b)::n,i,j,k
  real(dp)::rr,r,V
  real(dp),dimension(3)::vec
  real(dp),dimension(3,-1:1,-1:1,-1:1)::shift
  ! Make the shift vectors that point to the mirror atoms:
  do k = -1, 1; do j = -1, 1; do i = -1, 1
     shift(1,i,j,k) = i*L
     shift(2,i,j,k) = j*L
     shift(3,i,j,k) = k*L
  end do; end do; end do
  ! Find the energy:  
  Eatom = 0.
  if ( element == 'Ar' .or. element == 'ArGas' ) then
     if ( form == 'hybrid' ) then
        do n = 1, atom-1
           do k = -1, 1; do j = -1, 1; do i = -1, 1
              vec = theX - x(:,n) - shift(:,i,j,k)
              r = sqrt(dot_product(vec,vec))
              Eatom = Eatom + Var(r)
           end do; end do; end do
        end do
        do n = atom+1, Natoms
           do k = -1, 1; do j = -1, 1; do i = -1, 1
              vec = theX - x(:,n) - shift(:,i,j,k)
              r = sqrt(dot_product(vec,vec))
              Eatom = Eatom + Var(r)
           end do; end do; end do
        end do
     else if ( form == 'LJM' ) then
        do n = 1, atom-1
           do k = -1, 1; do j = -1, 1; do i = -1, 1
              vec = theX - x(:,n) - shift(:,i,j,k)
              r = sqrt(dot_product(vec,vec))
              Eatom = Eatom + VLJarM(r)
           end do; end do; end do
        end do
        do n = atom+1, Natoms
           do k = -1, 1; do j = -1, 1; do i = -1, 1
              vec = theX - x(:,n) - shift(:,i,j,k)
              r = sqrt(dot_product(vec,vec))
              Eatom = Eatom + VLJarM(r)
           end do; end do; end do
        end do
     else if ( form == 'LJH' ) then
        do n = 1, atom-1
           do k = -1, 1; do j = -1, 1; do i = -1, 1
              vec = theX - x(:,n) - shift(:,i,j,k)
              r = sqrt(dot_product(vec,vec))
              Eatom = Eatom + VLJarH(r)
           end do; end do; end do
        end do
        do n = atom+1, Natoms
           do k = -1, 1; do j = -1, 1; do i = -1, 1
              vec = theX - x(:,n) - shift(:,i,j,k)
              r = sqrt(dot_product(vec,vec))
              Eatom = Eatom + VLJarH(r)
           end do; end do; end do
        end do
     end if
  else if ( element == 'Kr' .or. element == 'KrGas' ) then
     if ( form == 'hybrid' ) then
        do n = 1, atom-1
           do k = -1, 1; do j = -1, 1; do i = -1, 1
              vec = theX - x(:,n) - shift(:,i,j,k)
              r = sqrt(dot_product(vec,vec))
              Eatom = Eatom + Vkr(r)
           end do; end do; end do
        end do
        do n = atom+1, Natoms
           do k = -1, 1; do j = -1, 1; do i = -1, 1
              vec = theX - x(:,n) - shift(:,i,j,k)
              r = sqrt(dot_product(vec,vec))
              Eatom = Eatom + Vkr(r)
           end do; end do; end do
        end do
     else if ( form == 'LJM' ) then
        do n = 1, atom-1
           do k = -1, 1; do j = -1, 1; do i = -1, 1
              vec = theX - x(:,n) - shift(:,i,j,k)
              r = sqrt(dot_product(vec,vec))
              Eatom = Eatom + VLJkrM(r)
           end do; end do; end do
        end do
        do n = atom+1, Natoms
           do k = -1, 1; do j = -1, 1; do i = -1, 1
              vec = theX - x(:,n) - shift(:,i,j,k)
              r = sqrt(dot_product(vec,vec))
              Eatom = Eatom + VLJkrM(r)
           end do; end do; end do
        end do
     else if ( form == 'LJH' ) then
        do n = 1, atom-1
           do k = -1, 1; do j = -1, 1; do i = -1, 1
              vec = theX - x(:,n) - shift(:,i,j,k)
              r = sqrt(dot_product(vec,vec))
              Eatom = Eatom + VLJkrH(r)
           end do; end do; end do
        end do
        do n = atom+1, Natoms
           do k = -1, 1; do j = -1, 1; do i = -1, 1
              vec = theX - x(:,n) - shift(:,i,j,k)
              r = sqrt(dot_product(vec,vec))
              Eatom = Eatom + VLJkrH(r)
           end do; end do; end do
        end do
     end if
  else
     open(999,file='badEnergyElement')
     write(999,*)'Element',element,' in energy subroutine',&
          'is inappropriate.'
     close(999); stop
  end if
contains
  real function Var(r) ! hybrid Ar
    real(dp)::r,r6
    real(dp),parameter::a=1720.,b=2.6920,c=0.2631
    real(dp),parameter::d=37.943,e=177588.
    r6 = r**6
    Var = a*exp(-b*r)*(1. - c*r) - d/(r6 + e/r6)
  end function Var
  real function Vkr(r) ! hybrid Kr
    real(dp)::r,r6
    real(dp),parameter::a=2499.,b=2.5249,c=0.2466
    real(dp),parameter::d=78.214,e=199064.
    r6 = r**6
    Vkr = a*exp(-b*r)*(1. - c*r) - d/(r6 + e/r6)
  end function Vkr
  real function VLJarM(r) ! LJ Ar, fitted to minimum
    real(dp)::r,r6,r12
    real(dp),parameter::r0 = 3.757,V0 = 0.01234
    real(dp),parameter::r06 = r0**6,r012 = r06*r06
    r6 = r**6; r12 = r6*r6
    VLJarM = V0*( r012/r12 - 2.*r06/r6 )
  end function VLJarM
  real function VLJkrM(r) ! LJ Kr fitted to minimum
    real(dp)::r,r6,r12
    real(dp),parameter::r0 = 4.008,V0 = 0.017338
    real(dp),parameter::r06 = r0**6,r012 = r06*r06
    r6 = r**6; r12 = r6*r6
    VLJkrM = V0*( r012/r12 - 2.*r06/r6 )
  end function VLJkrM
  real function VLJarH(r) ! LJ, Ar Hill's virials
    real(dp)::r,r6,r12
    real(dp),parameter::kB=8.617342D-05 ! eV/K
    real(dp),parameter::r0 = 3.822,V0 = 119.8*kB
    real(dp),parameter::r06 = r0**6,r012 = r06*r06
    r6 = r**6; r12 = r6*r6
    VLJarH = V0*( r012/r12 - 2.*r06/r6 )
  end function VLJarH
  real function VLJkrH(r) ! LJ, Kr  Hill's virials
    real(dp)::r,r6,r12
    real(dp),parameter::kB=8.617342D-05 ! eV/K
    real(dp),parameter::r0 = 4.04,V0 = 171.0*kB
    real(dp),parameter::r06 = r0**6,r012 = r06*r06
    r6 = r**6; r12 = r6*r6
    VLJkrH = V0*( r012/r12 - 2.*r06/r6 )
  end function VLJkrH
end subroutine energy