Coverage for tests/test_potential_functions.py: 100%

1import numpy as np 

2 

3def test_potential_functions() -> None: 

4 import gamdpy as gp 

5 import numba 

6 

7 # note: this example assumes these functions were implemented correctly in version bfa77f6e 

8 assert gp.LJ_12_6(1, [2, 3]) == (5.0, 42.0, 438.0), "Problem with gp.LJ_12_6" 

9 assert gp.LJ_12_6_sigma_epsilon(1, [2, 3]) == (48384.0, 585216.0, 7635456.0), "Problem with gp.LJ_12_6_sigma_epsilon" 

10 # gp.LJ_12_6_params_from_sigma_epsilon_cutoff seems not to be used 

11 #assert gp.LJ_12_6_params_from_sigma_epsilon_cutoff(1, [2, 3, 4]) == (5.0, 42.0, 438.0), "Problem with gp.LJ_12_6_params_from_sigma_epsilon_cutoff" 

12 # consider moving inner functions out for better testing 

13 

14 # Test make_IPL_n 

15 for n, r, a in [(12,1,1), (12,2**(1/6),3), (6,2**(1/6),1), (6,2,4), (1,2,3)]: 

16 ipl_n = gp.make_IPL_n(n) 

17 expected = (a*r**(-n), n*a*r**(-n-2), a*n*(n+1)*r**(-n-2)) 

18 assert np.all(np.isclose(ipl_n(r, (a,)), expected)), f'Problem with make_IPL_n, {(n,r,a)=}' 

19 

20 # Test add_potential_functions 

21 LJ = gp.add_potential_functions(gp.make_IPL_n(12), gp.make_IPL_n(6, first_parameter=1)) 

22 for r, a12, a6 in [(1,1,-1), (2**(1/6),3,-3), (2**(1/6),4,-4), (2,4,4), (2,4,-4)]: 

23 expected = gp.LJ_12_6(r, (a12, a6)) 

24 assert np.all(np.isclose(LJ(r, (a12,a6)), expected, atol=1e-5)), f'Problem with add_potential_functions, {(n,a12,a6)=}' 

25 

26 assert gp.harmonic_bond_function(2.5, [2, 100]) == (12.5, -20.0, 100.0), "Problem with gp.harmonic_bond_function" 

27 # seems correct way: https://stackoverflow.com/questions/624926/how-do-i-detect-whether-a-variable-is-a-function 

28 assert callable(gp.make_IPL_n(12)), "Problem with gp.make_IPL_n" 

29 from sympy.abc import r,s,e 

30 potLJ = 4*e*((s/r)**(12)-(s/r)**6) 

31 potLJ_gp = gp.make_potential_function_from_sympy(potLJ, (s, e)) 

32 assert potLJ_gp(1, (2,3)) == gp.LJ_12_6_sigma_epsilon(1, [2, 3]), "Problem with gp.make_potential_function_from_sympy" 

33 

34 # Test SAAP potential 

35 number_of_params = 8 

36 params = [1.0]*number_of_params 

37 dist = 1.0 

38 pot_SAAP = gp.SAAP(dist, params) 

39 assert len(pot_SAAP) == 3, "Problem with gp.SAAP" 

40 

41 # Test harmonic repulsion, here u=(1-r)² 

42 pair_pot = gp.PairPotential(gp.harmonic_repulsion, params=params, max_num_nbs=128) 

43 params = 2.0, 1.0 

44 dist = 0.5 

45 pot_harm_rep = gp.harmonic_repulsion(dist, params) 

46 assert np.isclose(pot_harm_rep[0],0.25), "Problem with gp.harmonic_repulsion" 

47 assert np.isclose(pot_harm_rep[1],2.0), "Problem with gp.harmonic_repulsion" 

48 assert np.isclose(pot_harm_rep[2],2.0), "Problem with gp.harmonic_repulsion" 

49 eps, sig = 1.43, 1.37 

50 r = 0.98 

51 pot_harm_rep_2 = gp.harmonic_repulsion(r, [eps, sig]) 

52 assert np.isclose(pot_harm_rep_2[0], np.float32(0.5*eps*(1.0-r/sig)**2)), f"Problem with gp.harmonic_repulsion" 

53 du_dr = -eps*(1.0-r/sig)/sig 

54 assert np.isclose(pot_harm_rep_2[1], -du_dr/r), "Problem with gp.harmonic_repulsion" 

55 assert np.isclose(pot_harm_rep_2[2], eps/sig**2), "Problem with gp.harmonic_repulsion" 

56 

57 # Test Hertzian pair potential, u=eps*(1-r/sig)**alpha 

58 params = 1.0, 2.0, 1.0 # Same as "harmonic repulsion" above 

59 dist = 0.5 

60 pot_hertzian = gp.hertzian(dist, params) 

61 assert np.isclose(pot_hertzian[0],0.25), "Problem with gp.hertzian" 

62 assert np.isclose(pot_hertzian[1],2.0), "Problem with gp.hertzian" 

63 assert np.isclose(pot_hertzian[2],2.0), "Problem with gp.hertzian" 

64 eps, alpha, sig = 1.43, 3.1, 1.24 

65 r = 0.98 

66 pot_hertzian_2 = gp.hertzian(r, [eps, alpha, sig]) 

67 assert np.isclose(pot_hertzian_2[0] , eps*(1.0-r/sig)**alpha ), "Problem with gp.hertzian" 

68 assert np.isclose(pot_hertzian_2[1] , alpha*eps*(1.0-r/sig)**(alpha-1)/sig/r ), "Problem with gp.hertzian" 

69 assert np.isclose(pot_hertzian_2[2] , eps*alpha*(alpha-1)*(1.0-r/sig)**(alpha-2)/sig/sig ), "Problem with gp.hertzian" 

70 

71 # needs to add test for apply_shifted_force_cutoff, apply_shifted_potential_cutoff 

72 

73if __name__ == '__main__': # pragma: no cover 

74 test_potential_functions()