1from ase.units import Bohr 

2 

3 

4def read_turbomole(fd): 

5 """Method to read turbomole coord file 

6 

7 coords in bohr, atom types in lowercase, format: 

8 $coord 

9 x y z atomtype 

10 x y z atomtype f 

11 $end 

12 Above 'f' means a fixed atom. 

13 """ 

14 from ase import Atoms 

15 from ase.constraints import FixAtoms 

16 

17 lines = fd.readlines() 

18 atoms_pos = [] 

19 atom_symbols = [] 

20 myconstraints = [] 

21 

22 # find $coord section; 

23 # does not necessarily have to be the first $<something> in file... 

24 for i, l in enumerate(lines): 

25 if l.strip().startswith('$coord'): 

26 start = i 

27 break 

28 for line in lines[start + 1:]: 

29 if line.startswith('$'): # start of new section 

30 break 

31 else: 

32 x, y, z, symbolraw = line.split()[:4] 

33 symbolshort = symbolraw.strip() 

34 symbol = symbolshort[0].upper() + symbolshort[1:].lower() 

35 # print(symbol) 

36 atom_symbols.append(symbol) 

37 atoms_pos.append( 

38 [float(x) * Bohr, float(y) * Bohr, float(z) * Bohr] 

39 ) 

40 cols = line.split() 

41 if (len(cols) == 5): 

42 fixedstr = line.split()[4].strip() 

43 if (fixedstr == "f"): 

44 myconstraints.append(True) 

45 else: 

46 myconstraints.append(False) 

47 else: 

48 myconstraints.append(False) 

49 

50 # convert Turbomole ghost atom Q to X 

51 atom_symbols = [element if element != 'Q' else 'X' for element in atom_symbols] 

52 atoms = Atoms(positions=atoms_pos, symbols=atom_symbols, pbc=False) 

53 c = FixAtoms(mask=myconstraints) 

54 atoms.set_constraint(c) 

55 return atoms 

56 

57 

58class TurbomoleFormatError(ValueError): 

59 default_message = ('Data format in file does not correspond to known ' 

60 'Turbomole gradient format') 

61 

62 def __init__(self, *args, **kwargs): 

63 if args or kwargs: 

64 ValueError.__init__(self, *args, **kwargs) 

65 else: 

66 ValueError.__init__(self, self.default_message) 

67 

68 

69def read_turbomole_gradient(fd, index=-1): 

70 """ Method to read turbomole gradient file """ 

71 

72 # read entire file 

73 lines = [x.strip() for x in fd.readlines()] 

74 

75 # find $grad section 

76 start = end = -1 

77 for i, line in enumerate(lines): 

78 if not line.startswith('$'): 

79 continue 

80 if line.split()[0] == '$grad': 

81 start = i 

82 elif start >= 0: 

83 end = i 

84 break 

85 

86 if end <= start: 

87 raise RuntimeError('File does not contain a valid \'$grad\' section') 

88 

89 # trim lines to $grad 

90 del lines[:start + 1] 

91 del lines[end - 1 - start:] 

92 

93 # Interpret $grad section 

94 from ase import Atoms, Atom 

95 from ase.calculators.singlepoint import SinglePointCalculator 

96 from ase.units import Bohr, Hartree 

97 images = [] 

98 while lines: # loop over optimization cycles 

99 # header line 

100 # cycle = 1 SCF energy = -267.6666811409 |dE/dxyz| = 0.157112 # noqa: E501 

101 fields = lines[0].split('=') 

102 try: 

103 # cycle = int(fields[1].split()[0]) 

104 energy = float(fields[2].split()[0]) * Hartree 

105 # gradient = float(fields[3].split()[0]) 

106 except (IndexError, ValueError) as e: 

107 raise TurbomoleFormatError() from e 

108 

109 # coordinates/gradient 

110 atoms = Atoms() 

111 forces = [] 

112 for line in lines[1:]: 

113 fields = line.split() 

114 if len(fields) == 4: # coordinates 

115 # 0.00000000000000 0.00000000000000 0.00000000000000 c # noqa: E501 

116 try: 

117 symbol = fields[3].lower().capitalize() 

118 # if dummy atom specified, substitute 'Q' with 'X' 

119 if symbol == 'Q': 

120 symbol = 'X' 

121 position = tuple([Bohr * float(x) for x in fields[0:3]]) 

122 except ValueError as e: 

123 raise TurbomoleFormatError() from e 

124 atoms.append(Atom(symbol, position)) 

125 elif len(fields) == 3: # gradients 

126 # -.51654903354681D-07 -.51654903206651D-07 0.51654903169644D-07 # noqa: E501 

127 grad = [] 

128 for val in fields[:3]: 

129 try: 

130 grad.append( 

131 -float(val.replace('D', 'E')) * Hartree / Bohr 

132 ) 

133 except ValueError as e: 

134 raise TurbomoleFormatError() from e 

135 forces.append(grad) 

136 else: # next cycle 

137 break 

138 

139 # calculator 

140 calc = SinglePointCalculator(atoms, energy=energy, forces=forces) 

141 atoms.calc = calc 

142 

143 # save frame 

144 images.append(atoms) 

145 

146 # delete this frame from data to be handled 

147 del lines[:2 * len(atoms) + 1] 

148 

149 return images[index] 

150 

151 

152def write_turbomole(fd, atoms): 

153 """ Method to write turbomole coord file 

154 """ 

155 from ase.constraints import FixAtoms 

156 

157 coord = atoms.get_positions() 

158 symbols = atoms.get_chemical_symbols() 

159 

160 # convert X to Q for Turbomole ghost atoms 

161 symbols = [element if element != 'X' else 'Q' for element in symbols] 

162 

163 fix_indices = set() 

164 if atoms.constraints: 

165 for constr in atoms.constraints: 

166 if isinstance(constr, FixAtoms): 

167 fix_indices.update(constr.get_indices()) 

168 

169 fix_str = [] 

170 for i in range(len(atoms)): 

171 if i in fix_indices: 

172 fix_str.append('f') 

173 else: 

174 fix_str.append('') 

175 

176 fd.write('$coord\n') 

177 for (x, y, z), s, fix in zip(coord, symbols, fix_str): 

178 fd.write('%20.14f %20.14f %20.14f %2s %2s \n' 

179 % (x / Bohr, y / Bohr, z / Bohr, s.lower(), fix)) 

180 

181 fd.write('$end\n')