Coverage for examples/ASD.py: 100%

1""" Simulate a system of asymmetric dumbbells (ASD) with a Lennard-Jones potential and a harmonic bond potential. """ 

2 

3import numpy as np 

4import gamdpy as gp 

5 

6# Specify state point 

7rho = 1.863 # Atomic density = Molecular density * 2 

8temperature = 0.465 

9filename = f'Data/ASD_rho{rho:.3f}_T{temperature:.3f}.h5' 

10 

11# Generate two-component configuration with a FCC lattice  

12configuration = gp.Configuration(D=3, compute_flags={'Fsq':True, 'lapU':True, 'Ptot':True, 'Vol':True}) 

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

14configuration['m'] = 1.0 

15B_particles = range(1, configuration.N, 2) 

16configuration.ptype[B_particles] = 1 # Setting particle type of B particles 

17configuration['m'][B_particles] = 0.195 # Setting masses of B particles 

18configuration.randomize_velocities(temperature=1.44) 

19 

20# Make bonds 

21bond_potential = gp.harmonic_bond_function 

22bond_params = [[0.584, 3000.], ] # Parameters for bond type 0, 1, 2 etc (here only 0) 

23bond_indices = [[i, i + 1, 0] for i in range(0, configuration.N - 1, 2)] # dumbells: i(even) and i+1 bonded with type 0 

24bonds = gp.Bonds(bond_potential, bond_params, bond_indices) 

25 

26# Make pair potential 

27# pair_func = gp.apply_shifted_potential_cutoff(gp.LJ_12_6_sigma_epsilon) 

28pair_func = gp.apply_shifted_force_cutoff(gp.LJ_12_6_sigma_epsilon) 

29sig = [[1.000, 0.894], 

30 [0.894, 0.788]] 

31eps = [[1.000, 0.342], 

32 [0.342, 0.117]] 

33cut = np.array(sig) * 2.5 

34exclusions = bonds.get_exclusions(configuration) 

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

36 

37# Make integrator 

38dt = 0.002 # timestep 

39num_blocks = 64 # Do simulation in this many 'blocks' 

40steps_per_block = 1*1024 # ... each of this many steps (increase for better statistics) 

41running_time = dt * num_blocks * steps_per_block 

42 

43Ttarget_function = gp.make_function_ramp(value0=10.000, x0=running_time * (1 / 8), 

44 value1=temperature, x1=running_time * (1 / 4)) 

45integrator0 = gp.integrators.NVT(Ttarget_function, tau=0.2, dt=dt) 

46 

47sim = gp.Simulation(configuration, [pair_pot, bonds], integrator0, 

48 runtime_actions=[gp.MomentumReset(100)], 

49 num_timeblocks=num_blocks, steps_per_timeblock=steps_per_block, 

50 storage='memory') 

51 

52print('High Temperature followed by cooling and equilibration:') 

53for block in sim.run_timeblocks(): 

54 if block % 10 == 0: 

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

56print(sim.summary()) 

57 

58print('Production:') 

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

60 

61runtime_actions = [gp.MomentumReset(100), 

62 gp.TrajectorySaver(), 

63 gp.ScalarSaver(32, {'Fsq':True, 'lapU':True, 'Ptot':True}), ] 

64 

65sim = gp.Simulation(configuration, [pair_pot, bonds], integrator, runtime_actions, 

66 num_timeblocks=num_blocks, steps_per_timeblock=steps_per_block, 

67 storage=filename) 

68 

69for block in sim.run_timeblocks(): 

70 if block % 10 == 0: 

71 print(f'{block=:4} {sim.status(per_particle=True)}') 

72print(sim.summary())