Running Polyhedral Template Matching

This is almost the same script as the previous example (Running Common Neighbor Analysis), but using the new Polyhedral Template Matching algorithm instead. If you run both scripts you will notice that PTM sees an abrupt crystallization of the cluster at a temperature where CNA still is not able to distiguish much structure in the cluster.

The final structure below looks very different from the one displayed with the previous example. This is because the two simulations produced different final clusters; the Langevin dynamics is not deterministic. At 300K the performance of CNA and PTM are very similar.

A cut through the cluster. FCC atoms are white, HCP atoms red and all other atoms are blue.

PTMdemo.py

"Melting a copper cluster."

from numpy import *
from asap3 import Atoms, EMT, units
from asap3.analysis import PTMobserver
from asap3.visualize.primiplotter import *
from ase.lattice.cubic import FaceCenteredCubic
from asap3.md.langevin import Langevin

# Color coding for PTM analysis
mycolors = {0: (0.3, 0.3, 1.0),  # Unknown is blue
    1: (0.2, 1.0, 0.2),  # SC is green
    2: (1.0, 1.0, 1.0),  # FCC is white
    3: (1.0, 0.0, 0.0),  # HCP is red
    4: (1.0, 1.0, 0.0),  # ICO is yellow
    5: (0.0, 0.0, 0.0)   # BCC is black
    }

# Create the atoms
atoms = FaceCenteredCubic(size=(6,6,6), symbol="Cu", pbc=False)

# Associate the EMT potential with the atoms
atoms.set_calculator(EMT())

# Temperature profile - brief run at 1750, then lower temperature gradually
t = linspace(1500, 250, 81)
temperatures = 1500 * ones(100)
temperatures[-len(t):] = t
print(temperatures)

# How many steps at each temperature
nsteps_total = 100000
nsteps = nsteps_total // len(temperatures)

# Interval between plots
plotinterval = 2000

# Make the Langevin dynamics module
dyn = Langevin(atoms, 5*units.fs, units.kB*temperatures[0], 0.002)

# The PTM analyser is called every plotinterval timesteps.
ptm = PTMobserver(atoms, rmsd_max=0.16)
dyn.attach(ptm.analyze, interval=plotinterval)

# The plotter
def invisible_atoms(a):
    """Return True for atoms that should be invisible."""
    r = atoms.get_positions()
    centerofmass = r.sum(axis=0) / len(atoms)
    return (r[:,2] < centerofmass[2])

plotter = PrimiPlotter(atoms)
plotter.set_invisibility_function(invisible_atoms)
plotter.set_colors(mycolors) # Map tags to colors
plotter.set_output(PngFile("ptm"))  # Save plots in files ptm0000.png ...
plotter.set_rotation((10.0, 5.0, 0))

# Attach the plotter to the PTMobserver object.  That guarantees
# that the plotter is called AFTER the PTM analysis has been done.
# Similarly, a Trajectory should be attached to the PTMobserver
# object.  By using interval=1 (the default), the plotter is called
# every time PTMobserver is called, i.e. every plotinterval
# timesteps.
ptm.attach(plotter.plot)

# The main loop
for t in temperatures:
    dyn.set_temperature(units.kB*t)
    for i in range(nsteps//100):
        dyn.run(100)
        print("E_total = %-10.5f  T = %.0f K  (goal: %.0f K, step %d of %d)" %\
          (atoms.get_total_energy()/len(atoms), atoms.get_temperature(), t, i, nsteps/100))