[gpaw-users] How to get the trajectory during TDDFT (ehrenfest md)

[Christian Glinsvad]

Hi Zhu Xi

It's a bit tricky because the Ehrenfest code is still in a refactoring state, so stuff like this is left up to the user to hack in.

Here's basically what I do:

import numpy as np
from ase.units import Hartree, Bohr
from ase.io.trajectory import PickleTrajectory
from ase.calculators.singlepoint import SinglePointCalculator
from gpaw.tddft import TDDFT
from gpaw.tddft.ehrenfest import EhrenfestVelocityVerlet

# Timestep and duration of simulation in attoseconds
timestep = 5.0
duration = 5e3

# Number of iterations and integer divisor to update stats every 0.1 fs
ndiv = int(np.ceil(0.1e3 / timestep))
niter = ndiv * int(np.ceil(duration / (ndiv * timestep)))

tdcalc = TDDFT('ground_state_file.gpw', solver='CSCG', propagator='EFSICN')
ehrenfest = EhrenfestVelocityVerlet(tdcalc)

traj = PickleTrajectory('ehrenfest.traj', 'w', tdcalc.get_atoms())

for i in range(niter):
    ehrenfest.propagate(timestep)

    if i % ndiv == 0:
        # Dummy atoms and calculator to prevent calls to GPAW::get_potential_energy
        spa = tdcalc.get_atoms()
        epot = tdcalc.get_td_energy() * Hartree
        F_av = ehrenfest.F * Hartree / Bohr
        spc = SinglePointCalculator(epot, F_av, None, None, spa)
        spa.set_calculator(spc)
        traj.write(spa)

traj.close()


Note that I'm using the Ehrenfest branch in GPAW for this...
https://trac.fysik.dtu.dk/projects/gpaw/browser/branches/ehrenfest

Regards
Christian Glinsvad


Zhu Xi wrote:
>Dear All;
>
>How can I  get the trajectory imformation during TDDFT (ehrenfest md)?
>
>do I need import the ase-dyn module?
>
>Thanks