Coverage for /builds/hweiske/ase/ase/calculators/cp2k.py: 83.38%

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.

104

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

109

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

125

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

167

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

169 command = None

170

171 default_parameters = dict(

172 auto_write=False,

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

174 basis_set_file='BASIS_MOLOPT',

175 charge=0,

176 cutoff=400 * Rydberg,

177 force_eval_method="Quickstep",

178 inp='',

179 max_scf=50,

180 multiplicity=None,

181 potential_file='POTENTIAL',

182 pseudo_potential='auto',

183 stress_tensor=True,

184 uks=False,

185 poisson_solver='auto',

186 xc='LDA',

187 print_level='LOW')

188

189 def __init__(self, restart=None,

190 ignore_bad_restart_file=Calculator._deprecated,

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

192 debug=False, **kwargs):

193 """Construct CP2K-calculator object."""

194

195 self._debug = debug

196 self._force_env_id = None

197 self._shell = None

198 self.label = None

199 self.parameters = None

200 self.results = None

201 self.atoms = None

202

203 # Several places are check to determine self.command

204 if command is not None:

205 self.command = command

206 elif CP2K.command is not None:

207 self.command = CP2K.command

208 else:

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

210

211 super().__init__(restart=restart,

212 ignore_bad_restart_file=ignore_bad_restart_file,

213 label=label, atoms=atoms, **kwargs)

214 if restart is not None:

215 self.read(restart)

216

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

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

219

220 def __del__(self):

221 """Terminate cp2k_shell child process"""

222 self.close()

223

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

225 self.close()

226

227 def close(self):

228 """Close the attached shell"""

229 if self._shell is not None:

230 self._shell.close()

231 self._shell = None

232 self._force_env_id = None # Force env must be recreated

233

234 def set(self, **kwargs):

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

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

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

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

239 for key in kwargs:

240 if key not in self.default_parameters:

241 raise CalculatorSetupError(msg % key)

242

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

244 if changed_parameters:

245 self.reset()

246

247 def write(self, label):

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

249 if self._debug:

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

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

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

253 from ase.io.jsonio import write_json

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

255 write_json(fd, self.results)

256

257 def read(self, label):

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

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

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

261 from ase.io.jsonio import read_json

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

263 self.results = read_json(fd)

264

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

266 system_changes=all_changes):

267 """Do the calculation."""

268

269 if not properties:

270 properties = ['energy']

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

272

273 # Start the shell if needed

274 if self._shell is None:

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

276

277 if self._debug:

278 print("system_changes:", system_changes)

279

280 if 'numbers' in system_changes:

281 self._release_force_env()

282

283 if self._force_env_id is None:

284 self._create_force_env()

285

286 # enable eV and Angstrom as units

287 self._shell.send('UNITS_EV_A')

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

289

290 n_atoms = len(self.atoms)

291 if 'cell' in system_changes:

292 cell = self.atoms.get_cell()

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

294 for i in range(3):

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

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

297

298 if 'positions' in system_changes:

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

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

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

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

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

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

305 assert max_change >= 0 # sanity check

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

307

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

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

310

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

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

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

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

315

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

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

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

319 assert nvals == 3 * n_atoms # sanity check

320 for i in range(n_atoms):

321 line = self._shell.recv()

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

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

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

325 self.results['forces'] = forces

326

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

328 line = self._shell.recv()

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

330

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

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

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

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

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

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

337

338 if self.parameters.auto_write:

339 self.write(self.label)

340

341 def _create_force_env(self):

342 """Instantiates a new force-environment"""

343 assert self._force_env_id is None

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

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

346 print('Creating directory: ' + label_dir)

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

348

349 inp = self._generate_input()

350 inp_fn = self.label + '.inp'

351 out_fn = self.label + '.out'

352 self._write_file(inp_fn, inp)

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

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

355 assert self._force_env_id > 0

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

357

358 def _write_file(self, fn, content):

359 """Write content to a file"""

360 if self._debug:

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

362 print(content)

363 if self._shell.version < 2.0:

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

365 fd.write(content)

366 else:

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

368 if self._shell.version < 2.1:

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

370 self._shell.send('WRITE_FILE')

371 self._shell.send(fn)

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

373 for line in lines:

374 self._shell.send(line)

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

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

377

378 def _release_force_env(self):

379 """Destroys the current force-environment"""

380 if self._force_env_id:

381 if self._shell.isready:

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

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

384 else:

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

386 warn(msg, RuntimeWarning)

387 self._force_env_id = None

388

389 def _generate_input(self):

390 """Generates a CP2K input file"""

391 p = self.parameters

392 root = parse_input(p.inp)

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

394 if p.print_level:

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

396 if p.force_eval_method:

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

398 if p.stress_tensor:

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

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

401 '_SECTION_PARAMETERS_ ON')

402 if p.basis_set_file:

403 root.add_keyword('FORCE_EVAL/DFT',

404 'BASIS_SET_FILE_NAME ' + p.basis_set_file)

405 if p.potential_file:

