Coverage for gpaw/test/gpwfile.py: 31%

1import functools

2from math import sqrt

3from pathlib import Path

4import numpy as np

6from ase import Atom, Atoms

7from ase.units import Bohr

8from ase.build import bulk, molecule

9from ase.lattice.compounds import L1_2

10from ase.lattice.hexagonal import Graphene

11from gpaw import Davidson, FermiDirac, GPAW, Mixer, PW, FD, LCAO

12from gpaw.directmin.etdm_fdpw import FDPWETDM

13from gpaw.directmin.etdm_lcao import LCAOETDM

14from gpaw.directmin.tools import excite

15from gpaw.directmin.derivatives import Davidson as SICDavidson

16from gpaw.mom import prepare_mom_calculation

17from gpaw.mpi import world, serial_comm

18from gpaw.new.ase_interface import GPAW as GPAWNew

19from gpaw.poisson import FDPoissonSolver, PoissonSolver

20from gpaw.test.cachedfilehandler import CachedFilesHandler

22response_band_cutoff = {}

23_all_gpw_methodnames = set()

26def with_band_cutoff(*, gpw, band_cutoff):

27 # Store the band cutoffs in a dictionary to aid response tests

28 response_band_cutoff[gpw] = band_cutoff

30 def decorator(func):

31 @functools.wraps(func)

32 def wrapper(*args, **kwargs):

33 return func(*args, band_cutoff=band_cutoff, **kwargs)

34 return wrapper

36 return decorator

39def partial(fun, **kwargs):

40 import copy

41 kwargs = copy.deepcopy(kwargs)

43 def f(self):

44 return fun(self, **kwargs)

45 return f

48def _si_gw(self, *, a, symm, name):

49 atoms = self.generate_si_systems()[a]

50 return self._si_gw(atoms=atoms,

51 symm=symm,

52 name=f'{name}.txt')

55def si_gpwfiles():

56 gpw_file_dict = {}

57 for a in [0, 1]:

58 for symm, name1 in [({}, 'all'), ('off', 'no'),

59 ({'point_group': False}, 'tr'),

60 ({'time_reversal': False}, 'pg')]:

61 name = f'si_gw_a{a}_{name1}'

62 """

63 In !2153, a bug related to late binding of local functions was

64 fixed, and the partial wrapper utilized here is a temporary fix.

65 Previously, the test would only test the last symmetry in the loop,

66 but four times.

67 """

68 fun = partial(_si_gw, a=a, symm=symm, name=name)

69 fun.__name__ = name

70 gpw_file_dict[name] = gpwfile(fun)

72 return gpw_file_dict

75def gpwfile(meth):

76 """Decorator to identify the methods that produce gpw files."""

77 _all_gpw_methodnames.add(meth.__name__)

78 return meth

81class GPWFiles(CachedFilesHandler):

82 """Create gpw-files."""

84 def __init__(self, folder: Path):

85 super().__init__(folder, '.gpw')

87 def _calculate_and_write(self, name, work_path):

88 calc = getattr(self, name)()

89 calc.write(work_path, mode='all')

91 @gpwfile

92 def bcc_li_pw(self):

93 return self.bcc_li({'name': 'pw', 'ecut': 200})

95 @gpwfile

96 def bcc_li_fd(self):

97 return self.bcc_li({'name': 'fd'})

99 @gpwfile

100 def bcc_li_lcao(self):

101 return self.bcc_li({'name': 'lcao'})

102

103 def bcc_li(self, mode):

104 li = bulk('Li', 'bcc', 3.49)

105 li.calc = GPAW(mode=mode,

106 kpts=(3, 3, 3),

107 txt=self.folder / f'bcc_li_{mode["name"]}.txt')

108 li.get_potential_energy()

109 return li.calc

110

111 @gpwfile

112 def be_atom_fd(self):

113 atoms = Atoms('Be', [(0, 0, 0)], pbc=False)

114 atoms.center(vacuum=6)

115 calc = GPAW(mode='fd', h=0.35, symmetry={'point_group': False},

116 txt=self.folder / 'be_atom_fd.txt')

117 atoms.calc = calc

118 atoms.get_potential_energy()

119 return atoms.calc

120

121 @gpwfile

122 def fcc_Ni_col(self):

123 return self._fcc_Ni('col')

124

125 @gpwfile

126 def fcc_Ni_ncol(self):

127 return self._fcc_Ni('ncol')

128

129 @gpwfile

130 def fcc_Ni_ncolsoc(self):

131 return self._fcc_Ni('ncolsoc')

132

133 def _fcc_Ni(self, calc_type):

134 Ni = bulk('Ni', 'fcc', 3.48)

135 Ni.center()

136

137 mm = 0.5

138 easy_axis = 1 / np.sqrt(3) * np.ones(3)

139 Ni.set_initial_magnetic_moments([mm])

140

141 symmetry = {'point_group': False, 'time_reversal': True} if \

142 calc_type == 'col' else 'off'

143 magmoms = None if calc_type == 'col' else [mm * easy_axis]

144 soc = True if calc_type == 'ncolsoc' else False

145

146 Ni.calc = GPAWNew(mode={'name': 'pw', 'ecut': 380},

147 xc='LDA',

148 nbands='200%',

149 kpts={'size': (4, 4, 4), 'gamma': True},

150 parallel={'domain': 1, 'band': 1},

151 mixer={'beta': 0.5},

152 symmetry=symmetry,

153 occupations={'name': 'fermi-dirac', 'width': 0.05},

154 convergence={'density': 1e-8,

155 'bands': 'CBM+10',

156 'eigenstates': 1e-12},

157 magmoms=magmoms,

158 soc=soc,

159 txt=self.folder / f'fcc_Ni_{calc_type}.txt')

160 Ni.get_potential_energy()

161 return Ni.calc

162

163 @gpwfile

164 def h2_pw(self):

165 return self.h2({'name': 'pw', 'ecut': 200})

166

167 @gpwfile

168 def h2_fd(self):

169 return self.h2({'name': 'fd'})

170

171 @gpwfile

172 def h2_lcao(self):

173 return self.h2({'name': 'lcao'})

174

175 def h2(self, mode):

176 h2 = Atoms('H2', positions=[[0, 0, 0], [0.74, 0, 0]])

177 h2.center(vacuum=2.5)

178 h2.calc = GPAW(mode=mode,

179 txt=self.folder / f'h2_{mode["name"]}.txt')

180 h2.get_potential_energy()

181 return h2.calc

182

183 @gpwfile

184 def h2_pw_0(self):

185 h2 = Atoms('H2',

186 positions=[[-0.37, 0, 0], [0.37, 0, 0]],

187 cell=[5.74, 5, 5],

188 pbc=True)

189 h2.calc = GPAW(mode={'name': 'pw', 'ecut': 200},

190 txt=self.folder / 'h2_pw_0.txt')

191 h2.get_potential_energy()

192 return h2.calc

193

194 @gpwfile

195 def h2_bcc_afm(self):

196 a = 2.75

197 atoms = bulk(name='H', crystalstructure='bcc', a=a, cubic=True)

198 atoms.set_initial_magnetic_moments([1., -1.])

199

200 atoms.calc = GPAW(xc='LDA',

201 txt=self.folder / 'h2_bcc_afm.txt',

202 mode=PW(250),

203 nbands=4,

204 convergence={'bands': 4},

205 kpts={'density': 2.0, 'gamma': True})

206 atoms.get_potential_energy()

207 return atoms.calc

208

209 @gpwfile

210 def c2h4_do_fd(self):

211 atm = Atoms('CCHHHH',

212 positions=[

213 [-0.66874198, -0.00001714, -0.00001504],

214 [0.66874210, 0.00001699, 0.00001504],

215 [-1.24409879, 0.00000108, -0.93244784],

216 [-1.24406253, 0.00000112, 0.93242153],

217 [1.24406282, -0.93242148, 0.00000108],

218 [1.24409838, 0.93244792, 0.00000112]

219 ]

220 )

221 atm.center(vacuum=4.0)

222 atm.set_pbc(False)

223 atm.calc = GPAW(mode=FD(), h=0.3,

224 xc='PBE',

225 occupations={'name': 'fixed-uniform'},

226 eigensolver={'name': 'etdm-fdpw',

227 'converge_unocc': True},

228 mixer={'backend': 'no-mixing'},

229 spinpol=True,

230 symmetry='off',

231 nbands=-5,

232 convergence={'eigenstates': 4.0e-6},

233 )

234 atm.get_potential_energy()

235 return atm.calc

236

237 @gpwfile

238 def c2h4_do_pw(self):

239 atm = Atoms(

240 "CCHHHH",

241 positions=[

242 [-0.66874198, -0.00001714, -0.00001504],

243 [0.66874210, 0.00001699, 0.00001504],

244 [-1.24409879, 0.00000108, -0.93244784],

245 [-1.24406253, 0.00000112, 0.93242153],

246 [1.24406282, -0.93242148, 0.00000108],

247 [1.24409838, 0.93244792, 0.00000112],

248 ],

249 )

250 atm.center(vacuum=4.0)

251 atm.set_pbc(False)

252 atm.calc = GPAW(

253 mode=PW(300, force_complex_dtype=True),

254 xc="PBE",

255 occupations={"name": "fixed-uniform"},

256 eigensolver={"name": "etdm-fdpw", "converge_unocc": True},

257 mixer={"backend": "no-mixing"},

258 spinpol=True,

259 symmetry="off",

260 nbands=-5,

261 convergence={"eigenstates": 4.0e-6},

262 )

263 atm.get_potential_energy()

264 return atm.calc

265

266 def h3_maker(self, vacuum=2.0, pbc=True, **kwargs):

267 atoms = Atoms('H3', positions=[(0, 0, 0),

268 (0.59, 0, 0),

269 (1.1, 0, 0)])

270 atoms.center(vacuum=vacuum)

271 atoms.set_pbc(pbc)

272 return atoms

273

274 @gpwfile

275 def h3_do_num_pw_complex(self):

276 atoms = self.h3_maker()

277 calc = GPAW(mode=PW(300, force_complex_dtype=True),

278 basis='sz(dzp)',

279 h=0.3,

280 spinpol=False,

281 convergence={'energy': np.inf,

282 'eigenstates': np.inf,

283 'density': np.inf,

284 'minimum iterations': 1},

285 eigensolver=FDPWETDM(converge_unocc=True),

286 occupations={'name': 'fixed-uniform'},

287 mixer={'backend': 'no-mixing'},

288 nbands='nao',

289 symmetry='off')

290 atoms.calc = calc

291 atoms.get_potential_energy()

292 return atoms.calc

293

294 @gpwfile

295 def h3_do_num_pw(self):

296 atoms = self.h3_maker()

297 calc = GPAW(mode=PW(300, force_complex_dtype=False),

298 basis='sz(dzp)',

299 h=0.3,

300 spinpol=False,

301 convergence={'energy': np.inf,

302 'eigenstates': np.inf,

303 'density': np.inf,

304 'minimum iterations': 1},

305 eigensolver=FDPWETDM(converge_unocc=True),

306 occupations={'name': 'fixed-uniform'},

307 mixer={'backend': 'no-mixing'},

308 nbands='nao',

309 symmetry='off')

310 atoms.calc = calc

311 atoms.get_potential_energy()

312 return atoms.calc

313

314 @gpwfile

315 def h3_do_sd_lcao(self):

316 atoms = self.h3_maker(pbc=False)

317 atoms.set_initial_magnetic_moments([1, 0, 0])

318 calc = GPAW(mode='lcao',

319 h=0.3,

320 spinpol=True,

321 convergence={'energy': 0.1,

322 'eigenstates': 1e-4,

323 'density': 1e-4},

324 eigensolver={'name': 'etdm-lcao'},

325 occupations={'name': 'fixed-uniform'},

326 mixer={'backend': 'no-mixing'},

327 nbands='nao',

328 symmetry='off')

