Coverage for /builds/kinetik161/ase/ase/calculators/cp2k.py: 84.62%

1"""This module defines an ASE interface to CP2K.

3https://www.cp2k.org/

4Author: Ole Schuett <ole.schuett@mat.ethz.ch>

5"""

8import os

9import os.path

10import subprocess

11from warnings import warn

13import numpy as np

15import ase.io

16from ase.calculators.calculator import (Calculator, CalculatorSetupError,

17 Parameters, all_changes)

18from ase.config import cfg

19from ase.units import Rydberg

22class CP2K(Calculator):

23 """ASE-Calculator for CP2K.

25 CP2K is a program to perform atomistic and molecular simulations of solid

26 state, liquid, molecular, and biological systems. It provides a general

27 framework for different methods such as e.g., density functional theory

28 (DFT) using a mixed Gaussian and plane waves approach (GPW) and classical

29 pair and many-body potentials.

31 CP2K is freely available under the GPL license.

32 It is written in Fortran 2003 and can be run efficiently in parallel.

34 Check https://www.cp2k.org about how to obtain and install CP2K.

35 Make sure that you also have the CP2K-shell available, since it is required

36 by the CP2K-calulator.

38 The CP2K-calculator relies on the CP2K-shell. The CP2K-shell was originally

39 designed for interactive sessions. When a calculator object is

40 instantiated, it launches a CP2K-shell as a subprocess in the background

41 and communications with it through stdin/stdout pipes. This has the

42 advantage that the CP2K process is kept alive for the whole lifetime of

43 the calculator object, i.e. there is no startup overhead for a sequence

44 of energy evaluations. Furthermore, the usage of pipes avoids slow file-

45 system I/O. This mechanism even works for MPI-parallelized runs, because

46 stdin/stdout of the first rank are forwarded by the MPI-environment to the

47 mpiexec-process.

49 The command used by the calculator to launch the CP2K-shell is

50 ``cp2k_shell``. To run a parallelized simulation use something like this::

52 CP2K.command="env OMP_NUM_THREADS=2 mpiexec -np 4 cp2k_shell.psmp"

54 Arguments:

56 auto_write: bool

57 Flag to enable the auto-write mode. If enabled the

58 ``write()`` routine is called after every

59 calculation, which mimics the behavior of the

60 ``FileIOCalculator``. Default is ``False``.

61 basis_set: str

62 Name of the basis set to be use.

63 The default is ``DZVP-MOLOPT-SR-GTH``.

64 basis_set_file: str

65 Filename of the basis set file.

66 Default is ``BASIS_MOLOPT``.

67 Set the environment variable $CP2K_DATA_DIR

68 to enabled automatic file discovered.

69 charge: float

70 The total charge of the system. Default is ``0``.

71 command: str

72 The command used to launch the CP2K-shell.

73 If ``command`` is not passed as an argument to the

74 constructor, the class-variable ``CP2K.command``,

75 and then the environment variable

76 ``$ASE_CP2K_COMMAND`` are checked.

77 Eventually, ``cp2k_shell`` is used as default.

78 cutoff: float

79 The cutoff of the finest grid level. Default is ``400 * Rydberg``.

80 debug: bool

81 Flag to enable debug mode. This will print all

82 communication between the CP2K-shell and the

83 CP2K-calculator. Default is ``False``.

84 force_eval_method: str

85 The method CP2K uses to evaluate energies and forces.

86 The default is ``Quickstep``, which is CP2K's

87 module for electronic structure methods like DFT.

88 inp: str

89 CP2K input template. If present, the calculator will

90 augment the template, e.g. with coordinates, and use

91 it to launch CP2K. Hence, this generic mechanism

92 gives access to all features of CP2K.

93 Note, that most keywords accept ``None`` to disable the generation

94 of the corresponding input section.

96 This input template is important for advanced CP2K

97 inputs, but is also needed for e.g. controlling the Brillouin

98 zone integration. The example below illustrates some common

99 options::

100

101 inp = '''&FORCE_EVAL

102 &DFT

103 &KPOINTS

104 SCHEME MONKHORST-PACK 12 12 8

105 &END KPOINTS

106 &SCF

107 ADDED_MOS 10

108 &SMEAR

109 METHOD FERMI_DIRAC

110 ELECTRONIC_TEMPERATURE [K] 500.0

111 &END SMEAR

112 &END SCF

113 &END DFT

114 &END FORCE_EVAL

115 '''

116

117 max_scf: int

118 Maximum number of SCF iteration to be performed for

119 one optimization. Default is ``50``.

120 multiplicity: int, default=None

121 Select the multiplicity of the system

122 (two times the total spin plus one).

123 If None, multiplicity is not explicitly given in the input file.

124 poisson_solver: str

125 The poisson solver to be used. Currently, the only supported

126 values are ``auto`` and ``None``. Default is ``auto``.

127 potential_file: str

128 Filename of the pseudo-potential file.

129 Default is ``POTENTIAL``.

130 Set the environment variable $CP2K_DATA_DIR

131 to enabled automatic file discovered.

132 pseudo_potential: str

133 Name of the pseudo-potential to be use.

134 Default is ``auto``. This tries to infer the

135 potential from the employed XC-functional,

136 otherwise it falls back to ``GTH-PBE``.

137 stress_tensor: bool

138 Indicates whether the analytic stress-tensor should be calculated.

139 Default is ``True``.

140 uks: bool

141 Requests an unrestricted Kohn-Sham calculations.

142 This is need for spin-polarized systems, ie. with an

143 odd number of electrons. Default is ``False``.

144 xc: str

145 Name of exchange and correlation functional.

146 Accepts all functions supported by CP2K itself or libxc.

147 Default is ``LDA``.

148 print_level: str

149 PRINT_LEVEL of global output.

150 Possible options are:

151 DEBUG Everything is written out, useful for debugging purposes only

152 HIGH Lots of output

153 LOW Little output

154 MEDIUM Quite some output

155 SILENT Almost no output

156 Default is 'LOW'

157 """

