Coverage for examples/yukawa.py: 98%

1""" Example of a user defined potential, example of a Yukawa potential. 

2 

3This script demonstrates how to define a user-defined potential function 

4and use it in a gamdpy simulation. The example uses the Yukawa potential 

5as an example, but the same approach can be used for other pair potentials. 

6 

7Comments 

8-------- 

9 

10The pair potential function is JIT compiled using numba.cuda, 

11and it is often the time-consuming part of the simulation. 

12Thus, you may want to optimize this function. For example, 

13it is the experience that numba.float32 is faster than numba.float64. 

14 

15For mathematical functions supported by numba.cuda use the math module, 

16as described in the numba documentation: 

17 

18 https://numba.readthedocs.io/en/stable/cuda/cudapysupported.html#math 

19 

20This example uses a syntax similar to the backend of gamdpy, making it easy to 

21include the code in the package, and making it available to the community. 

22 

23It is recommended ensuring that the analytical derivatives are correct. 

24 

25""" 

26 

27from math import exp # Note math.exp is supported by numba cuda 

28 

29import numpy as np 

30import matplotlib.pyplot as plt 

31import numba 

32 

33import gamdpy as gp 

34 

35 

36def yukawa(dist, params): 

37 """ The Yukawa potential: u(r) = A·exp(-κ·r)/r 

38 

39 parameters: κ, A (κ is the greek letter kappa) 

40 

41 The Yukawa potential is a simple screened Coulomb potential. 

42 The potential is given by: 

43 

44 u(r) = A·exp(-κ·r)/r 

45 

46 where A is the strength of the interaction, 

47 and kappa is the inverse of the screening length. 

48 

49 The s(r) function, used to compute pair forces (𝐅=s·𝐫), is defined as 

50 

51 s(r) = -u'(r)/r 

52 

53 and specifically for the Yukawa potential it is 

54 

55 s(r) = A·exp(-κ·r)·(κ·r + 1)/r³ 

56 

57 The second derivative (`d2u_dr2`) of the potential is given by 

58 

59 u''(r) = A·exp(-κ·r)*([κ·r]² + 2κ·r + 2)/r³ 

60 

61 """ 

62 

63 # Extract parameters 

64 kappa = numba.float32(params[0]) # κ 

65 prefactor = numba.float32(params[1]) # A 

66 

67 # Floats. Note: numba.float32's may make code faster 

68 one = numba.float32(1.0) 

69 two = numba.float32(2.0) 

70 

71 # Compute helper variables 

72 kappa_dist = kappa * dist # κ·r 

73 inv_dist = one / dist # 1/r 

74 inv_dist3 = inv_dist*inv_dist*inv_dist # 1/r³ 

75 exp_kappa_dist = prefactor * exp(-kappa_dist) # A·exp(-κ·r) 

76 

77 # Compute pair potential energy, pair force and pair curvature 

78 

79 # A·exp(-κ·r)/r 

80 u = exp_kappa_dist * inv_dist 

81 

82 # A·exp(-κ·r)·(κ·r + 1)/r³ 

83 s = exp_kappa_dist * (kappa_dist + one) * inv_dist3 

84 

85 # A·exp(-κ·r)*([κ·r]² + 2κ·r + 2)/r³ 

86 d2u_dr2 = exp_kappa_dist * (kappa_dist*kappa_dist + two * kappa_dist + two) * inv_dist3 

87 

88 return u, s, d2u_dr2 # u(r), -u'(r)/r, u''(r) 

89 

90 

91# Plot the Yukawa potential, and confirm the analytical derivatives 

92# are as expected from the numerical derivatives. 

93plt.figure() 

94r = np.linspace(0.8, 3, 200, dtype=np.float32) 

95params = [1.0, 1.0, 2.5] 

96u = [yukawa(rr, params)[0] for rr in r] 

97u_check = params[1] * np.exp(-params[0] * r) / r 

98s = [yukawa(rr, params)[1] for rr in r] 

99s_numerical = -np.gradient(u, r) / r 

100umm = [yukawa(rr, params)[2] for rr in r] 

101umm_numerical = np.gradient(np.gradient(u, r), r) 

102plt.plot(r, u, '-', label='u(r)') 

103plt.plot(r, u_check, '--', label='u(r), check') 

104plt.plot(r, s, '-', label='s(r)') 

105plt.plot(r, s_numerical, '--', label='s(r), numerical') 

106plt.plot(r, umm, label='u\'\'(r)') 

107plt.plot(r, umm_numerical, '--', label='u\'\'(r), numerical') 

108plt.xlabel('r') 

109plt.ylabel('u, s, u\'\'') 

110plt.legend() 

111if __name__ == "__main__": 

112 plt.show() 

113 

114# Setup configuration: FCC Lattice 

115configuration = gp.Configuration(D=3) 

116configuration.make_lattice(gp.unit_cells.FCC, cells=[8, 8, 8], rho=0.973) 

117configuration['m'] = 1.0 

118configuration.randomize_velocities(temperature=0.7) 

119 

120# Setup pair potential: Single component Yukawa system 

121pair_func = gp.apply_shifted_potential_cutoff(yukawa) # Note: We use the above yukawa function here 

122sig, eps, cut = 1.0, 1.0, 2.5 

123pair_pot = gp.PairPotential(pair_func, params=[sig, eps, cut], max_num_nbs=1000) 

124 

125# Setup integrator: NVT 

126integrator = gp.integrators.NVE(dt=0.005) 

127 

128runtime_actions = [gp.MomentumReset(100), 

129 gp.TrajectorySaver(), 

130 gp.ScalarSaver(), ] 

131 

132# Setup Simulation. 

133sim = gp.Simulation(configuration, pair_pot, integrator, runtime_actions, 

134 num_timeblocks=32, steps_per_timeblock=1024, 

135 storage='memory') 

136 

137# Run simulation 

138sim.run()