Applications

Once a sGDML model is trained (see force field reconstruction), it can be used like any other force field. Here, we show a few typical applications. Additional practical examples can be found in Chmiela et al. 2019 [1].

Warning

The examples below rely on tools from the Atomic Simulation Environment (ASE), which expects eV and Angstrom as energy and length units, respectively. However, sGDML models are unit agnostic, which is why a unit conversion might be necessary. By default, the included SGDMLCalculator ASE interface assumes that kcal/mol and Angstrom are used in the model and performs the appropriate conversion. Please specify your own conversion factors during initialization of the calculator, if this is not the case.

[1]Chmiela, S., Sauceda, H. E., Poltavsky, Igor, Müller, K.-R., Tkatchenko, A. (2019). sGDML: Constructing Accurate and Data Efficient Molecular Force Fields Using Machine Learning. ‎Comput. Phys. Commun., 240, 38-45.

MD simulations [ASE]

The following example (based on this tutorial) shows how to run a constant-energy molecular dynamics (MD) simulation of ethanol at a temperature of 300 K and a resolution of 0.2 fs for 10 ps. A pre-trained ethanol model can be downloaded as described here. As starting point, a geometry file ethanol.xyz in XYZ format is also needed (download).

from sgdml.intf.ase_calc import SGDMLCalculator

from ase.io import read
from ase.optimize import QuasiNewton
from ase.md.velocitydistribution import (MaxwellBoltzmannDistribution, Stationary, ZeroRotation)
from ase.md.verlet import VelocityVerlet
from ase import units

model_path = 'm_ethanol.npz'
calc = SGDMLCalculator(model_path)

mol = read('ethanol.xyz')
mol.set_calculator(calc)

# do a quick geometry relaxation
qn = QuasiNewton(mol)
qn.run(1e-4, 100)

# set the momenta corresponding to T=300K
MaxwellBoltzmannDistribution(mol, 300 * units.kB)
Stationary(mol) # zero linear momentum
ZeroRotation(mol) # zero angular momentum

# run MD with constant energy using the VelocityVerlet algorithm
dyn = VelocityVerlet(mol, 0.2 * units.fs, trajectory='md.traj')  # 0.2 fs time step.

def printenergy(a):
        # function to print the potential, kinetic and total energy
        epot = a.get_potential_energy() / len(a)
        ekin = a.get_kinetic_energy() / len(a)
        print('Energy per atom: Epot = %.3feV  Ekin = %.3feV (T=%3.0fK)  '
                'Etot = %.3feV' % (epot, ekin, ekin / (1.5 * units.kB), epot + ekin))

# now run the dynamics
printenergy(mol)
for i in range(500):
        dyn.run(10)
        printenergy(mol)

Visualizing results

The code above first relaxes the given geometry in the provided force field model, before running an MD simulation. It will then output a trajectory file md.traj that can be visualized by running

$ ase gui md.traj

The result will look something like this:

_images/ase_md.png

You can easily export each step of this trajectory as individual Extended XYZ files using the following code snippet:

import sys
import os.path
from ase.io.trajectory import Trajectory
from ase.io import write

our_dir = 'out'
if os.path.isdir(our_dir):
        sys.exit('Output directory already exists.')
else:
        os.mkdir(our_dir)

traj = Trajectory('md.traj')
for i,atoms in enumerate(traj):

atoms.wrap()
write(our_dir + '/step'+("{:05d}".format(i))+'.xyz', atoms, append=False)

Vibrational modes [ASE]

The following example (based the example here) shows how to obtain the vibrational modes for ethanol. Once again, a pre-trained ethanol model (here) and a geometry file ethanol.xyz in XYZ format are needed (download).

from sgdml.intf.ase_calc import SGDMLCalculator

from ase.io import read
from ase.optimize import QuasiNewton
from ase.vibrations import Vibrations

model_path = 'ethanol.npz'
calc = SGDMLCalculator(model_path)

