# How to apply periodic boundary conditions in molecular dynamics

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.

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.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)

``````

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!

Dear Ajjackson.

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