Dear all,
I encounter a problem which I think has already been mentionned in
previous messages: when simulating surfaces (i.e. slabs + vacuum), the
SCF loop goes completely crazy after some iterations. In my case, this
is always the same scenario:
- for small slabs (~ 4-5 layers) it works perfectly
- when the thickness increases, the SCF loop goes crazy, like that:
siesta: iscf Eharris(eV) E_KS(eV) FreeEng(eV) dDmax
Ef_up Ef_dn(eV)
siesta: 1 -33906.0903 -33846.2100 -33846.2100 0.9317
-4.0527 -4.0527
timer: Routine,Calls,Time,% = IterSCF 1 1060.402 74.54
elaps: Routine,Calls,Wall,% = IterSCF 1 2120.809 77.15
siesta: 2 -34774.8916 -33783.9534 -33783.9671 63.0993
-2.7994 -2.7994
siesta: 3 -59491.0873 -32509.4245 -32509.4388109.3144
-76.4386 -76.4386
siesta: 4 -111533.4330 -22127.9203
-22127.9366147.8978-277.2337-277.2337
siesta: 5 -77233.2656 -25713.9823
-25714.0108127.3315-162.8721-162.8721
siesta: 6 -61186.4413 -27094.5345 -27094.5430146.7737
-54.6958 -54.6958
siesta: 7 -119320.5730 -17529.8343 -17529.8563
63.4149-285.2653-285.2653
siesta: 8 -92386.7875 -18939.3725
-18939.3917193.6104-261.1572-261.1572
....
etc... until the last iteration of the loop:
siesta: 144 -49241.0248 -23355.0125 -23355.0307 89.2656
-36.4019 -36.4019
siesta: 145 -88556.6170 -21558.7038 -21558.7083
63.0390-173.2360-173.2360
siesta: 146 -92624.4787 -17754.6189 -17754.6299
75.6454-241.8001-241.8001
siesta: 147 -63806.0777 -23578.2794 -23578.3094
69.1123-101.6597-101.6597
siesta: 148 -87399.0418 -23078.6146
-23078.6443200.3330-112.9953-112.9953
siesta: 149 -103152.0975 -15736.6015 -15736.6163
87.2093-268.4937-268.4937
Of course the forces have very large values and it is impossible to
optimize the structure. The example above is taken from BiFeO3 slabs,
but I had exactly the same problem with
1) BaTiO3 slabs (so this is not related to magnetism)
2) Pt slabs
In this example, the 4-layer slab has been perfectly converged and optimized
(supercell height 5.517471 x 2.84 ang). The SCF loop goes crazy in the
case of 6 layers (supercell height 5.517471 x 3.20 ang).
This is not a question of compilation since I have the same problem on
different machines (PC linux and SGI altix supercomputer). According
to me this is specific to systems in which there is vacuum.
Does anybody have an idea to correct the problem?
Thanks in advance!
Gregory
Here follows my complete input file:
#
# General system descriptors
#
SystemName BFO
SystemLabel BFO
NumberOfAtoms 30
NumberOfSpecies 3
%block ChemicalSpeciesLabel
1 83 Bi
2 26 Fe
3 8 O
%endblock ChemicalSpeciesLabel
#
# Basis definition
#
PAO.EnergyShift 0.01 Ry
PAO.BasisType split
block PAO.Basis < basis.fdf
%block PS.lmax
Bi 3
Fe 3
O 3
%endblock PS.lmax
#
# Lattice, coordinates, k-sampling
#
LatticeConstant 5.517471 Ang
%block LatticeVectors
1.00 0.00 0.00
0.00 1.00 0.00
0.00 0.00 3.20
%endblock LatticeVectors
AtomicCoordinatesFormat Bohr
%block AtomicCoordinatesAndAtomicSpecies
0.00000 0.00000 0.00000 1
5.21326 5.21326 0.00000 1
0.00000 0.00000 7.40297 1
5.21326 5.21326 7.40297 1
0.00000 0.00000 14.80550 1
5.21326 5.21326 14.80550 1
0.00000 5.21326 3.70126 2
0.00000 5.21326 11.10424 2
5.21326 0.00000 3.70126 2
5.21326 0.00000 11.10424 2
0.00000 5.21326 18.50676 2
5.21326 0.00000 18.50676 2
0.00000 5.21326 0.00000 3
5.21326 0.00000 0.00000 3
0.00000 5.21326 7.40297 3
5.21326 0.00000 7.40297 3
2.60663 2.60663 3.70126 3
7.81988 7.81988 3.70126 3
2.60663 7.81988 3.70126 3
7.81988 2.60663 3.70126 3
2.60663 2.60663 11.10424 3
7.81988 7.81988 11.10424 3
2.60663 7.81988 11.10424 3
7.81988 2.60663 11.10424 3
0.00000 5.21326 14.80550 3
5.21326 0.00000 14.80550 3
2.60663 2.60663 18.50676 3
7.81988 7.81988 18.50676 3
2.60663 7.81988 18.50676 3
7.81988 2.60663 18.50676 3
%endblock AtomicCoordinatesAndAtomicSpecies
# Magnetism
FixSpin .true.
TotalSpin 0.0
NonCollinearSpin .false.
%block DM.InitSpin
7 + # Atom index, spin, theta, phi (deg)
8 -
9 -
10 +
11 +
12 -
%endblock DM.InitSpin
%block kgrid_Monkhorst_Pack
4 0 0 0.0
0 4 0 0.0
0 0 1 0.0
%endblock kgrid_Monkhorst_Pack
#
# DFT, Grid, SCF
#
XC.Functional LDA
XC.Authors CA # CA = Ceperley-Alder
SpinPolarized .true.
MeshCutoff 200 Ry
Diag.DivideAndConquer .false.
#note: this one has been added after the message:
#Failure to converge standard eigenproblem Stopping Program from Node: 0
DM.MixingWeight 0.1000 # New DM amount for next SCF cycle
DM.Tolerance 1.d-5 # Tolerance in maximum difference
# between input and output DM
DM.NumberPulay 3
#DM.NumberKick 21
#DM.KickMixingWeight 0.200
MaxSCFIterations 150
#
# Eigenvalue problem: order-N or diagonalization
#
SolutionMethod diagon
ElectronicTemperature 300 K
#
# Molecular dynamics and relaxations
#
MD.TypeOfRun CG # Type of dynamics:
# - CG
# - Verlet
# - Nose
# - Parrinello-Rahman
# - Nose-Parrinello-Rahman
# - Anneal
# - FC
MD.VariableCell .false. # The lattice is relaxed together with
# the atomic coordinates?
MD.NumCGsteps 300 # Number of CG steps for
# coordinate optimization
MD.MaxCGDispl 0.5 Bohr # Maximum atomic displacement
# in one CG step
MD.MaxForceTol 0.01 eV/Ang # Tolerance in the maximum
# atomic force
MD.MaxStressTol 0.0005 eV/Ang**3
#
# Output options
#
LongOutput .false.
WriteCoorStep .true.
WriteForces .true.
WriteMullikenPop 0 # Write Mulliken Population Analysis
# - 0 = None
# - 1 = atomic and orbital charges
# - 2 = 1 + atomic overlap population
# - 3 = 2 + orbital overlap population
WriteCoorXmol .true.
WriteMDCoorXmol .false.
WriteMDhistory .false.
WriteEigenvalues .true.
AllocReportLevel 2 # Sets the level of the allocation report
# - 0 = No report at all (default)
# - 1 = Only total memory peak
and where
# it ocurred
# - 2 = detailed report printed only
# at normal program termination
# - 3 = detailed report printed at
# every new memory peak
# - 4 = print every individual
# (re)allocation or deallocation
#
# Options for saving/reading information
#
DM.UseSaveDM .false. # Use DM Continuation files
MD.UseSaveXV .false. # Use stored positions and velocities
MD.UseSaveCG .false. # Use CG history information
SaveRho .false. # Write valence pseudocharge at the mesh
SaveDeltaRho .false. # Write RHOscf-RHOatm at the mesh
SaveElectrostaticPotential .false. # Write the total elect. pot. at the mesh
SaveTotalPotential .false. # Write the total pot. at the mesh
#WriteSiestaDim .false. # Write minimum dim to siesta.h and stop
WriteDenchar .false. # Write information for DENCHAR
My basis and pseudopotentials include semicore electrons for Bi and Fe
and perfectly work for bulk R3c BiFeO3:
PAO.BasisType split
%block PAO.Basis
Bi 4
n=5 2 2 P 1
0.000 0.000
1.000 1.000
n=6 0 2 P 1
0.000 0.000
1.000 1.000
n=6 1 2 P 1
0.000 0.000
1.000 1.000
n=7 0 1
0.000
1.000
Fe 4
n=3 1 2 P 1
0.000 0.000
1.000 1.000
n=4 0 2 P 1
0.000 0.000
1.000 1.000
n=3 2 2 P 1
0.000 0.000
1.000 1.000
n=4 1 1
0.000
1.000
O 2
n=2 0 2
0.000 0.000
1.000 1.000
n=2 1 2 P 1
0.000 0.000
1.000 1.000
%endblock PAO.Basis
