Coverage for examples/molecules.py: 100%

1import numpy as np 

2import gamdpy as gp 

3 

4gp.select_gpu() 

5 

6# Simulation params  

7rho, temperature = 0.85, 1.5 

8N_A, N_B, N_C = 8, 4, 4 # Number of atoms of each tyoe 

9particles_per_molecule = N_A + N_B + N_C 

10filename = 'Data/molecules' 

11num_timeblocks = 64 

12steps_per_timeblock = 1 * 1024 # 8 * 1024 to show reliable pattern formation 

13 

14positions = [] 

15particle_types = [] 

16masses = [] 

17 

18# A particles 

19for i in range(N_A): 

20 positions.append( [ i*1.0, (i%2)*.1, 0. ] ) # x, y, z for this particle 

21 particle_types.append( 0 ) 

22 masses.append( 1.0 ) 

23 

24# B particles 

25for i in range(N_B): 

26 positions.append( [ 0, (i+1)*1.0, ((i+1)%2)*.1 ] ) # x, y, z for this particle 

27 particle_types.append( 1 ) 

28 masses.append( 1.0 ) 

29 

30# C particles 

31for i in range(N_C): 

32 positions.append( [ ((i+1)%2)*.1, 0, (i+1)*1.0 ] ) # x, y, z for this particle 

33 particle_types.append( 2 ) 

34 masses.append( 1.0 ) 

35 

36# Setup configuration: Single molecule first, then duplicate 

37top = gp.Topology(['MyMolecule', ]) 

38top.bonds = gp.bonds_from_positions(positions, cut_off=1.1, bond_type=0) 

39top.angles = gp.angles_from_bonds(top.bonds, angle_type=0) 

40top.dihedrals = gp.dihedrals_from_angles(top.angles, dihedral_type=0) 

41top.molecules['MyMolecule'] = gp.molecules_from_bonds(top.bonds) 

42 

43 

44dict_this_mol = {"positions" : positions, 

45 "particle_types" : particle_types, 

46 "masses" : masses, 

47 "topology" : top} 

48 

49print('Initial Positions:') 

50for position in positions: 

51 print('\t\t', position) 

52print('Particle types:\t', particle_types) 

53print('Bonds: \t', top.bonds) 

54print('Angles: \t', top.angles) 

55print('Dihedrals: \t', top.dihedrals) 

56print() 

57 

58# This call creates the pdf "molecule.pdf" with a drawing of the molecule  

59# Use block=True to visualize the molecule before running the simulation 

60gp.plot_molecule(top, positions, particle_types, filename="molecule.pdf", block=False) 

61 

62configuration = gp.replicate_molecules([dict_this_mol], [216], safety_distance=2.0, compute_flags={"stresses":True}) 

63 

64configuration.randomize_velocities(temperature=temperature) 

65 

66print(f'Number of molecules: {len(configuration.topology.molecules["MyMolecule"])}') 

67print(f'Number of particles: {configuration.N}\n') 

68 

69# Make bond interactions 

70bond_potential = gp.harmonic_bond_function 

71bond_params = [[0.8, 1000.], ] 

72bonds = gp.Bonds(bond_potential, bond_params, configuration.topology.bonds) 

73 

74# Make angle interactions 

75angle_potential = gp.cos_angle_function 

76angle0, k = 2.0, 500.0 

77angles = gp.Angles(angle_potential, configuration.topology.angles, parameters=[[k, angle0],]) 

78 

79# Make dihedral interactions 

80dihedral_potential = gp.ryckbell_dihedral 

81rbcoef=[.0, 5.0, .0, .0, .0, .0] 

82dihedrals = gp.Dihedrals(dihedral_potential, configuration.topology.dihedrals, parameters=[rbcoef, ]) 

83 

84# Exlusion list 

85exclusions = dihedrals.get_exclusions(configuration) 

86#exclusions = bonds.get_exclusions(configuration) 

87 

88# Make pair potential 

89pair_func = gp.apply_shifted_force_cutoff(gp.LJ_12_6_sigma_epsilon) 

90sig = [[1.00, 1.00, 1.00], 

91 [1.00, 1.00, 1.00], 

92 [1.00, 1.00, 1.00],] 

93eps = [[1.00, 1.00, 1.00], 

94 [1.00, 1.00, 0.80], 

95 [1.00, 0.80, 1.00],] 

96cut = [[2.50, 1.12, 1.12], 

97 [1.12, 2.50, 2.50], 

98 [1.12, 2.50, 2.50]] 

99 

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

101 

102# Make integrator 

103integrator = gp.integrators.NVT(temperature=temperature, tau=0.1, dt=0.004) 

104 

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

106runtime_actions = [gp.TrajectorySaver(), 

107 gp.ScalarSaver(32), 

108 gp.StressSaver(32, compute_flags={'stresses':True}), 

109 gp.MomentumReset(100)] 

110 

111# Setup simulation 

112sim = gp.Simulation(configuration, [pair_pot, bonds, angles, dihedrals], integrator, runtime_actions, 

113 num_timeblocks=num_timeblocks, steps_per_timeblock=steps_per_timeblock, 

114 storage='memory') 

115 

116print('\nCompression and equilibration: ') 

117dump_filename = 'Data/dump_compress.lammps' 

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

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

120 

121initial_rho = configuration.N / configuration.get_volume() 

122for block in sim.run_timeblocks(): 

123 volume = configuration.get_volume() 

124 N = configuration.N 

125 print(sim.status(per_particle=True), f'rho= {N/volume:.3}', end='\t') 

126 print(f'P= {(N*temperature + np.sum(configuration["W"]))/volume:.3}') # pV = NkT + W 

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

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

129 

130 # Scale configuration to get closer to final density, rho 

131 if block<sim.num_blocks/2: 

132 desired_rho = (block+1)/(sim.num_blocks/2)*(rho - initial_rho) + initial_rho 

133 configuration.atomic_scale(density=desired_rho) 

134 configuration.copy_to_device() # Since we altered configuration, we need to copy it back to device 

135print(sim.summary()) 

136print(configuration) 

137 

138sim = gp.Simulation(configuration, [pair_pot, bonds, angles, dihedrals], integrator, runtime_actions, 

139 num_timeblocks=num_timeblocks, steps_per_timeblock=steps_per_timeblock, 

140 compute_plan=sim.compute_plan, storage=filename+'.h5') 

141 

142print('\nProduction: ') 

143dump_filename = 'Data/dump.lammps' 

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

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

146 

147for block in sim.run_timeblocks(): 

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

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

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

151 

152print(sim.summary()) 

153print(configuration) 

154 

155W = gp.extract_scalars(sim.output, 'W') 

156full_stress_tensor = gp.extract_stress_tensor(sim.output) 

157mean_diagonal_sts = (full_stress_tensor[:,0,0] + full_stress_tensor[:,1,1] + full_stress_tensor[:,2,2])/3 

158 

159print("Mean diagonal stress", np.mean(mean_diagonal_sts) ) 

160print("Pressure", np.mean(W)*rho/N + temperature*rho) 

161 

162print('\nAnalyse structure with:') 

163print(' python3 analyze_structure.py Data/molecules') 

164 

165print('\nAnalyze dynamics with:') 

166print(' python3 analyze_dynamics.py Data/molecules') 

167 

168print('\nVisualize simulation in ovito with:') 

169print(' ovito Data/dump.lammps') 

170 

171#print('\nVisualize simulation in VMD with:') 

172#print(' vmd -lammpstrj Data/dump.lammps')