Coverage for examples/D4.py: 100%

1import matplotlib.pyplot as plt 

2import numpy as np 

3import gamdpy as gp 

4import math 

5 

6temperature, density, npart, D = 3.40410425853, 0.34, 2401, 4 

7 

8configuration = gp.Configuration(D=D) 

9# Setup configuration. SC Lattice in 4D 

10configuration.make_positions(N=npart, rho=density) 

11# Setup masses and velocities 

12configuration['m'] = 1.0 # Set all masses to 1.0 

13configuration.randomize_velocities(temperature=temperature) 

14 

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

16pair_func = gp.apply_shifted_force_cutoff(gp.LJ_12_6_sigma_epsilon) 

17pair_func = gp.apply_shifted_potential_cutoff(gp.LJ_12_6_sigma_epsilon) 

18sig, eps, cut = 1.0, 1.0, 2.5 

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

20 

21# Setup integrator 

22integrator = gp.integrators.NVT(temperature=temperature, tau=0.08, dt=0.001) 

23 

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

25runtime_actions = [gp.TrajectorySaver(), 

26 gp.ScalarSaver(), 

27 gp.MomentumReset(100)] 

28 

29 

30# Setup Simulation 

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

32 num_timeblocks=16, steps_per_timeblock=4096, 

33 storage='memory') 

34 

35# Run simulation 

36sim.run(verbose=False) 

37sim.run(verbose=False) 

38 

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

40dU = U - np.mean(U) 

41dW = W - np.mean(W) 

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

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

44print(f'{density=} {temperature=} {R=:.3f} {gamma=:.3f}') 

45# Reference results are from http://glass.ruc.dk/pdf/articles/2016_JChemPhys_144_231101.pdf 

46assert math.isclose(0.80, R , rel_tol=0.1), f"{R=} but should be 0.80" 

47assert math.isclose(3.33, gamma, rel_tol=0.1), f"{gamma=} but should be 3.33" 

48