Coverage for examples/LJ.py: 84%

1import sys 

2 

3import numba 

4import numpy as np 

5import pandas as pd 

6 

7import gamdpy as gp 

8 

9integrator_name = 'NVE' 

10if 'NVT' in sys.argv: 

11 integrator_name = 'NVT' 

12if 'NVT_Langevin' in sys.argv: 

13 integrator_name = 'NVT_Langevin' 

14if 'NPT_Langevin' in sys.argv: # use with NoRDF since box size is varying 

15 integrator_name = 'NPT_Langevin' 

16 

17# Generate configuration with a FCC lattice 

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

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

20configuration['m'] = 1.0 

21configuration.randomize_velocities(temperature=1.44) 

22 

23# Make pair potential 

24pair_func = gp.apply_shifted_force_cutoff(gp.LJ_12_6_sigma_epsilon) 

25sig, eps, cut = 1.0, 1.0, 2.5 

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

27 

28# Make integrator 

29dt = 0.005 # timestep 

30num_blocks = 64 # Do simulation in this many 'blocks' 

31steps_per_block = 2*1024 # ... each of this many steps 

32running_time = dt*num_blocks*steps_per_block 

33temperature = 0.7 # Not used for NVE 

34pressure = 1.2 # Not used for NV* 

35 

36# Parameters, temperature and pressure, can be functions of time: 

37#temperature = gp.make_function_ramp(value0=0.7, x0=0.5*running_time, value1=1.7, x1=0.8*running_time) 

38#pressure = gp.make_function_ramp(value0=1.2, x0=0.3*running_time, value1=3.2, x1=0.6*running_time) 

39 

40if integrator_name=='NVE': 

41 integrator = gp.integrators.NVE(dt=dt) 

42if integrator_name=='NVT': 

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

44if integrator_name=='NVT_Langevin': 

45 integrator = gp.integrators.NVT_Langevin(temperature=temperature, alpha=0.2, dt=dt, seed=2023) 

46if integrator_name=='NPT_Langevin': 

47 integrator = gp.integrators.NPT_Langevin(temperature=temperature, pressure=pressure, 

48 alpha=0.1, alpha_baro=0.0001, mass_baro=0.0001, 

49 volume_velocity=0.0, barostatModeISO = True , boxFlucCoord = 2, 

50 dt=dt, seed=2023) 

51 

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

53compute_plan = gp.get_default_compute_plan(configuration) 

54print(compute_plan) 

55 

56runtime_actions = [gp.MomentumReset(100), 

57 gp.TrajectorySaver(), 

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

59 

60 

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

62 num_timeblocks=num_blocks, steps_per_timeblock=steps_per_block, 

63 compute_plan=compute_plan, storage='memory') 

64 

65for block in sim.run_timeblocks(): 

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

67print(sim.summary()) 

68 

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

70data = np.array(gp.extract_scalars(sim.output, columns, first_block=1)) 

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

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

73if integrator_name!='NVE' and callable(temperature): 

74 df['Ttarget'] = numba.vectorize(temperature)(np.array(df['t'])) 

75if integrator_name=='NPT_Langevin' and callable(pressure): 

76 df['Ptarget'] = numba.vectorize(pressure)(np.array(df['t'])) 

77 

78gp.plot_scalars(df, configuration.N, configuration.D, figsize=(10,8), block=True)