Using ASE in SAMSON

Using ASE in SAMSON

The Atomic Simulation Environment (ASE) is a set of tools and Python modules for setting up, manipulating, running, visualizing and analyzing atomistic simulations.

In this tutorial, we will show how to combine SAMSON and other molecular dynamics packages on the example of the ASE package.

To install ase package in conda environment:

conda install -c conda-forge ase

You may either copy-paste the examples below in your console inside SAMSON or create a py-files and import them as modules in the console. In the last case, let us change directory to some example directory and add this directory to the python system path:

# Change to the directory with the script
import os, sys
os.chdir('/path/to/examples')
sys.path.append('.')                     # add this path to the system path
# or, instead, just add the directory to the python system path without changing the directory
sys.path.append('/path/to/examples')

You can also find the code in SAMSON Python Scripting samples on github.

Content

Optimization of a water molecule using BFGS

In this example we will create a water molecule using ASE, copy it to SAMSON, do optimization in ASE, and update positions of the atoms in SAMSON accordingly.
For the sake of this example, we will run an optimization for several steps and we will be updating atoms’ positions in SAMSON according to their change in ASE after each step.

import numpy as np
from ase import Atoms
from ase.optimize import BFGS
from ase.calculators.emt import EMT
import samson as sam                                                            # import samson
 
