Coverage for examples/isomorph.py: 100%

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.py (generates: isomorph_dynamics.pdf) 

9python plot_isomorph_rdf.py (generates: isomorph_rdf.pdf) 

10 

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

12 

13""" 

14 

15import pickle 

16 

17import numpy as np 

18 

19import gamdpy as gp 

20 

21# Setup pair potential. 

22pair_func = gp.apply_shifted_force_cutoff(gp.LJ_12_6_sigma_epsilon) 

23sig, eps, cut = 1.0, 1.0, 2.5 

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

25 

26T = 2.00 

27rhos = [1.00, 1.05, 1.10, 1.15, 1.20, 1.20] 

28data = [] 

29 

30for index, rho in enumerate(rhos): 

31 print(f'\nRho = {rho}, Temperature = {T}') 

32 

33 # Setup fcc configuration 

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

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

36 configuration['m'] = 1.0 

37 configuration.randomize_velocities(temperature=2 * T) 

38 

39 # Setup integrator 

40 integrator = gp.integrators.NVT(temperature=T, tau=0.2, dt=0.0025) 

41 

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

43 runtime_actions = [gp.TrajectorySaver(), 

44 gp.ScalarSaver(16, {'W':True}), 

45 gp.MomentumReset(100)] 

46 

47 # Setup Simulation 

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

49 num_timeblocks=16, # try something like 128 for better statistics, 

50 steps_per_timeblock=512, 

51 storage='memory') 

52 

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

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

55 

56 print('Equilibration:', end='\t') 

57 for block in sim.run_timeblocks(): 

58 pass 

59 print(sim.status(per_particle=True)) 

60 

61 print('Production:', end='\t') 

62 for block in sim.run_timeblocks(): 

63 calc_rdf.update() 

64 print(sim.status(per_particle=True)) 

65 

66 # Do data analysis 

67 U, W = gp.extract_scalars(sim.output, ['U', 'W'], first_block=1) 

68 dU = U - np.mean(U) 

69 dW = W - np.mean(W) 

70 gamma = np.dot(dW,dU)/np.dot(dU,dU) 

71 R = np.dot(dW,dU)/(np.dot(dW,dW)*np.dot(dU,dU))**0.5 

72 print(f'Gamma = {gamma:.3f}, R = {R:.3f}') 

73 

74 dynamics = gp.tools.calc_dynamics(sim.output, 0, qvalues=7.5*rho**(1/3)) 

75 rdf = calc_rdf.read() 

76 data.append({'rho':rho, 'T':float(T), 'dynamics':dynamics, 'rdf':rdf}) 

77 

78 # Set temperature for next simulation 

79 if index+1<len(rhos)-1: 

80 T = round((rhos[index+1]/rho)**gamma*T, 3) # Isomorph theory 

81 else: 

82 T = 2.00 # Last simulation Isothermal to first simulation 

83 

84with open('Data/isomorph.pkl', 'wb') as f: 

85 pickle.dump(data, f) 

86 

87# To generat plots (isomorph_dynamics.pdf & isomorph_rdf.pdf) of data:  

88# python plot_isomorph_dynamics.py; python plot_isomorph_rdf.py