Coverage for /builds/kinetik161/ase/ase/calculators/nwchem.py: 68.42%

1"""This module defines an ASE interface to NWchem 

2 

3https://nwchemgit.github.io 

4""" 

5import os 

6 

7import numpy as np 

8 

9from ase import io 

10from ase.calculators.calculator import FileIOCalculator 

11from ase.spectrum.band_structure import BandStructure 

12from ase.units import Hartree 

13 

14 

15class NWChem(FileIOCalculator): 

16 implemented_properties = ['energy', 'free_energy', 

17 'forces', 'stress', 'dipole'] 

18 _legacy_default_command = 'nwchem PREFIX.nwi > PREFIX.nwo' 

19 accepts_bandpath_keyword = True 

20 discard_results_on_any_change = True 

21 

22 def __init__(self, restart=None, 

23 ignore_bad_restart_file=FileIOCalculator._deprecated, 

24 label='nwchem', atoms=None, command=None, **kwargs): 

25 """ 

26 NWChem keywords are specified using (potentially nested) 

27 dictionaries. Consider the following input file block:: 

28 

29 dft 

30 odft 

31 mult 2 

32 convergence energy 1e-9 density 1e-7 gradient 5e-6 

33 end 

34 

35 This can be generated by the NWChem calculator by using the 

36 following settings: 

37 >>> from ase.calculators.nwchem import NWChem 

38 >>> calc = NWChem(dft={'odft': None, 

39 ... 'mult': 2, 

40 ... 'convergence': {'energy': 1e-9, 

41 ... 'density': 1e-7, 

42 ... 'gradient': 5e-6, 

43 ... }, 

44 ... }, 

45 ... ) 

46 

47 In addition, the calculator supports several special keywords: 

48 

49 theory: str 

50 Which NWChem module should be used to calculate the 

51 energies and forces. Supported values are ``'dft'``, 

52 ``'scf'``, ``'mp2'``, ``'ccsd'``, ``'tce'``, ``'tddft'``, 

53 ``'pspw'``, ``'band'``, and ``'paw'``. If not provided, the 

54 calculator will attempt to guess which theory to use based 

55 on the keywords provided by the user. 

56 xc: str 

57 The exchange-correlation functional to use. Only relevant 

58 for DFT calculations. 

59 task: str 

60 What type of calculation is to be performed, e.g. 

61 ``'energy'``, ``'gradient'``, ``'optimize'``, etc. When 

62 using ``'SocketIOCalculator'``, ``task`` should be set 

63 to ``'optimize'``. In most other circumstances, ``task`` 

64 should not be set manually. 

65 basis: str or dict 

66 Which basis set to use for gaussian-type orbital 

67 calculations. Set to a string to use the same basis for all 

68 elements. To use a different basis for different elements, 

69 provide a dict of the form: 

70 

71 >>> calc = NWChem(..., 

72 ... basis={'O': '3-21G', 

73 ... 'Si': '6-31g'}) 

74 

75 basispar: str 

76 Additional keywords to put in the NWChem ``basis`` block, 

77 e.g. ``'rel'`` for relativistic bases. 

78 symmetry: int or str 

79 The point group (for gaussian-type orbital calculations) or 

80 space group (for plane-wave calculations) of the system. 

81 Supports both group names (e.g. ``'c2v'``, ``'Fm3m'``) and 

82 numbers (e.g. ``225``). 

83 autosym: bool 

84 Whether NWChem should automatically determine the symmetry 

85 of the structure (defaults to ``False``). 

86 center: bool 

87 Whether NWChem should automatically center the structure 

88 (defaults to ``False``). Enable at your own risk. 

89 autoz: bool 

90 Whether NWChem should automatically construct a Z-matrix 

91 for your molecular system (defaults to ``False``). 

92 geompar: str 

93 Additional keywords to put in the NWChem `geometry` block, 

94 e.g. ``'nucleus finite'`` for gaussian-shaped nuclear 

95 charges. Do not set ``'autosym'``, ``'center'``, or 

96 ``'autoz'`` in this way; instead, use the appropriate 

97 keyword described above for these settings. 

98 set: dict 

99 Used to manually create or modify entries in the NWChem 

100 rtdb. For example, the following settings enable 

101 pseudopotential filtering for plane-wave calculations:: 

102 

103 set nwpw:kbpp_ray .true. 

104 set nwpw:kbpp_filter .true. 

105 

106 These settings are generated by the NWChem calculator by 

107 passing the arguments: 

108 

109 >>> calc = NWChem(..., 

110 >>> set={'nwpw:kbpp_ray': True, 

111 >>> 'nwpw:kbpp_filter': True}) 

112 

113 kpts: (int, int, int), or dict 

114 Indicates which k-point mesh to use. Supported syntax is 

115 similar to that of GPAW. Implies ``theory='band'``. 

116 bandpath: BandPath object 

117 The band path to use for a band structure calculation. 

118 Implies ``theory='band'``. 

119 pretasks: list of dict 

120 Tasks used to produce a better initial guess 

121 for the wavefunction. 

122 These task typically use a cheaper level of theory 

123 or smaller basis set (but not both). 

124 The output energy and forces should remain unchanged 

125 regardless of the number of tasks or their parameters, 

126 but the runtime may be significantly improved. 

127 

128 For example, a MP2 calculation preceded by guesses at the 

129 DFT and HF levels would be 

130 

131 >>> calc = NWChem(theory='mp2', basis='aug-cc-pvdz', 

132 >>> pretasks=[ 

133 >>> {'dft': {'xc': 'hfexch'}, 

134 >>> 'set': {'lindep:n_dep': 0}}, 

135 >>> {'theory': 'scf', 'set': {'lindep:n_dep': 0}}, 

136 >>> ]) 

137 

138 Each dictionary could contain any of the other parameters, 

139 except those which pertain to global configurations 

140 (e.g., geometry details, scratch dir). 

141 

142 The default basis set is that of the final step in the calculation, 

143 or that of the previous step that which defines a basis set. 

144 For example, all steps in the example will use aug-cc-pvdz 

145 because the last step is the only one which defines a basis. 

146 

147 Steps which change basis set must use the same theory. 

148 The following specification would perform SCF using the 3-21G 

149 basis set first, then B3LYP//3-21g, and then B3LYP//6-31G(2df,p) 

150 

151 >>> calc = NWChem(theory='dft', xc='b3lyp', basis='6-31g(2df,p)', 

152 >>> pretasks=[ 

153 >>> {'theory': 'scf', 'basis': '3-21g', 

154 >>> 'set': {'lindep:n_dep': 0}}, 

155 >>> {'dft': {'xc': 'b3lyp'}}, 

156 >>> ]) 

157 

158 The :code:`'set': {'lindep:n_dep': 0}` option is highly suggested 

159 as a way to avoid errors relating to symmetry changes between tasks. 

160 

161 The calculator will configure appropriate options for saving 

162 and loading intermediate wavefunctions, and 

163 place an "ignore" task directive between each step so that 

164 convergence errors in intermediate steps do not halt execution. 

165 """ 

