I've been trying to calculate the polarization of perovskite BaTiO3 as a function of finite electric field. So, following the tutorial about polarization and electric filed, I calculated the polarization for various electric field (using a structure that I have optimized previously).
The problem I face is that the Polarization vs. Electric Filed curve comes out as a straight-line even for electric fields as high as 0.0001 a.u. Also there is no sign of hysteresis.
So I thought may be the problem is that the nuclear positions are fixed in the tutorial whereas the hysteresis / polarization switching in perovskites is due to the movement of central B atom and the O atoms under electric field.
Is there anyway to to calculate the polarization of a crystal without fixing the ions (i.e. with ionmov = 2 or 6)?
This is my input file:
Code: Select all
#Definition of the elementary cell
#*********************************
acell 7.5120216890E+00 7.5120216890E+00 7.6915729471E+00 Bohr
rprim 1 0 0
0 1 0
0 0 1
#Definition of the atoms
#***********************
natom 5
ntypat 3
znucl 56 22 08
typat 1 2 3 3 3
ixc 11
xred 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00
5.0000000000E-01 5.0000000000E-01 5.1430000000E-01
5.0000000000E-01 5.0000000000E-01 -3.0760000000E-02
5.0000000000E-01 0.0000000000E+00 4.8422377187E-01
0.0000000000E+00 5.0000000000E-01 4.8422377187E-01
#Definition of the SCF procedure
#*******************************
iscf 5
nstep 100
nband 12
nbdbuf 0
#Definition of the plane wave basis set
#**************************************
ecut 35
ecutsm 0.5
dilatmx 1.15
ngkpt 6 6 6
nshiftk 1
shiftk 0.5 0.5 0.5
toldfe 1.0d-12
tphysel 300 K
getwfk -1
ndtset 21
! ndtset 3
jdtset 11
21 22 23 24 25 # The additional 8 values of the field have been suppressed to spare CPU time
31 32 33 34 35
41 42 43 44 45
51 51 53 54 55
berryopt11 -1 rfdir11 1 1 1
berryopt21 4 efield21 0.000002 0.000002 0.000002
berryopt22 4 efield22 0.000004 0.000004 0.000004
berryopt23 4 efield23 0.000006 0.000006 0.000006
berryopt24 4 efield24 0.000008 0.000008 0.000008
berryopt25 4 efield25 0.000010 0.000010 0.000010
berryopt31 4 efield31 0.000008 0.000008 0.000008
berryopt32 4 efield32 0.000006 0.000006 0.000006
berryopt33 4 efield33 0.000004 0.000004 0.000004
berryopt34 4 efield34 0.000002 0.000002 0.000002
berryopt35 4 efield35 0.000 0.000 0.000
P.S. Another question has been bothering me. When I do the structural optimization run, the energy converges to a very high value (about 50 Ha). My input file for optimization:
Code: Select all
ndtset 2 # There are 2 datasets in this calculation
# Set 1 : Internal coordinate optimization
ionmov1 2 # Use BFGS algorithm for structural optimization
ntime1 7 # Maximum number of optimization steps
tolmxf1 1.0e-5 # Optimization is converged when maximum force
# (Hartree/Bohr) is less than this maximum
natfix1 3 # Fix the position of two symmetry-equivalent atoms
# in doing the structural optimization
iatfix1 1 2 3 # Choose atoms 1 and 2 as the fixed atoms (see discussion)
optcell1 0
# Set 2 : Lattice parameter relaxation (including re-optimization of
# internal coordinates)
dilatmx2 1.15 # Maximum scaling allowed for lattice parameters
getxred2 -1 # Start with relaxed coordinates from dataset 1
getwfk2 -1 # Start with wave functions from dataset 1
ionmov2 2 # Use BFGS algorithm
ntime2 10 # Maximum number of optimization steps
optcell2 3 # Fully optimize unit cell geometry, keeping symmetry
tolmxf2 1.0e-5 # Convergence limit for forces as above
strfact2 100 # Test convergence of stresses (Hartree/bohr^3) by
# multiplying by this factor and applying force
# convergence test
natfix2 3
iatfix2 1 2 3
#Common input data
acell 3.992 3.992 4.036 angstrom #this is a guess, with the c/a
#ratio based on ideal tetrahedral
#bond angles
rprim 1 0 0 # primitive vectors must be
0 1 0 #specified with high accuracy to be
0 0 1 #sure that the symmetry is recognized
#and preserved in the optimization
#process
#Definition of the atom types and atoms
ntypat 3
znucl 56 22 08
natom 5
typat 1 2 3 3 3
#Starting approximation for atomic positions in REDUCED coordinates
#based on ideal tetrahedral bond angles
xred 0.0 0.0 0.0
0.5 0.5 0.5143
0.5 0.5 -0.03076
0.5 0.0 0.4814
0.0 0.5 0.4814
#Gives the number of bands, explicitely (do not take the default)
nband 15 # For an insulator (if described correctly as an
# insulator by DFT), conduction bands should not
# be included in response-function calculations
#Definition of the plane wave basis set
ecut 50 # Maximum kinetic energy cutoff (Hartree)
ecutsm 0.5 # Smoothing energy needed for lattice paramete
# optimization. This will be retained for
# consistency throughout.
ixc = 11 !GGA, Perdew-Burke-Ernzerhof GGA functional
tphysel 300 K
#Definition of the k-point grid
ngkpt 6 6 6 # 6 6 6 grid
nshiftk 1 # Use one copy of grid only (default)
shiftk 0.5 0.5 0.5 # This choice of origin for the k point grid
# preserves the hexagonal symmetry of the grid,
# which would be broken by the default choice.
#Definition of the self-consistency procedure
diemac 4000 # Model dielectric preconditioner
iscf 7 # Pulay mixing of the potential
nstep 40 # Maxiumum number of SCF iterations
tolvrs 1.0d-18 # Strict tolerance on (squared) residual of the
# SCF potential needed for accurate forces and
# stresses in the structural optimization, and
# accurate wave functions in the RF calculations
The etotal and fcart are:
Code: Select all
etotal1 -5.3094090418E+01
etotal2 -5.3094307091E+01
fcart1 -0.0000000000E+00 -0.0000000000E+00 -6.9446291206E-03
-0.0000000000E+00 -0.0000000000E+00 -2.4101663245E-02
-0.0000000000E+00 -0.0000000000E+00 3.1042348651E-02
-0.0000000000E+00 -0.0000000000E+00 1.9718573189E-06
-0.0000000000E+00 -0.0000000000E+00 1.9718573189E-06
fcart2 -0.0000000000E+00 -0.0000000000E+00 -6.7464332558E-03
-0.0000000000E+00 -0.0000000000E+00 -1.9553866581E-02
-0.0000000000E+00 -0.0000000000E+00 2.6297646773E-02
-0.0000000000E+00 -0.0000000000E+00 1.3265321034E-06
-0.0000000000E+00 -0.0000000000E+00 1.3265321034E-06
. I'll add the appropriate commands to my input and let you know as soon as I get results. Since I'm not doing calculations in parallel, that may take a while.