How to apply periodic boundary conditions in molecular dynamics

ASE
1

Hello!
I’m tying to apply periodic boundary condition in molecular dynamics.
I set atoms.set_pbc((True, True, True)) , but atoms move out of the cell.
image

Could you please tell me how to solve it?

source code is here.

from ase.lattice.cubic import FaceCenteredCubic
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase.md.verlet import VelocityVerlet
from ase import units
from ase.io.trajectory import Trajectory
from ase.io import read, write

from ase.calculators.emt import EMT
size = 2

atoms = FaceCenteredCubic(directions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]],
                          symbol='Cu',
                          size=(size, size, size),
                          pbc=True)
atoms.cell=([10,10,10])
atoms.center()
atoms.set_pbc((True, True, True))
print(atoms)

atoms.set_calculator(EMT())

MaxwellBoltzmannDistribution(atoms, 5000 * units.kB)

aaa="x123"
test1= aaa+ ".traj"

dyn = VelocityVerlet(atoms, 5 * units.fs, trajectory = test1)  # 5 fs time step.

for i in range(1000):
    dyn.run(1)
    

ff= Trajectory("x123.traj")
write("h2otraj_22.xyz", ff)
1 Like
2

As long as the ASE calculator is implemented correctly, the PBC should be applied just fine in the energy/force calculation.

Wrapping the positions in the saved trajectory is often undesirable as it makes it more difficult to correctly calculate velocities, correlation functions, diffusion etc. I would recommend that if you want to wrap the positions afterwards (e.g. for visualisation) you could do something like

for atoms in traj:
    atoms.wrap()

If it is really important to wrap the positions during the simulation (perhaps to work around a bug in the atomistic code?) then you can use an ASE feature to attach a callback function to the dynamics driver.

...
dyn = VelocityVerlet(atoms, 5 * units.fs, trajectory=test1)

def wrap_atoms(a=atoms):
    a.wrap()
dyn.attach(wrap_atoms, interval=1)

dyn.run(1000)

But this will come with some computational overhead!

3

Dear Ajjackson.

Thank you for the advice !
I was able to solve the problem.

I would like to ask an additional question.
How can I convert a trajectory to a xyz file?