diff --git a/examples/EXX_example/README b/examples/EXX_example/README new file mode 100644 index 000000000..399765dc1 --- /dev/null +++ b/examples/EXX_example/README @@ -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] + + - 0.5* + 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.. + diff --git a/examples/EXX_example/run_example b/examples/EXX_example/run_example new file mode 100755 index 000000000..9798a16a5 --- /dev/null +++ b/examples/EXX_example/run_example @@ -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"