Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

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='engrad 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 Should in most cases contain "engrad" or method that 

39 writes the engrad file. If not (single point only), 

40 set the "task" parameter manually. 

41 Default is ``engrad tightscf PBE def2-SVP``. 

42 

43 orcablocks: str 

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

45 Default is ``%scf maxiter 200 end``. 

46 

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

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

49 available in ORCA instead of hardcoding stuff. 

50 

51 Default parameters are: 

52 

53 charge: 0 

54 

55 mult: 1 

56 

57 task: 'gradient' 

58 

59 Point Charge IO functionality added by A. Dohn. 

60 """ 

61 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file, 

62 label, atoms, **kwargs) 

63 

64 self.pcpot = None 

65 

66 def set(self, **kwargs): 

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

68 if changed_parameters: 

69 self.reset() 

70 

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

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

73 p = self.parameters 

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

75 p['label'] = self.label 

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

77 p['pcpot'] = self.pcpot 

78 

79 write_orca(atoms, **p) 

80 

81 def read(self, label): 

82 FileIOCalculator.read(self, label) 

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

84 raise ReadError 

85 

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

87 for line in fd: 

88 if line.startswith('geometry'): 

89 break 

90 symbols = [] 

91 positions = [] 

92 for line in fd: 

93 if line.startswith('end'): 

94 break 

95 words = line.split() 

96 symbols.append(words[0]) 

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

98 

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

100 self.read_results() 

101 

102 def read_results(self): 

103 self.read_energy() 

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

105 self.read_forces() 

106 

107 def read_energy(self): 

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

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

110 text = fd.read() 

111 # Energy: 

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

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

114 found_line = re_energy.finditer(text) 

115 for match in found_line: 

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

117 self.results['energy'] = float( 

118 match.group().split()[-1]) * Hartree 

119 

120 def read_forces(self): 

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

122 if not os.path.isfile(self.label + '.engrad'): 

123 raise ReadError("Engrad file missing.") 

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

125 lines = fd.readlines() 

126 getgrad = False 

127 gradients = [] 

128 tempgrad = [] 

129 for i, line in enumerate(lines): 

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

131 getgrad = True 

132 gradients = [] 

133 tempgrad = [] 

134 continue 

135 if getgrad and "#" not in line: 

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

137 tempgrad.append(float(grad)) 

138 if len(tempgrad) == 3: 

139 gradients.append(tempgrad) 

140 tempgrad = [] 

141 if '# The at' in line: 

142 getgrad = False 

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

144 

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

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

147 """ 

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

149 return self.pcpot 

150 

151 

152class PointChargePotential: 

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

154 """ Point Charge Potential Interface to ORCA """ 

155 if positions is not None: 

156 self.set_positions(positions) 

157 if directory is None: 

158 directory = os.getcwd() 

159 

160 self.directory = directory + os.sep 

161 self.mmcharges = mmcharges 

162 self.label = label 

163 

164 def set_positions(self, positions): 

165 self.positions = positions 

166 

167 def set_charges(self, mmcharges): 

168 self.q_p = mmcharges 

169 

170 def write_mmcharges(self, filename): 

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

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

173 

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

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

176 [x, y, z] = pos 

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

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

179 

180 pc_file.close() 

181 

182 def get_forces(self, calc): 

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

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

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

186 lines = fd.readlines() 

187 numpc = int(lines[0]) 

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

189 for i in range(numpc): 

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

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

192 

193 return -forces * Hartree / Bohr