Dipole-Dipole effect at second order

The detailed equations can be found in papers Phys. Rev. B 50, 13035(R) (1994) and Phys. Rev. B 55, 10355 (1997).

The code not-only provides the dipole-dipole effects on the dynamic matrices at second order, but also the derivative of the effect with respect to q-vector.

The code is developed quite recently, it is in pretty good shape and probably only requires minor variables renaming. One point of concern is that there is a Fortran version of the code developed to improve the speed of the calculation. It is written in Fortran for its easier handling of matrices than C/C++. But it should be discussed whether a C/C++ version is needed, and which is better.

class DipoleDipole(struct, zeu=None, epsilon=None, dataset=None, G_ran=5, Lambda=0.20, factor=1.0, tol=1.0E-6)

• struct: periodica compatible structure
• zeu: array of floats, Born effective charges
• epsilon: array of floats, dielectric constants
• dataset: string, path a yaml file containing data including zeu, epsilon, G_ran, Lambda and factor
• G_ran: integer, range of G vectors
• Lambda: float, $\lambda$ in the equation
• factor: float, unit conversion factor (for VASP the factor is ~14.40)
• tol: float, error tolerance.

There are 2 ways to instantiate the object:

dp = DipoleDipole(
    struct="path to struct file", # or a periodica object
    dataset="path to epsilon.yml",
)
dp = DipoleDipole(
    struct="path to struct file",
    zeu=zeu,
    epsilon=epsilon,
    factor=factor,
)

import numpy as np
np.set_printoptions(linewidth=180, suppress=True)
from dipole_dipole import DipoleDipole

poscar = """\
PuO2 Fluorite
1.0
0.00000000   2.76451179   2.76451179
2.76451179   0.00000000   2.76451179
2.76451179   2.76451179   0.00000000
1 2
d
0.00000000   0.00000000   0.00000000  Th:p
0.25000000   0.25000000   0.25000000  O:p
0.75000000   0.75000000   0.75000000  O:p"""
factor = 14.40
epsilon = [
    [     4.883362,     0.000000,     0.000000],
    [     0.000000,     4.883362,     0.000000],
    [     0.000000,     0.000000,     4.883362],
]
zeu = [
    [[     5.407510,     0.000000,     0.000000],
     [     0.000000,     5.407790,     0.000000],
     [     0.000000,     0.000000,     5.409020],],
    [[    -2.705480,     0.000000,     0.000000],
     [     0.000000,    -2.705730,     0.000000],
     [     0.000000,     0.000000,    -2.705850],],
    [[    -2.705480,     0.000000,     0.000000],
     [     0.000000,    -2.705730,     0.000000],
     [     0.000000,     0.000000,    -2.705850],],
    ]

dp = DipoleDipole(
    struct=poscar,
    zeu=zeu,
    epsilon=epsilon,
    factor=factor,
)

dd = np.round(dp.c_dd("1/2 0 0"), 6)
print "real part"
print dd.real
print "imaginary part"
print dd.imag

real part
[[ 6.607385 -6.697344 -6.698868 -1.830306  3.397245  3.397396  1.830306 -3.397245 -3.397396]
 [-6.697344  6.608032  6.699215  3.397107 -1.83057  -3.397572 -3.397107  1.83057   3.397572]
 [-6.698868  6.699215  6.611206  3.39788  -3.398194 -1.831067 -3.39788   3.398194  1.831067]
 [-1.830306  3.397107  3.39788   2.038075 -1.676541 -1.676616  0.       -0.       -0.      ]
 [ 3.397245 -1.83057  -3.398194 -1.676541  2.03847   1.676771 -0.       -0.        0.      ]
 [ 3.397396 -3.397572 -1.831067 -1.676616  1.676771  2.038567 -0.        0.        0.      ]
 [ 1.830306 -3.397107 -3.39788   0.       -0.       -0.        2.038075 -1.676541 -1.676616]
 [-3.397245  1.83057   3.398194 -0.       -0.        0.       -1.676541  2.03847   1.676771]
 [-3.397396  3.397572  1.831067 -0.        0.        0.       -1.676616  1.676771  2.038567]]
imaginary part
[[ 0. -0. -0. -0. -0.  0. -0. -0.  0.]
 [-0. -0. -0. -0. -0. -0. -0.  0. -0.]
 [-0. -0. -0.  0. -0.  0.  0. -0. -0.]
 [ 0.  0. -0. -0. -0.  0. -0. -0. -0.]
 [ 0.  0.  0. -0. -0. -0. -0. -0. -0.]
 [-0.  0. -0.  0. -0.  0. -0. -0.  0.]
 [ 0.  0. -0.  0.  0.  0.  0.  0. -0.]
 [ 0. -0.  0.  0.  0.  0.  0.  0.  0.]
 [-0.  0.  0.  0.  0. -0. -0.  0.  0.]]

