[QE-users] Constrained DFT (Hubbard alpha approach)

[Christoph Wolf]

Dear all,

Some time ago I read about the implementation of a constrained DFT approach
here: http://theossrv1.epfl.ch/Main/OxidationStates

I was naturally curious on how to simply test this and decided to use H-H
dissociation as outlined in the text. However I am quite skeptical if my
approach works and I was hoping someone would give it a look.

I repeat briefly the "algorithm" given in above link:

To the total energy calculated by QE add a penalty term -alpha*N_target
then the resulting "CDFT" energy term reads

E=E_KS+alpha[Tr[n]-N_target]
with "Tr[n]" the output of the lda_plus_u implementation, i.e. the
occupation of the target manifold.

find the "correct" alpha by maximizing partial E/partial alpha -> alpha
[Tr[n]-N_target]=0 (for alpha is not zero)

I attempted to do this as follows:

for alpha in `seq 1E-40 1.0 30` ;do
...
  lda_plus_u=.true.
  Hubbard_U(1)=1D-40
  Hubbard_alpha(1)=-$alpha
  Hubbard_U(2)=1D-40
  Hubbard_alpha(2)=$alpha
....
  H1   1.0  H.pbe-kjpaw_psl.1.0.0.UPF
  H2   1.0  H.pbe-kjpaw_psl.1.0.0.UPF
[...]

and I then collect for atom 1 and 2 Tr[n].

The main problem I encounter is "double counting" of electrons, i.e. the
total number of electrons is not 2 but over 3 and the "most integer"
occupation I can get for H1 is ~1.96 (ok, close enough to 2) but only for
fairly large distance H-H. I attached a plot and the run-script, if anyone
could give it a glance and tell me if this is totally wrong or a little bit
right I would be very grateful!

Thanks in advance for your help!

Chris
[image: ECDFT.png]





-- 
Postdoctoral Researcher
Center for Quantum Nanoscience, Institute for Basic Science
Ewha Womans University, Seoul, South Korea
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.quantum-espresso.org/pipermail/users/attachments/20181009/e886e430/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: ECDFT.png
Type: image/png
Size: 75450 bytes
Desc: not available
URL: <http://lists.quantum-espresso.org/pipermail/users/attachments/20181009/e886e430/attachment.png>
-------------- next part --------------
printf "#%-1s %8s %8s %8s %9s %9s %8s\n" "a" "occ1" "occ2" "n_tot" "e_TOT" "E_CDFT" "dist (A)">occ2.dat
printf "#%-1s %8s %8s %8s %9s %9s %8s\n" "1" "2" "3" "4" "5" "6" "7">>occ2.dat.dat


n_target=1.961980

for alpha in `seq 1E-40 1.0 30` ;do

cat >relax.in <<EOF
&control
 calculation='relax',
 prefix='FeO',
 pseudo_dir= './',
 outdir = './',
/
&SYSTEM
  ibrav = 1
  A =    20
  nat = 2
  ntyp = 2
  ecutwfc = 50,
  ecutrho =250
  occupations='smearing'
  degauss = 0.005
  lda_plus_u=.true.
  Hubbard_U(1)=1D-40
  Hubbard_alpha(1)=-$alpha
  Hubbard_U(2)=1D-40
  Hubbard_alpha(2)=$alpha
/
&electrons
 conv_thr=1.d-9,
 startingpot='file'
/
&IONS
/
ATOMIC_SPECIES
  H1   1.0  H.pbe-kjpaw_psl.1.0.0.UPF
  H2   1.0  H.pbe-kjpaw_psl.1.0.0.UPF
ATOMIC_POSITIONS {crystal}
H1       0.500000000   0.500000000   0.500000000    0   0   0
H2       0.500000000   0.500000000   0.537541992    0   0   1
K_POINTS gamma
EOF


mpirun -np 4 pw.x <relax.in >relax_$alpha.out

occ1=$(grep "atom    1   Tr" relax_$alpha.out |tail -1|awk '{print $5}')
occ2=$(grep "atom    2   Tr" relax_$alpha.out |tail -1|awk '{print $5}')
etot=$(grep "!" relax_$alpha.out |tail -1|awk '{print $5}')
pos=$(grep "H2" relax_$alpha.out |tail -1|awk '{print $4}')

awk -v alpha="$alpha" -v occ1="$occ1" -v occ2="$occ2" -v etot="$etot"  -v n_target="$n_target" -v pos="$pos" 'BEGIN{printf "%02i %1.6f %1.6f %1.6f %1.6f %1.6f %1.6f\n", alpha, occ1, occ2, occ1+occ2, etot, alpha*(occ1-n_target), (pos-0.5)*20 }' >>occ2.dat


done

cat >plot <<EOF
set terminal pngcairo enhanced size 1024,800
set output "ECDFT.png"


set multiplot layout 1,2 title "Hubbard-{/Symbol a} approach to H-H dissociation \n{/*0.5 number labels is the total electron count}"
set xrange [0:12]

set ylabel "E_{CDFT}"
set y2label "H-H length (A)"
set xlabel "{/Symbol a}"

set ytics nomirror
set y2tics

plot "occ2.dat" u 1:6 axes x1y1 w lp pt 7 lw 2 lt 1  title "E_{CDFT}","occ2.dat" u 1:6:4 axes x1y1 w labels  offset 1 rotate by 45 notitle, "occ2.dat" u 1:7 axes x1y2 w lp pt 8 lt 2 title "H-H", 0.74 axes x1y2 w l lt 2 dt 2 notitle

set ylabel "E_{CDFT}"
set xlabel "H-H length (A)"
unset y2label
set ytics mirror

set xrange [0.74:2.2]
plot "occ2.dat" u 7:6 axes x1y1 w lp pt 7 lw 2 lt 1  notitle

EOF

gnuplot plot