406 root.add_keyword('FORCE_EVAL/DFT',

407 'POTENTIAL_FILE_NAME ' + p.potential_file)

408 if p.cutoff:

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

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

411 if p.max_scf:

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

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

414

415 if p.xc:

416 legacy_libxc = ""

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

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

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

420 # libxc input section changed over time

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

422 legacy_libxc += " " + functional # handled later

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

424 s = InputSection(name='LIBXC')

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

426 xc_sec.subsections.append(s)

427 elif functional.startswith("XC_"):

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

429 xc_sec.subsections.append(s)

430 else:

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

432 xc_sec.subsections.append(s)

433 if legacy_libxc:

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

435 'FUNCTIONAL ' + legacy_libxc)

436

437 if p.uks:

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

439

440 if p.multiplicity:

441 root.add_keyword('FORCE_EVAL/DFT',

442 'MULTIPLICITY %d' % p.multiplicity)

443

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

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

446

447 # add Poisson solver if needed

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

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

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

451

452 # write coords

453 syms = self.atoms.get_chemical_symbols()

454 atoms = self.atoms.get_positions()

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

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

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

458

459 # write cell

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

461 if len(pbc) == 0:

462 pbc = 'NONE'

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

464 c = self.atoms.get_cell()

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

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

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

468

469 # determine pseudo-potential

470 potential = p.pseudo_potential

471 if p.pseudo_potential == 'auto':

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

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

474 else:

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

476 warn(msg, RuntimeWarning)

477 potential = 'GTH-PBE' # fall back

478

479 # write atomic kinds

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

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

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

483 if elem not in kinds.keys():

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

485 subsys.append(s)

486 kinds[elem] = s

487 if p.basis_set:

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

489 if potential:

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

491

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

493 return '\n'.join(output_lines)

494

495

496class Cp2kShell:

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

498

499 def __init__(self, command, debug):

500 """Construct CP2K-shell object"""

501

502 self.isready = False

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

504 self._debug = debug

505

506 # launch cp2k_shell child process

507 assert 'cp2k_shell' in command

508 if self._debug:

509 print(command)

510 self._child = subprocess.Popen(

511 command, shell=True, universal_newlines=True,

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

513 self.expect('* READY')

514

515 # check version of shell

516 self.send('VERSION')

517 line = self.recv()

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

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

520 'Probably the shell version is too old. '

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

522

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

524 self.version = float(shell_version)

525 assert self.version >= 1.0

526

527 self.expect('* READY')

528

529 # enable harsh mode, stops on any error

530 self.send('HARSH')

531 self.expect('* READY')

532

533 def __del__(self):

534 """Terminate cp2k_shell child process"""

535 self.close()

536

537 def close(self):

538 """Terminate cp2k_shell child process"""

539 if self.isready:

540 self.send('EXIT')

541 self._child.communicate()

542 rtncode = self._child.wait()

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

544 elif getattr(self, '_child', None) is not None:

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

546 self._child.terminate()

547 self._child.communicate()

548 self._child = None

549 self.version = None

550 self.isready = False

551

552 def send(self, line):

553 """Send a line to the cp2k_shell"""

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

555 if self._debug:

556 print('Sending: ' + line)

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

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

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

560 self.isready = False

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

562

563 def recv(self):

564 """Receive a line from the cp2k_shell"""

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

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

567 if self._debug:

568 print('Received: ' + line)

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

570 return line

571

572 def expect(self, line):

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

574 received = self.recv()

575 assert received == line

576

577

578class InputSection:

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

580

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

582 self.name = name.upper()

583 self.params = params

584 self.keywords = []

585 self.subsections = []

586

587 def write(self):

588 """Outputs input section as string"""

589 output = []

590 for k in self.keywords:

591 output.append(k)

592 for s in self.subsections:

593 if s.params:

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

595 else:

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

597 for l in s.write():

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

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

600 return output

601

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

603 """Adds a keyword to section."""

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

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

606 if len(candidates) == 0:

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

608 self.subsections.append(s)

609 candidates = [s]

610 elif len(candidates) != 1:

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

612

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

614 if len(parts) > 1:

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

616 elif key == '_SECTION_PARAMETERS_':

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

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

619 raise Exception(msg)

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

621 else:

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

623 if unique and key in old_keys:

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

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

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

627

628 def get_subsection(self, path):

629 """Finds a subsection"""

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

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

632 if len(candidates) > 1:

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

634 if len(candidates) == 0:

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

636 self.subsections.append(s)

637 candidates = [s]

638 if len(parts) == 1:

639 return candidates[0]

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

641

642

643def parse_input(inp):

644 """Parses the given CP2K input string"""

645 root_section = InputSection('CP2K_INPUT')

646 section_stack = [root_section]

647

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

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

650 if len(line) == 0:

651 continue

652

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

654 s = section_stack.pop()

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

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

657 name = parts[0][1:]

658 if len(parts) > 1:

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

660 else:

661 s = InputSection(name=name)

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

663 section_stack.append(s)

664 else:

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

666

667 return root_section