anchor.diff:  
     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
  
positive.code:   public :: lj_energy
    public :: lj_force_9_6
    public :: lj_energy_9_6
  
    type lj_params_t
       integer :: n1, n2
  
negative.code: 
    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
  