158

159 implemented_properties = ['energy', 'free_energy', 'forces', 'stress']

160 command = None

161

162 default_parameters = dict(

163 auto_write=False,

164 basis_set='DZVP-MOLOPT-SR-GTH',

165 basis_set_file='BASIS_MOLOPT',

166 charge=0,

167 cutoff=400 * Rydberg,

168 force_eval_method="Quickstep",

169 inp='',

170 max_scf=50,

171 multiplicity=None,

172 potential_file='POTENTIAL',

173 pseudo_potential='auto',

174 stress_tensor=True,

175 uks=False,

176 poisson_solver='auto',

177 xc='LDA',

178 print_level='LOW')

179

180 def __init__(self, restart=None,

181 ignore_bad_restart_file=Calculator._deprecated,

182 label='cp2k', atoms=None, command=None,

183 debug=False, **kwargs):

184 """Construct CP2K-calculator object."""

185

186 self._debug = debug

187 self._force_env_id = None

188 self._shell = None

189 self.label = None

190 self.parameters = None

191 self.results = None

192 self.atoms = None

193

194 # Several places are check to determine self.command

195 if command is not None:

196 self.command = command

197 elif CP2K.command is not None:

198 self.command = CP2K.command

199 else:

200 self.command = cfg.get('ASE_CP2K_COMMAND', 'cp2k_shell')

201

202 super().__init__(restart=restart,

203 ignore_bad_restart_file=ignore_bad_restart_file,

204 label=label, atoms=atoms, **kwargs)

205

206 self._shell = Cp2kShell(self.command, self._debug)

207

208 if restart is not None:

209 self.read(restart)

210

211 def __del__(self):

212 """Release force_env and terminate cp2k_shell child process"""

213 if self._shell:

214 self._release_force_env()

215 del self._shell

216

217 def set(self, **kwargs):

218 """Set parameters like set(key1=value1, key2=value2, ...)."""

219 msg = '"%s" is not a known keyword for the CP2K calculator. ' \

220 'To access all features of CP2K by means of an input ' \

221 'template, consider using the "inp" keyword instead.'

222 for key in kwargs:

223 if key not in self.default_parameters:

224 raise CalculatorSetupError(msg % key)

225

226 changed_parameters = Calculator.set(self, **kwargs)

227 if changed_parameters:

228 self.reset()

229

