Coverage for gamdpy/configuration/old_input_output.py: 92%

1import h5py 

2import gzip 

3import numpy as np 

4import numba 

5import math 

6from numba import cuda 

7 

8from .colarray import colarray 

9from ..simulation_boxes import Orthorhombic, LeesEdwards 

10from .topology import Topology, duplicate_topology, replicate_topologies 

11from ..simulation.get_default_compute_flags import get_default_compute_flags 

12from .Configuration import Configuration 

13 

14''' 

15def configuration_to_hdf5(configuration: Configuration, filename: str, meta_data=None) -> None: 

16 """ Write a configuration to a HDF5 file 

17 

18 Parameters 

19 ---------- 

20 

21 configuration : gamdpy.Configuration 

22 a gamdpy configuration object 

23 

24 filename : str 

25 filename of the output file .h5 

26 

27 meta_data : str 

28 not used in the function so far (default None) 

29 

30 Example 

31 ------- 

32 

33 >>> import os 

34 >>> import gamdpy as gp 

35 >>> conf = gp.Configuration(D=3) 

36 >>> conf.make_positions(N=10, rho=1.0) 

37 >>> gp.configuration_to_hdf5(configuration=conf, filename="final.h5") 

38 >>> os.remove("final.h5") # Removes file (for doctests) 

39 

40 """ 

41 

42 if not filename.endswith('.h5'): 

43 filename += '.h5' 

44 with h5py.File(filename, "w") as f: 

45 f.attrs['simbox'] = configuration.simbox.get_lengths() 

46 if meta_data is not None: 

47 for item in meta_data: 

48 f.attrs[item] = meta_data[item] 

49 

50 ds_r = f.create_dataset('r', shape=(configuration.N, configuration.D), dtype=np.float32) 

51 ds_v = f.create_dataset('v', shape=(configuration.N, configuration.D), dtype=np.float32) 

52 ds_p = f.create_dataset('ptype', shape=(configuration.N), dtype=np.int32) 

53 ds_m = f.create_dataset('m', shape=(configuration.N), dtype=np.float32) 

54 ds_r_im = f.create_dataset('r_im', shape=(configuration.N, configuration.D), dtype=np.int32) 

55 ds_r[:] = configuration['r'] 

56 ds_v[:] = configuration['v'] 

57 ds_p[:] = configuration.ptype 

58 ds_m[:] = configuration['m'] 

59 ds_r_im[:] = configuration.r_im 

60''' 

61 

62 

63def configuration_from_hdf5(filename: str, reset_images=False, compute_flags=None) -> Configuration: 

64 """ Read a configuration from a HDF5 file 

65 

66 Parameters 

67 ---------- 

68 

69 filename : str 

70 filename of the input file .h5 

71 

72 reset_images : bool 

73 if True set the images to zero (default False) 

74 

75 Returns 

76 ------- 

77 

78 configuration : gamdpy.Configuration 

79 a gamdpy configuration object 

80 

81 Example 

82 ------- 

83 

84 >>> import gamdpy as gp 

85 >>> conf = gp.configuration_from_hdf5("examples/Data/final.h5") 

86 >>> print(conf.D, conf.N, conf['r'][0]) # Print number of dimensions D, number of particles N and position of first particle 

87 3 10 [-0.7181449 -1.3644753 -1.5799187] 

88 

89 """ 

90 

91 if not filename.endswith('.h5'): 

92 raise ValueError('Filename not in HDF5 format') 

93 with h5py.File(filename, "r") as f: 

94 lengths = f.attrs['simbox'] 

95 r = f['r'][:] 

96 v = f['v'][:] 

97 ptype = f['ptype'][:] 

98 m = f['m'][:] 

99 r_im = f['r_im'][:] 

100 N, D = r.shape 

101 configuration = Configuration(D=D, compute_flags=compute_flags) 

102 configuration.simbox = Orthorhombic(D, lengths) 

103 configuration['r'] = r 

104 configuration['v'] = v 

