1""" Read/Write DL_POLY_4 CONFIG files """ 

2import re 

3 

4from numpy import zeros, isscalar 

5 

6from ase.atoms import Atoms 

7from ase.units import _auf, _amu, _auv 

8from ase.data import chemical_symbols 

9from ase.calculators.singlepoint import SinglePointCalculator 

10 

11__all__ = ['read_dlp4', 'write_dlp4'] 

12 

13# dlp4 labels will be registered in atoms.arrays[DLP4_LABELS_KEY] 

14DLP4_LABELS_KEY = 'dlp4_labels' 

15DLP4_DISP_KEY = 'dlp4_displacements' 

16dlp_f_si = 1.0e-10 * _amu / (1.0e-12 * 1.0e-12) # Å*Da/ps^2 

17dlp_f_ase = dlp_f_si / _auf 

18dlp_v_si = 1.0e-10 / 1.0e-12 # Å/ps 

19dlp_v_ase = dlp_v_si / _auv 

20 

21 

22def _get_frame_positions(f): 

23 """Get positions of frames in HISTORY file.""" 

24 # header line contains name of system 

25 init_pos = f.tell() 

26 f.seek(0) 

27 rl = len(f.readline()) # system name, and record size 

28 items = f.readline().strip().split() 

29 if len(items) == 5: 

30 classic = False 

31 elif len(items) == 3: 

32 classic = True 

33 else: 

34 raise RuntimeError("Cannot determine version of HISTORY file format.") 

35 

36 levcfg, imcon, natoms = [int(x) for x in items[0:3]] 

37 if classic: 

38 # we have to iterate over the entire file 

39 startpos = f.tell() 

40 pos = [] 

41 line = True 

42 while line: 

43 line = f.readline() 

44 if 'timestep' in line: 

45 pos.append(f.tell()) 

46 f.seek(startpos) 

47 else: 

48 nframes = int(items[3]) 

49 pos = [((natoms * (levcfg + 2) + 4) * i + 3) 

50 * rl for i in range(nframes)] 

51 

52 f.seek(init_pos) 

53 return levcfg, imcon, natoms, pos 

54 

55 

56def read_dlp_history(f, index=-1, symbols=None): 

57 """Read a HISTORY file. 

58 

59 Compatible with DLP4 and DLP_CLASSIC. 

60 

61 *Index* can be integer or slice. 

62 

63 Provide a list of element strings to symbols to ignore naming 

64 from the HISTORY file. 

65 """ 

66 levcfg, imcon, natoms, pos = _get_frame_positions(f) 

67 if isscalar(index): 

68 selected = [pos[index]] 

69 else: 

70 selected = pos[index] 

71 

72 images = [] 

73 for fpos in selected: 

74 f.seek(fpos + 1) 

75 images.append(read_single_image(f, levcfg, imcon, natoms, 

76 is_trajectory=True, symbols=symbols)) 

77 

78 return images 

79 

80 

81def iread_dlp_history(f, symbols=None): 

82 """Generator version of read_history""" 

83 levcfg, imcon, natoms, pos = _get_frame_positions(f) 

84 for p in pos: 

85 f.seek(p + 1) 

86 yield read_single_image(f, levcfg, imcon, natoms, is_trajectory=True, 

87 symbols=symbols) 

88 

89 

90def read_dlp4(f, symbols=None): 

91 """Read a DL_POLY_4 config/revcon file. 

92 

93 Typically used indirectly through read('filename', atoms, format='dlp4'). 

94 

95 Can be unforgiven with custom chemical element names. 

96 Please complain to alin@elena.space for bugs.""" 

97 

98 line = f.readline() 

99 line = f.readline().split() 

100 levcfg = int(line[0]) 

101 imcon = int(line[1]) 

102 

103 position = f.tell() 

104 line = f.readline() 

105 tokens = line.split() 

106 is_trajectory = tokens[0] == 'timestep' 

107 if not is_trajectory: 

108 # Difficult to distinguish between trajectory and non-trajectory 

109 # formats without reading one line too much. 

110 f.seek(position) 

111 

112 while line: 

113 if is_trajectory: 

114 tokens = line.split() 

115 natoms = int(tokens[2]) 

116 else: 

117 natoms = None 

118 yield read_single_image(f, levcfg, imcon, natoms, is_trajectory, 

119 symbols) 

120 line = f.readline() 

121 

122 

123def read_single_image(f, levcfg, imcon, natoms, is_trajectory, symbols=None): 

124 cell = zeros((3, 3)) 

125 ispbc = imcon > 0 

126 if ispbc or is_trajectory: 

127 for j in range(3): 

128 line = f.readline() 

129 line = line.split() 

130 for i in range(3): 

131 try: 

132 cell[j, i] = float(line[i]) 

133 except ValueError: 

134 raise RuntimeError("error reading cell") 

135 if symbols: 

136 sym = symbols 

137 else: 

138 sym = [] 

139 positions = [] 

140 velocities = [] 

141 forces = [] 

142 charges = [] 

143 masses = [] 

144 disp = [] 

145 

146 if is_trajectory: 

147 counter = range(natoms) 

148 else: 

149 from itertools import count 

150 counter = count() 

151 

152 labels = [] 

153 mass = None 

154 ch = None 

155 d = None 

156 

157 for a in counter: 

158 line = f.readline() 

159 if not line: 