230 def write(self, label):

231 'Write atoms, parameters and calculated results into restart files.'

232 if self._debug:

233 print("Writing restart to: ", label)

234 self.atoms.write(label + '_restart.traj')

235 self.parameters.write(label + '_params.ase')

236 from ase.io.jsonio import write_json

237 with open(label + '_results.json', 'w') as fd:

238 write_json(fd, self.results)

239

240 def read(self, label):

241 'Read atoms, parameters and calculated results from restart files.'

242 self.atoms = ase.io.read(label + '_restart.traj')

243 self.parameters = Parameters.read(label + '_params.ase')

244 from ase.io.jsonio import read_json

245 with open(label + '_results.json') as fd:

246 self.results = read_json(fd)

247

248 def calculate(self, atoms=None, properties=None,

249 system_changes=all_changes):

250 """Do the calculation."""

251

252 if not properties:

253 properties = ['energy']

254 Calculator.calculate(self, atoms, properties, system_changes)

255

256 if self._debug:

257 print("system_changes:", system_changes)

258

259 if 'numbers' in system_changes:

260 self._release_force_env()

261

262 if self._force_env_id is None:

263 self._create_force_env()

264

265 # enable eV and Angstrom as units

266 self._shell.send('UNITS_EV_A')

267 self._shell.expect('* READY')

268

269 n_atoms = len(self.atoms)

270 if 'cell' in system_changes:

271 cell = self.atoms.get_cell()

272 self._shell.send('SET_CELL %d' % self._force_env_id)

273 for i in range(3):

274 self._shell.send('%.18e %.18e %.18e' % tuple(cell[i, :]))

275 self._shell.expect('* READY')

276

277 if 'positions' in system_changes:

278 self._shell.send('SET_POS %d' % self._force_env_id)

279 self._shell.send('%d' % (3 * n_atoms))

280 for pos in self.atoms.get_positions():

281 self._shell.send('%.18e %.18e %.18e' % tuple(pos))

282 self._shell.send('*END')

283 max_change = float(self._shell.recv())

284 assert max_change >= 0 # sanity check

285 self._shell.expect('* READY')

286

287 self._shell.send('EVAL_EF %d' % self._force_env_id)

288 self._shell.expect('* READY')

289

290 self._shell.send('GET_E %d' % self._force_env_id)

291 self.results['energy'] = float(self._shell.recv())

292 self.results['free_energy'] = self.results['energy']

293 self._shell.expect('* READY')

294

295 forces = np.zeros(shape=(n_atoms, 3))

296 self._shell.send('GET_F %d' % self._force_env_id)

297 nvals = int(self._shell.recv())

298 assert nvals == 3 * n_atoms # sanity check

299 for i in range(n_atoms):

300 line = self._shell.recv()

301 forces[i, :] = [float(x) for x in line.split()]

302 self._shell.expect('* END')

303 self._shell.expect('* READY')

304 self.results['forces'] = forces

305

306 self._shell.send('GET_STRESS %d' % self._force_env_id)

307 line = self._shell.recv()

308 self._shell.expect('* READY')

309

310 stress = np.array([float(x) for x in line.split()]).reshape(3, 3)

311 assert np.all(stress == np.transpose(stress)) # should be symmetric

312 # Convert 3x3 stress tensor to Voigt form as required by ASE

313 stress = np.array([stress[0, 0], stress[1, 1], stress[2, 2],

314 stress[1, 2], stress[0, 2], stress[0, 1]])

315 self.results['stress'] = -1.0 * stress # cp2k uses the opposite sign

316

317 if self.parameters.auto_write:

318 self.write(self.label)

319

320 def _create_force_env(self):

321 """Instantiates a new force-environment"""

322 assert self._force_env_id is None

323 label_dir = os.path.dirname(self.label)

324 if len(label_dir) > 0 and not os.path.exists(label_dir):

325 print('Creating directory: ' + label_dir)

326 os.makedirs(label_dir) # cp2k expects dirs to exist

327

328 inp = self._generate_input()

329 inp_fn = self.label + '.inp'

330 out_fn = self.label + '.out'

331 self._write_file(inp_fn, inp)

332 self._shell.send(f'LOAD {inp_fn} {out_fn}')

