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

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

2 

3https://www.cp2k.org/ 

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

5""" 

6 

7import os 

8import os.path 

9import subprocess 

10from contextlib import AbstractContextManager 

11from warnings import warn 

12 

13import numpy as np 

14 

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 

24 

25 

26class CP2K(Calculator, AbstractContextManager): 

27 """ASE-Calculator for CP2K. 

28 

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. 

34 

35 CP2K is freely available under the GPL license. 

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

37 

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. 

41 

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. 

52 

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:: 

55 

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

57 

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:: 

61 

62 with CP2K() as calc: 

63 calc.get_potential_energy(atoms) 

64 

65 The shell will be restarted if you call the calculator object again. 

66 

67 Arguments: 

68 

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. 

108 

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:: 

113 

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 ''' 

129 

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 """ 

176 

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

178 command = None 

179 

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 ) 

199 

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.""" 

205 

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 

213 

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') 

221 

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) 

227 

228 # Start the shell by default, which is how SocketIOCalculator 

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

230 

231 def __del__(self): 

232 """Terminate cp2k_shell child process""" 

233 self.close() 

234 

235 def __exit__(self, __exc_type, __exc_value, __traceback): 

236 self.close() 

237 

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 

244 

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) 

253 

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

255 if changed_parameters: 

256 self.reset() 

257 

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) 

267 

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) 

275 

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

277 system_changes=all_changes): 

278 """Do the calculation.""" 

279 

280 if not properties: 

281 properties = ['energy'] 

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

283 

284 # Start the shell if needed 

285 if self._shell is None: 

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

287 

288 if self._debug: 

289 print("system_changes:", system_changes) 

290 

291 if 'numbers' in system_changes: 

292 self._release_force_env() 

293 

294 if self._force_env_id is None: 

295 self._create_force_env() 

296 

297 # enable eV and Angstrom as units 

298 self._shell.send('UNITS_EV_A') 

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

300 

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') 

308 

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') 

334 

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

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

337 

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') 

342 

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 

353 

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

355 line = self._shell.recv() 

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

357 

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 

364 

365 if self.parameters.auto_write: 

366 self.write(self.label) 

367 

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 

375 

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') 

384 

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') 

404 

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 

415 

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) 

441 

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) 

463 

464 if p.uks: 

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

466 

467 if p.multiplicity: 

468 root.add_keyword('FORCE_EVAL/DFT', 

469 'MULTIPLICITY %d' % p.multiplicity) 

470 

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

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

473 

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') 

478 

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) 

485 

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) 

495 

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 

505 

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) 

518 

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

520 return '\n'.join(output_lines) 

521 

522 

523class Cp2kShell: 

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

525 

526 def __init__(self, command, debug): 

527 """Construct CP2K-shell object""" 

528 

529 self.isready = False 

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

531 self._debug = debug 

532 

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') 

541 

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}') 

550 

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

552 self.version = float(shell_version) 

553 assert self.version >= 1.0 

554 

555 self.expect('* READY') 

556 

557 # enable harsh mode, stops on any error 

558 self.send('HARSH') 

559 self.expect('* READY') 

560 

561 def __del__(self): 

562 """Terminate cp2k_shell child process""" 

563 self.close() 

564 

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 

579 

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') 

590 

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 

599 

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 

604 

605 

606class InputSection: 

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

608 

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

610 self.name = name.upper() 

611 self.params = params 

612 self.keywords = [] 

613 self.subsections = [] 

614 

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 

629 

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 ') 

640 

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) 

655 

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]) 

669 

670 

671def parse_input(inp): 

672 """Parses the given CP2K input string""" 

673 root_section = InputSection('CP2K_INPUT') 

674 section_stack = [root_section] 

675 

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

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

678 if len(line) == 0: 

679 continue 

680 

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) 

694 

695 return root_section