1# Copyright (C) 2012, Jesper Friis, SINTEF 

2# (see accompanying license files for ASE). 

3"""Module to read and write atoms EON reactant.con files. 

4 

5See http://theory.cm.utexas.edu/eon/index.html for a description of EON. 

6""" 

7import os 

8from warnings import warn 

9from glob import glob 

10 

11import numpy as np 

12 

13from ase.atoms import Atoms 

14from ase.constraints import FixAtoms 

15from ase.geometry import cellpar_to_cell, cell_to_cellpar 

16from ase.utils import writer 

17 

18 

19def read_eon(fileobj, index=-1): 

20 """Reads an EON reactant.con file. If *fileobj* is the name of a 

21 "states" directory created by EON, all the structures will be read.""" 

22 if isinstance(fileobj, str): 

23 if (os.path.isdir(fileobj)): 

24 return read_states(fileobj) 

25 else: 

26 fd = open(fileobj) 

27 else: 

28 fd = fileobj 

29 

30 more_images_to_read = True 

31 images = [] 

32 

33 first_line = fd.readline() 

34 while more_images_to_read: 

35 

36 comment = first_line.strip() 

37 fd.readline() # 0.0000 TIME (??) 

38 cell_lengths = fd.readline().split() 

39 cell_angles = fd.readline().split() 

40 # Different order of angles in EON. 

41 cell_angles = [cell_angles[2], cell_angles[1], cell_angles[0]] 

42 cellpar = [float(x) for x in cell_lengths + cell_angles] 

43 fd.readline() # 0 0 (??) 

44 fd.readline() # 0 0 0 (??) 

45 ntypes = int(fd.readline()) # number of atom types 

46 natoms = [int(n) for n in fd.readline().split()] 

47 atommasses = [float(m) for m in fd.readline().split()] 

48 

49 symbols = [] 

50 coords = [] 

51 masses = [] 

52 fixed = [] 

53 for n in range(ntypes): 

54 symbol = fd.readline().strip() 

55 symbols.extend([symbol] * natoms[n]) 

56 masses.extend([atommasses[n]] * natoms[n]) 

57 fd.readline() # Coordinates of Component n 

58 for i in range(natoms[n]): 

59 row = fd.readline().split() 

60 coords.append([float(x) for x in row[:3]]) 

61 fixed.append(bool(int(row[3]))) 

62 

63 atoms = Atoms(symbols=symbols, 

64 positions=coords, 

65 masses=masses, 

66 cell=cellpar_to_cell(cellpar), 

67 constraint=FixAtoms(mask=fixed), 

68 info=dict(comment=comment)) 

69 

70 images.append(atoms) 

71 

72 first_line = fd.readline() 

73 if first_line == '': 

74 more_images_to_read = False 

75 

76 if isinstance(fileobj, str): 

77 fd.close() 

78 

79 if not index: 

80 return images 

81 else: 

82 return images[index] 

83 

84 

85def read_states(states_dir): 

86 """Read structures stored by EON in the states directory *states_dir*.""" 

87 subdirs = glob(os.path.join(states_dir, '[0123456789]*')) 

88 subdirs.sort(key=lambda d: int(os.path.basename(d))) 

89 images = [read_eon(os.path.join(subdir, 'reactant.con')) 

90 for subdir in subdirs] 

91 return images 

92 

93 

94@writer 

95def write_eon(fileobj, images): 

96 """Writes structure to EON reactant.con file 

97 Multiple snapshots are allowed.""" 

98 if isinstance(images, Atoms): 

99 atoms = images 

100 elif len(images) == 1: 

101 atoms = images[0] 

102 else: 

103 raise ValueError('Can only write one configuration to EON ' 

104 'reactant.con file') 

105 

106 out = [] 

107 out.append(atoms.info.get('comment', 'Generated by ASE')) 

108 out.append('0.0000 TIME') # ?? 

109 

110 a, b, c, alpha, beta, gamma = cell_to_cellpar(atoms.cell) 

111 out.append('%-10.6f %-10.6f %-10.6f' % (a, b, c)) 

112 out.append('%-10.6f %-10.6f %-10.6f' % (gamma, beta, alpha)) 

113 

114 out.append('0 0') # ?? 

115 out.append('0 0 0') # ?? 

116 

117 symbols = atoms.get_chemical_symbols() 

118 massdict = dict(list(zip(symbols, atoms.get_masses()))) 

119 atomtypes = sorted(massdict.keys()) 

120 atommasses = [massdict[at] for at in atomtypes] 

121 natoms = [symbols.count(at) for at in atomtypes] 

122 ntypes = len(atomtypes) 

123 

124 out.append(str(ntypes)) 

125 out.append(' '.join([str(n) for n in natoms])) 

126 out.append(' '.join([str(n) for n in atommasses])) 

127 

128 atom_id = 0 

129 for n in range(ntypes): 

130 fixed = np.array([False] * len(atoms)) 

131 out.append(atomtypes[n]) 

132 out.append('Coordinates of Component %d' % (n + 1)) 

133 indices = [i for i, at in enumerate(symbols) if at == atomtypes[n]] 

134 a = atoms[indices] 

135 coords = a.positions 

136 for c in a.constraints: 

137 if not isinstance(c, FixAtoms): 

138 warn('Only FixAtoms constraints are supported by con files. ' 

139 'Dropping %r', c) 

140 continue 

141 if c.index.dtype.kind == 'b': 

142 fixed = np.array(c.index, dtype=int) 

143 else: 

144 fixed = np.zeros((natoms[n], ), dtype=int) 

145 for i in c.index: 

146 fixed[i] = 1 

147 for xyz, fix in zip(coords, fixed): 

148 out.append('%22.17f %22.17f %22.17f %d %4d' % 

149 (tuple(xyz) + (fix, atom_id))) 

150 atom_id += 1 

151 fileobj.write('\n'.join(out)) 

152 fileobj.write('\n')