Coverage for /builds/ase/ase/ase/calculators/qchem.py: 54.22%

1# fmt: off 

2 

3import numpy as np 

4 

5import ase.units 

6from ase.calculators.calculator import FileIOCalculator, SCFError 

7 

8 

9class QChem(FileIOCalculator): 

10 """ 

11 QChem calculator 

12 """ 

13 name = 'QChem' 

14 

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

16 _legacy_default_command = 'qchem PREFIX.inp PREFIX.out' 

17 

18 # Following the minimal requirements given in 

19 # http://www.q-chem.com/qchem-website/manual/qchem43_manual/sect-METHOD.html 

20 default_parameters = {'method': 'hf', 

21 'basis': '6-31G*', 

22 'jobtype': None, 

23 'charge': 0} 

24 

25 def __init__(self, restart=None, 

26 ignore_bad_restart_file=FileIOCalculator._deprecated, 

27 label='qchem', scratch=None, np=1, nt=1, pbs=False, 

28 basisfile=None, ecpfile=None, atoms=None, **kwargs): 

29 """ 

30 The scratch directory, number of processor and threads as well as a few 

31 other command line options can be set using the arguments explained 

32 below. The remaining kwargs are copied as options to the input file. 

33 The calculator will convert these options to upper case 

34 (Q-Chem standard) when writing the input file. 

35 

36 scratch: str 

37 path of the scratch directory 

38 np: int 

39 number of processors for the -np command line flag 

40 nt: int 

41 number of threads for the -nt command line flag 

42 pbs: boolean 

43 command line flag for pbs scheduler (see Q-Chem manual) 

44 basisfile: str 

45 path to file containing the basis. Use in combination with 

46 basis='gen' keyword argument. 

47 ecpfile: str 

48 path to file containing the effective core potential. Use in 

49 combination with ecp='gen' keyword argument. 

50 """ 

51 

52 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file, 

53 label, atoms, **kwargs) 

54 

55 # Augment the command by various flags 

56 if pbs: 

57 self.command = 'qchem -pbs ' 

58 else: 

59 self.command = 'qchem ' 

60 if np != 1: 

61 self.command += '-np %d ' % np 

62 if nt != 1: 

63 self.command += '-nt %d ' % nt 

64 self.command += 'PREFIX.inp PREFIX.out' 

65 if scratch is not None: 

66 self.command += f' {scratch}' 

67 

68 self.basisfile = basisfile 

69 self.ecpfile = ecpfile 

70 

71 def read(self, label): 

72 raise NotImplementedError 

73 

74 def read_results(self): 

75 filename = self.label + '.out' 

76 

77 with open(filename) as fileobj: 

78 lineiter = iter(fileobj) 

79 for line in lineiter: 

80 if 'SCF failed to converge' in line: 

81 raise SCFError() 

82 elif 'ERROR: alpha_min' in line: 

83 # Even though it is not technically a SCFError: 

84 raise SCFError() 

85 elif ' Total energy in the final basis set =' in line: 

86 convert = ase.units.Hartree 

87 self.results['energy'] = float(line.split()[8]) * convert 

88 elif ' Gradient of SCF Energy' in line: 

89 # Read gradient as 3 by N array and transpose at the end 

90 gradient = [[] for _ in range(3)] 

91 # Skip first line containing atom numbering 

92 next(lineiter) 

93 while True: 

94 # Loop over the three Cartesian coordinates 

95 for i in range(3): 

96 # Cut off the component numbering and remove 

97 # trailing characters ('\n' and stuff) 

98 line = next(lineiter)[5:].rstrip() 

99 # Cut in chunks of 12 symbols and convert into 

100 # strings. This is preferred over string.split() as 

101 # the fields may overlap for large gradients 

102 gradient[i].extend(list(map( 

103 float, [line[i:i + 12] 

104 for i in range(0, len(line), 12)]))) 

105 

106 # After three force components we expect either a 

107 # separator line, which we want to skip, or the end of 

108 # the gradient matrix which is characterized by the 

109 # line ' Max gradient component'. 

110 # Maybe change stopping criterion to be independent of 

111 # next line. Eg. if not lineiter.next().startswith(' ') 

112 if ' Max gradient component' in next(lineiter): 

113 # Minus to convert from gradient to force 

114 self.results['forces'] = np.array(gradient).T * ( 

115 -ase.units.Hartree / ase.units.Bohr) 

116 break 

117 

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

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

120 filename = self.label + '.inp' 

121 

122 with open(filename, 'w') as fileobj: 

123 fileobj.write('$comment\n ASE generated input file\n$end\n\n') 

124 

125 fileobj.write('$rem\n') 

126 if self.parameters['jobtype'] is None: 

127 if 'forces' in properties: 

128 fileobj.write(' %-25s %s\n' % ('JOBTYPE', 'FORCE')) 

129 else: 

130 fileobj.write(' %-25s %s\n' % ('JOBTYPE', 'SP')) 

131 

132 for prm in self.parameters: 

133 if prm not in ['charge', 'multiplicity']: 

134 if self.parameters[prm] is not None: 

135 fileobj.write(' %-25s %s\n' % ( 

136 prm.upper(), self.parameters[prm].upper())) 

137 

138 # Not even a parameters as this is an absolute necessity 

139 fileobj.write(' %-25s %s\n' % ('SYM_IGNORE', 'TRUE')) 

140 fileobj.write('$end\n\n') 

141 

142 fileobj.write('$molecule\n') 

143 # Following the example set by the gaussian calculator 

144 if ('multiplicity' not in self.parameters): 

145 tot_magmom = atoms.get_initial_magnetic_moments().sum() 

146 mult = tot_magmom + 1 

147 else: 

148 mult = self.parameters['multiplicity'] 

149 # Default charge of 0 is defined in default_parameters 

150 fileobj.write(' %d %d\n' % (self.parameters['charge'], mult)) 

151 for a in atoms: 

152 fileobj.write(' {} {:f} {:f} {:f}\n'.format(a.symbol, 

153 a.x, a.y, a.z)) 

154 fileobj.write('$end\n\n') 

155 

156 if self.basisfile is not None: 

157 with open(self.basisfile) as f_in: 

158 basis = f_in.readlines() 

159 fileobj.write('$basis\n') 

160 fileobj.writelines(basis) 

161 fileobj.write('$end\n\n') 

162 

163 if self.ecpfile is not None: 

164 with open(self.ecpfile) as f_in: 

165 ecp = f_in.readlines() 

166 fileobj.write('$ecp\n') 

167 fileobj.writelines(ecp) 

168 fileobj.write('$end\n\n')