add EXX example

git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@2794 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
degironc 2006-02-09 13:44:51 +00:00
parent 56210fe18a
commit ed595a2422
2 changed files with 566 additions and 0 deletions

163
examples/EXX_example/README Normal file
View File

@ -0,0 +1,163 @@
Hybrid Hartree-Fock-DFT functionals are an experimental feature in PWscf.
Only few functionalities are implemented and it is possible (actually
desirable, see later) that the implementation evolves in the future
therefore its stability cannot be garanteed.
HOW TO COMPILE :
i) download the CVS version of QE and execute the ../configure step
ii) edit file make.sys and add to variable DFLAGS the -DEXX option
iii) issue "make pw" .
Pay attention that ./configure re-creates make.sys so if the QE distributions
is updated and ./configure run again also step ii) above has to be done again.
WHICH FUNCTIONALS ARE IMPLEMENTED ?
only a few hybrid functionals are implemented: Hartree-Fock, PBE0
and B3LYP . Their labels in QE are "HF", "PBE0" and "B3LP" .
This last one is not optimal but is needed to avoid conflict with LYP
correlation when parsing the dft string. It might change in the future.
Usually in PWscf the functional to be used is read from the
pseudopotential files but we do not have so far a pseudopotential
generator for hybrid functionals so one needs to use pseudopotentials
generated with some other functionals (eg. PBE or LDA) and force the
use of the hybrid functional in the PWscf input.
This is done using the variable input_dft in system namelist; for instace
input_dft="pbe0" will force the use of pbe0 unrespective of the functional
used in the pseudopotential generation.
HOW DOES THE SCF SEARCH PROCEED ?
use of hybrid-functionals is VERY expensive (see later), moreover in
more traditional DFT calculation the mixed quantity is the density while
for HF of hybrid-functionals the density-matrix should be mixed which is
not practical. The strategy used here is to consider (let's focus on HF
for simplicity) an auxiliary set of wavefunctions psi in addition to the
usual set phi and to minimize the auxiliary functional:
E[phi,psi] = T[phi] + E_ext[phi] + E_Hartree[phi] +
<phi|Vx[psi]|phi> - 0.5*<psi|Vx[psi]|psi>
where Vx[psi] is the fock operator defined with the auxiliary function psi.
Taking the functional derivatives w.r.t. phi it can be shown that the scf
condition for phi are the HF equation with fixed Fock operator, so Vx does
not enter in the scf procedure and one can mix density as usual. The
minimum condition w.r.t. psi is simply psi=phi so when both psi and phi
are minimized the standard HF energy is obtained.
Actually one can show that the functional E[phi,psi] above is
E[phi,psi] = E_HF[phi] + dexx[phi,psi]
where dexx is a positive definite addition to E_HF .
The scf procedure goes as follow.
1) a normal scf (with LDA or other scheme is performed) then
2) hybrid functional is switched on and psi = phi (the current best wfcs)
3) a new scf is performed w.r.t phi, keeping fixed Vx[psi]
4) dexx[phi,psi] is computed and if it exceeds the required tollerance
the proceedure is repeated from point 2)
HF may require several phi-scf cycles to reach full convergence. B3LYP
and PBE0, due to the smaller fraction of HF exchange included, require
usually a smaller number of phi-scf cycles
HOW MUCH EXPENSIVE IS THE CALCULATION ?
Very much.
Applying the Fock operator on a single vawefunction (phi_k,v) requires
the calculation of an integral over the whole BZ and all psi bands.
For each needed pair psi_k+q,v' and phi_k,v an auxiliary charge density
rho(-q+G) is built in real space and then FFT to reciprocal space where
the corresponding Poisson equation is solved. This auxiliary potential is
FFT back in real space where it is multiplied by psi_k+q,v' and added to
Vx[psi]phi...
The cost of the operation is therefore roughly NBND * NQS * ( 2 * FFT + ... )
where NQS is the number of q-points chosen to represent the BZ integration,
and depends in general on the localization of the Wannier functions of the
system.
For comparison non-local pseudopotentials in the KB formualtion (without
exploiting the locality of the KB projetors) cost NKB * (2 * NPW) where
NKB is tipically of the order of NBND but NPW cost at least an order
of magnitude less than an FFT.
Therefore even when one can take NQS=1 (for large non-metallic system
should be ok) hybrid-functionals will require at least an order of
magnitude more resources that a standard calculation.
HOW CAN I CHOSE NQS IN INPUT ?
In the system namelist there are three variables nqx1,nqx2,nqx3 that
define the regular q-grid in the BZ in a way similar to the automatic
k-points generation. Their vaule must be compatible with the k-points
used (that is k+q must be equivalent to some other k in the k-points list)
Their default value are nqx1=1,nqx2=1,nqx3=1 (BZ integration is
approximated by gamma point value only).
DIVERGENCE AT q->0
The BZ integral to be performed has a diverging kernel when (q+G)->0.
This is dealt with by adding and subtracting a term with the same
divergence that can be integrated analytically and performing
numerically the integration for the non divergent residue
[Gygi-Baldereschi, PRB 34, 4405 (1986)].
One problem is left: the now non divergent q=0 term is not easily determined
since it is a 0/0 (non analytic) limit. Several options have been considered:
1) just discard it ... this is not a good idea in general because it
induces an error proportional to 1/(NQS*Omega) in the total energy
where Omega is the Wigner-Seitze volume of the crystal. As one
wish to keep NQS as small as possible this may be large.
2) exploit the fact that the term has the above dependence and extract
it from a calculation with a given nqx1,nqx2,nqx3 and the one with
a grid twice as coarse in each direction. One does not really need to
perform two calculations but can do it internally (even when nqx? are
not even numbers...). This seems to work and it is set as the default.
In order to disable this feature [and get back to option 1)] set
x_gamma_extrapolation = .false.
3) perform calcualtions in q-grids that are shifted away from gamma so that
the 0/0 term is not needed. This create some extra complication in the
coding and cannot be used with Gamma-only k-point integration.
In some test does not seem to be superior to option 2) ... it has not
been fully implemented and now it has been removed.
4) use the value at small (q+G) to estimate the (q+G)->0 limit. This
again has been tried and found to offer, for low order numerical
differentiation, no better results that option 2). It is possible
than higher order formulas can be better but this has not been explored.
This option is currently not implemented but would be easy to resume it.
OTHER LIMITATIONS
So far only NORM-CONSERVING pseudopotentials are implemented.
there is no fundametal problem in defining HF for US pseudopotentials
but since some density-like object is required one would need to operate
on the dense charge-density FFT grid anyway with no computational gain.
Maybe this is not true and one can find ways to perform this integrals
more efficently. So far I did not think to much to this point.
PARALLEL IMPLEMENTATION ?
yes (and no).
At present only R-and-G parallelization has been implemented.
This is what is needed for large systems ...
For metals k-point parallelization is often useful but the need
in hybrid-functionals for the BZ integration prevent its simple
implementation, except when q-integration reduce to Gamma.
In this case it should work.
If the scratch area is common to all processors one could think about
implementing it in general (one would needs a VERYBIG direct-access
file that all processors can see) .. this feature is not implemented.
WHAT PROPERTIES CAN I COMPUTE ?
Energy and forces (thanks to Hellman-feynman theorem forces do not
require extra calculations). In principle also stresses but the
corresponding formulas have not been coded yet.
So structural optimization is OK if the cell shape is kept fixed.
Band structure ? yes and no. Obviously one compute wfc during the scf
cycle and their eigenvalues are printed in output.
This can be enough to draw a band structure or a DOS, however the
problem is when one wish non-scf calculations in k-points different
from those computed during the scf cycle. At present it is not possible
because it would require the knowledge of all bands at k+q that we do
not have. I do not know how to by-pass this problem.
ELECTRIC FIELD
I did not dig into this issue but Paolo Umari is using EXX with
electric field. For details it would be better to ask him directly.
AN EXAMPLE
run_example script in this directory performs two series of calculations:
1) total energy of Silicon using different values for nq,
2) calculation of binding energy of o2,co,n2 from calculations in a
12 au cubic box and gamma sampling.
Running it will generate directory "results" to be compared with directory
"reference"
Please report problems and suggestions to Stefano de Gironcoli
(degironc@sissa.it) and keep in mind that this feature is still
experimetal..

