Coverage for /builds/kinetik161/ase/ase/io/turbomole.py: 49.54%

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 != 

52 'Q' else 'X' for element in atom_symbols] 

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

54 c = FixAtoms(mask=myconstraints) 

55 atoms.set_constraint(c) 

56 return atoms 

57 

58 

59class TurbomoleFormatError(ValueError): 

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

61 'Turbomole gradient format') 

62 

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

64 if args or kwargs: 

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

66 else: 

67 ValueError.__init__(self, self.default_message) 

68 

69 

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

71 """ Method to read turbomole gradient file """ 

72 

73 # read entire file 

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

75 

76 # find $grad section 

77 start = end = -1 

78 for i, line in enumerate(lines): 

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

80 continue 

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

82 start = i 

83 elif start >= 0: 

84 end = i 

85 break 

86 

87 if end <= start: 

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

89 

90 # trim lines to $grad 

91 del lines[:start + 1] 

92 del lines[end - 1 - start:] 

93 

94 # Interpret $grad section 

95 from ase import Atom, Atoms 

96 from ase.calculators.singlepoint import SinglePointCalculator 

97 from ase.units import Bohr, Hartree 

98 images = [] 

99 while lines: # loop over optimization cycles 

100 # header line 

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

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

103 try: 

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

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

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

107 except (IndexError, ValueError) as e: 

108 raise TurbomoleFormatError() from e 

109 

110 # coordinates/gradient 

111 atoms = Atoms() 

112 forces = [] 

113 for line in lines[1:]: 

114 fields = line.split() 

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

116 # 0.00000000000000 0.00000000000000 0.00000000000000 c # noqa: E501 

117 try: 

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

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

120 if symbol == 'Q': 

121 symbol = 'X' 

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

123 except ValueError as e: 

124 raise TurbomoleFormatError() from e 

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

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

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

128 grad = [] 

129 for val in fields[:3]: 

130 try: 

131 grad.append( 

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

133 ) 

134 except ValueError as e: 

135 raise TurbomoleFormatError() from e 

136 forces.append(grad) 

137 else: # next cycle 

138 break 

139 

140 # calculator 

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

142 atoms.calc = calc 

143 

144 # save frame 

145 images.append(atoms) 

146 

147 # delete this frame from data to be handled 

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

149 

150 return images[index] 

151 

152 

153def write_turbomole(fd, atoms): 

154 """ Method to write turbomole coord file 

155 """ 

156 from ase.constraints import FixAtoms 

157 

158 coord = atoms.get_positions() 

159 symbols = atoms.get_chemical_symbols() 

160 

161 # convert X to Q for Turbomole ghost atoms 

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

163 

164 fix_indices = set() 

165 if atoms.constraints: 

166 for constr in atoms.constraints: 

167 if isinstance(constr, FixAtoms): 

168 fix_indices.update(constr.get_indices()) 

169 

170 fix_str = [] 

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

172 if i in fix_indices: 

173 fix_str.append('f') 

174 else: 

175 fix_str.append('') 

176 

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

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

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

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

181 

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