Coverage for examples/LJchain.py: 98%

1""" Example of a simulation of a Lennard-Jones chain with 10 beads per chain """ 

2 

3import h5py 

4import matplotlib.pyplot as plt 

5import numpy as np 

6import pandas as pd 

7 

8import gamdpy as gp 

9 

10# Generate configuration with a FCC lattice 

11rho = 1.0 

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

13configuration.make_positions(N=2000, rho=rho) 

14configuration['m'] = 1.0 

15configuration.randomize_velocities(temperature=1.44) 

16 

17# Make bonds 

18bond_potential = gp.harmonic_bond_function 

19bond_params = [[1.00, 3000.], ] 

20bond_indices = [[i, i + 1, 0] for i in range(0, configuration.N - 1, 1) if 

21 i % 10 != 9] # 10 bead chains (last molecule: N % 10) 

22bonds = gp.Bonds(bond_potential, bond_params, bond_indices) 

23 

24# Make pair potential 

25pair_func = gp.apply_shifted_force_cutoff(gp.LJ_12_6_sigma_epsilon) 

26sig, eps, cut = 1.0, 1.0, 2.5 

27exclusions = bonds.get_exclusions(configuration) 

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

29 

30# Define molecules (to be implemented) 

31molecules_A = [np.arange(0, 10) + j for j in range(0, configuration.N // 2, 10)] 

32molecules_B = [np.arange(0, 10) + j for j in range(configuration.N // 2, configuration.N, 10)] 

33# molecules = Molecules([molecules_A, molecules_B], names=['H2O', 'Hexaflourobenzene']) 

34# molecules.check_identical() # Maybe make check in __init__ and print warning if not all molecules of a type are identical 

35# molecules.get_numbe_particles() 

36# molecules.get_resulting_force() # Requires atomic forces to be calculated 

37#print(molecules_A[:3]) 

38 

39 

40# Make integrator 

41dt = 0.002 # timestep 

42num_blocks = 16 # Do simulation in this many 'blocks' 

43steps_per_block = 2048 # ... each of this many steps 

44running_time = dt * num_blocks * steps_per_block 

45temperature = 0.7 

46filename = 'Data/LJchain10_Rho1.00_T0.700.h5' 

47Ttarget_function = gp.make_function_ramp(value0=10.000, x0=running_time * (1 / 8), 

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

49integrator0 = gp.integrators.NVT(Ttarget_function, tau=0.2, dt=dt) 

50 

51compute_plan = gp.get_default_compute_plan(configuration) 

52print(compute_plan) 

53 

54runtime_actions = [gp.MomentumReset(100), 

55 gp.TrajectorySaver(), 

56 gp.ScalarSaver(32, {'Fsq':True, 'lapU':True}), ] 

57 

58 

59sim = gp.Simulation(configuration, [pair_pot, bonds], integrator0, runtime_actions, 

60 num_timeblocks=num_blocks, steps_per_timeblock=steps_per_block, 

61 compute_plan=compute_plan, storage=filename) 

62 

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

64for block in sim.run_timeblocks(): 

65 if block % 10 == 0: 

66 print(f'{block=:4} {sim.status(per_particle=True)}') 

67print(sim.summary()) 

68 

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

70sim = gp.Simulation(configuration, [pair_pot, bonds], integrator, runtime_actions, 

71 num_timeblocks=num_blocks, steps_per_timeblock=steps_per_block, 

72 compute_plan=compute_plan, storage=filename) 

73 

74# Setup on-the-fly calculation of Radial Distribution Function 

75calc_rdf = gp.CalculatorRadialDistribution(configuration, bins=1000) 

76 

77print('Production:') 

78for block in sim.run_timeblocks(): 

79 if block % 10 == 0: 

80 print(f'{block=:4} {sim.status(per_particle=True)}') 

81 calc_rdf.update() 

82print(sim.summary()) 

83 

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

85with h5py.File(filename, "r") as f: 

86 data = np.array(gp.extract_scalars(f, columns, first_block=1)) 

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

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

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

90 

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

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

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

94 

95print('gamdpy:') 

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

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

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

99 

100with h5py.File(filename, "r") as f: 

101 dynamics = gp.tools.calc_dynamics(f, first_block=0, qvalues=[7.118]) 

102fig, axs = plt.subplots(3, 1, figsize=(8, 9), sharex=True) 

103fig.subplots_adjust(hspace=0.00) # Remove vertical space between axes 

104axs[0].set_ylabel('MSD') 

105axs[1].set_ylabel('Non Gaussian parameter') 

106axs[2].set_ylabel('Intermediate scattering function') 

107axs[2].set_xlabel('Time') 

108axs[0].grid(linestyle='--', alpha=0.5) 

109axs[1].grid(linestyle='--', alpha=0.5) 

110axs[2].grid(linestyle='--', alpha=0.5) 

111 

112axs[0].loglog(dynamics['times'], dynamics['msd'], 'o--') 

113axs[1].semilogx(dynamics['times'], dynamics['alpha2'], 'o--') 

114axs[2].semilogx(dynamics['times'], dynamics['Fs'], 'o--') 

115if __name__ == "__main__": 

116 plt.show(block=False) 

117 

118rdf = calc_rdf.read() 

119rdf['rdf'] = np.mean(rdf['rdf'], axis=0) 

120fig, axs = plt.subplots(1, 1, figsize=(8, 4)) 

121axs.set_ylabel('RDF') 

122axs.set_xlabel('Distance') 

123axs.grid(linestyle='--', alpha=0.5) 

124axs.plot(rdf['distances'], rdf['rdf'], '-') 

125axs.set_xlim((0.5, 3.5)) 

126if __name__ == "__main__": 

127 plt.show(block=True)