# Preparing for an MD calculation¶

[Input: recipes/moleculardynamics/initialstructre/]

The initial structure for starting a molecular dynamics simulation should usually be fully structurally relaxed. This is done to remove any excess potential energy, which would otherwise rapidly convert into kinetic energy (so heating up the system).

This example relaxes the geometry of a molecule (PTCDA) to a local miniumum, prior to performing molecular dynamics.

After running the relaxation, the output file `geo_end.gen`

is now suitable to
use as a starting structure for MD (or for other property calculations).

# Vibrational modes¶

[Input: recipes/moleculardynamics/vibrations/]

Once at a structural minimum, the quasi-harmonic vibrational modes of the atoms in the system can be calculated. These can then be compated with the power spectrum of the system, as determined by molecular dynamics.

The vibrational modes can be obtained from the mass-weighted Hessian matrix of
the system. This can be calculated by `DFTB+`

from finite difference second
derivatives of the energy with respect to atom positions (technically, the first
derivatives of the forces are actually used).

This derivative calculation can be enabled as an alternative choice of the Driver:

```
Driver = SecondDerivatives {
Delta = 1E-4
}
```

The provided example also sets the size of the finite difference step used to differentiate the energy as 0.0001 atomic units.

The modes are known as quasi-harmonic, since depending on the numerical step size these derivatives will include a contribution from higher even derivatives of the function. Hence, for accurate evaluation of harmonic frequencies, the smallest stable choice of the step size should be used (however, a carefull choice of this value can sometimes be used to include anharonicity in the calculated vibrational energies, see for example Jones and Briddon 1998).

Running the code will produce a file called `hessian.out`

which contains the
derivative information.

## Calculating the modes¶

The `modes`

code can read the Hessian matrix and calculate the vibrational
modes. An example input file is included in this directory using the
`hessian.out`

file that is produced by the DFTB+ calculation using the
supplied `dftb_in.hsd`

.

The general structure of the input for modes is similar to DFTB+, also being in
the hsd format. The code expects a `modes_in.hsd`

file, which contains blocks
for the geometry, Hessian and the Slater-Koster files (used to get the mass of
the atoms). A typical input looks like:

```
# Needs the equilibrium geometry, at which the Hessian had been calculated
Geometry = GenFormat {
<<< geom.gen
}
DisplayModes = {
PlotModes = -20:-1 # Take the top 10 modes
Animate = Yes # make xyz files showing the atoms moving
}
# You need to specify the SK-files, as the mass of the elements is needed
SlaterKosterFiles = Type2FileNames {
Prefix = "../../slakos/mio-ext/"
Separator = "-"
Suffix = ".skf"
}
# Include the Hessian, which was calculated by DFTB+
Hessian = {
<<< "hessian.out"
}
# This file uses the 3rd input format of the modes code
InputVersion = 3
```

Here, the code produces animations of several modes, which can be viewed for example using the jmol cross platform viewer

```
jmol mode_114.xyz &
```

for the highest frequency stretch mode of this molecule.