329 atoms.calc = calc

330 atoms.get_potential_energy()

331 return atoms.calc

332

333 @gpwfile

334 def h3_do_num_lcao(self):

335 atm = self.h3_maker()

336 atm.calc = GPAW(mode=LCAO(force_complex_dtype=True),

337 basis='sz(dzp)',

338 h=0.3,

339 spinpol=False,

340 convergence={'eigenstates': 10.0,

341 'density': 10.0,

342 'energy': 10.0},

343 occupations={'name': 'fixed-uniform'},

344 eigensolver={'name': 'etdm-lcao',

345 'matrix_exp': 'egdecomp'},

346 mixer={'backend': 'no-mixing'},

347 nbands='nao',

348 symmetry='off',

349 txt=None)

350 atm.get_potential_energy()

351 return atm.calc

352

353 def h2o_maker(self, vacuum, t=np.pi / 180 * 104.51, eps=0, **kwargs):

354 d = 0.9575

355 H2O = Atoms('OH2',

356 positions=[(0, 0, 0),

357 (d + eps, 0, 0),

358 (d * np.cos(t), d * np.sin(t), 0)],

359 **kwargs)

360 H2O.center(vacuum=vacuum)

361 return H2O

362

363 @gpwfile

364 def h2o_do_gmf_lcao(self):

365 atm = self.h2o_maker(vacuum=4.0)

366 atm.calc = GPAW(

367 mode=LCAO(),

368 basis="dzp",

369 h=0.22,

370 occupations={"name": "fixed-uniform"},

371 eigensolver="etdm-lcao",

372 mixer={"backend": "no-mixing"},

373 nbands="nao",

374 symmetry="off",

375 spinpol=True,

376 convergence={"density": 1.0e-4, "eigenstates": 4.0e-8},

377 )

378 atm.get_potential_energy()

379 return atm.calc

380

381 @gpwfile

382 def h2o_do_lcao(self):

383 atm = self.h2o_maker(vacuum=5.0)

384 atm.calc = GPAW(

385 mode=LCAO(),

386 basis="dzp",

387 occupations={"name": "fixed-uniform"},

388 eigensolver="etdm",

389 mixer={"backend": "no-mixing"},

390 nbands="nao",

391 symmetry="off",

392 )

393 atm.get_potential_energy()

394 return atm.calc

395

396 @gpwfile

397 def h2o_cdo_lcao(self):

398 atm = self.h2o_maker(vacuum=4.0)

399 atm.calc = GPAW(

400 mode=LCAO(),

401 basis="dzp",

402 h=0.22,

403 occupations={"name": "fixed-uniform"},

404 eigensolver={"name": "etdm-lcao"},

405 mixer={"backend": "no-mixing"},

406 nbands="nao",

407 symmetry="off",

408 spinpol=True,

409 convergence={"density": 1.0e-4, "eigenstates": 4.0e-8},

410 )

411 atm.get_potential_energy()

412 return atm.calc

413

414 @gpwfile

415 def h2o_cdo_lcao_sic(self):

416 atm = self.h2o_maker(vacuum=4.0)

417 atm.calc = GPAW(

418 mode=LCAO(force_complex_dtype=True),

419 h=0.22,

420 occupations={"name": "fixed-uniform"},

421 eigensolver={

422 "name": "etdm-lcao",

423 "localizationtype": "PM_PZ",

424 "localizationseed": 42,

425 "subspace_convergence": 1e-3,

426 "functional": {"name": "PZ-SIC", "scaling_factor": (0.5, 0.5)},

427 },

428 convergence={"eigenstates": 1e-4},

429 mixer={"backend": "no-mixing"},

430 nbands="nao",

431 symmetry="off",

432 )

433 atm.get_potential_energy()

434 return atm.calc

435

436 @gpwfile

437 def h2o_fdsic(self):

438 atm = self.h2o_maker(vacuum=4.0,

439 t=np.pi / 180 * (104.51 + 2.0),

440 eps=0.02)

441 atm.calc = GPAW(

442 mode=FD(force_complex_dtype=True),

443 h=0.25,

444 occupations={"name": "fixed-uniform"},

445 eigensolver=FDPWETDM(

446 functional={"name": "PZ-SIC", "scaling_factor": (0.5, 0.5)},

447 localizationseed=42,

448 localizationtype="FB_ER",

449 grad_tol_pz_localization=1.0e-3,

450 maxiter_pz_localization=200,

451 converge_unocc=True,

452 ),

453 convergence={"eigenstates": 1e-4},

454 mixer={"backend": "no-mixing"},

455 symmetry="off",

456 spinpol=True,

457 )

458 atm.get_potential_energy()

459 return atm.calc

460

461 @gpwfile

462 def h2o_lcaosic(self):

463 atm = self.h2o_maker(vacuum=4.0)

464 atm.calc = GPAW(

465 mode=LCAO(force_complex_dtype=True),

466 h=0.22,

467 occupations={"name": "fixed-uniform"},

468 eigensolver=LCAOETDM(

469 localizationtype="PM_PZ",

470 localizationseed=42,

471 functional={"name": "PZ-SIC", "scaling_factor": (0.5, 0.5)},

472 ),

473 convergence={"eigenstates": 1e-4},

474 mixer={"backend": "no-mixing"},

475 nbands="nao",

476 symmetry="off",

477 )

478 atm.get_potential_energy()

479 return atm.calc

480

481 @gpwfile

482 def h2o_mom_lcaosic(self):

483 atm = self.h2o_maker(vacuum=3.0)

484 calc = GPAW(

485 mode=LCAO(force_complex_dtype=True),

486 h=0.24,

487 basis="sz(dzp)",

488 spinpol=True,

489 symmetry="off",

490 eigensolver="etdm-lcao",

491 mixer={"backend": "no-mixing"},

492 occupations={"name": "fixed-uniform"},

493 convergence={"eigenstates": 1e-4},

494 )

495 atm.calc = calc

496 atm.get_potential_energy()

497 atm.calc.set(eigensolver=LCAOETDM(excited_state=True))

498 f_sn = excite(atm.calc, 0, 0, (0, 0))

499 prepare_mom_calculation(atm.calc, atm, f_sn)

500 atm.get_potential_energy()

501 atm.calc.set(

502 eigensolver=LCAOETDM(

503 searchdir_algo={"name": "l-sr1p"},

504 linesearch_algo={"name": "max-step"},

505 need_init_orbs=False,

506 localizationtype="PM_PZ",

507 localizationseed=42,

508 functional={"name": "pz-sic", "scaling_factor": (0.5, 0.5)},

509 ),

510 convergence={"eigenstates": 1e-2},

511 )

512 atm.get_potential_energy()

513 return atm.calc

514

515 @gpwfile

516 def h2o_gmf_lcaosic(self):

517 atm = self.h2o_maker(vacuum=3.0)

518 calc = GPAW(

519 mode=LCAO(),

520 basis="sz(dzp)",

521 h=0.24,

522 occupations={"name": "fixed-uniform"},

523 eigensolver="etdm-lcao",

524 convergence={"eigenstates": 1e-4},

525 mixer={"backend": "no-mixing"},

526 nbands="nao",

527 spinpol=True,

528 symmetry="off",

529 )

530 atm.calc = calc

531 atm.get_potential_energy()

532 atm.calc.set(eigensolver=LCAOETDM(excited_state=True))

533 f_sn = excite(atm.calc, 0, 0, (0, 0))

534 prepare_mom_calculation(atm.calc, atm, f_sn)

535 atm.get_potential_energy()

536 dave = SICDavidson(atm.calc.wfs.eigensolver, None)

537 appr_sp_order = dave.estimate_sp_order(atm.calc)

538

539 for kpt in atm.calc.wfs.kpt_u:

540 f_sn[kpt.s] = kpt.f_n

541 atm.calc.set(

542 eigensolver=LCAOETDM(

543 partial_diagonalizer={

544 "name": "Davidson",

545 "seed": 42,

546 "m": 20,

547 "eps": 5e-3,

548 "remember_sp_order": True,

549 "sp_order": appr_sp_order,

550 },

551 linesearch_algo={"name": "max-step"},

552 searchdir_algo={"name": "LBFGS-P_GMF"},

553 localizationtype="PM",

554 functional={"name": "PZ-SIC", "scaling_factor": (0.5, 0.5)},

555 need_init_orbs=False,

556 ),

557 occupations={"name": "mom",

558 "numbers": f_sn,

559 "use_fixed_occupations": True},

560 )

561 atm.get_potential_energy()

562 return atm.calc

563

564 @gpwfile

565 def h2o_mom_pwsic(self):

566 atm = self.h2o_maker(vacuum=3.0)

567 calc = GPAW(

568 mode=PW(300, force_complex_dtype=True),

569 spinpol=True,

570 symmetry="off",

571 eigensolver=FDPWETDM(converge_unocc=True),

572 mixer={"backend": "no-mixing"},

573 occupations={"name": "fixed-uniform"},

574 convergence={"eigenstates": 1e-4},

575 )

576 atm.calc = calc

577 atm.get_potential_energy()

578 atm.calc.set(eigensolver=FDPWETDM(excited_state=True))

579 f_sn = excite(atm.calc, 0, 0, (0, 0))

580 prepare_mom_calculation(atm.calc, atm, f_sn)

581 atm.get_potential_energy()

582 atm.calc.set(

583 eigensolver=FDPWETDM(

584 excited_state=True,

585 need_init_orbs=False,

586 functional={"name": "PZ-SIC",

587 "scaling_factor": (0.5, 0.5)}, # SIC/2

588 localizationseed=42,

589 localizationtype="PM",

590 grad_tol_pz_localization=1.0e-2,

591 printinnerloop=False,

592 grad_tol_inner_loop=1.0e-2,

593 ),

594 convergence={"eigenstates": 1e-3, "density": 1e-3},

595 )

596 atm.get_potential_energy()

597 return atm.calc

598

599 @gpwfile

600 def h2o_pwsic(self):

601 atm = self.h2o_maker(vacuum=4.0,

602 t=np.pi / 180 * (104.51 + 2.0),

603 eps=0.02)

604 atm.calc = GPAW(

605 mode=PW(300, force_complex_dtype=True),

606 occupations={"name": "fixed-uniform"},

607 eigensolver=FDPWETDM(

608 functional={"name": "pz-sic", "scaling_factor": (0.5, 0.5)},

609 localizationseed=42,

610 localizationtype="FB_ER",

611 grad_tol_pz_localization=5.0e-3,

612 maxiter_pz_localization=200,

613 converge_unocc=True,

614 ),

615 convergence={"eigenstates": 1e-4},

616 mixer={"backend": "no-mixing"},

617 symmetry="off",

618 spinpol=True,

619 )

620 atm.get_potential_energy()

621 return atm.calc

622

623 @gpwfile

624 def h2o_mom_do_pw(self):

625 atm = self.h2o_maker(vacuum=4.0)

626 calc = GPAW(

627 mode=PW(300),

628 spinpol=True,

629 symmetry="off",

630 eigensolver={"name": "etdm-fdpw", "converge_unocc": True},

631 mixer={"backend": "no-mixing"},

632 occupations={"name": "fixed-uniform"},

633 convergence={"eigenstates": 1e-4},

634 txt=None,

635 )

636 atm.calc = calc

637 atm.get_potential_energy()

638 return atm.calc

639

640 @gpwfile

641 def co_mom_do_lcao_forces(self):

642 L = 4.0

643 d = 1.13

644 atoms = Atoms(

645 "CO",

646 [[0.5 * L, 0.5 * L, 0.5 * L - 0.5 * d],

647 [0.5 * L, 0.5 * L, 0.5 * L + 0.5 * d]],

648 )

649 atoms.set_cell([L, L, L])

650 atoms.set_pbc(True)

