# 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:

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

## 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):

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:

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