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