Coverage for tests/test_nvt_langevin.py: 68%

1def test_step_langevin(verbose=False, plot_figures=False) -> None: 

2 """ Test NVT langevin thermostat 

3 Test temperature T=1.2 (r_c=2.5) fcc-liquid coexistence state-point in https://doi.org/10.1063/1.4818747 """ 

4 import pandas as pd 

5 import numpy as np 

6 import matplotlib.pyplot as plt 

7 import time 

8 

9 

10 import gamdpy as gp 

11 

12 # State-point 

13 temperature = 1.2 

14 density = 1 / 0.9672 

15 

16 # Expected values 

17 expected_kinetic_energy = 3 / 2 * temperature 

18 expected_total_energy = -4.020 

19 expected_potential_energy = expected_total_energy - expected_kinetic_energy 

20 

21 # Setup configuration (give temperature kick to particles to get closer to equilibrium) 

22 configuration = gp.Configuration(D=3, compute_flags={'W':True, 'K':True}) 

23 configuration.make_lattice(gp.unit_cells.FCC, cells=[7, 7, 7], rho=density) 

24 configuration['m'] = 1.0 

25 configuration.randomize_velocities(temperature=2 * temperature, seed=0) 

26 

27 # Setup pair potential. 

28 pairfunc = gp.apply_shifted_potential_cutoff(gp.LJ_12_6_sigma_epsilon) 

29 sig, eps, cut = 1.0, 1.0, 2.5 

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

31 

32 # Setup integrator 

33 dt = 0.005 

34 alpha = 0.1 

35 integrator = gp.integrators.NVT_Langevin(temperature, alpha=alpha, dt=dt, seed=0) 

36 

37 runtime_actions = [gp.MomentumReset(100), 

38 gp.ScalarSaver(32, {'W':True, 'K':True}), ] 

39 

40 # Setup the Simulation 

41 num_blocks = 32 

42 steps_per_block = 512 

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

44 num_timeblocks=num_blocks, steps_per_timeblock=steps_per_block, 

45 storage='memory') 

46 

47 # Run simulation one block at a time 

48 for block in sim.run_timeblocks(): 

49 pass 

50 print(sim.summary()) 

51 

52 # Convert scalars to dataframe 

53 columns = ['U', 'W', 'K'] 

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

55 df = pd.DataFrame(data.T/configuration.N, columns=columns) 

56 

57 # Compute summary statistics 

58 summary_statistics = df.describe() 

59 potential_energy = summary_statistics['U']['mean'] # per particle 

60 kinetic_energy = summary_statistics['K']['mean'] # per particle 

61 if verbose: 

62 print( 

63 f'Potential energy (per particle): {potential_energy: 8.4f} (expected: {expected_potential_energy: 8.4f})') 

64 print(f'Kinetic energy (per particle): {kinetic_energy: 8.4f} (expected: {expected_kinetic_energy: 8.4f})') 

65 

66 # Assert that the energies are close to the expected values 

67 assert abs(potential_energy - expected_potential_energy) < 0.05 

68 assert abs(kinetic_energy - expected_kinetic_energy) < 0.05 

69 

70 if plot_figures: 

71 plt.figure(figsize=(6, 8)) 

72 plt.subplot(2, 1, 1) 

73 plt.plot(df['U'] , label='u') 

74 plt.subplot(2, 1, 2) 

75 plt.plot(df['K'] , label='k') 

76 df = df.iloc[len(df) // 2:, :] # last half 

77 if plot_figures: 

78 plt.subplot(2, 1, 1) 

79 plt.plot(df['U'] , label='u (last half)') 

80 plt.plot([0, 2*len(df)], [expected_potential_energy, expected_potential_energy], 'k--', label='expected') 

81 plt.ylabel(r'Potential energy, $u$') 

82 plt.legend() 

83 plt.subplot(2, 1, 2) 

84 plt.plot(df['K'] , label='k (last half)') 

85 plt.plot([0, 2*len(df)], [expected_kinetic_energy, expected_kinetic_energy], 'k--', label='expected') 

86 plt.ylabel(r'Kinetic energy, $k$') 

87 plt.xlabel('Outer loop step') 

88 plt.legend() 

89 plt.show() 

90 

91 return 

92 

93if __name__ == '__main__': 

94 test_step_langevin(verbose=True, plot_figures=True)