class ExampleASEOptimization:
 
        def __init__(self):
                '''
                Initialize ASE and SAMSON
                '''
                self.initialize_system_ase()                                    # initialize ASE
                self.initialize_system_samson()                                 # initialize SAMSON
 
        def compute(self):
                '''
                Do computations using ASE and update positions in SAMSON
                '''
                dyn = BFGS(self.molecule)                                       # set up computations in ASE
                nsteps = 10
                for i in range(nsteps):
                        dyn.run(fmax = 0.005, steps = 1)                        # run the computations in ASE
                        self.update_positions_samson()                          # update positions in SAMSON
 
        def initialize_system_ase(self):
                '''
                Create an H2O molecule using ASE package; set a calculator
                '''
                d = 0.9575                                                      # set parameters
                t = np.pi / 180 * 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 to the active document in SAMSON
                '''
                structural_model = sam.Modeling.StructuralModel.StructuralModel() # a new structural model
                structural_model.create()                                       # create this structural model in SAMSON
                structural_model.name = 'H20 opt'
                document = SAMSON.getActiveDocument()                           # get the active document
                document.addChild(structural_model)                             # add the structural model to the active document
 
                self.add_atoms_in_samson(structural_model, self.molecule)       # add atoms in SAMSON
 
                self.indexer = structural_model.getNodes('n.t a')               # get a list of all atoms in the structural model
 
                SAMSON.getActiveCamera().topView()                              # set a camera view in SAMSON
 
        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
                        sa = sam.Modeling.StructuralModel.Atom(SAMSON.getElementTypeBySymbol(a.symbol),
                                sam.DataModel.Quantity.angstrom(a.position[0]),
                                sam.DataModel.Quantity.angstrom(a.position[1]),
                                sam.DataModel.Quantity.angstrom(a.position[2])) # construct an atom in SAMSON with a given element type and position
                        sa.create()                                             # create an atom in SAMSON
                        structural_model.addChild( sa )                         # add an atom to the structural model
 
        def update_positions_samson(self):
                '''
                Update positions of atoms in SAMSON based on new positions of these atoms in ASE
                '''
                new_positions = self.molecule.get_positions()                   # get the new positions from ASE
 
                for i in range(self.indexer.size):
                        self.indexer[i].setX(sam.DataModel.Quantity.angstrom(new_positions[i, 0]))
                        self.indexer[i].setY(sam.DataModel.Quantity.angstrom(new_positions[i, 1]))
                        self.indexer[i].setZ(sam.DataModel.Quantity.angstrom(new_positions[i, 2]))
 
                SAMSON.processEvents()                                          # process all events in SAMSON (to update positions in the viewport)

You may save it as a file, let us call it ‘exaseopt.py’. Now, to use it in SAMSON, we need to import this module:

import exaseopt

Now, we can create an object of the ExampleASEOptimization class, which will initialize ASE and SAMSON

ex1 = exaseopt.ExampleASEOptimization()

You may see a water molecule in the SAMSON’s viewport.

Let us run the computations now:

ex1.compute()

In the console output you may get something like this:

      Step     Time          Energy         fmax
BFGS:    0 14:51:51        2.733104        8.4006
BFGS:    1 14:51:51        2.136563        4.2991
BFGS:    2 14:51:51        1.901701        1.1757
BFGS:    3 14:51:51        1.880064        0.1810
BFGS:    4 14:51:51        1.879461        0.0212
BFGS:    5 14:51:51        1.879452        0.0132
BFGS:    6 14:51:51        1.879436        0.0224
BFGS:    7 14:51:51        1.879399        0.0420
BFGS:    8 14:51:51        1.879322        0.0647
BFGS:    9 14:51:51        1.879193        0.0800

Nitrogen on copper

Let us for the tutorial reasons adopt a tutorial from the ASE: Nitrogen on copper. For that we will slightly modify our previous example.

import numpy as np
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
import samson as sam                                                            # import samson
 
class ExampleASEOptimization2:
 
        def __init__(self):
                '''
                Initialize ASE and SAMSON
                '''
                self.initialize_system_ase()                                    # initialize ASE
                self.initialize_system_samson()                                 # initialize SAMSON
 
        def compute_adsorption_energy(self):
                '''
                Do computations using ASE and 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, trajectory = 'N2Cu.traj')
                dyn.run(fmax=0.05)
 
                self.update_positions_samson()                                  # update positions in SAMSON
 
                print('Adsorption energy:', self.e_slab + self.e_N2 - self.slab.get_potential_energy())
 
        def simulate(self):
                '''
                Do some MD computations using ASE and update positions in SAMSON
                '''
                # set up computations in ASE
                dyn = VelocityVerlet(self.molecule, dt = 1.0 * units.fs)
                for i in range(10):
                        pot = self.molecule.get_potential_energy()
                        kin = self.molecule.get_kinetic_energy()
                        print('%2d: %.5f eV, %.5f eV, %.5f eV' % (i, pot + kin, pot, kin))
                        dyn.run(steps=20)                                       # run the computations in ASE
                        self.update_positions_samson()                          # update positions in SAMSON
 
        def initialize_system_ase(self):
                '''
                Create an H2O molecule using ASE package; set a calculator
                '''
                d = 1.10                                                        # set parameters
                # 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 to the active document in SAMSON
                '''
                structural_model = sam.Modeling.StructuralModel.StructuralModel() # a new structural model
                structural_model.create()                                       # create this structural model in SAMSON
                structural_model.name = 'ase slab opt'
                document = SAMSON.getActiveDocument()                           # get the active document
                document.addChild(structural_model)                             # add the structural model to the active document
 
                self.add_atoms_in_samson(structural_model, self.slab)           # add slab atoms in SAMSON
                self.add_atoms_in_samson(structural_model, self.molecule)       # add molecule atoms in SAMSON
 
                self.indexer = structural_model.getNodes('n.t a and N')         # get a list of all Nitrogen atoms in the structural model
 
                SAMSON.getActiveCamera().leftView()                             # set a camera view in SAMSON
 
        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
                        sa = sam.Modeling.StructuralModel.Atom(SAMSON.getElementTypeBySymbol(a.symbol),
                                sam.DataModel.Quantity.angstrom(a.position[0]),
                                sam.DataModel.Quantity.angstrom(a.position[1]),
                                sam.DataModel.Quantity.angstrom(a.position[2])) # construct an atom in SAMSON with a given element type and position
                        sa.create()                                             # create an atom in SAMSON
                        structural_model.addChild( sa )                         # add an atom to the structural model
 
        def update_positions_samson(self):
                '''
                Update positions of atoms in SAMSON based on new positions of these atoms in ASE
                '''
                new_positions = self.molecule.get_positions()                   # get the new positions from ASE
 
                for i in range(self.indexer.size):
                        self.indexer[i].setX(sam.DataModel.Quantity.angstrom(new_positions[i, 0]))
                        self.indexer[i].setY(sam.DataModel.Quantity.angstrom(new_positions[i, 1]))
                        self.indexer[i].setZ(sam.DataModel.Quantity.angstrom(new_positions[i, 2]))
 
                SAMSON.processEvents()                                          # process all events in SAMSON (to update positions in the viewport)

You may save it as a file, let us call it ‘exaseopt2.py’. Now, to use it in SAMSON, we need to import this module:

import exaseopt2

Now, we can create an object of the ExampleASEOptimization2 class, which will initialize ASE and SAMSON

ex2 = exaseopt2.ExampleASEOptimization2()

You may see Nitrogen and Copper molecules in the SAMSON’s viewport.

Let us compute the adsorption energy:

ex2.compute_adsorption_energy()
                Step[ FC]     Time          Energy          fmax
*Force-consistent energies used in optimization.
BFGSLineSearch:    0[  0] 16:49:46       11.689927*       1.0797
BFGSLineSearch:    1[  2] 16:49:46       11.670814*       0.4090
BFGSLineSearch:    2[  4] 16:49:47       11.625880*       0.0409
Adsorption energy: 0.3235194223181601

Let us perform MD simulation of Nitrogen atoms:

ex2.simulate()

During this short simulation, you may see in the SAMSON’s viewport how these Nitrogen atoms move atop of the Copper slab.

If you want to visualize the saved trajectory using ASE:

from ase.io import read
from ase.visualize import view
slab_mol = read('N2Cu.traj')
view(slab_mol)

Usage of the SAMSON’s Interaction Model

In this example, we will be using an Interaction Model provided by SAMSON together with ASE for minimization of the system. This 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).

from ase import Atoms, Atom
from ase.optimize import FIRE
from ase.calculators.calculator import Calculator, all_changes
import numpy as np
import samson as sam                                                    # import samson
 
class IMWrap(Calculator):
        ''' Calculator for ASE, see: https://wiki.fysik.dtu.dk/ase/development/calculators.html '''
 
        implemented_properties = ['energy', 'forces']                   # properties that calculator can handle 
 
        im = None
        molecule = None
        indexer = None
 
        def __init__(self, **kwargs):
                Calculator.__init__(self, **kwargs)
                if kwargs is not None:
                        print("set IMPS")
                        self.im = kwargs["im"]
                        self.molecule = kwargs["molecule"]
                        self.indexer = kwargs["indexer"]
                else:
                        print("Error: input parameters are not provided")
                        pass
 
        def calculate(self, atoms = None, properties = ['energy', 'forces'], system_changes = all_changes):
                ''' The main function for Calculator that does the calculations '''
 
                self.im.updateInteractions()                            # update interactions in SAMSON
 
                forces = np.zeros((len(self.molecule), 3))
                for i in range(self.indexer.size):                      # get forces from the SAMSON's InteractionModelParticleSystem
                        forces[i, 0] = self.im.getForce(i)[0].N
                        forces[i, 1] = self.im.getForce(i)[1].N
                        forces[i, 2] = self.im.getForce(i)[2].N
 
                force_conversion = 1. / 8.238722514e-8                  # 1 atomic unit force = Eh / a0 = 8.2387225(14) x 10^(-8) N
                # set calculated properties for ASE
                self.results['forces'] = forces * force_conversion      # convert forces from N to atomic unit force
                self.results['energy'] = self.im.energy.eV.value        # get energy from the SAMSON's InteractionModelParticleSystem
 
                self.update_positions_samson()                          # update positions of atoms in SAMSON
 
        def update_positions_samson(self):
                '''
                Update positions of atoms in SAMSON based on new positions of these atoms in ASE
                '''
                new_positions = self.molecule.get_positions()           # get the new positions from ASE
 
                for i in range(self.indexer.size):
                        self.indexer[i].setX(sam.DataModel.Quantity.angstrom(new_positions[i, 0]))
                        self.indexer[i].setY(sam.DataModel.Quantity.angstrom(new_positions[i, 1]))
                        self.indexer[i].setZ(sam.DataModel.Quantity.angstrom(new_positions[i, 2]))
 
                SAMSON.processEvents()                                  # process all events in SAMSON (to update positions in the viewport)
 
