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

2 

3https://www.cp2k.org/ 

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

5""" 

6 

7 

8import os 

9import os.path 

10from warnings import warn 

11from subprocess import Popen, PIPE 

12import numpy as np 

13import ase.io 

14from ase.units import Rydberg 

15from ase.calculators.calculator import (Calculator, all_changes, Parameters, 

16 CalculatorSetupError) 

17 

18 

19class CP2K(Calculator): 

20 """ASE-Calculator for CP2K. 

21 

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

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

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

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

26 pair and many-body potentials. 

27 

28 CP2K is freely available under the GPL license. 

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

30 

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

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

33 by the CP2K-calulator. 

34 

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

36 designed for interactive sessions. When a calculator object is 

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

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

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

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

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

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

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

44 mpiexec-process. 

45 

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

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

48 

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

50 

51 Arguments: 

52 

53 auto_write: bool 

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

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

56 calculation, which mimics the behavior of the 

57 ``FileIOCalculator``. Default is ``False``. 

58 basis_set: str 

59 Name of the basis set to be use. 

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

61 basis_set_file: str 

62 Filename of the basis set file. 

63 Default is ``BASIS_MOLOPT``. 

64 Set the environment variable $CP2K_DATA_DIR 

65 to enabled automatic file discovered. 

66 charge: float 

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

68 command: str 

69 The command used to launch the CP2K-shell. 

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

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

72 and then the environment variable 

73 ``$ASE_CP2K_COMMAND`` are checked. 

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

75 cutoff: float 

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

77 debug: bool 

78 Flag to enable debug mode. This will print all 

79 communication between the CP2K-shell and the 

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

81 force_eval_method: str 

82 The method CP2K uses to evaluate energies and forces. 

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

84 module for electronic structure methods like DFT. 

85 inp: str 

86 CP2K input template. If present, the calculator will 

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

88 it to launch CP2K. Hence, this generic mechanism 

89 gives access to all features of CP2K. 

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

91 of the corresponding input section. 

92 

93 This input template is important for advanced CP2K 

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

95 zone integration. The example below illustrates some common 

96 options:: 

97 

98 >>> inp = '''&FORCE_EVAL 

99 >>> &DFT 

100 >>> &KPOINTS 

101 >>> SCHEME MONKHORST-PACK 12 12 8 

102 >>> &END KPOINTS 

103 >>> &SCF 

104 >>> ADDED_MOS 10 

105 >>> &SMEAR 

106 >>> METHOD FERMI_DIRAC 

107 >>> ELECTRONIC_TEMPERATURE [K] 500.0 

108 >>> &END SMEAR 

109 >>> &END SCF 

110 >>> &END DFT 

111 >>> &END FORCE_EVAL 

112 >>> ''' 

113 

114 max_scf: int 

115 Maximum number of SCF iteration to be performed for 

116 one optimization. Default is ``50``. 

117 poisson_solver: str 

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

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

120 potential_file: str 

121 Filename of the pseudo-potential file. 

122 Default is ``POTENTIAL``. 

123 Set the environment variable $CP2K_DATA_DIR 

124 to enabled automatic file discovered. 

125 pseudo_potential: str 

126 Name of the pseudo-potential to be use. 

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

128 potential from the employed XC-functional, 

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

130 stress_tensor: bool 

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

132 Default is ``True``. 

133 uks: bool 

134 Requests an unrestricted Kohn-Sham calculations. 

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

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

137 xc: str 

138 Name of exchange and correlation functional. 

139 Accepts all functions supported by CP2K itself or libxc. 

140 Default is ``LDA``. 

141 print_level: str 

142 PRINT_LEVEL of global output. 

143 Possible options are: 

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

145 HIGH Lots of output 

146 LOW Little output 

147 MEDIUM Quite some output 

148 SILENT Almost no output 

149 Default is 'LOW' 

