[QE-users] strange and possibly wrong values of forces in assume_isolated='2D' and/or 'esm' calculations
Giuseppe Mattioli
giuseppe.mattioli at ism.cnr.it
Thu Feb 1 15:42:57 CET 2024
Dear all
I'm studying the absorption of molecules on metal surfaces. I've often
used without problems assume_isolated='esm' and esm_bc='bc1' to
decouple the interaction along z between periodically repeated images.
However, I plan to use tot_charge .neq. 0 this time, which seems to be
not compatible with 'esm', thus I also tried with assume_isolated='2D'.
I'm using version 7.2 together with well-tested gbrv or even older
ultrasoft pseudopotentials in a setup used dozens of times for this
kind of study. I've been able to reproduce part of the problem in the
case of a single, flat benzene molecule.
This is my input file. It is a hexagonal cell which I need for
Au(111), but there is no difference if I use a cubic one.
&control
calculation = 'relax'
restart_mode='from_scratch',
prefix='ben-test1',
pseudo_dir = '/home/giuseppe/PP_UPF//',
outdir='/mnt/giuseppe-jobid_27303/',
wf_collect=.true.,
max_seconds=7000,
/
&system
ibrav= 4, celldm(1)=38.0226, celldm(3)=1.40000,
nat=12, ntyp=2,
ecutwfc=40.0,
ecutrho=320.0,
occupations='smearing', smearing='gaussian', degauss=0.01,
nspin=1,
input_dft='vdw-df-c09'
assume_isolated = '2D',
/
&electrons
diagonalization='david',
mixing_mode='plain',
mixing_beta=0.1,
conv_thr=1.0d-6,
electron_maxstep=100
/
&ions
ion_dynamics='bfgs'
/
ATOMIC_SPECIES
C 12.011 c_pbe_v1.2.uspp.F.UPF
H 1.008 h_pbe_v1.4.uspp.F.UPF
ATOMIC_POSITIONS {angstrom}
C 6.03754185 4.09980430 0.00 1 1 0
C 6.61183619 5.36437480 0.00 1 1 0
C 4.65518189 3.96484180 0.00 1 1 0
C 3.84702389 5.09448512 0.00 1 1 0
C 4.42130535 6.35904350 0.00 1 1 0
C 5.80366482 6.49400044 0.00 1 1 0
H 7.68628666 5.46931959 0.00 1 1 0
H 6.25002662 7.47692529 0.00 1 1 0
H 6.66555853 3.22167059 0.00 1 1 0
H 4.20877138 2.98189546 0.00 1 1 0
H 2.77257143 4.98953422 0.00 1 1 0
H 3.79328781 7.23720529 0.00 1 1 0
K_POINTS {gamma}
After the end of scf, pw.x calculate strange forces
atom 1 type 1 force = 1.97313636 -2.75871856 0.00000000
atom 2 type 1 force = 3.37564882 0.32985501 0.00000000
atom 3 type 1 force = -1.40228454 -3.08800275 0.00000000
atom 4 type 1 force = -3.37564436 -0.32996097 0.00000000
atom 5 type 1 force = -1.97318664 2.75878311 0.00000000
atom 6 type 1 force = 1.40228523 3.08791716 0.00000000
atom 7 type 2 force = 1.35992166 0.13283850 0.00000000
atom 8 type 2 force = 0.56499090 1.24423036 0.00000000
atom 9 type 2 force = 0.79488126 -1.11132340 0.00000000
atom 10 type 2 force = -0.56496888 -1.24411482 0.00000000
atom 11 type 2 force = -1.35992315 -0.13282232 0.00000000
atom 12 type 2 force = -0.79485668 1.11131869 0.00000000
which are so large in the x,y components because ionic and non local
contributions are not almost compensated.
The ionic contribution to forces
atom 1 type 1 force = 4.08723253 -5.71437365 -0.00000000
atom 2 type 1 force = 6.99250543 0.68334134 -0.00000000
atom 3 type 1 force = -2.90464852 -6.39651738 -0.00000000
atom 4 type 1 force = -6.99247323 -0.68344291 -0.00000000
atom 5 type 1 force = -4.08729816 5.71451843 -0.00000000
atom 6 type 1 force = 2.90458905 6.39635110 -0.00000000
atom 7 type 2 force = 3.36898875 0.32905678 -0.00000000
atom 8 type 2 force = 1.39967299 3.08228687 -0.00000000
atom 9 type 2 force = 1.96916090 -2.75319510 -0.00000000
atom 10 type 2 force = -1.39966185 -3.08209645 -0.00000000
atom 11 type 2 force = -3.36898333 -0.32906411 -0.00000000
atom 12 type 2 force = -1.96908455 2.75313508 -0.00000000
The local contribution to forces
atom 1 type 1 force = -2.10775560 2.94695248 0.00000000
atom 2 type 1 force = -3.60605113 -0.35240198 0.00000001
atom 3 type 1 force = 1.49790156 3.29863450 0.00000001
atom 4 type 1 force = 3.60602686 0.35255311 -0.00000002
atom 5 type 1 force = 2.10781420 -2.94702214 0.00000002
atom 6 type 1 force = -1.49789744 -3.29849178 -0.00000002
atom 7 type 2 force = -1.81828506 -0.17758491 -0.00000000
atom 8 type 2 force = -0.75543049 -1.66354139 0.00000000
atom 9 type 2 force = -1.06277449 1.48594231 0.00000000
atom 10 type 2 force = 0.75542343 1.66344016 -0.00000000
atom 11 type 2 force = 1.81828464 0.17760125 0.00000000
atom 12 type 2 force = 1.06275927 -1.48594251 0.00000000
In this conditions, bfgs goes bananas and after 50 scf steps I have
this unpleasant situation.
Total force = 8.956657 Total SCF correction = 0.001270
Energy error = 2.8E-03 Ry
Gradient error = 3.4E+00 Ry/Bohr
number of scf cycles = 50
number of bfgs steps = 1
energy old = -256.3549703404 Ry
energy new = -256.3521766047 Ry
CASE: energy _new > energy _old
new trust radius = 0.0000268942 bohr
trust_radius < trust_radius_min
resetting bfgs history
The maximum number of steps has been reached.
End of BFGS Geometry Optimization
If I only drop assume_isolated = '2D' from the input above, everything
goes as expected
Forces acting on atoms (cartesian axes, Ry/au):
atom 1 type 1 force = 0.01500130 -0.02089768 0.00000000
atom 2 type 1 force = 0.02582976 0.00235533 0.00000000
atom 3 type 1 force = -0.01069085 -0.02372350 0.00000000
atom 4 type 1 force = -0.02581131 -0.00235619 0.00000000
atom 5 type 1 force = -0.01500318 0.02083674 0.00000000
atom 6 type 1 force = 0.01071567 0.02379058 0.00000000
atom 7 type 2 force = -0.00097453 -0.00008554 0.00000000
atom 8 type 2 force = -0.00040726 -0.00089695 0.00000000
atom 9 type 2 force = -0.00047082 0.00066260 0.00000000
atom 10 type 2 force = 0.00040426 0.00088135 0.00000000
atom 11 type 2 force = 0.00094869 0.00008311 0.00000000
atom 12 type 2 force = 0.00045827 -0.00064987 0.00000000
...
The ionic contribution to forces
atom 1 type 1 force = 4.16313923 -5.81989640 0.00000000
atom 2 type 1 force = 7.12215767 0.69533929 0.00000000
atom 3 type 1 force = -2.95838629 -6.51519380 0.00000000
atom 4 type 1 force = -7.12195513 -0.69531315 0.00000000
atom 5 type 1 force = -4.16328737 5.81985352 0.00000000
atom 6 type 1 force = 2.95844582 6.51515895 0.00000000
atom 7 type 2 force = 3.29424118 0.32170505 0.00000000
atom 8 type 2 force = 1.36853432 3.01372692 0.00000000
atom 9 type 2 force = 1.92589235 -2.69253439 0.00000000
atom 10 type 2 force = -1.36858820 -3.01373442 0.00000000
atom 11 type 2 force = -3.29428594 -0.32169572 0.00000000
atom 12 type 2 force = -1.92590764 2.69258415 0.00000000
The local contribution to forces
atom 1 type 1 force = -4.15700443 5.81134910 -0.00000041
atom 2 type 1 force = -7.11159870 -0.69439936 -0.00000386
atom 3 type 1 force = 2.95401211 6.50546039 0.00000456
atom 4 type 1 force = 7.11140469 0.69436550 0.00000203
atom 5 type 1 force = 4.15710464 -5.81134937 -0.00000062
atom 6 type 1 force = -2.95410437 -6.50540079 0.00000027
atom 7 type 2 force = -3.11354868 -0.30405716 -0.00000088
atom 8 type 2 force = -1.29348330 -2.84844469 -0.00000016
atom 9 type 2 force = -1.82014254 2.54468681 -0.00000009
atom 10 type 2 force = 1.29353174 2.84844194 0.00000054
atom 11 type 2 force = 3.11358189 0.30404786 0.00000091
atom 12 type 2 force = 1.82014521 -2.54472427 0.00000036
and BFGS converges after only 5 steps
Total force = 0.000993 Total SCF correction = 0.000134
SCF correction compared to forces is large: reduce conv_thr to
get better values
Energy error = 7.3E-06 Ry
Gradient error = 4.5E-04 Ry/Bohr
bfgs converged in 5 scf cycles and 4 bfgs steps
(criteria: energy < 1.0E-04 Ry, force < 1.0E-03 Ry/Bohr)
End of BFGS Geometry Optimization
Final energy = -75.8258627406 Ry
File /mnt/giuseppe-jobid_27305/ben-test3.bfgs deleted, as requested
BFGS converges without problems *in this case* also if I use
assume_isolated = 'esm' instead of '2D', possibly because even if I
ask for a vacuum-slab-vacuum geometry, there is no dipole in this cell
and forces along z are killed by symmetry (am I right?). HOWEVER
(sorry, I prefer to put together all the problems in a single report
:-D ), when I have a surface+molecule system, which is not anymore
symmetric along z like benzene (I can share all the details, if
needed; it is a phthalocyanine on Au(111)), both '2D' and 'esm' seem
to fail in the calculation of forces, I report here only a few atoms
to show where the calculations fail
This is assume_isolated='2D'
Forces acting on atoms (cartesian axes, Ry/au):
atom 1 type 1 force = -2.46581549 0.05555305 -44.68953990
atom 2 type 1 force = 1.33203881 -0.12969950 -39.29480550
atom 3 type 1 force = -0.24014829 0.13983948 -38.22826314
atom 4 type 1 force = -1.62848060 -1.75347368 -39.49598715
...
The ionic contribution to forces
atom 1 type 1 force = -4.92903372 0.11027640 -89.37571261
atom 2 type 1 force = 2.66753532 -0.25975231 -78.57707321
atom 3 type 1 force = -0.48180291 0.28021951 -76.44140193
atom 4 type 1 force = -3.25281219 -3.50837127 -78.97964210
...
The local contribution to forces
atom 1 type 1 force = 2.46715770 -0.05676727 44.60116355
atom 2 type 1 force = -1.33007965 0.12925384 39.20659054
atom 3 type 1 force = 0.23794546 -0.13850138 38.14136450
atom 4 type 1 force = 1.63063742 1.75189936 39.40956985
and this is assume_isolated='esm' with esm_bc='bc1', where at least
forces are compensated in x,y, but wild in z.
Forces acting on atoms (cartesian axes, Ry/au):
atom 1 type 1 force = 0.00087886 -0.00095173 29.72776373
atom 2 type 1 force = 0.00134883 0.00103179 29.72611292
atom 3 type 1 force = -0.00054959 0.00029250 29.75075633
atom 4 type 1 force = 0.00146407 0.00086933 29.74575523
...
The ionic contribution to forces
atom 1 type 1 force = -4.86679572 0.10850005 -89.30125557
atom 2 type 1 force = 2.64875995 -0.25726046 -78.60005420
atom 3 type 1 force = -0.47928395 0.27875618 -76.47410067
atom 4 type 1 force = -3.23123811 -3.48390938 -78.99960387
...
The local contribution to forces
atom 1 type 1 force = 4.87145206 -0.11217577 118.94587481
atom 2 type 1 force = -2.64125008 0.25728307 108.25252541
atom 3 type 1 force = 0.47492849 -0.27616873 106.15561996
atom 4 type 1 force = 3.24009805 3.48198049 108.67291065
For the sake of completeness, I also attach a few details of the input
file, where you can see that all atoms are placed quite close to z=0,
as requested by assume_isolated='esm'
&control
calculation = 'relax'
restart_mode='from_scratch',
prefix='au+fepc-test1',
pseudo_dir = '/home/giuseppe/PP_UPF//',
outdir='/mnt/giuseppe-jobid_27219/',
wf_collect=.true.,
max_seconds=160000,
/
&system
ibrav= 4, celldm(1)=38.0226, celldm(3)=1.40000,
nat=155, ntyp=5,
ecutwfc=40.0,
ecutrho=320.0,
occupations='smearing', smearing='gaussian', degauss=0.01,
nspin=2,
tot_magnetization=2.0,
input_dft='vdw-df-c09'
assume_isolated = '2D',
/
&electrons
diagonalization='david',
mixing_mode='plain',
mixing_beta=0.1,
conv_thr=1.0d-6,
electron_maxstep=30,
scf_must_converge=.false.,
/
&ions
ion_dynamics='bfgs'
/
ATOMIC_SPECIES
Au 196.967 Au-pbe-n-van.UPF
Fe 55.845 fe_pbe_v1.5.uspp.F.UPF
N 14.007 N.pbe-van_bm.UPF
C 12.011 C_pbe.van.UPF
H 1.008 H_pbe.van.UPF
ATOMIC_POSITIONS {angstrom}
Au 17.24592053 1.65974626 -3.36715305 0 0 0
Au 8.62333050 1.65969366 -3.36703580 0 0 0
Au 2.87432523 11.61662117 -3.36637204 0 0 0
Au 4.31134743 14.10584071 -3.36658733 0 0 0
...
Fe 0.00000000 1.65973416 1.50000000 0 0 1
...
H -7.55645579 0.91853478 1.50000000
H -6.91357881 -1.47304547 1.50000000
H -4.49348975 -2.15506125 1.50000000
H -5.80254887 2.72071813 1.50000000
K_POINTS {gamma}
HUBBARD {atomic}
U Fe-3d 4.0
As you have been so patient that you arrived all the way here, I hope
that you can also provide suggestions to solve my problem. I'm
available to provide further details, if needed.
Thank you in advance
Giuseppe
GIUSEPPE MATTIOLI
CNR - ISTITUTO DI STRUTTURA DELLA MATERIA
Via Salaria Km 29,300 - C.P. 10
I-00015 - Monterotondo Scalo (RM)
Mob (*preferred*) +39 373 7305625
Tel + 39 06 90672342 - Fax +39 06 90672316
E-mail: <giuseppe.mattioli at ism.cnr.it>
More information about the users
mailing list