Coverage for gpaw/xc/libvdwxc.py: 57%

1import numpy as np 

2 

3from gpaw.mpi import have_mpi 

4from gpaw.utilities import compiled_with_libvdwxc 

5from gpaw.utilities.grid_redistribute import Domains, general_redistribute 

6from gpaw.utilities.timing import nulltimer 

7from gpaw.xc.functional import XCFunctional 

8from gpaw.xc.lda import stress_lda_term 

9from gpaw.xc.gga import GGA, gga_vars, add_gradient_correction, stress_gga_term 

10from gpaw.xc.libxc import LibXC 

11from gpaw.xc.mgga import MGGA 

12import gpaw 

13import gpaw.cgpaw as cgpaw 

14 

15 

16def libvdwxc_has_mpi(): 

17 return have_mpi and cgpaw.libvdwxc_has('mpi') 

18 

19 

20def libvdwxc_has_pfft(): 

21 return have_mpi and cgpaw.libvdwxc_has('pfft') 

22 

23 

24def libvdwxc_has_spin(): 

25 try: 

26 return cgpaw.libvdwxc_has("spin") 

27 except SystemError as err: 

28 raise SystemError(interface_error) from err 

29 

30 

31def libvdwxc_has_stress(): 

32 try: 

33 val = cgpaw.libvdwxc_has("stress") 

34 except SystemError as err: 

35 raise SystemError(interface_error) from err 

36 return val 

37 

38 

39def check_grid_descriptor(gd): 

40 assert gd.parsize_c[1] == 1 and gd.parsize_c[2] == 1 

41 nxpts_p = gd.n_cp[0][1:] - gd.n_cp[0][:-1] 

42 nxpts0 = nxpts_p[0] 

43 for nxpts in nxpts_p[1:-1]: 

44 assert nxpts == nxpts0 

45 assert nxpts_p[-1] <= nxpts0 

46 

47 

48def get_domains(N_c, parsize_c): 

49 # We want a distribution like this: 

50 # [B, B, ..., B, remainder, 0, 0, ..., 0]. 

51 # with blocksize B chosen as large as possible for better load balance. 

52 # This function returns the blocksize and the cumulative sum of indices 

53 # starting with 0. 