105 configuration.ptype = ptype 

106 configuration['m'] = m 

107 if reset_images: 

108 configuration.r_im = np.zeros((N, D), dtype=np.int32) 

109 else: 

110 configuration.r_im = r_im 

111 return configuration 

112 

113 

114''' 

115def configuration_from_hdf5_group(f, group_name, reset_images=False, compute_flags=None) -> Configuration: 

116 """ Read a configuration from an open HDF5 file identified by group-name 

117 

118 Parameters 

119 ---------- 

120 f : HDF5 File 

121 open HDF5 open, as returned by h5py.File() 

122 

123 reset_images : bool 

124 if True set the images to zero (default False) 

125 

126 Returns 

127 ------- 

128 

129 configuration : gamdpy.Configuration 

130 a gamdpy configuration object 

131 

132 

133 Example: 

134 -------- 

135 

136 >>> import gamdpy as gp 

137 >>> output_file = h5py.File('examples/Data/LJ_r0.973_T0.70_toread.h5') 

138 >>> conf = gp.configuration_from_hdf5_group(output_file, 'restarts/restart0000') 

139 >>> print(conf.D, conf.N, conf['r'][0]) # Print number of dimensions D, number of particles N and position of first particle 

140 3 2048 [-6.384221 -6.3622074 -6.3125153] 

141 

142 """ 

143 

144 

145 vectors_array = f[group_name]['vectors'][:] 

146 _, N, D = vectors_array.shape 

147 configuration = Configuration(D=D, N=N, compute_flags=compute_flags) 

148 

149 

150 configuration.vector_columns = f[group_name]['vectors'].attrs['vector_columns'] 

151 configuration.scalar_columns = f[group_name]['scalars'].attrs['scalar_columns'] 

152 

153 configuration.ptype = f[group_name]['ptype'][:] 

154 

155 

156 scalars_array = f[group_name]['scalars'][:] 

157 configuration.scalars = scalars_array 

158 configuration.vectors.array = vectors_array 

159 

160 

161 simbox_name = f[group_name].attrs['simbox_name'] 

162 simbox_data = f[group_name].attrs['simbox_data'] 

163 

164 if simbox_name == 'Orthorhombic': 

165 configuration.simbox = Orthorhombic(D, simbox_data) 

166 elif simbox_name == 'LeesEdwards': 

167 box_shift_image = {True:0, False: int(simbox_data[D+1])} [reset_images] 

168 configuration.simbox = LeesEdwards(D, simbox_data[:D], simbox_data[D], box_shift_image) 

169 else: 

170 raise ValueError('simbox_name %s not recognized in group %s' % (simbox_name, group_name)) 

171 

172 if reset_images: 

173 configuration.r_im = np.zeros((N, D), dtype=np.int32) 

174 else: 

175 configuration.r_im = f[group_name]['r_im'][:] 

176 return configuration 

177''' 

178 

179def configuration_to_rumd3(configuration: Configuration, filename: str) -> None: 

180 """ Write a configuration to a RUMD3 file  

181 

182 Parameters 

183 ---------- 

184 

185 configuration : gamdpy.Configuration 

186 a gamdpy configuration object 

187 

188 filename : str 

189 filename of the output file .xyz.gz 

190 

191 Example 

192 ------- 

193 

194 >>> import os 

195 >>> import gamdpy as gp 

196 >>> conf = gp.Configuration(D=3) 

197 >>> conf.make_positions(N=10, rho=1.0) 

198 >>> gp.configuration_to_rumd3(configuration=conf, filename="restart.xyz.gz") 

199 >>> os.remove("restart.xyz.gz") # Removes file (for doctests) 

200 

201 """ 

202 N = configuration.N 

203 if configuration.D != 3: 

204 raise ValueError("Only D==3 is compatibale with RUMD-3") 

205 

206 r = configuration['r'] 

207 v = configuration['v'] 

208 ptype = configuration.ptype 

209 m = configuration['m'] 

