Coverage for examples/evaluator_einstein_crystal.py: 100%

1""" 

2Simulate Lennard-Jones crystal and evaluate the harmonic potential.  

3 

4This script sets up a configuration with a face-centered cubic (FCC) lattice  

5and randomizes the velocities. The pair potential is the single component 12-6  

6Lennard-Jones potential. The integrator is NVT.  

7An evaluator is created to calculate the potential energy of the system  

8using a none-interacting potential (eps=0.0) and harmonic springs.  

9The simulation is then performed, and the mean potential energy of the reference 

10einstein crystal (harmonic springs) is calculated and printed. 

11The mean displacement from the ideal lattice, sqrt(2*u_spring), is calculated and printed. 

12 

13""" 

14 

15import numpy as np 

16 

17import gamdpy as gp 

18 

19# Setup configuration: FCC Lattice 

20configuration = gp.Configuration(D=3) 

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

22configuration['m'] = 1.0 

23configuration.randomize_velocities(temperature=0.7) 

24# anchor_points = np.array(configuration['r']).copy() 

25 

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

27pair_func = gp.apply_shifted_potential_cutoff(gp.LJ_12_6_sigma_epsilon) 

28eps, sig, cut = 1.0, 1.0, 2.5 

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

30 

31# Setup integrator: NVT 

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

33 

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

35runtime_actions = [gp.TrajectorySaver(), 

36 gp.ScalarSaver(32), 

37 gp.MomentumReset(100)] 

38 

39# Setup Simulation. 

40sim = gp.Simulation( 

41 configuration, pair_pot, integrator, runtime_actions, 

42 num_timeblocks=16, 

43 steps_per_timeblock=1024, 

44 storage='memory' 

45) 

46 

47# Create evaluator for einstein crystal 

48# (replace with your potential of interest) 

49none_interacting = pair_pot = gp.PairPotential(pair_func, params=[0.0, 1.0, 2.5], max_num_nbs=1000) 

50harmonic_springs = gp.Tether() # U = 0.5*k*(r-r0)^2 

51harmonic_springs.set_anchor_points_from_lists( 

52 particle_indices=list(range(configuration.N)), 

53 spring_constants=[1.0]*configuration.N, 

54 configuration=configuration 

55) 

56evaluator = gp.Evaluator(sim.configuration, [none_interacting, harmonic_springs]) 

57 

58# Run simulation 

59u_spring = [] 

60displacements = [] 

61for block in sim.run_timeblocks(): 

62 evaluator.evaluate(sim.configuration) 

63 this_u_spring = evaluator.configuration['U'] 

64 u_spring.append(np.sum(this_u_spring)) 

65 

66 # Use the harmonic potential energy to calculate 

67 # the displacement relative to the ideal lattice 

68 dr = np.sqrt(2*this_u_spring) 

69 displacements.append(dr) 

70 

71print(f'Mean harmonic potential energy: {np.mean(u_spring)}') 

72print(f'Mean displacement from ideal lattice: {np.mean(displacements)}') 

73