7.6. Si : Anharmonic interatomic force constants (IFCs)

This page explains another method to calculate anharmonic interatomic force constants (IFCs) using ALAMODE.

The example input files are provided in example/Si/anharm_IFCs.

Let us move to the example directory.

$ cd ${ALAMODE_ROOT}/example/Si/anharm_IFCs

1. Calculate the harmonic IFCs

First, we calculate the harmonic IFCs in the considering supercell.

The procedure is the same as explained in Tutorial 7.1, so we briefly describe the outline here.

We first calculate the displacement patterns by alm.

$ cd ${ALAMODE_ROOT}/example/Si/anharm_IFCs/1_harmonic
$ ${ALAMODE_ROOT}/alm/alm si_alm_sug.in > si_alm_sug.log

Then, generate the VASP input files using the generated si222_harmonic.pattern_HARMONIC and the input files in the VASP_input directory. After running the VASP calculation and obtaining DFSET_harmonic (You can skip this part because DFSET_harmonic is provided in this tutorial), calculate the harmonic IFCs with

$ ${ALAMODE_ROOT}/alm/alm si_alm_opt.in > si_alm_opt.log

2. Generate the displacement-force data

In Tutorial 7.5, we used the AIMD calculation to generate random configurations, from which we calculated the anharmonic IFCs.

However, the low-energy region of the potential energy surface (PES) can be sampled more efficiently for weakly anharmonic materials or materials without imaginary harmonic frequencies. Here, we use the harmonic IFCs and perform the independent random sampling from the harmonic PES at a given temperature.

In preparation, we calculate the phonon frequencies and the polarization vectors at commensurate \(q\)-points.

$ cd ${ALAMODE_ROOT}/example/Si/anharm_IFCs/2_generate_config
$ python3 ${ALAMODE_ROOT}/tools/displace.py --VASP POSCAR_supercell --prim POSCAR_primitive_cell --random_normalcoord >> "si_anphon.in"

The information on the \(q\)-points are written out in si_anphon.in as follows.

&general
  PREFIX = si222_harmonic
  MODE   = phonons
  FCSXML = si222_harmonic.xml

  NKD = 1; KD = Si
/

*****************************************************************
    displace.py --  Generator of displaced configurations
                      Version. 1.2.1
*****************************************************************

Output format                  : VASP POSCAR
Structure before displacements : POSCAR_supercell
Output file names              : disp{counter}.POSCAR
Magnitude of displacements     : 0.02 Angstrom
Number of atoms                : 64

The --evec option is necessary when '--random_normalcoord'
option is used.
Please generate a PREFIX.evec file by using the ANPHON code
with the following inputs and then run displace.py again with
--evec=PREFIX.evec option:

&cell
1.0
  0.000000000000000   5.131551292420093   5.131551292420093
  5.131551292420093   0.000000000000000   5.131551292420093
  5.131551292420093   5.131551292420093   0.000000000000000
/
&kpoint
0
  0.000000000000000    0.000000000000000    0.000000000000000
  ...

Now, delete the unnecessary part of the output and run the anphon calculation.

$ ${ALAMODE_ROOT}/anphon/anphon si_anphon.in > si_anphon.log

The calculated phonon frequencies and the polarization vectors are stored in si222_harmonic.evec.

With these preparations, we can generate supercells with random displacements by

$ mkdir configurations
$ cd configurations
$ cp ../POSCAR_primitive_cell  ../POSCAR_supercell ../si222_harmonic.evec ./
$ python3 ${ALAMODE_ROOT}/tools/displace.py --VASP POSCAR_supercell --prim POSCAR_primitive_cell --random_normalcoord --evec si222_harmonic.evec --temp 300 --prefix randomQ_ -nd 100

Here, we generated -nd 100 configurations by randomly sampling from the distribution at -temp 300 K in the harmonic PES.

Please run the DFT calculation for each generated supercell using the VASP input in example/Si/anharm_IFCs/1_harmonic/VASP_input. Then, use extract.py to obtain DFSET_randomQ using the procedure explained in Tutorial 7.1.

Note

The imaginary frequencies are replaced by their absolute values in the random sampling. Thus, the procedure can be performed for the strongly anharmonic materials as well.

However, the user must be careful whether the generated set of random configurations is a good dataset for calculating IFCs.

3. Cross validation (CV)

In this step, we explain how to run the different sets of CV calculations separately.

The calculation of different sets can be executed in parallel because they are independent of each other. So, if you have a cluster computer with multiple cores, you can run the calculations of each CV set in separate jobs. The preparation of the input files is slightly complicated, but it will be time-saving when the computational cost of the CV calculation is significant.

The input files are si_alm_cvset1.in to si_alm_cvset4.in in example/Si/anharm_IFCs/3_cv.

The essential parts of the input file si_alm_cvset1.in are as follows.

We have NDATA = 100 displacement-force data, and we will perform CV with 4 sets. Thus, we want to use the first 25 data (NSTART_CV = 1, NEND_CV = 25) in the validation process in the calculation of the first CV set (set1). Note that these 25 sets have to be excluded in the training process (SKIP = 1-25)

The input files of the other CV sets are set accordingly. It is important that we use different PREFIX for each set because the result of another CV set will overwrite the output file otherwise.

&general
  PREFIX = si222_cvset1
  ...
/

...
&optimize
  ...
  NDATA = 100
...
  SKIP = 1-25
  NSTART_CV = 1
  NEND_CV = 25
/

...

Run the calculation with

$ cd ${ALAMODE_ROOT}/example/Si/anharm_IFCs/3_cv
$ ${ALAMODE_ROOT}/alm/alm si_alm_cvset1.in > si_alm_cvset1.log
$ ${ALAMODE_ROOT}/alm/alm si_alm_cvset2.in > si_alm_cvset2.log
$ ${ALAMODE_ROOT}/alm/alm si_alm_cvset3.in > si_alm_cvset3.log
$ ${ALAMODE_ROOT}/alm/alm si_alm_cvset4.in > si_alm_cvset4.log

After all the calculations are finished, collect the cvscore data with

$ python3 cvscore.py *cvset > si222.cvscore

Note

The number of \(\alpha\) for which the calculation is performed can differ depending on the CV sets because the calculation stops in the middle due to the STOP_CRITERION-tag. If the calculations stop at different steps, the Python script stops with an error of “Inconsistent number of entries”.

In that case, please manually adjust the cvset files so that the number of entries is consistent.

The optimal amplitude of regularization (\(\alpha\)) can be read from the last line of si222.cvscore.

#Minimum cvscore at  2.25633e-07

4. Calculation of IFCs

Finally, we calculate the IFCs of silicon in example/Si/anharm_IFCs/4_optimize.

The input file is si_alm_opt.in. Set CV = 0 and set the optimal \(\alpha\) with L1_ALPHA = 2.25633e-07 in &optimize-field.

Run the calculation with

$ cd ${ALAMODE_ROOT}/example/Si/anharm_IFCs/4_optimize
$ ${ALAMODE_ROOT}/alm/alm si_alm_opt.in > si_alm_opt.log

The calculated IFCs are written out in si222.xml and si222.fcs. The fitting error is

RESIDUAL (%): 0.524303