Coverage for examples/rubber_cube.py: 83%

1""" A rubber cube modeled as particles connected by springs """ 

2 

3from itertools import product 

4from math import pi, sin, cos 

5 

6import numpy as np 

7import matplotlib.pyplot as plt 

8 

9import gamdpy as gp 

10 

11# A rotated cube of particles 

12lattice_length = 8 # Number of particles in each direction 

13cm = np.array([0.0, 0.0, 0.0]) # Initial position of center of mass 

14theta_z = pi/12 # 1st rotation around z-axis 

15theta_x = pi/24 # 2nd rotation around x-axis 

16cube = np.array([(x, y, z) for x, y, z in product(range(lattice_length), repeat=3)]) 

17center = cube.mean(axis=0) 

18rot_matrix_z = np.array([ 

19 [cos(theta_z), sin(theta_z), 0], 

20 [-sin(theta_z), cos(theta_z), 0], 

21 [0, 0, 1] 

22]) 

23rot_matrix_x = np.array([ 

24 [1, 0, 0], 

25 [0, cos(theta_z), sin(theta_z)], 

26 [0, -sin(theta_z), cos(theta_z)], 

27]) 

28N = lattice_length**3 

29cube = (cube - center) @ rot_matrix_z @ rot_matrix_x + cm 

30configuration = gp.Configuration(D=3, N=N) 

31configuration['r'] = cube 

32configuration['m'] = 1.0 

33box_length = 64 

34configuration.simbox = gp.Orthorhombic(3, [box_length, box_length, box_length]) 

35configuration.randomize_velocities(temperature=0.001) 

36# Compute plan 

37compute_plan = {'pb': 16, 'tp': 4, 'skin': 2.0, 'UtilizeNIII': False, 'gridsync': False, 'nblist': 'N squared'} 

38 

39# Create bonds 

40bond_potential = gp.harmonic_bond_function 

41bond_params = [ 

42 [1.0, 100.0], 

43 [2**0.5, 10.0] 

44] 

45bond_directions =[ 

46 # (dx, dy, dz), bond_type 

47 ((1, 0, 0), 0), 

48 ((0, 1, 0), 0), 

49 ((0, 0, 1), 0), 

50 ((0, 1, 1), 1), 

51 ((1, 0, 1), 1), 

52 ((1, 1, 0), 1) 

53] 

54L = lattice_length 

55def xyz2idx(x, y, z): 

56 return x + y*L + z*L**2 

57neighbour_bonds = [] 

58for x, y, z in product(range(L), repeat=3): 

59 this = xyz2idx(x, y, z) 

60 for (dx, dy, dz), bond_type in bond_directions: 

61 if x + dx < L and y + dy < L and z + dz < L: 

62 that = xyz2idx(x+dx, y+dy, z+dz) 

63 neighbour_bonds.append([this, that, bond_type]) 

64bonds = gp.Bonds(bond_potential, bond_params, neighbour_bonds) 

65 

66# Create a 3D figure showing bonds 

67plot = False 

68if plot: 

69 fig = plt.figure() 

70 ax = fig.add_subplot(111, projection='3d') 

71 pos = configuration['r'] 

72 for bond in neighbour_bonds: 

73 n, m, _ = bond 

74 xs = [pos[n, 0], pos[m, 0]] 

75 ys = [pos[n, 1], pos[m, 1]] 

76 zs = [pos[n, 2], pos[m, 2]] 

77 ax.plot(xs, ys, zs, 'k-', linewidth=1) 

78 ax.scatter(cube[:, 0], cube[:, 1], cube[:, 2], c=cube[:, 2]) 

79 plt.show() 

80 

81# Add two smooth walls 

82wall_distance = box_length/2 

83walls = gp.interactions.Planar( 

84 potential=gp.harmonic_repulsion, 

85 params=[[100.0, 1.0], [100.0, 1.0]], 

86 indices=[[n, 0] for n in range(N)] + [[n, 1] for n in range(N)], # All particles feel both walls 

87 normal_vectors=[[0.0, 1.0, 0.0], [0.0, -1.0, 0.0]], 

88 points=[[0.0, -wall_distance/2, 0], [0.0, wall_distance/2, 0]] 

89) 

90 

91# Add gravity 

92mg = 0.0005 

93potential_gravity = gp.make_IPL_n(-1) 

94gravity = gp.interactions.Planar( 

95 potential=potential_gravity, 

96 params= [[mg, 10*wall_distance]], 

97 indices= [[n, 0] for n in range(N)], # All particles feel the gravity 

98 normal_vectors= [[0,1,0], ], 

99 points= [[0, -wall_distance/2.0, 0] ] # Defining 0 for potential energy on lower wall 

100) 

101 

102# Setup simulation 

103integrator = gp.integrators.NVE(dt=0.01) 

104runtime_actions = [gp.TrajectorySaver(), 

105 gp.ScalarSaver(steps_between_output=1)] 

106interactions = [bonds, walls, gravity] 

107sim = gp.Simulation(configuration, interactions, integrator, runtime_actions, 

108 num_timeblocks=64, steps_per_timeblock=1024, 

109 storage='memory') 

110 

111# Run simulation and save as lammps trajectory 

112dump_filename = 'Data/rubber_cube.lammps' 

113with open(dump_filename, 'w') as f: 

114 print(gp.configuration_to_lammps(sim.configuration, timestep=0), file=f) 

115 

116for block in sim.run_timeblocks(): 

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

118 with open(dump_filename, 'a') as f: 

119 print(gp.configuration_to_lammps(sim.configuration, timestep=sim.steps_per_block*(block+1)), file=f) 

120 

121print(sim.summary())