651 calc = GPAW(

652 mode="lcao",

653 basis="dzp",

654 h=0.22,

655 xc="PBE",

656 spinpol=True,

657 symmetry="off",

658 occupations={"name": "fixed-uniform"},

659 eigensolver={"name": "etdm-lcao", "linesearch_algo": "max-step"},

660 mixer={"backend": "no-mixing"},

661 nbands="nao",

662 convergence={"density": 1.0e-4, "eigenstates": 4.0e-8},

663 )

664 atoms.calc = calc

665 atoms.get_potential_energy()

666 return atoms.calc

667

668 @gpwfile

669 def h2o_mom_do_lcao(self):

670 atm = self.h2o_maker(vacuum=4.0)

671 atm.calc = GPAW(

672 mode=LCAO(),

673 basis="dzp",

674 h=0.22,

675 occupations={"name": "fixed-uniform"},

676 eigensolver="etdm-lcao",

677 mixer={"backend": "no-mixing"},

678 nbands="nao",

679 symmetry="off",

680 spinpol=True,

681 convergence={"density": 1.0e-4, "eigenstates": 4.0e-8},

682 )

683 atm.get_potential_energy()

684 return atm.calc

685

686 @gpwfile

687 def h2o_pz_localization_pw(self):

688 atm = self.h2o_maker(vacuum=3.0,

689 t=np.pi / 180 * (104.51 + 2.0),

690 eps=0.02)

691 atm.calc = GPAW(

692 mode=PW(300, force_complex_dtype=True),

693 occupations={"name": "fixed-uniform"},

694 convergence={

695 "energy": np.inf,

696 "eigenstates": np.inf,

697 "density": np.inf,

698 "minimum iterations": 0,

699 },

700 eigensolver=FDPWETDM(converge_unocc=False),

701 mixer={"backend": "no-mixing"},

702 symmetry="off",

703 spinpol=True,

704 )

705 atm.get_potential_energy()

706 atm.calc.set(

707 eigensolver=FDPWETDM(

708 functional={"name": "PZ-SIC", "scaling_factor": (0.5, 0.5)},

709 localizationseed=42,

710 localizationtype="KS_PZ",

711 localization_tol=5.0e-2,

712 converge_unocc=False,

713 )

714 )

715 atm.get_potential_energy()

716 return atm.calc

717

718 @gpwfile

719 def c2h4_do_lcao(self):

720 atoms = Atoms(

721 "C2H4",

722 [

723 [6.68748500e-01, 2.00680000e-04, 5.55800000e-05],

724 [-6.68748570e-01, -2.00860000e-04, -5.51500000e-05],

725 [4.48890600e-01, -5.30146300e-01, 9.32670330e-01],

726 [4.48878120e-01, -5.30176640e-01, -9.32674730e-01],

727 [-1.24289513e00, 1.46164400e-02, 9.32559990e-01],

728 [-1.24286000e00, -1.46832100e-02, -9.32554970e-01],

729 ],

730 )

731 atoms.center(vacuum=4)

732 eigensolver = LCAOETDM(

733 searchdir_algo={"name": "l-sr1p"},

734 linesearch_algo={"name": "max-step"}

735 )

736 calc = GPAW(

737 mode="lcao",

738 basis="dzp",

739 h=0.24,

740 xc="PBE",

741 symmetry="off",

742 occupations={"name": "fixed-uniform"},

743 eigensolver=eigensolver,

744 mixer={"backend": "no-mixing"},

745 nbands="nao",

746 convergence={"density": 1.0e-4, "eigenstates": 4.0e-8},

747 )

748 atoms.calc = calc

749 atoms.get_potential_energy()

750 return atoms.calc

751

752 @gpwfile

753 def h3_orthonorm_lcao(self):

754 atoms = Atoms("H3", positions=[(0, 0, 0), (0.59, 0, 0), (1.1, 0, 0)])

755 atoms.set_initial_magnetic_moments([1, 0, 0])

756

757 atoms.center(vacuum=2.0)

758 atoms.set_pbc(False)

759 calc = GPAW(

760 mode="lcao",

761 basis="sz(dzp)",

762 h=0.3,

763 spinpol=True,

764 convergence={

765 "energy": np.inf,

766 "eigenstates": np.inf,

767 "density": np.inf,

768 "minimum iterations": 1,

769 },

770 eigensolver={"name": "etdm-lcao"},

771 occupations={"name": "fixed-uniform"},

772 mixer={"backend": "no-mixing"},

773 nbands="nao",

774 symmetry="off",

775 txt=None,

776 )

777 atoms.calc = calc

778 atoms.get_potential_energy()

779 return atoms.calc

780

781 @gpwfile

782 def h2_sic_scfsic(self):

783 a = 6.0

784 atm = Atoms("H2", positions=[(0, 0, 0), (0, 0, 0.737)], cell=(a, a, a))

785 atm.center()

786 calc = GPAW(mode="fd", xc="LDA-PZ-SIC",

787 eigensolver="rmm-diis", setups="hgh")

788 atm.calc = calc

789 atm.get_potential_energy()

790 return atm.calc

791

792 @gpwfile

793 def h_magmom(self):

794 a = 6.0

795 atm = Atoms("H", magmoms=[1.0], cell=(a, a, a))

796 atm.center()

797 calc = GPAW(mode="fd", xc="LDA-PZ-SIC",

798 eigensolver="rmm-diis", setups="hgh")

799 atm.calc = calc

800 atm.get_potential_energy()

801 return atm.calc

802

803 @gpwfile

804 def h_hess_num_pw(self):

805 calc = GPAW(

806 xc="PBE",

807 mode=PW(300, force_complex_dtype=False),

808 h=0.25,

809 convergence={

810 "energy": np.inf,

811 "eigenstates": np.inf,

812 "density": np.inf,

813 "minimum iterations": 1,

814 },

815 spinpol=False,

816 eigensolver=FDPWETDM(converge_unocc=True),

817 occupations={"name": "fixed-uniform"},

818 mixer={"backend": "no-mixing"},

819 nbands=2,

820 symmetry="off",

821 )

822 atoms = Atoms("H", positions=[[0, 0, 0]])

823 atoms.center(vacuum=5.0)

824 atoms.set_pbc(False)

825 atoms.calc = calc

826 atoms.get_potential_energy()

827 return atoms.calc

828

829 @gpwfile

830 def h2_break_ilcao(self):

831 atoms = Atoms('H2', positions=[(0, 0, 0), (0, 0, 2.0)])

832 atoms.center(vacuum=2.0)

833 atoms.set_pbc(False)

834 calc = GPAW(xc='PBE',

835 mode='lcao',

836 h=0.24,

837 basis='sz(dzp)',

838 spinpol=True,

839 eigensolver='etdm-lcao',

840 convergence={'density': 1.0e-2,

841 'eigenstates': 1.0e-2},

842 occupations={'name': 'fixed-uniform'},

843 mixer={'backend': 'no-mixing'},

844 nbands='nao',

845 symmetry='off')

846 atoms.calc = calc

847 atoms.get_potential_energy()

848 return atoms.calc

849

850 @gpwfile

851 def h_do_gdavid_lcao(self):

852 calc = GPAW(xc='PBE',

853 mode=LCAO(force_complex_dtype=True),

854 h=0.25,

855 basis='dz(dzp)',

856 spinpol=False,

857 eigensolver={'name': 'etdm-lcao',

858 'representation': 'u-invar'},

859 occupations={'name': 'fixed-uniform'},

860 mixer={'backend': 'no-mixing'},

861 nbands='nao',

862 symmetry='off',

863 )

864

865 atoms = Atoms('H', positions=[[0, 0, 0]])

866 atoms.center(vacuum=5.0)

867 atoms.set_pbc(False)

868 atoms.calc = calc

869 atoms.get_potential_energy()

870 return atoms.calc

871

872 @gpwfile

873 def h2_mom_do_pwh(self):

874 d = 1.4 * Bohr

875 h2 = Atoms('H2',

876 positions=[[-d / 2, 0, 0],

877 [d / 2, 0, 0]])

878 h2.center(vacuum=3)

879 calc = GPAW(mode=PW(300),

880 # h=0.3,

881 xc={'name': 'HSE06', 'backend': 'pw'},

882 eigensolver={'name': 'etdm-fdpw',

883 'converge_unocc': True},

884 mixer={'backend': 'no-mixing'},

885 occupations={'name': 'fixed-uniform'},

886 symmetry='off',

887 nbands=2,

888 convergence={'eigenstates': 4.0e-6})

889 h2.calc = calc

890 h2.get_potential_energy()

891 return h2.calc

892

893 @gpwfile

894 def h_hess_num_lcao(self):

895 calc = GPAW(xc='PBE',

896 mode=LCAO(force_complex_dtype=True),

897 h=0.25,

898 basis='dz(dzp)',

899 spinpol=False,

900 eigensolver={'name': 'etdm-lcao',

901 'representation': 'u-invar'},

902 occupations={'name': 'fixed-uniform'},

903 mixer={'backend': 'no-mixing'},

904 nbands='nao',

905 symmetry='off',

906 )

907 atoms = Atoms('H', positions=[[0, 0, 0]])

908 atoms.center(vacuum=5.0)

909 atoms.set_pbc(False)

910 atoms.calc = calc

911 atoms.get_potential_energy()

912 return atoms.calc

913

914 @gpwfile

915 def silicon_pdens_tool(self):

916 # used by response code's pdens tool test

917 pw = 200

918 kpts = 3

919 nbands = 8

920

921 a = 5.431

922 atoms = bulk('Si', 'diamond', a=a)

923

924 calc = GPAW(mode=PW(pw),

925 kpts=(kpts, kpts, kpts),

926 nbands=nbands,

927 convergence={'bands': -1},

928 xc='LDA',

929 occupations=FermiDirac(0.001))

930

931 atoms.calc = calc

932 atoms.get_potential_energy()

933 return calc

934

935 @gpwfile

936 def h_pw(self):

937 h = Atoms('H', magmoms=[1])

938 h.center(vacuum=4.0)

939 h.calc = GPAW(mode={'name': 'pw', 'ecut': 500},

940 txt=self.folder / 'h_pw.txt')

941 h.get_potential_energy()

942 return h.calc

943

944 @gpwfile

945 def h_chain(self):

946 from gpaw.new.ase_interface import GPAW

947 a = 2.5

948 k = 4

949 """Compare 2*H AFM cell with 1*H q=1/2 spin-spiral cell."""

950 h = Atoms('H',

951 magmoms=[1],

952 cell=[a, 0, 0],

953 pbc=[1, 0, 0])

954 h.center(vacuum=2.0, axis=(1, 2))

955 h.calc = GPAW(mode={'name': 'pw',

956 'ecut': 400,

957 'qspiral': [0.5, 0, 0]},

958 magmoms=[[1, 0, 0]],

959 symmetry='off',

960 kpts=(2 * k, 1, 1),

961 txt=self.folder / 'h_chain.txt')

962 h.get_potential_energy()

963 return h.calc

964

965 @gpwfile

966 def h_shift(self):

967 # Check for Hydrogen atom

968 atoms = Atoms('H', cell=(3 * np.eye(3)), pbc=True)

969

970 # Do a GS and save it

971 calc = GPAW(

972 mode=PW(600), symmetry={'point_group': False},

973 kpts={'size': (2, 2, 2)}, nbands=5, txt=None)

974 atoms.calc = calc

975 atoms.get_potential_energy()

976

977 return atoms.calc

978

979 @gpwfile

980 def h2_chain(self):

981 a = 2.5

982 k = 4

983 h2 = Atoms('H2',

984 [(0, 0, 0), (a, 0, 0)],

985 magmoms=[1, -1],

986 cell=[2 * a, 0, 0],

987 pbc=[1, 0, 0])

988 h2.center(vacuum=2.0, axis=(1, 2))

989 h2.calc = GPAW(mode={'name': 'pw',

990 'ecut': 400},

991 kpts=(k, 1, 1),

992 txt=self.folder / 'h2_chain.txt')

