Coverage for /builds/debichem-team/python-ase/ase/io/orca.py: 92.31%

1import re 

2from io import StringIO 

3from pathlib import Path 

4from typing import List, Optional 

5 

6import numpy as np 

7 

8from ase.io import read 

9from ase.units import Bohr, Hartree 

10from ase.utils import reader, writer 

11 

12# Made from NWChem interface 

13 

14 

15@reader 

16def read_geom_orcainp(fd): 

17 """Method to read geometry from an ORCA input file.""" 

18 lines = fd.readlines() 

19 

20 # Find geometry region of input file. 

21 stopline = 0 

22 for index, line in enumerate(lines): 

23 if line[1:].startswith('xyz '): 

24 startline = index + 1 

25 stopline = -1 

26 elif (line.startswith('end') and stopline == -1): 

27 stopline = index 

28 elif (line.startswith('*') and stopline == -1): 

29 stopline = index 

30 # Format and send to read_xyz. 

31 xyz_text = '%i\n' % (stopline - startline) 

32 xyz_text += ' geometry\n' 

33 for line in lines[startline:stopline]: 

34 xyz_text += line 

35 atoms = read(StringIO(xyz_text), format='xyz') 

36 atoms.set_cell((0., 0., 0.)) # no unit cell defined 

37 

38 return atoms 

39 

40 

41@writer 

42def write_orca(fd, atoms, params): 

43 # conventional filename: '<name>.inp' 

44 fd.write(f"! {params['orcasimpleinput']} \n") 

45 fd.write(f"{params['orcablocks']} \n") 

46 

47 if 'coords' not in params['orcablocks']: 

48 fd.write('*xyz') 

49 fd.write(" %d" % params['charge']) 

50 fd.write(" %d \n" % params['mult']) 

51 for atom in atoms: 

52 if atom.tag == 71: # 71 is ascii G (Ghost) 

53 symbol = atom.symbol + ' : ' 

54 else: 

55 symbol = atom.symbol + ' ' 

56 fd.write( 

57 symbol 

58 + str(atom.position[0]) 

59 + " " 

60 + str(atom.position[1]) 

61 + " " 

62 + str(atom.position[2]) 

63 + "\n" 

64 ) 

65 fd.write('*\n') 

66 

67 

68def read_charge(lines: List[str]) -> Optional[float]: 

69 """Read sum of atomic charges.""" 

70 charge = None 

71 for line in lines: 

72 if 'Sum of atomic charges' in line: 

73 charge = float(line.split()[-1]) 

74 return charge 

75 

76 

77def read_energy(lines: List[str]) -> Optional[float]: 

78 """Read energy.""" 

79 energy = None 

80 for line in lines: 

81 if 'FINAL SINGLE POINT ENERGY' in line: 

82 if "Wavefunction not fully converged" in line: 

83 energy = float('nan') 

84 else: 

85 energy = float(line.split()[-1]) 

86 if energy is not None: 

87 return energy * Hartree 

88 return energy 

89 

90 

91def read_center_of_mass(lines: List[str]) -> Optional[np.ndarray]: 

92 """ Scan through text for the center of mass """ 

93 # Example: 

94 # 'The origin for moment calculation is the CENTER OF MASS = 

95 # ( 0.002150, -0.296255 0.086315)' 

96 # Note the missing comma in the output 

97 com = None 

98 for line in lines: 

99 if 'The origin for moment calculation is the CENTER OF MASS' in line: 

100 line = re.sub(r'[(),]', '', line) 

101 com = np.array([float(_) for _ in line.split()[-3:]]) 

102 if com is not None: 

103 return com * Bohr # return the last match 

104 return com 

105 

106 

107def read_dipole(lines: List[str]) -> Optional[np.ndarray]: 

108 """Read dipole moment. 

109 

110 Note that the read dipole moment is for the COM frame of reference. 

111 """ 

112 dipole = None 

113 for line in lines: 

114 if 'Total Dipole Moment' in line: 

115 dipole = np.array([float(_) for _ in line.split()[-3:]]) 

116 if dipole is not None: 

117 return dipole * Bohr # Return the last match 

118 return dipole 

119 

120 

121@reader 

122def read_orca_output(fd): 

123 """ From the ORCA output file: Read Energy and dipole moment 

124 in the frame of reference of the center of mass " 

125 """ 

126 lines = fd.readlines() 

127 

128 energy = read_energy(lines) 

129 charge = read_charge(lines) 

130 com = read_center_of_mass(lines) 

131 dipole = read_dipole(lines) 

132 

133 results = {} 

134 results['energy'] = energy 

135 results['free_energy'] = energy 

136 

137 if com is not None and dipole is not None: 

138 dipole = dipole + com * charge 

139 results['dipole'] = dipole 

140 

141 return results 

142 

143 

144@reader 

145def read_orca_engrad(fd): 

146 """Read Forces from ORCA .engrad file.""" 

147 getgrad = False 

148 gradients = [] 

149 tempgrad = [] 

150 for _, line in enumerate(fd): 

151 if line.find('# The current gradient') >= 0: 

152 getgrad = True 

153 gradients = [] 

154 tempgrad = [] 

155 continue 

156 if getgrad and "#" not in line: 

157 grad = line.split()[-1] 

158 tempgrad.append(float(grad)) 

159 if len(tempgrad) == 3: 

160 gradients.append(tempgrad) 

161 tempgrad = [] 

162 if '# The at' in line: 

163 getgrad = False 

164 

165 forces = -np.array(gradients) * Hartree / Bohr 

166 return forces 

167 

168 

169def read_orca_outputs(directory, stdout_path): 

170 stdout_path = Path(stdout_path) 

171 results = {} 

172 results.update(read_orca_output(stdout_path)) 

173 

174 # Does engrad always exist? - No! 

175 # Will there be other files -No -> We should just take engrad 

176 # as a direct argument. Or maybe this function does not even need to 

177 # exist. 

178 engrad_path = stdout_path.with_suffix('.engrad') 

179 if engrad_path.is_file(): 

180 results['forces'] = read_orca_engrad(engrad_path) 

181 return results