Coverage for examples/shear_SLLOD.py: 97%

1""" Example of a Simulation using gamdpy, using explicit blocks. 

2 

3Simulation of a Lennard-Jones crystal in the NVT ensemble followed by shearing with SLLOD  

4and Lees-Edwards boundary conditions. Runs one shear rate but easy to make a loop over shear rates. 

5 

6""" 

7import os 

8import h5py 

9import numpy as np 

10import gamdpy as gp 

11import matplotlib.pyplot as plt 

12 

13run_NVT = True # False #  

14 

15# Setup pair potential: Single component 12-6 Lennard-Jones 

16pairfunc = gp.apply_shifted_force_cutoff(gp.LJ_12_6_sigma_epsilon) 

17sig, eps, cut = 1.0, 1.0, 2.5 

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

19 

20temperature_low = 0.700 

21gridsync = True 

22 

23 

24if run_NVT: 

25 # Setup configuration: FCC Lattice 

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

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

28 configuration['m'] = 1.0 

29 configuration.randomize_velocities(temperature=2.0) 

30 

31 

32 # Setup integrator to melt the crystal 

33 dt = 0.005 

34 num_blocks = 10 

35 steps_per_block = 2048 

36 running_time = dt*num_blocks*steps_per_block 

37 

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

39 value1=temperature_low, x1=running_time*(7/8)) 

40 integrator_NVT = gp.integrators.NVT(Ttarget_function, tau=0.2, dt=dt) 

41 

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

43 runtime_actions = [gp.TrajectorySaver(), 

44 gp.ScalarSaver(), 

45 gp.MomentumReset(100)] 

46 

47 

48 

49 # Set simulation up. Total number of timesteps: num_blocks * steps_per_block 

50 sim_NVT = gp.Simulation(configuration, pairpot, integrator_NVT, runtime_actions, 

51 num_timeblocks=num_blocks, steps_per_timeblock=steps_per_block, 

52 storage='memory') 

53 

54 

55 

56 for block in sim_NVT.run_timeblocks(): 

57 print(block) 

58 print(sim_NVT.status(per_particle=True)) 

59 

60 # save both in hdf5 and rumd-3 formats 

61 with h5py.File('LJ_cooled_0.70.h5', 'w') as fout: 

62 configuration.save(fout, "configuration") 

63 

64else: 

65 with h5py.File('LJ_cooled_0.70.h5', 'r') as fin: 

66 configuration = Configuration.from_h5(fin, "configuration", compute_flags={'stresses':True}) 

67 

68compute_plan = gp.get_default_compute_plan(configuration) 

69compute_plan['gridsync'] = gridsync 

70print('compute_plan') 

71print(compute_plan) 

72print("Now run SLLOD simulation on what should now be a glass or polycrystal") 

73 

74sc_output = 8 

75 

76 

77dt = 0.01 

78sr = 0.02 # restuls for different values shown in comments below. This value only takes 4 seconds to run so good for running as a test 

79 

80configuration.simbox = gp.LeesEdwards(configuration.D, configuration.simbox.get_lengths()) 

81 

82integrator_SLLOD = gp.integrators.SLLOD(shear_rate=sr, dt=dt) 

83 

84# set the kinetic temperature to the exact value associated with the desired 

85# temperature since SLLOD uses an isokinetic thermostat 

86configuration.set_kinetic_temperature(temperature_low, ndofs=configuration.N*3-4) # remove one DOF due to constraint on total KE 

87 

88# Setup Simulation. Total number of timesteps: num_blocks * steps_per_block 

89totalStrain = 10.0 

90steps_per_block = 4096 

91total_steps = int(totalStrain / (sr*dt)) + 1 

92num_blocks = total_steps // steps_per_block + 1 

93strain_transient = 1.0 # how much of the output to ignore 

94num_steps_transient = int(strain_transient / (sr*dt) ) + 1 

95 

96 

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

98runtime_actions = [gp.TrajectorySaver(include_simbox=True), 

99 gp.MomentumReset(100), 

100 gp.StressSaver(sc_output, compute_flags={'stresses':True}), 

101 gp.ScalarSaver(sc_output)] 

102 

103 

104print(f'num_blocks={num_blocks}') 

105sim_SLLOD = gp.Simulation(configuration, pairpot, integrator_SLLOD, runtime_actions, 

106 num_timeblocks=num_blocks, steps_per_timeblock=steps_per_block, 

107 storage='memory', compute_plan=compute_plan) 

108 

109# Run simulation one block at a time 

110for block in sim_SLLOD.run_timeblocks(): 

111 print(sim_SLLOD.status(per_particle=True)) 

112 configuration.simbox.copy_to_host() 

113 box_shift = configuration.simbox.box_shift 

114 lengths = configuration.simbox.get_lengths() 

115 print(f'box-shift={box_shift:.4f}, strain = {box_shift/lengths[1]:.4f}') 

116print(sim_SLLOD.summary()) 

117 

118 

119U, K, W, V_sxy = gp.extract_scalars(sim_SLLOD.output, ['U', 'K', 'W', 'Sxy']) 

120N = configuration.N 

121u, k, sxy = U/N,K/N, V_sxy / configuration.get_volume() 

122 

123# alternative (newer way) to get the shear stress 

124full_stress_tensor = gp.StressSaver.extract(sim_SLLOD.output) 

125sxy_alt = full_stress_tensor[:,0,1] 

126 

127times = np.arange(len(u)) * sc_output * dt 

128stacked_output = np.column_stack((times, u, k, sxy, sxy_alt)) 

129np.savetxt('shear_run.txt', stacked_output, delimiter=' ', fmt='%f') 

130 

131 

132 

133strains = times * sr 

134 

135num_items_transient = num_steps_transient // sc_output 

136print(f'num_items_transient={num_items_transient}') 

137sxy_SS = sxy[num_items_transient:] 

138 

139sxy_mean = np.mean(sxy_SS) 

140print(f'{sr:.2g} {sxy_mean:.6f}') 

141 

142#plt.figure(1) 

143#plt.plot(strains, k) 

144#plt.plot(time, u) 

145#plt.figure(2) 

146#plt.plot(strains, sxy) 

147#plt.show() 

148 

149# STRAINRATE VS MEAN STRESS 

150 

151 

152# 0.001 0.014687 

153# 0.0025 0.037265 

154# 0.005 0.071615 

155# 0.0075 0.111890 

156# 0.01 0.140723 

157# 0.02 0.264056 

158# 0.03 0.363014 

159 

160 

161# quadratic fit gives the following 

162# -0.00083871 + 15.436 * x - 110.19*x^2 

163 

164# The small value of the stress at zero is consistent with zero 

165# which is promising and we can read the Newtonian viscosity off as 15.436