993 h2.get_potential_energy()

994 return h2.calc

995

996 @gpwfile

997 def h2_lcao_pair(self):

998 atoms = molecule('H2')

999 atoms.set_cell([6.4, 6.4, 6.4])

1000 atoms.center()

1001

1002 atoms.calc = GPAW(mode='lcao', occupations=FermiDirac(0.1),

1003 poissonsolver={'name': 'fd'},

1004 txt=self.folder / 'h2_lcao_pair.txt')

1005 atoms.get_potential_energy()

1006 return atoms.calc

1007

1008 @gpwfile

1009 def na_chain(self):

1010 a = 3

1011 atom = Atoms('Na',

1012 cell=[a, 0, 0],

1013 pbc=[1, 1, 1])

1014 atom.center(vacuum=1 * a, axis=(1, 2))

1015 atom.center()

1016 atoms = atom.repeat((2, 1, 1))

1017

1018 calc = GPAW(mode=PW(200),

1019 nbands=4,

1020 setups={'Na': '1'},

1021 txt=self.folder / 'na_chain.txt',

1022 kpts=(16, 2, 2))

1023

1024 atoms.calc = calc

1025 atoms.get_potential_energy()

1026 return calc

1027

1028 @gpwfile

1029 def n2_pw(self):

1030 N2 = molecule('N2')

1031 N2.center(vacuum=2.0)

1032

1033 N2.calc = GPAW(mode=PW(force_complex_dtype=True),

1034 xc='PBE',

1035 parallel={'domain': 1},

1036 eigensolver='rmm-diis',

1037 txt=self.folder / 'n2_pw.txt')

1038

1039 N2.get_potential_energy()

1040 N2.calc.diagonalize_full_hamiltonian(nbands=104, scalapack=True)

1041 return N2.calc

1042

1043 @gpwfile

1044 def n_pw(self):

1045 N2 = molecule('N2')

1046 N2.center(vacuum=2.0)

1047

1048 N = molecule('N')

1049 N.set_cell(N2.cell)

1050 N.center()

1051

1052 N.calc = GPAW(mode=PW(force_complex_dtype=True),

1053 xc='PBE',

1054 parallel={'domain': 1},

1055 eigensolver='rmm-diis',

1056 txt=self.folder / 'n_pw.txt')

1057 N.get_potential_energy()

1058 N.calc.diagonalize_full_hamiltonian(nbands=104, scalapack=True)

1059 return N.calc

1060

1061 @gpwfile

1062 def o2_pw(self):

1063 d = 1.1

1064 a = Atoms('O2', positions=[[0, 0, 0], [d, 0, 0]], magmoms=[1, 1])

1065 a.center(vacuum=4.0)

1066 a.calc = GPAW(mode={'name': 'pw', 'ecut': 800},

1067 txt=self.folder / 'o2_pw.txt')

1068 a.get_potential_energy()

1069 return a.calc

1070

1071 @gpwfile

1072 def Cu3Au_qna(self):

1073 ecut = 300

1074 kpts = (1, 1, 1)

1075

1076 QNA = {'alpha': 2.0,

1077 'name': 'QNA',

1078 'stencil': 1,

1079 'orbital_dependent': False,

1080 'parameters': {'Au': (0.125, 0.1), 'Cu': (0.0795, 0.005)},

1081 'setup_name': 'PBE',

1082 'type': 'qna-gga'}

1083

1084 atoms = L1_2(['Au', 'Cu'], latticeconstant=3.7)

1085 atoms[0].position[0] += 0.01 # Break symmetry already here

1086 calc = GPAW(mode=PW(ecut),

1087 eigensolver=Davidson(2),

1088 nbands='120%',

1089 mixer=Mixer(0.4, 7, 50.0),

1090 parallel=dict(domain=1),

1091 convergence={'density': 1e-4},

1092 xc=QNA,

1093 kpts=kpts,

1094 txt=self.folder / 'Cu3Au_qna.txt')

1095 atoms.calc = calc

1096 atoms.get_potential_energy()

1097 return atoms.calc

1098

1099 @gpwfile

1100 def co_mom(self):

1101 atoms = molecule('CO')

1102 atoms.center(vacuum=2)

1103

1104 calc = GPAW(mode='lcao',

1105 basis='dzp',

1106 nbands=7,

1107 h=0.24,

1108 xc='PBE',

1109 spinpol=True,

1110 symmetry='off',

1111 convergence={'energy': 100,

1112 'density': 1e-3},

1113 txt=self.folder / 'co_mom.txt')

1114

1115 atoms.calc = calc

1116 # Ground-state calculation

1117 atoms.get_potential_energy()

1118 return atoms.calc

1119

1120 @gpwfile

1121 def co_lcao(self):

1122 d = 1.1

1123 co = Atoms('CO', positions=[[0, 0, 0], [d, 0, 0]])

1124 co.center(vacuum=4.0)

1125 co.calc = GPAW(mode='lcao',

1126 txt=self.folder / 'co_lcao.txt')

1127 co.get_potential_energy()

1128 return co.calc

1129

1130 def _c2h4(self):

1131 d = 1.54

1132 h = 1.1

1133 x = d * (2 / 3)**0.5

1134 z = d / 3**0.5

1135 pe = Atoms('C2H4',

1136 positions=[[0, 0, 0],

1137 [x, 0, z],

1138 [0, -h * (2 / 3)**0.5, -h / 3**0.5],

1139 [0, h * (2 / 3)**0.5, -h / 3**0.5],

1140 [x, -h * (2 / 3)**0.5, z + h / 3**0.5],

1141 [x, h * (2 / 3)**0.5, z + h / 3**0.5]],

1142 cell=[2 * x, 0, 0],

1143 pbc=(1, 0, 0))

1144 pe.center(vacuum=2.0, axis=(1, 2))

1145 return pe

1146

1147 def c2h4_pw_nosym(self):

1148 pe = self._c2h4()

1149 pe.calc = GPAW(mode='pw',

1150 kpts=(3, 1, 1),

1151 symmetry='off',

1152 txt=self.folder / 'c2h4_pw_nosym.txt')

1153 pe.get_potential_energy()

1154 return pe.calc

1155

1156 @gpwfile

1157 def c6h12_pw(self):

1158 pe = self._c2h4()

1159 pe = pe.repeat((3, 1, 1))

1160 pe.calc = GPAW(mode='pw', txt=self.folder / 'c6h12_pw.txt')

1161 pe.get_potential_energy()

1162 return pe.calc

1163

1164 @gpwfile

1165 def h2o_lcao(self):

1166 atoms = molecule('H2O', cell=[8, 8, 8], pbc=1)

1167 atoms.center()

1168 atoms.calc = GPAW(mode='lcao', txt=self.folder / 'h2o_lcao.txt')

1169 atoms.get_potential_energy()

1170 return atoms.calc

1171

1172 @gpwfile

1173 def h2o_xas(self):

1174 setupname = 'h2o_xas_hch1s'

1175 self.generator2_setup(

1176 'O', 8, '2s,s,2p,p,d', [1.2], 1.0, None, 2,

1177 core_hole='1s,0.5',

1178 name=setupname)

1179

1180 a = 5.0

1181 H2O = self.h2o_maker(vacuum=None, cell=[a, a, a], pbc=False)

1182 calc = GPAW(

1183 txt=self.folder / 'h2o_xas.txt',

1184 mode='fd',

1185 nbands=10,

1186 h=0.2,

1187 setups={'O': 'h2o_xas_hch1s'},

1188 poissonsolver=FDPoissonSolver(use_charge_center=True),

1189 )

1190 H2O.calc = calc

1191 _ = H2O.get_potential_energy()

1192 return calc

1193

1194 @gpwfile

1195 def h20_lr2_nbands8(self):

1196 return self.h2o_nbands(

1197 {'poissonsolver': {'name': 'fd'}, 'nbands': 8, 'basis': 'dzp'})

1198

1199 @gpwfile

1200 def h20_lr2_nbands6(self):

1201 return self.h2o_nbands({'nbands': 6, 'basis': 'sz(dzp)'})

1202

1203 def h2o_nbands(self, params):

1204 atoms = molecule('H2O')

1205 atoms.center(vacuum=4)

1206

1207 calc = GPAW(h=0.4, **params, xc='LDA', mode='lcao',

1208 txt=self.folder / f'h2o_lr2_nbands{params["nbands"]}.out')

1209 atoms.calc = calc

1210 atoms.get_potential_energy()

1211

1212 return atoms.calc

1213

1214 @gpwfile

1215 def si_fd_ibz(self):

1216 si = bulk('Si', 'diamond', a=5.43)

1217 k = 3

1218 si.calc = GPAW(mode='fd', kpts=(k, k, k),

1219 txt=self.folder / 'si_fd_ibz.txt')

1220 si.get_potential_energy()

1221 return si.calc

1222

1223 @gpwfile

1224 def si_fd_bz(self):

1225 si = bulk('Si', 'diamond', a=5.43)

1226 k = 3

1227 si.calc = GPAW(mode='fd', kpts=(k, k, k,),

1228 symmetry={'point_group': False,

1229 'time_reversal': False},

1230 txt=self.folder / 'si_fd_bz.txt')

1231 si.get_potential_energy()

1232 return si.calc

1233

1234 @gpwfile

1235 def si_pw(self):

1236 si = bulk('Si')

1237 calc = GPAW(mode='pw',

1238 xc='LDA',

1239 occupations=FermiDirac(width=0.001),

1240 kpts={'size': (2, 2, 2), 'gamma': True},

1241 txt=self.folder / 'si_pw.txt')

1242 si.calc = calc

1243 si.get_potential_energy()

1244 return si.calc

1245

1246 @gpwfile

1247 def si_noisy_kpoints(self):

1248 # Test system for guarding against inconsistent kpoints as in #1178.

1249 from ase.calculators.calculator import kpts2kpts

1250 atoms = bulk('Si')

1251 kpts = kpts2kpts(kpts={'size': (2, 2, 2), 'gamma': True}, atoms=atoms)

1252

1253 # Error happened when qpoint was ~1e-17 yet was not considered gamma.

1254 # Add a bit of noise on purpose so we are sure to hit such a case,

1255 # even if the underlying implementation changes:

1256 kpts.kpts += np.linspace(1e-16, 1e-15, 24).reshape(8, 3)

1257

1258 calc = GPAW(mode='pw',

1259 xc='LDA',

1260 occupations=FermiDirac(width=0.001),

1261 kpts=kpts.kpts,

1262 txt=self.folder / 'si_noisy_kpoints.txt')

1263 atoms.calc = calc

1264 atoms.get_potential_energy()

1265 return calc

1266

1267 @gpwfile

1268 def si_pw_nbands10_converged(self):

1269 calc = GPAW(mode='pw',

1270 kpts={'size': (2, 2, 2), 'gamma': True},

1271 occupations=FermiDirac(0.01),

1272 nbands=10,

1273 symmetry='off',

1274 convergence={'bands': -4, 'density': 1e-7,

1275 'eigenstates': 1e-10})

1276

1277 atoms = bulk('Si', 'diamond', a=5.431)

1278 atoms.calc = calc

1279 atoms.get_potential_energy()

1280 return calc

1281

1282 @property

1283 def testing_setup_path(self):

1284 # Some calculations in gpwfile fixture like to use funny setups.

1285 # This is not so robust since the setups will be all jumbled.

1286 # We could improve the mechanism by programmatic naming/subfolders.

1287 return self.folder / 'setups'

1288

1289 def save_setup(self, setup):

1290 self.testing_setup_path.mkdir(parents=True, exist_ok=True)

1291 setup_file = self.testing_setup_path / setup.stdfilename

1292 if world.rank == 0:

1293 setup.write_xml(setup_file)

1294 world.barrier()

1295 return setup

1296

