1import os 

2import re 

3 

4import numpy as np 

5 

6from ase.units import Hartree, Bohr 

7from ase.io.orca import write_orca 

8from ase.calculators.calculator import FileIOCalculator, Parameters, ReadError 

9 

10 

11class ORCA(FileIOCalculator): 

12 implemented_properties = ['energy', 'forces'] 

13 

14 if 'ORCA_COMMAND' in os.environ: 

15 command = os.environ['ORCA_COMMAND'] + ' PREFIX.inp > PREFIX.out' 

16 else: 

17 command = 'orca PREFIX.inp > PREFIX.out' 

18 

19 default_parameters = dict( 

20 charge=0, mult=1, 

21 task='gradient', 

22 orcasimpleinput='tightscf PBE def2-SVP', 

23 orcablocks='%scf maxiter 200 end') 

24 

25 def __init__(self, restart=None, 

26 ignore_bad_restart_file=FileIOCalculator._deprecated, 

27 label='orca', atoms=None, **kwargs): 

28 """ ASE interface to ORCA 4 

29 by Ragnar Bjornsson, Based on NWchem interface but simplified. 

30 Only supports energies and gradients (no dipole moments, 

31 orbital energies etc.) for now. 

32 

33 For more ORCA-keyword flexibility, method/xc/basis etc. 

34 keywords are not used. Instead, two keywords: 

35 

36 orcasimpleinput: str 

37 What you'd put after the "!" in an orca input file. 

38 

39 orcablock: str 

40 What you'd put in the "% ... end"-blocks. 

41 

42 are used to define the ORCA simple-inputline and the ORCA-block input. 

43 This allows for more flexible use of any ORCA method or keyword 

44 available in ORCA instead of hardcoding stuff. 

45 

46 Point Charge IO functionality added by A. Dohn. 

47 """ 

48 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file, 

49 label, atoms, **kwargs) 

50 

51 self.pcpot = None 

52 

53 def set(self, **kwargs): 

54 changed_parameters = FileIOCalculator.set(self, **kwargs) 

55 if changed_parameters: 

56 self.reset() 

57 

58 def write_input(self, atoms, properties=None, system_changes=None): 

59 FileIOCalculator.write_input(self, atoms, properties, system_changes) 

60 p = self.parameters 

61 p.write(self.label + '.ase') 

62 p['label'] = self.label 

63 if self.pcpot: # also write point charge file and add things to input 

64 p['pcpot'] = self.pcpot 

65 

66 write_orca(atoms, **p) 

67 

68 def read(self, label): 

69 FileIOCalculator.read(self, label) 

70 if not os.path.isfile(self.label + '.out'): 

71 raise ReadError 

72 

73 with open(self.label + '.inp') as fd: 

74 for line in fd: 

75 if line.startswith('geometry'): 

76 break 

77 symbols = [] 

78 positions = [] 

79 for line in fd: 

80 if line.startswith('end'): 

81 break 

82 words = line.split() 

83 symbols.append(words[0]) 

84 positions.append([float(word) for word in words[1:]]) 

85 

86 self.parameters = Parameters.read(self.label + '.ase') 

87 self.read_results() 

88 

89 def read_results(self): 

90 self.read_energy() 

91 if self.parameters.task.find('gradient') > -1: 

92 self.read_forces() 

93 

94 def read_energy(self): 

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

96 with open(self.label + '.out', mode='r', encoding='utf-8') as fd: 

97 text = fd.read() 

98 # Energy: 

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

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

101 found_line = re_energy.search(text) 

102 if found_line and not re_not_converged.search(found_line.group()): 

103 self.results['energy'] = float(found_line.group().split()[-1]) * Hartree 

104 

105 def read_forces(self): 

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

107 with open(f'{self.label}.engrad', 'r') as fd: 

108 lines = fd.readlines() 

109 getgrad = False 

110 gradients = [] 

111 tempgrad = [] 

112 for i, line in enumerate(lines): 

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

114 getgrad = True 

115 gradients = [] 

116 tempgrad = [] 

117 continue 

118 if getgrad and "#" not in line: 

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

120 tempgrad.append(float(grad)) 

121 if len(tempgrad) == 3: 

122 gradients.append(tempgrad) 

123 tempgrad = [] 

124 if '# The at' in line: 

125 getgrad = False 

126 self.results['forces'] = -np.array(gradients) * Hartree / Bohr 

127 

128 def embed(self, mmcharges=None, **parameters): 

129 """Embed atoms in point-charges (mmcharges) 

130 """ 

131 self.pcpot = PointChargePotential(mmcharges, label=self.label) 

132 return self.pcpot 

133 

134 

135class PointChargePotential: 

136 def __init__(self, mmcharges, label=None, positions=None, directory=None): 

137 """ Point Charge Potential Interface to ORCA """ 

138 if positions is not None: 

139 self.set_positions(positions) 

140 if directory is None: 

141 directory = os.getcwd() 

142 

143 self.directory = directory + os.sep 

144 self.mmcharges = mmcharges 

145 self.label = label 

146 

147 def set_positions(self, positions): 

148 self.positions = positions 

149 

150 def set_charges(self, mmcharges): 

151 self.q_p = mmcharges 

152 

153 def write_mmcharges(self, filename): 

154 pc_file = open(os.path.join(self.directory, 

155 filename + '.pc'), 'w') 

156 

157 pc_file.write('{0:d}\n'.format(len(self.mmcharges))) 

158 for [pos, pc] in zip(self.positions, self.mmcharges): 

159 [x, y, z] = pos 

160 pc_file.write('{0:12.6f} {1:12.6f} {2:12.6f} {3:12.6f}\n' 

161 .format(pc, x, y, z)) 

162 

163 pc_file.close() 

164 

165 def get_forces(self, calc): 

166 ''' reads forces on point charges from .pcgrad file ''' 

167 with open(os.path.join(self.directory, self.label + '.pcgrad'), 

168 'r', encoding='utf-8') as fd: 

169 lines = fd.readlines() 

170 numpc = int(lines[0]) 

171 forces = np.zeros((numpc, 3)) 

172 for i in range(numpc): 

173 [fx, fy, fz] = [float(f) for f in lines[i + 1].split()] 

174 forces[i, :] = fx, fy, fz 

175 

176 return -forces * Hartree / Bohr