Combined dynamics module in ASE

ASE
1

Hi there,

I wonder whether it’s possible to create a new dynamic module combined the VelocityVerlet (all atoms) and the Langevin (only to a group of atoms, with indices given, serves as the thermostat layer) in ASE?

Many thanks,

Catherine

2

This seems the sort of thing that a user might write for a particular research application, and share the code with their publication. Is it a routine method in some field?

Perhaps it could be implemented in a general sense by creating a MolecularDynamics class that contains instances of other MolecularDynamics classes, wrapping copies of the Atoms. When its step() method is called, it gets forces from the “real” Atoms, applies them to the copies (via SinglePointCalculator), calls each of their step() methods, and then recombines the consequences as desired onto the “real” Atoms.

3

Yes this approach can help control the system temperature (from a thermostat layer), at the same time prevent the research of interest region out of thermal fluctuation.

This sounds doable, and the structure probably will be as follows?

class VerletLangevin(VelocityVerlet):
    def __init__(self, atoms, timestep, thermostat_indices, temperature_K, damping, **kwargs):
   
    def apply_langevin(self):

    def step(self):
        super().step()
        self.apply_langevin()

dyn = VerletLangevin(...)

Cheers,

Catherine

4

Perhaps! I had in mind something more like

class MultiThermostat(MolecularDynamics):
    def __init__(self, atoms, timestep,
        thermostat_2_indices: list[int],
        thermostat_1_class: type[MolecularDynamics],
        thermostat_1_kwargs: dict[str, Any],
        thermostat_2_class: type[MolecularDynamics],
        thermostat_2_kwargs: dict[str, Any]
    ):
        self.thermostat_1 = thermostat_1_class(atoms.copy(), timestep=timestep, **thermostat_1_kwargs)
        self.thermostat_2 = thermostat_2_class(atoms.copy(), timestep=timestep, **thermostat_2_kwargs)
        ...

    def step(self):
        self.thermostat_1.atoms.set_positions(self.atoms.positions)
        self.thermostat_1.atoms.set_velocities(
            self.atoms.get_velocities())
        self.thermostat_1.atoms.calc = SinglePointCalculator(
            atoms=self.thermostat_1.atoms,
            forces=self.atoms.get_forces())

        self.thermostat_2.atoms.set_positions(self.atoms.positions)
        self.thermostat_2.atoms.set_velocities(
            self.atoms.get_velocities())
        self.thermostat_2.atoms.calc = SinglePointCalculator(
            atoms=self.thermostat_2.atoms,
            forces=self.atoms.get_forces())

        self.thermostat_1.step()
        self.atoms.positions = self.thermostat_1.atoms.positions
        self.atoms.set_velocities(
            self.thermostat_1.atoms.get_velocities())
        self.thermostat_2.step()
        self.atoms.positions[self._thermostat_2_indices] = self.thermostat_2.atoms.positions[self.thermostat_2_indices]
        self.atoms[self._thermostat_2_indices].set_velocities(
            self.thermostat_2.atoms.get_velocities()
            )[self._thermostat_2_indices]

Probably I am missing some other essential things :sweat_smile:

5

Thanks Adam, I will give it a try. Cheers, Catherine