quantum-espresso/examples/example31

This example shows how to perform electronic structure calculations
using pw.x for a system undergoing the presence of a static homogeneous
finite electric field. The method is explained in:

P. Umari and A. Pasquarello, PRL 89,157602 (2002)
I. Souza, J.Iniguez, and D.Vanderbilt, PRL 89, 117602 (2002)


The concerned parameters are:

In namelist &CONTROL

 
lelfield       LOGICAL ( default = .FALSE. )
               If .TRUE. a homogeneous finite electric field
               described through the modern theory of the polarization
               is applied.
 
 
gdir           INTEGER
               For Berry phase calculation: direction of the k-point
               strings in reciprocal space. Allowed values: 1, 2, 3
               1=first, 2=second, 3=third reciprocal lattice vector
               For calculations with finite electric fields
               (lelfield==.true.), gdir is the direction of the field
 
nppstr         INTEGER
               For Berry phase calculation: number of k-points to be
               calculated along each symmetry-reduced string
               The same for calculation with finite electric fields
               (lelfield==.true.)
 
nberrycyc      INTEGER ( default = 1 )
               In the case of a finite electric field  (lelfield==.true.)
               it defines the number of iterations for converging the
               wavefunctions in the electric field Hamiltonian, for each
               external iteration on the charge density
 

In namelist &ELECTRONS

 
efield         REAL ( default = 0.D0 )
               For finite electric field calculations (lelfield == .true.),
               it defines the intensity of the field in a.u.


To perform a calculations with an electric field, an estimate of
the optimized wavefunctions is needed to build the electric field
operator (See: I. Souza, J.Iniguez and D. Vanderbilt, PRB 69, 085106,
2004). Therefore when lelfield ==.true. a copy of the wavefunctions
is read from disk (i.e. restart_mode must be 'restart').
The parameters GDIR defines the direction of the electric field.
The k_points must be given as a series of k-points-strings.
A k-points-string is a series of NPPSTR uniform spaced k-points
along the direction gdir. All the k-points in a string must have
the same weight. 
PAY ATTENTION: in pw.x the default units for k-points coordinates
is 2pi/alat and NOT  crystalline units.

Example of k-strings:

nppstr=4
gdir=1

0.0   KY  KZ  1.
0.25  KY  KZ  1.
0.50  KY  KZ  1.
0.75  KY  KZ  1.

nppstr=4
gdir=3

KX  KY  0.0  1.
KX  KY  0.25 1.
KX  KY  0.50 1.
KX  KY  0.75 1.
 

NOTE:It works fine in parallel with parallelization over the G points.
    However, for PARALLEL calculations gdir MUST BE 3.

For every usual iteration of pw.x when the Hartree and exchange-correlation 
potentials are kept fixed, when lelfield==.true. there are NBERRYCYC
iterations. During each of these iterations,  the electric field operator 
(which depends on the wave-functions) is kept fixed; then
the new electric field operator is built from the eigen-wavefunctions,
and a new iteration starts. This has been introduced because
the electric field Hamiltonian depends selfconsistently on the
wavefunctions.

NOTE:The option lelfield==.true. has been tested SO FAR only for
      orthorhombic unit cells

For every iteration on the charge (usual pw.x iterations), the code
reports the Electronic and Ionic Dipole in a.u. per unit cell and
the expectation value of the operator e^{+iGz}. The letter
is given for the corresponding supercell  containing N_kx*N_ky*N_kz
unit cells (N_kx,N_ky,N_kz are the number of k-points along x,y,z)


Example:

With this example, we show how to calculate the dielectric constant
of bulk silicon. The system is described by a 8-atom cubic unit cell.
We use a regular mesh of 3X3X7 k-points, where we have 7 k-points
along the directions of the electric field: gdir=3,nppstr=7
The first calculation just calculates the electronic structure
without electric field. The second calculation turns on the field
but with 0 a.u. intensity. The third calculation applies a field
of 0.001 a.u..

After the second calculation the electronic dipole D[0.a.u.] at 0 field reads:
Electronic Dipole per cell (a.u.) -2.502553344652915E-004
After the third calculation the electronic dipole D[0.001 a.u.] at 0.001 a.u. field reads:
Electronic Dipole per cell (a.u.)  0.934898995836340

The high-frequency dielectric constant eps_inf is then given by

eps_inf=4*pi*(D[0.001 a.u.]-D[0.0 a.u.])/(0.001 a.u. * Omega) + 1

where Omega is the volume of the unit cell (1054.9778 (a.u.)^3). 
We obtain:

eps_inf=12.14

(Compare: other DFT calculations, 12.7-13.1 , exp. 11.4 )
The result 12.14 is not fully converged on nppstr, for a discussion
on convergence see:
P.Umari and A. Pasquarello, PRB 68, 085114 (2003).