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

2 

3https://nwchemgit.github.io 

4""" 

5import os 

6import numpy as np 

7 

8from ase import io 

9from ase.units import Hartree 

10from ase.calculators.calculator import FileIOCalculator 

11from ase.spectrum.band_structure import BandStructure 

12 

13 

14class NWChem(FileIOCalculator): 

15 implemented_properties = ['energy', 'free_energy', 

16 'forces', 'stress', 'dipole'] 

17 command = 'nwchem PREFIX.nwi > PREFIX.nwo' 

18 accepts_bandpath_keyword = True 

19 discard_results_on_any_change = True 

20 

21 def __init__(self, restart=None, 

22 ignore_bad_restart_file=FileIOCalculator._deprecated, 

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

24 """ 

25 NWChem keywords are specified using (potentially nested) 

26 dictionaries. Consider the following input file block: 

27 

28 >>> dft 

29 >>> odft 

30 >>> mult 2 

31 >>> convergence energy 1e-9 density 1e-7 gradient 5e-6 

32 >>> end 

33 

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

35 following settings: 

36 

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

38 >>> 'mult': 2, 

39 >>> 'convergence': {'energy': 1e-9, 

40 >>> 'density': 1e-7, 

41 >>> 'gradient': 5e-6, 

42 >>> }, 

43 >>> }, 

44 >>> ) 

45 

46 In addition, the calculator supports several special keywords: 

47 

48 theory: str 

49 Which NWChem module should be used to calculate the 

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

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

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

53 calculator will attempt to guess which theory to use based 

54 on the keywords provided by the user. 

55 xc: str 

56 The exchange-correlation functional to use. Only relevant 

57 for DFT calculations. 

58 task: str 

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

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

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

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

63 should not be set manually. 

64 basis: str or dict 

65 Which basis set to use for gaussian-type orbital 

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

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

68 provide a dict of the form: 

69 

70 >>> calc = NWChem(..., 

71 >>> basis={'O': '3-21G', 

72 >>> 'Si': '6-31g'}) 

73 

74 basispar: str 

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

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

77 symmetry: int or str 

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

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

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

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

82 autosym: bool 

83 Whether NWChem should automatically determine the symmetry 

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

85 center: bool 

86 Whether NWChem should automatically center the structure 

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

88 autoz: bool 

89 Whether NWChem should automatically construct a Z-matrix 

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

91 geompar: str 

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

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

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

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

96 keyword described above for these settings. 

97 set: dict 

98 Used to manually create or modify entries in the NWChem 

99 rtdb. For example, the following settings enable 

100 pseudopotential filtering for plane-wave calculations: 

101 

102 >>> set nwpw:kbpp_ray .true. 

103 >>> set nwpw:kbpp_filter .true. 

104 

105 These settings are generated by the NWChem calculator by 

106 passing the arguments: 

107 

108 >>> calc = NWChem(..., 

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

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

111 

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

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

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

115 bandpath: BandPath object 

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

117 Implies ``theory='band'``. 

118 """ 

119 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file, 

120 label, atoms, command, **kwargs) 

121 self.calc = None 

122 

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

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

125 

126 # Prepare perm and scratch directories 

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

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

129 os.makedirs(perm, exist_ok=True) 

130 os.makedirs(scratch, exist_ok=True) 

131 

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

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

134 

135 def read_results(self): 

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

137 self.calc = output.calc 

138 self.results = output.calc.results 

139 

140 def band_structure(self): 

141 self.calculate() 

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

143 if self.calc.get_spin_polarized(): 

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

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

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

147 else: 

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

149 self.label + '.restricted_band')) 

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

151 eref = self.calc.get_fermi_level() 

152 if eref is None: 

153 eref = 0. 

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