Examples of usage of ASE package#
Below you can find examples of how to perform automated simulations using the Atomic Simulation Environment (ASE) Python package together with SAMSON.
The Atomic Simulation Environment (ASE) is a set of tools and Python modules for setting up, manipulating, running, visualizing and analyzing atomistic simulations. [ASE docs]
Requirements
First, you need to install the ASE package in the integrated Python in SAMSON. Please see Installing Python packages - type ase and click Install.
Optimization of a water molecule using ASE#
In this example, we create a water molecule using the ASE package, copy it to SAMSON, optimize it using ASE (with BFGS), and update atoms positions in SAMSON accordingly.
For the sake of this example, we will run the optimization for several steps and we will be updating atoms’ positions in SAMSON according to their change in ASE after each step.
Optimize a water molecule using ASE#
import numpy as np
from ase import Atoms
from ase.optimize import BFGS
from ase.calculators.emt import EMT
class ExampleASEOptimizationH20:
molecule = None
atom_indexer = None
def __init__(self):
"""Initialize ASE and SAMSON"""
self.initialize_system_ase() # initialize the system in ASE
self.initialize_system_samson() # initialize the system in SAMSON
def compute(self):
"""Do computations using ASE and update positions in SAMSON"""
dyn = BFGS(self.molecule) # set up computations in ASE
nsteps = 100
for i in range(nsteps):
dyn.run(fmax = 0.005, steps = 20) # run the computations in ASE
self.update_positions_in_samson() # update positions in SAMSON
def initialize_system_ase(self):
"""Create an H2O molecule using ASE package; set a calculator"""
# the H2O molecule parameters: distance and angle
d = 0.9575
t = np.pi / 180.0 * 104.51
# create a molecule in ASE (positions in ASE are in Angstrom); set calculator
self.molecule = Atoms('H2O',
positions = [(d, 0, 0.1),
(d * np.cos(t), d * np.sin(t), 0),
(0, 0, 0)],
calculator = EMT())
def initialize_system_samson(self):
"""Adds atoms created by ASE into the active document in SAMSON"""
structural_model = SBStructuralModel() # a new structural model
structural_model.name = 'H20 optimization example'
# add atoms in the structural model
self.add_atoms_in_samson(structural_model, self.molecule)
# add covalent bonds
structural_model.createCovalentBonds()
# enable undo/redo for the creation of the structural model
if SAMSON.isHolding(): SAMSON.hold(structural_model)
# create this structural model in SAMSON
structural_model.create()
# add the structural model into the active document
SAMSON.getActiveDocument().addChild(structural_model)
# get an indexer of all atoms in the structural_model
self.atom_indexer = structural_model.getNodes('n.t a')
def add_atoms_in_samson(self, structural_model, molecule):
"""Add atoms in SAMSON according to their type and positions"""
for a in molecule: # loop through a list of atoms created by ase
# construct an atom in SAMSON with a given element type and position
sa = SBAtom(SAMSON.getElementTypeBySymbol(a.symbol),
SBQuantity.angstrom(a.position[0]),
SBQuantity.angstrom(a.position[1]),
SBQuantity.angstrom(a.position[2]))
# add an atom to the root of the structural model
structural_model.addChild( sa )
def update_positions_in_samson(self):
"""Update positions of atoms in SAMSON based on new positions of these atoms in ASE"""
# get the new positions from ASE
new_positions = self.molecule.get_positions()
for i in range(len(self.atom_indexer)):
atom = self.atom_indexer[i]
atom.setX(SBQuantity.angstrom(new_positions[i, 0]))
atom.setY(SBQuantity.angstrom(new_positions[i, 1]))
atom.setZ(SBQuantity.angstrom(new_positions[i, 2]))
# process all events in SAMSON (to update positions in the viewport)
SAMSON.processEvents()
def store_current_conformation(self, name):
"""Create a conformation in SAMSON based
on the current state of the structural model
and add it to the active document"""
if len(self.atom_indexer):
conformation = SBConformation(name, self.atom_indexer)
if SAMSON.isHolding(): SAMSON.hold(conformation)
conformation.create()
SAMSON.getActiveDocument().addChild(conformation)
# make the operation undoable
SAMSON.beginHolding('Water optimization example')
# Create an object of the ExampleASEOptimizationH20 class, which will initialize ASE and SAMSON.
# After this step, you will see a water molecule in the Viewport and in the Document View.
ase_h20_opt = ExampleASEOptimizationH20()
# store the initial conformation in the active document
ase_h20_opt.store_current_conformation('Initial conformation')
# set a camera view in SAMSON
SAMSON.getActiveCamera().topView()
# process events: updates viewport
SAMSON.processEvents()
# run the computations
ase_h20_opt.compute()
# process events: updates viewport
SAMSON.processEvents()
# store the optimized conformation in the active document
ase_h20_opt.store_current_conformation('Optimized conformation')
# stop the undoable operation
SAMSON.endHolding()
In the console output you may get something like this:
Step Time Energy fmax
BFGS: 0 15:02:48 2.733104 8.4006
BFGS: 1 15:02:48 1.925778 1.5675
...
BFGS: 10 15:02:48 1.878888 0.0139
BFGS: 11 15:02:48 1.878884 0.0018
The script also creates the initial and final conformations of the system - you can switch between them by double-clicking on them in the Document view or in Python via SBConformation.restore().
Computing adsorption energy for Nitrogen on Copper system#
In this example, we create using ASE a system consisting of Nitrogen molecule and a slab of Copper atoms using ASE package, copy it to SAMSON, compute adsorption energy, simulate it in ASE using velocity Verlet, and update atoms positions in SAMSON accordingly.
This example is based on the ASE tutorial: Surface adsorption study using the ASE database.
Optimize a system using ASE#
from ase import Atoms
from ase.calculators.emt import EMT
from ase.constraints import FixAtoms
from ase.optimize import QuasiNewton
from ase.build import fcc111, add_adsorbate
from ase.md.verlet import VelocityVerlet
from ase import units
class ExampleASEOptimizationN2Cu:
slab = None
molecule = None
atom_indexer = None
def __init__(self):
"""Initialize ASE and SAMSON"""
self.initialize_system_ase() # initialize the system in ASE
self.initialize_system_samson() # initialize the system in SAMSON
def compute_adsorption_energy(self):
"""Computations using ASE; update positions in SAMSON
Optimize the structure of the N2 molecule adsorbed on the Cu surface"""
# set up computations in ASE
h = 1.85
add_adsorbate(self.slab, self.molecule, h, 'ontop')
constraint = FixAtoms(mask = [a.symbol != 'N' for a in self.slab])
self.slab.set_constraint(constraint)
dyn = QuasiNewton(self.slab
# if you would like to save the trajectory,
# uncomment and modify the line below
#, trajectory = '/path/to/N2Cu.traj'
)
dyn.run(fmax = 0.05)
# update positions in SAMSON
self.update_positions_in_samson()
print('Adsorption energy:', self.e_slab + self.e_N2 - self.slab.get_potential_energy())
def simulate(self):
"""MD computations using ASE; update positions in SAMSON"""
# set up computations in ASE
dyn = VelocityVerlet(self.molecule, timestep = 1.0 * units.fs)
nsteps = 100
for i in range(nsteps):
pot = self.molecule.get_potential_energy()
kin = self.molecule.get_kinetic_energy()
print('%2d: Etotal = %.5f eV, Epot = %.5f eV, Ekin = %.5f eV' %
(i, pot + kin, pot, kin))
dyn.run(steps = 20) # run the computations in ASE
self.update_positions_in_samson() # update positions in SAMSON
def initialize_system_ase(self):
"""Create molecules using ASE package; set a calculator"""
# the slab parameters
d = 1.10
# create atoms using ASE (positions in ASE are in Angstrom); set calculator
self.slab = fcc111('Cu', size = (4, 4, 2), vacuum = 10.0) # build Cu crystal
self.slab.set_calculator(EMT())
self.e_slab = self.slab.get_potential_energy()
self.molecule = Atoms('2N', positions = [(0., 0., 0.), (0., 0., d)])
self.molecule.set_calculator(EMT())
self.e_N2 = self.molecule.get_potential_energy()
def initialize_system_samson(self):
"""Adds atoms created by ASE into a new layer in SAMSON"""
structural_model = SBStructuralModel() # a new structural model
structural_model.name = 'ASE slab optimization'
# add slab atoms in SAMSON
self.add_atoms_in_samson(structural_model, self.slab)
# add molecule atoms in SAMSON
self.add_atoms_in_samson(structural_model, self.molecule)
# enable undo/redo for the creation of the structural model
if SAMSON.isHolding(): SAMSON.hold(structural_model)
# create this structural model in SAMSON
structural_model.create()
# add the structural model into the active document
SAMSON.getActiveDocument().addChild(structural_model)
# get an indexer of all Nitrogen atoms in the structural_model
self.atom_indexer = structural_model.getNodes('n.t a and N')
def add_atoms_in_samson(self, structural_model, molecule):
"""Add atoms in SAMSON according to their type and positions"""
for a in molecule: # loop through a list of atoms created by ase
# construct an atom in SAMSON with a given element type and position
sa = SBAtom(SAMSON.getElementTypeBySymbol(a.symbol),
SBQuantity.angstrom(a.position[0]),
SBQuantity.angstrom(a.position[1]),
SBQuantity.angstrom(a.position[2]))
# add an atom to the root of the structural model
structural_model.addChild( sa )
def update_positions_in_samson(self):
"""Update positions of atoms in SAMSON based on new positions of these atoms in ASE"""
# get the new positions from ASE
new_positions = self.molecule.get_positions()
for i in range(len(self.atom_indexer)):
atom = self.atom_indexer[i]
atom.setX(SBQuantity.angstrom(new_positions[i, 0]))
atom.setY(SBQuantity.angstrom(new_positions[i, 1]))
atom.setZ(SBQuantity.angstrom(new_positions[i, 2]))
# process all events in SAMSON (to update positions in the viewport)
SAMSON.processEvents()
def store_current_conformation(self, name):
"""Create a conformation in SAMSON based
on the current state of the structural model
and add it to the active document"""
if len(self.atom_indexer):
conformation = SBConformation(name, self.atom_indexer)
if SAMSON.isHolding(): SAMSON.hold(conformation)
conformation.create()
SAMSON.getActiveDocument().addChild(conformation)
# make the operation undoable
SAMSON.beginHolding('Water optimization example')
# Create an object of the ExampleASEOptimizationN2Cu class, which will initialize ASE and SAMSON.
# After this step, you will see Nitrogen and Copper slab in the Viewport and in the Document View.
ase_N2Cu_opt = ExampleASEOptimizationN2Cu()
# store the initial conformation in the active document
ase_N2Cu_opt.store_current_conformation('Initial conformation')
# set a camera view in SAMSON
SAMSON.getActiveCamera().rightView()
# process events: updates viewport
SAMSON.processEvents()
# compute the adsorption energy
ase_N2Cu_opt.compute_adsorption_energy()
# perform MD simulation
ase_N2Cu_opt.simulate()
# process events: updates viewport
SAMSON.processEvents()
# store the optimized conformation in the active document
ase_N2Cu_opt.store_current_conformation('Optimized conformation')
# stop the undoable operation
SAMSON.endHolding()
In the console output you may get something like this:
Step[ FC] Time Energy fmax
*Force-consistent energies used in optimization.
BFGSLineSearch: 0[ 0] 15:01:59 11.689927* 1.0797
BFGSLineSearch: 1[ 2] 15:01:59 11.670814* 0.4090
BFGSLineSearch: 2[ 4] 15:01:59 11.625880* 0.0409
Adsorption energy: 0.32351942231800734
0: Etotal = 0.44034 eV, Epot = 0.44034 eV, Ekin = 0.00000 eV
1: Etotal = 0.43816 eV, Epot = 0.26289 eV, Ekin = 0.17527 eV
...
98: Etotal = 0.43999 eV, Epot = 0.40155 eV, Ekin = 0.03843 eV
99: Etotal = 0.43858 eV, Epot = 0.31076 eV, Ekin = 0.12782 eV
During this short simulation, you may see in the SAMSON’s viewport how these Nitrogen atoms move atop of the Copper slab.
The script also creates the initial and final conformations of the system - you can switch between them by double-clicking on them in the Document view or in Python via SBConformation.restore().
If you want to visualize the saved trajectory using ASE, you can save it in the compute_adsorption_energy function (uncomment the corresponding line) and then you can use the following code:
Visualize saved trajectory using ASE#
from ase.io import read
from ase.visualize import view
slab_mol = read('/path/to/N2Cu.traj')
view(slab_mol)
Combining SAMSON simulation with ASE calculators#
This example demonstrates how you can combine a simulation in SAMSON with ASE calculators. We use an interaction model provided by SAMSON together with ASE calculators for minimization of a system. The interaction model will be used to compute forces in SAMSON and ASE will be used for minimization of the system by using FIRE. For that, we implement our own calculator in ASE (see ASE calculators).
ASE provides interfaces to different codes through Calculators which are used together with the central Atoms object and the many available algorithms in ASE. [ASE docs]
Note
SAMSON also provides the FIRE state update.
Open some hydrocarbon from the Samples folder (e.g. cyclopropane.pdb) or create one from assets. If you want, you can slightly move some atoms in the molecule using SAMSON editors, e.g., choose the Displacer move editor (shortcut: D) and move slightly some atoms.
In the example below we apply to models in the active document a simulator with the Interaction model set to “Brenner” and the State Updater set to “Interactive modeling”. To interact between SAMSON interaction model and ASE calculator we create the InteractionModelWrapper class based on ASE calculator.
Optimize a system using ASE#
import os
import numpy as np
from ase import Atoms, Atom
from ase.optimize import FIRE
from ase.calculators.calculator import Calculator, all_changes
class InteractionModelWrapper(Calculator):
"""Calculator for ASE, see: https://wiki.fysik.dtu.dk/ase/development/calculators.html"""
# properties that calculator can handle
implemented_properties = ['energy', 'forces']
interaction_model = None
molecule = None
atom_indexer = None
def __init__(self, **kwargs):
Calculator.__init__(self, **kwargs)
#print("Setting Interaction Model wrapper...")
if kwargs is not None:
self.interaction_model = kwargs["interaction_model"]
self.molecule = kwargs["molecule"]
self.atom_indexer = kwargs["atom_indexer"]
else:
print("Error: input parameters are not provided")
pass
if self.interaction_model == None or self.molecule == None or self.atom_indexer == None:
print("Error: improper input parameters")
pass
def calculate(self, atoms = None,
properties = ['energy', 'forces'], system_changes = all_changes):
"""The main function for Calculator that does the calculations"""
# update interactions in SAMSON
self.interaction_model.updateInteractions()
forces = np.zeros((len(self.molecule), 3))
# get forces from the SAMSON's interaction model
for i in range(len(self.atom_indexer)):
# get force in Newtons
force = self.interaction_model.getForce(i)
forces[i, 0] = force[0].N.value
forces[i, 1] = force[1].N.value
forces[i, 2] = force[2].N.value
# set calculated properties for ASE
# convert forces from N to atomic unit force, 1 atomic unit force = Eh / a0
force_convertion = 1. / (SBQuantity.Eh(1).J / SBQuantity.a0(1).m).value
self.results['forces'] = forces * force_convertion
# get energy from the SAMSON's interaction model
self.results['energy'] = self.interaction_model.energy.eV.value
# update positions of atoms in SAMSON
self.update_positions_in_samson()
def update_positions_in_samson(self):
"""Update positions of atoms in SAMSON based on new positions of these atoms in ASE"""
# get the new positions from ASE
new_positions = self.molecule.get_positions()
for i in range(len(self.atom_indexer)):
atom = self.atom_indexer[i]
atom.setX(SBQuantity.angstrom(new_positions[i, 0]))
atom.setY(SBQuantity.angstrom(new_positions[i, 1]))
atom.setZ(SBQuantity.angstrom(new_positions[i, 2]))
# process all events in SAMSON (to update positions in the viewport)
SAMSON.processEvents()
class ExampleASECalculator:
"""An example of the usage of a SAMSON's InteractionModelParticleSystem together with ASE.
We assume, that the InteractionModelParticleSystem has already been created in SAMSON
for the atoms in the SAMSON's data graph.
"""
molecule = None
atom_indexer = None
def __init__(self):
# get an indexer of all atoms in SAMSON's active document
self.atom_indexer = SAMSON.getNodes('node.type atom')
# create atoms in ASE
self.create_atoms_in_ase(self.atom_indexer)
# get all interaction models from the active document in SAMSON
imps_indexer = SAMSON.getNodes('n.t imps')
# get the first interaction model
interaction_model = imps_indexer[0].castToInteractionModelParticleSystem()
# set an ASE calculator
self.molecule.set_calculator(
InteractionModelWrapper(interaction_model = interaction_model,
molecule = self.molecule,
atom_indexer = self.atom_indexer)
)
def create_atoms_in_ase(self, atom_indexer):
"""Create atoms in ASE based on atoms in SAMSON"""
# create an empty molecule in ASE
self.molecule = Atoms()
# add atoms from SAMSON to the molecule in ASE
for a in atom_indexer:
self.molecule.append(
Atom(a.elementSymbol,
(a.getX().angstrom.value, a.getY().angstrom.value, a.getZ().angstrom.value)
) )
def compute(self):
# set up the computation
dyn = FIRE(self.molecule)
# run the computation
dyn.run(fmax = 0.005, steps = 100)
def store_current_conformation(self, name):
"""Create a conformation in SAMSON based
on the current state of the structural model
and add it to the active document"""
if len(self.atom_indexer):
conformation = SBConformation(name, self.atom_indexer)
if SAMSON.isHolding(): SAMSON.hold(conformation)
conformation.create()
SAMSON.getActiveDocument().addChild(conformation)
def samson_add_simulator(simulator):
"""Initialize and add simulator in the active document in SAMSON.
To add a simulator in SAMSON we need to create and add
a dynamical model, an interaction model, and a state updater.
"""
if simulator == None: return
# get the dynamical model from the simulator
dynamicalModel = simulator.getDynamicalModel()
# get the interaction model from the simulator
interactionModel = simulator.getInteractionModel()
# get the state updater from the simulator
stateUpdater = simulator.getStateUpdater()
# make the operations in SAMSON undoable
SAMSON.beginHolding('Add simulator')
if SAMSON.isHolding():
# hold objects for undo/redo in SAMSON
SAMSON.hold(dynamicalModel)
SAMSON.hold(interactionModel)
SAMSON.hold(simulator)
# create nodes
dynamicalModel.create()
interactionModel.create()
simulator.create()
stateUpdater.create()
# add nodes to the active document in SAMSON
document = SAMSON.getActiveDocument()
document.addChild(dynamicalModel)
document.addChild(interactionModel)
document.addChild(simulator)
SAMSON.endHolding()
# show a window with interaction model properties
SAMSON.showProperties(interactionModel)
# initialize interaction model
interactionModel.initializeInteractions()
# show a window with state updator properties
SAMSON.showProperties(stateUpdater)
# Apply a simulator
# with Brenner interaction model and Interactive modeling state updater
# to a molecule in the active document
# get an indexer of all structural models in the active document
indexer = SAMSON.getNodes('n.t sm')
# make an instance of a simulator
simulator = SAMSON.makeSimulator(
nodeIndexer = indexer,
# the name of the interaction model class
interactionModelClassName = 'SMMBrennerInteractionModel',
# the UUID of the Brenner SAMSON Extension
interactionModelExtensionUUID = SBUUID('AD608CB6-6971-7CD4-6FCC-34531998E743'),
# the name of the state updater class
stateUpdaterClassName = 'SESInteractiveModelingUpdater',
# the UUID of the State Updater pack SAMSON Extension
stateUpdaterExtensionUUID = SBUUID('F912F119-7CBB-B5BD-972A-0A02DFCF683D')
)
if simulator:
# add the simulator in the active document in SAMSON
samson_add_simulator(simulator)
# make the operations in SAMSON undoable
SAMSON.beginHolding('Simulate with ASE')
# Run ASE calculations
# Create an object of the ExampleASECalculator class,
# which will initialize system in ASE and in SAMSON
ase_calc = ExampleASECalculator()
# store the initial conformation in the active document
ase_calc.store_current_conformation('Initial conformation')
# Run computations: force computation using SAMSON and minimization using ASE
ase_calc.compute()
# store the final conformation in the active document
ase_calc.store_current_conformation('Final conformation')
SAMSON.endHolding()
During these computations, you may see in the SAMSON’s viewport how atoms in the system are moving towards the minima.
In the console output you may get something like this:
Step Time Energy fmax
FIRE: 0 14:33:42 -35.994858 0.1356
FIRE: 1 14:33:42 -36.023361 0.1356
...
FIRE: 34 14:33:43 -36.881455 0.0077
FIRE: 35 14:33:43 -36.887776 0.0081
The script also creates the initial and final conformations of the system - you can switch between them by double-clicking on them in the Document view or in Python via SBConformation.restore().