Coverage for tests/test_LJ_nxyz.py: 95%

1import sys 

2import numpy as np 

3import gamdpy as gp 

4from numba import cuda, config 

5import pandas as pd 

6import pytest 

7 

8from hypothesis import given, strategies as st, settings, Verbosity, example 

9 

10def LJ(nx, ny, nz, rho=0.8442, pb=None, tp=None, skin=None, gridsync=None, UtilizeNIII=None, cut=2.5, integrator='NVE', verbose=True): 

11 

12 # Generate configuration with a FCC lattice 

13 configuration = gp.Configuration(D=3, compute_flags={'Fsq':True, 'lapU':True}) 

14 configuration.make_lattice(gp.unit_cells.FCC, cells=[nx, ny, nz], rho=rho) 

15 configuration['m'] = 1.0 

16 configuration.randomize_velocities(temperature=1.44, seed=0) 

17 assert configuration.N==nx*ny*nz*4, f'Wrong number particles (FCC), {configuration.N} <> {nx*ny*nz*4}' 

18 assert configuration.D==3, f'Wrong dimension (FCC), {configuration.D} <> {3}' 

19 

20 # Allow for overwritting of the default compute_plan 

21 compute_plan = gp.get_default_compute_plan(configuration) 

22 if pb!=None: 

23 compute_plan['pb'] = pb 

24 if tp!=None: 

25 compute_plan['tp'] = tp 

26 if skin!=None: 

27 compute_plan['skin'] = np.float32(skin) 

28 if gridsync!=None: 

29 compute_plan['gridsync'] = gridsync 

30 if UtilizeNIII!=None: 

31 compute_plan['UtilizeNb'] = UtilizeNIII 

32 if verbose: 

33 print('simbox lengths:', configuration.simbox.get_lengths()) 

34 print('compute_plan: ', compute_plan) 

35 

36 # Make the pair potential. 

37 pairfunc = gp.apply_shifted_force_cutoff(gp.LJ_12_6_sigma_epsilon) 

38 sig, eps, cut = 1.0, 1.0, 2.5 

39 pairpot = gp.PairPotential(pairfunc, params=[sig, eps, cut], max_num_nbs=1000) 

40 

41 # Setup the integrator 

42 dt = 0.005 

43 

44 if integrator=='NVE': 

45 integrator = gp.integrators.NVE(dt=dt) 

46 

47 if integrator=='NVT': 

48 integrator = gp.integrators.NVT(temperature=0.70, tau=0.2, dt=dt) 

49 

50 if integrator=='NVT_Langevin': 

51 integrator = gp.integrators.NVT_Langevin(temperature=0.70, alpha=0.1, dt=dt, seed=213) 

52 

53 runtime_actions = [gp.MomentumReset(100), 

54 gp.TrajectorySaver(), 

55 gp.ScalarSaver(8, {'Fsq':True, 'lapU':True, 'stresses':False}), ] 

56 

57 

58 # Setup the Simulation 

59 sim = gp.Simulation(configuration, pairpot, integrator, runtime_actions, 

60 num_timeblocks=2, steps_per_timeblock=1024 * 4, 

61 storage='memory', verbose=False) 

62 

63 # Run simulation one block at a time 

64 for block in sim.run_timeblocks(): 

65 pass 

66 

67 # Make conversion to dataframe a method at some point... 

68 columns = ['U', 'W', 'lapU', 'Fsq', 'K'] 

69 data = np.array(gp.extract_scalars(sim.output, columns, first_block=1)) 

70 df = pd.DataFrame(data.T, columns=columns) 

71 return df 

72 

73def get_results_from_df(df, N, D): 

74 df['E'] = df['U'] + df['K'] # Total energy 

75 df['Tkin'] =2*df['K']/D/(N-1) 

76 df['Tconf'] = df['Fsq']/df['lapU'] 

77 df['dU'] = df['U'] - np.mean(df['U']) 

78 df['dE'] = df['E'] - np.mean(df['E']) 

79 df['dW'] = df['W'] - np.mean(df['W']) 

80 