class ExampleASECalculator:
        '''
        An example of the usage of a SAMSON's InteractionModelParticleSystem together with ASE
        We assume, that the InteractionModelParticleSystem has been already created in SAMSON for the atoms in the SAMSON's data graph
        '''
 
        molecule = None
 
        def __init__(self):
 
                indexer = SAMSON.getNodes('n.t a')                      # get list of all atoms in SAMSON
 
                self.create_atoms_in_ase(indexer)                       # create atoms in ASE
 
                imps = SAMSON.getNodes('n.t imps')[0].castToInteractionModelParticleSystem()  # get a SAMSON's InteractionModelParticleSystem
 
                self.molecule.set_calculator(IMWrap(im = imps, molecule = self.molecule, indexer = indexer))  # set an ASE calculator
 
        def create_atoms_in_ase(self, indexer):
                '''create atoms in ASE'''
                self.molecule = Atoms()                                 # create an empty molecule in ASE
 
                for a in indexer:                                       # add atoms in the molecule
                        self.molecule.append( Atom( a.elementSymbol, (a.getX().angstrom, a.getY().angstrom, a.getZ().angstrom) ) )
 
        def compute(self):
                dyn = FIRE(self.molecule)                               # set the computation
                dyn.run(fmax = 0.005, steps = 100)                      # run the computation

Let us open some hydrocarbon from the Samples folder (e.g. Cyclopropane.pdb) and apply a simulator to it with the “Interaction Model” set to “Brenner” and “State Updater” to “Interactive modeling”. Now you may slightly move some atoms in the molecule using SAMSON Editors (e.g. choose the Point selection and move slightly some atoms).

You may save it as a file, let us call it ‘exaseim.py’. Now, to use it in SAMSON, we need to import this module:

import exaseim

Now, we can create an object of the ExampleASECalculator class, which will initialize ASE and SAMSON

ex = exaseim.ExampleASECalculator()

Now we can run computations (force computation using SAMSON and minimization using ASE):

ex.compute()

During these computations, you may see in the SAMSON’s viewport how atoms in the system are moving towards the minima.

Python Script, tutorial: ASE + Interaction Model

Comments are closed.