210 r_im = configuration.r_im 

211 

212 num_types = max(ptype) + 1 # assumes consecutive types starting from zero 

213 # find corresponding masses assuming unique mass for each type as required by RUMD-3 

214 masses = np.ones(num_types, dtype=np.float32) 

215 for type in range(num_types): 

216 type_first_idx = np.where(ptype == type)[0][0] 

217 masses[type] = m[type_first_idx] 

218 

219 sim_box = configuration.simbox.get_lengths() 

220 if not filename.endswith('.gz'): 

221 filename += '.gz' 

222 

223 with gzip.open(filename, 'wt') as f: 

224 f.write('%d\n' % N) 

225 comment_line = 'ioformat=2 numTypes=%d' % (num_types) 

226 comment_line += ' sim_box=RectangularSimulationBox,%f,%f,%f' % (sim_box[0], sim_box[1], sim_box[2]) 

227 comment_line += ' mass=%f' % (masses[0]) 

228 for mass in masses[1:]: 

229 comment_line += ',%f' % mass 

230 comment_line += ' columns=type,x,y,z,imx,imy,imz,vx,vy,vz' 

231 comment_line += '\n' 

232 f.write(comment_line) 

233 for idx in range(N): 

234 line_out = '%d %.9f %.9f %.9f %d %d %d %f %f %f\n' % ( 

235 ptype[idx], r[idx, 0], r[idx, 1], r[idx, 2], r_im[idx, 0], r_im[idx, 1], r_im[idx, 2], v[idx, 0], 

236 v[idx, 1], 

237 v[idx, 2]) 

238 f.write(line_out) 

239 

240 

241def configuration_from_rumd3(filename: str, reset_images=False, compute_flags=None) -> Configuration: 

242 """ Read a configuration from a RUMD3 file  

243 

244 Parameters 

245 ---------- 

246 

247 filename : str 

248 filename of the output file .xyz.gz 

249 

250 Returns 

251 ------- 

252 

253 configuration : gamdpy.Configuration 

254 a gamdpy configuration object 

255 

256 Example 

257 ------- 

258 

259 >>> import gamdpy as gp 

260 >>> conf = gp.configuration_from_rumd3("examples/Data/NVT_N4000_T2.0_rho1.2_KABLJ_rumd3/TrajectoryFiles/restart0000.xyz.gz") 

261 >>> print(conf.D, conf.N, conf['r'][0]) # Print number of dimensions D, number of particles N and position of first particle 

262 3 4000 [ 7.197245 6.610052 -4.7467813] 

263 

264 """ 

265 with gzip.open(filename) as f: 

266 line1 = f.readline().decode() 

267 N = int(line1) 

268 

269 line2 = f.readline().decode() 

270 meta_data_items = line2.split() 

271 meta_data = {} 

272 for item in meta_data_items: 

273 key, val = item.split("=") 

274 meta_data[key] = val 

275 

276 num_types = int(meta_data['numTypes']) 

277 masses = [float(x) for x in meta_data['mass'].split(',')] 

278 assert len(masses) == num_types 

279 if meta_data['ioformat'] == '1': 

280 lengths = np.array([float(x) for x in meta_data['boxLengths'].split(',')], dtype=np.float32) 

281 else: 

282 assert meta_data['ioformat'] == '2' 

283 sim_box_data = meta_data['sim_box'].split(',') 

284 sim_box_type = sim_box_data[0] 

285 sim_box_params = [float(x) for x in sim_box_data[1:]] 

286 assert sim_box_type == 'RectangularSimulationBox' 

287 lengths = np.array(sim_box_params) 

288 # TO DO: handle LeesEdwards sim box 

289 assert meta_data['columns'].startswith('type,x,y,z,imx,imy,imz') 

290 has_velocities = (meta_data['columns'].startswith('type,x,y,z,imx,imy,imz,vx,vy,vz')) 

291 type_array = np.zeros(N, dtype=np.int32) 

292 r_array = np.zeros((N, 3), dtype=np.float32) 

