Chemical potential and molecular dynamics

[ivca]

 Hi!

I would like to calculate and print the chemical potential for every step of a molecular dynamics (MD) simulation.
I have used one of the tutorials to generate the MD script and to print the chemical potential. Here the last lines of the script:
"
nlsave('output.nc', bulk_configuration)

# -------------------------------------------------------------
# Molecular Dynamics
# -------------------------------------------------------------

initial_velocity = None

method = NVTBerendsen(
    time_step=0.5*femtoSecond,
    reservoir_temperature=300*Kelvin,
    thermal_coupling_constant=500*femtoSecond,
    initial_velocity=initial_velocity
)

md_trajectory = MolecularDynamics(
    bulk_configuration,
    trajectory_filename='trajectory.nc',
    steps=100000,
    log_interval=2,
    method=method
)

bulk_configuration = md_trajectory.lastImage()
nlsave('output.nc', md_trajectory)

# -------------------------------------------------------------
# Chemical Potential
# -------------------------------------------------------------
chemical_potential = ChemicalPotential(bulk_configuration)
nlsave('output.nc', chemical_potential)
nlprint(chemical_potential)
"

In this way, the chemical potential is not printed out though.
Is it possible to print the chemical potential at each step? If yes, how?

Thank you in advance,
best regards,
Ivano

[Julian Schneider]

Yes, that is possible using a hook function during Molecular Dynamics. To do so, you first have to define a suitable hook function in your script, e.g.

Code

def chemical_potential_hook(step, time, configuration):
    """ 
    Hook function to calculate and print the chemical potential at 
    every MD step.
    """
    chemical_potential = ChemicalPotential(configuration)
    nlprint(chemical_potential)

Then, you define your MD simulation as usual, but with the additional argument

Code

post_step_hook=chemical_potential_hook

e.g.

Code

md_trajectory = MolecularDynamics(
    bulk_configuration,
    trajectory_filename='trajectory.nc',
    steps=100000,
    log_interval=2,
    post_step_hook=chemical_potential_hook,
    method=method
)

The hook function will now be called after each MD step and the specified operation, in this case the calculation and printing of the chemical potential, will be performed. Of course, the calculator that you use for running molecular dynamics must support the calculation of the chemical potential (i.e. it must be an electronic-structure-based calculator such as the LCAOCalculator, as this is the electronic chemical potential).