160 a -= 1 

161 break 

162 

163 m = re.match(r'\s*([A-Za-z][a-z]?)(\S*)', line) 

164 assert m is not None, line 

165 symbol, label = m.group(1, 2) 

166 symbol = symbol.capitalize() 

167 if is_trajectory: 

168 ll = line.split() 

169 if len(ll) == 5: 

170 mass, ch, d = [float(i) for i in line.split()[2:5]] 

171 else: 

172 mass, ch = [float(i) for i in line.split()[2:4]] 

173 

174 charges.append(ch) 

175 masses.append(mass) 

176 disp.append(d) 

177 if not symbols: 

178 assert symbol in chemical_symbols 

179 sym.append(symbol) 

180 # make sure label is not empty 

181 if label: 

182 labels.append(label) 

183 else: 

184 labels.append('') 

185 

186 x, y, z = f.readline().split()[:3] 

187 positions.append([float(x), float(y), float(z)]) 

188 if levcfg > 0: 

189 vx, vy, vz = f.readline().split()[:3] 

190 velocities.append([float(vx) * dlp_v_ase, float(vy) 

191 * dlp_v_ase, float(vz) * dlp_v_ase]) 

192 if levcfg > 1: 

193 fx, fy, fz = f.readline().split()[:3] 

194 forces.append([float(fx) * dlp_f_ase, float(fy) 

195 * dlp_f_ase, float(fz) * dlp_f_ase]) 

196 

197 if symbols and (a + 1 != len(symbols)): 

198 raise ValueError("Counter is at {} but you gave {} symbols" 

199 .format(a + 1, len(symbols))) 

200 

201 if imcon == 0: 

202 pbc = False 

203 elif imcon == 6: 

204 pbc = [True, True, False] 

205 else: 

206 assert imcon in [1, 2, 3] 

207 pbc = True 

208 

209 atoms = Atoms(positions=positions, 

210 symbols=sym, 

211 cell=cell, 

212 pbc=pbc, 

213 # Cell is centered around (0, 0, 0) in dlp4: 

214 celldisp=-cell.sum(axis=0) / 2) 

215 

216 if is_trajectory: 

217 atoms.set_masses(masses) 

218 atoms.set_array(DLP4_DISP_KEY, disp, float) 

219 atoms.set_initial_charges(charges) 

220 

221 atoms.set_array(DLP4_LABELS_KEY, labels, str) 

222 if levcfg > 0: 

223 atoms.set_velocities(velocities) 

224 if levcfg > 1: 

225 atoms.calc = SinglePointCalculator(atoms, forces=forces) 

226 return atoms 

227 

228 

229def write_dlp4(fd, atoms, levcfg=0, title='CONFIG generated by ASE'): 

230 """Write a DL_POLY_4 config file. 

231 

232 Typically used indirectly through write('filename', atoms, format='dlp4'). 

233 

234 Can be unforgiven with custom chemical element names. 

235 Please complain to alin@elena.space in case of bugs""" 

236 

237 fd.write('{0:72s}\n'.format(title)) 

238 natoms = len(atoms) 

239 npbc = sum(atoms.pbc) 

240 if npbc == 0: 

241 imcon = 0 

242 elif npbc == 3: 

243 imcon = 3 

244 elif tuple(atoms.pbc) == (True, True, False): 

245 imcon = 6 

246 else: 

247 raise ValueError('Unsupported pbc: {}. ' 

248 'Supported pbc are 000, 110, and 000.' 

249 .format(atoms.pbc)) 

250 

251 fd.write('{0:10d}{1:10d}{2:10d}\n'.format(levcfg, imcon, natoms)) 

252 if imcon > 0: 

253 cell = atoms.get_cell() 

254 for j in range(3): 

255 fd.write('{0:20.10f}{1:20.10f}{2:20.10f}\n'.format( 

256 cell[j, 0], cell[j, 1], cell[j, 2])) 

257 vels = [] 

258 forces = [] 

259 if levcfg > 0: 

260 vels = atoms.get_velocities() / dlp_v_ase 

261 if levcfg > 1: 

262 forces = atoms.get_forces() / dlp_f_ase 

263 

264 labels = atoms.arrays.get(DLP4_LABELS_KEY) 

265 

266 for i, a in enumerate(atoms): 

267 sym = a.symbol 

268 if labels is not None: 

269 sym += labels[i] 

270 fd.write("{0:8s}{1:10d}\n{2:20.10f}{3:20.10f}{4:20.10f}\n".format( 

271 sym, a.index + 1, a.x, a.y, a.z)) 

272 if levcfg > 0: 

273 if vels is None: 

274 fd.write("{0:20.10f}{1:20.10f}{2:20.10f}\n".format( 

275 0.0, 0.0, 0.0)) 

276 else: 

277 fd.write("{0:20.10f}{1:20.10f}{2:20.10f}\n".format( 

278 vels[a.index, 0], vels[a.index, 1], vels[a.index, 2])) 

279 if levcfg > 1: 

280 if forces is None: 

281 fd.write("{0:20.10f}{1:20.10f}{2:20.10f}\n".format( 

282 0.0, 0.0, 0.0)) 

283 else: 

284 fd.write("{0:20.10f}{1:20.10f}{2:20.10f}\n".format( 

285 forces[a.index, 0], forces[a.index, 1], forces[a.index, 2]))