1297 def generate_setup(self, *args, **kwargs):

1298 from gpaw.test import gen

1299 setup = gen(*args, **kwargs, write_xml=False)

1300 self.save_setup(setup)

1301 return setup

1302

1303 def generator2_setup(self, *args, name, **kwargs):

1304 from gpaw.atom.generator2 import generate

1305 gen = generate(*args, **kwargs)

1306 setup = gen.make_paw_setup(name)

1307 self.save_setup(setup)

1308 return setup

1309

1310 @gpwfile

1311 def si_corehole_pw(self):

1312 # Generate setup for oxygen with half a core-hole:

1313 setupname = 'si_corehole_pw_hch1s'

1314 self.generate_setup('Si', name=setupname,

1315 corehole=(1, 0, 0.5), gpernode=30)

1316

1317 a = 2.6

1318 si = Atoms('Si', cell=(a, a, a), pbc=True)

1319

1320 calc = GPAW(mode='fd',

1321 txt=self.folder / 'si_corehole_pw.txt',

1322 h=0.25,

1323 occupations=FermiDirac(width=0.05),

1324 setups='si_corehole_pw_hch1s',

1325 convergence={'maximum iterations': 1})

1326 si.calc = calc

1327 _ = si.get_potential_energy()

1328 return si.calc

1329

1330 @gpwfile

1331 def si_corehole_sym_pw(self):

1332 setupname = 'si_corehole_sym_pw_hch1s'

1333 self.generate_setup('Si', name=setupname, corehole=(1, 0, 0.5),

1334 gpernode=30)

1335 return self.si_corehole_sym(sym={}, setupname=setupname)

1336

1337 @gpwfile

1338 def si_corehole_nosym_pw(self):

1339 setupname = 'si_corehole_sym_pw_hch1s'

1340 # XXX same setup as above, but we have it twice since caching

1341 # works per gpw file and not per setup

1342 self.generate_setup('Si', name=setupname, corehole=(1, 0, 0.5),

1343 gpernode=30)

1344 return self.si_corehole_sym(sym='off', setupname=setupname)

1345

1346 def si_corehole_sym(self, sym, setupname):

1347 tag = 'nosym' if sym == 'off' else 'sym'

1348

1349 a = 5.43095

1350 si_nonortho = Atoms(

1351 [Atom('Si', (0, 0, 0)), Atom('Si', (a / 4, a / 4, a / 4))],

1352 cell=[(a / 2, a / 2, 0), (a / 2, 0, a / 2), (0, a / 2, a / 2)],

1353 pbc=True,

1354 )

1355 # calculation with full symmetry

1356 calc = GPAW(

1357 mode='fd',

1358 txt=self.folder / f'si_corehole_{tag}_hch1s.txt',

1359 nbands=-10,

1360 h=0.25,

1361 kpts=(2, 2, 2),

1362 occupations=FermiDirac(width=0.05),

1363 setups={0: setupname},

1364 symmetry=sym

1365 )

1366 si_nonortho.calc = calc

1367 _ = si_nonortho.get_potential_energy()

1368 return calc

1369

1370 @gpwfile

1371 @with_band_cutoff(gpw='fancy_si_pw',

1372 band_cutoff=8) # 2 * (3s, 3p)

1373 def _fancy_si(self, *, band_cutoff, symmetry=None):

1374 if symmetry is None:

1375 symmetry = {}

1376 xc = 'LDA'

1377 kpts = 4

1378 pw = 300

1379 occw = 0.01

1380 conv = {'bands': band_cutoff + 1,

1381 'density': 1.e-8}

1382 atoms = bulk('Si')

1383 atoms.center()

1384

1385 tag = '_nosym' if symmetry == 'off' else ''

1386 atoms.calc = GPAW(

1387 xc=xc,

1388 mode=PW(pw),

1389 kpts={'size': (kpts, kpts, kpts), 'gamma': True},

1390 nbands=band_cutoff + 12, # + 2 * (3s, 3p),

1391 occupations=FermiDirac(occw),

1392 convergence=conv,

1393 txt=self.folder / f'fancy_si_pw{tag}.txt',

1394 symmetry=symmetry)

1395

1396 atoms.get_potential_energy()

1397 return atoms.calc

1398

1399 @gpwfile

1400 def fancy_si_pw(self):

1401 return self._fancy_si()

1402

1403 @gpwfile

1404 def fancy_si_pw_nosym(self):

1405 return self._fancy_si(symmetry='off')

1406

1407 @with_band_cutoff(gpw='sic_pw',

1408 band_cutoff=8) # (3s, 3p) + (2s, 2p)

1409 def _sic_pw(self, *, band_cutoff, spinpol=False):

1410 """Simple semi-conductor with broken inversion symmetry."""

1411 # Use the diamond crystal structure as blue print

1412 diamond = bulk('C', 'diamond')

1413 si = bulk('Si', 'diamond')

1414 # Break inversion symmetry by substituting one Si for C

1415 atoms = si.copy()

1416 atoms.symbols = 'CSi'

1417 # Scale the cell to the diamond/Si average

1418 cell_cv = (diamond.get_cell() + si.get_cell()) / 2.

1419 atoms.set_cell(cell_cv)

1420

1421 # Set up calculator

1422 tag = '_spinpol' if spinpol else ''

1423 atoms.calc = GPAW(

1424 mode=PW(400),

1425 xc='LDA',

1426 kpts={'size': (4, 4, 4)},

1427 symmetry={'point_group': False,

1428 'time_reversal': True},

1429 nbands=band_cutoff + 6,

1430 occupations=FermiDirac(0.001),

1431 convergence={'bands': band_cutoff + 1,

1432 'density': 1e-8},

1433 spinpol=spinpol,

1434 txt=self.folder / f'sic_pw{tag}.txt'

1435 )

1436

1437 atoms.get_potential_energy()

1438 return atoms.calc

1439

1440 @gpwfile

1441 def sic_pw(self):

1442 return self._sic_pw()

1443

1444 @gpwfile

1445 def sic_pw_spinpol(self):

1446 return self._sic_pw(spinpol=True)

1447

1448 @staticmethod

1449 def generate_si_systems():

1450 a = 5.43

1451 si1 = bulk('Si', 'diamond', a=a)

1452 si2 = si1.copy()

1453 si2.positions -= a / 8

1454 return [si1, si2]

1455

1456 def _si_gw(self, atoms, symm, name):

1457 atoms.calc = GPAW(mode=PW(250),

1458 eigensolver='rmm-diis',

1459 occupations=FermiDirac(0.01),

1460 symmetry=symm,

1461 kpts={'size': (2, 2, 2), 'gamma': True},

1462 convergence={'density': 1e-7},

1463 parallel={'domain': 1},

1464 txt=self.folder / name)

1465 atoms.get_potential_energy()

1466 scalapack = atoms.calc.wfs.bd.comm.size

1467 atoms.calc.diagonalize_full_hamiltonian(nbands=8, scalapack=scalapack)

1468 return atoms.calc

1469

1470 @gpwfile

1471 def c2_gw_more_bands(self):

1472 a = 3.567

1473 atoms = bulk('C', 'diamond', a=a)

1474 atoms.calc = GPAW(mode=PW(400),

1475 parallel={'domain': 1},

1476 kpts={'size': (2, 2, 2), 'gamma': True},

1477 xc='LDA',

1478 occupations=FermiDirac(0.001))

1479 atoms.get_potential_energy()

1480 atoms.calc.diagonalize_full_hamiltonian(nbands=128)

1481 return atoms.calc

1482

1483 @gpwfile

1484 def na_pw(self):

1485 from ase.build import bulk

1486

1487 blk = bulk('Na', 'bcc', a=4.23)

1488

1489 ecut = 350

1490 blk.calc = GPAW(mode=PW(ecut),

1491 basis='dzp',

1492 kpts={'size': (4, 4, 4), 'gamma': True},

1493 parallel={'domain': 1, 'band': 1},

1494 txt=self.folder / 'na_pw.txt',

1495 nbands=4,

1496 occupations=FermiDirac(0.01),

1497 setups={'Na': '1'})

1498 blk.get_potential_energy()

1499 blk.calc.diagonalize_full_hamiltonian(nbands=520)

1500 return blk.calc

1501

1502 @gpwfile

1503 def na2_tddft_poisson_sym(self):

1504 return self._na2_tddft(name='poisson_sym', basis='dzp', xc='LDA')

1505

1506 @gpwfile

1507 def na2_tddft_poisson(self):

1508 return self._na2_tddft(name='poisson', basis='dzp', xc='LDA')

1509

1510 @gpwfile

1511 def na2_tddft_dzp(self):

1512 return self._na2_tddft(name='dzp', basis='dzp', xc='LDA')

1513

1514 @gpwfile

1515 def na2_tddft_sz(self):

1516 return self._na2_tddft(name='sz', basis='sz(dzp)', xc='oldLDA')

1517

1518 def _na2_tddft(self, name, basis, xc):

1519 atoms = molecule('Na2')

1520 atoms.center(vacuum=4.0)

1521

1522 names = ['poisson_sym', 'poisson']

1523 poisson = PoissonSolver('fd', eps=1e-16) if name in names \

1524 else None

1525 symmetry = {'point_group': True if name == 'poisson_sym' else False}

1526

1527 calc = GPAW(nbands=2, h=0.4, setups=dict(Na='1'),

1528 basis=basis, mode='lcao', xc=xc,

1529 convergence={'density': 1e-8},

1530 poissonsolver=poisson,

1531 communicator=serial_comm if xc == 'oldLDA' else world,

1532 symmetry=symmetry,

1533 txt=self.folder / f'na2_tddft_{name}.out')

1534 atoms.calc = calc

1535 atoms.get_potential_energy()

1536

1537 return calc

1538

1539 @gpwfile

1540 def na2_fd(self):

1541 """Sodium dimer, Na2."""

1542 d = 1.5

1543 atoms = Atoms(symbols='Na2',

1544 positions=[(0, 0, d),

1545 (0, 0, -d)],

1546 pbc=False)

1547

1548 atoms.center(vacuum=6.0)

1549 # Larger grid spacing, LDA is ok

1550 gs_calc = GPAW(mode='fd',

1551 txt=self.folder / 'na2_fd.txt',

1552 nbands=1, h=0.35, xc='LDA',

1553 setups={'Na': '1'},

1554 symmetry={'point_group': False})

1555 atoms.calc = gs_calc

1556 atoms.get_potential_energy()

1557 return atoms.calc

1558

1559 @gpwfile

1560 def na2_fd_with_sym(self):

1561 """Sodium dimer, Na2."""

1562 d = 1.5

1563 atoms = Atoms(symbols='Na2',

1564 positions=[(0, 0, d),

1565 (0, 0, -d)],

1566 pbc=False)

1567

1568 atoms.center(vacuum=6.0)

1569 # Larger grid spacing, LDA is ok

1570 gs_calc = GPAW(mode='fd', nbands=1, h=0.35, xc='LDA',

1571 txt=self.folder / 'na2_fd_with_sym.txt',

1572 setups={'Na': '1'})

1573 atoms.calc = gs_calc

1574 atoms.get_potential_energy()

1575 return atoms.calc

1576

1577 @gpwfile

1578 def na2_isolated(self):

1579 # Permittivity file

1580 if world.rank == 0:

1581 fo = open(self.folder / 'ed.txt', 'w')

1582 fo.writelines(['1.20 0.20 25.0'])

1583 fo.close()

1584 world.barrier()

1585

1586 from gpaw.fdtd.poisson_fdtd import FDTDPoissonSolver

1587 from gpaw.fdtd.polarizable_material import (

1588 PermittivityPlus,

1589 PolarizableMaterial,

1590 PolarizableSphere,

1591 )

1592

1593 # Whole simulation cell (Angstroms)

1594 large_cell = [20, 20, 30]

1595

1596 # Quantum subsystem

