Coverage for examples/isochore.py: 100%

24 statements  

« prev     ^ index     » next       coverage.py v7.9.1, created at 2025-06-14 15:55 +0200

1""" Simple example of performing several simulation in one go using gamdpy. 

2 

3Simulation of heating a Lennard-Jones crystal on an isochore in the NVT ensemble. 

4For an even simpler script, see minimal.py 

5 

6""" 

7 

8import gamdpy as gp 

9 

10# Setup fcc configuration 

11configuration = gp.Configuration(D=3) 

12configuration.make_lattice(gp.unit_cells.FCC, cells=[6, 6, 6], rho=0.973) 

13configuration['m'] = 1.0 

14configuration.randomize_velocities(temperature=1.6) 

15 

16# Setup pair potential. 

17pair_func = gp.apply_shifted_force_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# Specify duration of simulations.  

22# Increase 'num_timelocks' for longer runs, better statistics, AND larger storage consumption 

23# Increase 'steps_per_block' for longer runs 

24dt = 0.004 # timestep 

25num_timeblocks = 8 # Do simulation in this many 'timeblocks'.  

26steps_per_timeblock = 1*1024 # ... each of this many steps 

27 

28 

29 

30for temperature in ['0.70', '1.10', '1.50']: 

31 print('\n\nTemperature: ' + temperature) 

32 

33 # Setup integrator 

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

35 

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

37 runtime_actions = [gp.TrajectorySaver(), 

38 gp.ScalarSaver(), 

39 gp.MomentumReset(100)] 

40 

41 # Setup Simulation 

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

43 num_timeblocks=num_timeblocks, steps_per_timeblock=steps_per_timeblock, 

44 storage='Data/LJ_r0.973_T'+temperature+'.h5') 

45 

46 print('Equilibration:') 

47 for block in sim.run_timeblocks(): 

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

49 print(sim.summary()) 

50 

51 print('Production:') 

52 for block in sim.run_timeblocks(): 

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

54 print(sim.summary()) 

55 

56# To get a plot of the MSD do something like this: 

57# python3 -m gamdpy.tools.calc_dynamics -o Data/msd_r0.973.pdf Data/LJ_r0.973_T*.h5