Coverage for tests/test_SLLOD.py: 95%

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 

5 

6""" 

7 

8 

9 

10def test_SLLOD(run_NVT=False): 

11 from pathlib import Path 

12 

13 import numpy as np 

14 import matplotlib.pyplot as plt 

15 

16 import gamdpy as gp 

17 

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

19 pairfunc = gp.apply_shifted_force_cutoff(gp.LJ_12_6_sigma_epsilon) 

20 sig, eps, cut = 1.0, 1.0, 2.5 

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

22 

23 temperature = 0.700 

24 gridsync = True 

25 

26 # read reference configuration 

27 configuration = None 

28 possible_file_paths = ['reference_data/conf_LJ_N2048_rho0.973_T0.700.h5', 'tests/reference_data/conf_LJ_N2048_rho0.973_T0.700.h5'] 

29 for path in possible_file_paths: 

30 if Path(path).is_file(): 

31 configuration = gp.configuration_from_hdf5(path, compute_flags={'stresses':True}) 

32 break 

33 if configuration is None: 

34 raise FileNotFoundError(f'Could not find configuration file in {possible_file_paths}') 

35 

36 compute_plan = gp.get_default_compute_plan(configuration) 

37 compute_plan['gridsync'] = gridsync 

38 sc_output = 1 

39 sr = 0.1 

40 dt = 0.01 

41 

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

43 

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

45 

46 # Test get_kernel 

47 integrator_SLLOD.get_kernel(configuration=configuration, 

48 compute_plan = gp.get_default_compute_plan(configuration), 

49 compute_flags = gp.get_default_compute_flags(), 

50 interactions_kernel=None, 

51 verbose=True) 

52 

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

54 # temperature since SLLOD uses an isokinetic thermostat 

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

56 

57 runtime_actions = [gp.MomentumReset(100), 

58 gp.TrajectorySaver(include_simbox=True), 

59 gp.StressSaver(sc_output), 

60 gp.ScalarSaver(sc_output, {'stresses':True}), ] 

61 

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

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

64 num_timeblocks=3, steps_per_timeblock=128, 

65 storage='memory', compute_plan=compute_plan) 

66 

67 # Run simulation one block at a time 

68 for block in sim_SLLOD.run_timeblocks(): 

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

70 configuration.simbox.copy_to_host() 

71 box_shift = configuration.simbox.box_shift 

72 lengths = configuration.simbox.get_lengths() 

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

74 print(sim_SLLOD.summary()) 

75 

76 sxy = gp.StressSaver.extract(sim_SLLOD.output)[:,0,1] 

77 sxy_mean = np.mean(sxy) 

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

79 assert (np.isclose(sxy_mean, 2.71, atol=0.005 )) 

80 assert(np.isclose(pairpot.nblist.d_nbflag[2], 49, atol=1)) 

81 

82 #sxy_sc = gp.extract_scalars(sim_SLLOD.output, ['Sxy'])/configuration.get_volume() 

83 #sxy_mean_sc = np.mean(sxy) 

84 #assert (np.isclose(sxy_mean_sc, 2.71, atol=0.005 )) 

85 

86 

87 

88if __name__ == '__main__': 

89 test_SLLOD()