54 blocksize_c = -(-N_c // parsize_c) 

55 return (np.arange(1 + parsize_c) * blocksize_c).clip(0, N_c) 

56 

57 

58def get_auto_pfft_grid(size): 

59 nproc1 = size 

60 nproc2 = 1 

61 while nproc1 > nproc2 and nproc1 % 2 == 0: 

62 nproc1 //= 2 

63 nproc2 *= 2 

64 return nproc1, nproc2 

65 

66 

67interface_error = ("The libvdwxc interface has changed. You might need to " 

68 "recompile the GPAW C extensions.") 

69spinwarning = """\ 

70GPAW uses the total density to evaluate the van der Waals functional 

71for a spin-polarized system. This is not entirely rigorous, so the 

72calculation cannot be considered a true vdW-DF-family calculation.""" 

73_VDW_NUMERICAL_CODES = {'vdW-DF': 1, 

74 'vdW-DF2': 2, 

75 'vdW-DF-cx': 3} 

76 

77 

78class LibVDWXC: 

79 """Minimum-tomfoolery object-oriented interface to libvdwxc.""" 

80 def __init__(self, funcname, N_c, cell_cv, comm, mode='auto', 

81 pfft_grid=None, nspins=1): 

82 self.initialized = False 

83 if not compiled_with_libvdwxc(): 

84 raise ImportError('libvdwxc not compiled into GPAW') 

85 

86 self.vdw_functional_name = funcname 

87 code = _VDW_NUMERICAL_CODES[funcname] 

88 self.shape = tuple(N_c) 

89 ptr = np.empty(1, np.intp) 

90 cgpaw.libvdwxc_create(ptr, code, nspins, self.shape, 

91 tuple(np.ravel(cell_cv))) 

92 # assign ptr only now that it is initialized (so __del__ always works) 

93 self._ptr = ptr 

94 

95 # Choose mode automatically if necessary: 

96 if mode == 'auto': 

97 if pfft_grid is not None: 

98 mode = 'pfft' 

99 elif comm.size > 1: 

100 mode = 'mpi' 

101 else: 

102 mode = 'serial' 

103 assert mode in ['serial', 'mpi', 'pfft'] 

104 

105 if mode != 'serial' and not have_mpi: 

106 raise ImportError('MPI not available for libvdwxc-%s ' 

107 'because GPAW is serial' % mode) 

108 

109 if mode == 'pfft': 

110 if pfft_grid is None: 

111 pfft_grid = get_auto_pfft_grid(comm.size) 

112 nx, ny = pfft_grid 

113 assert nx * ny == comm.size 

114 # User might have passed a list, but we make sure to store a tuple: 

115 self.pfft_grid = (nx, ny) 

116 elif pfft_grid is not None: 

117 raise ValueError('pfft_grid specified with mode %s' % mode) 

118 

119 self.mode = mode 

120 self.comm = comm 

121 

122 def initialize_backend(self): 

123 assert not self.initialized 

124 

125 mode = self.mode 

126 comm = self.comm 

127 if mode == 'serial': 

128 assert comm.size == 1, ('You cannot run in serial with %d cores' 

129 % comm.size) 

130 cgpaw.libvdwxc_init_serial(self._ptr) 

131 elif gpaw.dry_run: 

132 pass 

133 # XXXXXX liable to cause really nasty errors. 

134 elif mode == 'mpi': 

135 if not libvdwxc_has_mpi(): 

136 raise ImportError('libvdwxc not compiled with MPI') 

137 cgpaw.libvdwxc_init_mpi(self._ptr, comm.get_c_object()) 

138 elif mode == 'pfft': 

139 nx, ny = self.pfft_grid 

140 cgpaw.libvdwxc_init_pfft(self._ptr, comm.get_c_object(), nx, ny) 

141 

142 self.initialized = True 

143 

144 def calculate(self, n_sg, sigma_xg, dedn_sg, dedsigma_xg, 

145 get_stress=False): 

146 """Calculate energy and add partial derivatives to arrays.""" 

147 if not self.initialized: 

148 self.initialize_backend() 

149 

150 for arr in [n_sg, sigma_xg, dedn_sg, dedsigma_xg]: 

151 assert arr.flags.contiguous 

152 assert arr.dtype == float 

153 # XXX We cannot actually ask libvdwxc about its expected shape 

154 # assert arr.shape == self.shape, [arr.shape, self.shape] 

155 if get_stress: 

156 stress_vv = np.zeros((3, 3)) 

157 else: 

158 stress_vv = None 

159 try: 

160 energy = cgpaw.libvdwxc_calculate(self._ptr, stress_vv, n_sg, 

161 sigma_xg, dedn_sg, dedsigma_xg) 

162 except TypeError as err: 

163 raise TypeError(interface_error) from err 

164 return energy, stress_vv 

165 

166 def get_description(self): 

167 if self.mode == 'serial': 

168 pardesc = 'fftw3 serial' 

169 elif self.mode == 'mpi': 

170 size = self.comm.size 

171 cores = '%d cores' % size if size != 1 else 'one core' 

172 pardesc = 'fftw3-mpi with %s' % cores 

173 else: 

174 assert self.mode == 'pfft' 

175 pardesc = 'pfft with %d x %d CPU grid' % self.pfft_grid 

176 return f'{self.vdw_functional_name} [libvdwxc/{pardesc}]' 

177 

178 def tostring(self): 

179 return cgpaw.libvdwxc_tostring(self._ptr) 

180 

181 def __del__(self): 

182 if hasattr(self, '_ptr'): 

183 cgpaw.libvdwxc_free(self._ptr) 

184 

185 

186class RedistWrapper: 

187 """Call libvdwxc redistributing automatically from and to GPAW grid.""" 

188 def __init__(self, libvdwxc, distribution, timer=nulltimer, 

189 vdwcoef=1.0): 

190 # It is hacky for the RedistWrapper to apply the vdwcoef, but this 

191 # is the only accessible place where we take copies of the arrays, 

192 # and therefore the only 'good' place to apply a factor without 

193 # applying it to the existing contents of the array. 

194 self.libvdwxc = libvdwxc 

195 self.distribution = distribution 

196 self.timer = timer 

197 self.vdwcoef = vdwcoef 

198 

199 def calculate(self, n_sg, sigma_xg, v_sg, dedsigma_xg, get_stress=False): 

200 zeros = self.distribution.block_zeros 

201 sshape = (len(n_sg),) 

202 xshape = (len(sigma_xg),) 

203 nblock_sg = zeros(sshape) 

204 sigmablock_xg = zeros(xshape) 

205 vblock_sg = zeros(sshape) 

206 dedsigmablock_xg = zeros(xshape) 

207 

208 self.timer.start('redistribute') 

209 self.distribution.gd2block(n_sg, nblock_sg) 

210 self.distribution.gd2block(sigma_xg, sigmablock_xg) 

211 self.timer.stop('redistribute') 

212 

213 self.timer.start('libvdwxc nonlocal') 

214 energy, stress = self.libvdwxc.calculate(nblock_sg, sigmablock_xg, 

215 vblock_sg, dedsigmablock_xg, 

216 get_stress=get_stress) 

217 self.timer.stop('libvdwxc nonlocal') 

218 energy *= self.vdwcoef 

219 if get_stress: 

220 stress[:] *= self.vdwcoef 

221 for arr in vblock_sg, dedsigmablock_xg: 

222 arr *= self.vdwcoef 

223 

224 self.timer.start('redistribute') 

225 self.distribution.block2gd_add(vblock_sg, v_sg) 

226 self.distribution.block2gd_add(dedsigmablock_xg, dedsigma_xg) 

227 self.timer.stop('redistribute') 

228 return energy, stress 

229 

230 

231class FFTDistribution: 

232 def __init__(self, gd, parsize_c): 

233 self.input_gd = gd 

234 assert np.prod(parsize_c) == gd.comm.size 

235 self.local_input_size_c = gd.n_c 

236 self.domains_in = Domains(gd.n_cp) 

237 N_c = gd.get_size_of_global_array(pad=True) 

238 self.domains_out = Domains([get_domains(N_c[i], parsize_c[i]) 

239 for i in range(3)]) 

240 

241 if parsize_c[1] == 1 and parsize_c[2] == 1: 

242 def aux_rank_to_parpos(rank=gd.comm.rank): 

243 return np.array([rank, 0, 0], int) 

244 else: 

245 # For 2D distributions we use a grid descriptor. We could 

246 # actually use a grid descriptor always, but that causes trouble 

247 # when there are more cores than grid points, which could well 

248 # be the case when the distribution is only 1D. 

249 # 

250 # The auxiliary gd actually is used *only* for the rank/parpos 

251 # correspondence. The actual domains it defines are unused!! 

252 aux_gd = gd.new_descriptor(comm=gd.comm, parsize_c=parsize_c) 

253 aux_rank_to_parpos = aux_gd.get_processor_position_from_rank 

254 

255 parpos_c = aux_rank_to_parpos() 

256 self.aux_rank_to_parpos = aux_rank_to_parpos 

257 self.local_output_size_c = tuple(self.domains_out.get_box(parpos_c)[1]) 

258 

259 def block_zeros(self, shape=()): 

260 return np.zeros(shape + self.local_output_size_c) 

261 

262 def gd2block(self, a_xg, b_xg): 

263 general_redistribute(self.input_gd.comm, 

264 self.domains_in, self.domains_out, 

265 self.input_gd.get_processor_position_from_rank, 

266 self.aux_rank_to_parpos, 

267 a_xg, b_xg, behavior='overwrite') 

268 

269 def block2gd_add(self, a_xg, b_xg): 

270 general_redistribute(self.input_gd.comm, 

271 self.domains_out, self.domains_in, 

272 self.aux_rank_to_parpos, 

273 self.input_gd.get_processor_position_from_rank, 

274 a_xg, b_xg, behavior='add') 

275 

276 

277class VDWXC(XCFunctional): 

278 def __init__(self, semilocal_xc, name, mode='auto', 

279 pfft_grid=None, libvdwxc_name=None, 

280 setup_name='revPBE', vdwcoef=1.0, 

281 accept_partial_decomposition=False): 

282 """Initialize VDWXC object (further initialization required). 

283 

284 mode can be 'auto', 'serial', 'mpi', or 'pfft'. 

285 

286 * 'serial' uses FFTW and only works with serial decompositions. 

287 

288 * 'mpi' uses FFTW-MPI with communicator of the grid 

289 descriptor, parallelizing along the x axis. 

290 

291 * 'pfft' uses PFFT and works with any decomposition, 

292 parallelizing along two directions for best scalability. 

293 

294 * 'auto' uses PFFT if pfft_grid is given, else FFTW-MPI if the 

295 calculation uses more than one core, else serial FFTW. 

296 

297 pfft_grid is the 2D CPU grid used by PFFT and can be a tuple 

298 (nproc1, nproc2) that multiplies to total communicator size, 

299 or None. It is an error to specify pfft_grid unless using 

300 PFFT. If left unspecified, a hopefully reasonable automatic 

301 choice will be made. 

302 """ 

303 

304 # XXX We should probably not initialize with the same data as the 

305 # semilocal XC kernel. 

306 XCFunctional.__init__(self, semilocal_xc.kernel.name, 

307 semilocal_xc.kernel.type) 

308 # Really, 'type' should be something along the lines of vdw-df. 

309 self.semilocal_xc = semilocal_xc 

310 

311 # We set these in the initialize later (ugly). 

312 self.libvdwxc = None 

313 self.distribution = None 

314 self.redist_wrapper = None 

315 self.timer = nulltimer 

316 self.setup_name = setup_name 

317 

318 # These parameters are simply stored for later forwarding 

319 self.name = name 

320 if libvdwxc_name is None: 

321 libvdwxc_name = name 

322 self._libvdwxc_name = libvdwxc_name 

323 self._mode = mode 

324 self._pfft_grid = pfft_grid 

325 self._vdwcoef = vdwcoef 

326 self.accept_partial_decomposition = accept_partial_decomposition 

327 self._nspins = 1 

328 

329 self.last_nonlocal_energy = None 

330 self.last_semilocal_energy = None 

331 

332 # XXXXXXXXXXXXXXXXX 

333 self.calculate_paw_correction = semilocal_xc.calculate_paw_correction 

334 self.calculate_spherical = semilocal_xc.calculate_spherical 

335 self.apply_orbital_dependent_hamiltonian = \ 

336 semilocal_xc.apply_orbital_dependent_hamiltonian 

337 self.add_forces = semilocal_xc.add_forces 

338 self.get_kinetic_energy_correction = \ 

339 semilocal_xc.get_kinetic_energy_correction 

340 self.rotate = semilocal_xc.rotate 

341 

342 def __str__(self): 

343 tokens = [self._mode] 

344 if self._libvdwxc_name != self.name: 

345 tokens.append(f'nonlocal-name={self._libvdwxc_name}') 

346 tokens.append('gga-kernel={}' 

347 .format(self.semilocal_xc.kernel.name)) 

348 if self._pfft_grid is not None: 

349 tokens.append(f'pfft={self._pfft_grid}') 

350 if self._vdwcoef != 1.0: 

351 tokens.append(f'vdwcoef={self._vdwcoef}') 

352 

353 qualifier = ', '.join(tokens) 

354 return f'{self.name} [libvdwxc/{qualifier}]' 

355 

356 def todict(self): 

357 dct = dict(backend='libvdwxc', 

358 semilocal_xc=self.semilocal_xc.name, 

359 name=self.name, 

360 # mode=self._mode, 

361 # pfft_grid=self._pfft_grid, 

362 libvdwxc_name=self._libvdwxc_name, 

363 setup_name=self.setup_name, 

364 vdwcoef=self._vdwcoef) 

365 return dct 

366 

367 def set_grid_descriptor(self, gd): 

368 XCFunctional.set_grid_descriptor(self, gd) 

369 self.semilocal_xc.set_grid_descriptor(gd) 

370 

371 def get_description(self): 

372 lines = [] 

373 app = lines.append 

374 app(f'{self.name} with libvdwxc') 

375 mode = self.libvdwxc.mode 

376 ncores = self.libvdwxc.comm.size 

377 cores = 'core' if ncores == 1 else 'cores' 

378 if mode == 'mpi': 

379 mode = f'mpi with {ncores} {cores}' 

380 elif mode == 'pfft': 

381 nx, ny = self.libvdwxc.pfft_grid 

382 mode = f'pfft with {nx} x {ny} {cores}' 

383 app(f'Mode: {mode}') 

384 app(f'Semilocal: {self.semilocal_xc.get_description()}') 

385 if self.libvdwxc.vdw_functional_name != self.name: 

386 app('Corresponding non-local functional: {}' 

387 .format(self.libvdwxc.vdw_functional_name)) 

388 app('Local blocksize: {} x {} x {}' 

389 .format(*self.distribution.local_output_size_c)) 

390 app(f'PAW datasets: {self.get_setup_name()}') 

391 return '\n'.join(lines) 

392 

393 def summary(self, log): 

394 from ase.units import Hartree as Ha 

395 enl = self.libvdwxc.comm.sum_scalar(self.last_nonlocal_energy) 

396 esl = self.gd.comm.sum_scalar(self.last_semilocal_energy) 

397 # In the current implementation these communicators have the same 

398 # processes always: 

399 assert self.libvdwxc.comm.size == self.gd.comm.size 

400 log(f'Non-local {self.name} correlation energy: {(enl * Ha):.6f}') 

401 log(f'Semilocal {self.semilocal_xc.kernel.name} ' 

402 f'energy: {esl * Ha:.6f}') 

403 log('(Not including atomic contributions)') 

404 

405 def get_setup_name(self): 

406 return self.setup_name 

407 

408 def initialize_backend(self, gd): 

409 N_c = gd.get_size_of_global_array(pad=True) 

410 # garbage-collect before allocating next one, if an object was 

411 # already allocated: 

412 self.libvdwxc = None 

413 self.libvdwxc = LibVDWXC(self._libvdwxc_name, N_c, gd.cell_cv, 

414 gd.comm, mode=self._mode, 

415 pfft_grid=self._pfft_grid, 

416 nspins=self._nspins) 

417 cpugrid = [1, 1, 1] 

418 if self.libvdwxc.mode == 'mpi': 

419 cpugrid[0] = gd.comm.size 

420 elif self.libvdwxc.mode == 'pfft': 

421 cpugrid[0], cpugrid[1] = self.libvdwxc.pfft_grid 

422 self.distribution = FFTDistribution(gd, cpugrid) 

423 self.redist_wrapper = RedistWrapper(self.libvdwxc, 

424 self.distribution, 

425 self.timer, 

426 self._vdwcoef) 

427 

428 def set_positions(self, spos_ac): 

429 self.semilocal_xc.set_positions(spos_ac) 

430 

431 def initialize(self, density, hamiltonian, wfs): 

432 self.timer = hamiltonian.timer # fragile object robbery 

433 self.semilocal_xc.initialize(density, hamiltonian, wfs) 

434 gd = self.gd 

435 self._nspins = density.nspins 

436 

437 self.initialize_backend(self.gd) 

438 

439 if (wfs.world.size > gd.comm.size and np.prod(gd.N_c) > 64**3 

440 and not self.accept_partial_decomposition): 

441 # We could issue a warning if an excuse turns out to exist some day 

442 raise ValueError('You are using libvdwxc with only ' 

443 '%d out of %d available cores in a non-small ' 

444 'calculation (%s points). This is not ' 

445 'a crime but is likely silly and therefore ' 

446 'triggers an error. Please use ' 

447 'parallel={\'augment_grids\': True} ' 

448 'or complain to the developers. ' 

449 'You can also disable this error ' 

450 'by passing accept_partial_decomposition=True ' 

451 'to the libvdwxc functional object, ' 

452 'but be careful to use good domain ' 

453 'decomposition.' % 

454 (gd.comm.size, wfs.world.size, 

455 ' x '.join(str(N) for N in gd.N_c))) 

456 # TODO Here we could decide FFT padding. 

457 

458 def calculate_vdw(self, e_g, n_sg, v_sg, sigma_xg, dedsigma_xg): 

459 assert self.libvdwxc is not None 

460 # Note: Redistwrapper handles vdwcoef. For now 

461 energy_nonlocal = self.redist_wrapper.calculate(n_sg, sigma_xg, 

462 v_sg, dedsigma_xg)[0] 

463 self.last_semilocal_energy = e_g.sum() * self.gd.dv 

464 self.last_nonlocal_energy = energy_nonlocal 

465 e_g[0, 0, 0] += energy_nonlocal / self.gd.dv 

466 

467 def calculate_impl(self, gd, n_sg, v_sg, e_g): 

468 """Calculate energy and potential. 

469 

470 gd may be non-periodic. To be distinguished from self.gd 

471 which is always periodic due to priminess of FFT dimensions. 

472 (To do: proper padded FFTs.)""" 

473 assert gd == self.distribution.input_gd 

474 assert self.libvdwxc is not None 

475 semiloc = self.semilocal_xc 

476 

477 self.timer.start('van der Waals') 

478 

479 self.timer.start('semilocal') 

480 # XXXXXXX taken from GGA 

481 grad_v = self.grad_v 

482 sigma_xg, dedsigma_xg, gradn_svg = gga_vars(gd, grad_v, n_sg) 

483 n_sg[:] = np.abs(n_sg) # XXXX What to do about this? 

484 sigma_xg[:] = np.abs(sigma_xg) 

485 

486 # Grrr, interface still sucks 

487 if hasattr(semiloc, 'process_mgga'): 

488 semiloc.process_mgga(e_g, n_sg, v_sg, sigma_xg, dedsigma_xg) 

489 else: 

490 semiloc.kernel.calculate(e_g, n_sg, v_sg, sigma_xg, dedsigma_xg) 

491 self.timer.stop('semilocal') 

492 

493 self.calculate_vdw(e_g, n_sg, v_sg, sigma_xg, dedsigma_xg) 

494 add_gradient_correction(grad_v, gradn_svg, sigma_xg, dedsigma_xg, v_sg) 

495 self.timer.stop('van der Waals') 

496 

497 def stress_tensor_contribution(self, n_sg, skip_sum=False): 

498 assert self.libvdwxc is not None 

499 

500 stress0_vv = self.semilocal_xc.stress_tensor_contribution( 

501 n_sg, skip_sum=skip_sum 

502 ) 

503 

504 sigma_xg, dedsigma_xg, gradn_svg = gga_vars(self.gd, self.grad_v, n_sg) 

505 n_sg[:] = np.abs(n_sg) # XXXX What to do about this? 

506 sigma_xg[:] = np.abs(sigma_xg) 

507 v_sg = np.zeros_like(n_sg) 

508 dedsigma_xg = np.zeros_like(sigma_xg) 

509 

510 stress_vv = self.redist_wrapper.calculate(n_sg, sigma_xg, 

511 v_sg, dedsigma_xg, 

512 get_stress=True)[1] 

513 

514 stress_vv += stress_lda_term(self.gd, None, n_sg, v_sg) 

515 stress_vv += stress_gga_term(self.gd, sigma_xg, gradn_svg, dedsigma_xg) 

516 

517 if not skip_sum: 

518 self.gd.comm.sum(stress_vv) 

519 return stress_vv + stress0_vv 

520 

521 @property 

522 def grad_v(self): 

523 return self.semilocal_xc.grad_v 

524 

525 @property 

526 def dedtaut_sG(self): 

527 return self.semilocal_xc.dedtaut_sG 

528 

529 @property 

530 def tauct(self): 

531 if hasattr(self.semilocal_xc, "tauct"): 

532 return self.semilocal_xc.tauct 

533 else: 

534 raise ValueError("This vdW functional does not have tauct") 

535 

536 

537def vdw_df(**kwargs): 

538 kwargs1 = dict(name='vdW-DF', setup_name='revPBE', 

539 semilocal_xc=GGA(LibXC('GGA_X_PBE_R+LDA_C_PW'), 

540 stencil=kwargs.pop('stencil', 2))) 

541 kwargs1.update(kwargs) 

542 return VDWXC(**kwargs1) 

543 

544 

545def vdw_df2(**kwargs): 

546 kwargs1 = dict(name='vdW-DF2', setup_name='PBE', 

547 semilocal_xc=GGA(LibXC('GGA_X_RPW86+LDA_C_PW'), 

548 stencil=kwargs.pop('stencil', 2))) 

549 kwargs1.update(kwargs) 

550 return VDWXC(**kwargs1) 

551 

552 

553def vdw_df_cx(**kwargs): 

554 # cx semilocal exchange is in libxc 2.2.2 or newer (or maybe from older) 

555 kernel = kwargs.get('kernel') 

556 if kernel is None: 

557 kernel = LibXC('GGA_X_LV_RPW86+LDA_C_PW') 

558 

559 kwargs1 = dict(name='vdW-DF-cx', setup_name='PBE', 

560 # PBEsol is most correct but not distributed by default. 

561 semilocal_xc=GGA(kernel, stencil=kwargs.pop('stencil', 2))) 

562 kwargs1.update(kwargs) 

563 return VDWXC(**kwargs1) 

564 

565 

566def vdw_optPBE(**kwargs): 

567 kwargs1 = dict(name='vdW-optPBE', libvdwxc_name='vdW-DF', setup_name='PBE', 

568 semilocal_xc=GGA(LibXC('GGA_X_OPTPBE_VDW+LDA_C_PW'), 

569 stencil=kwargs.pop('stencil', 2))) 

570 kwargs1.update(kwargs) 

571 return VDWXC(**kwargs1) 

572 

573 

574def vdw_optB88(**kwargs): 

575 kwargs1 = dict(name='optB88', libvdwxc_name='vdW-DF', setup_name='PBE', 

576 semilocal_xc=GGA(LibXC('GGA_X_OPTB88_VDW+LDA_C_PW'), 

577 stencil=kwargs.pop('stencil', 2))) 

578 kwargs1.update(kwargs) 

579 return VDWXC(**kwargs1) 

580 

581 

582def vdw_C09(**kwargs): 

583 kwargs1 = dict(name='vdW-C09', libvdwxc_name='vdW-DF', setup_name='PBE', 

584 semilocal_xc=GGA(LibXC('GGA_X_C09X+LDA_C_PW'), 

585 stencil=kwargs.pop('stencil', 2))) 

586 kwargs1.update(kwargs) 

587 return VDWXC(**kwargs1) 

588 

589 

590def vdw_beef(**kwargs): 

591 # Kernel parameters stolen from vdw.py 

592 from gpaw.xc.bee import BEEVDWKernel 

593 kernel = BEEVDWKernel('BEE2', None, 

594 0.600166476948828631066, 

595 0.399833523051171368934) 

596 kwargs1 = dict(name='vdW-BEEF', libvdwxc_name='vdW-DF2', setup_name='PBE', 

597 semilocal_xc=GGA(kernel, stencil=kwargs.pop('stencil', 2))) 

598 kwargs1.update(kwargs) 

599 return VDWXC(**kwargs1) 

600 

601 

602def vdw_mbeef(**kwargs): 

603 # Note: Parameters taken from vdw.py 

604 from gpaw.xc.bee import BEEVDWKernel 

605 kernel = BEEVDWKernel('BEE3', None, 0.405258352, 0.356642240) 

606 kwargs1 = dict(name='vdW-mBEEF', setup_name='PBEsol', 

607 libvdwxc_name='vdW-DF2', vdwcoef=0.886774972, 

608 semilocal_xc=MGGA(kernel, stencil=kwargs.pop('stencil', 2))) 

609 kwargs1.update(kwargs) 

610 return VDWXC(**kwargs1) 

611 

612 

613# String to functional mapping 

614def get_libvdwxc_functional(name, **kwargs): 

615 if 'name' in kwargs: 

616 name2 = kwargs.pop('name') 

617 assert name == name2 

618 funcs = {'vdW-DF': vdw_df, 

619 'vdW-DF2': vdw_df2, 

620 'vdW-DF-cx': vdw_df_cx, 

621 'optPBE-vdW': vdw_optPBE, 

622 'optB88-vdW': vdw_optB88, 

623 'C09-vdW': vdw_C09, 

624 'BEEF-vdW': vdw_beef, 

625 'mBEEF-vdW': vdw_mbeef} 

626 

627 semilocal_xc = kwargs.pop('semilocal_xc', None) 

628 if semilocal_xc is not None: 

629 from gpaw.xc import XC 

630 semilocal_xc = XC(semilocal_xc) 

631 

632 func = funcs[name](**kwargs) 

633 if semilocal_xc is not None: 

634 assert semilocal_xc.name == func.semilocal_xc.name 

635 return func 

636 

637 

638class CXGGAKernel: 

639 def __init__(self, just_kidding=False): 

640 self.just_kidding = just_kidding 

641 self.type = 'GGA' 

642 self.lda_c = LibXC('LDA_C_PW') 

643 if self.just_kidding: 

644 self.name = 'purepython rPW86_with_%s' % self.lda_c.name 

645 else: 

646 self.name = 'purepython cx' 

647 

648 def calculate(self, e_g, n_sg, v_sg, sigma_xg, dedsigma_xg): 

649 e_g[:] = 0.0 

650 dedsigma_xg[:] = 0.0 

651 

652 self.lda_c.calculate(e_g, n_sg, v_sg, sigma_xg, dedsigma_xg) 

653 

654 for arr in [n_sg, v_sg, sigma_xg, dedsigma_xg]: 

655 assert len(arr) == 1 

656 self._exchange(n_sg[0], sigma_xg[0], e_g, v_sg[0], dedsigma_xg[0]) 

657 

658 def _exchange(self, rho, grho, sx, v1x, v2x): 

659 """Calculate cx local exchange. 

660 

661 Note that this *adds* to the energy density sx so that it can 

662 be called after LDA correlation part without ruining anything. 

663 Also it adds to v1x and v2x as is normal in GPAW.""" 

664 tol = 1e-20 

665 rho[rho < tol] = tol 

666 grho[grho < tol] = tol 

667 alp = 0.021789 

668 beta = 1.15 

669 a = 1.851 

670 b = 17.33 

671 c = 0.163 

672 mu_LM = 0.09434 

673 s_prefactor = 6.18733545256027 

674 Ax = -0.738558766382022 # = -3./4. * (3./pi)**(1./3) 

675 four_thirds = 4. / 3. 

676 

677 grad_rho = np.sqrt(grho) 

678 

679 # eventually we need s to power 12. Clip to avoid overflow 

680 # (We have individual tolerances on both rho and grho, but 

681 # they are not sufficient to guarantee this) 

682 s_1 = (grad_rho / (s_prefactor * rho**four_thirds)).clip(0.0, 1e20) 

683 s_2 = s_1 * s_1 

684 s_3 = s_2 * s_1 

685 s_4 = s_3 * s_1 

686 s_5 = s_4 * s_1 

687 s_6 = s_5 * s_1 

688 

689 fs_rPW86 = (1.0 + a * s_2 + b * s_4 + c * s_6)**(1. / 15.) 

690 

691 if self.just_kidding: 

692 fs = fs_rPW86 

693 else: 

694 fs = (1.0 + mu_LM * s_2) / (1.0 + alp * s_6) \ 

695 + alp * s_6 / (beta + alp * s_6) * fs_rPW86 

696 

697 # the energy density for the exchange. 

698 sx[:] += Ax * rho**four_thirds * fs 

699 

700 df_rPW86_ds = (1. / (15. * fs_rPW86**14.0)) * \ 

701 (2 * a * s_1 + 4 * b * s_3 + 6 * c * s_5) 

702 

703 if self.just_kidding: 

704 df_ds = df_rPW86_ds # XXXXXXXXXXXXXXXXXXXX 

705 else: 

706 df_ds = 1. / (1. + alp * s_6)**2 \ 

707 * (2.0 * mu_LM * s_1 * (1. + alp * s_6) - 

708 6.0 * alp * s_5 * (1. + mu_LM * s_2)) \ 

709 + alp * s_6 / (beta + alp * s_6) * df_rPW86_ds \ 

710 + 6.0 * alp * s_5 * fs_rPW86 / (beta + alp * s_6) \ 

711 * (1. - alp * s_6 / (beta + alp * s_6)) 

712 

713 # de/dn. This is the partial derivative of sx wrt. n, for s constant 

714 v1x[:] += Ax * four_thirds * (rho**(1. / 3.) * fs - 

715 grad_rho / (s_prefactor * rho) * df_ds) 

716 # de/d(nabla n). The other partial derivative 

717 v2x[:] += 0.5 * Ax * df_ds / (s_prefactor * grad_rho) 

718 # (We may or may not understand what that grad_rho is doing here.) 

719 

720 

721def test_derivatives(): 

722 gen = np.random.RandomState(1) 

723 shape = (1, 20, 20, 20) 

724 ngpts = np.prod(shape) 

725 n_sg = gen.rand(*shape) 

726 sigma_xg = np.zeros(shape) 

727 sigma_xg[:] = gen.rand(*shape) 

728 

729 qe_kernel = CXGGAKernel(just_kidding=True) 

730 libxc_kernel = LibXC('GGA_X_RPW86+LDA_C_PW') 

731 

732 cx_kernel = CXGGAKernel(just_kidding=False) 

733 

734 def check(kernel, n_sg, sigma_xg): 

735 e_g = np.zeros(shape[1:]) 

736 dedn_sg = np.zeros(shape) 

737 dedsigma_xg = np.zeros(shape) 

738 kernel.calculate(e_g, n_sg, dedn_sg, sigma_xg, dedsigma_xg) 

739 return e_g, dedn_sg, dedsigma_xg 

740 

741 def check_and_write(kernel): 

742 n1_sg = n_sg.copy() 

743 e_g, dedn_sg, dedsigma_xg = check(kernel, n_sg, sigma_xg) 

744 dedn = dedn_sg[0, 0, 0, 0] 

745 dedsigma = dedsigma_xg[0, 0, 0, 0] 

746 

747 dn = 1e-6 

748 n1_sg = n_sg.copy() 

749 n1_sg[0, 0, 0, 0] -= dn / 2. 

750 e1_g, _, _ = check(kernel, n1_sg, sigma_xg) 

751 

752 n1_sg[0, 0, 0, 0] += dn 

753 e2_g, _, _ = check(kernel, n1_sg, sigma_xg) 

754 

755 dedn_fd = (e2_g[0, 0, 0] - e1_g[0, 0, 0]) / dn