mirror of https://gitlab.com/QEF/q-e.git
*** empty log message ***
git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@226 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
parent
a85db0682b
commit
d6193c5491
10
PW/Makefile
10
PW/Makefile
|
@ -21,6 +21,16 @@ atomic_wfc.o \
|
||||||
bachel.o \
|
bachel.o \
|
||||||
becmod.o \
|
becmod.o \
|
||||||
bfgs.o \
|
bfgs.o \
|
||||||
|
bp_bess.o \
|
||||||
|
bp_calc_btq.o \
|
||||||
|
bp_c_phase.o \
|
||||||
|
bp_dbess.o \
|
||||||
|
bp_qvan3.o \
|
||||||
|
bp_radin.o \
|
||||||
|
bp_strings.o \
|
||||||
|
bp_ylm_q.o \
|
||||||
|
bp_zgedi.o \
|
||||||
|
bp_zgefa.o \
|
||||||
broadcast.o \
|
broadcast.o \
|
||||||
c_bands.o \
|
c_bands.o \
|
||||||
ccalbec.o \
|
ccalbec.o \
|
||||||
|
|
|
@ -0,0 +1,81 @@
|
||||||
|
C-------------------------------------------------------------------------
|
||||||
|
SUBROUTINE BESS(XG,L,MMAX,R,JL)
|
||||||
|
C-------------------------------------------------------------------------
|
||||||
|
C CALCULATES SPHERICAL BESSEL FUNCTIONS j_l(Gr)
|
||||||
|
C
|
||||||
|
IMPLICIT REAL*8 (A-H,O-Z)
|
||||||
|
PARAMETER(EPS=1.E-8)
|
||||||
|
REAL*8 JL(MMAX),R(MMAX)
|
||||||
|
|
||||||
|
IF(L.EQ.1) THEN ! S PART
|
||||||
|
IF(XG.LT.EPS) THEN
|
||||||
|
DO 41 IR=1,MMAX
|
||||||
|
41 JL(IR)=1.D0
|
||||||
|
ELSE
|
||||||
|
JL(1)=1.D0
|
||||||
|
DO 42 IR=2,MMAX
|
||||||
|
XRG=R(IR)*XG
|
||||||
|
JL(IR)=SIN(XRG)/XRG
|
||||||
|
42 CONTINUE
|
||||||
|
ENDIF
|
||||||
|
ENDIF
|
||||||
|
|
||||||
|
IF(L.EQ.2) THEN ! P PART
|
||||||
|
IF(XG.LT.EPS) THEN
|
||||||
|
DO 43 IR=1,MMAX
|
||||||
|
43 JL(IR)=0.D0
|
||||||
|
ELSE
|
||||||
|
JL(1)=0.
|
||||||
|
DO 44 IR=2,MMAX
|
||||||
|
XRG=R(IR)*XG
|
||||||
|
JL(IR)=(SIN(XRG)/XRG-COS(XRG))/XRG
|
||||||
|
44 CONTINUE
|
||||||
|
ENDIF
|
||||||
|
ENDIF
|
||||||
|
|
||||||
|
IF(L.EQ.3) THEN ! D PART
|
||||||
|
IF(XG.LT.EPS) THEN
|
||||||
|
DO 45 IR=1,MMAX
|
||||||
|
45 JL(IR)=0.D0
|
||||||
|
ELSE
|
||||||
|
JL(1)=0.D0
|
||||||
|
DO 46 IR=2,MMAX
|
||||||
|
XRG=R(IR)*XG
|
||||||
|
JL(IR)=(SIN(XRG)*(3./(XRG*XRG)-1.)
|
||||||
|
+ -3.*COS(XRG)/XRG) /XRG
|
||||||
|
46 CONTINUE
|
||||||
|
ENDIF
|
||||||
|
ENDIF
|
||||||
|
|
||||||
|
IF(L.EQ.4) THEN ! F PART
|
||||||
|
IF(XG.LT.EPS) THEN
|
||||||
|
DO 47 IR=1,MMAX
|
||||||
|
47 JL(IR)=0.D0
|
||||||
|
ELSE
|
||||||
|
JL(1)=0.D0
|
||||||
|
DO 48 IR=2,MMAX
|
||||||
|
XRG=R(IR)*XG
|
||||||
|
XRG2=XRG*XRG
|
||||||
|
JL(IR)=( SIN(XRG)*(15./(XRG2*XRG)-6./XRG)
|
||||||
|
+ +COS(XRG)*(1.-15./XRG2) )/XRG
|
||||||
|
48 CONTINUE
|
||||||
|
ENDIF
|
||||||
|
ENDIF
|
||||||
|
|
||||||
|
IF(L.EQ.5) THEN ! G PART
|
||||||
|
IF(XG.LT.EPS) THEN
|
||||||
|
DO 49 IR=1,MMAX
|
||||||
|
49 JL(IR)=0.D0
|
||||||
|
ELSE
|
||||||
|
JL(1)=0.D0
|
||||||
|
DO 50 IR=2,MMAX
|
||||||
|
XRG=R(IR)*XG
|
||||||
|
XRG2=XRG*XRG
|
||||||
|
JL(IR)=( SIN(XRG)*(105./(XRG2*XRG2)-45./XRG2+1.)
|
||||||
|
+ +COS(XRG)*(10./XRG-105./(XRG2*XRG)) )/XRG
|
||||||
|
50 CONTINUE
|
||||||
|
ENDIF
|
||||||
|
ENDIF
|
||||||
|
|
||||||
|
RETURN
|
||||||
|
END
|
|
@ -0,0 +1,823 @@
|
||||||
|
!##############################################################################!
|
||||||
|
!# #!
|
||||||
|
!# #!
|
||||||
|
!# This is the main one of a set of Fortran 90 files designed to compute #!
|
||||||
|
!# the electrical polarization in a crystaline solid. #!
|
||||||
|
!# #!
|
||||||
|
!# #!
|
||||||
|
!# AUTHORS #!
|
||||||
|
!# ~~~~~~~ #!
|
||||||
|
!# This set of subprograms is based on code written in an early Fortran #!
|
||||||
|
!# 77 version of PWSCF by Alessio Filippetti. These were later ported #!
|
||||||
|
!# into another version by Lixin He. Oswaldo Dieguez, in collaboration #!
|
||||||
|
!# with Lixin He and Jeff Neaton, ported these routines into Fortran 90 #!
|
||||||
|
!# version 1.2.1 of PWSCF. He, Dieguez, and Neaton were working at the #!
|
||||||
|
!# time in David Vanderbilt's group at Rutgers, The State University of #!
|
||||||
|
!# New Jersey, USA. #!
|
||||||
|
!# #!
|
||||||
|
!# #!
|
||||||
|
!# LIST OF FILES #!
|
||||||
|
!# ~~~~~~~~~~~~~ #!
|
||||||
|
!# The complete list of files added to the PWSCF distribution is: #!
|
||||||
|
!# * ../PW/bp_bess.f #!
|
||||||
|
!# * ../PW/bp_calc_btq.f90 #!
|
||||||
|
!# * ../PW/bp_c_phase.f90 #!
|
||||||
|
!# * ../PW/bp_dbess.f #!
|
||||||
|
!# * ../PW/bp_qvan3.f90 #!
|
||||||
|
!# * ../PW/bp_radin.f #!
|
||||||
|
!# * ../PW/bp_strings.f90 #!
|
||||||
|
!# * ../PW/bp_ylm_q.f #!
|
||||||
|
!# * ../PW/bp_zgedi.f #!
|
||||||
|
!# * ../PW/bp_zgefa.f #!
|
||||||
|
!# #!
|
||||||
|
!# The PWSCF files that needed (minor) modifications were: #!
|
||||||
|
!# * ../PW/electrons.f90 #!
|
||||||
|
!# * ../PW/input.f90 #!
|
||||||
|
!# * ../PW/pwcom.f90 #!
|
||||||
|
!# * ../PW/setup.f90 #!
|
||||||
|
!# #!
|
||||||
|
!# #!
|
||||||
|
!# BRIEF SUMMARY OF THE METHODOLOGY #!
|
||||||
|
!# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #!
|
||||||
|
!# The spontaneous polarization has two contibutions, electronic #!
|
||||||
|
!# and ionic. With these additional routines, PWSCF will output both. #!
|
||||||
|
!# #!
|
||||||
|
!# The ionic contribution is relatively trivial to compute, requiring #!
|
||||||
|
!# knowledge only of the atomic positions and core charges. The new #!
|
||||||
|
!# subroutines focus mainly on evaluating the electronic contribution, #!
|
||||||
|
!# computed as a Berry phase, i.e., a global phase property that can #!
|
||||||
|
!# be computed from inner products of Bloch states at neighboring #!
|
||||||
|
!# points in k-space. #!
|
||||||
|
!# #!
|
||||||
|
!# The standard procedure would be for the user to first perform a #!
|
||||||
|
!# self-consistent (sc) calculation to obtain a converged charge density. #!
|
||||||
|
!# With well-converged sc charge density, the user would then run one #!
|
||||||
|
!# or more non-self consistent (or "band structure") calculations, #!
|
||||||
|
!# using the same main code, but with a flag to ask for the polarization. #!
|
||||||
|
!# Each such run would calculate the projection of the polarization #!
|
||||||
|
!# onto one of the three primitive reciprocal lattice vectors. In #!
|
||||||
|
!# cases of high symmetry (e.g. a tetragonal ferroelectric phase), one #!
|
||||||
|
!# such run would suffice. In the general case of low symmetry, the #!
|
||||||
|
!# user would have to submit up to three jobs to compute the three #!
|
||||||
|
!# components of polarization, and would have to obtain the total #!
|
||||||
|
!# polarization "by hand" by summing these contributions. #!
|
||||||
|
!# #!
|
||||||
|
!# Accurate calculation of the electronic or "Berry-phase" polarization #!
|
||||||
|
!# requires overlaps between wavefunctions along fairly dense lines (or #!
|
||||||
|
!# "strings") in k-space in the direction of the primitive G-vector for #!
|
||||||
|
!# which one is calculating the projection of the polarization. The #!
|
||||||
|
!# code would use a higher-density k-mesh in this direction, and a #!
|
||||||
|
!# standard-density mesh in the two other directions. See below for #!
|
||||||
|
!# details. #!
|
||||||
|
!# #!
|
||||||
|
!# #!
|
||||||
|
!# FUNCTIONALITY/COMPATIBILITY #!
|
||||||
|
!# ~~~~~~~~~~~~~~~~~~~~~~~~~~~ #!
|
||||||
|
!# * Berry phases for a given G-vector. #!
|
||||||
|
!# #!
|
||||||
|
!# * Contribution to the polarization (in relevant units) for a given #!
|
||||||
|
!# G-vector. #!
|
||||||
|
!# #!
|
||||||
|
!# * Spin-polarized systems supported. #!
|
||||||
|
!# #!
|
||||||
|
!# * Ultrasoft and norm-conserving pseudopotentials supported. #!
|
||||||
|
!# #!
|
||||||
|
!# * The value of the "polarization quantum" and the ionic contribution #!
|
||||||
|
!# to the polarization are reported. #!
|
||||||
|
!# #!
|
||||||
|
!# #!
|
||||||
|
!# NEW INPUT PARAMETERS #!
|
||||||
|
!# ~~~~~~~~~~~~~~~~~~~~ #!
|
||||||
|
!# * lberry (.TRUE. or .FALSE.) #!
|
||||||
|
!# Tells PWSCF that a Berry phase calcultion is desired. #!
|
||||||
|
!# #!
|
||||||
|
!# * gdir (1, 2, or 3) #!
|
||||||
|
!# Specifies the direction of the k-point strings in reciprocal space. #!
|
||||||
|
!# '1' refers to the first reciprocal lattice vector, '2' to the #!
|
||||||
|
!# second, and '3' to the third. #!
|
||||||
|
!# #!
|
||||||
|
!# * nppstr (integer) #!
|
||||||
|
!# Specifies the number of k-points to be calculated along each #!
|
||||||
|
!# symmetry-reduced string. #!
|
||||||
|
!# #!
|
||||||
|
!# #!
|
||||||
|
!# EXPLANATION OF K-POINT MESH #!
|
||||||
|
!# ~~~~~~~~~~~~~~~~~~~~~~~~~~~ #!
|
||||||
|
!# If gdir=1, the program takes the standard input specification of the #!
|
||||||
|
!# k-point mesh (nk1 x nk2 x nk3) and stops if the k-points in dimension #!
|
||||||
|
!# 1 are not equally spaced or if its number is not equal to nppstr, #!
|
||||||
|
!# working with a mesh of dimensions (nppstr x nk2 x nk3). That is, for #!
|
||||||
|
!# each point of the (nk2 x nk3) two-dimensional mesh, it works with a #!
|
||||||
|
!# string of nppstr k-points extending in the third direction. Symmetry #!
|
||||||
|
!# will be used to reduce the number of strings (and assign them weights) #!
|
||||||
|
!# if possible. Of course, if gdir=2 or 3, the variables nk2 or nk3 will #!
|
||||||
|
!# be overridden instead, and the strings constructed in those #!
|
||||||
|
!# directions, respectively. #!
|
||||||
|
!# #!
|
||||||
|
!# #!
|
||||||
|
!# BIBLIOGRAPHY #!
|
||||||
|
!# ~~~~~~~~~~~~ #!
|
||||||
|
!# The theory behind this implementation is described in: #!
|
||||||
|
!# [1] R D King-Smith and D Vanderbilt, "Theory of polarization of #!
|
||||||
|
!# crystaline solids", Phys Rev B 47, 1651 (1993). #!
|
||||||
|
!# #!
|
||||||
|
!# Other relevant sources of information are: #!
|
||||||
|
!# [2] D Vanderbilt and R D King-Smith, "Electronic polarization in the #!
|
||||||
|
!# ultrasoft pseudopotential formalism", internal report (1998), #!
|
||||||
|
!# [3] D Vanderbilt, "Berry phase theory of proper piezoelectric #!
|
||||||
|
!# response", J Phys Chem Solids 61, 147 (2000). #!
|
||||||
|
!# #!
|
||||||
|
!# #!
|
||||||
|
!# dieguez@physics.rutgers.edu #!
|
||||||
|
!# 09 June 2003 #!
|
||||||
|
!# #!
|
||||||
|
!# #!
|
||||||
|
!##############################################################################!
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
!======================================================================!
|
||||||
|
|
||||||
|
SUBROUTINE c_phase
|
||||||
|
|
||||||
|
!----------------------------------------------------------------------!
|
||||||
|
|
||||||
|
! Geometric phase calculation along a strip of nppstr k-points
|
||||||
|
! averaged over a 2D grid of nkort k-points ortogonal to nppstr
|
||||||
|
|
||||||
|
#include "machine.h"
|
||||||
|
|
||||||
|
! --- Make use of the module with common information ---
|
||||||
|
USE pwcom
|
||||||
|
|
||||||
|
! --- Avoid implicit definitions ---
|
||||||
|
IMPLICIT NONE
|
||||||
|
|
||||||
|
! --- Internal definitions ---
|
||||||
|
INTEGER :: i
|
||||||
|
INTEGER :: igk1(npwx)
|
||||||
|
INTEGER :: igk0(npwx)
|
||||||
|
INTEGER :: ik
|
||||||
|
INTEGER :: ind1
|
||||||
|
INTEGER :: info
|
||||||
|
INTEGER :: is
|
||||||
|
INTEGER :: istring
|
||||||
|
INTEGER :: iv
|
||||||
|
INTEGER :: ivpt(nbnd)
|
||||||
|
INTEGER :: j
|
||||||
|
INTEGER :: jkb
|
||||||
|
INTEGER :: jkb_bp
|
||||||
|
INTEGER :: jkb1
|
||||||
|
INTEGER :: job
|
||||||
|
INTEGER :: jv
|
||||||
|
INTEGER :: kindex
|
||||||
|
INTEGER :: kort
|
||||||
|
INTEGER :: kpar
|
||||||
|
INTEGER :: kpoint
|
||||||
|
INTEGER :: kstart
|
||||||
|
INTEGER :: mb
|
||||||
|
INTEGER :: mk1
|
||||||
|
INTEGER :: mk2
|
||||||
|
INTEGER :: mk3
|
||||||
|
INTEGER , ALLOCATABLE :: mod_elec(:)
|
||||||
|
INTEGER :: mod_elec_dw
|
||||||
|
INTEGER :: mod_elec_tot
|
||||||
|
INTEGER :: mod_elec_up
|
||||||
|
INTEGER :: mod_ion(nat)
|
||||||
|
INTEGER :: mod_ion_dw
|
||||||
|
INTEGER :: mod_ion_tot
|
||||||
|
INTEGER :: mod_ion_up
|
||||||
|
INTEGER :: mod_tot
|
||||||
|
INTEGER :: n1
|
||||||
|
INTEGER :: n2
|
||||||
|
INTEGER :: n3
|
||||||
|
INTEGER :: na
|
||||||
|
INTEGER :: nb
|
||||||
|
INTEGER :: ng
|
||||||
|
INTEGER :: nhjkb
|
||||||
|
INTEGER :: nhjkbm
|
||||||
|
INTEGER :: nkbtona(nkb)
|
||||||
|
INTEGER :: nkbtonh(nkb)
|
||||||
|
INTEGER :: nkort
|
||||||
|
INTEGER :: np
|
||||||
|
INTEGER :: npw1
|
||||||
|
INTEGER :: npw0
|
||||||
|
INTEGER :: nstring
|
||||||
|
INTEGER :: nt
|
||||||
|
LOGICAL :: lodd
|
||||||
|
REAL(dp) :: dk(3)
|
||||||
|
REAL(dp) :: dkmod
|
||||||
|
REAL(dp) :: el_loc
|
||||||
|
REAL(dp) :: eps
|
||||||
|
REAL(dp) :: fac
|
||||||
|
REAL(dp) :: g2kin_bp(npwx)
|
||||||
|
REAL(dp) :: gpar(3)
|
||||||
|
REAL(dp) :: gtr(3)
|
||||||
|
REAL(dp) :: gvec
|
||||||
|
REAL(dp) :: ln(-nr1:nr1,-nr2:nr2,-nr3:nr3)
|
||||||
|
REAL(dp), ALLOCATABLE :: loc_k(:)
|
||||||
|
REAL(dp) , ALLOCATABLE :: pdl_elec(:)
|
||||||
|
REAL(dp) :: pdl_elec_dw
|
||||||
|
REAL(dp) :: pdl_elec_tot
|
||||||
|
REAL(dp) :: pdl_elec_up
|
||||||
|
REAL(dp) :: pdl_ion(nat)
|
||||||
|
REAL(dp) :: pdl_ion_dw
|
||||||
|
REAL(dp) :: pdl_ion_tot
|
||||||
|
REAL(dp) :: pdl_ion_up
|
||||||
|
REAL(dp) :: pdl_tot
|
||||||
|
REAL(dp) , ALLOCATABLE :: phik(:)
|
||||||
|
REAL(dp) :: phidw
|
||||||
|
REAL(dp) :: phiup
|
||||||
|
REAL(dp) :: rmod
|
||||||
|
REAL(dp) :: qrad_dk(nbrx,nbrx,lqx,ntyp)
|
||||||
|
REAL(dp) :: upol(3)
|
||||||
|
REAL(dp) :: weight
|
||||||
|
REAL(dp), ALLOCATABLE :: wstring(:)
|
||||||
|
REAL(dp) :: ylm_dk(lqx*lqx)
|
||||||
|
REAL(dp) :: zeta_mod
|
||||||
|
COMPLEX(dp) :: aux(ngm)
|
||||||
|
COMPLEX(dp) :: aux0(ngm)
|
||||||
|
COMPLEX(dp) :: becp0(nkb,nbnd)
|
||||||
|
COMPLEX(dp) :: becp_bp(nkb,nbnd)
|
||||||
|
COMPLEX(dp) :: cdet(2)
|
||||||
|
COMPLEX(dp) :: cdwork(nbnd)
|
||||||
|
COMPLEX(dp) :: cave
|
||||||
|
COMPLEX(dp) :: cave_dw
|
||||||
|
COMPLEX(dp) :: cave_up
|
||||||
|
COMPLEX(dp) , ALLOCATABLE :: cphik(:)
|
||||||
|
COMPLEX(dp) :: det
|
||||||
|
COMPLEX(dp) :: dtheta
|
||||||
|
COMPLEX(dp) :: mat(nbnd,nbnd)
|
||||||
|
COMPLEX(dp) :: pref
|
||||||
|
COMPLEX(dp) :: psi(npwx,nbnd)
|
||||||
|
COMPLEX(dp) :: q_dk(nhm,nhm,ntyp)
|
||||||
|
COMPLEX(dp) :: struc(nat)
|
||||||
|
COMPLEX(dp) :: theta0
|
||||||
|
COMPLEX(dp) :: zdotc
|
||||||
|
COMPLEX(dp) :: zeta
|
||||||
|
|
||||||
|
|
||||||
|
! ------------------------------------------------------------------------- !
|
||||||
|
! INITIALIZATIONS
|
||||||
|
! ------------------------------------------------------------------------- !
|
||||||
|
|
||||||
|
! --- Write header ---
|
||||||
|
WRITE(6,"(/,/,/,15X,50('='))")
|
||||||
|
WRITE(6,"(28X,'POLARIZATION CALCULATION')")
|
||||||
|
WRITE(6,"(15X,50('-'),/)")
|
||||||
|
|
||||||
|
! --- Check that we are working with an insulator with no empty bands ---
|
||||||
|
IF ((degauss > 0.01) .OR. (nbnd /= nelec/2)) CALL errore('c_phase', &
|
||||||
|
'Polarization only for insulators and no empty bands',1)
|
||||||
|
|
||||||
|
! --- Define a small number ---
|
||||||
|
eps=1.0E-6_dp
|
||||||
|
|
||||||
|
! --- Recalculate FFT correspondence (see ggen.f90) ---
|
||||||
|
DO ng=1,ngm
|
||||||
|
mk1=nint(g(1,ng)*at(1,1)+g(2,ng)*at(2,1)+g(3,ng)*at(3,1))
|
||||||
|
mk2=nint(g(1,ng)*at(1,2)+g(2,ng)*at(2,2)+g(3,ng)*at(3,2))
|
||||||
|
mk3=nint(g(1,ng)*at(1,3)+g(2,ng)*at(2,3)+g(3,ng)*at(3,3))
|
||||||
|
ln(mk1,mk2,mk3) = ng
|
||||||
|
END DO
|
||||||
|
|
||||||
|
! --- Initialize arrays ---
|
||||||
|
jkb_bp=0
|
||||||
|
DO nt=1,ntyp
|
||||||
|
DO na=1,nat
|
||||||
|
IF (ityp(na).eq.nt) THEN
|
||||||
|
DO i=1, nh(nt)
|
||||||
|
jkb_bp=jkb_bp+1
|
||||||
|
nkbtona(jkb_bp) = na
|
||||||
|
nkbtonh(jkb_bp) = i
|
||||||
|
END DO
|
||||||
|
END IF
|
||||||
|
END DO
|
||||||
|
END DO
|
||||||
|
|
||||||
|
! --- Get the number of strings ---
|
||||||
|
nstring=nks/nppstr
|
||||||
|
nkort=nstring/(nspin)
|
||||||
|
|
||||||
|
! --- Allocate memory for arrays ---
|
||||||
|
ALLOCATE(phik(nstring))
|
||||||
|
ALLOCATE(loc_k(nstring))
|
||||||
|
ALLOCATE(cphik(nstring))
|
||||||
|
ALLOCATE(wstring(nstring))
|
||||||
|
ALLOCATE(pdl_elec(nstring))
|
||||||
|
ALLOCATE(mod_elec(nstring))
|
||||||
|
|
||||||
|
|
||||||
|
! ------------------------------------------------------------------------- !
|
||||||
|
! electronic polarization: set values for k-points strings !
|
||||||
|
! ------------------------------------------------------------------------- !
|
||||||
|
|
||||||
|
! --- Find vector along strings ---
|
||||||
|
gpar(1)=xk(1,nppstr)-xk(1,1)
|
||||||
|
gpar(2)=xk(2,nppstr)-xk(2,1)
|
||||||
|
gpar(3)=xk(3,nppstr)-xk(3,1)
|
||||||
|
gvec=dsqrt(gpar(1)**2+gpar(2)**2+gpar(3)**2)*tpiba
|
||||||
|
|
||||||
|
! --- Find vector between consecutive points in strings ---
|
||||||
|
dk(1)=xk(1,2)-xk(1,1)
|
||||||
|
dk(2)=xk(2,2)-xk(2,1)
|
||||||
|
dk(3)=xk(3,2)-xk(3,1)
|
||||||
|
dkmod=SQRT(dk(1)**2+dk(2)**2+dk(3)**2)*tpiba
|
||||||
|
IF (ABS(dkmod-gvec/(nppstr-1)) > eps) &
|
||||||
|
CALL errore('c_phase','Wrong k-strings?',1)
|
||||||
|
|
||||||
|
! --- Check that k-points form strings ---
|
||||||
|
DO i=1,nspin*nkort
|
||||||
|
DO j=2,nppstr
|
||||||
|
kindex=j+(i-1)*nppstr
|
||||||
|
IF (ABS(xk(1,kindex)-xk(1,kindex-1)-dk(1)) > eps) &
|
||||||
|
CALL errore('c_phase','Wrong k-strings?',1)
|
||||||
|
IF (ABS(xk(2,kindex)-xk(2,kindex-1)-dk(2)) > eps) &
|
||||||
|
CALL errore('c_phase','Wrong k-strings?',1)
|
||||||
|
IF (ABS(xk(3,kindex)-xk(3,kindex-1)-dk(3)) > eps) &
|
||||||
|
CALL errore('c_phase','Wrong k-strings?',1)
|
||||||
|
IF (ABS(wk(kindex)-wk(kindex-1)-dk(1)) > eps) &
|
||||||
|
CALL errore('c_phase','Wrong k-strings weights?',1)
|
||||||
|
END DO
|
||||||
|
END DO
|
||||||
|
|
||||||
|
|
||||||
|
! ------------------------------------------------------------------------- !
|
||||||
|
! electronic polarization: weight strings !
|
||||||
|
! ------------------------------------------------------------------------- !
|
||||||
|
|
||||||
|
! --- Calculate string weights, normalizing to 1 (no spin) or 1+1 (spin) ---
|
||||||
|
DO is=1,nspin
|
||||||
|
weight=0.0_dp
|
||||||
|
DO kort=1,nkort
|
||||||
|
istring=kort+(is-1)*nkort
|
||||||
|
wstring(istring)=wk(nppstr*istring)
|
||||||
|
weight=weight+wstring(istring)
|
||||||
|
END DO
|
||||||
|
DO kort=1,nkort
|
||||||
|
istring=kort+(is-1)*nkort
|
||||||
|
wstring(istring)=wstring(istring)/weight
|
||||||
|
END DO
|
||||||
|
END DO
|
||||||
|
|
||||||
|
|
||||||
|
! ------------------------------------------------------------------------- !
|
||||||
|
! electronic polarization: structure factor !
|
||||||
|
! ------------------------------------------------------------------------- !
|
||||||
|
|
||||||
|
! --- Calculate structure factor e^{-i dk*R} ---
|
||||||
|
DO na=1,nat
|
||||||
|
fac=(dk(1)*tau(1,na)+dk(2)*tau(2,na)+dk(3)*tau(3,na))*tpi
|
||||||
|
struc(na)=CMPLX(cos(fac),-sin(fac))
|
||||||
|
END DO
|
||||||
|
|
||||||
|
|
||||||
|
! ------------------------------------------------------------------------- !
|
||||||
|
! electronic polarization: form factor !
|
||||||
|
! ------------------------------------------------------------------------- !
|
||||||
|
|
||||||
|
! --- Calculate Bessel transform of Q_ij(|r|) at dk [Q_ij^L(|r|)] ---
|
||||||
|
CALL calc_btq(dkmod,qrad_dk,0)
|
||||||
|
|
||||||
|
! --- Calculate the q-space real spherical harmonics at dk [Y_LM] ---
|
||||||
|
CALL ylm_q(lqx*lqx,dk,dkmod,ylm_dk)
|
||||||
|
|
||||||
|
! --- Form factor: 4 pi sum_LM c_ij^LM Y_LM(Omega) Q_ij^L(|r|) ---
|
||||||
|
CALL setv(nhm*nhm*ntyp,0.d0,q_dk,1)
|
||||||
|
DO np =1, ntyp
|
||||||
|
DO iv = 1, nh(np)
|
||||||
|
DO jv = iv, nh(np)
|
||||||
|
call qvan3(iv,jv,np,pref,ylm_dk,qrad_dk)
|
||||||
|
q_dk(iv,jv,np) = omega*pref
|
||||||
|
q_dk(jv,iv,np) = omega*pref
|
||||||
|
ENDDO
|
||||||
|
ENDDO
|
||||||
|
ENDDO
|
||||||
|
|
||||||
|
|
||||||
|
! ------------------------------------------------------------------------- !
|
||||||
|
! electronic polarization: strings phases !
|
||||||
|
! ------------------------------------------------------------------------- !
|
||||||
|
|
||||||
|
el_loc = 0.d0
|
||||||
|
kpoint=0
|
||||||
|
|
||||||
|
! --- Start loop over spin ---
|
||||||
|
DO is=1,nspin
|
||||||
|
|
||||||
|
! --- Start loop over orthogonal k-points ---
|
||||||
|
DO kort=1,nkort
|
||||||
|
|
||||||
|
! --- Index for this string ---
|
||||||
|
istring=kort+(is-1)*nkort
|
||||||
|
|
||||||
|
! --- Initialize expectation value of the phase operator ---
|
||||||
|
zeta=(1.d0,0.d0)
|
||||||
|
zeta_mod = 1.d0
|
||||||
|
|
||||||
|
! --- Start loop over parallel k-points ---
|
||||||
|
DO kpar = 1,nppstr
|
||||||
|
|
||||||
|
! --- Set index of k-point ---
|
||||||
|
kpoint = kpoint + 1
|
||||||
|
|
||||||
|
! --- Calculate dot products between wavefunctions and betas ---
|
||||||
|
IF (kpar /= 1) THEN
|
||||||
|
|
||||||
|
! --- Dot wavefunctions and betas for PREVIOUS k-point ---
|
||||||
|
CALL gk_sort(xk(1,kpoint-1),ngm,g,ecutwfc/tpiba2, &
|
||||||
|
npw0,igk0,g2kin_bp)
|
||||||
|
CALL davcio(psi,nwordwfc,iunwfc,kpoint-1,-1)
|
||||||
|
CALL init_us_2 (npw0,igk0,xk(1,kpoint-1),vkb)
|
||||||
|
CALL ccalbec(nkb, npwx, npw, nbnd, becp0, vkb, psi)
|
||||||
|
|
||||||
|
! --- Dot wavefunctions and betas for CURRENT k-point ---
|
||||||
|
IF (kpar /= nppstr) THEN
|
||||||
|
CALL gk_sort(xk(1,kpoint),ngm,g,ecutwfc/tpiba2, &
|
||||||
|
npw1,igk1,g2kin_bp)
|
||||||
|
CALL davcio(evc,nwordwfc,iunwfc,kpoint,-1)
|
||||||
|
CALL init_us_2 (npw1,igk1,xk(1,kpoint),vkb)
|
||||||
|
CALL ccalbec(nkb,npwx,npw,nbnd,becp_bp,vkb,evc)
|
||||||
|
ELSE
|
||||||
|
kstart = kpoint-nppstr+1
|
||||||
|
CALL gk_sort(xk(1,kstart),ngm,g,ecutwfc/tpiba2, &
|
||||||
|
npw1,igk1,g2kin_bp)
|
||||||
|
CALL davcio(evc,nwordwfc,iunwfc,kstart,-1)
|
||||||
|
CALL init_us_2 (npw1,igk1,xk(1,kstart),vkb)
|
||||||
|
CALL ccalbec(nkb,npwx,npw,nbnd,becp_bp,vkb,evc)
|
||||||
|
ENDIF
|
||||||
|
|
||||||
|
! --- Matrix elements calculation ---
|
||||||
|
CALL setv(2*nbnd*nbnd,0.d0,mat,1)
|
||||||
|
DO nb=1,nbnd
|
||||||
|
DO mb=1,nbnd
|
||||||
|
CALL setv(2*ngm,0.d0,aux,1)
|
||||||
|
CALL setv(2*ngm,0.d0,aux0,1)
|
||||||
|
DO ik=1,npw0
|
||||||
|
aux0(igk0(ik))=psi(ik,nb)
|
||||||
|
END DO
|
||||||
|
DO ik=1,npw1
|
||||||
|
IF (kpar /= nppstr) THEN
|
||||||
|
aux(igk1(ik))=evc(ik,mb)
|
||||||
|
ELSE
|
||||||
|
! --- If k'=k+G_o, the relation psi_k+G_o (G-G_o) ---
|
||||||
|
! --- = psi_k(G) is used, gpar=G_o, gtr = G-G_o ---
|
||||||
|
gtr(1)=g(1,igk1(ik))-gpar(1)
|
||||||
|
gtr(2)=g(2,igk1(ik))-gpar(2)
|
||||||
|
gtr(3)=g(3,igk1(ik))-gpar(3)
|
||||||
|
! --- Find crystal coordinates of gtr, n1,n2,n3 ---
|
||||||
|
! --- and the position ng in the ngm array ---
|
||||||
|
IF (gtr(1)**2+gtr(2)**2+gtr(3)**2 <= gcutm) THEN
|
||||||
|
n1=NINT(gtr(1)*at(1,1)+gtr(2)*at(2,1) &
|
||||||
|
+gtr(3)*at(3,1))
|
||||||
|
n2=NINT(gtr(1)*at(1,2)+gtr(2)*at(2,2) &
|
||||||
|
+gtr(3)*at(3,2))
|
||||||
|
n3=NINT(gtr(1)*at(1,3)+gtr(2)*at(2,3) &
|
||||||
|
+gtr(3)*at(3,3))
|
||||||
|
ng=ln(n1,n2,n3)
|
||||||
|
IF ((ABS(g(1,ng)-gtr(1)) > eps) .OR. &
|
||||||
|
(ABS(g(2,ng)-gtr(2)) > eps) .OR. &
|
||||||
|
(ABS(g(3,ng)-gtr(3)) > eps)) THEN
|
||||||
|
WRITE(6,*) ' error: translated G=', &
|
||||||
|
gtr(1),gtr(2),gtr(3), &
|
||||||
|
' with crystal coordinates',n1,n2,n3, &
|
||||||
|
' corresponds to ng=',ng,' but G(ng)=', &
|
||||||
|
g(1,ng),g(2,ng),g(3,ng)
|
||||||
|
WRITE(6,*) ' probably because G_par is NOT', &
|
||||||
|
' a reciprocal lattice vector '
|
||||||
|
WRITE(6,*) ' Possible choices as smallest ', &
|
||||||
|
' G_par:'
|
||||||
|
DO i=1,50
|
||||||
|
WRITE(6,*) ' i=',i,' G=', &
|
||||||
|
g(1,i),g(2,i),g(3,i)
|
||||||
|
ENDDO
|
||||||
|
STOP
|
||||||
|
ENDIF
|
||||||
|
ELSE
|
||||||
|
WRITE(6,*) ' |gtr| > gcutm for gtr=', &
|
||||||
|
gtr(1),gtr(2),gtr(3)
|
||||||
|
STOP
|
||||||
|
END IF
|
||||||
|
aux(ng)=evc(ik,mb)
|
||||||
|
ENDIF
|
||||||
|
END DO
|
||||||
|
mat(nb,mb) = zdotc(ngm,aux0,1,aux,1)
|
||||||
|
! --- Calculate the augmented part: ij=KB projectors, ---
|
||||||
|
! --- R=atom index: SUM_{ijR} q(ijR) <u_nk|beta_iR> ---
|
||||||
|
! --- <beta_jR|u_mk'> e^i(k-k')*R = ---
|
||||||
|
! --- also <u_nk|beta_iR>=<psi_nk|beta_iR> = becp^* ---
|
||||||
|
pref = (0.d0,0.d0)
|
||||||
|
DO jkb=1,nkb
|
||||||
|
nhjkb = nkbtonh(jkb)
|
||||||
|
na = nkbtona(jkb)
|
||||||
|
np = ityp(na)
|
||||||
|
nhjkbm = nh(np)
|
||||||
|
jkb1 = jkb - nhjkb
|
||||||
|
DO j = 1,nhjkbm
|
||||||
|
pref = pref+conjg(becp0(jkb,nb))*becp_bp(jkb1+j,mb) &
|
||||||
|
*q_dk(nhjkb,j,np)*struc(na)
|
||||||
|
ENDDO
|
||||||
|
ENDDO
|
||||||
|
mat(nb,mb) = mat(nb,mb) + pref
|
||||||
|
ENDDO
|
||||||
|
ENDDO
|
||||||
|
|
||||||
|
! --- Calculate matrix determinant ---
|
||||||
|
CALL zgefa(mat,nbnd,nbnd,ivpt,info)
|
||||||
|
CALL errore('c_phase','error in zgefa',abs(info))
|
||||||
|
job=10
|
||||||
|
CALL zgedi(mat,nbnd,nbnd,ivpt,cdet,cdwork,job)
|
||||||
|
det=cdet(1)*10.d0**cdet(2)
|
||||||
|
|
||||||
|
! --- Multiply by the already calculated determinants ---
|
||||||
|
zeta=zeta*det
|
||||||
|
|
||||||
|
! --- End of dot products between wavefunctions and betas ---
|
||||||
|
ENDIF
|
||||||
|
|
||||||
|
! --- End loop over parallel k-points ---
|
||||||
|
END DO
|
||||||
|
|
||||||
|
! --- Calculate the phase for this string ---
|
||||||
|
phik(istring)=IMAG(LOG(zeta))
|
||||||
|
cphik(istring)=COS(phik(istring))*(1.0_dp,0.0_dp) &
|
||||||
|
+SIN(phik(istring))*(0.0_dp,1.0_dp)
|
||||||
|
|
||||||
|
! --- Calculate the localization for current kort ---
|
||||||
|
zeta_mod=dreal(conjg(zeta)*zeta)
|
||||||
|
loc_k(istring)= - (nppstr-1) / gvec**2 / nbnd *log(zeta_mod)
|
||||||
|
|
||||||
|
! --- End loop over orthogonal k-points ---
|
||||||
|
END DO
|
||||||
|
|
||||||
|
! --- End loop over spin ---
|
||||||
|
END DO
|
||||||
|
|
||||||
|
|
||||||
|
! ------------------------------------------------------------------------- !
|
||||||
|
! electronic polarization: phase average !
|
||||||
|
! ------------------------------------------------------------------------- !
|
||||||
|
|
||||||
|
! --- Initializations ---
|
||||||
|
cave_up=(0.0_dp,0.0_dp)
|
||||||
|
cave_dw=(0.0_dp,0.0_dp)
|
||||||
|
|
||||||
|
! --- Start loop over spins ---
|
||||||
|
DO is=1,nspin
|
||||||
|
|
||||||
|
! --- Initialize average of phases as complex numbers ---
|
||||||
|
cave=(0.0_dp,0.0_dp)
|
||||||
|
|
||||||
|
! --- Start loop over strings with same spin ---
|
||||||
|
DO kort=1,nkort
|
||||||
|
|
||||||
|
! --- Calculate string index ---
|
||||||
|
istring=kort+(is-1)*nkort
|
||||||
|
|
||||||
|
! --- Average phases as complex numbers ---
|
||||||
|
cave=cave+wstring(istring)*cphik(istring)
|
||||||
|
|
||||||
|
! --- End loop over strings with same spin ---
|
||||||
|
END DO
|
||||||
|
|
||||||
|
! --- Get the angle corresponding to the complex numbers average ---
|
||||||
|
theta0=atan2(IMAG(cave),REAL(cave))
|
||||||
|
|
||||||
|
! --- Assign this angle to the corresponding spin phase average ---
|
||||||
|
IF (nspin == 1) THEN
|
||||||
|
phiup=theta0
|
||||||
|
phidw=theta0
|
||||||
|
ELSE IF (nspin == 2) THEN
|
||||||
|
IF (is == 1) THEN
|
||||||
|
phiup=theta0
|
||||||
|
ELSE IF (is == 2) THEN
|
||||||
|
phidw=theta0
|
||||||
|
END IF
|
||||||
|
END IF
|
||||||
|
|
||||||
|
! --- Put the phases in an around theta0 ---
|
||||||
|
cphik(istring)=cphik(istring)/cave
|
||||||
|
dtheta=atan2(IMAG(cphik(istring)),REAL(cphik(istring)))
|
||||||
|
phik(istring)=theta0+dtheta
|
||||||
|
|
||||||
|
! --- End loop over spins
|
||||||
|
END DO
|
||||||
|
|
||||||
|
|
||||||
|
! ------------------------------------------------------------------------- !
|
||||||
|
! electronic polarization: remap phases !
|
||||||
|
! ------------------------------------------------------------------------- !
|
||||||
|
|
||||||
|
! --- Remap string phases to interval [-0.5,0.5) ---
|
||||||
|
pdl_elec=phik/(2.0_dp*pi)
|
||||||
|
mod_elec=1
|
||||||
|
|
||||||
|
! --- Remap spin average phases to interval [-0.5,0.5) ---
|
||||||
|
pdl_elec_up=phiup/(2.0_dp*pi)
|
||||||
|
mod_elec_up=1
|
||||||
|
pdl_elec_dw=phidw/(2.0_dp*pi)
|
||||||
|
mod_elec_dw=1
|
||||||
|
|
||||||
|
! --- Depending on nspin, remap total phase to [-1,1) or [-0.5,0.5) ---
|
||||||
|
pdl_elec_tot=pdl_elec_up+pdl_elec_dw
|
||||||
|
IF (nspin == 1) THEN
|
||||||
|
pdl_elec_tot=pdl_elec_tot-2.0_dp*NINT(pdl_elec_tot/2.0_dp)
|
||||||
|
mod_elec_tot=2
|
||||||
|
ELSE IF (nspin == 2) THEN
|
||||||
|
pdl_elec_tot=pdl_elec_tot-1.0_dp*NINT(pdl_elec_tot/1.0_dp)
|
||||||
|
mod_elec_tot=1
|
||||||
|
END IF
|
||||||
|
|
||||||
|
|
||||||
|
! ------------------------------------------------------------------------- !
|
||||||
|
! ionic polarization !
|
||||||
|
! ------------------------------------------------------------------------- !
|
||||||
|
|
||||||
|
! --- Look for ions with odd number of charges ---
|
||||||
|
mod_ion=2
|
||||||
|
lodd=.FALSE.
|
||||||
|
DO na=1,nat
|
||||||
|
IF (MOD(NINT(zv(ityp(na))),2) == 1) THEN
|
||||||
|
mod_ion(na)=1
|
||||||
|
lodd=.TRUE.
|
||||||
|
END IF
|
||||||
|
END DO
|
||||||
|
|
||||||
|
! --- Calculate ionic polarization phase for every ion ---
|
||||||
|
pdl_ion=0.0_dp
|
||||||
|
DO na=1,nat
|
||||||
|
DO i=1,3
|
||||||
|
pdl_ion(na)=pdl_ion(na)+zv(ityp(na))*tau(i,na)*gpar(i)
|
||||||
|
ENDDO
|
||||||
|
IF (mod_ion(na) == 1) THEN
|
||||||
|
pdl_ion(na)=pdl_ion(na)-1.0_dp*nint(pdl_ion(na)/1.0_dp)
|
||||||
|
ELSE IF (mod_ion(na) == 2) THEN
|
||||||
|
pdl_ion(na)=pdl_ion(na)-2.0_dp*nint(pdl_ion(na)/2.0_dp)
|
||||||
|
END IF
|
||||||
|
ENDDO
|
||||||
|
|
||||||
|
! --- Add up the phases modulo 2 iff the ionic charges are even numbers ---
|
||||||
|
pdl_ion_tot=SUM(pdl_ion(1:nat))
|
||||||
|
IF (lodd) THEN
|
||||||
|
pdl_ion_tot=pdl_ion_tot-1.d0*nint(pdl_ion_tot/1.d0)
|
||||||
|
mod_ion_tot=1
|
||||||
|
ELSE
|
||||||
|
pdl_ion_tot=pdl_ion_tot-2.d0*nint(pdl_ion_tot/2.d0)
|
||||||
|
mod_ion_tot=2
|
||||||
|
END IF
|
||||||
|
|
||||||
|
|
||||||
|
! ------------------------------------------------------------------------- !
|
||||||
|
! total polarization !
|
||||||
|
! ------------------------------------------------------------------------- !
|
||||||
|
|
||||||
|
! --- Add electronic and ionic contributions to total phase ---
|
||||||
|
pdl_tot=pdl_elec_tot+pdl_ion_tot
|
||||||
|
IF ((.NOT.lodd).AND.(nspin == 1)) THEN
|
||||||
|
mod_tot=2
|
||||||
|
ELSE
|
||||||
|
mod_tot=1
|
||||||
|
END IF
|
||||||
|
|
||||||
|
|
||||||
|
! ------------------------------------------------------------------------- !
|
||||||
|
! write output information !
|
||||||
|
! ------------------------------------------------------------------------- !
|
||||||
|
|
||||||
|
! --- Information about the k-points string used ---
|
||||||
|
WRITE(6,"(/,21X,'K-POINTS STRINGS USED IN CALCULATIONS')")
|
||||||
|
WRITE(6,"(21X,37('~'),/)")
|
||||||
|
WRITE(6,"(7X,'G-vector along string (2 pi/a):',3F9.5)") &
|
||||||
|
gpar(1),gpar(2),gpar(3)
|
||||||
|
WRITE(6,"(7X,'Modulus of the vector (1/bohr):',F9.5)") &
|
||||||
|
gvec
|
||||||
|
WRITE(6,"(7X,'Number of k-points per string:',I4)") nppstr
|
||||||
|
WRITE(6,"(7X,'Number of different strings :',I4)") nkort
|
||||||
|
|
||||||
|
! --- Information about ionic polarization phases ---
|
||||||
|
WRITE(6,"(2/,31X,'IONIC POLARIZATION')")
|
||||||
|
WRITE(6,"(31X,18('~'),/)")
|
||||||
|
WRITE(6,"(8X,'Note: (mod 1) means that the phases (angles ranging from' &
|
||||||
|
/,8X,'-pi to pi) have been mapped to the interval [-1/2,+1/2) by',&
|
||||||
|
/,8X,'dividing by 2*pi; (mod 2) refers to the interval [-1,+1)',&
|
||||||
|
/)")
|
||||||
|
WRITE(6,"(2X,76('='))")
|
||||||
|
WRITE(6,"(4X,'Ion',4X,'Species',4X,'Charge',14X, &
|
||||||
|
'Position',16X,'Phase')")
|
||||||
|
WRITE(6,"(2X,76('-'))")
|
||||||
|
DO na=1,nat
|
||||||
|
WRITE(6,"(3X,I3,8X,A2,F12.3,5X,3F8.4,F12.5,' (mod ',I1,')')") &
|
||||||
|
na,atm(ityp(na)),zv(ityp(na)), &
|
||||||
|
tau(1,na),tau(2,na),tau(3,na),pdl_ion(na),mod_ion(na)
|
||||||
|
END DO
|
||||||
|
WRITE(6,"(2X,76('-'))")
|
||||||
|
WRITE(6,"(47X,'IONIC PHASE: ',F9.5,' (mod ',I1,')')") pdl_ion_tot,mod_ion_tot
|
||||||
|
WRITE(6,"(2X,76('='))")
|
||||||
|
|
||||||
|
! --- Information about electronic polarization phases ---
|
||||||
|
WRITE(6,"(2/,28X,'ELECTRONIC POLARIZATION')")
|
||||||
|
WRITE(6,"(28X,23('~'),/)")
|
||||||
|
WRITE(6,"(8X,'Note: (mod 1) means that the phases (angles ranging from' &
|
||||||
|
/,8X,'-pi to pi) have been mapped to the interval [-1/2,+1/2) by',&
|
||||||
|
/,8X,'dividing by 2*pi; (mod 2) refers to the interval [-1,+1)',&
|
||||||
|
/)")
|
||||||
|
WRITE(6,"(2X,76('='))")
|
||||||
|
WRITE(6,"(3X,'Spin',4X,'String',5X,'Weight',6X, &
|
||||||
|
'First k-point in string',9X,'Phase')")
|
||||||
|
WRITE(6,"(2X,76('-'))")
|
||||||
|
DO istring=1,nstring/nspin
|
||||||
|
ind1=1+(istring-1)*nppstr
|
||||||
|
WRITE(6,"(3X,' up ',3X,I5,F14.6,4X,3(F8.4),F12.5' (mod ',I1,')')") &
|
||||||
|
istring,wstring(istring), &
|
||||||
|
xk(1,ind1),xk(2,ind1),xk(3,ind1),pdl_elec(istring),mod_elec(istring)
|
||||||
|
END DO
|
||||||
|
WRITE(6,"(2X,76('-'))")
|
||||||
|
! --- Treat unpolarized/polarized spin cases ---
|
||||||
|
IF (nspin == 1) THEN
|
||||||
|
! --- In unpolarized spin, just copy again the same data ---
|
||||||
|
DO istring=1,nstring
|
||||||
|
ind1=1+(istring-1)*nppstr
|
||||||
|
WRITE(6,"(3X,'down',3X,I5,F14.6,4X,3(F8.4),F12.5' (mod ',I1,')')") &
|
||||||
|
istring,wstring(istring), xk(1,ind1),xk(2,ind1),xk(3,ind1), &
|
||||||
|
pdl_elec(istring),mod_elec(istring)
|
||||||
|
END DO
|
||||||
|
ELSE IF (nspin == 2) THEN
|
||||||
|
! --- If there is spin polarization, write information for new strings ---
|
||||||
|
DO istring=nstring/2+1,nstring
|
||||||
|
ind1=1+(istring-1)*nppstr
|
||||||
|
WRITE(6,"(3X,'down',3X,I4,F15.6,4X,3(F8.4),F12.5' (mod ',I1,')')") &
|
||||||
|
istring,wstring(istring), xk(1,ind1),xk(2,ind1),xk(3,ind1), &
|
||||||
|
pdl_elec(istring),mod_elec(istring)
|
||||||
|
END DO
|
||||||
|
END IF
|
||||||
|
WRITE(6,"(2X,76('-'))")
|
||||||
|
WRITE(6,"(40X,'Average phase (up): ',F9.5,' (mod ',I1,')')") &
|
||||||
|
pdl_elec_up,mod_elec_up
|
||||||
|
WRITE(6,"(38X,'Average phase (down): ',F9.5,' (mod ',I1,')')")&
|
||||||
|
pdl_elec_dw,mod_elec_dw
|
||||||
|
WRITE(6,"(42X,'ELECTRONIC PHASE: ',F9.5,' (mod ',I1,')')") &
|
||||||
|
pdl_elec_tot,mod_elec_tot
|
||||||
|
WRITE(6,"(2X,76('='))")
|
||||||
|
|
||||||
|
! --- Information about total phase ---
|
||||||
|
WRITE(6,"(2/,31X,'SUMMARY OF PHASES')")
|
||||||
|
WRITE(6,"(31X,17('~'),/)")
|
||||||
|
WRITE(6,"(26X,'Ionic Phase:',F9.5,' (mod ',I1,')')") &
|
||||||
|
pdl_ion_tot,mod_ion_tot
|
||||||
|
WRITE(6,"(21X,'Electronic Phase:',F9.5,' (mod ',I1,')')") &
|
||||||
|
pdl_elec_tot,mod_elec_tot
|
||||||
|
WRITE(6,"(26X,'TOTAL PHASE:',F9.5,' (mod ',I1,')')") &
|
||||||
|
pdl_tot,mod_tot
|
||||||
|
|
||||||
|
! --- Information about the value of polarization ---
|
||||||
|
WRITE(6,"(2/,29X,'VALUES OF POLARIZATION')")
|
||||||
|
WRITE(6,"(29X,22('~'),/)")
|
||||||
|
WRITE(6,"( &
|
||||||
|
8X,'The calculation of phases done along the direction of vector ',I1, &
|
||||||
|
/,8X,'of the reciprocal lattice gives the following contribution to', &
|
||||||
|
/,8X,'the polarization vector (in different units, and being Omega', &
|
||||||
|
/,8X,'the volume of the unit cell):')") &
|
||||||
|
gdir
|
||||||
|
! --- Calculate direction of polarization and modulus of lattice vector ---
|
||||||
|
rmod=SQRT(at(1,gdir)*at(1,gdir)+at(2,gdir)*at(2,gdir) &
|
||||||
|
+at(3,gdir)*at(3,gdir))
|
||||||
|
upol(:)=at(:,gdir)/rmod
|
||||||
|
rmod=alat*rmod
|
||||||
|
! --- Give polarization in units of (e/Omega).bohr ---
|
||||||
|
fac=rmod
|
||||||
|
WRITE(6,"(/,11X'P = ',F11.7,' (mod ',F11.7,') (e/Omega).bohr')") &
|
||||||
|
fac*pdl_tot,fac*REAL(mod_tot)
|
||||||
|
! --- Give polarization in units of e.bohr ---
|
||||||
|
fac=rmod/omega
|
||||||
|
WRITE(6,"(/,11X'P = ',F11.7,' (mod ',F11.7,') e/bohr^2')") &
|
||||||
|
fac*pdl_tot,fac*REAL(mod_tot)
|
||||||
|
! --- Give polarization in SI units (C/m^2) ---
|
||||||
|
fac=(rmod/omega)*(1.60097E-19_dp/5.29177E-11_dp**2)
|
||||||
|
WRITE(6,"(/,11X'P = ',F11.7,' (mod ',F11.7,') C/m^2')") &
|
||||||
|
fac*pdl_tot,fac*REAL(mod_tot)
|
||||||
|
! --- Write polarization direction ---
|
||||||
|
WRITE(6,"(/,8X,'The polarization direction is: ( ', &
|
||||||
|
F7.5,' , ',F7.5,' , ',F7.5,' )'))") upol(1),upol(2),upol(3)
|
||||||
|
|
||||||
|
! --- End of information relative to polarization calculation ---
|
||||||
|
WRITE(6,"(/,/,15X,50('=')/,/,)")
|
||||||
|
|
||||||
|
|
||||||
|
! ------------------------------------------------------------------------- !
|
||||||
|
! finalization !
|
||||||
|
! ------------------------------------------------------------------------- !
|
||||||
|
|
||||||
|
! --- Free memory ---
|
||||||
|
DEALLOCATE(pdl_elec)
|
||||||
|
DEALLOCATE(mod_elec)
|
||||||
|
DEALLOCATE(wstring)
|
||||||
|
DEALLOCATE(loc_k)
|
||||||
|
DEALLOCATE(phik)
|
||||||
|
DEALLOCATE(cphik)
|
||||||
|
|
||||||
|
|
||||||
|
!------------------------------------------------------------------------------!
|
||||||
|
|
||||||
|
END SUBROUTINE c_phase
|
||||||
|
|
||||||
|
!==============================================================================!
|
|
@ -0,0 +1,86 @@
|
||||||
|
!----------------------------------------------------------------------
|
||||||
|
subroutine calc_btq(ql,qr_k,idbes)
|
||||||
|
!----------------------------------------------------------------------
|
||||||
|
!
|
||||||
|
! Calculates the Bessel-transform (or its derivative if idbes=1)
|
||||||
|
! of the augmented qrad charges at a given ql point.
|
||||||
|
! Rydberg atomic units are used.
|
||||||
|
!
|
||||||
|
use pwcom
|
||||||
|
!
|
||||||
|
integer :: ik, msh_bp, i, np, m, k, l
|
||||||
|
integer :: n,idbes,ilmin,ilmax,iv,jv
|
||||||
|
real(DP) :: jl(ndm), ql, sum, jlp1(ndm), aux(ndm), &
|
||||||
|
qr_k(nbrx,nbrx,lqx,ntyp)
|
||||||
|
|
||||||
|
! declaration readvan quantities
|
||||||
|
! integer NBETA,KKBETA,iver,nqf,ifqopt,nqlc,lll
|
||||||
|
! real*8 DION,BETAR,QQQ,QFUNC,qfcoef,rinner
|
||||||
|
! COMMON/NCPRM/DION(NBRX,NBRX,NPSX),
|
||||||
|
! C BETAR(0:ndm,NBRX,NPSX),QQQ(NBRX,NBRX,NPSX),
|
||||||
|
! C QFUNC(0:ndm,NBRX,NBRX,NPSX),
|
||||||
|
! C NBETA(NPSX),KKBETA(NPSX),NVALES(NPSX),lll(nbrx,npsx),
|
||||||
|
! C iver(3,npsx),nqf(npsx),ifqopt(npsx),nqlc(npsx),
|
||||||
|
! C qfcoef(nqfm,lqx,NBRX,NBRX,npsx),rinner(lqx,npsx)
|
||||||
|
! common/ncprm/dion(nbrx,nbrx,npsx),
|
||||||
|
! + betar(0:ndm,nbrx,npsx), qqq(nbrx,nbrx,npsx),
|
||||||
|
! + qfunc(0:ndm,nbrx,nbrx,npsx),
|
||||||
|
! + qfcoef(nqfm,lqx,nbrx,nbrx,npsx), rinner(lqx,npsx),
|
||||||
|
! + nbeta(npsx), kkbeta(npsx),
|
||||||
|
! + nqf(npsx), nqlc(npsx), ifqopt(npsx), lll(nbrx,npsx),
|
||||||
|
! + iver(3,npsx)
|
||||||
|
|
||||||
|
|
||||||
|
!
|
||||||
|
do np=1,ntyp
|
||||||
|
msh_bp=kkbeta(np)
|
||||||
|
if (tvanp(np)) then
|
||||||
|
do iv =1, nbeta(np)
|
||||||
|
do jv =iv, nbeta(np)
|
||||||
|
ilmin = iabs(lll(iv,np)-lll(jv,np))
|
||||||
|
ilmax = iabs(lll(iv,np)+lll(jv,np))
|
||||||
|
! only need to calculate for for lmin,lmin+2 ...lmax-2,lmax
|
||||||
|
do l = ilmin,ilmax,2
|
||||||
|
do i = msh_bp,2,-1
|
||||||
|
if (r(i,np) .lt. rinner(l+1,np)) goto 100
|
||||||
|
aux(i) = qfunc(i,iv,jv,np)
|
||||||
|
enddo
|
||||||
|
100 call setqf(qfcoef(1,l+1,iv,jv,np),aux(1),r(1,np) &
|
||||||
|
,nqf(np),l,i)
|
||||||
|
|
||||||
|
if (idbes .eq. 1) then
|
||||||
|
call dbess(ql,l+1,msh_bp,r(1,np), &
|
||||||
|
jl)
|
||||||
|
else
|
||||||
|
call bess(ql,l+1,msh_bp,r(1,np), &
|
||||||
|
jl)
|
||||||
|
endif
|
||||||
|
|
||||||
|
! jl is the Bessel function (or its derivative) calculated at ql
|
||||||
|
! now integrate qfunc*jl*r^2 = Bessel transform of qfunc
|
||||||
|
|
||||||
|
do i=1, msh_bp
|
||||||
|
jlp1(i) = jl(i)*aux(i)
|
||||||
|
enddo
|
||||||
|
! if (tlog(np)) then
|
||||||
|
if(tvanp(np)) then
|
||||||
|
call radlg1(msh_bp,jlp1,rab(1,np),sum)
|
||||||
|
else
|
||||||
|
call radlg(msh_bp,jlp1,r(1,np),dx(np),sum)
|
||||||
|
endif
|
||||||
|
! else
|
||||||
|
! call radin(msh_bp,dx(np),jlp1,sum)
|
||||||
|
! endif
|
||||||
|
qr_k(iv,jv,l+1,np) = sum*fpi/omega
|
||||||
|
qr_k(jv,iv,l+1,np) = qr_k(iv,jv,l+1,np)
|
||||||
|
|
||||||
|
!c write(6,*) 'qr_k=',qr_k(iv,jv,l+1,np)
|
||||||
|
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
enddo
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
!
|
||||||
|
return
|
||||||
|
end
|
|
@ -0,0 +1,92 @@
|
||||||
|
C
|
||||||
|
C-------------------------------------------------------------------------
|
||||||
|
SUBROUTINE DBESS(XG,L,MMAX,R,DJL)
|
||||||
|
C-------------------------------------------------------------------------
|
||||||
|
C CALCULATES DERIVATIVES OF SPHERICAL BESSEL FUNCTIONS j_l(Gr)
|
||||||
|
C WITH RESPECT TO h_alpha,beta (WITHOUT THE FACTOR GAGK(KK,IG)*HTM1)
|
||||||
|
C I.E. -x * D(jl(x))/dx
|
||||||
|
IMPLICIT REAL*8 (A-H,O-Z)
|
||||||
|
PARAMETER(EPS=1.E-8)
|
||||||
|
REAL*8 DJL(MMAX),R(MMAX)
|
||||||
|
|
||||||
|
IF(L.EQ.1) THEN ! S PART
|
||||||
|
IF(XG.LT.EPS) THEN
|
||||||
|
DO IR=1,MMAX
|
||||||
|
DJL(IR) = 0.D0
|
||||||
|
END DO
|
||||||
|
ELSE
|
||||||
|
DJL(1) = 0.D0
|
||||||
|
DO IR=2,MMAX
|
||||||
|
XRG=R(IR)*XG
|
||||||
|
DJL(IR) = SIN(XRG)/XRG-COS(XRG)
|
||||||
|
END DO
|
||||||
|
ENDIF
|
||||||
|
ENDIF
|
||||||
|
|
||||||
|
IF(L.EQ.2) THEN ! P PART
|
||||||
|
IF(XG.LT.EPS) THEN
|
||||||
|
DO IR=1,MMAX
|
||||||
|
DJL(IR) = 0.D0
|
||||||
|
END DO
|
||||||
|
ELSE
|
||||||
|
DJL(1) = 0.D0
|
||||||
|
DO IR=2,MMAX
|
||||||
|
XRG=R(IR)*XG
|
||||||
|
DJL(IR) = 2.D0*(SIN(XRG)/XRG-COS(XRG))/XRG - SIN(XRG)
|
||||||
|
END DO
|
||||||
|
ENDIF
|
||||||
|
ENDIF
|
||||||
|
|
||||||
|
IF(L.EQ.3) THEN ! D PART
|
||||||
|
IF(XG.LT.EPS) THEN
|
||||||
|
DO IR=1,MMAX
|
||||||
|
DJL(IR) = 0.D0
|
||||||
|
END DO
|
||||||
|
ELSE
|
||||||
|
DJL(1) = 0.D0
|
||||||
|
DO IR=2,MMAX
|
||||||
|
XRG=R(IR)*XG
|
||||||
|
DJL(IR) = ( SIN(XRG)*(9.D0/(XRG*XRG)-4.D0) -
|
||||||
|
- 9.D0*COS(XRG)/XRG ) /XRG + COS(XRG)
|
||||||
|
END DO
|
||||||
|
ENDIF
|
||||||
|
ENDIF
|
||||||
|
|
||||||
|
IF(L.EQ.4) THEN ! F PART
|
||||||
|
IF(XG.LT.EPS) THEN
|
||||||
|
DO IR=1,MMAX
|
||||||
|
DJL(IR) = 0.D0
|
||||||
|
END DO
|
||||||
|
ELSE
|
||||||
|
DJL(1) = 0.D0
|
||||||
|
DO IR=2,MMAX
|
||||||
|
XRG=R(IR)*XG
|
||||||
|
XRG2=XRG*XRG
|
||||||
|
DJL(IR)=SIN(XRG)*(60.D0/(XRG2*XRG2)-27.D0/XRG2+1.d0)
|
||||||
|
$ -COS(XRG)*(60.D0/XRG2-7.D0)/XRG
|
||||||
|
END DO
|
||||||
|
ENDIF
|
||||||
|
ENDIF
|
||||||
|
|
||||||
|
IF(L.EQ.5) THEN ! G PART
|
||||||
|
IF(XG.LT.EPS) THEN
|
||||||
|
DO IR=1,MMAX
|
||||||
|
DJL(IR) = 0.D0
|
||||||
|
END DO
|
||||||
|
ELSE
|
||||||
|
DJL(1) = 0.D0
|
||||||
|
DO IR=2,MMAX
|
||||||
|
XRG=R(IR)*XG
|
||||||
|
XRG2=XRG*XRG
|
||||||
|
DJL(IR)=SIN(XRG)*(525.D0/(XRG2*XRG2)-240.D0/XRG2+11.D0)/XRG
|
||||||
|
$ - COS(XRG)*(525.D0/(XRG2*XRG2)-65.D0/XRG2+1.D0)
|
||||||
|
END DO
|
||||||
|
ENDIF
|
||||||
|
ENDIF
|
||||||
|
|
||||||
|
IF(L.LE.0 .OR. L.GE.6) THEN
|
||||||
|
CALL ERRORE('DBESS',' L NOT PROGRAMMED, L= ',L)
|
||||||
|
END IF
|
||||||
|
|
||||||
|
RETURN
|
||||||
|
END
|
|
@ -0,0 +1,68 @@
|
||||||
|
!
|
||||||
|
!--------------------------------------------------------------------------
|
||||||
|
subroutine qvan3(iv,jv,is,qg,ylm_k,qr)
|
||||||
|
!--------------------------------------------------------------------------
|
||||||
|
!
|
||||||
|
! calculate qg = SUM_LM (-I)^L AP(LM,iv,jv) YR_LM QRAD(iv,jv,L,is)
|
||||||
|
|
||||||
|
use pwcom
|
||||||
|
integer :: iv,jv,is
|
||||||
|
complex(DP) :: qg,sig
|
||||||
|
real(DP) :: ylm_k(lqx*lqx)
|
||||||
|
real(DP) :: qr(nbrx,nbrx,lqx,ntyp)
|
||||||
|
|
||||||
|
integer ivs,jvs,ivl,jvl,ig,lp,l,i
|
||||||
|
! IV = 1..8 ! s_1 p_x1 p_y1 p_z1 s_2 p_x2 p_z2 p_y2
|
||||||
|
! IVS = 1..4 ! s_1 s_2 p_1 p_2 d_1 d_2
|
||||||
|
! IVL = 1..4 ! s p_x p_y p_z
|
||||||
|
!
|
||||||
|
! NOTE : IV = 1..8 (sppp sppp) IVS = 1..4 (sspp) OR 1..2 (sp)
|
||||||
|
! IVL = 1..4 (sppp)
|
||||||
|
!
|
||||||
|
ivs = indv(iv,is)
|
||||||
|
jvs = indv(jv,is)
|
||||||
|
ivl = nhtol(iv,is)*nhtol(iv,is)+nhtom(iv,is)
|
||||||
|
jvl = nhtol(jv,is)*nhtol(jv,is)+nhtom(jv,is)
|
||||||
|
|
||||||
|
IF(IVL.GT.NLX) CALL ERRORE(' QVAN ',' IVL.GT.NLX ',IVL)
|
||||||
|
IF(JVL.GT.NLX) CALL ERRORE(' QVAN ',' JVL.GT.NLX ',JVL)
|
||||||
|
IF(IVS.GT.NBRX) CALL ERRORE(' QVAN ',' IVS.GT.NBRX ',IVS)
|
||||||
|
IF(JVS.GT.NBRX) CALL ERRORE(' QVAN ',' JVS.GT.NBRX ',JVS)
|
||||||
|
|
||||||
|
qg = (0.0d0,0.0d0)
|
||||||
|
|
||||||
|
!odl Write(*,*) 'QVAN3 -- ivs jvs = ',ivs,jvs
|
||||||
|
!odl Write(*,*) 'QVAN3 -- ivl jvl = ',ivl,jvl
|
||||||
|
do i=1,lpx(ivl,jvl)
|
||||||
|
!odl Write(*,*) 'QVAN3 -- i = ',i
|
||||||
|
lp = lpl(ivl,jvl,i)
|
||||||
|
!odl Write(*,*) 'QVAN3 -- lp = ',lp
|
||||||
|
|
||||||
|
! EXTRACTION OF ANGULAR MOMENT L FROM LP:
|
||||||
|
|
||||||
|
if (lp.eq.1) then
|
||||||
|
l = 1
|
||||||
|
else if ((lp.ge.2) .and. (lp.le.4)) then
|
||||||
|
l = 2
|
||||||
|
else if ((lp.ge.5) .and. (lp.le.9)) then
|
||||||
|
l = 3
|
||||||
|
else if ((lp.ge.10).and.(lp.le.16)) then
|
||||||
|
l = 4
|
||||||
|
else if ((lp.ge.17).and.(lp.le.25)) then
|
||||||
|
l = 5
|
||||||
|
else if (lp.ge.26) then
|
||||||
|
call errore(' qvan3 ',' lp.ge.26 ',lp)
|
||||||
|
end if
|
||||||
|
|
||||||
|
sig = (0.d0,-1.d0)**(l-1)
|
||||||
|
sig = sig * ap(lp,ivl,jvl)
|
||||||
|
|
||||||
|
!odl Write(*,*) 'QVAN3 -- sig = ',sig
|
||||||
|
|
||||||
|
! write(*,*) 'qvan3',ng1,LP,L,ivs,jvs
|
||||||
|
|
||||||
|
qg = qg + sig * ylm_k(lp) * qr(ivs,jvs,l,is)
|
||||||
|
|
||||||
|
end do
|
||||||
|
return
|
||||||
|
end
|
|
@ -0,0 +1,107 @@
|
||||||
|
C
|
||||||
|
C ---------------------------------------------------------------
|
||||||
|
SUBROUTINE RADIN(MESH,C,FUNC,ASUM)
|
||||||
|
C ---------------------------------------------------------------
|
||||||
|
C SIMPSONS RULE INTEGRATION FOR HERMAN SKILLMAN MESH
|
||||||
|
C MESH - # OF MESH POINTS
|
||||||
|
C C - 0.8853418/Z**(1/3.)
|
||||||
|
C
|
||||||
|
IMPLICIT REAL*8 (A-H,O-Z)
|
||||||
|
DIMENSION FUNC(mesh)
|
||||||
|
A1=0.0
|
||||||
|
A2E=0.0
|
||||||
|
ASUM=0.0
|
||||||
|
H=0.0025*C
|
||||||
|
NBLOCK=MESH/40
|
||||||
|
I=1
|
||||||
|
c FUNC(1)=0.0
|
||||||
|
DO 39 J=1,NBLOCK
|
||||||
|
DO 38 K=1,20
|
||||||
|
I=I+2
|
||||||
|
I1=I-1
|
||||||
|
A2ES=A2E
|
||||||
|
A2O=FUNC(I1)/12.0
|
||||||
|
A2E=FUNC(I)/12.0
|
||||||
|
A1=A1+5.0*A2ES+8.0*A2O-A2E
|
||||||
|
c FUNC(I1)=ASUM+A1*H
|
||||||
|
A1=A1-A2ES+8.0*A2O+5.0*A2E
|
||||||
|
c FUNC(I)=ASUM+A1*H
|
||||||
|
fi = ASUM+A1*H
|
||||||
|
38 CONTINUE
|
||||||
|
c ASUM=FUNC(I)
|
||||||
|
asum = fi
|
||||||
|
A1=0.0
|
||||||
|
39 H=H+H
|
||||||
|
C
|
||||||
|
RETURN
|
||||||
|
END
|
||||||
|
C
|
||||||
|
c-----------------------------------------------------------------------
|
||||||
|
subroutine radlg(mesh,func,r,dx,asum)
|
||||||
|
c-----------------------------------------------------------------------
|
||||||
|
c
|
||||||
|
c simpson's rule integrator for function stored on the
|
||||||
|
c radial logarithmic mesh
|
||||||
|
c
|
||||||
|
c.....logarithmic radial mesh information
|
||||||
|
IMPLICIT REAL*8 (A-H,O-Z)
|
||||||
|
dimension r(mesh)
|
||||||
|
c.....function to be integrated
|
||||||
|
dimension func(mesh)
|
||||||
|
c
|
||||||
|
c.....variable for file = 0
|
||||||
|
c
|
||||||
|
c routine assumes that mesh is an odd number so run check
|
||||||
|
if ( mesh - ( mesh / 2 ) * 2 .ne. 1 ) then
|
||||||
|
write(*,*) '***error in subroutine radlg'
|
||||||
|
write(*,*) 'routine assumes mesh is odd but mesh =',mesh
|
||||||
|
stop
|
||||||
|
endif
|
||||||
|
|
||||||
|
asum = func(1)*r(1)+func(mesh)*r(mesh)
|
||||||
|
do i = 2,mesh-1,2
|
||||||
|
asum = asum + 4.0d0*func(i)*r(i)+2.0d0*func(i+1)*r(i+1)
|
||||||
|
enddo
|
||||||
|
asum = asum*dx/3.0d0
|
||||||
|
return
|
||||||
|
end
|
||||||
|
|
||||||
|
C
|
||||||
|
c-----------------------------------------------------------------------
|
||||||
|
subroutine radlg1(mesh,func,rab,asum)
|
||||||
|
c-----------------------------------------------------------------------
|
||||||
|
c
|
||||||
|
c simpson's rule integrator for function stored on the
|
||||||
|
c radial logarithmic mesh
|
||||||
|
c
|
||||||
|
c.....logarithmic radial mesh information
|
||||||
|
IMPLICIT REAL*8 (A-H,O-Z)
|
||||||
|
dimension rab(mesh)
|
||||||
|
c.....function to be integrated
|
||||||
|
dimension func(mesh)
|
||||||
|
c
|
||||||
|
c.....variable for file = 0
|
||||||
|
c
|
||||||
|
c routine assumes that mesh is an odd number so run check
|
||||||
|
if ( mesh - ( mesh / 2 ) * 2 .ne. 1 ) then
|
||||||
|
write(*,*) '***error in subroutine radlg'
|
||||||
|
write(*,*) 'routine assumes mesh is odd but mesh =',mesh
|
||||||
|
stop
|
||||||
|
endif
|
||||||
|
|
||||||
|
asum = 0.0d0
|
||||||
|
r12 = 1.0d0 / 12.0d0
|
||||||
|
f3 = func(1) * rab(1) * r12
|
||||||
|
c func(1) = 0.0d0
|
||||||
|
|
||||||
|
do 100 i = 2,mesh-1,2
|
||||||
|
f1 = f3
|
||||||
|
f2 = func(i) * rab(i) * r12
|
||||||
|
f3 = func(i+1) * rab(i+1) * r12
|
||||||
|
asum = asum + 5.0d0*f1 + 8.0d0*f2 - 1.0d0*f3
|
||||||
|
c func(i) = asum
|
||||||
|
asum = asum - 1.0d0*f1 + 8.0d0*f2 + 5.0d0*f3
|
||||||
|
c func(i+1) = asum
|
||||||
|
100 continue
|
||||||
|
return
|
||||||
|
end
|
|
@ -0,0 +1,113 @@
|
||||||
|
c
|
||||||
|
c-----------------------------------------------------------------------
|
||||||
|
subroutine ylm_q(lmax,gx,g,ylm)
|
||||||
|
c-----------------------------------------------------------------------
|
||||||
|
c REAL SPHERIcAL HARMONIcS, L IS cOMBINED INDEX FOR LM (L=1,2...25)
|
||||||
|
c ORDER: S, P_X, P_Y, P_Z, D_XY, D_XZ, D_Z^2, D_YZ, D_X^2-Y^2 ....
|
||||||
|
c THE REAL SPHERIcAL HARMONIcS USED HERE FORM BASES FOR THE
|
||||||
|
c IRRIDUcBLE REPRESENTATIONS OF THE gROUP O
|
||||||
|
c
|
||||||
|
c SEE WIESSBLUTH 'ATOMS AND MOLEcULES' PAgES 128-130
|
||||||
|
c ERRORS IN WEISSBLUTH HAVE BEEN cORREcTED:
|
||||||
|
c 1.) ELIMINATION OF THE 7'S FROM L=20
|
||||||
|
c 2.) ADDITION OF THE FAcTOR 1./sqrt(12.) TO L=25
|
||||||
|
c
|
||||||
|
implicit real*8 (A-H,O-Z)
|
||||||
|
dimension ylm(lmax),gx(3)
|
||||||
|
PI=4.D0*DATAN(1.D0)
|
||||||
|
fpi=4.D0*PI
|
||||||
|
eps=1e-9
|
||||||
|
|
||||||
|
if (lmax.ge.26) call errore
|
||||||
|
& (' ylm_q',' not programmed for L>',L)
|
||||||
|
|
||||||
|
|
||||||
|
c note : ylm(q=0) = 1/sqrt(fpi) WHEN L=0 AND = 0 WHEN L>0
|
||||||
|
|
||||||
|
|
||||||
|
ylm(1) = sqrt(1./fpi)
|
||||||
|
|
||||||
|
if(g.lt.eps) then
|
||||||
|
do l=2,lmax
|
||||||
|
ylm(l) = 0.0
|
||||||
|
enddo
|
||||||
|
return
|
||||||
|
endif
|
||||||
|
|
||||||
|
c=sqrt(3./fpi)
|
||||||
|
|
||||||
|
c p_x p_y p_z
|
||||||
|
|
||||||
|
ylm(2) = c*gx(1)/sqrt(g) ! X
|
||||||
|
ylm(3) = c*gx(2)/sqrt(g) ! Y
|
||||||
|
ylm(4) = c*gx(3)/sqrt(g) ! Z
|
||||||
|
|
||||||
|
|
||||||
|
c d_xy d_xz d_yz
|
||||||
|
|
||||||
|
c=sqrt(15./fpi)
|
||||||
|
|
||||||
|
ylm(5) = c*gx(1)*gx(2)/g ! X*Y
|
||||||
|
ylm(6) = c*gx(1)*gx(3)/g ! X*Z
|
||||||
|
ylm(8) = c*gx(2)*gx(3)/g ! Y*Z
|
||||||
|
|
||||||
|
c=sqrt(5.0/fpi/4.0)
|
||||||
|
|
||||||
|
ylm(7) = c*(3.*gx(3)**2/g-1.) ! (3.*Z*Z-1.0)
|
||||||
|
|
||||||
|
c=sqrt(15./fpi/4.)
|
||||||
|
|
||||||
|
ylm(9) = c*(gx(1)**2-gx(2)**2)/g ! X*X-Y*Y
|
||||||
|
|
||||||
|
c=sqrt(7./fpi)*5./2.
|
||||||
|
|
||||||
|
ylm(10) = c*gx(1)*(gx(1)**2-0.6*g)/(g*sqrt(g)) ! X(X^2-3R^2/5)
|
||||||
|
ylm(11) = c*gx(2)*(gx(2)**2-0.6*g)/(g*sqrt(g))
|
||||||
|
|
||||||
|
c=sqrt(7.*15./fpi)
|
||||||
|
|
||||||
|
ylm(12) = c*gx(1)*gx(2)*gx(3)/(g*sqrt(g))
|
||||||
|
|
||||||
|
c=sqrt(7./fpi)*5./2.
|
||||||
|
|
||||||
|
ylm(13) = c*gx(3)*(gx(3)**2-0.6*g)/(g*sqrt(g))
|
||||||
|
|
||||||
|
c=sqrt(7.*15./fpi)/2.
|
||||||
|
|
||||||
|
ylm(14) = c*gx(3)*(gx(1)**2-gx(2)**2)/(g*sqrt(g))
|
||||||
|
ylm(15) = c*gx(2)*(gx(3)**2-gx(1)**2)/(g*sqrt(g))
|
||||||
|
ylm(16) = c*gx(1)*(gx(2)**2-gx(3)**2)/(g*sqrt(g))
|
||||||
|
|
||||||
|
c=sqrt(3.*7./fpi)*5./4.
|
||||||
|
|
||||||
|
ylm(17) = c*((gx(1)**4+gx(2)**4+gx(3)**4)/(g*g)-0.6)
|
||||||
|
|
||||||
|
c=sqrt(9.*35./fpi)/2.
|
||||||
|
|
||||||
|
ylm(18) = c*gx(2)*gx(3)*(gx(2)**2-gx(3)**2)/g**2
|
||||||
|
ylm(19) = c*gx(1)*gx(3)*(gx(3)**2-gx(1)**2)/g**2
|
||||||
|
|
||||||
|
c=sqrt(9.*5./fpi)/4.
|
||||||
|
|
||||||
|
ylm(20) = c*((gx(1)**4-gx(2)**4)-
|
||||||
|
+ 6.*gx(3)**2*(gx(1)**2-gx(2)**2))/(g*g)
|
||||||
|
|
||||||
|
c=sqrt(9.*35./fpi)/2.
|
||||||
|
|
||||||
|
ylm(21) = c*gx(1)*gx(2)*(gx(1)**2-gx(2)**2)/g**2
|
||||||
|
|
||||||
|
|
||||||
|
c=sqrt(9.*5./fpi)*7./2.
|
||||||
|
|
||||||
|
ylm(22) = c*gx(1)*gx(2)*(gx(3)**2-g/7.)/g**2
|
||||||
|
ylm(23) = c*gx(1)*gx(3)*(gx(2)**2-g/7.)/g**2
|
||||||
|
ylm(24) = c*gx(2)*gx(3)*(gx(1)**2-g/7.)/g**2
|
||||||
|
|
||||||
|
c=sqrt(9.*5./fpi/3.)*7./2.
|
||||||
|
|
||||||
|
ylm(25) = c*( gx(3)**4-0.5*(gx(1)**4+gx(2)**4)-
|
||||||
|
+ 6./7.*g*(gx(3)**2-0.5*(gx(1)**2+gx(2)**2) ))/ g**2
|
||||||
|
|
||||||
|
|
||||||
|
return
|
||||||
|
end
|
|
@ -0,0 +1,138 @@
|
||||||
|
c
|
||||||
|
subroutine zgedi(a,lda,n,ipvt,det,work,job)
|
||||||
|
integer lda,n,ipvt(1),job
|
||||||
|
complex*16 a(lda,1),det(2),work(1)
|
||||||
|
c
|
||||||
|
c zgedi computes the determinant and inverse of a matrix
|
||||||
|
c using the factors computed by zgeco or zgefa.
|
||||||
|
c
|
||||||
|
c on entry
|
||||||
|
c
|
||||||
|
c a complex*16(lda, n)
|
||||||
|
c the output from zgeco or zgefa.
|
||||||
|
c
|
||||||
|
c lda integer
|
||||||
|
c the leading dimension of the array a .
|
||||||
|
c
|
||||||
|
c n integer
|
||||||
|
c the order of the matrix a .
|
||||||
|
c
|
||||||
|
c ipvt integer(n)
|
||||||
|
c the pivot vector from zgeco or zgefa.
|
||||||
|
c
|
||||||
|
c work complex*16(n)
|
||||||
|
c work vector. contents destroyed.
|
||||||
|
c
|
||||||
|
c job integer
|
||||||
|
c = 11 both determinant and inverse.
|
||||||
|
c = 01 inverse only.
|
||||||
|
c = 10 determinant only.
|
||||||
|
c
|
||||||
|
c on return
|
||||||
|
c
|
||||||
|
c a inverse of original matrix if requested.
|
||||||
|
c otherwise unchanged.
|
||||||
|
c
|
||||||
|
c det complex*16(2)
|
||||||
|
c determinant of original matrix if requested.
|
||||||
|
c otherwise not referenced.
|
||||||
|
c determinant = det(1) * 10.0**det(2)
|
||||||
|
c with 1.0 .le. cabs1(det(1)) .lt. 10.0
|
||||||
|
c or det(1) .eq. 0.0 .
|
||||||
|
c
|
||||||
|
c error condition
|
||||||
|
c
|
||||||
|
c a division by zero will occur if the input factor contains
|
||||||
|
c a zero on the diagonal and the inverse is requested.
|
||||||
|
c it will not occur if the subroutines are called correctly
|
||||||
|
c and if zgeco has set rcond .gt. 0.0 or zgefa has set
|
||||||
|
c info .eq. 0 .
|
||||||
|
c
|
||||||
|
c linpack. this version dated 08/14/78 .
|
||||||
|
c cleve moler, university of new mexico, argonne national lab.
|
||||||
|
c
|
||||||
|
c subroutines and functions
|
||||||
|
c
|
||||||
|
c blas zaxpy,zscal,zswap
|
||||||
|
c fortran dabs,dcmplx,mod
|
||||||
|
c
|
||||||
|
c internal variables
|
||||||
|
c
|
||||||
|
c
|
||||||
|
complex*16 t
|
||||||
|
double precision ten
|
||||||
|
integer i,j,k,kb,kp1,l,nm1
|
||||||
|
c
|
||||||
|
complex*16 zdum
|
||||||
|
double precision cabs1
|
||||||
|
double precision dreal,dimag
|
||||||
|
complex*16 zdumr,zdumi
|
||||||
|
dreal(zdumr) = zdumr
|
||||||
|
dimag(zdumi) = (0.0d0,-1.0d0)*zdumi
|
||||||
|
cabs1(zdum) = dabs(dreal(zdum)) + dabs(dimag(zdum))
|
||||||
|
c
|
||||||
|
c compute determinant
|
||||||
|
c
|
||||||
|
if (job/10 .eq. 0) go to 70
|
||||||
|
det(1) = (1.0d0,0.0d0)
|
||||||
|
det(2) = (0.0d0,0.0d0)
|
||||||
|
ten = 10.0d0
|
||||||
|
do 50 i = 1, n
|
||||||
|
if (ipvt(i) .ne. i) det(1) = -det(1)
|
||||||
|
det(1) = a(i,i)*det(1)
|
||||||
|
c ...exit
|
||||||
|
if (cabs1(det(1)) .eq. 0.0d0) go to 60
|
||||||
|
10 if (cabs1(det(1)) .ge. 1.0d0) go to 20
|
||||||
|
det(1) = dcmplx(ten,0.0d0)*det(1)
|
||||||
|
det(2) = det(2) - (1.0d0,0.0d0)
|
||||||
|
go to 10
|
||||||
|
20 continue
|
||||||
|
30 if (cabs1(det(1)) .lt. ten) go to 40
|
||||||
|
det(1) = det(1)/dcmplx(ten,0.0d0)
|
||||||
|
det(2) = det(2) + (1.0d0,0.0d0)
|
||||||
|
go to 30
|
||||||
|
40 continue
|
||||||
|
50 continue
|
||||||
|
60 continue
|
||||||
|
70 continue
|
||||||
|
c
|
||||||
|
c compute inverse(u)
|
||||||
|
c
|
||||||
|
if (mod(job,10) .eq. 0) go to 150
|
||||||
|
do 100 k = 1, n
|
||||||
|
a(k,k) = (1.0d0,0.0d0)/a(k,k)
|
||||||
|
t = -a(k,k)
|
||||||
|
call zscal(k-1,t,a(1,k),1)
|
||||||
|
kp1 = k + 1
|
||||||
|
if (n .lt. kp1) go to 90
|
||||||
|
do 80 j = kp1, n
|
||||||
|
t = a(k,j)
|
||||||
|
a(k,j) = (0.0d0,0.0d0)
|
||||||
|
call zaxpy(k,t,a(1,k),1,a(1,j),1)
|
||||||
|
80 continue
|
||||||
|
90 continue
|
||||||
|
100 continue
|
||||||
|
c
|
||||||
|
c form inverse(u)*inverse(l)
|
||||||
|
c
|
||||||
|
nm1 = n - 1
|
||||||
|
if (nm1 .lt. 1) go to 140
|
||||||
|
do 130 kb = 1, nm1
|
||||||
|
k = n - kb
|
||||||
|
kp1 = k + 1
|
||||||
|
do 110 i = kp1, n
|
||||||
|
work(i) = a(i,k)
|
||||||
|
a(i,k) = (0.0d0,0.0d0)
|
||||||
|
110 continue
|
||||||
|
do 120 j = kp1, n
|
||||||
|
t = work(j)
|
||||||
|
call zaxpy(n,t,a(1,j),1,a(1,k),1)
|
||||||
|
120 continue
|
||||||
|
l = ipvt(k)
|
||||||
|
if (l .ne. k) call zswap(n,a(1,k),1,a(1,l),1)
|
||||||
|
130 continue
|
||||||
|
140 continue
|
||||||
|
150 continue
|
||||||
|
return
|
||||||
|
end
|
||||||
|
|
|
@ -0,0 +1,114 @@
|
||||||
|
c
|
||||||
|
subroutine zgefa(a,lda,n,ipvt,info)
|
||||||
|
integer lda,n,ipvt(1),info
|
||||||
|
complex*16 a(lda,1)
|
||||||
|
c
|
||||||
|
c zgefa factors a complex*16 matrix by gaussian elimination.
|
||||||
|
c
|
||||||
|
c zgefa is usually called by zgeco, but it can be called
|
||||||
|
c directly with a saving in time if rcond is not needed.
|
||||||
|
c (time for zgeco) = (1 + 9/n)*(time for zgefa) .
|
||||||
|
c
|
||||||
|
c on entry
|
||||||
|
c
|
||||||
|
c a complex*16(lda, n)
|
||||||
|
c the matrix to be factored.
|
||||||
|
c
|
||||||
|
c lda integer
|
||||||
|
c the leading dimension of the array a .
|
||||||
|
c
|
||||||
|
c n integer
|
||||||
|
c the order of the matrix a .
|
||||||
|
c
|
||||||
|
c on return
|
||||||
|
c
|
||||||
|
c a an upper triangular matrix and the multipliers
|
||||||
|
c which were used to obtain it.
|
||||||
|
c the factorization can be written a = l*u where
|
||||||
|
c l is a product of permutation and unit lower
|
||||||
|
c triangular matrices and u is upper triangular.
|
||||||
|
c
|
||||||
|
c ipvt integer(n)
|
||||||
|
c an integer vector of pivot indices.
|
||||||
|
c
|
||||||
|
c info integer
|
||||||
|
c = 0 normal value.
|
||||||
|
c = k if u(k,k) .eq. 0.0 . this is not an error
|
||||||
|
c condition for this subroutine, but it does
|
||||||
|
c indicate that zgesl or zgedi will divide by zero
|
||||||
|
c if called. use rcond in zgeco for a reliable
|
||||||
|
c indication of singularity.
|
||||||
|
c
|
||||||
|
c linpack. this version dated 08/14/78 .
|
||||||
|
c cleve moler, university of new mexico, argonne national lab.
|
||||||
|
c
|
||||||
|
c subroutines and functions
|
||||||
|
c
|
||||||
|
c blas zaxpy,zscal,izamax
|
||||||
|
c fortran dabs
|
||||||
|
c
|
||||||
|
c internal variables
|
||||||
|
c
|
||||||
|
complex*16 t
|
||||||
|
integer izamax,j,k,kp1,l,nm1
|
||||||
|
c
|
||||||
|
complex*16 zdum
|
||||||
|
double precision cabs1
|
||||||
|
double precision dreal,dimag
|
||||||
|
complex*16 zdumr,zdumi
|
||||||
|
dreal(zdumr) = zdumr
|
||||||
|
dimag(zdumi) = (0.0d0,-1.0d0)*zdumi
|
||||||
|
cabs1(zdum) = dabs(dreal(zdum)) + dabs(dimag(zdum))
|
||||||
|
c
|
||||||
|
c gaussian elimination with partial pivoting
|
||||||
|
c
|
||||||
|
info = 0
|
||||||
|
nm1 = n - 1
|
||||||
|
if (nm1 .lt. 1) go to 70
|
||||||
|
do 60 k = 1, nm1
|
||||||
|
kp1 = k + 1
|
||||||
|
c
|
||||||
|
c find l = pivot index
|
||||||
|
c
|
||||||
|
l = izamax(n-k+1,a(k,k),1) + k - 1
|
||||||
|
ipvt(k) = l
|
||||||
|
c
|
||||||
|
c zero pivot implies this column already triangularized
|
||||||
|
c
|
||||||
|
if (cabs1(a(l,k)) .eq. 0.0d0) go to 40
|
||||||
|
c
|
||||||
|
c interchange if necessary
|
||||||
|
c
|
||||||
|
if (l .eq. k) go to 10
|
||||||
|
t = a(l,k)
|
||||||
|
a(l,k) = a(k,k)
|
||||||
|
a(k,k) = t
|
||||||
|
10 continue
|
||||||
|
c
|
||||||
|
c compute multipliers
|
||||||
|
c
|
||||||
|
t = -(1.0d0,0.0d0)/a(k,k)
|
||||||
|
call zscal(n-k,t,a(k+1,k),1)
|
||||||
|
c
|
||||||
|
c row elimination with column indexing
|
||||||
|
c
|
||||||
|
do 30 j = kp1, n
|
||||||
|
t = a(l,j)
|
||||||
|
if (l .eq. k) go to 20
|
||||||
|
a(l,j) = a(k,j)
|
||||||
|
a(k,j) = t
|
||||||
|
20 continue
|
||||||
|
call zaxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
|
||||||
|
30 continue
|
||||||
|
go to 50
|
||||||
|
40 continue
|
||||||
|
info = k
|
||||||
|
50 continue
|
||||||
|
60 continue
|
||||||
|
70 continue
|
||||||
|
ipvt(n) = n
|
||||||
|
if (cabs1(a(n,n)) .eq. 0.0d0) info = n
|
||||||
|
return
|
||||||
|
end
|
||||||
|
|
||||||
|
|
|
@ -126,6 +126,7 @@ implicit none
|
||||||
|
|
||||||
!
|
!
|
||||||
call c_bands (iter, ik_, dr2)
|
call c_bands (iter, ik_, dr2)
|
||||||
|
|
||||||
!
|
!
|
||||||
!! skip all the rest if not lscf
|
!! skip all the rest if not lscf
|
||||||
if (.not.lscf) then
|
if (.not.lscf) then
|
||||||
|
@ -143,6 +144,12 @@ implicit none
|
||||||
write (6, 9020) (xk (i, ik), i = 1, 3)
|
write (6, 9020) (xk (i, ik), i = 1, 3)
|
||||||
write (6, 9030) (et (ibnd, ik) * 13.6058, ibnd = 1, nbnd)
|
write (6, 9030) (et (ibnd, ik) * 13.6058, ibnd = 1, nbnd)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
|
! Do a Berry phase polarization calculation if required
|
||||||
|
if ((lberry).and.(iswitch /= -1)) call errore('electrons', &
|
||||||
|
'calculation=''nscf'' is mandatory if lberry=''.true.''',1)
|
||||||
|
if (lberry) call c_phase
|
||||||
|
|
||||||
! jump to the end
|
! jump to the end
|
||||||
goto 999
|
goto 999
|
||||||
|
|
||||||
|
|
12
PW/input.f90
12
PW/input.f90
|
@ -27,7 +27,8 @@ subroutine iosys
|
||||||
lmovecell, imix, at, omega, ityp, tau, nks, xk, wk, uakbar, amconv, &
|
lmovecell, imix, at, omega, ityp, tau, nks, xk, wk, uakbar, amconv, &
|
||||||
force, at_old, omega_old, starting_scf_threshold, title, crystal, &
|
force, at_old, omega_old, starting_scf_threshold, title, crystal, &
|
||||||
atm, nk1, nk2, nk3, k1, k2, k3, &
|
atm, nk1, nk2, nk3, k1, k2, k3, &
|
||||||
tefield, edir, emaxpos, eopreg, eamp
|
tefield, edir, emaxpos, eopreg, eamp, &
|
||||||
|
lberry, gdir, nppstr
|
||||||
use io, only : tmp_dir, prefix, pseudo_dir, pseudop
|
use io, only : tmp_dir, prefix, pseudo_dir, pseudop
|
||||||
use constants, only: pi
|
use constants, only: pi
|
||||||
#ifdef __PARA
|
#ifdef __PARA
|
||||||
|
@ -54,7 +55,8 @@ subroutine iosys
|
||||||
NAMELIST / control / title, calculation, verbosity, &
|
NAMELIST / control / title, calculation, verbosity, &
|
||||||
restart_mode, nstep, iprint, isave, tstress, tprnfor, &
|
restart_mode, nstep, iprint, isave, tstress, tprnfor, &
|
||||||
dt, ndr, ndw, outdir, prefix, max_seconds, ekin_conv_thr,&
|
dt, ndr, ndw, outdir, prefix, max_seconds, ekin_conv_thr,&
|
||||||
etot_conv_thr, forc_conv_thr, pseudo_dir, disk_io, tefield
|
etot_conv_thr, forc_conv_thr, pseudo_dir, disk_io, tefield, &
|
||||||
|
lberry, gdir, nppstr
|
||||||
|
|
||||||
! SYSTEM namelist
|
! SYSTEM namelist
|
||||||
|
|
||||||
|
@ -159,6 +161,9 @@ subroutine iosys
|
||||||
disk_io = 'default'
|
disk_io = 'default'
|
||||||
tefield=.false.
|
tefield=.false.
|
||||||
noinv = .false. ! not actually used
|
noinv = .false. ! not actually used
|
||||||
|
lberry=.false.
|
||||||
|
gdir=0
|
||||||
|
nppstr=0
|
||||||
!
|
!
|
||||||
#ifdef __T3E
|
#ifdef __T3E
|
||||||
call pxfgetenv('HOME',0,pseudo_dir,i,ios)
|
call pxfgetenv('HOME',0,pseudo_dir,i,ios)
|
||||||
|
@ -390,6 +395,9 @@ subroutine iosys
|
||||||
CALL mp_bcast( pseudo_dir, ionode_id )
|
CALL mp_bcast( pseudo_dir, ionode_id )
|
||||||
CALL mp_bcast( disk_io, ionode_id )
|
CALL mp_bcast( disk_io, ionode_id )
|
||||||
CALL mp_bcast( tefield, ionode_id )
|
CALL mp_bcast( tefield, ionode_id )
|
||||||
|
CALL mp_bcast( lberry, ionode_id )
|
||||||
|
CALL mp_bcast( gdir, ionode_id )
|
||||||
|
CALL mp_bcast( nppstr, ionode_id )
|
||||||
!
|
!
|
||||||
! ... SYSTEM Variables Broadcast
|
! ... SYSTEM Variables Broadcast
|
||||||
!
|
!
|
||||||
|
|
13
PW/pwcom.f90
13
PW/pwcom.f90
|
@ -604,6 +604,18 @@ module sticks
|
||||||
! potential grid, and its wave functions sub-grid.
|
! potential grid, and its wave functions sub-grid.
|
||||||
end module
|
end module
|
||||||
|
|
||||||
|
module bp
|
||||||
|
use parameters
|
||||||
|
!
|
||||||
|
! The variables needed for the Berry phase polarization calculation
|
||||||
|
!
|
||||||
|
logical :: &
|
||||||
|
lberry ! if true, calculate polarization
|
||||||
|
integer :: &
|
||||||
|
gdir, & ! G-vector for polarization calculation
|
||||||
|
nppstr ! number of k-points (parallel vector)
|
||||||
|
end module bp
|
||||||
|
!
|
||||||
|
|
||||||
|
|
||||||
module pwcom
|
module pwcom
|
||||||
|
@ -636,5 +648,6 @@ module pwcom
|
||||||
use ldaU
|
use ldaU
|
||||||
use extfield
|
use extfield
|
||||||
use sticks
|
use sticks
|
||||||
|
use bp
|
||||||
end module pwcom
|
end module pwcom
|
||||||
!
|
!
|
||||||
|
|
|
@ -237,7 +237,12 @@ subroutine setup
|
||||||
call setupkpoint (s, nrot, xk, wk, nks, npk, nk1, &
|
call setupkpoint (s, nrot, xk, wk, nks, npk, nk1, &
|
||||||
nk2, nk3, k1, k2, k3, at, bg, tipo)
|
nk2, nk3, k1, k2, k3, at, bg, tipo)
|
||||||
else if (nks == 0) then
|
else if (nks == 0) then
|
||||||
call kpoint_grid ( nrot, s, bg, npk, k1,k2,k3, nk1,nk2,nk3, nks, xk, wk)
|
if (lberry) then
|
||||||
|
call kp_strings &
|
||||||
|
( nppstr, gdir, nrot, s, bg, npk, k1,k2,k3, nk1,nk2,nk3, nks, xk, wk)
|
||||||
|
else
|
||||||
|
call kpoint_grid ( nrot, s, bg, npk, k1,k2,k3, nk1,nk2,nk3, nks, xk, wk)
|
||||||
|
end if
|
||||||
end if
|
end if
|
||||||
!
|
!
|
||||||
input_nks = nks
|
input_nks = nks
|
||||||
|
|
Loading…
Reference in New Issue