# Copyright (c) 2026, PyRETIS Development Team.
# Distributed under the LGPLv2.1+ License. See LICENSE for more info.
"""A CP2K external MD integrator interface (relaunch / restart-chain).
This module defines the relaunch CP2K engine, selected by
``class = "cp2k_steps"``. It advances CP2K one order-parameter point at a
time, relaunching the executable for each step and chaining the run
through ``previous.restart`` / ``previous.wfn`` files (the slower,
start-and-stop implementation; the canonical continuous engine, which
runs CP2K as a single persistent process and reads its trajectory
on-the-fly, lives in :py:mod:`pyretis.engines.cp2k`).
Important classes defined here
------------------------------
CP2KEngineSteps (:py:class:`.CP2KEngineSteps`)
A class responsible for interfacing CP2K via relaunching.
"""
import logging
import os
import re
import shlex
import numpy as np
from pyretis.engines.external import ExternalMDEngine
from pyretis.core.forcefieldhelp import rescale_system_velocities
from pyretis.core.systemhelp import update_system_box
from pyretis.core.random_gen import create_random_generator
from pyretis.inout.settings import look_for_input_files
from pyretis.core.velocity import is_aimless_velocity_setting
from pyretis.inout.formats.xyz import (
read_xyz_file,
write_xyz_trajectory,
convert_snapshot
)
from pyretis.inout.formats.cp2k import (
update_cp2k_input,
read_cp2k_restart,
read_cp2k_box,
read_cp2k_energy,
read_cp2k_input,
set_parents,
)
logger = logging.getLogger(__name__) # pylint: disable=invalid-name
logger.addHandler(logging.NullHandler())
[docs]
def guess_particle_mass(particle_no, particle_type):
"""Guess a particle mass from its type and convert to CP2K units.
Parameters
----------
particle_no : integer
The particle number. Only used for the logger message.
particle_type : string
The particle type as a string; an element of the periodic table.
Returns
-------
out : float
The guessed particle mass in the engine's internal units.
"""
# Lazy import to avoid the _parts -> setup circular import at load.
from pyretis.engines._parts import PERIODIC_TABLE
logger.info(
'Mass not specified for particle no. %i\n'
'Will guess from particle type "%s"',
particle_no, particle_type,
)
mass = PERIODIC_TABLE.get(particle_type, None)
if mass is None:
msg = (
f"Could not infer mass for particle number {particle_no} with "
f"element '{particle_type}' as it is not in our periodic table."
)
raise ValueError(msg)
particle_mass = 1822.8884858012982 * mass
logger.info(
'-> Using a mass of %f g/mol (%f in internal units)',
mass, particle_mass,
)
return particle_mass
OUTPUT_FILES = {
'energy': '{}-1.ener',
'restart': '{}-1.restart',
'pos': '{}-pos-1.xyz',
'vel': '{}-vel-1.xyz',
'wfn': '{}-RESTART.wfn',
'wfn-bak': '{}-RESTART.wfn.bak-'
}
REGEXP_BACKUP = re.compile(r'\.bak-\d$')
[docs]
def write_for_step_vel(infile, outfile, timestep, subcycles, posfile, vel,
name='md_step', print_freq=None):
"""Create input file for a single step.
Note, the single step actually consists of a number of subcycles.
But from PyRETIS' point of view, this is a single step.
Further, we here assume that we start from a given xyz file and
we also explicitly give the velocities here.
Parameters
----------
infile : string
The input template to use.
outfile : string
The file to create.
timestep : float
The time-step to use for the simulation.
subcycles : integer
The number of sub-cycles to perform.
posfile : string
The (base)name for the input file to read positions from.
vel : numpy.array
The velocities to set in the input.
name : string, optional
A name for the CP2K project.
print_freq : integer, optional
How often we should print to the trajectory file.
"""
if print_freq is None:
print_freq = subcycles
to_update = {
'GLOBAL': {
'data': [f'PROJECT {name}',
'RUN_TYPE MD',
'PRINT_LEVEL LOW'],
'replace': True,
},
'MOTION->MD': {
'data': {'STEPS': subcycles,
'TIMESTEP': timestep}
},
'MOTION->PRINT->RESTART': {
'data': ['BACKUP_COPIES 0'],
'replace': True,
},
'MOTION->PRINT->RESTART->EACH': {
'data': {'MD': print_freq}
},
'MOTION->PRINT->VELOCITIES->EACH': {
'data': {'MD': print_freq}
},
'MOTION->PRINT->TRAJECTORY->EACH': {
'data': {'MD': print_freq}
},
'FORCE_EVAL->SUBSYS->TOPOLOGY': {
'data': {'COORD_FILE_NAME': posfile,
'COORD_FILE_FORMAT': 'xyz'}
},
'FORCE_EVAL->SUBSYS->VELOCITY': {
'data': [],
'replace': True,
},
'FORCE_EVAL->DFT->SCF->PRINT->RESTART': {
'data': ['BACKUP_COPIES 0'],
'replace': True,
},
}
for veli in vel:
to_update['FORCE_EVAL->SUBSYS->VELOCITY']['data'].append(
f'{veli[0]} {veli[1]} {veli[2]}'
)
remove = [
'EXT_RESTART',
'FORCE_EVAL->SUBSYS->COORD'
]
update_cp2k_input(infile, outfile, update=to_update, remove=remove)
[docs]
def write_for_integrate(infile, outfile, timestep, subcycles, posfile,
name='md_step', print_freq=None):
"""Create input file for a single step for the integrate method.
Here, we do minimal changes and just set the time step and subcycles
and the starting configuration.
Parameters
----------
infile : string
The input template to use.
outfile : string
The file to create.
timestep : float
The time-step to use for the simulation.
subcycles : integer
The number of sub-cycles to perform.
posfile : string
The (base)name for the input file to read positions from.
name : string, optional
A name for the CP2K project.
print_freq : integer, optional
How often we should print to the trajectory file.
"""
if print_freq is None:
print_freq = subcycles
to_update = {
'GLOBAL': {
'data': [f'PROJECT {name}',
'RUN_TYPE MD',
'PRINT_LEVEL LOW'],
'replace': True,
},
'MOTION->MD': {
'data': {'STEPS': subcycles,
'TIMESTEP': timestep}
},
'MOTION->PRINT->RESTART': {
'data': ['BACKUP_COPIES 0'],
'replace': True,
},
'MOTION->PRINT->RESTART->EACH': {
'data': {'MD': print_freq}
},
'MOTION->PRINT->VELOCITIES->EACH': {
'data': {'MD': print_freq}
},
'MOTION->PRINT->TRAJECTORY->EACH': {
'data': {'MD': print_freq}
},
'FORCE_EVAL->SUBSYS->TOPOLOGY': {
'data': {'COORD_FILE_NAME': posfile,
'COORD_FILE_FORMAT': 'xyz'}
},
'FORCE_EVAL->DFT->SCF->PRINT->RESTART': {
'data': ['BACKUP_COPIES 0'],
'replace': True,
},
}
remove = [
'EXT_RESTART',
'FORCE_EVAL->SUBSYS->COORD'
]
update_cp2k_input(infile, outfile, update=to_update, remove=remove)
[docs]
def write_for_continue(infile, outfile, timestep, subcycles,
name='md_continue'):
"""
Create input file for a single step.
Note, the single step actually consists of a number of subcycles.
But from PyRETIS' point of view, this is a single step.
Here, we make use of restart files named ``previous.restart``
and ``previous.wfn`` to continue a run.
Parameters
----------
infile : string
The input template to use.
outfile : string
The file to create.
timestep : float
The time-step to use for the simulation.
subcycles : integer
The number of sub-cycles to perform.
name : string, optional
A name for the CP2K project.
"""
to_update = {
'GLOBAL': {
'data': [f'PROJECT {name}',
'RUN_TYPE MD',
'PRINT_LEVEL LOW'],
'replace': True,
},
'MOTION->MD': {
'data': {'STEPS': subcycles,
'TIMESTEP': timestep}
},
'MOTION->PRINT->RESTART': {
'data': ['BACKUP_COPIES 0'],
'replace': True,
},
'MOTION->PRINT->RESTART->EACH': {
'data': {'MD': subcycles}
},
'MOTION->PRINT->VELOCITIES->EACH': {
'data': {'MD': subcycles}
},
'MOTION->PRINT->TRAJECTORY->EACH': {
'data': {'MD': subcycles}
},
'EXT_RESTART': {
'data': ['RESTART_VEL',
'RESTART_POS',
'RESTART_FILE_NAME previous.restart'],
'replace': True
},
'FORCE_EVAL->DFT': {
'data': {'WFN_RESTART_FILE_NAME': 'previous.wfn'},
},
'FORCE_EVAL->DFT->SCF->PRINT->RESTART': {
'data': ['BACKUP_COPIES 0'],
'replace': True,
},
}
remove = [
'FORCE_EVAL->SUBSYS->TOPOLOGY',
'FORCE_EVAL->SUBSYS->VELOCITY',
'FORCE_EVAL->SUBSYS->COORD'
'FORCE_EVAL->DFT->RESTART_FILE_NAME',
]
update_cp2k_input(infile, outfile, update=to_update, remove=remove)
[docs]
def write_for_genvel(infile, outfile, posfile, seed,
name='genvel'): # pragma: no cover
"""Create input file for velocity generation.
2022.05.25: Note: This function is no longer in use
as we now generate velocities internally. However we keep
this function in case of future need.
Parameters
----------
infile : string
The input template to use.
outfile : string
The file to create.
posfile : string
The (base)name for the input file to read positions from.
seed : integer
A seed for generating velocities.
name : string, optional
A name for the CP2K project.
"""
to_update = {
'GLOBAL': {
'data': [f'PROJECT {name}',
f'SEED {seed}',
'RUN_TYPE MD',
'PRINT_LEVEL LOW'],
'replace': True,
},
'FORCE_EVAL->DFT->SCF': {
'data': {'SCF_GUESS': 'ATOMIC'}
},
'MOTION->MD': {
'data': {'STEPS': 1,
'TIMESTEP': 0}
},
'MOTION->PRINT->RESTART': {
'data': ['BACKUP_COPIES 0'],
'replace': True,
},
'MOTION->PRINT->RESTART->EACH': {
'data': {'MD': 1}
},
'MOTION->PRINT->VELOCITIES->EACH': {
'data': {'MD': 1}
},
'MOTION->PRINT->TRAJECTORY->EACH': {
'data': {'MD': 1}
},
'FORCE_EVAL->SUBSYS->TOPOLOGY': {
'data': {'COORD_FILE_NAME': posfile,
'COORD_FILE_FORMAT': 'xyz'}
},
'FORCE_EVAL->DFT->SCF->PRINT->RESTART': {
'data': ['BACKUP_COPIES 0'],
'replace': True,
},
}
remove = [
'EXT_RESTART',
'FORCE_EVAL->SUBSYS->VELOCITY',
'FORCE_EVAL->DFT->RESTART_FILE_NAME',
]
update_cp2k_input(infile, outfile, update=to_update, remove=remove)
[docs]
class CP2KEngineSteps(ExternalMDEngine):
"""
A class for interfacing CP2K.
This class defines the interface to CP2K.
Attributes
----------
cp2k : string
The command for executing CP2K.
input_path : string
The directory where the input files are stored.
timestep : float
The time step used in the CP2K MD simulation.
subcycles : integer
The number of steps each CP2K run is composed of.
rgen : object like :py:class:`.RandomGenerator`
An object we use to set seeds for velocity generation.
extra_files : list
List of extra files which may be required to run CP2K.
"""
default_units = 'cp2k'
[docs]
@classmethod
def get_default_units(cls, settings=None):
"""Return the default unit system for the CP2K engine.
Parameters
----------
settings : dict, optional
Full simulation settings or engine settings. This argument
is ignored for the CP2K engine.
Returns
-------
out : string
The default unit system for the CP2K engine.
"""
_ = settings
return 'cp2k'
[docs]
def __init__(self, cp2k, input_path, timestep, subcycles,
extra_files=None, exe_path=os.path.abspath('.'), seed=0,
conf='initial.xyz', template='cp2k.inp',
velocity_generation='maxwell', temperature=None):
"""Set up the CP2K engine.
Parameters
----------
cp2k : string
The CP2K executable.
input_path : string
The path to where the input files are stored.
timestep : float
The time step used in the CP2K simulation.
subcycles : integer
The number of steps each CP2K run is composed of.
extra_files : list
List of extra files which may be required to run CP2K.
exe_path: string, optional
The path on which the engine is executed.
seed : integer, optional
A seed for the random number generator.
conf : string, optional
The default configuration file (default is 'initial.xyz').
template : string, optional
The CP2K input template (default is 'cp2k.inp').
velocity_generation : string, optional
How shooting velocities are drawn. ``"maxwell"`` (default)
draws them in Python from a Maxwell-Boltzmann distribution
via the injected random generator, using engine-level masses
and temperature (so the draw does not require the system to
carry masses). This mirrors the ``velocity_generation`` knob
of the GROMACS/LAMMPS engines.
temperature : float, optional
The MD temperature for the Maxwell draw. When ``None`` it is
read from the ``MOTION->MD`` section of the CP2K input
template, so the same input directory drives velocity
generation without an extra ``[engine]`` setting.
"""
super().__init__('CP2K external engine', timestep,
subcycles)
self.rgen = create_random_generator({'seed': seed})
self.ext = os.path.splitext(conf)[1][1:]
self.cp2k = shlex.split(cp2k)
logger.info('Command for execution of CP2K: %s', ' '.join(self.cp2k))
# Store input path:
self.input_path = os.path.join(exe_path, input_path)
# Set the defaults input files:
default_files = {
'conf': conf,
'template': template,
}
# Check the presence of the defaults input files or, if absent,
# try to find then by extension.
self.input_files = look_for_input_files(self.input_path, default_files)
# todo, these info can be processed by look_for_input_files using
# the extra_files option.
self.extra_files = []
if extra_files is not None:
for key in extra_files:
fname = os.path.join(self.input_path, key)
if not os.path.isfile(fname):
logger.critical('Extra CP2K input file "%s" not found!',
fname)
else:
self.extra_files.append(fname)
# Velocity-generation policy (engine parity with GROMACS/LAMMPS).
if velocity_generation not in ('maxwell',):
msg = ('cp2k velocity_generation only supports "maxwell" '
'(engine-native regeneration is not implemented); '
'got "%s"')
logger.error(msg, velocity_generation)
raise ValueError(msg % velocity_generation)
self.velocity_generation = velocity_generation
# Engine-level masses + inverse thermal energy for the Maxwell
# draw. Computed lazily on first velocity modification (so simply
# constructing the engine does not require the input files to be
# readable on disk yet); see ``_ensure_velocity_setup``.
self.kb = 3.16681534e-6 # hartree
self._temperature = temperature
self.mass = None
self.temperature = None
self._beta = None
[docs]
def _ensure_velocity_setup(self):
"""Compute engine-level masses + beta on first use.
The masses are guessed from the initial-config atom types and the
temperature is taken from the ``[engine]`` setting or, when that
is absent, the ``MOTION->MD`` section of the CP2K input template.
Deferred from ``__init__`` so the engine can be constructed before
the input files exist on disk.
"""
if self.mass is not None:
return
_, _, _, atoms = self._read_configuration(self.input_files['conf'])
mass = []
if atoms is not None:
mass = [guess_particle_mass(i, name)
for i, name in enumerate(atoms)]
self.mass = np.reshape(mass, (len(mass), 1))
cp2k_temp = None
nodes = read_cp2k_input(self.input_files['template'])
md_settings = set_parents(nodes).get('MOTION->MD')
if md_settings is not None:
for data in md_settings.data:
if 'temperature' in data.lower():
cp2k_temp = float(data.split()[1])
temperature = (
self._temperature if self._temperature is not None else cp2k_temp
)
self.temperature = temperature
self._beta = (
None if temperature is None else 1.0 / (temperature * self.kb)
)
[docs]
def run_cp2k(self, input_file, proj_name):
"""
Run the CP2K executable.
Returns
-------
out : dict
The files created by the run.
"""
cmd = self.cp2k + ['-i', input_file]
logger.debug('Executing CP2K %s: %s', proj_name, input_file)
self.execute_command(cmd, cwd=self.exe_dir, inputs=None)
out = {}
for key, name in OUTPUT_FILES.items():
out[key] = os.path.join(self.exe_dir, name.format(proj_name))
return out
[docs]
def _propagate_from(self, name, path, system, ens_set, msg_file,
reverse=False):
"""
Propagate with CP2K from the current system configuration.
Here, we assume that this method is called after the propagate()
has been called in the parent. The parent is then responsible
for reversing the velocities and also for setting the initial
state of the system.
Parameters
----------
name : string
A name to use for the trajectory we are generating.
path : object like :py:class:`.PathBase`
This is the path we use to fill in phase-space points.
system : object like :py:class:`.System`
The system to act on.
ens_set : dict
It contains the simulations info:
* `interfaces` : list of floats
These defines the interfaces for which we will check the
crossing(s).
msg_file : object like :py:class:`.FileIO`
An object we use for writing out messages that are useful
for inspecting the status of the current propagation.
reverse : boolean, optional
If True, the system will be propagated backward in time.
Returns
-------
success : boolean
This is True if we generated an acceptable path.
status : string
A text description of the current status of the
propagation.
"""
status = f'propagating with CP2K (reverse = {reverse})'
interfaces = ens_set['interfaces']
logger.debug(status)
success = False
left, _, right = interfaces
logger.debug('Adding input files for CP2K')
# First, copy the required input files:
self.add_input_files(self.exe_dir)
# Get positions and velocities from the input file.
initial_conf = system.config[0]
box, xyz, vel, atoms = self._read_configuration(initial_conf)
if box is None:
box, _ = read_cp2k_box(self.input_files['template'])
# Add CP2K input for a single step:
step_input = os.path.join(self.exe_dir, 'step.inp')
write_for_step_vel(self.input_files['template'], step_input,
self.timestep, self.subcycles,
os.path.basename(initial_conf),
vel, name=name)
# And create the input file for continuing:
continue_input = os.path.join(self.exe_dir, 'continue.inp')
write_for_continue(self.input_files['template'], continue_input,
self.timestep, self.subcycles, name=name)
# Get the order parameter before the run:
order = self.calculate_order(system, xyz=xyz, vel=vel, box=box)
traj_file = os.path.join(self.exe_dir, f'{name}.{self.ext}')
# Create a message file with some info about this run:
msg_file.write(
f'# Initial order parameter: {" ".join([str(i) for i in order])}'
)
msg_file.write(f'# Trajectory file is: {traj_file}')
# Run the first step:
msg_file.write('# Running first CP2k step.')
out_files = self.run_cp2k('step.inp', name)
restart_file = os.path.join(self.exe_dir, out_files['restart'])
prestart_file = os.path.join(self.exe_dir, 'previous.restart')
wave_file = os.path.join(self.exe_dir, out_files['wfn'])
pwave_file = os.path.join(self.exe_dir, 'previous.wfn')
# Note: Order is calculated at the END of each iteration!
i = 0
# Write the config so we have a non-empty file:
write_xyz_trajectory(traj_file, xyz, vel, atoms, box, step=i,
append=False)
msg_file.write('# Running main CP2k propagation loop.')
msg_file.write('# Step order parameter cv1 cv2 ...')
for i in range(path.maxlen):
msg_file.write(f'{i} {" ".join([str(j) for j in order])}')
snapshot = {'order': order, 'config': (traj_file, i),
'vel_rev': reverse}
phase_point = self.snapshot_to_system(system, snapshot)
status, success, stop, add = self.add_to_path(path, phase_point,
left, right)
if add and i > 0:
# Write the previous configuration:
write_xyz_trajectory(traj_file, xyz, vel, atoms, box,
step=i)
if stop:
logger.debug('CP2K propagation ended at %i. Reason: %s',
i, status)
break
if i == 0:
pass
elif i > 0:
self._movefile(restart_file, prestart_file)
self._movefile(wave_file, pwave_file)
if i < path.maxlen - 1:
out_files = self.run_cp2k('continue.inp', name)
self._remove_files(self.exe_dir,
self._find_backup_files(self.exe_dir))
# Read config after the step
if i < path.maxlen - 1:
atoms, xyz, vel, box, _ = read_cp2k_restart(restart_file)
order = self.calculate_order(system,
xyz=xyz, vel=vel, box=box)
msg_file.write('# Propagation done.')
energy_file = out_files['energy']
msg_file.write(f'# Reading energies from: {energy_file}')
energy = read_cp2k_energy(energy_file)
end = (i + 1) * self.subcycles
ekin = energy.get('ekin', [])
vpot = energy.get('vpot', [])
etot = energy.get('etot', [])
temp = energy.get('temp', [])
path.update_energies(ekin[:end:self.subcycles],
vpot[:end:self.subcycles],
etot[:end:self.subcycles],
temp[:end:self.subcycles])
for _, files in out_files.items():
self._removefile(files)
self._removefile(prestart_file)
self._removefile(pwave_file)
self._removefile(continue_input)
self._removefile(step_input)
return success, status
[docs]
def integrate(self, system, order_function, steps, thermo='full'):
"""
Propagate several integration steps.
This method will perform several integration steps using
CP2K. It will also calculate the order parameter(s) and
energy terms if requested.
Parameters
----------
system : object like :py:class:`.System`
The system object gives the initial state for the
integration.
order_function : object like :py:class:`.OrderParameter` or None
An order function can be specified if we want to calculate
the order parameter along with the simulation.
steps : integer
The number of steps we are going to perform. Note that we
do not integrate on the first step (e.g. step 0) but we do
obtain the other properties. This is to output the starting
configuration.
thermo : string, optional
Select the thermodynamic properties we are to calculate.
Yields
------
results : dict
The results from a MD step. This contains the state of the system
and order parameter(s) and energies (if calculated).
"""
logger.debug('Integrating with CP2K')
# First, copy the required input files:
logger.debug('Adding input files for CP2K')
self.add_input_files(self.exe_dir)
self.order_function = order_function
# Get positions and velocities from the input file.
initial_conf = self.dump_frame(system)
box, xyz, vel, atoms = self._read_configuration(initial_conf)
if box is None:
box, _ = read_cp2k_box(self.input_files['template'])
name = 'pyretis-cp2k'
logger.debug('Thermo was set to: %s', thermo)
# Add CP2K input for a single step:
step_input = os.path.join(self.exe_dir, 'step.inp')
write_for_integrate(self.input_files['template'], step_input,
self.timestep, self.subcycles,
os.path.basename(initial_conf),
name=name)
# And create the input file for continuing:
continue_input = os.path.join(self.exe_dir, 'continue.inp')
write_for_continue(self.input_files['template'], continue_input,
self.timestep, self.subcycles, name=name)
# Get the order parameter before the run:
if order_function:
order = self.calculate_order(system,
xyz=xyz, vel=vel, box=box)
else:
order = None
traj_file = os.path.join(self.exe_dir, f'{name}.{self.ext}')
# Run the first step:
out_files = self.run_cp2k('step.inp', name)
restart_file = os.path.join(self.exe_dir, out_files['restart'])
prestart_file = os.path.join(self.exe_dir, 'previous.restart')
wave_file = os.path.join(self.exe_dir, out_files['wfn'])
pwave_file = os.path.join(self.exe_dir, 'previous.wfn')
# Note: Order is calculated at the END of each iteration!
i = 0
# Write the config so we have a non-empty file:
write_xyz_trajectory(traj_file, xyz, vel, atoms, box, step=i,
append=False)
for i in range(steps):
if i > 0:
# Write the previous configuration:
write_xyz_trajectory(traj_file, xyz, vel, atoms, box,
step=i)
if i == 0:
pass
elif i > 0:
self._movefile(restart_file, prestart_file)
self._movefile(wave_file, pwave_file)
if i < steps - 1: # Do not integrate the last step.
out_files = self.run_cp2k('continue.inp', name)
self._remove_files(self.exe_dir,
self._find_backup_files(self.exe_dir))
results = {}
if order:
results['order'] = order
# Read config after the step
if i < steps - 1:
atoms, xyz, vel, box, _ = read_cp2k_restart(restart_file)
if order_function:
order = self.calculate_order(
system, xyz=xyz, vel=vel, box=box
)
energy = read_cp2k_energy(out_files['energy'])
end = (i + 1) * self.subcycles
results['thermo'] = {}
for key, val in energy.items():
results['thermo'][key] = val[:end:self.subcycles][-1]
yield results
# Note we do not remove the run-files here as they might be
# useful for the user.
[docs]
def step(self, system, name):
"""Perform a single step with CP2K.
Parameters
----------
system : object like :py:class:`.System`
The system we are integrating.
name : string
To name the output files from the CP2K step.
Returns
-------
out : string
The name of the output configuration, obtained after
completing the step.
"""
initial_conf = self.dump_frame(system)
# Prepare input files etc.:
self.add_input_files(self.exe_dir)
box, xyz, vel, atoms = self._read_configuration(initial_conf)
if box is None:
box, _ = read_cp2k_box(self.input_files['template'])
# Add CP2K input for a single step:
step_input = os.path.join(self.exe_dir, 'step.inp')
write_for_step_vel(self.input_files['template'], step_input,
self.timestep, self.subcycles,
os.path.basename(initial_conf),
vel, name=name)
# Execute single step CP2K...
logger.debug('Executing CP2K')
out_files = self.run_cp2k('step.inp', name)
energy = read_cp2k_energy(out_files['energy'])
# Get the output configuration:
atoms, xyz, vel, box, _ = read_cp2k_restart(out_files['restart'])
conf_out = os.path.join(self.exe_dir, f'{name}.{self.ext}')
write_xyz_trajectory(conf_out, xyz, vel, atoms, box, append=False)
system.config = (conf_out, None)
system.vel_rev = False
system.ekin = energy['ekin'][-1]
system.vpot = energy['vpot'][-1]
update_system_box(system, box)
# Prepare input files etc.:
logger.debug('Removing CP2K output after single step.')
# Remove run-files etc:
for _, files in out_files.items():
self._removefile(files)
self._remove_files(
self.exe_dir,
self._find_backup_files(self.exe_dir)
)
return conf_out
[docs]
@staticmethod
def _find_backup_files(dirname):
"""Return backup-files in the given directory."""
out = []
for entry in os.scandir(dirname):
if entry.is_file():
match = REGEXP_BACKUP.search(entry.name)
if match is not None:
out.append(entry.name)
return out
[docs]
@staticmethod
def _read_configuration(filename):
"""
Read CP2K output configuration.
This method is used when we calculate the order parameter.
Parameters
----------
filename : string
The file to read the configuration from.
Returns
-------
box : numpy.array
The box dimensions if we manage to read it.
xyz : numpy.array
The positions.
vel : numpy.array
The velocities.
names : list of strings
The atom names found in the file.
"""
xyz, vel, box, names = None, None, None, None
for snapshot in read_xyz_file(filename):
box, xyz, vel, names = convert_snapshot(snapshot)
break # Stop after the first snapshot.
return box, xyz, vel, names
[docs]
def _reverse_velocities(self, filename, outfile):
"""Reverse velocity in a given snapshot.
Parameters
----------
filename : string
The configuration to reverse velocities in.
outfile : string
The output file for storing the configuration with
reversed velocities.
"""
box, xyz, vel, names = self._read_configuration(filename)
write_xyz_trajectory(outfile, xyz, -1.0*vel, names, box, append=False)
[docs]
def _prepare_shooting_point(self, input_file): # pragma: no cover
"""
Create initial configuration for a shooting move.
This creates a new initial configuration with random velocities.
2022.05.25: Note: This function is no longer in use
as we now generate velocities internally. However we keep
this function in case of future need.
Parameters
----------
input_file : string
The input configuration to generate velocities for.
Returns
-------
output_file : string
The name of the file created.
energy : dict
The energy terms read from the CP2K energy file.
"""
box, xyz, vel, atoms = self._read_configuration(input_file)
if box is None:
box, _ = read_cp2k_box(self.input_files['template'])
input_config = os.path.join(self.exe_dir, 'genvel_input.xyz')
write_xyz_trajectory(input_config, xyz, vel, atoms, box, append=False)
# Create input file for CP2K:
run_file = os.path.join(self.exe_dir, 'run.inp')
# draw_integer_inclusive handles both the native RandomGenerator
# and the scheduler's bare numpy Generator (no random_integers).
from pyretis.core.common import draw_integer_inclusive
write_for_genvel(self.input_files['template'],
run_file,
'genvel_input.xyz',
draw_integer_inclusive(self.rgen, 1, 999999999),
name='genvel')
# Prepare to run it:
self.add_input_files(self.exe_dir)
out_files = self.run_cp2k(run_file, 'genvel')
energy = read_cp2k_energy(out_files['energy'])
# Get the output configuration:
atoms, xyz, vel, box, _ = read_cp2k_restart(out_files['restart'])
conf_out = os.path.join(self.exe_dir, f'genvel.{self.ext}')
write_xyz_trajectory(conf_out, xyz, vel, atoms, box, append=False)
# Remove run-files etc:
for _, files in out_files.items():
self._removefile(files)
self._remove_files(
self.exe_dir,
self._find_backup_files(self.exe_dir)
)
return conf_out, energy
[docs]
def modify_velocities(self, system, vel_settings=None, rgen=None):
"""
Modify the velocities of the current state.
This method will modify the velocities of a time slice.
Parameters
----------
system : object like :py:class:`.System`
System is used here since we need access to the particle list.
vel_settings: dict, optional
It contains info about the velocity settings:
* `sigma_v` : numpy.array
Negative values select aimless shooting. Zero or positive
values set the standard deviation, one for each particle,
for soft velocity perturbations.
* `momentum` : boolean
If True, we reset the linear momentum to zero after
generating.
* `rescale` : float
In some NVE simulations, we may wish to re-scale the
energy to a fixed value. If `rescale` is a float > 0,
we will re-scale the energy (after modification of
the velocities) to match the given float.
rgen : object like :py:class:`.RandomGenerator`
This is the random generator that will be used.
Returns
-------
dek : float
The change in the kinetic energy.
kin_new : float
The new kinetic energy.
"""
# Engine-level masses + inverse thermal energy drive the Maxwell
# draw and the kinetic energy, so the draw does not depend on the
# cp2k system carrying masses/temperature (the GROMACS/LAMMPS
# engines work the same way). Lazy import to avoid the
# _parts -> setup circular import at module load.
from pyretis.engines._parts import (
kinetic_energy as _kinetic_energy,
reset_momentum as _reset_momentum,
)
self._ensure_velocity_setup()
mass = self.mass
beta = self._beta
rescale = vel_settings.get('rescale')
pos = self.dump_frame(system)
box, xyz, vel, atoms = self._read_configuration(pos)
system.vel = vel
if box is None:
box, _ = read_cp2k_box(self.input_files['template'])
# to-do: retrieve system.vpot from previous energy file.
if None not in ((rescale, system.vpot)) and rescale is not False:
if rescale > 0:
kin_old = rescale - system.vpot
do_rescale = True
else:
logger.warning('Ignored re-scale 6.2%f < 0.0.', rescale)
return 0.0, _kinetic_energy(vel, mass)[0]
else:
kin_old = _kinetic_energy(vel, mass)[0]
do_rescale = False
if is_aimless_velocity_setting(vel_settings):
vel, _ = self.draw_maxwellian_velocities(vel, mass, beta, rgen)
else:
dvel, _ = self.draw_maxwellian_velocities(
vel, mass, beta, rgen, sigma_v=vel_settings['sigma_v'])
vel = vel + dvel
if vel_settings.get('momentum', False):
vel = _reset_momentum(vel, mass)
system.vel = vel
if do_rescale:
# rescale_system_velocities reads the system's masses; the cp2k
# system carries none, so populate the engine-level masses here
# (same source as the draw above) to keep the energy rescale
# consistent with the Maxwell draw.
system.particles.mass = mass
system.particles.imass = np.reciprocal(mass * 1.0)
rescale_system_velocities(system, rescale, external=True)
vel = system.vel
conf_out = os.path.join(self.exe_dir, f'genvel.{self.ext}')
write_xyz_trajectory(conf_out, xyz, vel,
atoms, box, append=False)
kin_new = _kinetic_energy(vel, mass)[0]
system.config = (conf_out, None)
system.vel_rev = False
system.ekin = kin_new
if kin_old == 0.0:
dek = float('inf')
logger.debug(('Kinetic energy not found for previous point.'
'\n(This happens when the initial configuration '
'does not contain energies.)'))
else:
dek = kin_new - kin_old
return dek, kin_new