Coverage for /builds/debichem-team/python-ase/ase/calculators/qchem.py: 54.22%

1import numpy as np 

2 

3import ase.units 

4from ase.calculators.calculator import FileIOCalculator, SCFError 

5 

6 

7class QChem(FileIOCalculator): 

8 """ 

9 QChem calculator 

10 """ 

11 name = 'QChem' 

12 

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

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

15 

16 # Following the minimal requirements given in 

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

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

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

20 'jobtype': None, 

21 'charge': 0} 

22 

23 def __init__(self, restart=None, 

24 ignore_bad_restart_file=FileIOCalculator._deprecated, 

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

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

27 """ 

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

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

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

31 The calculator will convert these options to upper case 

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

33 

34 scratch: str 

35 path of the scratch directory 

36 np: int 

37 number of processors for the -np command line flag 

38 nt: int 

39 number of threads for the -nt command line flag 

40 pbs: boolean 

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

42 basisfile: str 

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

44 basis='gen' keyword argument. 

45 ecpfile: str 

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

47 combination with ecp='gen' keyword argument. 

48 """ 

49 

50 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file, 

51 label, atoms, **kwargs) 

52 

53 # Augment the command by various flags 

54 if pbs: 

55 self.command = 'qchem -pbs ' 

56 else: 

57 self.command = 'qchem ' 

58 if np != 1: 

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

60 if nt != 1: 

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

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

63 if scratch is not None: 

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

65 

66 self.basisfile = basisfile 

67 self.ecpfile = ecpfile 

68 

69 def read(self, label): 

70 raise NotImplementedError 

71 

72 def read_results(self): 

73 filename = self.label + '.out' 

74 

75 with open(filename) as fileobj: 

76 lineiter = iter(fileobj) 

77 for line in lineiter: 

78 if 'SCF failed to converge' in line: 

79 raise SCFError() 

80 elif 'ERROR: alpha_min' in line: 

81 # Even though it is not technically a SCFError: 

82 raise SCFError() 

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

84 convert = ase.units.Hartree 

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

86 elif ' Gradient of SCF Energy' in line: 

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

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

89 # Skip first line containing atom numbering 

90 next(lineiter) 

91 while True: 

92 # Loop over the three Cartesian coordinates 

93 for i in range(3): 

94 # Cut off the component numbering and remove 

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

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

97 # Cut in chunks of 12 symbols and convert into 

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

99 # the fields may overlap for large gradients 

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

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

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

103 

104 # After three force components we expect either a 

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

106 # the gradient matrix which is characterized by the 

107 # line ' Max gradient component'. 

108 # Maybe change stopping criterion to be independent of 

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

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

111 # Minus to convert from gradient to force 

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

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

114 break 

115 

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

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

118 filename = self.label + '.inp' 

119 

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

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

122 

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

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

125 if 'forces' in properties: 

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

127 else: 

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

129 

130 for prm in self.parameters: 

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

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

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

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

135 

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

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

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

139 

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

141 # Following the example set by the gaussian calculator 

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

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

144 mult = tot_magmom + 1 

145 else: 

146 mult = self.parameters['multiplicity'] 

147 # Default charge of 0 is defined in default_parameters 

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

149 for a in atoms: 

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

151 a.x, a.y, a.z)) 

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

153 

154 if self.basisfile is not None: 

155 with open(self.basisfile) as f_in: 

156 basis = f_in.readlines() 

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

158 fileobj.writelines(basis) 

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

160 

161 if self.ecpfile is not None: 

162 with open(self.ecpfile) as f_in: 

163 ecp = f_in.readlines() 

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

165 fileobj.writelines(ecp) 

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