81 df2 = df.drop(range(len(df)//2)) 

82 

83 df2['dU'] = df2['U'] - np.mean(df2['U']) 

84 df2['dE'] = df2['E'] - np.mean(df2['E']) 

85 df2['dW'] = df2['W'] - np.mean(df2['W']) 

86 

87 var_e = np.var(df['E'])/N 

88 Tkin = np.mean(df2['Tkin']) 

89 Tconf = np.mean(df2['Tconf']) 

90 R = np.dot(df2['dW'], df2['dU'])/(np.dot(df2['dW'], df2['dW'])*np.dot(df2['dU'], df2['dU']))**0.5 

91 Gamma = np.dot(df2['dW'], df2['dU'])/(np.dot(df2['dU'], df2['dU'])) 

92 

93 return var_e, Tkin, Tconf, R, Gamma 

94 

95@pytest.mark.slow 

96@settings(deadline=200_000, max_examples = 8) 

97@given(nx=st.integers(min_value=4, max_value=16), ny=st.integers(min_value=4, max_value=16), nz=st.integers(min_value=4, max_value=16)) 

98@example(nx=4, ny=4, nz=4) 

99@example(nx=16, ny=16, nz=32) 

100@example(nx=5, ny=5, nz=13) 

101def test_nve(nx, ny, nz): 

102 N = nx*ny*nz*4 

103 D = 3 

104 df = LJ(nx, ny, nz, integrator='NVE', cut=2.5, verbose=False) 

105 var_e, Tkin, Tconf, R, Gamma = get_results_from_df(df, N, D) 

106 print(N, '\t', nx, '\t', ny, '\t', nz, '\t', var_e, '\t', Tkin, '\t',Tconf, '\t',R, '\t',Gamma) 

107 assert var_e < 1e-4 

108 assert 0.67 < Tkin < 0.72, print(f'{Tkin=}') 

109 assert 0.67 < Tconf < 0.72, print(f'{Tkin=}') 

110 assert 0.86 < R < 0.99, print(f'{R=}') 

111 assert 5.1 < Gamma < 6.6, print(f'{Gamma=}') 

112 

113 return 

114@pytest.mark.slow 

115@settings(deadline=200_000, max_examples = 8) 

116@given(nx=st.integers(min_value=4, max_value=16), ny=st.integers(min_value=4, max_value=16), nz=st.integers(min_value=4, max_value=16)) 

117@example(nx=4, ny=4, nz=4) 

118@example(nx=16, ny=16, nz=32) 

119@example(nx=10, ny=13, nz=15) 

120def test_nvt(nx, ny, nz): 

121 N = nx*ny*nz*4 

122 D = 3 

123 df = LJ(nx, ny, nz, integrator='NVT', verbose=False) 

124 var_e, Tkin, Tconf, R, Gamma = get_results_from_df(df, N, D) 

125 print(N, '\t', nx, '\t', ny, '\t', nz, '\t', var_e, '\t', Tkin, '\t',Tconf, '\t',R, '\t',Gamma) 

126 # assert var_e < 0.001 

127 assert 0.67 < Tkin < 0.73, print(f'{Tkin=}') 

128 assert 0.67 < Tconf < 0.73, print(f'{Tkin=}') 

129 assert 0.91 < R < 1.00, print(f'{R=}') 

130 assert 4.9 < Gamma < 7.2, print(f'{Gamma=}') 

131 

132 return 

133 

134 

135@pytest.mark.slow 

136@settings(deadline=200_000, max_examples = 8) 

137@given(nx=st.integers(min_value=4, max_value=16), ny=st.integers(min_value=4, max_value=16), nz=st.integers(min_value=4, max_value=16)) 

138@example(nx=4, ny=4, nz=4) 

139@example(nx=16, ny=16, nz=32) 

140def test_nvt_langevin(nx, ny, nz): 

141 N = nx*ny*nz*4 

142 D = 3 

143 df = LJ(nx, ny, nz, integrator='NVT_Langevin', verbose=False) 

144 var_e, Tkin, Tconf, R, Gamma = get_results_from_df(df, N, D) 

145 print(N, '\t', nx, '\t', ny, '\t', nz, '\t', var_e, '\t', Tkin, '\t',Tconf, '\t',R, '\t',Gamma) 

146 # assert var_e < 0.001 

147 assert 0.62 < Tkin < 0.78, print(f'{Tkin=}') 

148 assert 0.62 < Tconf < 0.78, print(f'{Tconf=}') 

149 assert 0.85 < R < 1.00, print(f'{R=}') 

150 assert 4.8 < Gamma < 6.8, print(f'{Gamma=}') 

151 

152 return 

153 

154if __name__ == "__main__": 

155 config.CUDA_LOW_OCCUPANCY_WARNINGS = False 

156 if len(sys.argv)==1 or 'NVE' in sys.argv: 

157 print('Testing LJ NVE:') 

158 test_nve() 

159 print('Passed: LJ NVE!') 

160 if len(sys.argv)==1 or 'NVT' in sys.argv: 

161 print('Testing LJ NVT:') 

162 test_nvt() 

163 print('Passed: LJ NVT!') 

164 if len(sys.argv)==1 or 'NVT_Langevin' in sys.argv: 

165 print('Testing LJ NVT Langevin:') 

166 test_nvt_langevin() 

167 print('Passed: LJ NVT Langevin!') 

168