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

1# fmt: off 

2 

3from ase.units import Bohr 

4 

5 

6def read_turbomole(fd): 

7 """Method to read turbomole coord file 

8 

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

10 $coord 

11 x y z atomtype 

12 x y z atomtype f 

13 $end 

14 Above 'f' means a fixed atom. 

15 """ 

16 from ase import Atoms 

17 from ase.constraints import FixAtoms 

18 

19 lines = fd.readlines() 

20 atoms_pos = [] 

21 atom_symbols = [] 

22 myconstraints = [] 

23 

24 # find $coord section; 

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

26 for i, l in enumerate(lines): 

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

28 start = i 

29 break 

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

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

32 break 

33 else: 

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

35 symbolshort = symbolraw.strip() 

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

37 # print(symbol) 

38 atom_symbols.append(symbol) 

39 atoms_pos.append( 

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

41 ) 

42 cols = line.split() 

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

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

45 if (fixedstr == "f"): 

46 myconstraints.append(True) 

47 else: 

48 myconstraints.append(False) 

49 else: 

50 myconstraints.append(False) 

51 

52 # convert Turbomole ghost atom Q to X 

53 atom_symbols = [element if element != 

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

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

56 c = FixAtoms(mask=myconstraints) 

57 atoms.set_constraint(c) 

58 return atoms 

59 

60 

61class TurbomoleFormatError(ValueError): 

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

63 'Turbomole gradient format') 

64 

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

66 if args or kwargs: 

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

68 else: 

69 ValueError.__init__(self, self.default_message) 

70 

71 

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

73 """ Method to read turbomole gradient file """ 

74 

75 # read entire file 

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

77 

78 # find $grad section 

79 start = end = -1 

80 for i, line in enumerate(lines): 

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

82 continue 

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

84 start = i 

85 elif start >= 0: 

86 end = i 

87 break 

88 

89 if end <= start: 

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

91 

92 # trim lines to $grad 

93 del lines[:start + 1] 

94 del lines[end - 1 - start:] 

95 

96 # Interpret $grad section 

97 from ase import Atom, Atoms 

98 from ase.calculators.singlepoint import SinglePointCalculator 

99 from ase.units import Bohr, Hartree 

100 images = [] 

101 while lines: # loop over optimization cycles 

102 # header line 

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

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

105 try: 

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

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

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

109 except (IndexError, ValueError) as e: 

110 raise TurbomoleFormatError from e 

111 

112 # coordinates/gradient 

113 atoms = Atoms() 

114 forces = [] 

115 for line in lines[1:]: 

116 fields = line.split() 

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

118 # 0.00000000000000 0.00000000000000 0.00000000000000 c # noqa: E501 

119 try: 

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

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

122 if symbol == 'Q': 

123 symbol = 'X' 

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

125 except ValueError as e: 

126 raise TurbomoleFormatError from e 

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

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

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

130 grad = [] 

131 for val in fields[:3]: 

132 try: 

133 grad.append( 

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

135 ) 

136 except ValueError as e: 

137 raise TurbomoleFormatError from e 

138 forces.append(grad) 

139 else: # next cycle 

140 break 

141 

142 # calculator 

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

144 atoms.calc = calc 

145 

146 # save frame 

147 images.append(atoms) 

148 

149 # delete this frame from data to be handled 

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

151 

152 return images[index] 

153 

154 

155def write_turbomole(fd, atoms): 

156 """ Method to write turbomole coord file 

157 """ 

158 from ase.constraints import FixAtoms 

159 

160 coord = atoms.get_positions() 

161 symbols = atoms.get_chemical_symbols() 

162 

163 # convert X to Q for Turbomole ghost atoms 

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

165 

166 fix_indices = set() 

167 if atoms.constraints: 

168 for constr in atoms.constraints: 

169 if isinstance(constr, FixAtoms): 

170 fix_indices.update(constr.get_indices()) 

171 

172 fix_str = [] 

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

174 if i in fix_indices: 

175 fix_str.append('f') 

176 else: 

177 fix_str.append('') 

178 

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

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

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

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

183 

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