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].

[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 a starting point, a geometry file ethanol.xyz in XYZ format is also needed (download).

from sgdml.intf.ase_calc import SGDMLCalculator

from ase.io.xyz import read_xyz
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 = 'ethanol.npz'
calc = SGDMLCalculator(model_path)

mol = read_xyz('ethanol.xyz').next()
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

This 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

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.xyz import read_xyz
from ase.optimize import QuasiNewton
from ase.vibrations import Vibrations

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

mol = read_xyz('ethanol.xyz').next()
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 perfectly relaxed yet.

Visualizing results

The code above also outputs a file vib.xyz for viewing of 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.xyz import read_xyz
from ase.optimize import QuasiNewton
from ase.optimize.minimahopping import (MinimaHopping, MHPlot)

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

mol = read_xyz('ethanol.xyz').next()
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.”