dd = np.round(dp.c_dd("0 0 0"), 6)
print "real part"
print dd.real
print "imaginary part"
print dd.imag

real part
[[-0.911526  0.        0.        0.455763  0.        0.        0.455763 -0.       -0.      ]
 [ 0.       -0.911658  0.        0.        0.455829  0.       -0.        0.455829 -0.      ]
 [ 0.        0.       -0.911905  0.        0.        0.455953 -0.       -0.        0.455953]
 [ 0.455763  0.        0.        0.155948 -0.       -0.       -0.611711 -0.        0.      ]
 [ 0.        0.455829  0.       -0.        0.155995 -0.       -0.       -0.611824 -0.      ]
 [ 0.        0.        0.455953 -0.       -0.        0.155925  0.       -0.       -0.611878]
 [ 0.455763 -0.       -0.       -0.611711 -0.        0.        0.155948  0.        0.      ]
 [-0.        0.455829 -0.       -0.       -0.611824 -0.        0.        0.155995  0.      ]
 [-0.       -0.        0.455953  0.       -0.       -0.611878  0.        0.        0.155925]]
imaginary part
[[ 0. -0. -0. -0. -0.  0.  0.  0. -0.]
 [-0. -0. -0. -0. -0. -0.  0.  0.  0.]
 [-0. -0. -0.  0. -0.  0. -0.  0.  0.]
 [ 0.  0. -0. -0. -0.  0. -0. -0.  0.]
 [ 0.  0.  0. -0. -0. -0. -0. -0. -0.]
 [-0.  0. -0.  0. -0.  0.  0. -0.  0.]
 [-0. -0.  0.  0.  0. -0.  0.  0. -0.]
 [-0. -0. -0.  0.  0.  0.  0.  0.  0.]
 [ 0. -0. -0. -0.  0. -0. -0.  0.  0.]]

# The limit of the dipole dipole effect with q->0 from (1, 0, 0) direction (primitive convention)
dd = np.round(dp.c_dd("0 0 0", q_direction="1 0 0"), 6)
print "real part"
print dd.real
print "imaginary part"
print dd.imag

real part
[[ 7.636021 -8.54799  -8.549934 -3.820737  4.276895  4.277085 -3.820737  4.276895  4.277085]
 [-8.54799   7.636775  8.550377  4.276722 -3.821288 -4.277307  4.276722 -3.821288 -4.277307]
 [-8.549934  8.550377  7.640416  4.277694 -4.27809  -3.822327  4.277694 -4.27809  -3.822327]
 [-3.820737  4.276722  4.277694  2.295562 -2.139812 -2.139907  1.527904 -2.139812 -2.139907]
 [ 4.276895 -3.821288 -4.27809  -2.139812  2.296005  2.140105 -2.139812  1.528186  2.140105]
 [ 4.277085 -4.277307 -3.822327 -2.139907  2.140105  2.296125 -2.139907  2.140105  1.528322]
 [-3.820737  4.276722  4.277694  1.527904 -2.139812 -2.139907  2.295562 -2.139812 -2.139907]
 [ 4.276895 -3.821288 -4.27809  -2.139812  1.528186  2.140105 -2.139812  2.296005  2.140105]
 [ 4.277085 -4.277307 -3.822327 -2.139907  2.140105  1.528322 -2.139907  2.140105  2.296125]]
imaginary part
[[ 0. -0. -0. -0. -0.  0.  0.  0. -0.]
 [-0. -0. -0. -0. -0. -0.  0.  0.  0.]
 [-0. -0. -0.  0. -0.  0. -0.  0.  0.]
 [ 0.  0. -0. -0. -0.  0. -0. -0.  0.]
 [ 0.  0.  0. -0. -0. -0. -0. -0. -0.]
 [-0.  0. -0.  0. -0.  0.  0. -0.  0.]
 [-0. -0.  0.  0.  0. -0.  0.  0. -0.]
 [-0. -0. -0.  0.  0.  0.  0.  0.  0.]
 [ 0. -0. -0. -0.  0. -0. -0.  0.  0.]]

Including LO-TO from command line

In order to execute LO-TO splitting on some calculation from the command line, do the following:

1. you can make the epsilon (dielectric constants + born effective charges) file from the OUTCAR using the following command:

~/git_repos/software/phonons/examples/fourier_interp_Dqn/parse_epsilon.py struct=poscar > epsilon.yml

2. Run the Fourier interpolation using the above generated epsilon file:

fourier_interp_Dqn.py epsilon=<path to epsilon file>

3. Regenerate the band structure along the desired path throught k-space using the new force constants and the epsilon file:

band_plot_phonons.py epsilon=<path to epsilon file> pos=poscar units=THz band=kinput phi=PhiN_WS.pkl