1597 atom_center = np.array([10.0, 10.0, 20.0])

1598 atoms = Atoms('Na2', [atom_center + [0.0, 0.0, -1.50],

1599 atom_center + [0.0, 0.0, +1.50]])

1600

1601 # Classical subsystem

1602 classical_material = PolarizableMaterial()

1603 sphere_center = np.array([10.0, 10.0, 10.0])

1604 classical_material.add_component(

1605 PolarizableSphere(

1606 permittivity=PermittivityPlus(self.folder / 'ed.txt'),

1607 center=sphere_center,

1608 radius=5.0

1609 )

1610 )

1611

1612 # Accuracy

1613 energy_eps = 0.0005

1614 density_eps = 1e-6

1615 poisson_eps = 1e-12

1616

1617 # Combined Poisson solver

1618 poissonsolver = FDTDPoissonSolver(

1619 classical_material=classical_material,

1620 eps=poisson_eps,

1621 qm_spacing=0.40,

1622 cl_spacing=0.40 * 4,

1623 cell=large_cell,

1624 remove_moments=(1, 4),

1625 communicator=world,

1626 potential_coupler='Refiner',

1627 )

1628 poissonsolver.set_calculation_mode('iterate')

1629

1630 # Combined system

1631 atoms.set_cell(large_cell)

1632 atoms, qm_spacing, gpts = poissonsolver.cut_cell(atoms, vacuum=2.50)

1633

1634 # Initialize GPAW

1635 gs_calc = GPAW(mode='fd',

1636 txt=self.folder / 'na2_isolated.txt',

1637 gpts=gpts,

1638 eigensolver='cg',

1639 nbands=-1,

1640 poissonsolver=poissonsolver,

1641 symmetry={'point_group': False},

1642 convergence={'energy': energy_eps,

1643 'density': density_eps})

1644 atoms.calc = gs_calc

1645

1646 # Ground state

1647 atoms.get_potential_energy()

1648

1649 return gs_calc

1650

1651 @gpwfile

1652 def na3_pw_restart(self):

1653 params = dict(mode=PW(200), convergence={

1654 'eigenstates': 1.24, 'energy': 2e-1, 'density': 1e-1})

1655 return self._na3_restart(params=params)

1656

1657 @gpwfile

1658 def na3_fd_restart(self):

1659 params = dict(mode='fd', h=0.30, convergence={

1660 'eigenstates': 1.24, 'energy': 2e-1, 'density': 1e-1})

1661 return self._na3_restart(params=params)

1662

1663 @gpwfile

1664 def na3_fd_kp_restart(self):

1665 params = dict(mode='fd', h=0.30, kpts=(1, 1, 3), convergence={

1666 'eigenstates': 1.24, 'energy': 2e-1, 'density': 1e-1})

1667 return self._na3_restart(params=params)

1668

1669 @gpwfile

1670 def na3_fd_density_restart(self):

1671 params = dict(mode='fd', h=0.30, convergence={

1672 'eigenstates': 1.e-3, 'energy': 2e-1, 'density': 1e-1})

1673 return self._na3_restart(params=params)

1674

1675 def _na3_restart(self, params):

1676 d = 3.0

1677 atoms = Atoms('Na3',

1678 positions=[(0, 0, 0),

1679 (0, 0, d),

1680 (0, d * sqrt(3 / 4), d / 2)],

1681 magmoms=[1.0, 1.0, 1.0],

1682 cell=(3.5, 3.5, 4 + 2 / 3),

1683 pbc=True)

1684

1685 atoms.calc = GPAW(nbands=3,

1686 setups={'Na': '1'},

1687 **params)

1688 atoms.get_potential_energy()

1689

1690 return atoms.calc

1691

1692 @gpwfile

1693 def sih4_xc_gllbsc_lcao(self):

1694 return self._sih4_gllbsc(mode='lcao', basis='dzp')

1695

1696 @gpwfile

1697 def sih4_xc_gllbsc_fd(self):

1698 return self._sih4_gllbsc(mode='fd', basis={})

1699

1700 def _sih4_gllbsc(self, mode, basis):

1701 atoms = molecule('SiH4')

1702 atoms.center(vacuum=4.0)

1703

1704 # Ground-state calculation

1705 calc = GPAW(mode=mode, nbands=7, h=0.4,

1706 convergence={'density': 1e-8},

1707 xc='GLLBSC',

1708 basis=basis,

1709 symmetry={'point_group': False},

1710 txt=self.folder / f'sih4_xc_gllbsc_{mode}.txt')

1711 atoms.calc = calc

1712 atoms.get_potential_energy()

1713 return atoms.calc

1714

1715 @gpwfile

1716 def nacl_nospin(self):

1717 return self._nacl_mol(spinpol=False)

1718

1719 @gpwfile

1720 def nacl_spin(self):

1721 return self._nacl_mol(spinpol=True)

1722

1723 def _nacl_mol(self, spinpol):

1724 atoms = molecule('NaCl')

1725 atoms.center(vacuum=4.0)

1726 calc = GPAW(nbands=6,

1727 h=0.4,

1728 setups=dict(Na='1'),

1729 basis='dzp',

1730 mode='lcao',

1731 convergence={'density': 1e-8},

1732 spinpol=spinpol,

1733 communicator=world,

1734 symmetry={'point_group': False},

1735 txt=self.folder / 'gs.out')

1736 atoms.calc = calc

1737 atoms.get_potential_energy()

1738 return atoms.calc

1739

1740 @gpwfile

1741 def nacl_fd(self):

1742 d = 4.0

1743 atoms = Atoms('NaCl', [(0, 0, 0), (0, 0, d)])

1744 atoms.center(vacuum=4.5)

1745

1746 gs_calc = GPAW(

1747 txt=self.folder / 'nacl_fd.txt',

1748 mode='fd', nbands=4, # eigensolver='cg',

1749 gpts=(32, 32, 44), xc='LDA', symmetry={'point_group': False},

1750 setups={'Na': '1'})

1751 atoms.calc = gs_calc

1752 atoms.get_potential_energy()

1753 return atoms.calc

1754

1755 @gpwfile

1756 def bn_pw(self):

1757 atoms = bulk('BN', 'zincblende', a=3.615)

1758 atoms.calc = GPAW(mode=PW(400),

1759 kpts={'size': (2, 2, 2), 'gamma': True},

1760 nbands=12,

1761 convergence={'bands': 9},

1762 occupations=FermiDirac(0.001),

1763 txt=self.folder / 'bn_pw.txt')

1764 atoms.get_potential_energy()

1765 return atoms.calc

1766

1767 def _hbn_pw(self, symmetry={}):

1768 atoms = Graphene(symbol='B',

1769 latticeconstant={'a': 2.5, 'c': 1.0},

1770 size=(1, 1, 1))

1771 atoms[0].symbol = 'N'

1772 atoms.pbc = (1, 1, 0)

1773 atoms.center(axis=2, vacuum=3.0)

1774 atoms.calc = GPAW(txt=self.folder / 'hbn_pw.txt',

1775 mode=PW(400),

1776 xc='LDA',

1777 nbands=50,

1778 occupations=FermiDirac(0.001),

1779 parallel={'domain': 1},

1780 convergence={'bands': 26},

1781 kpts={'size': (3, 3, 1), 'gamma': True},

1782 symmetry=symmetry)

1783 atoms.get_potential_energy()

1784 return atoms.calc

1785

1786 @gpwfile

1787 def hbn_pw(self):

1788 return self._hbn_pw()

1789

1790 @gpwfile

1791 def hbn_pw_nosym(self):

1792 return self._hbn_pw(symmetry='off')

1793

1794 @gpwfile

1795 def graphene_pw(self):

1796 from ase.lattice.hexagonal import Graphene

1797 atoms = Graphene(symbol='C',

1798 latticeconstant={'a': 2.45, 'c': 1.0},

1799 size=(1, 1, 1))

1800 atoms.pbc = (1, 1, 0)

1801 atoms.center(axis=2, vacuum=4.0)

1802 ecut = 250

1803 nkpts = 6

1804 atoms.calc = GPAW(mode=PW(ecut),

1805 kpts={'size': (nkpts, nkpts, 1), 'gamma': True},

1806 nbands=len(atoms) * 6,

1807 txt=self.folder / 'graphene_pw.txt')

1808 atoms.get_potential_energy()

1809 return atoms.calc

1810

1811 @gpwfile

1812 def i2sb2_pw_nosym(self):

1813 # Structure from c2db

1814 atoms = Atoms('I2Sb2',

1815 positions=[[0.02437357, 0.05048655, 6.11612164],

1816 [0.02524896, 3.07135573, 11.64646853],

1817 [0.02717742, 0.01556495, 8.89278807],

1818 [0.02841809, 3.10675382, 8.86983839]],

1819 cell=[[5.055642258802973, -9.89475498615942e-15, 0.0],

1820 [-2.5278211265136266, 4.731999711338355, 0.0],

1821 [3.38028806436979e-15, 0.0, 18.85580293064]],

1822 pbc=(1, 1, 0))

1823 atoms.calc = GPAW(mode=PW(250),

1824 xc='PBE',

1825 kpts={'size': (6, 6, 1), 'gamma': True},

1826 txt=self.folder / 'i2sb2_pw_nosym.txt',

1827 symmetry='off')

1828

1829 atoms.get_potential_energy()

1830 return atoms.calc

1831

1832 @with_band_cutoff(gpw='bi2i6_pw',

1833 band_cutoff=36)

1834 def _bi2i6(self, *, band_cutoff, symmetry=None):

1835 if symmetry is None:

1836 symmetry = {}

1837 positions = [[4.13843656, 2.38932746, 9.36037077],

1838 [0.00000000, 4.77865492, 9.36034750],

1839 [3.89827619, 0.00000000, 7.33713295],

1840 [2.18929748, 3.79197674, 7.33713295],

1841 [-1.94913711, 3.37600678, 7.33713295],

1842 [3.89827619, 0.00000000, 11.3835853],

1843 [2.18929961, 3.79197551, 11.3835853],

1844 [-1.94913924, 3.37600555, 11.3835853]]

1845 cell = [[8.276873113486648, 0.0, 0.0],

1846 [-4.138436556743325, 7.167982380179831, 0.0],

1847 [0.0, 0.0, 0.0]]

1848 pbc = [True, True, False]

1849 atoms = Atoms('Bi2I6',

1850 positions=positions,

1851 cell=cell,

1852 pbc=pbc)

1853 atoms.center(vacuum=3.0, axis=2)

1854

1855 ecut = 120

1856 nkpts = 4

1857 conv = {'bands': band_cutoff + 1,

1858 'density': 1.e-4}

1859

1860 tag = '_nosym' if symmetry == 'off' else ''

1861 atoms.calc = GPAW(mode=PW(ecut),

1862 xc='LDA',

1863 kpts={'size': (nkpts, nkpts, 1), 'gamma': True},

1864 occupations=FermiDirac(0.01),

1865 mixer={'beta': 0.5},

1866 convergence=conv,

1867 nbands=band_cutoff + 9,

1868 txt=self.folder / f'bi2i6_pw{tag}.txt',

1869 symmetry=symmetry)

1870

1871 atoms.get_potential_energy()

1872 return atoms.calc

1873

1874 @gpwfile

1875 def bi2i6_pw(self):

1876 return self._bi2i6()

1877

1878 @gpwfile

1879 def bi2i6_pw_nosym(self):

1880 return self._bi2i6(symmetry='off')

1881

1882 def _mos2(self, symmetry=None):

1883 if symmetry is None:

1884 symmetry = {}

1885 from ase.build import mx2

1886 atoms = mx2(formula='MoS2', kind='2H', a=3.184, thickness=3.127,

1887 size=(1, 1, 1), vacuum=5)

1888 atoms.pbc = (1, 1, 0)

1889 ecut = 250

1890 nkpts = 6

1891 tag = '_nosym' if symmetry == 'off' else ''

