Coverage for /builds/debichem-team/python-ase/ase/calculators/cp2k.py : 84.35%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""This module defines an ASE interface to CP2K.
3https://www.cp2k.org/
4Author: Ole Schuett <ole.schuett@mat.ethz.ch>
5"""
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)
19class CP2K(Calculator):
20 """ASE-Calculator for CP2K.
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.
28 CP2K is freely available under the GPL license.
29 It is written in Fortran 2003 and can be run efficiently in parallel.
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.
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.
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:
49 >>> CP2K.command="env OMP_NUM_THREADS=2 mpiexec -np 4 cp2k_shell.psmp"
51 Arguments:
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.
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::
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 >>> '''
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 """
152 implemented_properties = ['energy', 'free_energy', 'forces', 'stress']
153 command = None
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')
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."""
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
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
196 Calculator.__init__(self, restart=restart,
197 ignore_bad_restart_file=ignore_bad_restart_file,
198 label=label, atoms=atoms, **kwargs)
200 self._shell = Cp2kShell(self.command, self._debug)
202 if restart is not None:
203 self.read(restart)
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)
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)
220 changed_parameters = Calculator.set(self, **kwargs)
221 if changed_parameters:
222 self.reset()
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)
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)
242 def calculate(self, atoms=None, properties=None,
243 system_changes=all_changes):
244 """Do the calculation."""
246 if not properties:
247 properties = ['energy']
248 Calculator.calculate(self, atoms, properties, system_changes)
250 if self._debug:
251 print("system_changes:", system_changes)
253 if 'numbers' in system_changes:
254 self._release_force_env()
256 if self._force_env_id is None:
257 self._create_force_env()
259 # enable eV and Angstrom as units
260 self._shell.send('UNITS_EV_A')
261 self._shell.expect('* READY')
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')
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')
281 self._shell.send('EVAL_EF %d' % self._force_env_id)
282 self._shell.expect('* READY')
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')
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
300 self._shell.send('GET_STRESS %d' % self._force_env_id)
301 line = self._shell.recv()
302 self._shell.expect('* READY')
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
311 if self.parameters.auto_write:
312 self.write(self.label)
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
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')
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')
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
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)
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_"):
397 s = InputSection(name='LIBXC')
398 s.keywords.append('FUNCTIONAL ' + functional)
399 xc_sec.subsections.append(s)
400 else:
401 s = InputSection(name=functional.upper())
402 xc_sec.subsections.append(s)
403 if legacy_libxc:
404 root.add_keyword('FORCE_EVAL/DFT/XC/XC_FUNCTIONAL/LIBXC',
405 'FUNCTIONAL ' + legacy_libxc)
407 if p.uks:
408 root.add_keyword('FORCE_EVAL/DFT', 'UNRESTRICTED_KOHN_SHAM ON')
410 if p.charge and p.charge != 0:
411 root.add_keyword('FORCE_EVAL/DFT', 'CHARGE %d' % p.charge)
413 # add Poisson solver if needed
414 if p.poisson_solver == 'auto' and not any(self.atoms.get_pbc()):
415 root.add_keyword('FORCE_EVAL/DFT/POISSON', 'PERIODIC NONE')
416 root.add_keyword('FORCE_EVAL/DFT/POISSON', 'PSOLVER MT')
418 # write coords
419 syms = self.atoms.get_chemical_symbols()
420 atoms = self.atoms.get_positions()
421 for elm, pos in zip(syms, atoms):
422 line = '%s %.18e %.18e %.18e' % (elm, pos[0], pos[1], pos[2])
423 root.add_keyword('FORCE_EVAL/SUBSYS/COORD', line, unique=False)
425 # write cell
426 pbc = ''.join([a for a, b in zip('XYZ', self.atoms.get_pbc()) if b])
427 if len(pbc) == 0:
428 pbc = 'NONE'
429 root.add_keyword('FORCE_EVAL/SUBSYS/CELL', 'PERIODIC ' + pbc)
430 c = self.atoms.get_cell()
431 for i, a in enumerate('ABC'):
432 line = '%s %.18e %.18e %.18e' % (a, c[i, 0], c[i, 1], c[i, 2])
433 root.add_keyword('FORCE_EVAL/SUBSYS/CELL', line)
435 # determine pseudo-potential
436 potential = p.pseudo_potential
437 if p.pseudo_potential == 'auto':
438 if p.xc and p.xc.upper() in ('LDA', 'PADE', 'BP', 'BLYP', 'PBE',):
439 potential = 'GTH-' + p.xc.upper()
440 else:
441 msg = 'No matching pseudo potential found, using GTH-PBE'
442 warn(msg, RuntimeWarning)
443 potential = 'GTH-PBE' # fall back
445 # write atomic kinds
446 subsys = root.get_subsection('FORCE_EVAL/SUBSYS').subsections
447 kinds = dict([(s.params, s) for s in subsys if s.name == "KIND"])
448 for elem in set(self.atoms.get_chemical_symbols()):
449 if elem not in kinds.keys():
450 s = InputSection(name='KIND', params=elem)
451 subsys.append(s)
452 kinds[elem] = s
453 if p.basis_set:
454 kinds[elem].keywords.append('BASIS_SET ' + p.basis_set)
455 if potential:
456 kinds[elem].keywords.append('POTENTIAL ' + potential)
458 output_lines = ['!!! Generated by ASE !!!'] + root.write()
459 return '\n'.join(output_lines)
462class Cp2kShell:
463 """Wrapper for CP2K-shell child-process"""
465 def __init__(self, command, debug):
466 """Construct CP2K-shell object"""
468 self.isready = False
469 self.version = 1.0 # assume oldest possible version until verified
470 self._child = None
471 self._debug = debug
473 # launch cp2k_shell child process
474 assert 'cp2k_shell' in command
475 if self._debug:
476 print(command)
477 self._child = Popen(command, shell=True, universal_newlines=True,
478 stdin=PIPE, stdout=PIPE, bufsize=1)
479 self.expect('* READY')
481 # check version of shell
482 self.send('VERSION')
483 line = self.recv()
484 if not line.startswith('CP2K Shell Version:'):
485 raise RuntimeError('Cannot determine version of CP2K shell. '
486 'Probably the shell version is too old. '
487 'Please update to CP2K 3.0 or newer.')
489 shell_version = line.rsplit(":", 1)[1]
490 self.version = float(shell_version)
491 assert self.version >= 1.0
493 self.expect('* READY')
495 # enable harsh mode, stops on any error
496 self.send('HARSH')
497 self.expect('* READY')
499 def __del__(self):
500 """Terminate cp2k_shell child process"""
501 if self.isready:
502 self.send('EXIT')
503 rtncode = self._child.wait()
504 assert rtncode == 0 # child process exited properly?
505 else:
506 warn("CP2K-shell not ready, sending SIGTERM.", RuntimeWarning)
507 self._child.terminate()
508 self._child = None
509 self.version = None
510 self.isready = False
512 def send(self, line):
513 """Send a line to the cp2k_shell"""
514 assert self._child.poll() is None # child process still alive?
515 if self._debug:
516 print('Sending: ' + line)
517 if self.version < 2.1 and len(line) >= 80:
518 raise Exception('Buffer overflow, upgrade CP2K to r16779 or later')
519 assert(len(line) < 800) # new input buffer size
520 self.isready = False
521 self._child.stdin.write(line + '\n')
523 def recv(self):
524 """Receive a line from the cp2k_shell"""
525 assert self._child.poll() is None # child process still alive?
526 line = self._child.stdout.readline().strip()
527 if self._debug:
528 print('Received: ' + line)
529 self.isready = line == '* READY'
530 return line
532 def expect(self, line):
533 """Receive a line and asserts that it matches the expected one"""
534 received = self.recv()
535 assert received == line
538class InputSection:
539 """Represents a section of a CP2K input file"""
540 def __init__(self, name, params=None):
541 self.name = name.upper()
542 self.params = params
543 self.keywords = []
544 self.subsections = []
546 def write(self):
547 """Outputs input section as string"""
548 output = []
549 for k in self.keywords:
550 output.append(k)
551 for s in self.subsections:
552 if s.params:
553 output.append('&%s %s' % (s.name, s.params))
554 else:
555 output.append('&%s' % s.name)
556 for l in s.write():
557 output.append(' %s' % l)
558 output.append('&END %s' % s.name)
559 return output
561 def add_keyword(self, path, line, unique=True):
562 """Adds a keyword to section."""
563 parts = path.upper().split('/', 1)
564 candidates = [s for s in self.subsections if s.name == parts[0]]
565 if len(candidates) == 0:
566 s = InputSection(name=parts[0])
567 self.subsections.append(s)
568 candidates = [s]
569 elif len(candidates) != 1:
570 raise Exception('Multiple %s sections found ' % parts[0])
572 key = line.split()[0].upper()
573 if len(parts) > 1:
574 candidates[0].add_keyword(parts[1], line, unique)
575 elif key == '_SECTION_PARAMETERS_':
576 if candidates[0].params is not None:
577 msg = 'Section parameter of section %s already set' % parts[0]
578 raise Exception(msg)
579 candidates[0].params = line.split(' ', 1)[1].strip()
580 else:
581 old_keys = [k.split()[0].upper() for k in candidates[0].keywords]
582 if unique and key in old_keys:
583 msg = 'Keyword %s already present in section %s'
584 raise Exception(msg % (key, parts[0]))
585 candidates[0].keywords.append(line)
587 def get_subsection(self, path):
588 """Finds a subsection"""
589 parts = path.upper().split('/', 1)
590 candidates = [s for s in self.subsections if s.name == parts[0]]
591 if len(candidates) > 1:
592 raise Exception('Multiple %s sections found ' % parts[0])
593 if len(candidates) == 0:
594 s = InputSection(name=parts[0])
595 self.subsections.append(s)
596 candidates = [s]
597 if len(parts) == 1:
598 return candidates[0]
599 return candidates[0].get_subsection(parts[1])
602def parse_input(inp):
603 """Parses the given CP2K input string"""
604 root_section = InputSection('CP2K_INPUT')
605 section_stack = [root_section]
607 for line in inp.split('\n'):
608 line = line.split('!', 1)[0].strip()
609 if len(line) == 0:
610 continue
612 if line.upper().startswith('&END'):
613 s = section_stack.pop()
614 elif line[0] == '&':
615 parts = line.split(' ', 1)
616 name = parts[0][1:]
617 if len(parts) > 1:
618 s = InputSection(name=name, params=parts[1].strip())
619 else:
620 s = InputSection(name=name)
621 section_stack[-1].subsections.append(s)
622 section_stack.append(s)
623 else:
624 section_stack[-1].keywords.append(line)
626 return root_section