Simulation backends for atooms¶

This is a demo of running a molecular dynamics simulation in the NVT ensemble for the 80-20 Kob-Ansersen binary mixture using the backend shipped with atooms-dynamics. We use the symplectic Nosé-Poincaré algorithm to equilibrate the system at normal liquid conditions. The interactions are computed using the built-in Fortran 90 code shipped with atooms-models.

This is explicitly written down as a monolitic piece of notebook.

Setup¶

Import the necessary modules and turn off verbosity

import numpy
import matplotlib
import matplotlib.pyplot as plt
import atooms.core
from atooms.system import Thermostat
from atooms.simulation import Simulation, Scheduler, write_trajectory, store
from atooms.dynamics import NosePoincare, VelocityVerlet
from atooms.trajectory import TrajectoryXYZ

atooms.core.progress.active = False
atooms.core.utils.setup_logging(level=40, update=True)

Create the starting configuration¶

We create a small system of \(N=100\) particles at the canonical density \(\rho=1.2\) and composition 80-20

from atooms.system import System
system = System(N=100)
system.density = 1.2
system.composition = {'A': 80, 'B': 20}
print(system)

system composed by N=100 particles
with chemical composition C={'A': 80, 'B': 20}
with chemical concentration x={'A': 0.8, 'B': 0.2}
enclosed in a cubic box at number density rho=1.200000

Define the interaction model¶

We look for the interaction model of the Kob-Andersen mixture in the atooms.models database.

from atooms import models
from pprint import pprint
pprint(models.get('kob_andersen'))

{'cutoff': [{'parameters': {'rcut': [[2.5, 2.0], [2.0, 2.2]]},
             'type': 'cut_shift'}],
 'potential': [{'parameters': {'epsilon': [[1.0, 1.5], [1.5, 0.5]],
                               'sigma': [[1.0, 0.8], [0.8, 0.88]]},
                'type': 'lennard_jones'}],
 'reference': 'W. Kob and H. C. Andersen, Phys. Rev. Lett. 73, 1376 (1994)'}

We define a new model, changing the cutoff to a smoother one, and set up the interaction

model = {'cutoff': [{'parameters': {'rcut': [[2.5, 2.0], [2.0, 2.2]]},
                     'type': 'linear_cut_shift'}],
         'potential': [{'parameters': {'epsilon': [[1.0, 1.5], [1.5, 0.5]],
                                       'sigma': [[1.0, 0.8], [0.8, 0.88]]},
                        'type': 'lennard_jones'}],
}
system.species_layout = 'F'  # This ensures species are Fortran-styled
system.interaction = models.f90.Interaction(model)
# system.interaction.neighbor_list = f90.VerletList(skin=0.3)

>>>

Run the simulation¶

We create a simulation backend that integrates the equation of motions using the Nosé-Poincaré algorithm at a temperature \(T=2.0\). The simulation s fairly short because at this temperature the system is a normal liquid and thermalies rapidly. We also add

system.set_temperature(2.0)
system.thermostat = Thermostat(temperature=2.0, mass=5.0)
bck = NosePoincare(system, timestep=0.002)
sim = Simulation(bck)
data = {}
# This adds a callback to store every 100 steps several system attributes in the data dictionary
# Check the atooms.simulation.observers.store() doc for more info.
sim.add(store, 100, ['steps', 'potential energy per particle',
                     'kinetic energy per particle', 'temperature'], data)

Run the simulation and show a snapshot of the final state of the system (this requires ovito)

sim.run(int(1e4))

# Temporarily set alphabetical species for ovito
system.species_layout = 'A'
# system.show('ovito', outfile='_images/config.png')
system.species_layout = 'F'

Check that the kinetic temperature looks good

plt.plot(data['steps'], data['temperature'], 'o', label='Kinetic temperature')
plt.plot(data['steps'], [system.thermostat.temperature]*len(data['steps']), '-', label='Thermostat temperature')
plt.xlabel('Steps')

We now do a production run in which we store configurations in a trajectory file, written in xyz format. This time we store the configurations using an exponential sampling, which allows to resolve both short and long time dynamics. This will take about a minute.

th = TrajectoryXYZ('/tmp/config.xyz', 'w')
bck = VelocityVerlet(system, timestep=0.006)
sim = Simulation(bck)
# This adds a callback to write the system in the th trajectory using an exponential schedule
sim.add(write_trajectory, Scheduler(block=[2**n for n in range(10)]),
                                    variables=['species', 'pos', 'vel'], trajectory=th)
data = {}
sim.add(store, 100, ['steps', 'potential energy per particle',
                     'total energy per particle', 'temperature'], data)
sim.run(int(1e5))
th.close()

Chech that the total energy is conserved

plt.plot(data['steps'], data['total energy per particle'], '-', label='Total')
plt.plot(data['steps'], data['potential energy per particle'], '-', label='Potential')
plt.xlabel('Steps')

Analyze the trajectory file¶

The trajectory file can be analyzed using the atooms-postprocessing package. We compute the velocity auto-correlation function \(Z(t)\) over a short time grid

from atooms.postprocessing import VelocityAutocorrelation

with TrajectoryXYZ('/tmp/config.xyz') as th:
    cf = VelocityAutocorrelation(th, tgrid=th.times[:9])
    cf.compute()
plt.plot(cf.grid, cf.value, '-o')
plt.xlabel('t')
plt.ylabel('Z(t)')

Here we compute instead the mean square displacement

from atooms.postprocessing import MeanSquareDisplacement

with TrajectoryXYZ('/tmp/config.xyz') as th:
    cf = MeanSquareDisplacement(th)
    cf.compute()
plt.plot(cf.grid[:20], cf.value[:20], '-o')
plt.xlabel('t')
plt.ylabel('MSD')

Dynamics