sarkas.utilities.maths.force_error_analytic_lcl#

sarkas.utilities.maths.force_error_analytic_lcl(potential_type, cutoff_length, potential_matrix, rescaling_const)[source]#

Calculate the force error from its analytic formula.

Parameters

  • potential_type (str) – Type of potential used.

  • potential_matrix (numpy.ndarray) – Potential parameters.

  • cutoff_length (float) – Cutoff radius of the potential.

  • rescaling_const (float) – Constant by which to rescale the force error.

    In case of electric forces = \(Q^2/(4 \pi \epsilon_0) 1/a^2\).

Returns

force_error (float) – Force error in units of rescaling_const.

Examples

Yukawa, QSP, EGS force errors are calculated the same way.

>>> import numpy as np
>>> # Look more about potential matrix
>>> potential_matrix = np.zeros((2, 2, 2))
>>> kappa = 2.0 # in units of a_ws
>>> potential_matrix[1,:,:] = kappa
>>> rc = 6.0 # in units of a_ws
>>> const = 1.0 # Rescaling const
>>> force_error_analytic_lcl("yukawa", rc, potential_matrix, const)
2.1780665692875655e-05

Lennard jones potential

>>> import numpy as np
>>> potential_matrix = np.zeros( (5,2,2))
>>> sigma = 3.4e-10
>>> pot_const = 4.0 * 1.656e-21    # 4*epsilon
>>> high_pow, low_pow = 12, 6
>>> potential_matrix[0] = pot_const
>>> potential_matrix[1] = sigma
>>> potential_matrix[2] = high_pow
>>> potential_matrix[3] = low_pow
>>> rc = 10 * sigma
>>> force_error_analytic_lcl("lj", rc, potential_matrix, 1.0)
1.4590050212983888e-16

Moliere potential

>>> import numpy as np
>>> from scipy.constants import epsilon_0, pi, elementary_charge
>>> charge = 4.0 * elementary_charge  # = 4e [C] mks units
>>> coul_const = 1.0 / (4.0 * pi * epsilon_0)
>>> screening_charges = np.array([0.5, -0.5, 1.0])
>>> screening_lengths = np.array([5.99988000e-11, 1.47732309e-11, 1.47732309e-11])  # [m]
>>> params_len = len(screening_lengths)
>>> pot_mat = np.zeros((2 * params_len + 1, 2, 2))
>>> pot_mat[0] = coul_const * charge ** 2
>>> pot_mat[1: params_len + 1] = screening_charges.reshape((3, 1, 1))
>>> pot_mat[params_len + 1:] = 1. / screening_lengths.reshape((3, 1, 1))
>>> rc = 6.629e-10
>>> force_error_analytic_lcl("moliere", rc, pot_mat, 1.0)
2.1223648580087958e-14