Coverage for examples/isochore.py: 100%
24 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
« 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.
3Simulation of heating a Lennard-Jones crystal on an isochore in the NVT ensemble.
4For an even simpler script, see minimal.py
6"""
8import gamdpy as gp
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)
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)
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
30for temperature in ['0.70', '1.10', '1.50']:
31 print('\n\nTemperature: ' + temperature)
33 # Setup integrator
34 integrator = gp.integrators.NVT(temperature=temperature, tau=0.2, dt=0.005)
36 # Setup runtime actions, i.e. actions performed during simulation of timeblocks
37 runtime_actions = [gp.TrajectorySaver(),
38 gp.ScalarSaver(),
39 gp.MomentumReset(100)]
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')
46 print('Equilibration:')
47 for block in sim.run_timeblocks():
48 print(sim.status(per_particle=True))
49 print(sim.summary())
51 print('Production:')
52 for block in sim.run_timeblocks():
53 print(sim.status(per_particle=True))
54 print(sim.summary())
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