166 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file, 

167 label, atoms, command, **kwargs) 

168 self.calc = None 

169 

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

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

172 

173 # Prepare perm and scratch directories 

174 perm = os.path.abspath(self.parameters.get('perm', self.label)) 

175 scratch = os.path.abspath(self.parameters.get('scratch', self.label)) 

176 os.makedirs(perm, exist_ok=True) 

177 os.makedirs(scratch, exist_ok=True) 

178 

179 io.write(self.label + '.nwi', atoms, properties=properties, 

180 label=self.label, **self.parameters) 

181 

182 def read_results(self): 

183 output = io.read(self.label + '.nwo') 

184 self.calc = output.calc 

185 self.results = output.calc.results 

186 

187 def band_structure(self): 

188 self.calculate() 

189 perm = self.parameters.get('perm', self.label) 

190 if self.calc.get_spin_polarized(): 

191 alpha = np.loadtxt(os.path.join(perm, self.label + '.alpha_band')) 

192 beta = np.loadtxt(os.path.join(perm, self.label + '.beta_band')) 

193 energies = np.array([alpha[:, 1:], beta[:, 1:]]) * Hartree 

194 else: 

195 data = np.loadtxt(os.path.join(perm, 

196 self.label + '.restricted_band')) 

197 energies = data[np.newaxis, :, 1:] * Hartree 

198 eref = self.calc.get_fermi_level() 

199 if eref is None: 

200 eref = 0. 

201 return BandStructure(self.parameters.bandpath, energies, eref)