Coverage for examples/switching_integrator.py: 98%

1""" Example of performing several simulation in one go using gamdpy. 

2 

3An isomorph is traced out using the gamma method. The script demomstrates 

4the possibility of keeping the output of the simulation in memory (storage='memory'). 

5This is usefull when a lot of short simulations are performed. 

6 

7To plot the results do:  

8python plot_isomorph_dynamics.pdf 

9python plot_isomorph_rdf.pdf 

10 

11For a simpler script performing multiple simulations, see isochore.py 

12 

13""" 

14 

15import matplotlib.pyplot as plt 

16 

17import gamdpy as gp 

18 

19# Setup pair potential. 

20pair_func = gp.apply_shifted_force_cutoff(gp.LJ_12_6_sigma_epsilon) 

21sig, eps, cut = 1.0, 1.0, 2.5 

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

23 

24T = 0.8 

25 

26# Setup fcc configuration 

27configuration = gp.Configuration(D=3) 

28configuration.make_lattice(gp.unit_cells.FCC, cells=[8, 8, 8], rho=0.84) 

29configuration['m'] = 1.0 

30configuration.randomize_velocities(temperature=0.7) 

31 

32# Setup integrators 

33integrator1 = gp.integrators.NVT(temperature=T, tau=0.2, dt=0.0025) 

34integrator2 = gp.integrators.NVT(temperature=T, tau=0.2, dt=0.0025) 

35 

36runtime_actions = [gp.TrajectorySaver(), 

37 gp.ScalarSaver(1), 

38 gp.MomentumReset(100)] 

39 

40# Setup Simulations 

41sim1 = gp.Simulation(configuration, pair_pot, integrator1, runtime_actions, 

42 num_timeblocks=4, steps_per_timeblock=512, 

43 storage='memory') 

44 

45print(configuration['r'][1]) 

46print('Integrator1, Equilibration:', end='\t') 

47for block in sim1.run_timeblocks(): 

48 pass 

49print(sim1.status(per_particle=True)) 

50U1, K1 = gp.extract_scalars(sim1.output, ['U', 'K'], first_block=0) 

51E1 = U1 + K1 

52 

53print('Integrator1, Production:', end='\t') 

54for block in sim1.run_timeblocks(): 

55 pass 

56print(sim1.status(per_particle=True)) 

57U2, K2 = gp.extract_scalars(sim1.output, ['U', 'K'], first_block=0) 

58E2 = U2 + K2 

59 

60sim2 = gp.Simulation(configuration, pair_pot, integrator2, runtime_actions, 

61 num_timeblocks=4, steps_per_timeblock=512, 

62 storage='memory') 

63 

64print(configuration['r'][1]) 

65 

66print('Integrator2, Production:', end='\t') 

67for block in sim2.run_timeblocks(): 

68 pass 

69print(sim2.status(per_particle=True)) 

70 

71U3, K3 = gp.extract_scalars(sim2.output, ['U', 'K'], first_block=0) 

72E3 = U3 + K3 

73 

74 

75plt.plot(U1, '.-', label='Integrator1, Equilibration') 

76plt.plot(U2, '.-', label='Integrator1, Production') 

77plt.plot(U3, '.-', label='Integrator2, Production') 

78plt.legend() 

79if __name__ == "__main__": 

80 plt.show()