Molecular Dynamics Simulation

DFTpy performs molecular dynamics (MD) simulations with ASE. This is one example to run NVT (canonical ensemble) simulation:

import os
import numpy as np
from ase.lattice.cubic import FaceCenteredCubic
from ase.md.langevin import Langevin
from ase.md.verlet import VelocityVerlet
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase.io.trajectory import Trajectory
from ase import units
from ase.md.npt import NPT

from dftpy.config import DefaultOption, OptionFormat
from dftpy.interface import OptimizeDensityConf
from dftpy.api.api4ase import DFTpyCalculator

############################## initial config ##############################
conf = DefaultOption()
conf['PATH']['pppath'] = os.environ.get('DFTPY_DATA_PATH') 
conf['PP']['Al'] = 'Al_lda.oe01.recpot'
conf['OPT']['method'] = 'TN'
conf['KEDF']['kedf'] = 'WT'
conf['JOB']['calctype'] = 'Energy Force'
conf = OptionFormat(conf)
#-----------------------------------------------------------------------

size = 3
a = 4.24068463425528
T = 1023  # Kelvin
T *= units.kB
atoms = FaceCenteredCubic(directions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]],
                          latticeconstant = a,
                          symbol="Al",
                          size=(size, size, size),
                          pbc=True)

calc = DFTpyCalculator(config = conf)
atoms.set_calculator(calc)

MaxwellBoltzmannDistribution(atoms, T, force_temp = True)

dyn = Langevin(atoms, 2 * units.fs, T, 0.1)

step = 0
interval = 1
def printenergy(a=atoms):  
    global step, interval
    epot = a.get_potential_energy() / len(a)
    ekin = a.get_kinetic_energy() / len(a)
    print('Step={:<8d} Epot={:.5f} Ekin={:.5f} T={:.3f} Etot={:.5f}'.format(step, epot, ekin, ekin / (1.5 * units.kB), epot + ekin))
    step += interval
    
dyn.attach(printenergy, interval=1)

traj = Trajectory('md.traj', 'w', atoms)
dyn.attach(traj.write, interval=5)

dyn.run(100)

The output md.traj can be converted with ASE:

#!/usr/bin/env python3

from ase.io.trajectory import Trajectory
from ase.io import write

infile = sys.argv[1]
outfile = sys.argv[2]

traj = Trajectory(infile)
write(outfile, traj)