1892 atoms.calc = GPAW(mode=PW(ecut),

1893 xc='LDA',

1894 kpts={'size': (nkpts, nkpts, 1), 'gamma': True},

1895 occupations=FermiDirac(0.01),

1896 txt=self.folder / f'mos2_pw{tag}.txt',

1897 symmetry=symmetry)

1898

1899 atoms.get_potential_energy()

1900 return atoms.calc

1901

1902 @gpwfile

1903 def mos2_pw(self):

1904 return self._mos2()

1905

1906 @gpwfile

1907 def mos2_pw_nosym(self):

1908 return self._mos2(symmetry='off')

1909

1910 @gpwfile

1911 def mos2_5x5_pw(self):

1912 calc = GPAW(mode=PW(180),

1913 xc='PBE',

1914 nbands='nao',

1915 setups={'Mo': '6'},

1916 occupations=FermiDirac(0.001),

1917 convergence={'bands': -5},

1918 kpts=(5, 5, 1))

1919

1920 from ase.build import mx2

1921 layer = mx2(formula='MoS2', kind='2H', a=3.1604, thickness=3.172,

1922 size=(1, 1, 1), vacuum=3.414)

1923 layer.pbc = (1, 1, 0)

1924 layer.calc = calc

1925 layer.get_potential_energy()

1926 return layer.calc

1927

1928 @with_band_cutoff(gpw='p4_pw',

1929 band_cutoff=40)

1930 def _p4(self, band_cutoff, spinpol=False):

1931 atoms = Atoms('P4', positions=[[0.03948480, -0.00027057, 7.49990646],

1932 [0.86217564, -0.00026338, 9.60988536],

1933 [2.35547782, 1.65277230, 9.60988532],

1934 [3.17816857, 1.65277948, 7.49990643]],

1935 cell=[4.63138807675, 3.306178252090, 17.10979291],

1936 pbc=[True, True, False])

1937 atoms.center(vacuum=1.5, axis=2)

1938 tag = '_spinpol' if spinpol else ''

1939 nkpts = 2

1940 atoms.calc = GPAW(mode=PW(250),

1941 xc='LDA', spinpol=spinpol,

1942 kpts={'size': (nkpts, nkpts, 1), 'gamma': True},

1943 occupations={'width': 0},

1944 nbands=band_cutoff + 10,

1945 convergence={'bands': band_cutoff + 1},

1946 txt=self.folder / f'p4_pw{tag}.txt')

1947 atoms.get_potential_energy()

1948 return atoms.calc

1949

1950 @gpwfile

1951 def p4_pw(self):

1952 return self._p4()

1953

1954 @gpwfile

1955 def p4_pw_spinpol(self):

1956 return self._p4(spinpol=True)

1957

1958 def _ni_pw_kpts333(self, setups={'Ni': '10'}):

1959 from ase.dft.kpoints import monkhorst_pack

1960 # from gpaw.mpi import serial_comm

1961 Ni = bulk('Ni', 'fcc')

1962 Ni.set_initial_magnetic_moments([0.7])

1963

1964 kpts = monkhorst_pack((3, 3, 3))

1965

1966 calc = GPAW(mode='pw',

1967 txt=self.folder / 'ni_pw_kpts333.txt',

1968 kpts=kpts,

1969 occupations=FermiDirac(0.001),

1970 setups=setups,

1971 parallel=dict(domain=1), # >1 fails on 8 cores

1972 # communicator=serial_comm

1973 )

1974

1975 Ni.calc = calc

1976 Ni.get_potential_energy()

1977 calc.diagonalize_full_hamiltonian()

1978 return calc

1979

1980 @gpwfile

1981 def ni_pw(self):

1982 return self._ni_pw_kpts333(setups={})

1983

1984 @gpwfile

1985 def ni_pw_kpts333(self):

1986 return self._ni_pw_kpts333()

1987

1988 @gpwfile

1989 def c_pw(self):

1990 atoms = bulk('C')

1991 atoms.center()

1992 calc = GPAW(mode=PW(150),

1993 txt=self.folder / 'c_pw.txt',

1994 convergence={'bands': 6},

1995 nbands=12,

1996 kpts={'gamma': True, 'size': (2, 2, 2)},

1997 xc='LDA')

1998

1999 atoms.calc = calc

2000 atoms.get_potential_energy()

2001 return atoms.calc

2002

2003 def _nicl2_pw(self, vacuum=3.0, identifier='', ecut=285, **kwargs):

2004 from ase.build import mx2

2005

2006 # Define input parameters

2007 kpts = 6

2008 occw = 0.01

2009 conv = {'density': 1.e-3}

2010

2011 a = 3.502

2012 thickness = 2.617

2013 mm = 2.0

2014

2015 # Set up atoms

2016 atoms = mx2(formula='NiCl2', kind='1T', a=a,

2017 thickness=thickness, vacuum=vacuum)

2018 atoms.set_initial_magnetic_moments([mm, 0.0, 0.0])

2019 # Use pbc to allow for real-space density interpolation

2020 atoms.pbc = True

2021

2022 dct = dict(

2023 mixer={'beta': 0.75, 'nmaxold': 8, 'weight': 100.0},

2024 mode=PW(ecut,

2025 # Interpolate the density in real-space

2026 interpolation=3),

2027 kpts={'size': (kpts, kpts, 1), 'gamma': True},

2028 occupations=FermiDirac(occw),

2029 convergence=conv,

2030 txt=self.folder / f'nicl2_pw{identifier}.txt')

2031 dct.update(kwargs)

2032 atoms.calc = GPAW(**dct)

2033

2034 atoms.get_potential_energy()

2035

2036 return atoms.calc

2037

2038 @gpwfile

2039 def nicl2_pw(self):

2040 return self._nicl2_pw(vacuum=3.0, ecut=300.0)

2041

2042 @gpwfile

2043 def nicl2_pw_evac(self):

2044 return self._nicl2_pw(vacuum=10.0, identifier='_evac')

2045

2046 @gpwfile

2047 def Tl_box_pw(self):

2048 Tl = Atoms('Tl',

2049 cell=[5, 5, 5],

2050 scaled_positions=[[0.5, 0.5, 0.5]],

2051 pbc=False)

2052

2053 Tl.calc = GPAWNew(

2054 mode={'name': 'pw', 'ecut': 300},

2055 xc='LDA',

2056 occupations={'name': 'fermi-dirac', 'width': 0.01},

2057 symmetry='off',

2058 convergence={'density': 1e-6},

2059 parallel={'domain': 1, 'band': 1},

2060 magmoms=[[0, 0, 0.5]],

2061 soc=True,

2062 txt=self.folder / 'Tl_box_pw.txt')

2063 Tl.get_potential_energy()

2064 return Tl.calc

2065

2066 @with_band_cutoff(gpw='v2br4_pw',

2067 band_cutoff=28) # V(4s,3d) = 6, Br(4s,4p) = 4

2068 def _v2br4(self, *, band_cutoff, symmetry=None):

2069 from ase.build import mx2

2070

2071 if symmetry is None:

2072 symmetry = {}

2073

2074 # Define input parameters

2075 xc = 'LDA'

2076 kpts = 4

2077 pw = 200

2078 occw = 0.01

2079 conv = {'density': 1.e-4,

2080 'bands': band_cutoff + 1}

2081

2082 a = 3.840

2083 thickness = 2.897

2084 vacuum = 3.0

2085 mm = 3.0

2086

2087 # Set up atoms

2088 atoms = mx2(formula='VBr2', kind='1T', a=a,

2089 thickness=thickness, vacuum=vacuum)

2090 atoms = atoms.repeat((1, 2, 1))

2091 atoms.set_initial_magnetic_moments([mm, 0.0, 0.0, -mm, 0.0, 0.0])

2092 # Use pbc to allow for real-space density interpolation

2093 atoms.pbc = True

2094

2095 # Set up calculator

2096 tag = '_nosym' if symmetry == 'off' else ''

