anchor.diff:  
     end function lj_energy_9_6
   
  +  pure function lj_energy_9_3(z_sq, epsilon, sigma) result(e)
  +     double precision, intent(in) :: z_sq, epsilon, sigma
  +     double precision :: e
  +
  +     double precision :: sig3_o_r3
  +     
  +     sig3_o_r3 = sigma**3/z_sq**(3.d0/2.d0)
  +
  +     e = 3*sqrt(3.d0)/2*epsilon*(sig3_o_r3**3 - sig3_o_r3) + epsilon
  +  end function lj_energy_9_3
  +
   end module interaction
  
positive.code:     allocate(this% sigma_sq(n1, n2))
      allocate(this% cut(n1, n2))
      allocate(this% cut_sq(n1, n2))
  
      this% epsilon = epsilon
      this% sigma = sigma
      this% sigma_sq = sigma**2
      this% cut = cut
      this% cut_sq = cut**2
  
    end subroutine lj_params_init
  
  
  
negative.code: 
    end function lj_force_9_6
  
    pure function lj_force_9_3(z, z_sq, epsilon, sigma) result(e)
       double precision, intent(in) ::z, z_sq, epsilon, sigma
       double precision :: e
  
       double precision :: sig3_o_r3, force_o_z
       
       sig3_o_r3 = sigma**3/z_sq**(3.d0/2.d0)
  
       force_o_z = 9*sqrt(3.d0)/2*epsilon/z_sq* (3.d0*sig3_o_r3**3 - sig3_o_r3)
       e = force_o_z*z
    end function lj_force_9_3
  
    pure function lj_energy(r_sq, epsilon, sigma) result(e)
      double precision, intent(in) :: r_sq, epsilon, sigma
      double precision :: e
  