Coverage for /builds/debichem-team/python-ase/ase/calculators/cp2k.py: 81.23%
373 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-03-06 04:00 +0000
« prev ^ index » next coverage.py v7.5.3, created at 2025-03-06 04:00 +0000
1"""This module defines an ASE interface to CP2K.
3https://www.cp2k.org/
4Author: Ole Schuett <ole.schuett@mat.ethz.ch>
5"""
7import os
8import os.path
9import subprocess
10from contextlib import AbstractContextManager
11from warnings import warn
13import numpy as np
15import ase.io
16from ase.calculators.calculator import (
17 Calculator,
18 CalculatorSetupError,
19 Parameters,
20 all_changes,
21)
22from ase.config import cfg
23from ase.units import Rydberg
26class CP2K(Calculator, AbstractContextManager):
27 """ASE-Calculator for CP2K.
29 CP2K is a program to perform atomistic and molecular simulations of solid
30 state, liquid, molecular, and biological systems. It provides a general
31 framework for different methods such as e.g., density functional theory
32 (DFT) using a mixed Gaussian and plane waves approach (GPW) and classical
33 pair and many-body potentials.
35 CP2K is freely available under the GPL license.
36 It is written in Fortran 2003 and can be run efficiently in parallel.
38 Check https://www.cp2k.org about how to obtain and install CP2K.
39 Make sure that you also have the CP2K-shell available, since it is required
40 by the CP2K-calulator.
42 The CP2K-calculator relies on the CP2K-shell. The CP2K-shell was originally
43 designed for interactive sessions. When a calculator object is
44 instantiated, it launches a CP2K-shell as a subprocess in the background
45 and communications with it through stdin/stdout pipes. This has the
46 advantage that the CP2K process is kept alive for the whole lifetime of
47 the calculator object, i.e. there is no startup overhead for a sequence
48 of energy evaluations. Furthermore, the usage of pipes avoids slow file-
49 system I/O. This mechanism even works for MPI-parallelized runs, because
50 stdin/stdout of the first rank are forwarded by the MPI-environment to the
51 mpiexec-process.
53 The command used by the calculator to launch the CP2K-shell is
54 ``cp2k_shell``. To run a parallelized simulation use something like this::
56 CP2K.command="env OMP_NUM_THREADS=2 mpiexec -np 4 cp2k_shell.psmp"
58 The CP2K-shell can be shut down by calling :meth:`close`.
59 The close method will be called automatically when using the calculator as
60 part of a with statement::
62 with CP2K() as calc:
63 calc.get_potential_energy(atoms)
65 The shell will be restarted if you call the calculator object again.
67 Arguments:
69 auto_write: bool
70 Flag to enable the auto-write mode. If enabled the
71 ``write()`` routine is called after every
72 calculation, which mimics the behavior of the
73 ``FileIOCalculator``. Default is ``False``.
74 basis_set: str
75 Name of the basis set to be use.
76 The default is ``DZVP-MOLOPT-SR-GTH``.
77 basis_set_file: str
78 Filename of the basis set file.
79 Default is ``BASIS_MOLOPT``.
80 Set the environment variable $CP2K_DATA_DIR
81 to enabled automatic file discovered.
82 charge: float
83 The total charge of the system. Default is ``0``.
84 command: str
85 The command used to launch the CP2K-shell.
86 If ``command`` is not passed as an argument to the
87 constructor, the class-variable ``CP2K.command``,
88 and then the environment variable
89 ``$ASE_CP2K_COMMAND`` are checked.
90 Eventually, ``cp2k_shell`` is used as default.
91 cutoff: float
92 The cutoff of the finest grid level. Default is ``400 * Rydberg``.
93 debug: bool
94 Flag to enable debug mode. This will print all
95 communication between the CP2K-shell and the
96 CP2K-calculator. Default is ``False``.
97 force_eval_method: str
98 The method CP2K uses to evaluate energies and forces.
99 The default is ``Quickstep``, which is CP2K's
100 module for electronic structure methods like DFT.
101 inp: str
102 CP2K input template. If present, the calculator will
103 augment the template, e.g. with coordinates, and use
104 it to launch CP2K. Hence, this generic mechanism
105 gives access to all features of CP2K.
106 Note, that most keywords accept ``None`` to disable the generation
107 of the corresponding input section.
109 This input template is important for advanced CP2K
110 inputs, but is also needed for e.g. controlling the Brillouin
111 zone integration. The example below illustrates some common
112 options::
114 inp = '''&FORCE_EVAL
115 &DFT
116 &KPOINTS
117 SCHEME MONKHORST-PACK 12 12 8
118 &END KPOINTS
119 &SCF
120 ADDED_MOS 10
121 &SMEAR
122 METHOD FERMI_DIRAC
123 ELECTRONIC_TEMPERATURE [K] 500.0
124 &END SMEAR
125 &END SCF
126 &END DFT
127 &END FORCE_EVAL
128 '''
130 max_scf: int
131 Maximum number of SCF iteration to be performed for
132 one optimization. Default is ``50``.
133 multiplicity: int, default=None
134 Select the multiplicity of the system
135 (two times the total spin plus one).
136 If None, multiplicity is not explicitly given in the input file.
137 poisson_solver: str
138 The poisson solver to be used. Currently, the only supported
139 values are ``auto`` and ``None``. Default is ``auto``.
140 potential_file: str
141 Filename of the pseudo-potential file.
142 Default is ``POTENTIAL``.
143 Set the environment variable $CP2K_DATA_DIR
144 to enabled automatic file discovered.
145 pseudo_potential: str
146 Name of the pseudo-potential to be use.
147 Default is ``auto``. This tries to infer the
148 potential from the employed XC-functional,
149 otherwise it falls back to ``GTH-PBE``.
150 stress_tensor: bool
151 Indicates whether the analytic stress-tensor should be calculated.
152 Default is ``True``.
153 uks: bool
154 Requests an unrestricted Kohn-Sham calculations.
155 This is need for spin-polarized systems, ie. with an
156 odd number of electrons. Default is ``False``.
157 xc: str
158 Name of exchange and correlation functional.
159 Accepts all functions supported by CP2K itself or libxc.
160 Default is ``LDA``.
161 print_level: str
162 PRINT_LEVEL of global output.
163 Possible options are:
164 DEBUG Everything is written out, useful for debugging purposes only
165 HIGH Lots of output
166 LOW Little output
167 MEDIUM Quite some output
168 SILENT Almost no output
169 Default is 'LOW'
170 set_pos_file: bool
171 Send updated positions to the CP2K shell via file instead of
172 via stdin, which can bypass limitations for sending large
173 structures via stdin for CP2K built with some MPI libraries.
174 Requires CP2K 2024.2
175 """
177 implemented_properties = ['energy', 'free_energy', 'forces', 'stress']
178 command = None
180 default_parameters = dict(
181 auto_write=False,
182 basis_set='DZVP-MOLOPT-SR-GTH',
183 basis_set_file='BASIS_MOLOPT',
184 charge=0,
185 cutoff=400 * Rydberg,
186 force_eval_method="Quickstep",
187 inp='',
188 max_scf=50,
189 multiplicity=None,
190 potential_file='POTENTIAL',
191 pseudo_potential='auto',
192 stress_tensor=True,
193 uks=False,
194 poisson_solver='auto',
195 xc='LDA',
196 print_level='LOW',
197 set_pos_file=False,
198 )
200 def __init__(self, restart=None,
201 ignore_bad_restart_file=Calculator._deprecated,
202 label='cp2k', atoms=None, command=None,
203 debug=False, **kwargs):
204 """Construct CP2K-calculator object."""
206 self._debug = debug
207 self._force_env_id = None
208 self._shell = None
209 self.label = None
210 self.parameters = None
211 self.results = None
212 self.atoms = None
214 # Several places are check to determine self.command
215 if command is not None:
216 self.command = command
217 elif CP2K.command is not None:
218 self.command = CP2K.command
219 else:
220 self.command = cfg.get('ASE_CP2K_COMMAND', 'cp2k_shell')
222 super().__init__(restart=restart,
223 ignore_bad_restart_file=ignore_bad_restart_file,
224 label=label, atoms=atoms, **kwargs)
225 if restart is not None:
226 self.read(restart)
228 # Start the shell by default, which is how SocketIOCalculator
229 self._shell = Cp2kShell(self.command, self._debug)
231 def __del__(self):
232 """Terminate cp2k_shell child process"""
233 self.close()
235 def __exit__(self, __exc_type, __exc_value, __traceback):
236 self.close()
238 def close(self):
239 """Close the attached shell"""
240 if self._shell is not None:
241 self._shell.close()
242 self._shell = None
243 self._force_env_id = None # Force env must be recreated
245 def set(self, **kwargs):
246 """Set parameters like set(key1=value1, key2=value2, ...)."""
247 msg = '"%s" is not a known keyword for the CP2K calculator. ' \
248 'To access all features of CP2K by means of an input ' \
249 'template, consider using the "inp" keyword instead.'
250 for key in kwargs:
251 if key not in self.default_parameters:
252 raise CalculatorSetupError(msg % key)
254 changed_parameters = Calculator.set(self, **kwargs)
255 if changed_parameters:
256 self.reset()
258 def write(self, label):
259 'Write atoms, parameters and calculated results into restart files.'
260 if self._debug:
261 print("Writing restart to: ", label)
262 self.atoms.write(label + '_restart.traj')
263 self.parameters.write(label + '_params.ase')
264 from ase.io.jsonio import write_json
265 with open(label + '_results.json', 'w') as fd:
266 write_json(fd, self.results)
268 def read(self, label):
269 'Read atoms, parameters and calculated results from restart files.'
270 self.atoms = ase.io.read(label + '_restart.traj')
271 self.parameters = Parameters.read(label + '_params.ase')
272 from ase.io.jsonio import read_json
273 with open(label + '_results.json') as fd:
274 self.results = read_json(fd)
276 def calculate(self, atoms=None, properties=None,
277 system_changes=all_changes):
278 """Do the calculation."""
280 if not properties:
281 properties = ['energy']
282 Calculator.calculate(self, atoms, properties, system_changes)
284 # Start the shell if needed
285 if self._shell is None:
286 self._shell = Cp2kShell(self.command, self._debug)
288 if self._debug:
289 print("system_changes:", system_changes)
291 if 'numbers' in system_changes:
292 self._release_force_env()
294 if self._force_env_id is None:
295 self._create_force_env()
297 # enable eV and Angstrom as units
298 self._shell.send('UNITS_EV_A')
299 self._shell.expect('* READY')
301 n_atoms = len(self.atoms)
302 if 'cell' in system_changes:
303 cell = self.atoms.get_cell()
304 self._shell.send('SET_CELL %d' % self._force_env_id)
305 for i in range(3):
306 self._shell.send('%.18e %.18e %.18e' % tuple(cell[i, :]))
307 self._shell.expect('* READY')
309 if 'positions' in system_changes:
310 if self.parameters.set_pos_file:
311 # TODO: Update version number when released
312 if self._shell.version < 7:
313 raise ValueError('SET_POS_FILE requires > CP2K 2024.2')
314 pos: np.ndarray = self.atoms.get_positions()
315 fn = self.label + '.pos'
316 with open(fn, 'w') as fp:
317 print(3 * n_atoms, file=fp)
318 for pos in self.atoms.get_positions():
319 print('%.18e %.18e %.18e' % tuple(pos), file=fp)
320 self._shell.send(f'SET_POS_FILE {fn} {self._force_env_id}')
321 else:
322 if len(atoms) > 100 and 'psmp' in self.command:
323 warn('ASE may stall when passing large structures'
324 ' to MPI versions of CP2K.'
325 ' Consider using `set_pos_file=True`.')
326 self._shell.send('SET_POS %d' % self._force_env_id)
327 self._shell.send('%d' % (3 * n_atoms))
328 for pos in self.atoms.get_positions():
329 self._shell.send('%.18e %.18e %.18e' % tuple(pos))
330 self._shell.send('*END')
331 max_change = float(self._shell.recv())
332 assert max_change >= 0 # sanity check
333 self._shell.expect('* READY')
335 self._shell.send('EVAL_EF %d' % self._force_env_id)
336 self._shell.expect('* READY')
338 self._shell.send('GET_E %d' % self._force_env_id)
339 self.results['energy'] = float(self._shell.recv())
340 self.results['free_energy'] = self.results['energy']
341 self._shell.expect('* READY')
343 forces = np.zeros(shape=(n_atoms, 3))
344 self._shell.send('GET_F %d' % self._force_env_id)
345 nvals = int(self._shell.recv())
346 assert nvals == 3 * n_atoms # sanity check
347 for i in range(n_atoms):
348 line = self._shell.recv()
349 forces[i, :] = [float(x) for x in line.split()]
350 self._shell.expect('* END')
351 self._shell.expect('* READY')
352 self.results['forces'] = forces
354 self._shell.send('GET_STRESS %d' % self._force_env_id)
355 line = self._shell.recv()
356 self._shell.expect('* READY')
358 stress = np.array([float(x) for x in line.split()]).reshape(3, 3)
359 assert np.all(stress == np.transpose(stress)) # should be symmetric
360 # Convert 3x3 stress tensor to Voigt form as required by ASE
361 stress = np.array([stress[0, 0], stress[1, 1], stress[2, 2],
362 stress[1, 2], stress[0, 2], stress[0, 1]])
363 self.results['stress'] = -1.0 * stress # cp2k uses the opposite sign
365 if self.parameters.auto_write:
366 self.write(self.label)
368 def _create_force_env(self):
369 """Instantiates a new force-environment"""
370 assert self._force_env_id is None
371 label_dir = os.path.dirname(self.label)
372 if len(label_dir) > 0 and not os.path.exists(label_dir):
373 print('Creating directory: ' + label_dir)
374 os.makedirs(label_dir) # cp2k expects dirs to exist
376 inp = self._generate_input()
377 inp_fn = self.label + '.inp'
378 out_fn = self.label + '.out'
379 self._write_file(inp_fn, inp)
380 self._shell.send(f'LOAD {inp_fn} {out_fn}')
381 self._force_env_id = int(self._shell.recv())
382 assert self._force_env_id > 0
383 self._shell.expect('* READY')
385 def _write_file(self, fn, content):
386 """Write content to a file"""
387 if self._debug:
388 print('Writting to file: ' + fn)
389 print(content)
390 if self._shell.version < 2.0:
391 with open(fn, 'w') as fd:
392 fd.write(content)
393 else:
394 lines = content.split('\n')
395 if self._shell.version < 2.1:
396 lines = [l.strip() for l in lines] # save chars
397 self._shell.send('WRITE_FILE')
398 self._shell.send(fn)
399 self._shell.send('%d' % len(lines))
400 for line in lines:
401 self._shell.send(line)
402 self._shell.send('*END')
403 self._shell.expect('* READY')
405 def _release_force_env(self):
406 """Destroys the current force-environment"""
407 if self._force_env_id:
408 if self._shell.isready:
409 self._shell.send('DESTROY %d' % self._force_env_id)
410 self._shell.expect('* READY')
411 else:
412 msg = "CP2K-shell not ready, could not release force_env."
413 warn(msg, RuntimeWarning)
414 self._force_env_id = None
416 def _generate_input(self):
417 """Generates a CP2K input file"""
418 p = self.parameters
419 root = parse_input(p.inp)
420 root.add_keyword('GLOBAL', 'PROJECT ' + self.label)
421 if p.print_level:
422 root.add_keyword('GLOBAL', 'PRINT_LEVEL ' + p.print_level)
423 if p.force_eval_method:
424 root.add_keyword('FORCE_EVAL', 'METHOD ' + p.force_eval_method)
425 if p.stress_tensor:
426 root.add_keyword('FORCE_EVAL', 'STRESS_TENSOR ANALYTICAL')
427 root.add_keyword('FORCE_EVAL/PRINT/STRESS_TENSOR',
428 '_SECTION_PARAMETERS_ ON')
429 if p.basis_set_file:
430 root.add_keyword('FORCE_EVAL/DFT',
431 'BASIS_SET_FILE_NAME ' + p.basis_set_file)
432 if p.potential_file:
433 root.add_keyword('FORCE_EVAL/DFT',
434 'POTENTIAL_FILE_NAME ' + p.potential_file)
435 if p.cutoff:
436 root.add_keyword('FORCE_EVAL/DFT/MGRID',
437 'CUTOFF [eV] %.18e' % p.cutoff)
438 if p.max_scf:
439 root.add_keyword('FORCE_EVAL/DFT/SCF', 'MAX_SCF %d' % p.max_scf)
440 root.add_keyword('FORCE_EVAL/DFT/LS_SCF', 'MAX_SCF %d' % p.max_scf)
442 if p.xc:
443 legacy_libxc = ""
444 for functional in p.xc.split():
445 functional = functional.replace("LDA", "PADE") # resolve alias
446 xc_sec = root.get_subsection('FORCE_EVAL/DFT/XC/XC_FUNCTIONAL')
447 # libxc input section changed over time
448 if functional.startswith("XC_") and self._shell.version < 3.0:
449 legacy_libxc += " " + functional # handled later
450 elif functional.startswith("XC_") and self._shell.version < 5.0:
451 s = InputSection(name='LIBXC')
452 s.keywords.append('FUNCTIONAL ' + functional)
453 xc_sec.subsections.append(s)
454 elif functional.startswith("XC_"):
455 s = InputSection(name=functional[3:])
456 xc_sec.subsections.append(s)
457 else:
458 s = InputSection(name=functional.upper())
459 xc_sec.subsections.append(s)
460 if legacy_libxc:
461 root.add_keyword('FORCE_EVAL/DFT/XC/XC_FUNCTIONAL/LIBXC',
462 'FUNCTIONAL ' + legacy_libxc)
464 if p.uks:
465 root.add_keyword('FORCE_EVAL/DFT', 'UNRESTRICTED_KOHN_SHAM ON')
467 if p.multiplicity:
468 root.add_keyword('FORCE_EVAL/DFT',
469 'MULTIPLICITY %d' % p.multiplicity)
471 if p.charge and p.charge != 0:
472 root.add_keyword('FORCE_EVAL/DFT', 'CHARGE %d' % p.charge)
474 # add Poisson solver if needed
475 if p.poisson_solver == 'auto' and not any(self.atoms.get_pbc()):
476 root.add_keyword('FORCE_EVAL/DFT/POISSON', 'PERIODIC NONE')
477 root.add_keyword('FORCE_EVAL/DFT/POISSON', 'PSOLVER MT')
479 # write coords
480 syms = self.atoms.get_chemical_symbols()
481 atoms = self.atoms.get_positions()
482 for elm, pos in zip(syms, atoms):
483 line = f'{elm} {pos[0]:.18e} {pos[1]:.18e} {pos[2]:.18e}'
484 root.add_keyword('FORCE_EVAL/SUBSYS/COORD', line, unique=False)
486 # write cell
487 pbc = ''.join([a for a, b in zip('XYZ', self.atoms.get_pbc()) if b])
488 if len(pbc) == 0:
489 pbc = 'NONE'
490 root.add_keyword('FORCE_EVAL/SUBSYS/CELL', 'PERIODIC ' + pbc)
491 c = self.atoms.get_cell()
492 for i, a in enumerate('ABC'):
493 line = f'{a} {c[i, 0]:.18e} {c[i, 1]:.18e} {c[i, 2]:.18e}'
494 root.add_keyword('FORCE_EVAL/SUBSYS/CELL', line)
496 # determine pseudo-potential
497 potential = p.pseudo_potential
498 if p.pseudo_potential == 'auto':
499 if p.xc and p.xc.upper() in ('LDA', 'PADE', 'BP', 'BLYP', 'PBE',):
500 potential = 'GTH-' + p.xc.upper()
501 else:
502 msg = 'No matching pseudo potential found, using GTH-PBE'
503 warn(msg, RuntimeWarning)
504 potential = 'GTH-PBE' # fall back
506 # write atomic kinds
507 subsys = root.get_subsection('FORCE_EVAL/SUBSYS').subsections
508 kinds = {s.params: s for s in subsys if s.name == "KIND"}
509 for elem in set(self.atoms.get_chemical_symbols()):
510 if elem not in kinds.keys():
511 s = InputSection(name='KIND', params=elem)
512 subsys.append(s)
513 kinds[elem] = s
514 if p.basis_set:
515 kinds[elem].keywords.append('BASIS_SET ' + p.basis_set)
516 if potential:
517 kinds[elem].keywords.append('POTENTIAL ' + potential)
519 output_lines = ['!!! Generated by ASE !!!'] + root.write()
520 return '\n'.join(output_lines)
523class Cp2kShell:
524 """Wrapper for CP2K-shell child-process"""
526 def __init__(self, command, debug):
527 """Construct CP2K-shell object"""
529 self.isready = False
530 self.version = 1.0 # assume oldest possible version until verified
531 self._debug = debug
533 # launch cp2k_shell child process
534 assert 'cp2k_shell' in command
535 if self._debug:
536 print(command)
537 self._child = subprocess.Popen(
538 command, shell=True, universal_newlines=True,
539 stdin=subprocess.PIPE, stdout=subprocess.PIPE, bufsize=1)
540 self.expect('* READY')
542 # check version of shell
543 self.send('VERSION')
544 line = self.recv()
545 if not line.startswith('CP2K Shell Version:'):
546 raise RuntimeError('Cannot determine version of CP2K shell. '
547 'Probably the shell version is too old. '
548 'Please update to CP2K 3.0 or newer. '
549 f'Received: {line}')
551 shell_version = line.rsplit(":", 1)[1]
552 self.version = float(shell_version)
553 assert self.version >= 1.0
555 self.expect('* READY')
557 # enable harsh mode, stops on any error
558 self.send('HARSH')
559 self.expect('* READY')
561 def __del__(self):
562 """Terminate cp2k_shell child process"""
563 self.close()
565 def close(self):
566 """Terminate cp2k_shell child process"""
567 if self.isready:
568 self.send('EXIT')
569 self._child.communicate()
570 rtncode = self._child.wait()
571 assert rtncode == 0 # child process exited properly?
572 elif getattr(self, '_child', None) is not None:
573 warn('CP2K-shell not ready, sending SIGTERM.', RuntimeWarning)
574 self._child.terminate()
575 self._child.communicate()
576 self._child = None
577 self.version = None
578 self.isready = False
580 def send(self, line):
581 """Send a line to the cp2k_shell"""
582 assert self._child.poll() is None # child process still alive?
583 if self._debug:
584 print('Sending: ' + line)
585 if self.version < 2.1 and len(line) >= 80:
586 raise Exception('Buffer overflow, upgrade CP2K to r16779 or later')
587 assert len(line) < 800 # new input buffer size
588 self.isready = False
589 self._child.stdin.write(line + '\n')
591 def recv(self):
592 """Receive a line from the cp2k_shell"""
593 assert self._child.poll() is None # child process still alive?
594 line = self._child.stdout.readline().strip()
595 if self._debug:
596 print('Received: ' + line)
597 self.isready = line == '* READY'
598 return line
600 def expect(self, line):
601 """Receive a line and asserts that it matches the expected one"""
602 received = self.recv()
603 assert received == line
606class InputSection:
607 """Represents a section of a CP2K input file"""
609 def __init__(self, name, params=None):
610 self.name = name.upper()
611 self.params = params
612 self.keywords = []
613 self.subsections = []
615 def write(self):
616 """Outputs input section as string"""
617 output = []
618 for k in self.keywords:
619 output.append(k)
620 for s in self.subsections:
621 if s.params:
622 output.append(f'&{s.name} {s.params}')
623 else:
624 output.append(f'&{s.name}')
625 for l in s.write():
626 output.append(f' {l}')
627 output.append(f'&END {s.name}')
628 return output
630 def add_keyword(self, path, line, unique=True):
631 """Adds a keyword to section."""
632 parts = path.upper().split('/', 1)
633 candidates = [s for s in self.subsections if s.name == parts[0]]
634 if len(candidates) == 0:
635 s = InputSection(name=parts[0])
636 self.subsections.append(s)
637 candidates = [s]
638 elif len(candidates) != 1:
639 raise Exception(f'Multiple {parts[0]} sections found ')
641 key = line.split()[0].upper()
642 if len(parts) > 1:
643 candidates[0].add_keyword(parts[1], line, unique)
644 elif key == '_SECTION_PARAMETERS_':
645 if candidates[0].params is not None:
646 msg = f'Section parameter of section {parts[0]} already set'
647 raise Exception(msg)
648 candidates[0].params = line.split(' ', 1)[1].strip()
649 else:
650 old_keys = [k.split()[0].upper() for k in candidates[0].keywords]
651 if unique and key in old_keys:
652 msg = 'Keyword %s already present in section %s'
653 raise Exception(msg % (key, parts[0]))
654 candidates[0].keywords.append(line)
656 def get_subsection(self, path):
657 """Finds a subsection"""
658 parts = path.upper().split('/', 1)
659 candidates = [s for s in self.subsections if s.name == parts[0]]
660 if len(candidates) > 1:
661 raise Exception(f'Multiple {parts[0]} sections found ')
662 if len(candidates) == 0:
663 s = InputSection(name=parts[0])
664 self.subsections.append(s)
665 candidates = [s]
666 if len(parts) == 1:
667 return candidates[0]
668 return candidates[0].get_subsection(parts[1])
671def parse_input(inp):
672 """Parses the given CP2K input string"""
673 root_section = InputSection('CP2K_INPUT')
674 section_stack = [root_section]
676 for line in inp.split('\n'):
677 line = line.split('!', 1)[0].strip()
678 if len(line) == 0:
679 continue
681 if line.upper().startswith('&END'):
682 s = section_stack.pop()
683 elif line[0] == '&':
684 parts = line.split(' ', 1)
685 name = parts[0][1:]
686 if len(parts) > 1:
687 s = InputSection(name=name, params=parts[1].strip())
688 else:
689 s = InputSection(name=name)
690 section_stack[-1].subsections.append(s)
691 section_stack.append(s)
692 else:
693 section_stack[-1].keywords.append(line)
695 return root_section