2097 atoms.calc = GPAW(

2098 xc=xc,

2099 mode=PW(pw,

2100 # Interpolate the density in real-space

2101 interpolation=3),

2102 kpts={'size': (kpts, kpts // 2, 1), 'gamma': True},

2103 mixer={'beta': 0.5},

2104 setups={'V': '5'},

2105 nbands=band_cutoff + 12,

2106 occupations=FermiDirac(occw),

2107 convergence=conv,

2108 symmetry=symmetry,

2109 txt=self.folder / f'v2br4_pw{tag}.txt')

2110

2111 atoms.get_potential_energy()

2112

2113 return atoms.calc

2114

2115 @gpwfile

2116 def v2br4_pw(self):

2117 return self._v2br4()

2118

2119 @gpwfile

2120 def v2br4_pw_nosym(self):

2121 return self._v2br4(symmetry='off')

2122

2123 @with_band_cutoff(gpw='fe_pw',

2124 band_cutoff=9) # 4s, 4p, 3d = 9

2125 def _fe(self, *, band_cutoff, symmetry=None):

2126 if symmetry is None:

2127 symmetry = {}

2128 """See also the fe_fixture_test.py test."""

2129 xc = 'LDA'

2130 kpts = 4

2131 pw = 300

2132 occw = 0.01

2133 conv = {'bands': band_cutoff + 1,

2134 'density': 1.e-8}

2135 a = 2.867

2136 mm = 2.21

2137 atoms = bulk('Fe', 'bcc', a=a)

2138 # It is necessary to rattle the atoms to make sure that all tests pass

2139 # on all machines - see https://gitlab.com/gpaw/gpaw/-/issues/1397

2140 atoms.rattle(0.01, seed=42)

2141 atoms.set_initial_magnetic_moments([mm])

2142 atoms.center()

2143 tag = '_nosym' if symmetry == 'off' else ''

2144

2145 atoms.calc = GPAW(

2146 xc=xc,

2147 mode=PW(pw),

2148 kpts={'size': (kpts, kpts, kpts)},

2149 nbands=band_cutoff + 9,

2150 occupations=FermiDirac(occw),

2151 convergence=conv,

2152 txt=self.folder / f'fe_pw{tag}.txt',

2153 symmetry=symmetry)

2154

2155 atoms.get_potential_energy()

2156 return atoms.calc

2157

2158 @gpwfile

2159 def fe_pw(self):

2160 return self._fe()

2161

2162 @gpwfile

2163 def fe_pw_nosym(self):

2164 return self._fe(symmetry='off')

2165

2166 @with_band_cutoff(gpw='co_pw',

2167 band_cutoff=14) # 2 * (4s + 3d)

2168 def _co(self, *, band_cutoff, symmetry=None):

2169 if symmetry is None:

2170 symmetry = {}

2171 # ---------- Inputs ---------- #

2172

2173 # Atomic configuration

2174 a = 2.5071

2175 c = 4.0695

2176 mm = 1.6

2177 atoms = bulk('Co', 'hcp', a=a, c=c)

2178 atoms.set_initial_magnetic_moments([mm, mm])

2179 atoms.center()

2180

2181 # Ground state parameters

2182 xc = 'LDA'

2183 occw = 0.01

2184 ebands = 2 * 2 # extra bands for ground state calculation

2185 pw = 200

2186 conv = {'density': 1e-8,

2187 'forces': 1e-8,

2188 'bands': band_cutoff + 1}

2189

2190 # ---------- Calculation ---------- #

2191

2192 tag = '_nosym' if symmetry == 'off' else ''

2193 atoms.calc = GPAW(xc=xc,

2194 mode=PW(pw),

2195 kpts={'size': (4, 4, 4), 'gamma': True},

2196 occupations=FermiDirac(occw),

2197 convergence=conv,

2198 nbands=band_cutoff + ebands,

2199 symmetry=symmetry,

2200 txt=self.folder / f'co_pw{tag}.txt')

2201

2202 atoms.get_potential_energy()

2203 return atoms.calc

2204

2205 @gpwfile

2206 def co_pw(self):

2207 return self._co()

2208

2209 @gpwfile

2210 def co_pw_nosym(self):

2211 return self._co(symmetry='off')

2212

2213 @with_band_cutoff(gpw='srvo3_pw',

2214 band_cutoff=20)

2215 def _srvo3(self, *, band_cutoff, symmetry=None):

2216 if symmetry is None:

2217 symmetry = {}

2218

2219 nk = 3

2220 cell = bulk('V', 'sc', a=3.901).cell

2221 atoms = Atoms('SrVO3', cell=cell, pbc=True,

2222 scaled_positions=((0.5, 0.5, 0.5),

2223 (0, 0, 0),

2224 (0, 0.5, 0),

2225 (0, 0, 0.5),

2226 (0.5, 0, 0)))

2227 # Ground state parameters

2228 xc = 'LDA'

2229 occw = 0.01

2230 ebands = 10 # extra bands for ground state calculation

2231 pw = 200

2232 conv = {'density': 1e-8,

2233 'bands': band_cutoff + 1}

2234

2235 # ---------- Calculation ---------- #

2236

2237 tag = '_nosym' if symmetry == 'off' else ''

2238 atoms.calc = GPAW(xc=xc,

2239 mode=PW(pw),

2240 kpts={'size': (nk, nk, nk), 'gamma': True},

2241 occupations=FermiDirac(occw),

2242 convergence=conv,

2243 nbands=band_cutoff + ebands,

2244 symmetry=symmetry,

2245 txt=self.folder / f'srvo3_pw{tag}.txt')

2246

2247 atoms.get_potential_energy()

2248 return atoms.calc

2249

2250 @gpwfile

2251 def srvo3_pw(self):

2252 return self._srvo3()

2253

2254 @gpwfile

2255 def srvo3_pw_nosym(self):

2256 return self._srvo3(symmetry='off')

2257

2258 @with_band_cutoff(gpw='al_pw',

2259 band_cutoff=10) # 3s, 3p, 4s, 3d

2260 def _al(self, *, band_cutoff, symmetry=None):

2261 if symmetry is None:

2262 symmetry = {}

2263 xc = 'LDA'

2264 kpts = 4

2265 pw = 300

2266 occw = 0.01

2267 conv = {'bands': band_cutoff + 1,

2268 'density': 1.e-8}

2269 a = 4.043

2270 atoms = bulk('Al', 'fcc', a=a)

2271 atoms.center()

2272 tag = '_nosym' if symmetry == 'off' else ''

2273

2274 atoms.calc = GPAW(

2275 xc=xc,

2276 mode=PW(pw),

2277 kpts={'size': (kpts, kpts, kpts), 'gamma': True},

2278 nbands=band_cutoff + 4, # + 4p, 5s

2279 occupations=FermiDirac(occw),

2280 convergence=conv,

2281 txt=self.folder / f'al_pw{tag}.txt',

2282 symmetry=symmetry)

2283

2284 atoms.get_potential_energy()

2285 return atoms.calc

2286

2287 @gpwfile

2288 def al_pw(self):

2289 return self._al()

2290

2291 @gpwfile

2292 def al_pw_nosym(self):

2293 return self._al(symmetry='off')

2294

2295 @gpwfile

2296 def bse_al(self):

2297 a = 4.043

2298 atoms = bulk('Al', 'fcc', a=a)

2299 calc = GPAW(mode='pw',

2300 txt=self.folder / 'bse_al.txt',

2301 kpts={'size': (4, 4, 4), 'gamma': True},

2302 xc='LDA',

2303 nbands=4,

2304 convergence={'bands': 'all'})

2305

2306 atoms.calc = calc

2307 atoms.get_potential_energy()

2308 return atoms.calc

2309

2310 @gpwfile

2311 def ag_plusU_pw(self):

2312 xc = 'LDA'

2313 kpts = 2

2314 nbands = 6

2315 pw = 300

2316 occw = 0.01

2317 conv = {'bands': nbands,

2318 'density': 1e-12}

2319 a = 4.07

2320 atoms = bulk('Ag', 'fcc', a=a)

2321 atoms.center()

2322

2323 atoms.calc = GPAW(

2324 xc=xc,

2325 mode=PW(pw),

2326 kpts={'size': (kpts, kpts, kpts), 'gamma': True},

2327 setups={'Ag': '11:d,2.0,0'},

2328 nbands=nbands,

2329 occupations=FermiDirac(occw),

2330 convergence=conv,

2331 parallel={'domain': 1},

2332 txt=self.folder / 'ag_plusU_pw.txt')

2333

2334 atoms.get_potential_energy()

2335

2336 atoms.calc.diagonalize_full_hamiltonian()

2337

2338 return atoms.calc

2339

2340 @gpwfile

2341 def gaas_pw_nosym(self):

2342 return self._gaas(symmetry='off')

2343

2344 @gpwfile

2345 def gaas_pw(self):

2346 return self._gaas()

2347

2348 @with_band_cutoff(gpw='gaas_pw',

2349 band_cutoff=8)

2350 def _gaas(self, *, band_cutoff, symmetry=None):

2351 if symmetry is None:

2352 symmetry = {}

2353 nk = 4

2354 cell = bulk('Ga', 'fcc', a=5.68).cell

2355 atoms = Atoms('GaAs', cell=cell, pbc=True,

2356 scaled_positions=((0, 0, 0), (0.25, 0.25, 0.25)))

2357 tag = '_nosym' if symmetry == 'off' else ''

2358 conv = {'bands': band_cutoff + 1,

2359 'density': 1.e-8}

2360

2361 calc = GPAW(mode=PW(400),

2362 xc='LDA',

2363 occupations=FermiDirac(width=0.01),

2364 convergence=conv,

2365 nbands=band_cutoff + 1,

2366 kpts={'size': (nk, nk, nk), 'gamma': True},

2367 txt=self.folder / f'gs_GaAs{tag}.txt',

2368 symmetry=symmetry)

2369

2370 atoms.calc = calc

2371 atoms.get_potential_energy()

2372 return atoms.calc

2373

2374 @gpwfile

2375 def h_pw280_fulldiag(self):

2376 return self._pw_280_fulldiag(Atoms('H'), hund=True, nbands=4)

2377

2378 @gpwfile

2379 def h2_pw280_fulldiag(self):

2380 return self._pw_280_fulldiag(

2381 Atoms('H2', [(0, 0, 0), (0, 0, 0.7413)]), nbands=8)

2382

2383 def _pw_280_fulldiag(self, atoms, **kwargs):

2384 atoms.set_pbc(True)

2385 atoms.set_cell((2., 2., 3.))

2386 atoms.center()

2387 calc = GPAW(mode=PW(280, force_complex_dtype=True),

2388 txt=self.folder / f'{atoms.symbols}_pw_280_fulldiag.txt',

2389 xc='LDA',

2390 basis='dzp',

2391 parallel={'domain': 1},

2392 convergence={'density': 1.e-6},

2393 **kwargs)

2394 atoms.calc = calc

2395 atoms.get_potential_energy()

2396 calc.diagonalize_full_hamiltonian(nbands=80)

2397 return calc

2398

2399 @gpwfile

2400 def fe_pw_distorted(self):

2401 xc = 'revTPSS'

2402 m = [2.9]

2403 fe = bulk('Fe')

2404 fe.set_initial_magnetic_moments(m)

2405 k = 3

2406 fe.calc = GPAW(mode=PW(800),

2407 h=0.15,

2408 occupations=FermiDirac(width=0.03),

2409 xc=xc,

2410 kpts=(k, k, k),

2411 convergence={'energy': 1e-8},

2412 parallel={'domain': 1, 'augment_grids': True},

2413 txt=self.folder / 'fe_pw_distorted.txt')

2414 fe.set_cell(np.dot(fe.cell,

2415 [[1.02, 0, 0.03],

2416 [0, 0.99, -0.02],

2417 [0.2, -0.01, 1.03]]),

2418 scale_atoms=True)

2419 fe.get_potential_energy()

2420 return fe.calc

2421

2422 @gpwfile

2423 def si8_fd(self):

2424 a = 5.404

2425 bulk = Atoms(symbols='Si8',

2426 scaled_positions=[(0, 0, 0),

2427 (0, 0.5, 0.5),

2428 (0.5, 0, 0.5),

2429 (0.5, 0.5, 0),

2430 (0.25, 0.25, 0.25),

2431 (0.25, 0.75, 0.75),

2432 (0.75, 0.25, 0.75),

2433 (0.75, 0.75, 0.25)],

2434 pbc=True, cell=(a, a, a))

2435 n = 20

2436 calc = GPAW(mode='fd',

2437 gpts=(n, n, n),

2438 nbands=8 * 3,

2439 occupations=FermiDirac(width=0.01),

2440 kpts=(1, 1, 1))

2441 bulk.calc = calc

2442 bulk.get_potential_energy()

2443

2444 return bulk.calc

2445

2446 @gpwfile

2447 def si_pw_distorted(self):

2448 xc = 'TPSS'

2449 si = bulk('Si')

2450 k = 3

2451 si.calc = GPAW(mode=PW(250),

2452 mixer=Mixer(0.7, 5, 50.0),

2453 xc=xc,

2454 occupations=FermiDirac(0.01),

2455 kpts=(k, k, k),

2456 convergence={'energy': 1e-8},

2457 parallel={'domain': min(2, world.size)},

2458 txt=self.folder / 'si_pw_distorted.txt')

2459 si.set_cell(np.dot(si.cell,

2460 [[1.02, 0, 0.03],

2461 [0, 0.99, -0.02],

2462 [0.2, -0.01, 1.03]]),

2463 scale_atoms=True)

2464 si.get_potential_energy()

2465 return si.calc

2466

2467 @gpwfile

2468 def IBiTe_pw_monolayer(self):

2469 # janus material. material parameters obtained from c2db.

2470 from ase.atoms import Atoms

2471 IBiTe_positions = np.array([[0, 2.552, 7.802],

2472 [0, 0, 9.872],

2473 [2.210, 1.276, 11.575]])

2474 IBiTe = Atoms('IBiTe', positions=IBiTe_positions)

2475 IBiTe.pbc = [True, True, False]

2476 cell = np.array([[4.4219, 0, 0.0, ],

2477 [-2.211, 3.829, 0.0],

2478 [0.0, 0.0, 19.5]])

2479 IBiTe.cell = cell

2480 calc = GPAW(mode=PW(200),

2481 xc='LDA',

2482 occupations=FermiDirac(0.01),

2483 kpts={'size': (6, 6, 1), 'gamma': True},

2484 txt=None)

2485 IBiTe.calc = calc

2486 IBiTe.get_potential_energy()

2487 return IBiTe.calc

2488

2489 def _intraband(self, spinpol: bool):

2490 atoms = bulk('Na')

2491 if spinpol:

2492 atoms.set_initial_magnetic_moments([[0.1]])

2493 atoms.calc = GPAW(mode=PW(300),

2494 kpts={'size': (8, 8, 8), 'gamma': True},

2495 parallel={'band': 1},

2496 txt=None)

2497 atoms.get_potential_energy()

2498 atoms.calc.diagonalize_full_hamiltonian(nbands=20)

2499 return atoms.calc

2500

2501 @gpwfile

2502 def intraband_spinpaired_fulldiag(self):

2503 return self._intraband(False)

2504

2505 @gpwfile

2506 def intraband_spinpolarized_fulldiag(self):

2507 return self._intraband(True)

2508

2509

2510# We add Si fixtures with various symmetries to the GPWFiles namespace

2511for name, method in si_gpwfiles().items():

2512 setattr(GPWFiles, name, method)

2513

2514

2515if __name__ == '__main__':

2516 import sys

2517 name = sys.argv[1]

2518 calc = getattr(GPWFiles(Path()), name)()

2519 calc.write(name + '.gpw', mode='all')