Coverage for examples/minimal_NPT.py: 100%

1""" Minimal example of a constant NpT simulation 

2 

3Simulation of a Lennard-Jones liquid in the NPT ensemble. 

4After equilibration, the simulation runs and calculates the mean potential energy,  

5pressure, density and isothermal compressibility. 

6 

7""" 

8 

9import numpy as np 

10 

11import gamdpy as gp 

12 

13# Setup configuration: FCC Lattice 

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

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

16configuration['m'] = 1.0 

17configuration.randomize_velocities(temperature=2.0) 

18 

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

20pair_func = gp.apply_shifted_potential_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 

24# Setup NPT integrator 

25target_temperature = 2.0 # target temperature for barostat 

26target_pressure = 4.7 # target pressure for barostat 

27integrator = gp.integrators.NPT_Atomic(temperature=target_temperature, 

28 tau=0.4, 

29 pressure=target_pressure, 

30 tau_p=20, 

31 dt=0.001) 

32 

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

34runtime_actions = [gp.TrajectorySaver(), 

35 gp.ScalarSaver(32), 

36 gp.MomentumReset(100)] 

37 

38 

39# NPT Simulation  

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

41 num_timeblocks=16, steps_per_timeblock=2048, 

42 storage='memory') 

43 

44sim.run() # Equilibration run 

45sim.run() # Production run 

46 

47# Thermodynamic properties 

48U, W, K, V = gp.extract_scalars(sim.output, ['U', 'W', 'K', 'Vol'], first_block=1) 

49print(f'Mean U: {np.mean(U)/configuration.N}') 

50print(f'Kinetic temperature (consistency check): {2*np.mean(K)/3/(configuration.N-1)}') 

51print(f'Pressure (consistency check): {(2*np.mean(K)/3+np.mean(W))/np.mean(V)}') 

52rho = configuration.N/np.mean(V) # Average density 

53print(f"Density: {rho}") 

54compressibility = np.var(V)/np.mean(V)/target_temperature 

55print(f'Isothermal compressibility: {compressibility}') 

56print(f'Isothermal bulk modulus: {1/compressibility}')