333 self._force_env_id = int(self._shell.recv())

334 assert self._force_env_id > 0

335 self._shell.expect('* READY')

336

337 def _write_file(self, fn, content):

338 """Write content to a file"""

339 if self._debug:

340 print('Writting to file: ' + fn)

341 print(content)

342 if self._shell.version < 2.0:

343 with open(fn, 'w') as fd:

344 fd.write(content)

345 else:

346 lines = content.split('\n')

347 if self._shell.version < 2.1:

348 lines = [l.strip() for l in lines] # save chars

349 self._shell.send('WRITE_FILE')

350 self._shell.send(fn)

351 self._shell.send('%d' % len(lines))

352 for line in lines:

353 self._shell.send(line)

354 self._shell.send('*END')

355 self._shell.expect('* READY')

356

357 def _release_force_env(self):

358 """Destroys the current force-environment"""

359 if self._force_env_id:

360 if self._shell.isready:

361 self._shell.send('DESTROY %d' % self._force_env_id)

362 self._shell.expect('* READY')

363 else:

364 msg = "CP2K-shell not ready, could not release force_env."

365 warn(msg, RuntimeWarning)

366 self._force_env_id = None

367

368 def _generate_input(self):

369 """Generates a CP2K input file"""

370 p = self.parameters

371 root = parse_input(p.inp)

372 root.add_keyword('GLOBAL', 'PROJECT ' + self.label)

373 if p.print_level:

374 root.add_keyword('GLOBAL', 'PRINT_LEVEL ' + p.print_level)

375 if p.force_eval_method:

376 root.add_keyword('FORCE_EVAL', 'METHOD ' + p.force_eval_method)

377 if p.stress_tensor:

378 root.add_keyword('FORCE_EVAL', 'STRESS_TENSOR ANALYTICAL')

379 root.add_keyword('FORCE_EVAL/PRINT/STRESS_TENSOR',

380 '_SECTION_PARAMETERS_ ON')

381 if p.basis_set_file:

382 root.add_keyword('FORCE_EVAL/DFT',

383 'BASIS_SET_FILE_NAME ' + p.basis_set_file)

384 if p.potential_file:

385 root.add_keyword('FORCE_EVAL/DFT',

386 'POTENTIAL_FILE_NAME ' + p.potential_file)

387 if p.cutoff:

388 root.add_keyword('FORCE_EVAL/DFT/MGRID',

389 'CUTOFF [eV] %.18e' % p.cutoff)

390 if p.max_scf:

391 root.add_keyword('FORCE_EVAL/DFT/SCF', 'MAX_SCF %d' % p.max_scf)

392 root.add_keyword('FORCE_EVAL/DFT/LS_SCF', 'MAX_SCF %d' % p.max_scf)

393

394 if p.xc:

395 legacy_libxc = ""

396 for functional in p.xc.split():

397 functional = functional.replace("LDA", "PADE") # resolve alias

398 xc_sec = root.get_subsection('FORCE_EVAL/DFT/XC/XC_FUNCTIONAL')

399 # libxc input section changed over time

400 if functional.startswith("XC_") and self._shell.version < 3.0:

401 legacy_libxc += " " + functional # handled later

402 elif functional.startswith("XC_") and self._shell.version < 5.0:

403 s = InputSection(name='LIBXC')

404 s.keywords.append('FUNCTIONAL ' + functional)

405 xc_sec.subsections.append(s)

406 elif functional.startswith("XC_"):

407 s = InputSection(name=functional[3:])

408 xc_sec.subsections.append(s)

409 else:

410 s = InputSection(name=functional.upper())

411 xc_sec.subsections.append(s)

412 if legacy_libxc:

413 root.add_keyword('FORCE_EVAL/DFT/XC/XC_FUNCTIONAL/LIBXC',

414 'FUNCTIONAL ' + legacy_libxc)

415

416 if p.uks:

417 root.add_keyword('FORCE_EVAL/DFT', 'UNRESTRICTED_KOHN_SHAM ON')

418

419 if p.multiplicity:

420 root.add_keyword('FORCE_EVAL/DFT',

421 'MULTIPLICITY %d' % p.multiplicity)

422

