Skip to content

Vibrational analysis from beginning to end

Timing: \(\sim\) 45 minutes total

Here, we will investigate the vibrational modes of the zundel cation. For this, we are going to use the built-in phonon calculator of the i-PI code.1

Steps of the procedure:

  1. geometry optimization of the zundel cation using purely FHIaims.
  2. building the i-PI input.xml file
  3. i-PI + FHI-aims vibrational analysis

1. Geometry optimization

We continue using LIGHT settings for the O and H atoms.

Convergence of BFGS trm

In general, using a convergence tolerance of 1E-3 eV/Ang is enough to converge most of the systems at hand, however, it is important to note that there exist special cases, complex structures, where occasionally, a much tighter (1E-4 eV/Ang) criterium is required in order to move away from saddle point structures and obtain correct vibrational frequencies.

For geometry optimization, please add the following keywords

relax_geometry trm 1E-3

With these settings we choose the BFGS optimization2 algorithm where the trust region method is applied to determine the optimization step.3

2. Building the i-PI input file

We take the optimized structure from the FHIaims (geometry.in.next_step) and create the appropriate .xyz file that i-PI will use as a starting geometry. The .xyz input file for i-PI should read

7
# positions{angstrom}
O    x  y  z
H    x  y  z
H    x  y  z
H    x  y  z
O    x  y  z
H    x  y  z
H    x  y  z

where you have to substitute the x, y and z values by the relaxed positions obtained in FHI-aims.

Reminder:

The important bit, as we already discussed, is that the order of the atomic species in this file is the same as in geometry.in.

Be aware that as a convention i-PI is using atomic units in its input and output files unless specified otherwise. The block in the i-PI input file input.xml that controls the options to perform vibrations is the following

<motion mode="vibrations">
  <vibrations mode='fd'>
    <pos_shift> 0.001 </pos_shift>
  </vibrations>
</motion>

3. Running the vibrational analysis

Please run i-PI and FHI-aims to calculate the vibrational eigenvalues by submitting

i-pi input.xml > log.ipi &

and

mpirun -np 4 aims.x > aims_vib.out &

4. Analyzing the output files

Take a look at the ouput files that i-PI has generated at the end of the run. You can see the collection of eigenvalues in simulation.phonons.eigval and the matrix of eigenvectors corresponding to the eigenvalues in simulation.phonons.eigvec.

What can you say about the stationary point judging by the number of zero eigenvalues? How many of these do you expect to see for an N-atom system?

Answer

For a non-linear N-atom system we can expect to see 6 eigenvalues very close to zero. They correspond to translational and rotational molecular motions. For a linear system the number of zero eigenvalues is 5. Therefore, one expects to see 3N-6 (3N-5) nonzero eigenvalues at the minimum geometry of the system for non-linear and linear molecules, respectively.

What would the presence of a negative eigenvalue mean? What can you suggest to an inexperienced user to do in such a situation?

Answer

A negative eigenvalue means that along that direction the potential energy is not a minimum but rather the potential energy is a concave function along the corresponding direction. This indicates that the system has not been properly minimized and is stuck at a saddle point of the 3N-dimensional potential energy surface. Although we expect relaxation to end in the stationary point, negative eigenvalues may also indicate that relaxation was aborted, or accuracy of relaxation was not enough for tiny finite displacements, or something else.

Possible suggestions for such a scenario could involve

1) reoptimizing the geometry with different convergence settings

2) using different grid settings (LIGHT/INTERMEDIATE/TIGHT) to enhance the precision of the DFT energy calculations

3) breaking the symmetry of the system by randomly changing some of the atomic positions before reoptimization of the geometry. This may help in finding a global, rather than local minimum structure, too.

Calculate the normal mode frequencies of the system in cm-1 units. In what spectral region do these vibrations give signal?

Answer

The normal mode frequencies for each vibrational motion can be calculated using the following formula. $$ \omega_i = \sqrt{\lambda_i}*219474.63,$$ where \(\lambda_i\) are the eigenvalues which are given in atomic units (atomic-time-2) in i-PI.

Visualize the motion along the first 3 vibrational normal modes of the system using a graphical visualizer of your choice. Note that in order to visualize the modes, one should be in Cartesian space, while the eigenvectors are in mass-weighted space. The vibrational modes to be visualized are stored in the simulation.phonons.mode file.

176 cm-1 380 cm-1 420 cm-1

Moreover, we receive as an output the dynamical matrix, which is the matrix that is further diagonalized by i-PI for one to obtain the vibrational modes. Check your notes and write out the relation between the Hessian and the dynmical matrix.

Answer

The elemens of the dynamical matrix relate to the Hessian matrix by mass-weighting that is achieved by: $$ K_{ij} = \dfrac{1}{\sqrt{m_i m_j}} H_{ij}$$


  1. There are many ways one can calculate vibrations for isolated and periodic systems with FHI-aims, this is by no means the only one. We recommend that the interested user has a look at the respective section in the FHI-aims manual and also looks at FHI-vibes

  2. Numerical Optimization, Jorge Nocedal and Stephen J. Wright, Springer, 2006. 

  3. It is also possible to perform the geometry relaxation in i-PI itself. That is triggered by the motion mode minimize. The same trm algorithm as implemented in FHI-aims is available in i-PI.