-------- Forwarded Message -------- Subject: Problem with electric polarization calculation in quantum espresso Date: Thu, 14 Mar 2024 01:05:32 +0000 From: Akshay Mahajan <akshaymah...@iisc.ac.in> To: users@lists.quantum-espresso.org <users@lists.quantum-espresso.org>
Dear QE Users,I have been trying to calculate electric polarization for a 2D Bi monolayer (along the in-plane x-direction). I am using QE version 6.4.1.There are two issues that I am facing,
1.
The magnitude of the polarization changes significantly when the k_y
values change (I am keeping k_z to 1 only since it is a 2D system).
I expected it to reach some converged value but rather it keeps
increasing on increasing the k_y value and then drops and then again
starts increasing. The bash file is included in the email for this
calculation. In the bash file (polarization_conv.sh) I am changing
both k_x and k_y in the scf calculation and changing nppstr to 2*k_x
and keeping k_x =nppstr in the nscf calculation. The result from
running this bash script run is also included(pol_conv_Bi_mono). One
can modify the script and check that only changing the k_y is
causing this huge change in the polarization values.
2.
As per this paper
(https://onlinelibrary.wiley.com/doi/full/10.1002/adfm.201707383
<https://onlinelibrary.wiley.com/doi/full/10.1002/adfm.201707383>)
the positive delta_h (refer to the paper for definition) should
result in a positive polarization. While trying to calculate the
variation of electric polarization with delta_h, I am getting
negative polarization values for positive delta_h near the
centrosymmetric structure (delta_h=0) and positive values for
negative delta_h (again near the centrosymmetric structure
(delta_h=0)). The bash script (polarization.sh) for this calculation
and the result(pol_data_Bi_mono) are also attached.
I am also attaching the pseudopotential file for reference. Please let
me know if anyone can determine the issue behind these incorrect results.
Thank you. Regards, *Akshay Mahajan* *Ph.D. Student * *Prime Minister Research Fellow (PMRF) * *Solid State and Structural Chemistry Unit (SSCU) * *Indian Institute of Science (IISc)* *Bengaluru, India *
pol_conv_Bi_mono
Description: Binary data
pol_data_Bi_mono
Description: Binary data
#!/bin/bash
# Script to get the energy versus polarization and energy versus displacement plots
# Exclusive for Bi project (uses verticle dispalcement for polarization)
######### Energy conversion from Ry to eV parameter #########
Q=13.605703976
######### Intial values from relaxed FE structure #########
h_u_pol=18.419321249 # polar h displacement in upper layer
h_l_pol=15.420625873 # polar h displacement in lower layer
u_ref=17.790480876 # ref for polar h displacement in upper layer
l_ref=14.791753001 # ref for polar h displacement in lower layer
z_vacc=33.211090088 # Lattice parameter along the vaccum
######### Initializing the delta_h values for polar structure #########
delta_h_u_pol=`echo "scale=9; $h_u_pol - $u_ref" | bc`
delta_h_l_pol=`echo "scale=9; $h_l_pol - $l_ref" | bc`
######### Initializing the delta_h value of equally stable other polar structure #########
neg=-1
neg_delta_h_u_pol=`echo "scale=9; $delta_h_u_pol * $neg" | bc`
neg_delta_h_l_pol=`echo "scale=9; $delta_h_l_pol * $neg" | bc`
######### Intializing the outskirt of polar displacement for plotting W curves #########
## positive outskirt
#u_pos_max=`echo "scale=9; 2.5*$delta_h_u_pol" | bc`
#l_pos_max=`echo "scale=9; 2.5*$delta_h_l_pol" | bc`
## negative outskirt
#u_neg_max=`echo "scale=9; 2.5*$neg_delta_h_u_pol" | bc`
#l_neg_max=`echo "scale=9; 2.5*$neg_delta_h_l_pol" | bc`
######### Creating delta_h values for upper layer #########
#delta_h_u_list=($(seq $u_pos_max -0.15 $delta_h_u_pol))
#delta_h_u_list+=($(seq $delta_h_u_pol -0.05 0))
delta_u=`echo "scale=9; $delta_h_u_pol/10" | bc`
delta_u_lst=`echo "scale=9; $delta_u * $neg" | bc`
delta_h_u_list=($(seq $delta_h_u_pol $delta_u_lst 0))
## getting the negative displacements
u_neg_delta_list=()
for (( idx=${#delta_h_u_list[@]}-1 ; idx>=0 ; idx-- )); do
u_neg_delta_list+=(`echo "scale=9; ${delta_h_u_list[$idx]} * $neg" | bc`)
done
## appending the negative displacements
delta_h_u_list+=(0)
delta_h_u_list+=(${u_neg_delta_list[@]})
######### Creating delta_h values for lower layer #########
#delta_h_l_list=($(seq $l_pos_max -0.15 $delta_h_l_pol))
#delta_h_l_list+=($(seq $delta_h_l_pol -0.05 0))
delta_l=`echo "scale=9; $delta_h_l_pol/10" | bc`
delta_l_lst=`echo "scale=9; $delta_l * $neg" | bc`
delta_h_l_list=($(seq $delta_h_l_pol $delta_l_lst 0))
## getting the negative displacements
l_neg_delta_list=()
for (( idx=${#delta_h_l_list[@]}-1 ; idx>=0 ; idx-- )); do
l_neg_delta_list+=(`echo "scale=9; ${delta_h_l_list[$idx]} * $neg" | bc`)
done
## appending the negative dispalcements
delta_h_l_list+=(0)
delta_h_l_list+=(${l_neg_delta_list[@]})
######### Normalized Displacement upper part #########
normalized_disp=()
for (( idx=0 ; idx<=${#delta_h_u_list[@]}-1 ; idx++ )); do
normalized_disp+=(`echo "scale=9; ${delta_h_u_list[$idx]} / $delta_h_u_pol" | bc`)
done
######### Normalized Displacement lower part #########
normalized_disp_2=()
for (( idx=0 ; idx<=${#delta_h_l_list[@]}-1 ; idx++ )); do
normalized_disp_2+=(`echo "scale=9; ${delta_h_l_list[$idx]} / $delta_h_l_pol" | bc`)
done
######### Sweeping across all values of delta_h for W curve calcualtions #########
#echo "Normalized Displacement Energy(meV/atom) Polarization (10^(-10) C/m)" >> pol_data_Bi_mono
echo "Delta_h_u = $delta_h_u_pol Delta_h_l = $delta_h_l_pol" >> normalized_disp_data
echo "Normalized displacment for both layers of monolayer are provided below" >> normalized_disp_data
echo "Upper (layer-1) Lower (layer-1)" >> normalized_disp_data
polarization=()
tot_en=()
for (( idx=0 ; idx<=${#normalized_disp[@]}-1 ; idx++ )); do
scf_name="Bi_BP_mono"
nscf_name="Bi_BP_mono"
#mkdir folder-${normalized_disp[$idx]}
mkdir folder-$idx
cp main_job.sh folder-$idx
cd folder-$idx
coor_u_pol=`echo "scale=9; $u_ref + ${delta_h_u_list[$idx]}" | bc`
coor_l_pol=`echo "scale=9; $l_ref + ${delta_h_l_list[$idx]}" | bc`
# Create the input file XXXX.scf.in.
cat > $scf_name.scf.in << EOF
&control
calculation = 'scf'
restart_mode = 'from_scratch'
prefix = 'bi_bp'
pseudo_dir = '/home/akshay/pseudo/pbe.0.3.1/PSEUDOPOTENTIALS'
outdir = './tmp/'
/
&system
ibrav = 0
nat = 4
ntyp = 1
ecutwfc = 74.0
ecutrho = 592.0
vdw_corr = 'dft-d'
/
&electrons
diagonalization = 'david'
mixing_mode = 'plain'
mixing_beta = 0.7
conv_thr = 7.35d-10
/
ATOMIC_SPECIES
Bi 208.9804 Bi.pbe-dn-rrkjus_psl.0.2.2.UPF
CELL_PARAMETERS {angstrom}
4.898252066 0.000000000 0.000000000
0.000000000 4.443566326 0.000000000
0.000000000 0.000000000 33.211090088
ATOMIC_POSITIONS {angstrom}
Bi 0.625057731 3.332674564 14.791753001
Bi 3.074182872 1.110891844 $coor_u_pol
Bi 0.218643999 3.332674564 17.790480876
Bi 2.667804736 1.110891844 $coor_l_pol
K_POINTS {automatic}
16 16 1 0 0 0
EOF
# Create the input file XXXX.scf.in.
cat > $nscf_name.pol.in << EOF
&control
calculation = 'nscf'
restart_mode = 'from_scratch'
prefix = 'bi_bp'
pseudo_dir = '/home/akshay/pseudo/pbe.0.3.1/PSEUDOPOTENTIALS'
outdir = './tmp/'
lberry = .true.
gdir = 1
nppstr = 32
/
&system
ibrav = 0
nat = 4
ntyp = 1
ecutwfc = 74.0
ecutrho = 592.0
vdw_corr = 'dft-d'
/
&electrons
diagonalization = 'david'
mixing_mode = 'plain'
mixing_beta = 0.7
conv_thr = 7.35d-10
/
ATOMIC_SPECIES
Bi 208.9804 Bi.pbe-dn-rrkjus_psl.0.2.2.UPF
CELL_PARAMETERS {angstrom}
4.898252066 0.000000000 0.000000000
0.000000000 4.443566326 0.000000000
0.000000000 0.000000000 33.211090088
ATOMIC_POSITIONS {angstrom}
Bi 0.625057731 3.332674564 14.791753001
Bi 3.074182872 1.110891844 $coor_u_pol
Bi 0.218643999 3.332674564 17.790480876
Bi 2.667804736 1.110891844 $coor_l_pol
K_POINTS {automatic}
32 16 1 0 0 0
EOF
# Running the files
sed -i -e "2s/Bi_BP/Bi_pol_$idx/g" main_job.sh
qsub main_job.sh
stop=0
while [ $stop -eq 0 ]
do
grep -q 'JOB DONE.' $nscf_name.pol.out && stop=1
done
#rm -r ./tmp
cd ../ # Go back to the directory from where this script is executed.
ene=`grep "! total energy" folder-$idx/${scf_name}.scf.out | awk '{printf "%15.9f\n",$5}'`
tot_en+=(`grep "! total energy" folder-$idx/${scf_name}.scf.out | awk '{printf "%15.9f\n",$5}'`)
pol=`grep "C/m^2" folder-$idx/${nscf_name}.pol.out | awk '{printf "%5.7f\n",$3}'`
polarization+=(`echo "scale=15; $pol * $z_vacc" | bc`)
echo "${normalized_disp[$idx]} $ene $pol" >> pol_data_Bi_mono
# The >> symbols indicates that the data in each run must be appended to the same file and not overwrite the file.
echo "${normalized_disp[$idx]} ${tot_en[$idx]} ${polarization[$idx]}" >> run_check
echo "${normalized_disp[$idx]} ${normalized_disp_2[$idx]}" >> normalized_disp_data
done
######### Post-processing polarization values and storing the values #########
#para_idx=`echo "scale=0; ( (${#normalized_disp[@]}-1) / 2 )" | bc`
#for (( idx=0 ; idx<=${#normalized_disp[@]}-1 ; idx++ )); do
# energy_final=`echo "scale=9; ( (${tot_en[$idx]} * $Q * 1000) / 4 )" | bc`
# polarization_final=`echo "scale=9; ${polarization[$idx]} - ${polarization[$para_idx]}" | bc`
# echo "${normalized_disp[$idx]} $energy_final $polarization_final" >> pol_data_Bi_mono
#done
#!/bin/bash
# Script to get the energy versus polarization and energy versus displacement plots
# Exclusive for Bi project (uses verticle dispalcement for polarization)
z_vacc=33.211090088
#k=16
#pol_polar=`grep "C/m^2" folder-16/Bi_BP_polar.pol.out | awk '{printf "%5.7f\n",$3}'`
#pol_non_polar=`grep "C/m^2" folder-16/Bi_BP_non_polar.pol.out | awk '{printf "%5.7f\n",$3}'`
#polar_diff=`echo "scale=15; $pol_polar - $pol_non_polar" | bc`
#polarization=`echo "scale=15; $polar_diff * $z_vacc" | bc`
#echo "$k $polarization" >> pol_k_conv_Bi_mono
echo "Checking polarization convergence with k (npppstr=2*k)" >> pol_conv_Bi_mono
for k in $(seq 33 1 50)
do
nscf_pol_name="Bi_BP_mono.pol"
scf_pol_name="Bi_BP_mono.scf"
mkdir folder-$k
cp main_job.sh folder-$k
cd folder-$k
nppstr=`echo "scale=15; $k * 2" | bc`
# Create the input file XXXX.scf.in.
cat > $scf_pol_name.in << EOF
&control
calculation = 'scf'
restart_mode = 'from_scratch'
prefix = 'bi_bp'
pseudo_dir = '/home/akshay/pseudo/pbe.0.3.1/PSEUDOPOTENTIALS'
outdir = './tmp/'
/
&system
ibrav = 0
nat = 4
ntyp = 1
ecutwfc = 74.0
ecutrho = 592.0
vdw_corr = 'dft-d'
/
&electrons
diagonalization = 'david'
mixing_mode = 'plain'
mixing_beta = 0.7
conv_thr = 7.35d-10
/
ATOMIC_SPECIES
Bi 208.9804 Bi.pbe-dn-rrkjus_psl.0.2.2.UPF
CELL_PARAMETERS {angstrom}
4.898252066 0.000000000 0.000000000
0.000000000 4.443566326 0.000000000
0.000000000 0.000000000 33.211090088
ATOMIC_POSITIONS {angstrom}
Bi 0.625057731 3.332674564 14.791753001
Bi 3.074182872 1.110891844 18.419321249
Bi 0.218643999 3.332674564 17.790480876
Bi 2.667804736 1.110891844 15.420625873
K_POINTS {automatic}
$k $k 1 0 0 0
EOF
# Create the input file XXXX.nscf.in.
cat > $nscf_pol_name.in << EOF
&control
calculation = 'nscf'
restart_mode = 'from_scratch'
prefix = 'bi_bp'
pseudo_dir = '/home/akshay/pseudo/pbe.0.3.1/PSEUDOPOTENTIALS'
outdir = './tmp/'
lberry = .true.
gdir = 1
nppstr = $nppstr
/
&system
ibrav = 0
nat = 4
ntyp = 1
ecutwfc = 74.0
ecutrho = 592.0
vdw_corr = 'dft-d'
/
&electrons
diagonalization = 'david'
mixing_mode = 'plain'
mixing_beta = 0.7
conv_thr = 7.35d-10
/
ATOMIC_SPECIES
Bi 208.9804 Bi.pbe-dn-rrkjus_psl.0.2.2.UPF
CELL_PARAMETERS {angstrom}
4.898252066 0.000000000 0.000000000
0.000000000 4.443566326 0.000000000
0.000000000 0.000000000 33.211090088
ATOMIC_POSITIONS {angstrom}
Bi 0.625057731 3.332674564 14.791753001
Bi 3.074182872 1.110891844 18.419321249
Bi 0.218643999 3.332674564 17.790480876
Bi 2.667804736 1.110891844 15.420625873
K_POINTS {automatic}
$nppstr $k 1 0 0 0
EOF
# Running the files
sed -i -e "2s/Bi_BP/pol_$k/g" main_job.sh
# Starting both runs
qsub main_job.sh
# waiting for polar run to finish
stop=0
while [ $stop -eq 0 ]
do
grep -q 'JOB DONE.' $nscf_pol_name.out && stop=1
done
cd ../ # Go back to the directory from where this script is executed.
# Getting polarization values
pol_val=`grep "C/m^2" folder-${k}/${nscf_pol_name}.out | awk '{printf "%5.7f\n",$3}'`
polarization=`echo "scale=15; $pol_val * $z_vacc" | bc`
echo "$k $polarization" >> pol_conv_Bi_mono
done
_______________________________________________ The Quantum ESPRESSO community stands by the Ukrainian people and expresses its concerns about the devastating effects that the Russian military offensive has on their country and on the free and peaceful scientific, cultural, and economic cooperation amongst peoples _______________________________________________ Quantum ESPRESSO is supported by MaX (www.max-centre.eu) users mailing list users@lists.quantum-espresso.org https://lists.quantum-espresso.org/mailman/listinfo/users
