Coverage for tests/test_stress.py: 100%

1""" Example of a Simulation using gamdpy, including the stress tensor (configurational part only). 

2 

3""" 

4 

5import gamdpy as gp 

6import numpy as np 

7 

8# Setup configuration: FCC Lattice 

9rho = 0.973 

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

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

12configuration['m'] = 1.0 

13configuration.randomize_velocities(temperature=0.8 * 2) 

14 

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

16pairfunc = gp.apply_shifted_force_cutoff(gp.LJ_12_6_sigma_epsilon) 

17sig, eps, cut = 1.0, 1.0, 2.5 

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

19 

20# Setup integrator: NVT 

21integrator = gp.integrators.NVT(temperature=0.70, tau=0.2, dt=0.005) 

22 

23num_timeblocks = 4 

24 

25runtime_actions = [gp.MomentumReset(100), 

26 gp.TrajectorySaver(), 

27 gp.ScalarSaver(32, {'stresses':True}), ] 

28 

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

30sim = gp.Simulation(configuration, pairpot, integrator, runtime_actions, 

31 num_timeblocks=num_timeblocks, steps_per_timeblock=128, 

32 storage='memory',) 

33 

34p_conf_array = np.zeros(num_timeblocks) 

35scalar_stress_array = np.zeros(num_timeblocks) 

36 

37# Run simulation one block at a time 

38for block in sim.run_timeblocks(): 

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

40 vol = configuration.get_volume() 

41 

42 sts_x_row = (np.sum(sim.configuration['sx'], axis=0)/vol) 

43 sts_y_row = (np.sum(sim.configuration['sy'], axis=0)/vol) 

44 sts_z_row = (np.sum(sim.configuration['sz'], axis=0)/vol) 

45 #print('Stress tensor (configurational part)') 

46 mean_diagonal = (sts_x_row[0] + sts_y_row[1] + sts_z_row[2])/3 

47 virial = np.mean(configuration['W']) 

48 scalar_stress_array[block] = mean_diagonal 

49 p_conf_array[block] = virial * rho 

50 

51assert np.allclose(p_conf_array, -scalar_stress_array)