[ase-users] Save results of dynamic simulation of atoms position, energy, forces

[Dong gi Kang]

Hi all,

I would like to store as an extended xyz file the results of dynamic simulations (atoms position, energy, forces).
This is the script that I ran MD calcs for Si slab. The end of the script (where is marked) only stores atom position and momenta instead of position/energy/force (atoms.copy() does not copy force/energy)

Could you please give any advice?

` ` `
import os
import sys
import numpy as np
import matplotlib.pylab as plt
from quippy.potential import Potential
from quippy.descriptors import Descriptor
from ase import Atoms, units
from ase.build import add_vacuum
from ase.lattice.cubic import Diamond
from ase.io import write
from ase.constraints import FixAtoms
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase.md.verlet import VelocityVerlet
from ase.md.langevin import Langevin
from ase.optimize.precon import PreconLBFGS, Exp
sys.path.append('/home/{user}/software/QUIP/src/GAP/doc_src/gap_si_surface/')
from visualise import ViewStructure

T = 1000.0 # Temperature [Kelvin]
timestep = 1.0 * units.fs

##############
# building struc #
##############
def build_slab(size=(1,2,2), vacuum=10.):
    # Build Si lattice.
    # lattice = Diamond('Si', directions=([1, 0, 0], [0, 1, 0], [0, 0, 1]), latticeconstant=5.44, size=size)
    lattice = Diamond('Si', latticeconstant=5.44, size=size)
    atoms = Atoms(lattice)

    # Fixing the bottom layer
    bottom = atoms.positions[:,2].min()
    fixed_mask = (abs(atoms.positions[:,2] - bottom) < 2.0)
    atoms.set_constraint(FixAtoms(mask=fixed_mask))

    # build surface by adding space to z direction
    add_vacuum(atoms, vacuum)
    # atoms.center(vacuum=10.0, axis=2)

    return atoms

atoms = build_slab()
print('Number of atoms:', len(atoms))

############
# Setting Calc #
############
# attach tight binding calculator
qm_pot = Potential('TB DFTB', param_filename='/home/uccatka/software/QUIP/share/Parameters/tightbind.parms.DFTB.pbc-0-1.xml')
atoms.set_calculator(qm_pot)  # Set the potential

# Thermalize atoms
MaxwellBoltzmannDistribution(atoms, 2.0 * T * units.kB)
# dynamics = VelocityVerlet(atoms, timestep)
dynamics = Langevin(atoms, timestep, T * units.kB, 0.002)

#attach observer to dynamics to track quantity of interest (temperature)
def print_status():
    print('Step = {}, time = {} [fs], T = {} [K]'.format(
        dynamics.nsteps,
        dynamics.nsteps * dynamics.dt / units.fs,
        atoms.get_kinetic_energy() / (1.5 * units.kB * len(atoms))
    ))

def print_energy(a=atoms):  # store a reference to atoms in the definition.
    """Function to print the potential, kinetic and total energy."""
    epot = a.get_potential_energy() / len(a)
    ekin = a.get_kinetic_energy() / len(a)
    force = a.get_forces()
    print('Energy per atom: Epot = %.3feV  Ekin = %.3feV (T=%3.0fK)  '
          'Etot = %.3feV force = %.3f'% (epot, ekin, ekin / (1.5 * units.kB), epot + ekin, force))

dynamics.attach(print_status, interval=10)
dynamics.attach(print_energy, interval=1)
dynamics.run(steps=100)

##################
# store position/E/F  #
##################
db = []
def collect_data():
    db.append(atoms.copy())

dynamics.attach(collect_data, interval=10)
dynamics.run(steps=500)

write(‘./tmp/atoms_db.xyz’, db, write_results=True)
` ` `
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://listserv.fysik.dtu.dk/pipermail/ase-users/attachments/20201210/47f6c914/attachment.html>