293 im_array = np.zeros((N, 3), dtype=np.int32) 

294 v_array = np.zeros((N, 3), dtype=np.float32) 

295 m_array = np.ones(N, dtype=np.float32) 

296 

297 for idx in range(N): 

298 p_data = f.readline().decode().split() 

299 ptype = int(p_data[0]) 

300 type_array[idx] = ptype 

301 r_array[idx, :] = [float(x) for x in p_data[1:4]] 

302 if not reset_images: 

303 im_array[idx, :] = [int(x) for x in p_data[4:7]] 

304 if has_velocities: 

305 v_array[idx, :] = [float(x) for x in p_data[7:10]] 

306 m_array[idx] = masses[ptype] 

307 

308 configuration = Configuration(D=3, compute_flags=compute_flags) 

309 configuration.simbox = Orthorhombic(3, lengths) 

310 configuration['r'] = r_array 

311 configuration['v'] = v_array 

312 configuration.r_im = im_array 

313 configuration.ptype = type_array 

314 configuration['m'] = m_array 

315 

316 return configuration 

317 

318 

319def configuration_to_lammps(configuration, timestep=0) -> str: 

320 """ Convert a configuration to a string formatted as LAMMPS dump file  

321 

322 Parameters 

323 ---------- 

324 

325 configuration : gamdpy.Configuration 

326 a gamdpy configuration object 

327 

328 timestep : float 

329 time at which the configuration is saved 

330 

331 Returns 

332 ------- 

333 

334 str 

335 string formatted as LAMMPS dump file 

336 

337 Example 

338 ------- 

339 

340 >>> import gamdpy as gp 

341 >>> conf = gp.Configuration(D=3) 

342 >>> conf.make_positions(N=10, rho=1.0) 

343 >>> lmp_dump = gp.configuration_to_lammps(configuration=conf) 

344 

345 """ 

346 D = configuration.D 

347 if D != 3 and D!=2: 

348 raise ValueError('Only 3D and 2D configurations are supported') 

349 masses = configuration['m'] 

350 positions = configuration['r'] 

351 image_coordinates = configuration.r_im 

352 forces = configuration['f'] 

353 velocities = configuration['v'] 

354 ptypes = configuration.ptype 

355 simulation_box = configuration.simbox.get_lengths() 

356 

357 # Header 

358 header = f'ITEM: TIMESTEP\n{timestep:d}\n' 

359 number_of_atoms = positions.shape[0] 

360 header += f'ITEM: NUMBER OF ATOMS\n{number_of_atoms:d}\n' 

361 header += f'ITEM: BOX BOUNDS pp pp pp\n' 

362 for k in range(D): 

363 header += f'{-simulation_box[k] / 2:e} {simulation_box[k] / 2:e}\n' 

364 if D==2: 

365 header += f'{-1 / 2:e} {1 / 2:e}\n' 

366 # Atoms 

367 atom_data = 'ITEM: ATOMS id type mass x y z ix iy iz vx vy vz fx fy fz' 

368 for i in range(number_of_atoms): 

369 atom_data += f'\n{i + 1:d} {ptypes[i] + 1:d} {masses[i]:f} ' 

370 for k in range(D): 

371 atom_data += f'{positions[i, k]:f} ' 

372 if D==2: 

373 atom_data += f'{0.0:f} ' 

374 for k in range(D): 

375 atom_data += f'{image_coordinates[i, k]:d} ' 

376 if D==2: 

377 atom_data += f'{0.0:f} ' 

378 for k in range(D): 

379 atom_data += f'{velocities[i, k]:f} ' 

380 if D==2: 

381 atom_data += f'{0.0:f} ' 

382 for k in range(D): 

383 atom_data += f'{forces[i, k]:f} ' 

384 if D==2: 

385 atom_data += f'{0.0:f} ' 

386 #atom_data += '\n' 

387 # Combine header and atom lengths 

388 lammps_dump = header + atom_data 

389 return lammps_dump