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.
Parameter |
Default Value |
Description |
|---|---|---|
ICHAIN |
1 |
Use the Dynamical Matrix method |
IBRION |
3 |
Specify that VASP do MD with a zero time step |
POTIM |
0 |
Zero time step so that VASP does not move the ions |
EDIFF |
1E-8 |
Strict force criteria to accurately find force differences (curvature) |
EDIFFG |
-1E-8 |
Make sure that vasp does not quit if the forces are low |
NSW |
DOF+1 |
Do a force call at equilibrium and along each degree of freedom |
Compile a version of VASP with the vtst source.
Generate a DISPLACECAR file from either:
- dymselsph.pl (POSCAR) (Central atom) (radius) (displacement)
Selects atoms within a radius around the central atom(s) when calculating harmonic frequencies of a molecule.
or
- dymseldsp.pl (POSCAR 1) (POSCAR 2) (atoms to include) (displacement)
Selects the atoms with most displacement between two POSCAR files.
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.
Skip this step if you choose to run the calculation with a single image where a single displacement is calculated at a time.
If you want to run in parallel (multiple images where many displacements are calculated simultaneously), you need to choose a number of images, M, which evenly divides the number of degrees of freedom (DOF) you displaced. Similarly to how a nudged elastic band calculation is set up, create M subdirectories numbered 01 to M, using two digits for each directory name. Place the POSCAR file in each subdirectories.
Set up the INCAR file. You need to have the following variables set:
ICHAIN=1
POTIM=0.0
IMAGES=M (only if running in multiple geometries at once mode from step 2)
NSW=(DOF/M)+1
Note about setting the number of ionic steps, NSW The reason for adding 1 to the NSW (number of ionic steps) is that each calculation first computes the forces on the initial structure before any displacements occur. In a single image calculation, the NSW should be one more than the total number of degrees of freedom. In a parallel calculation, NSW is calculated as one more than the number of degrees of freedom divided by the number of images, M. For example, in a single image calculation for the vibrational modes of ethane, with 3 degrees of freedom per atom for 8 atoms, the NSW would be 25. For a parallel calculation with M=4 images, set IMAGES=4 and NSW would be calculated as NSW=(24/4)+1=7.
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.
Extract the force constants from the output using the dymmatrix.pl script specified in one of the following formats:
dymmatrix.pl (DISPLACECAR) (OUTCAR)
dymmatrix.pl (#DISPLACECAR) (DISPLACECAR1) (DISPLACECAR2) …
dymmatrix.pl (OUTCAR1) (OUTCAR2) (OUTCAR3) …
The output of this script is the dynamical matrix, its eigenvalues, the normal mode frequencies and the normal mode eigenvectors. Some examples are:
- Default, single calculation and single image:
dymmatrix.pl or dymmatrix.pl DISPLACECAR OUTCAR
- Single calculation with multiple images:
dymmatrix.pl 1 DISPLACECAR ??/OUTCAR
- Multiple calculations with n DISPLACECAR and OUTCAR files:
dymmatrix.pl n DISPLACECAR1 .. DISPLACECAR[n] OUTCAR1 .. OUTCAR[n]
If all calculated frequencies are 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.
Check for convergence of frequencies with respect to all parameters, and particularly the finite difference magnitude in the DISPLACECAR file, the accuracy of the forces (EDIFF) and the magnitude of the smearing. A small smearing value is particularly important for calculating the normal modes of molecules with large homo-lumo gaps.
Calculate the prefactor (not for simple frequency analysis) with dymprefactor.pl.
dymprefactor.pl (freq.dat of minimum) (freq.dat of transition state)
The input is the eigenvalue files of both the minimum and the transition state. The output is the prefactor in inverse centimeters.
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:
Save the OUTCAR files from the previous calculation.
Create the new DISPLACECAR as above.
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.
Modify NSW in your INCAR file as before.
Run Vasp.
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.