Coverage for /builds/kinetik161/ase/ase/io/orca.py: 96.39%

1import os 

2import re 

3from io import StringIO 

4from pathlib import Path 

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 fd.write('*xyz') 

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

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

50 for atom in atoms: 

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

52 symbol = atom.symbol + ' : ' 

53 else: 

54 symbol = atom.symbol + ' ' 

55 fd.write(symbol + 

56 str(atom.position[0]) + ' ' + 

57 str(atom.position[1]) + ' ' + 

58 str(atom.position[2]) + '\n') 

59 fd.write('*\n') 

60 

61 

62@reader 

63def read_orca_energy(fd): 

64 """Read Energy from ORCA output file.""" 

65 text = fd.read() 

66 re_energy = re.compile(r"FINAL SINGLE POINT ENERGY.*\n") 

67 re_not_converged = re.compile(r"Wavefunction not fully converged") 

68 

69 found_line = re_energy.finditer(text) 

70 energy = float('nan') 

71 for match in found_line: 

72 if not re_not_converged.search(match.group()): 

73 energy = float(match.group().split()[-1]) * Hartree 

74 if np.isnan(energy): 

75 raise RuntimeError('No energy') 

76 else: 

77 return energy 

78 

79 

80@reader 

81def read_orca_forces(fd): 

82 """Read Forces from ORCA output file.""" 

83 getgrad = False 

84 gradients = [] 

85 tempgrad = [] 

86 for i, line in enumerate(fd): 

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

88 getgrad = True 

89 gradients = [] 

90 tempgrad = [] 

91 continue 

92 if getgrad and "#" not in line: 

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

94 tempgrad.append(float(grad)) 

95 if len(tempgrad) == 3: 

96 gradients.append(tempgrad) 

97 tempgrad = [] 

98 if '# The at' in line: 

99 getgrad = False 

100 

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

102 return forces 

103 

104 

105def read_orca_outputs(directory, stdout_path): 

106 results = {} 

107 energy = read_orca_energy(Path(stdout_path)) 

108 results['energy'] = energy 

109 results['free_energy'] = energy 

110 

111 # Does engrad always exist? - No! 

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

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

114 # exist. 

115 engrad_path = Path(stdout_path).with_suffix('.engrad') 

116 if os.path.isfile(engrad_path): 

117 results['forces'] = read_orca_forces(engrad_path) 

118 return results