Coverage for examples/D2.py: 96%

1import matplotlib.pyplot as plt 

2 

3import gamdpy as gp 

4 

5# Setup configuration. BCC Lattice 

6configuration = gp.Configuration(D=2) 

7configuration.make_lattice(unit_cell=gp.unit_cells.HEXAGONAL, cells=[16, 10], rho=1.0) 

8 

9# Setup masses and velocities 

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

11configuration.randomize_velocities(temperature=0.7 * 2) 

12 

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

14pair_func = gp.apply_shifted_force_cutoff(gp.LJ_12_6_sigma_epsilon) 

15sig, eps, cut = 1.0, 1.0, 2.5 

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

17 

18# Setup integrator 

19integrator = gp.integrators.NVT(temperature=0.7, tau=0.2, dt=0.005) 

20 

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

22runtime_actions = [gp.TrajectorySaver(), 

23 gp.ScalarSaver(), 

24 gp.MomentumReset(100)] 

25 

26 

27# Setup Simulation. 

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

29 num_timeblocks=32, steps_per_timeblock=32*1024, 

30 storage='memory') 

31 

32# Plot initial configuration 

33plt.figure() 

34 

35# Run simulation 

36sim.run() 

37 

38# Plot final configuration and box 

39plt.plot(configuration['r'][:, 0], configuration['r'][:, 1], 'o', 

40 color='blue', label='Final configuration') 

41box_lengths = configuration.simbox.get_lengths() 

42x_min, x_max = -box_lengths[0] / 2, box_lengths[0] / 2 

43y_min, y_max = -box_lengths[1] / 2, box_lengths[1] / 2 

44x_vals = [x_min, x_max, x_max, x_min, x_min] 

45y_vals = [y_min, y_min, y_max, y_max, y_min] 

46plt.plot(x_vals, y_vals, 'k--', label='Box') 

47plt.axis('equal') 

48plt.xlabel(r'$x$') 

49plt.ylabel(r'$y$') 

50if __name__ == "__main__": 

51 plt.show()