Vasp TST Tools
     
DYNAMICAL MATRIX

The dynamical matrix code in VASP allows for the calculation of harmonic frequencies of and the prefactor of a reaction. The formalism was developed by Vineyard, and the Arrhenius rate includes the Vineyard prefactor. The basic idea is that the potential is assumed to be harmonic for both the initial state and the transition state. With this assumption, one can formulate the migration prefactor (the product of the attempt frequency and the exponential of the migration entropy over k) as the ratio of the products of the frequencies of the normal modes of the initial state over the transition state. Another way of saying this is that the entropy of the initial and transition state is evaluated using a harmonic approximation.

To be able to do a dynamical matrix calculation for a prefactor, one needs to know beforehand both the initial state and the transition state. The nudged elastic band (NEB) method and the dimer method are two ways of locating transition states.

The dynamical matrix calculation may also be used to calculate the harmonic frequencies and modes of a given molecule whose geometry has already been optimized in VASP. This dynamical matrix calculation requires only beforehand knowledge of the equilibrium optimized geometry of the initial state.

In implementing the dynamical matrix calculation in VASP, several support programs have been written, both to help set up the calculation as well as process the output. A page describing these programs, most of which are Perl scripts, can be found here.


Setting up a dynamical matrix calculation

This is a flow chart of the steps that need to be done in order to do a dynamical matrix calculation.

    CALCULATING THE FIRST DYNAMICAL MATRIX

  1. Compile a version of VASP with the vtst source.
  2. if you want to run VASP, set up the directory structure in a multiple-geometries-at-once mode in the same way that the nudged elastic band is run, then create the number of directories you want to run, from 01 to NN, where NN is the number of directories you want. To simply calculate vibrational frequencies of a molecule, this step is not necessary and may be ignored.
  3. Select the degrees of freedom to calculate in a DISPLACECAR file. This file is a 3 by N list of displacements where N is the number of atoms. A typical displacement and entry in this list is on the order of 0.001. This file may also be generated with either the dymselsph.pl script (to select atoms in a radius about central atom(s) when calculating harmonic frequencies of a molecule) or the dymseldsp.pl script (to select the atoms that are displaced the most between two POSCARs).
  4. Set up the INCAR file. You need to have the following variables set: ICHAIN=1, POTIM=0.0, IMAGES=NN (only if running in multiple geometries at once mode from step 2), NSW=(Degrees of freedom)/N+1. The reason for the +1 is because each directory will calculate the forces on the initial structure before doing any displacements.

    For instance, in the calculation of ethane modes, the NSW tag would be 25 (3 degrees of freedom for each of 8 atoms, plus 1). For ethylene, the NSW tag would be 19. No IMAGES tag would be added unless running in parallel as specified in step 2.

    Parameter
    Value
    Description
    ICHAIN1Run the dynamical matrix code
    IMAGESNNNumber of parallel images to run, if desired as in step 2 above. Otherwise, do not add.
    NSW(DOF)/N+1Number of ionic steps
    IBRION3Tell VASP to run dynamics
    POTIM0.0with a timestep of zero (ie, do nothing)

    Note: Choosing the number of time steps. The dynamical matrix calculation may be set up run in parallel, both with multiple processors working on one degree of freedom, and with different degrees of freedom in different directories. This parallelization is identical to that of the nudged elastic band. If there are N degrees of freedom and directories in which they are being calculated, then each directory needs to do N/M time steps. However, each directory does one time step without a displacement, one time step on the initial structure, and so each directory needs to do one more time step. Thus, the total number of time steps needed is: N/M + 1. Also note that M = IMAGES of the INCAR file.

  5. Run VASP after verifying that POTCAR, POSCAR, INCAR, DISPLACECAR, and KPOINTS files all exist and are properly configured. To ensure accuracy of calculated modes, the geometry in POSCAR should have first been optimized in VASP. All INCAR settings such as ISMEAR, SIGMA, PREC, etc used in the geometry optimization should be present in the dynamical matrix INCAR as well.
  6. Extract the force constants from the output. You can use the dymmatrix.pl script to do this. Please see the documentation for this script for its arguments. The output of this script is the dynamical matrix, its eigenvalues, the normal mode frequencies and the normal mode eigenvectors.
  7. If all calculated frequencies are all zero, check the EDIFFG tag in dynamical matrix INCAR. AFTER geometry minimization when running the dynamical matrix code, this INCAR tag may be set arbitrarily small to prevent premature termination of the job (as if minimum structural accuracy had been obtained). The number of force calls made should be equal to the number specified in the NSW tag.
  8. Check for convergence of frequencies with respect to all parameters.
  9. Calculate the prefactor (not for simple frequency analysis). The script dymprefactor.pl does this. It takes as input the eigenvalue files of both the minimum and the transition state. The output is the prefactor in inverse cm.

CALCULATING ADDITIONAL DEGREES OF FREEDOM FOR THE SAME DISPLACEMENT

This completes one prefactor calculation. If you want to do another calculation with the same system, but including more degrees of freedom, then you need to:

  1. Save the OUTCAR files from the previous calculation.
  2. Create the new DISPLACECAR as above.
  3. Generate the difference between the new and old DISPLACECAR with the dymcmpdisp.pl program. You will run the dynamical matrix code with the output of this program.
  4. Modify NSW in your INCAR file as before.
  5. Run Vasp.
  6. Extract the forces from the OUTCAR files. Use dymmatrix.pl to do this. To do this, you will need the DISPLACECAR from the original calculation, the DISPLACECAR from this calculation, and all of the corresponding OUTCARs. Again, the second and subsequent DISPLACECARs should only contain the new degrees of freedom not calculated before.