mol = read('ethanol.xyz')
mol.set_calculator(calc)

# do a quick geometry relaxation
qn = QuasiNewton(mol)
qn.run(1e-4, 100)

# run the vibration calculations
vib = Vibrations(mol)
vib.run()

vib.summary() # print a summary of the vibrational frequencies
vib.write_jmol() # write file for viewing of the modes with jmol

vib.clean() # remove pickle-files

The result will look something like this:

---------------------
  #    meV     cm^-1
---------------------
  0    1.4i     11.1i
  1    0.5i      4.1i
  2    0.1i      0.9i
  3    0.0i      0.0i
  4    0.0       0.4
  5    0.1       0.9
  6   27.6     222.8
  7   33.8     273.0
  8   50.2     405.0
  9   97.9     789.9
 10  108.0     870.7
 11  124.0    1000.3
 12  133.1    1073.9
 13  141.0    1137.5
 14  153.4    1237.6
 15  155.8    1256.7
 16  166.6    1344.0
 17  173.5    1399.0
 18  177.7    1433.0
 19  180.2    1453.8
 20  181.8    1466.6
 21  358.9    2894.5
 22  364.9    2943.4
 23  368.4    2971.1
 24  377.9    3047.8
 25  379.0    3056.6
 26  461.2    3719.5
---------------------
Zero-point energy: 2.108 eV

Imaginary frequencies indicate that the geometry is not yet perfectly relaxed.

Visualizing results

The code above also outputs a file vib.xyz that can be used to visualize the modes, e.g. using Jmol:

_images/vib.gif

This is vibrational mode #8 from the table above.

Minima hopping [ASE]

This example (based on this tutorial) shows how to search for a globally optimal geometric configuration using the minima hopping algorithm. Once again, a pre-trained ethanol model and a geometry file ethanol.xyz in XYZ format are needed (download):

from sgdml.intf.ase_calc import SGDMLCalculator

from ase.io import read
from ase.optimize import QuasiNewton
from ase.optimize.minimahopping import (MinimaHopping, MHPlot)

model_path = 'ethanol.npz'
calc = SGDMLCalculator(model_path)

mol = read('ethanol.xyz')
mol.set_calculator(calc)

# do a quick geometry relaxation
qn = QuasiNewton(mol)
qn.run(1e-4, 100)

# instantiate and run the minima hopping algorithm
hop = MinimaHopping(mol, T0=4000., fmax=1e-4)
hop(totalsteps=5)

plot = MHPlot()
plot.save_figure('ase_mh.png')

Visualizing results

The script above will run five molecular dynamics trajectories, each generating its own output file. It will also produce a file called minima.traj which contains all accepted minima. Those can be easily viewed with ASE, using:

$ ase gui minima.traj

We can also easily convert it to other formats, such as Extended XYZ:

$ ase convert minima.traj minima.xyz

This command renders the first found minimum as PNG file:

$ ase convert -n 0 minima.traj minimum_0.png

Here are renderings of the structures found in each of the five runs (some minimas are visited multiple times):

_images/minimas.png

Additionally, the minima hopping process is visualized in a summary figure ase_mh.png. As the search is inherently random, your results will look different. In this example, the two distinct minima of ethanol have been successfully found:

_images/ase_mh.png

This type of plot is described here:

“In this figure, you will see on the \(E_{\text{pot}}\) axes the energy levels of the conformers found. The flat bars represent the energy at the end of each local optimization step. The checkmark indicates the local minimum was accepted; red arrows indicate it was rejected for the three possible reasons. The black path between steps is the potential energy during the MD portion of the step; the dashed line is the local optimization on termination of the MD step. Note the y axis is broken to allow different energy scales between the local minima and the space explored in the MD simulations. The \(T\) and \(E_{\text{diff}}\) plots show the values of the self-adjusting parameters as the algorithm progresses.”