I'm trying to extract the first order Hamiltonian from a response function calculation in which the perturbation is a displacement of an atom. The code produces a POT file for each perturbation specified in the input file (reduced by symmetry if possible). Using cut3d one can read the files easily enough, and, for qpt 0 0 0, the output looks convincing and makes a nice, smooth function that meets itself at the cell boundary. I am guessing that this quantity is H^(1) defined by Eq. (17) of PRB 55 10337 (Gonze, First-principles responses of solids to atomic displacements and homogeneous electric fields: Implementation of a conjugate-gradient algorithm). However at non-trivial wave vectors, the contents of the corresponding POT files (while still having the correct length) contain data that appear to be two smooth functions interleaved (each pair of successive numbers in the file has one point belonging to each of the two smooth functions). Also, the large values for these non-zero wave-vector perturbations do not seem to be clustered close to the atom corresponding to the perturbation as the gamma point ones are.
So I surmise that I (or rather cut3d) am not interpreting the POT file correctly. Can somebody point me some documentation or comments in the code from which the structure of the POT file can be deduced in this case (RF calculation at non-zero wave vector)? I would be very grateful for any hint.
I have attached a plot of the same perturbation at vanishing wave vector (ds21.ps) and at q=(0.5,0.5,0.5) (ds31.ps).
Here is the input file used to generate the data for the above plots:
Code: Select all
ndtset 3
getwfk 1 # Use GS wave functions from dataset1
nqpt 1 # read a q for the perturbation
kptopt 3 # Need full k-point set for finite-Q response
rfphon 1 # Do phonon response
rfatpol 1 2 # Treat displacements of all atoms
rfdir 1 1 1 # Do all directions (symmetry will be used)
tolvrs 1.0d-8
prtpot 1
#first, calculate ground state
getwfk1 0 # Cancel default
kptopt1 1 # Automatic generation of k points, taking
# into account the symmetry
nqpt1 0 # Cancel default
tolvrs1 1.0d-18 # SCF stopping criterion (modify default)
rfphon1 0 # Cancel default
# Definition of the unit cell: fcc
acell 3*8.64
rprim 1.0 0.0 0.0 # primitive vectors (to be scaled by acell)
0.0 1.0 0.0
0.0 0.0 1.0
# Definition of the atom types
ntypat 2
znucl 55 53
# Definition of the atoms
natom 2 # There are two atoms
typat 1 2
xred # Reduced coordinate of atoms
0.0 0.0 0.0
0.5 0.5 0.5
# Definition of the k-point grid
prtvol 1
ngkpt 2 2 2
nshiftk 1
shiftk 0.0 0.0 0.0
# Definition of the planewave basis set
ecut 20.0 # Maximal kinetic energy cut-off, in Hartree
# Definition of the SCF procedure
diemac 5.65 # Although this is not mandatory, it is worth to
# precondition the SCF cycle. The model dielectric
# function used as the standard preconditioner
# is described in the "dielng" input variable section.
#list of q's
qpt2 0.0 0.0 0.0
qpt3 0.5 0.5 0.5
Here are the options I used to run cut3d:
Code: Select all
cut3d <<EOF > cut3d.out
phonono_DS2_POT1
1
8
pot.txt
0
EOF
The attached plots were made by using the last (36th in this case) plane of constant z in the abinit grid (this plane contains the perturbed atom). These results are from abinit 6.8.2/gfortran. If you have read this far thank you very much for your attention.
Happy computing,
Micah Prange