403
examples/EXX_example/run_example Executable file
View File

@ -0,0 +1,403 @@
#!/bin/sh
# run from directory where this script is
cd `echo $0 | sed 's/\(.*\)\/.*/\1/'` # extract pathname
EXAMPLE_DIR=`pwd`
# check whether echo has the -e option
if test "`echo -e`" = "-e" ; then ECHO=echo ; else ECHO="echo -e" ; fi
$ECHO
$ECHO "$EXAMPLE_DIR : starting"
$ECHO
$ECHO "This example shows how to use pw.x to calculate the total energy and"
$ECHO "of silicon and a few small molecules using hybrid functionals."
# set the needed environment variables
. ../environment_variables
# required executables and pseudopotentials
BIN_LIST="pw.x"
PSEUDO_LIST="Si.vbc.UPF"
$ECHO
$ECHO " executables directory: $BIN_DIR"
$ECHO " pseudo directory: $PSEUDO_DIR"
$ECHO " temporary directory: $TMP_DIR"
$ECHO " checking that needed directories and files exist...\c"
# check for directories
for DIR in "$BIN_DIR" "$PSEUDO_DIR" ; do
if test ! -d $DIR ; then
$ECHO
$ECHO "ERROR: $DIR not existent or not a directory"
$ECHO "Aborting"
exit 1
fi
done
for DIR in "$TMP_DIR" "$EXAMPLE_DIR/results" ; do
if test ! -d $DIR ; then
mkdir $DIR
fi
done
cd $EXAMPLE_DIR/results
# check for executables
for FILE in $BIN_LIST ; do
if test ! -x $BIN_DIR/$FILE ; then
$ECHO
$ECHO "ERROR: $BIN_DIR/$FILE not existent or not executable"
$ECHO "Aborting"
exit 1
fi
done
# check for pseudopotentials
for FILE in $PSEUDO_LIST ; do
if test ! -r $PSEUDO_DIR/$FILE ; then
$ECHO
$ECHO "ERROR: $PSEUDO_DIR/$FILE not existent or not readable"
$ECHO "Aborting"
exit 1
fi
done
$ECHO " done"
# how to run executables
PW_COMMAND="$PARA_PREFIX $BIN_DIR/pw.x $PARA_POSTFIX"
$ECHO
$ECHO " running pw.x as: $PW_COMMAND"
$ECHO
# clean TMP_DIR
$ECHO " cleaning $TMP_DIR...\c"
rm -rf $TMP_DIR/*
$ECHO " done"
$ECHO
$ECHO " running PBE0 calculation for Si with nq=1,2,4 \c"
$ECHO
for nq in 1 2 4 ; do
# self-consistent calculation
cat > si.in << EOF
&control
calculation = 'scf'
restart_mode='from_scratch',
prefix='silicon',
pseudo_dir = '$PSEUDO_DIR/',
outdir='$TMP_DIR/'
/
&system
ibrav= 2, celldm(1) =10.20, nat= 2, ntyp= 1,
ecutwfc =12.0, nbnd = 8,
input_dft='pbe0', nqx1 = $nq, nqx2 = $nq, nqx3 = $nq,
/
&electrons
mixing_beta = 0.7
/
ATOMIC_SPECIES
Si 28.086 Si.vbc.UPF
ATOMIC_POSITIONS
Si 0.00 0.00 0.00
Si 0.25 0.25 0.25
K_POINTS
10
0.1250000 0.1250000 0.1250000 1.00
0.1250000 0.1250000 0.3750000 3.00
0.1250000 0.1250000 0.6250000 3.00
0.1250000 0.1250000 0.8750000 3.00
0.1250000 0.3750000 0.3750000 3.00
0.1250000 0.3750000 0.6250000 6.00
0.1250000 0.3750000 0.8750000 6.00
0.1250000 0.6250000 0.6250000 3.00
0.3750000 0.3750000 0.3750000 1.00
0.3750000 0.3750000 0.6250000 3.00
EOF
$ECHO " running the scf calculation for Si with nq = $nq ...\c"
$PW_COMMAND < si.in > si.PBE0_nq=${nq}.out
$ECHO " done"
grep -e ! si.PBE0_nq=${nq}.out | tail -1
rm -f si.in
done
$ECHO
$ECHO " running now a few molecules with Gamma sampling ...\c"
$ECHO
PSEUDO_DIR=../Pseudo
$ECHO " pseudo directory cahnged to: $PSEUDO_DIR"
xc=pbe0
ps=1nlcc
ecut=80
rm -fr $TMP_DIR/*
cat > o.inp << EOF
&CONTROL
calculation = 'scf' ,
restart_mode = 'from_scratch' ,
outdir = '$TMP_DIR/' ,
pseudo_dir = '$PSEUDO_DIR/',
prefix = 'o',
disk_io = 'minimal' ,
iprint = 1
tprnfor = .true.
/
&SYSTEM
ibrav = 1,
celldm(1) = 12.0,
nat = 1,
ntyp = 1,
ecutwfc = $ecut ,
input_dft='$xc'
nspin = 2
starting_magnetization(1) = 0.2,
nbnd = 4
nelec = 6.0
nelup = 4.0
neldw = 2.0
/
&ELECTRONS
conv_thr = 0.5d-3
/
ATOMIC_SPECIES
O 16.0 OPBE$ps.RRKJ3
ATOMIC_POSITIONS angstrom
O 0.1 0.2 0.3
K_POINTS gamma
#automatic
#1 1 1 0 0 0
EOF
$ECHO " running oxygen atom..\c"
$PW_COMMAND < o.inp > o.$xc.$ps.out-$ecut
$ECHO " done"
rm -fr $TMP_DIR/*
cat > c.inp << EOF
&CONTROL
calculation = 'scf' ,
restart_mode = 'from_scratch' ,
outdir = '$TMP_DIR/' ,
pseudo_dir = '$PSEUDO_DIR/',
prefix = 'o2',
disk_io = 'minimal' ,
iprint = 1
tprnfor = .true.
/
&SYSTEM
ibrav = 1,
celldm(1) = 12.0,
nat = 1,
ntyp = 1,
ecutwfc = $ecut ,
input_dft='$xc'
nspin = 2
starting_magnetization(1) = 0.2,
nbnd = 4
nelec = 4.0
nelup = 3.0
neldw = 1.0
/
&ELECTRONS
conv_thr = 0.5d-3
/
ATOMIC_SPECIES
C 16.0 CPBE$ps.RRKJ3
ATOMIC_POSITIONS angstrom
C 0.1 0.2 0.3
K_POINTS gamma
#automatic
#1 1 1 0 0 0
EOF
$ECHO " running carbon atom..\c"
$PW_COMMAND < c.inp > c.$xc.$ps.out-$ecut
$ECHO " done"
rm -fr $TMP_DIR/*
cat > n.inp << EOF
&CONTROL
calculation = 'scf' ,
restart_mode = 'from_scratch' ,
outdir = '$TMP_DIR/' ,
pseudo_dir = '$PSEUDO_DIR/',
prefix = 'o2',
disk_io = 'minimal' ,
iprint = 1
tprnfor = .true.
/
&SYSTEM
ibrav = 1,
celldm(1) = 12.0,
nat = 1,
ntyp = 1,
ecutwfc = $ecut ,
input_dft='$xc'
nspin = 2
starting_magnetization(1) = 0.2,
nbnd = 4
nelec = 5.0
nelup = 4.0
neldw = 1.0
/
&ELECTRONS
conv_thr = 0.5d-4
/
ATOMIC_SPECIES
N 16.0 NPBE$ps.RRKJ3
ATOMIC_POSITIONS angstrom
N 0.1 0.2 0.3
K_POINTS gamma
#automatic
#1 1 1 0 0 0
EOF
$ECHO " running nitrogen atom..\c"
$PW_COMMAND < n.inp > n.$xc.$ps.out-$ecut
$ECHO " done"
rm -fr $TMP_DIR/*
b=0.3169
cat > n2.inp << EOF
&CONTROL
calculation = 'scf' ,
restart_mode = 'from_scratch' ,
outdir = '$TMP_DIR/' ,
pseudo_dir = '$PSEUDO_DIR/',
prefix = 'o2',
disk_io = 'minimal' ,
iprint = 1
tprnfor = .true.
/
&SYSTEM
ibrav = 1,
celldm(1) = 12.0,
nat = 2,
ntyp = 1,
ecutwfc = $ecut ,
input_dft='$xc'
nbnd = 8
/
&ELECTRONS
conv_thr = 1.d-4
/
&IONS
/
ATOMIC_SPECIES
N 16.0 NPBE$ps.RRKJ3
ATOMIC_POSITIONS angstrom
N $b $b $b
N -$b -$b -$b
K_POINTS gamma
#automatic
#1 1 1 0 0 0
EOF
$ECHO " running n2 molecule..\c"
$PW_COMMAND < n2.inp > n2.$xc.$ps.out-$ecut
$ECHO " done"
rm -fr $TMP_DIR/*
b=0.3256
cat > co.inp << EOF
&CONTROL
calculation = 'scf' ,
restart_mode = 'from_scratch' ,
outdir = '$TMP_DIR/' ,
pseudo_dir = '$PSEUDO_DIR/',
prefix = 'o2',
disk_io = 'minimal' ,
iprint = 1
tprnfor = .true.
/
&SYSTEM
ibrav = 1,
celldm(1) = 12.0,
nat = 2,
ntyp = 2,
ecutwfc = $ecut ,
input_dft='$xc'
nbnd = 8
/
&ELECTRONS
conv_thr = 0.5d-3
/
&IONS
/
ATOMIC_SPECIES
C 16.0 CPBE$ps.RRKJ3
O 16.0 OPBE$ps.RRKJ3
ATOMIC_POSITIONS angstrom
C $b $b $b
O -$b -$b -$b
K_POINTS gamma
#automatic
#1 1 1 0 0 0
EOF
$ECHO " running co molecule..\c"
$PW_COMMAND < co.inp > co.$xc.$ps.out-$ecut
$ECHO " done"
rm -fr $TMP_DIR/*
b=0.3478
cat > o2.inp << EOF
&CONTROL
calculation = 'scf' ,
restart_mode = 'from_scratch' ,
outdir = '$TMP_DIR/' ,
pseudo_dir = '$PSEUDO_DIR/',
prefix = 'o2',
disk_io = 'minimal' ,
iprint = 1
tprnfor = .true.
/
&SYSTEM
ibrav = 1,
celldm(1) = 12.0,
nat = 2,
ntyp = 1,
ecutwfc = $ecut ,
input_dft='$xc'
nspin = 2
starting_magnetization(1) = 0.2,
nbnd = 8
nelec = 12.0
nelup = 7.0
neldw = 5.0
/
&ELECTRONS
conv_thr = 0.5d-3
/
&IONS
/
ATOMIC_SPECIES
O 16.0 OPBE$ps.RRKJ3
ATOMIC_POSITIONS angstrom
O $b $b $b
O -$b -$b -$b
K_POINTS gamma
#automatic
#1 1 1 0 0 0
EOF
$ECHO " running o2 molecule..\c"
$PW_COMMAND < o2.inp > o2.$xc.$ps.out-$ecut
$ECHO " done"
$ECHO " summarize"
cat > summarize << EOF
grep -e ! n.$xc.$ps.out-$ecut | tail -1 | awk '{print \$5}' > N
grep -e ! n2.$xc.$ps.out-$ecut | tail -1 | awk '{print \$5}' > N2
paste N2 N | awk '{be= (\$1-\$2*2.0) * 13.6058 * 23.06; print "N2 : ",be}'
grep -e ! o.$xc.$ps.out-$ecut | tail -1 | awk '{print \$5}' > O
grep -e ! o2.$xc.$ps.out-$ecut | tail -1 | awk '{print \$5}' > O2
paste O2 O | awk '{be= (\$1-\$2*2.0) * 13.6058 * 23.06 ; print "O2 : ",be}'
grep -e ! c.$xc.$ps.out-$ecut | tail -1 | awk '{print \$5}' > C
grep -e ! co.$xc.$ps.out-$ecut | tail -1 | awk '{print \$5}' > CO
paste CO O C | awk '{be= (\$1-\$2-\$3) * 13.6058 * 23.06; print "CO : ",be}'
rm C N O CO O2 N2
EOF
sh summarize
rm -f *.inp
rm -fr $TMP_DIR/*
$ECHO "$EXAMPLE_DIR : done"