423 if p.charge and p.charge != 0:

424 root.add_keyword('FORCE_EVAL/DFT', 'CHARGE %d' % p.charge)

425

426 # add Poisson solver if needed

427 if p.poisson_solver == 'auto' and not any(self.atoms.get_pbc()):

428 root.add_keyword('FORCE_EVAL/DFT/POISSON', 'PERIODIC NONE')

429 root.add_keyword('FORCE_EVAL/DFT/POISSON', 'PSOLVER MT')

430

431 # write coords

432 syms = self.atoms.get_chemical_symbols()

433 atoms = self.atoms.get_positions()

434 for elm, pos in zip(syms, atoms):

435 line = f'{elm} {pos[0]:.18e} {pos[1]:.18e} {pos[2]:.18e}'

436 root.add_keyword('FORCE_EVAL/SUBSYS/COORD', line, unique=False)

437

438 # write cell

439 pbc = ''.join([a for a, b in zip('XYZ', self.atoms.get_pbc()) if b])

440 if len(pbc) == 0:

441 pbc = 'NONE'

442 root.add_keyword('FORCE_EVAL/SUBSYS/CELL', 'PERIODIC ' + pbc)

443 c = self.atoms.get_cell()

444 for i, a in enumerate('ABC'):

445 line = f'{a} {c[i, 0]:.18e} {c[i, 1]:.18e} {c[i, 2]:.18e}'

446 root.add_keyword('FORCE_EVAL/SUBSYS/CELL', line)

447

448 # determine pseudo-potential

449 potential = p.pseudo_potential

450 if p.pseudo_potential == 'auto':

451 if p.xc and p.xc.upper() in ('LDA', 'PADE', 'BP', 'BLYP', 'PBE',):

452 potential = 'GTH-' + p.xc.upper()

453 else:

454 msg = 'No matching pseudo potential found, using GTH-PBE'

455 warn(msg, RuntimeWarning)

456 potential = 'GTH-PBE' # fall back

457

458 # write atomic kinds

459 subsys = root.get_subsection('FORCE_EVAL/SUBSYS').subsections

460 kinds = {s.params: s for s in subsys if s.name == "KIND"}

461 for elem in set(self.atoms.get_chemical_symbols()):

462 if elem not in kinds.keys():

463 s = InputSection(name='KIND', params=elem)

464 subsys.append(s)

465 kinds[elem] = s

466 if p.basis_set:

467 kinds[elem].keywords.append('BASIS_SET ' + p.basis_set)

468 if potential:

469 kinds[elem].keywords.append('POTENTIAL ' + potential)

470

471 output_lines = ['!!! Generated by ASE !!!'] + root.write()

472 return '\n'.join(output_lines)

473

474

475class Cp2kShell:

476 """Wrapper for CP2K-shell child-process"""

477

478 def __init__(self, command, debug):

479 """Construct CP2K-shell object"""

480

481 self.isready = False

482 self.version = 1.0 # assume oldest possible version until verified

483 self._debug = debug

484

485 # launch cp2k_shell child process

486 assert 'cp2k_shell' in command

487 if self._debug:

488 print(command)

489 self._child = subprocess.Popen(

490 command, shell=True, universal_newlines=True,

491 stdin=subprocess.PIPE, stdout=subprocess.PIPE, bufsize=1)

492 self.expect('* READY')

493

494 # check version of shell

495 self.send('VERSION')

496 line = self.recv()

497 if not line.startswith('CP2K Shell Version:'):

498 raise RuntimeError('Cannot determine version of CP2K shell. '

499 'Probably the shell version is too old. '

500 'Please update to CP2K 3.0 or newer.')

501

502 shell_version = line.rsplit(":", 1)[1]

503 self.version = float(shell_version)

504 assert self.version >= 1.0

505

506 self.expect('* READY')

507

508 # enable harsh mode, stops on any error

509 self.send('HARSH')

510 self.expect('* READY')

511

512 def __del__(self):

513 """Terminate cp2k_shell child process"""

514 if self.isready:

515 self.send('EXIT')

516 self._child.communicate()

517 rtncode = self._child.wait()

518 assert rtncode == 0 # child process exited properly?

519 else:

