Coverage for examples/kablj.py: 98%

1""" Example of a binary LJ simulation using gamdpy. 

2 

3NVT simulation of the Kob-Andersen mixture, and compare results with Rumd3 (rumd.org) 

4""" 

5 

6import gamdpy as gp 

7 

8import matplotlib.pyplot as plt 

9import numpy as np 

10import pandas as pd 

11import os.path 

12 

13# Specify statepoint 

14num_part = 2000 

15rho = 1.200 

16temperature = 0.80 

17 

18# Setup configuration:  

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

20configuration.make_positions(N=num_part, rho=rho) 

21configuration['m'] = 1.0 # Specify all masses to unity  

22configuration.randomize_velocities(temperature=2.0) # Initial high temperature for randomizing 

23configuration.ptype[::5] = 1 # Every fifth particle set to type 1 (4:1 mixture) 

24 

25# Setup pair potential: Binary Kob-Andersen LJ mixture. 

26pair_func = gp.apply_shifted_potential_cutoff(gp.LJ_12_6_sigma_epsilon) 

27sig = [[1.00, 0.80], 

28 [0.80, 0.88]] 

29eps = [[1.00, 1.50], 

30 [1.50, 0.50]] 

31cut = np.array(sig)*2.5 

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

33 

34# Setup integrator.  

35# Increase 'num_blocks' for longer runs, better statistics, AND bigger storage consumption 

36# Increase 'steps_per_block' for longer runs 

37dt = 0.004 # timestep 

38num_timeblocks = 32 # Do simulation in this many 'blocks'.  

39steps_per_timeblock = 2*1024 # ... each of this many steps 

40running_time = dt*num_timeblocks*steps_per_timeblock 

41filename = f'Data/KABLJ_Rho{rho:.3f}_T{temperature:.3f}.h5' 

42 

43print('High Temperature followed by cooling and equilibration:') 

44Ttarget_function = gp.make_function_ramp(value0=2.000, x0=running_time*(1/8), 

45 value1=temperature, x1=running_time*(1/4)) 

46integrator = gp.integrators.NVT(temperature=Ttarget_function, tau=0.2, dt=dt) 

47 

48# Setup runtime actions, i.e. actions performed during simulation of timeblocks 

49runtime_actions = [gp.MomentumReset(100),] 

50 

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

52 num_timeblocks=num_timeblocks, steps_per_timeblock=steps_per_timeblock, 

53 storage="memory") 

54 

55for block in sim.run_timeblocks(): 

56 print(f'{sim.status(per_particle=True)}') 

57print(sim.summary()) 

58 

59# Print current status of configuration 

60print(configuration) 

61 

62print('\nProduction:') 

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

64 

65# Setup runtime actions, i.e. actions performed during simulation of timeblocks 

66#runtime_actions = [gp.TrajectorySaver(), 

67# gp.ScalarSaver(16, {'Fsq':True, 'lapU':True}), 

68# gp.MomentumReset(100)] 

69 

70runtime_actions = [gp.MomentumReset(100), 

71 gp.TrajectorySaver(), 

72 gp.ScalarSaver(16, {'Fsq':True, 'lapU':True}), ] 

73 

74 

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

76 num_timeblocks=num_timeblocks, steps_per_timeblock=steps_per_timeblock, 

77 storage=filename) 

78for block in sim.run_timeblocks(): 

79 print(f'{sim.status(per_particle=True)}') 

80print(sim.summary()) 

81 

82# Print current status of configuration 

83print(configuration) 

84 

85 

86 

87columns = ['U', 'W', 'K', 'Fsq', 'lapU', 'Vol'] 

88data = np.array(gp.extract_scalars(sim.output, columns, first_block=0)) 

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

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

91df['t'] = np.arange(len(df['U'])) * dt * sim.output['scalar_saver'].attrs["steps_between_output"] 

92 

93mu = np.mean(df['U'])/configuration.N 

94mw = np.mean(df['W'])/configuration.N 

95cvex = np.var(df['U'])/temperature**2/configuration.N 

96 

97print('\ngamdpy:') 

98print(f'Potential energy: {mu:.4f}') 

99print(f'Excess heat capacity: {cvex:.3f}') 

100print(f'Virial {mw:.4f}') 

101 

102if rho==1.200 and temperature==0.800: 

103 mu3 = -6.346 

104 mw3 = 5.534 

105 cvex3 = 0.0001089086505*10000/0.8**2 

106 print('\nRumd3:') 

107 print(f'Potential energy: {mu3:.4f}') 

108 print(f'Excess heat capacity: {cvex3:.3f}') 

109 print(f'Virial {mw3:.4f}') 

110 

111if __name__ == "__main__": 

112 gp.plot_scalars(df, configuration.N, configuration.D, figsize=(10,8), block=False) 

113 

114dyn = gp.tools.calc_dynamics(sim.output, first_block=0, qvalues=[7.5, 5.5]) 

115fig, axs = plt.subplots(1, 1, figsize=(6,4)) 

116axs.loglog(dyn['times'], dyn['msd'], '.-', label=['A (gamdpy)', 'B (gamdpy)']) 

117axs.set_xlabel('Time') 

118axs.set_ylabel('MSD') 

119 

120rumd3filename = f'Data/KABLJ_msd_R{rho:.3f}_T{temperature:.3f}_rumd3.dat' 

121if os.path.isfile(rumd3filename): 

122 msd3 = np.loadtxt(rumd3filename) 

123 axs.loglog(msd3[:21,0], msd3[:21,1:], '--', label=['A (rumd3)', 'B (rumd3)']) 

124 

125axs.legend() 

126if __name__ == "__main__": 

127 plt.show(block=True) 

128 

129if rho==1.200 and temperature==0.800: 

130 print('\nTesting complience with Rumd3:') 

131 assert abs(mu - mu3) < 0.01, f"{mu=} but in rumd3 is {mu3=}" 

132 assert abs(cvex - cvex3) < 0.2 , f"{cvex=} but in rumd3 is {cvex3=}" 

133 assert abs(mw - mw3) < 0.03, f"{mw=} but in rumd3 is {mw3=}" 

134 print('Passed')