150 """ 

151 

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

153 command = None 

154 

155 default_parameters = dict( 

156 auto_write=False, 

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

158 basis_set_file='BASIS_MOLOPT', 

159 charge=0, 

160 cutoff=400 * Rydberg, 

161 force_eval_method="Quickstep", 

162 inp='', 

163 max_scf=50, 

164 potential_file='POTENTIAL', 

165 pseudo_potential='auto', 

166 stress_tensor=True, 

167 uks=False, 

168 poisson_solver='auto', 

169 xc='LDA', 

170 print_level='LOW') 

171 

172 def __init__(self, restart=None, 

173 ignore_bad_restart_file=Calculator._deprecated, 

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

175 debug=False, **kwargs): 

176 """Construct CP2K-calculator object.""" 

177 

178 self._debug = debug 

179 self._force_env_id = None 

180 self._shell = None 

181 self.label = None 

182 self.parameters = None 

183 self.results = None 

184 self.atoms = None 

185 

186 # Several places are check to determine self.command 

187 if command is not None: 

188 self.command = command 

189 elif CP2K.command is not None: 

190 self.command = CP2K.command 

191 elif 'ASE_CP2K_COMMAND' in os.environ: 

192 self.command = os.environ['ASE_CP2K_COMMAND'] 

193 else: 

194 self.command = 'cp2k_shell' # default 

195 

196 Calculator.__init__(self, restart=restart, 

197 ignore_bad_restart_file=ignore_bad_restart_file, 

198 label=label, atoms=atoms, **kwargs) 

199 

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

201 

202 if restart is not None: 

203 self.read(restart) 

204 

205 def __del__(self): 

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

207 if self._shell: 

208 self._release_force_env() 

209 del self._shell 

210 

211 def set(self, **kwargs): 

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

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

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

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

216 for key in kwargs: 

217 if key not in self.default_parameters: 

218 raise CalculatorSetupError(msg % key) 

219 

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

221 if changed_parameters: 

222 self.reset() 

223 

224 def write(self, label): 

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

226 if self._debug: 

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

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

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

230 from ase.io.jsonio import write_json 

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

232 write_json(fd, self.results) 

233 

234 def read(self, label): 

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

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

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

238 from ase.io.jsonio import read_json 

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

240 self.results = read_json(fd) 

241 

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

243 system_changes=all_changes): 

244 """Do the calculation.""" 

245 

246 if not properties: 

247 properties = ['energy'] 

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

249 

250 if self._debug: 

251 print("system_changes:", system_changes) 

252 

253 if 'numbers' in system_changes: 

254 self._release_force_env() 

255 

256 if self._force_env_id is None: 

257 self._create_force_env() 

258 

259 # enable eV and Angstrom as units 

260 self._shell.send('UNITS_EV_A') 

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

262 

263 n_atoms = len(self.atoms) 

264 if 'cell' in system_changes: 

265 cell = self.atoms.get_cell() 

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

267 for i in range(3): 

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

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

270 

271 if 'positions' in system_changes: 

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

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

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

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

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

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

278 assert max_change >= 0 # sanity check 

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

280 

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

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

283 

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

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

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

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

288 

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

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

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

292 assert nvals == 3 * n_atoms # sanity check 

293 for i in range(n_atoms): 

294 line = self._shell.recv() 

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

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

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

298 self.results['forces'] = forces 

299 

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

301 line = self._shell.recv() 

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

303 

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

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

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

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

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

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

310 

311 if self.parameters.auto_write: 

312 self.write(self.label) 

313 

314 def _create_force_env(self): 

315 """Instantiates a new force-environment""" 

316 assert self._force_env_id is None 

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

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

319 print('Creating directory: ' + label_dir) 

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

321 

322 inp = self._generate_input() 

323 inp_fn = self.label + '.inp' 

324 out_fn = self.label + '.out' 

325 self._write_file(inp_fn, inp) 

326 self._shell.send('LOAD %s %s' % (inp_fn, out_fn)) 

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

328 assert self._force_env_id > 0 

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

330 

331 def _write_file(self, fn, content): 

332 """Write content to a file""" 

333 if self._debug: 

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

335 print(content) 

336 if self._shell.version < 2.0: 

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

338 fd.write(content) 

339 else: 

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

341 if self._shell.version < 2.1: 

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

343 self._shell.send('WRITE_FILE') 

344 self._shell.send(fn) 

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

346 for line in lines: 

347 self._shell.send(line) 

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

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

350 

351 def _release_force_env(self): 

352 """Destroys the current force-environment""" 

353 if self._force_env_id: 

354 if self._shell.isready: 

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

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

357 else: 

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

359 warn(msg, RuntimeWarning) 

360 self._force_env_id = None 

361 

362 def _generate_input(self): 

363 """Generates a CP2K input file""" 

364 p = self.parameters 

365 root = parse_input(p.inp) 

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

367 if p.print_level: 

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

369 if p.force_eval_method: 

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

371 if p.stress_tensor: 

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

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

374 '_SECTION_PARAMETERS_ ON') 

375 if p.basis_set_file: 

376 root.add_keyword('FORCE_EVAL/DFT', 

377 'BASIS_SET_FILE_NAME ' + p.basis_set_file) 

378 if p.potential_file: 

379 root.add_keyword('FORCE_EVAL/DFT', 

380 'POTENTIAL_FILE_NAME ' + p.potential_file) 

381 if p.cutoff: 

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

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

384 if p.max_scf: 

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

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

387 

388 if p.xc: 

389 legacy_libxc = "" 

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

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

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

393 # libxc input section changed over time 

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

395 legacy_libxc += " " + functional # handled later 

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

397 s = InputSection(name='LIBXC') 

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

399 xc_sec.subsections.append(s) 

400 elif functional.startswith("XC_"): 

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

402 xc_sec.subsections.append(s) 

403 else: 

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

405 xc_sec.subsections.append(s) 

406 if legacy_libxc: 

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

408 'FUNCTIONAL ' + legacy_libxc) 

409 

410 if p.uks: 

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

412 

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

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

415 

416 # add Poisson solver if needed 

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

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

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

420 

421 # write coords 

422 syms = self.atoms.get_chemical_symbols() 

423 atoms = self.atoms.get_positions() 

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

425 line = '%s %.18e %.18e %.18e' % (elm, pos[0], pos[1], pos[2]) 

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

427 

428 # write cell 

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

430 if len(pbc) == 0: 

431 pbc = 'NONE' 

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

433 c = self.atoms.get_cell() 

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

435 line = '%s %.18e %.18e %.18e' % (a, c[i, 0], c[i, 1], c[i, 2]) 

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

437 

438 # determine pseudo-potential 

439 potential = p.pseudo_potential 

440 if p.pseudo_potential == 'auto': 

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

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

443 else: 

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

445 warn(msg, RuntimeWarning) 

446 potential = 'GTH-PBE' # fall back 

447 

448 # write atomic kinds 

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

450 kinds = dict([(s.params, s) for s in subsys if s.name == "KIND"]) 

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

452 if elem not in kinds.keys(): 

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

454 subsys.append(s) 

455 kinds[elem] = s 

456 if p.basis_set: 

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

458 if potential: 

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

460 

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

462 return '\n'.join(output_lines) 

463 

464 

465class Cp2kShell: 

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

467 

468 def __init__(self, command, debug): 

469 """Construct CP2K-shell object""" 

470 

471 self.isready = False 

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

473 self._debug = debug 

474 

475 # launch cp2k_shell child process 

476 assert 'cp2k_shell' in command 

477 if self._debug: 

478 print(command) 

479 self._child = Popen(command, shell=True, universal_newlines=True, 

480 stdin=PIPE, stdout=PIPE, bufsize=1) 

481 self.expect('* READY') 

482 

483 # check version of shell 

484 self.send('VERSION') 

485 line = self.recv() 

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

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

488 'Probably the shell version is too old. ' 

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

490 

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

492 self.version = float(shell_version) 

493 assert self.version >= 1.0 

494 

495 self.expect('* READY') 

496 

497 # enable harsh mode, stops on any error 

498 self.send('HARSH') 

499 self.expect('* READY') 

500 

501 def __del__(self): 

502 """Terminate cp2k_shell child process""" 

503 if self.isready: 

504 self.send('EXIT') 

505 self._child.communicate() 

506 rtncode = self._child.wait() 

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

508 else: 

509 warn("CP2K-shell not ready, sending SIGTERM.", RuntimeWarning) 

510 self._child.terminate() 

511 self._child.communicate() 

512 self._child = None 

513 self.version = None 

514 self.isready = False 

515 

516 def send(self, line): 

517 """Send a line to the cp2k_shell""" 

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

519 if self._debug: 

520 print('Sending: ' + line) 

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

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

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

524 self.isready = False 

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

526 

527 def recv(self): 

528 """Receive a line from the cp2k_shell""" 

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

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

531 if self._debug: 

532 print('Received: ' + line) 

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

534 return line 

535 

536 def expect(self, line): 

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

538 received = self.recv() 

539 assert received == line 

540 

541 

542class InputSection: 

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

544 

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

546 self.name = name.upper() 

547 self.params = params 

548 self.keywords = [] 

549 self.subsections = [] 

550 

551 def write(self): 

552 """Outputs input section as string""" 

553 output = [] 

554 for k in self.keywords: 

555 output.append(k) 

556 for s in self.subsections: 

557 if s.params: 

558 output.append('&%s %s' % (s.name, s.params)) 

559 else: 

560 output.append('&%s' % s.name) 

561 for l in s.write(): 

562 output.append(' %s' % l) 

563 output.append('&END %s' % s.name) 

564 return output 

565 

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

567 """Adds a keyword to section.""" 

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

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

570 if len(candidates) == 0: 

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

572 self.subsections.append(s) 

573 candidates = [s] 

574 elif len(candidates) != 1: 

575 raise Exception('Multiple %s sections found ' % parts[0]) 

576 

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

578 if len(parts) > 1: 

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

580 elif key == '_SECTION_PARAMETERS_': 

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

582 msg = 'Section parameter of section %s already set' % parts[0] 

583 raise Exception(msg) 

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

585 else: 

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

587 if unique and key in old_keys: 

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

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

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

591 

592 def get_subsection(self, path): 

593 """Finds a subsection""" 

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

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

596 if len(candidates) > 1: 

597 raise Exception('Multiple %s sections found ' % parts[0]) 

598 if len(candidates) == 0: 

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

600 self.subsections.append(s) 

601 candidates = [s] 

602 if len(parts) == 1: 

603 return candidates[0] 

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

605 

606 

607def parse_input(inp): 

608 """Parses the given CP2K input string""" 

609 root_section = InputSection('CP2K_INPUT') 

610 section_stack = [root_section] 

611 

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

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

614 if len(line) == 0: 

615 continue 

616 

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

618 s = section_stack.pop() 

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

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

621 name = parts[0][1:] 

622 if len(parts) > 1: 

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

624 else: 

625 s = InputSection(name=name) 

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

627 section_stack.append(s) 

628 else: 

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

630 

631 return root_section