520 if hasattr(self, '_child'):

521 warn('CP2K-shell not ready, sending SIGTERM.', RuntimeWarning)

522 self._child.terminate()

523 self._child.communicate()

524 self._child = None

525 self.version = None

526 self.isready = False

527

528 def send(self, line):

529 """Send a line to the cp2k_shell"""

530 assert self._child.poll() is None # child process still alive?

531 if self._debug:

532 print('Sending: ' + line)

533 if self.version < 2.1 and len(line) >= 80:

534 raise Exception('Buffer overflow, upgrade CP2K to r16779 or later')

535 assert len(line) < 800 # new input buffer size

536 self.isready = False

537 self._child.stdin.write(line + '\n')

538

539 def recv(self):

540 """Receive a line from the cp2k_shell"""

541 assert self._child.poll() is None # child process still alive?

542 line = self._child.stdout.readline().strip()

543 if self._debug:

544 print('Received: ' + line)

545 self.isready = line == '* READY'

546 return line

547

548 def expect(self, line):

549 """Receive a line and asserts that it matches the expected one"""

550 received = self.recv()

551 assert received == line

552

553

554class InputSection:

555 """Represents a section of a CP2K input file"""

556

557 def __init__(self, name, params=None):

558 self.name = name.upper()

559 self.params = params

560 self.keywords = []

561 self.subsections = []

562

563 def write(self):

564 """Outputs input section as string"""

565 output = []

566 for k in self.keywords:

567 output.append(k)

568 for s in self.subsections:

569 if s.params:

570 output.append(f'&{s.name} {s.params}')

571 else:

572 output.append(f'&{s.name}')

573 for l in s.write():

574 output.append(f' {l}')

575 output.append(f'&END {s.name}')

576 return output

577

578 def add_keyword(self, path, line, unique=True):

579 """Adds a keyword to section."""

580 parts = path.upper().split('/', 1)

581 candidates = [s for s in self.subsections if s.name == parts[0]]

582 if len(candidates) == 0:

583 s = InputSection(name=parts[0])

584 self.subsections.append(s)

585 candidates = [s]

586 elif len(candidates) != 1:

587 raise Exception(f'Multiple {parts[0]} sections found ')

588

589 key = line.split()[0].upper()

590 if len(parts) > 1:

591 candidates[0].add_keyword(parts[1], line, unique)

592 elif key == '_SECTION_PARAMETERS_':

593 if candidates[0].params is not None:

594 msg = f'Section parameter of section {parts[0]} already set'

595 raise Exception(msg)

596 candidates[0].params = line.split(' ', 1)[1].strip()

597 else:

598 old_keys = [k.split()[0].upper() for k in candidates[0].keywords]

599 if unique and key in old_keys:

600 msg = 'Keyword %s already present in section %s'

601 raise Exception(msg % (key, parts[0]))

602 candidates[0].keywords.append(line)

603

604 def get_subsection(self, path):

605 """Finds a subsection"""

606 parts = path.upper().split('/', 1)

607 candidates = [s for s in self.subsections if s.name == parts[0]]

608 if len(candidates) > 1:

609 raise Exception(f'Multiple {parts[0]} sections found ')

610 if len(candidates) == 0:

611 s = InputSection(name=parts[0])

612 self.subsections.append(s)

613 candidates = [s]

614 if len(parts) == 1:

615 return candidates[0]

616 return candidates[0].get_subsection(parts[1])

617

618

619def parse_input(inp):

620 """Parses the given CP2K input string"""

621 root_section = InputSection('CP2K_INPUT')

622 section_stack = [root_section]

623

624 for line in inp.split('\n'):

625 line = line.split('!', 1)[0].strip()

626 if len(line) == 0:

627 continue

628

629 if line.upper().startswith('&END'):

630 s = section_stack.pop()

631 elif line[0] == '&':

632 parts = line.split(' ', 1)

633 name = parts[0][1:]

634 if len(parts) > 1:

635 s = InputSection(name=name, params=parts[1].strip())

636 else:

637 s = InputSection(name=name)

638 section_stack[-1].subsections.append(s)

639 section_stack.append(s)

640 else:

641 section_stack[-1].keywords.append(line)

642

643 return root_section