Skip to content
GitLab
Projects
Groups
Snippets
Help
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
fleur
Project overview
Project overview
Details
Activity
Releases
Repository
Repository
Files
Commits
Branches
Tags
Contributors
Graph
Compare
Issues
52
Issues
52
List
Boards
Labels
Service Desk
Milestones
Operations
Operations
Incidents
Packages & Registries
Packages & Registries
Container Registry
Analytics
Analytics
Repository
Value Stream
Wiki
Wiki
Members
Members
Collapse sidebar
Close sidebar
Activity
Graph
Create a new issue
Commits
Issue Boards
Open sidebar
fleur
fleur
Commits
567d369d
Commit
567d369d
authored
Jul 15, 2019
by
Matthias Redies
Browse files
Options
Browse Files
Download
Email Patches
Plain Diff
try to remove every 1d0-like expression
parent
bc50dedf
Changes
18
Expand all
Hide whitespace changes
Inline
Side-by-side
Showing
18 changed files
with
237 additions
and
237 deletions
+237
-237
hybrid/HF_init.F90
hybrid/HF_init.F90
+1
-1
hybrid/checkolap.F90
hybrid/checkolap.F90
+2
-2
hybrid/coulombmatrix.F90
hybrid/coulombmatrix.F90
+24
-24
hybrid/exchange_core.F90
hybrid/exchange_core.F90
+3
-3
hybrid/exchange_val_hf.F90
hybrid/exchange_val_hf.F90
+4
-4
hybrid/gen_wavf.F90
hybrid/gen_wavf.F90
+1
-1
hybrid/hf_setup.F90
hybrid/hf_setup.F90
+1
-1
hybrid/hsefunctional.F90
hybrid/hsefunctional.F90
+3
-3
hybrid/kp_perturbation.F90
hybrid/kp_perturbation.F90
+16
-16
hybrid/mixedbasis.F90
hybrid/mixedbasis.F90
+1
-1
hybrid/olap.F90
hybrid/olap.F90
+6
-6
hybrid/spmvec.F90
hybrid/spmvec.F90
+11
-11
hybrid/subvxc.F90
hybrid/subvxc.F90
+3
-3
hybrid/symm_hf.F90
hybrid/symm_hf.F90
+1
-1
hybrid/symmetrizeh.F90
hybrid/symmetrizeh.F90
+11
-11
hybrid/trafo.F90
hybrid/trafo.F90
+26
-26
hybrid/wavefproducts.F90
hybrid/wavefproducts.F90
+21
-21
hybrid/wrapper.F
hybrid/wrapper.F
+102
-102
No files found.
hybrid/HF_init.F90
View file @
567d369d
...
...
@@ -58,7 +58,7 @@ CONTAINS
hybdat
%
sfac
(
0
)
=
1
DO
i
=
1
,
hybdat
%
maxfac
hybdat
%
fac
(
i
)
=
hybdat
%
fac
(
i
-1
)
*
i
! hybdat%fac(i) = i!
hybdat
%
sfac
(
i
)
=
hybdat
%
sfac
(
i
-1
)
*
sqrt
(
i
*
1
d
0
)
! hybdat%sfac(i) = sqrt(i!)
hybdat
%
sfac
(
i
)
=
hybdat
%
sfac
(
i
-1
)
*
sqrt
(
i
*
1
.
0
)
! hybdat%sfac(i) = sqrt(i!)
END
DO
...
...
hybrid/checkolap.F90
View file @
567d369d
...
...
@@ -47,7 +47,7 @@
REAL
::
qnorm
COMPLEX
::
cexp
,
cdum
COMPLEX
,
PARAMETER
::
img
=
(
0
d0
,
1d
0
)
COMPLEX
,
PARAMETER
::
img
=
(
0
.0
,
1.
0
)
! -local arrays -
INTEGER
::
iarr
(
2
),
gpt
(
3
)
...
...
@@ -315,7 +315,7 @@
! & ikpt,sum(rarr(:1)**2/nbands(ikpt)),maxval(rarr(:1))
! CALL writeout(outtext,mpi%irank)
! IF( iatom .eq. 6 ) THEN
! cdum = exp(2*pi*img*dot_product(bkf(:,ikpt),(/0
d0,0d0,1d
0/) ))
! cdum = exp(2*pi*img*dot_product(bkf(:,ikpt),(/0
.0,0.0,1.
0/) ))
! lm = 0
! DO l = 0,lmax(itype)
! DO m = -l,l
...
...
hybrid/coulombmatrix.F90
View file @
567d369d
...
...
@@ -146,7 +146,7 @@ CONTAINS
facC
(
-1
:
0
)
=
1
! facA(i) = i!
DO
i
=
1
,
maxfac
! facB(i) = sqrt(i!)
facA
(
i
)
=
facA
(
i
-
1
)
*
i
! facC(i) = (2i+1)!!
facB
(
i
)
=
facB
(
i
-
1
)
*
SQRT
(
i
*
1
d
0
)
!
facB
(
i
)
=
facB
(
i
-
1
)
*
SQRT
(
i
*
1
.
0
)
!
facC
(
i
)
=
facC
(
i
-
1
)
*
(
2
*
i
+
1
)
!
END
DO
...
...
@@ -172,7 +172,7 @@ CONTAINS
coulomb
=
0
call
timestop
(
"coulomb allocation"
)
IF
(
mpi
%
irank
==
0
)
WRITE
(
6
,
'(/A,F6.1," MB")'
)
'Size of coulomb matrix:'
,
16
d
0
/
1048576
*
SIZE
(
coulomb
)
IF
(
mpi
%
irank
==
0
)
WRITE
(
6
,
'(/A,F6.1," MB")'
)
'Size of coulomb matrix:'
,
16
.
0
/
1048576
*
SIZE
(
coulomb
)
! Generate Symmetry:
! Reduce list of g-Points so that only one of each symm-equivalent is calculated
...
...
@@ -223,7 +223,7 @@ CONTAINS
END
IF
IF
(
ANY
(
sym
%
tau
(:,
isym2
)
/
=
0
))
CYCLE
IF
(
ALL
(
ABS
(
MATMUL
(
rrot
(:,
:,
isym
),
kpts
%
bk
(:,
ikpt
))
-
kpts
%
bk
(:,
ikpt
))
.LT.
1
d
-12
))
THEN
IF
(
ALL
(
ABS
(
MATMUL
(
rrot
(:,
:,
isym
),
kpts
%
bk
(:,
ikpt
))
-
kpts
%
bk
(:,
ikpt
))
.LT.
1
0.0
**
-12
))
THEN
isym1
=
isym1
+
1
sym1
(
isym1
,
ikpt
)
=
isym
END
IF
...
...
@@ -300,7 +300,7 @@ CONTAINS
IF
(
lm2
.GT.
lm1
)
EXIT
lp1
! Don't cross the diagonal!
gmat
(
lm1
,
lm2
)
=
facB
(
l1
+
l2
+
m2
-
m1
)
*
facB
(
l1
+
l2
+
m1
-
m2
)/
&
(
facB
(
l1
+
m1
)
*
facB
(
l1
-
m1
)
*
facB
(
l2
+
m2
)
*
facB
(
l2
-
m2
))/
&
SQRT
(
1
d0
*
(
2
*
l1
+
1
)
*
(
2
*
l2
+
1
)
*
(
2
*
(
l1
+
l2
)
+
1
))
*
(
4
*
pi_const
)
**
1.5d0
SQRT
(
1
.0
*
(
2
*
l1
+
1
)
*
(
2
*
l2
+
1
)
*
(
2
*
(
l1
+
l2
)
+
1
))
*
(
4
*
pi_const
)
**
1.5
gmat
(
lm2
,
lm1
)
=
gmat
(
lm1
,
lm2
)
END
DO
END
DO
LP1
...
...
@@ -580,7 +580,7 @@ CONTAINS
q
=
MATMUL
(
kpts
%
bk
(:,
ikpt
)
+
hybrid
%
gptm
(:,
igptp
),
cell
%
bmat
)
qnorm
=
SQRT
(
SUM
(
q
**
2
))
iqnrm
=
pqnrm
(
igpt
,
ikpt
)
IF
(
ABS
(
qnrm
(
iqnrm
)
-
qnorm
)
.GT.
1
d
-12
)
then
IF
(
ABS
(
qnrm
(
iqnrm
)
-
qnorm
)
.GT.
1
0.0
**
-12
)
then
STOP
'coulombmatrix: qnorm does not equal corresponding & element in qnrm (bug?)'
! We shouldn't stop here!
endif
...
...
@@ -896,7 +896,7 @@ CONTAINS
! Add corrections from higher orders in (3b) to coulomb(:,1)
! (1) igpt1 > 1 , igpt2 > 1 (finite G vectors)
call
timestart
(
"add corrections from higher orders"
)
rdum
=
(
4
*
pi_const
)
**
(
1.5
d0
)/
cell
%
vol
**
2
*
gmat
(
1
,
1
)
rdum
=
(
4
*
pi_const
)
**
(
1.5
)/
cell
%
vol
**
2
*
gmat
(
1
,
1
)
DO
igpt0
=
1
,
hybrid
%
ngptm1
(
1
)
igpt2
=
hybrid
%
pgptm1
(
igpt0
,
1
);
IF
(
igpt2
==
1
)
CYCLE
ix
=
hybrid
%
nbasp
+
igpt2
...
...
@@ -1692,8 +1692,8 @@ CONTAINS
+
CONJG
(
claplace
(
i
))
*
coeff
(
j
))/
2
)
END
DO
END
DO
coeff
(
hybrid
%
nbasp
+
1
)
=
1
d
0
coeff
(
hybrid
%
nbasp
+
2
:)
=
0
d
0
coeff
(
hybrid
%
nbasp
+
1
)
=
1
.
0
coeff
(
hybrid
%
nbasp
+
2
:)
=
0
.
0
IF
(
sym
%
invs
)
THEN
CALL
desymmetrize
(
coeff
,
1
,
nbasm1
(
1
),
2
,
&
...
...
@@ -1724,7 +1724,7 @@ CONTAINS
!
! Convergence parameter
#define CONVPARAM 1
d
-18
#define CONVPARAM 1
0.0**
-18
! Do some additional shells ( real-space and Fourier-space sum )
#define ADDSHELL1 40
#define ADDSHELL2 0
...
...
@@ -1780,7 +1780,7 @@ CONTAINS
IF
(
mpi
%
irank
/
=
0
)
first
=
.FALSE.
rdum
=
cell
%
vol
**
(
1
d
0
/
3
)
! define "average lattice parameter"
rdum
=
cell
%
vol
**
(
1
.
0
/
3
)
! define "average lattice parameter"
! ewaldlambda = ewaldscale
scale
=
hybrid
%
ewaldlambda
/
rdum
...
...
@@ -1877,7 +1877,7 @@ CONTAINS
radsh
(
i
)
=
a
END
DO
call
timestart
(
"rorderpf"
)
CALL
rorderpf
(
pnt
,
radsh
,
nptsh
,
MAX
(
0
,
INT
(
LOG
(
nptsh
*
0.001
d0
)/
LOG
(
2d
0
))))
CALL
rorderpf
(
pnt
,
radsh
,
nptsh
,
MAX
(
0
,
INT
(
LOG
(
nptsh
*
0.001
)/
LOG
(
2.
0
))))
call
timestop
(
"rorderpf"
)
ptsh
=
ptsh
(:,
pnt
)
radsh
=
radsh
(
pnt
)
...
...
@@ -1889,13 +1889,13 @@ CONTAINS
DO
i
=
1
,
nptsh
IF
(
ALL
(
conv
.NE.
HUGE
(
i
)))
EXIT
IF
(
i
.NE.
1
)
THEN
IF
(
ABS
(
radsh
(
i
)
-
radsh
(
i
-
1
))
.GT.
1
d
-10
)
ishell
=
ishell
+
1
IF
(
ABS
(
radsh
(
i
)
-
radsh
(
i
-
1
))
.GT.
1
0.0
**
-10
)
ishell
=
ishell
+
1
ENDIF
ra
=
MATMUL
(
cell
%
amat
,
ptsh
(:,
i
))
+
rc
a
=
scale
*
SQRT
(
SUM
(
ra
**
2
))
IF
(
a
.EQ.
0
)
THEN
CYCLE
ELSE
IF
(
ABS
(
a
-
a1
)
.GT.
1
d
-10
)
THEN
ELSE
IF
(
ABS
(
a
-
a1
)
.GT.
1
0.0
**
-10
)
THEN
a1
=
a
rexp
=
EXP
(
-
a
)
IF
(
ishell
.LE.
conv
(
0
))
g
(
0
)
=
rexp
/
a
&
...
...
@@ -1981,13 +1981,13 @@ CONTAINS
conv
=
HUGE
(
i
)
DO
i
=
1
,
nptsh
IF
(
i
.GT.
1
)
THEN
IF
(
ABS
(
radsh
(
i
)
-
radsh
(
i
-
1
))
.GT.
1
d
-10
)
ishell
=
ishell
+
1
IF
(
ABS
(
radsh
(
i
)
-
radsh
(
i
-
1
))
.GT.
1
0.0
**
-10
)
ishell
=
ishell
+
1
ENDIF
ki
=
ptsh
(:,
i
)
+
k
-
NINT
(
k
)
! -nint(...) transforms to Wigner-Seitz cell ( i.e. -0.5 <= x,y,z < 0.5 )
ka
=
MATMUL
(
ki
,
cell
%
bmat
)
a
=
SQRT
(
SUM
(
ka
**
2
))/
scale
aa
=
(
1
+
a
**
2
)
**
(
-1
)
IF
(
ABS
(
a
-
a1
)
.GT.
1
d
-10
)
THEN
IF
(
ABS
(
a
-
a1
)
.GT.
1
0.0
**
-10
)
THEN
a1
=
a
IF
(
a
.EQ.
0
)
THEN
g
(
0
)
=
pref
*
(
-4
)
...
...
@@ -2013,7 +2013,7 @@ CONTAINS
call
timestart
(
"harmonics"
)
call
ylm4
(
maxl
,
ka
,
y
)
call
timestop
(
"harmonics"
)
cdum
=
1
d
0
cdum
=
1
.
0
lm
=
0
DO
l
=
0
,
maxl
IF
(
ishell
.LE.
conv
(
l
))
THEN
...
...
@@ -2046,7 +2046,7 @@ CONTAINS
!
! Add contribution for l=0 to diagonal elements and rescale structure constants
!
structconst
(
1
,
1
,
1
,
:)
=
structconst
(
1
,
1
,
1
,
:)
-
5
d
0
/
16
/
SQRT
(
4
*
pi_const
)
structconst
(
1
,
1
,
1
,
:)
=
structconst
(
1
,
1
,
1
,
:)
-
5
.
0
/
16
/
SQRT
(
4
*
pi_const
)
DO
i
=
2
,
atoms
%
nat
structconst
(:,
i
,
i
,
:)
=
structconst
(:,
1
,
1
,
:)
END
DO
...
...
@@ -2054,11 +2054,11 @@ CONTAINS
structconst
(
l
**
2
+
1
:(
l
+
1
)
**
2
,
:,
:,
:)
=
structconst
(
l
**
2
+
1
:(
l
+
1
)
**
2
,
:,
:,
:)
*
scale
**
(
l
+
1
)
END
DO
rad
=
(
cell
%
vol
*
3
/
4
/
pi_const
)
**
(
1
d
0
/
3
)
! Wigner-Seitz radius (rad is recycled)
rad
=
(
cell
%
vol
*
3
/
4
/
pi_const
)
**
(
1
.
0
/
3
)
! Wigner-Seitz radius (rad is recycled)
! Calculate accuracy of Gamma-decomposition
IF
(
ALL
(
kpts
%
bk
.EQ.
0
))
GOTO
4
a
=
1
d
30
! ikpt = index of shortest non-zero k-point
a
=
1
0.0
**
30
! ikpt = index of shortest non-zero k-point
DO
i
=
2
,
kpts
%
nkpt
rdum
=
SUM
(
MATMUL
(
kpts
%
bk
(:,
i
),
cell
%
bmat
)
**
2
)
IF
(
rdum
.LT.
a
)
THEN
...
...
@@ -2163,12 +2163,12 @@ CONTAINS
ptsh
=
ihelp
DEALLOCATE
(
rhelp
,
ihelp
)
CALL
rorderpf
(
pnt
,
radsh
,
nptsh
,
MAX
(
0
,
INT
(
LOG
(
nptsh
*
0.001
d0
)/
LOG
(
2d
0
))))
CALL
rorderpf
(
pnt
,
radsh
,
nptsh
,
MAX
(
0
,
INT
(
LOG
(
nptsh
*
0.001
)/
LOG
(
2.
0
))))
radsh
=
radsh
(
pnt
)
ptsh
=
ptsh
(:,
pnt
)
nshell
=
1
DO
i
=
2
,
nptsh
IF
(
radsh
(
i
)
-
radsh
(
i
-
1
)
.GT.
1
d
-10
)
nshell
=
nshell
+
1
IF
(
radsh
(
i
)
-
radsh
(
i
-
1
)
.GT.
1
0.0
**
-10
)
nshell
=
nshell
+
1
END
DO
IF
(
lwrite
)
&
...
...
@@ -2203,7 +2203,7 @@ CONTAINS
q
=
MATMUL
(
kpts
%
bk
(:,
ikpt
)
+
gpt
(:,
igptp
),
cell
%
bmat
)
qnorm
=
SQRT
(
SUM
(
q
**
2
))
DO
j
=
1
,
i
IF
(
ABS
(
qnrm
(
j
)
-
qnorm
)
.LT.
1
d
-12
)
THEN
IF
(
ABS
(
qnrm
(
j
)
-
qnorm
)
.LT.
1
0.0
**
-12
)
THEN
pqnrm
(
igpt
,
ikpt
)
=
j
CYCLE
igptloop
END
IF
...
...
@@ -2292,7 +2292,7 @@ CONTAINS
r2
=
MIN
(
ABS
(
db
/
b1
),
ABS
(
dc
/
c1
))
! Ensure numerical stability. If both formulas are not sufficiently stable, the program stops.
IF
(
r1
.GT.
r2
)
THEN
IF
(
r1
.LT.
1
d
-6
.AND.
l_warn
)
THEN
IF
(
r1
.LT.
1
0.0
**
-6
.AND.
l_warn
)
THEN
WRITE
(
6
,
'(A,E12.5,A,E12.5,A)'
)
'sphbessel_integral: Warning! Formula One possibly unstable. Ratios:'
,
&
r1
,
'('
,
r2
,
')'
WRITE
(
6
,
'(A,2F15.10,I4)'
)
' Current qnorms and atom type:'
,
q1
,
q2
,
itype
...
...
@@ -2300,7 +2300,7 @@ CONTAINS
END
IF
sphbessel_integral
=
s
**
3
/
dq
*
da
ELSE
IF
(
r2
.LT.
1
d
-6
.AND.
l_warn
)
THEN
IF
(
r2
.LT.
1
0.0
**
-6
.AND.
l_warn
)
THEN
WRITE
(
6
,
'(A,E13.5,A,E13.5,A)'
)
'sphbessel_integral: Warning! Formula Two possibly unstable. Ratios:'
,
&
r2
,
'('
,
r1
,
')'
WRITE
(
6
,
'(A,2F15.10,I4)'
)
' Current qnorms and atom type:'
,
&
...
...
hybrid/exchange_core.F90
View file @
567d369d
...
...
@@ -203,7 +203,7 @@ CONTAINS
IF
(
l_real
)
THEN
IF
(
ANY
(
ABS
(
AIMAG
(
exchange
))
.GT.
1
d
-10
)
)
THEN
IF
(
ANY
(
ABS
(
AIMAG
(
exchange
))
.GT.
1
0.0
**
-10
)
)
THEN
IF
(
mpi
%
irank
==
0
)
WRITE
(
6
,
'(A)'
)
'exchangeCore: Warning! Unusually large imaginary component.'
WRITE
(
*
,
*
)
MAXVAL
(
ABS
(
AIMAG
(
exchange
)))
STOP
'exchangeCore: Unusually large imaginary component.'
...
...
@@ -394,7 +394,7 @@ CONTAINS
IF
(
mat_ex
%
l_real
)
THEN
IF
(
ANY
(
ABS
(
AIMAG
(
exchange
))
.GT.
1
d
-10
)
)
THEN
IF
(
ANY
(
ABS
(
AIMAG
(
exchange
))
.GT.
1
0.0
**
-10
)
)
THEN
IF
(
mpi
%
irank
==
0
)
WRITE
(
6
,
'(A)'
)
'exchangeCore: Warning! Unusually large imaginary component.'
WRITE
(
*
,
*
)
MAXVAL
(
ABS
(
AIMAG
(
exchange
)))
STOP
'exchangeCore: Unusually large imaginary component.'
...
...
@@ -600,7 +600,7 @@ CONTAINS
REAL
::
rdum0
,
rdum1
,
rdum2
,
rdum3
,
rdum4
COMPLEX
::
cdum
COMPLEX
,
PARAMETER
::
img
=
(
0
d0
,
1d
0
)
COMPLEX
,
PARAMETER
::
img
=
(
0
.0
,
1.
0
)
! - local arrays -
INTEGER
::
point
(
hybdat
%
maxindxc
,
-
hybdat
%
lmaxcd
:
hybdat
%
lmaxcd
,
0
:
hybdat
%
lmaxcd
,
atoms
%
nat
)
INTEGER
::
lmstart
(
0
:
atoms
%
lmaxd
,
atoms
%
ntype
)
...
...
hybrid/exchange_val_hf.F90
View file @
567d369d
...
...
@@ -44,7 +44,7 @@
! the "tail" is
! vol/(8*pi**3) INT F(k) d^3k - P SUM(k) F(k) ( P = principal value ) .
! For expo->0 the two terms diverge. Therefore a cutoff radius q0 is introduced and related to
! expo by exp(-expo*q0**2)=delta ( delta = small value, e.g., delta = 1
d
-10 ) .
! expo by exp(-expo*q0**2)=delta ( delta = small value, e.g., delta = 1
*10.0**
-10 ) .
! The resulting formula
! vol/(4*pi**1.5*sqrt(expo)) * erf(sqrt(a)*q0) - sum(q,0<q<q0) exp(-expo*q**2)/q**2
! converges well with q0. (Should be the default.)
...
...
@@ -320,7 +320,7 @@ SUBROUTINE exchange_valence_hf(nk,kpts,nkpt_EIBZ,sym,atoms,hybrid,cell,dimension
cprod_vv_c
(:
hybrid
%
nbasm
(
ikpt0
),:,:)
=
carr3_vv_c
(:
hybrid
%
nbasm
(
ikpt0
),:,:)
ENDIF
ELSE
phase_vv
(:,:)
=
(
1
d0
,
0d
0
)
phase_vv
(:,:)
=
(
1
.0
,
0.
0
)
END
IF
! calculate exchange matrix at ikpt0
...
...
@@ -519,8 +519,8 @@ SUBROUTINE calc_divergence(cell,kpts,divergence)
REAL
::
expo
,
rrad
,
k
(
3
),
kv1
(
3
),
kv2
(
3
),
kv3
(
3
),
knorm2
COMPLEX
::
cdum
expo
=
5
d
-3
rrad
=
sqrt
(
-
log
(
5
d
-3
)/
expo
)
expo
=
5
*
10.0
**
-3
rrad
=
sqrt
(
-
log
(
5
*
10.0
**
-3
)/
expo
)
cdum
=
sqrt
(
expo
)
*
rrad
divergence
=
cell
%
omtil
/
(
tpi_const
**
2
)
*
sqrt
(
pi_const
/
expo
)
*
cerf
(
cdum
)
rrad
=
rrad
**
2
...
...
hybrid/gen_wavf.F90
View file @
567d369d
...
...
@@ -59,7 +59,7 @@ CONTAINS
TYPE
(
t_mat
)
::
zhlp
INTEGER
::
ikpt0
,
ikpt
,
itype
,
iop
,
ispin
,
ieq
,
indx
,
iatom
INTEGER
::
i
,
j
,
l
,
ll
,
lm
,
ng
,
ok
COMPLEX
::
img
=
(
0
d0
,
1d
0
)
COMPLEX
::
img
=
(
0
.0
,
1.
0
)
INTEGER
::
nodem
,
noded
REAL
::
wronk
...
...
hybrid/hf_setup.F90
View file @
567d369d
...
...
@@ -150,7 +150,7 @@ SUBROUTINE hf_setup(hybrid,input,sym,kpts,DIMENSION,atoms,mpi,noco,cell,oneD,res
END
DO
DO
i
=
1
,
hybrid
%
ne_eig
(
nk
)
IF
(
results
%
w_iks
(
i
,
nk
,
jsp
)
.GT.
0
d
0
)
hybrid
%
nobd
(
nk
)
=
hybrid
%
nobd
(
nk
)
+
1
IF
(
results
%
w_iks
(
i
,
nk
,
jsp
)
.GT.
0
.
0
)
hybrid
%
nobd
(
nk
)
=
hybrid
%
nobd
(
nk
)
+
1
END
DO
IF
(
hybrid
%
nobd
(
nk
)
>
hybrid
%
nbands
(
nk
))
THEN
...
...
hybrid/hsefunctional.F90
View file @
567d369d
...
...
@@ -1013,7 +1013,7 @@ CONTAINS
! REAL, PARAMETER :: omega = omega_VHSE() ! omega of the HSE functional
! REAL, PARAMETER :: r1_omega2 = 1.0 / omega**2 ! 1/omega^2
! REAL, PARAMETER :: r1_4omega2 = 0.25 * r1_omega2 ! 1/(4 omega^2)
COMPLEX
,
PARAMETER
::
img
=
(
0
d0
,
1d
0
)
! i
COMPLEX
,
PARAMETER
::
img
=
(
0
.0
,
1.
0
)
! i
! private arrays
INTEGER
::
gPts
(
3
,
noGPts
)
! g vectors (internal units)
...
...
@@ -1349,7 +1349,7 @@ CONTAINS
! REAL, PARAMETER :: omega = omega_VHSE() ! omega of the HSE functional
! REAL, PARAMETER :: r1_omega2 = 1.0 / omega**2 ! 1/omega^2
! REAL, PARAMETER :: r1_4omega2 = 0.25 * r1_omega2 ! 1/(4 omega^2)
COMPLEX
,
PARAMETER
::
img
=
(
0
d0
,
1d
0
)
! i
COMPLEX
,
PARAMETER
::
img
=
(
0
.0
,
1.
0
)
! i
! private arrays
INTEGER
::
gPts
(
3
,
noGPts
)
! g vectors (internal units)
...
...
@@ -2326,7 +2326,7 @@ CONTAINS
END
DO
#ifdef CPP_INVERSION
IF
(
ANY
(
ABS
(
aimag
(
exchange
))
.GT.
1
d
-10
)
)
THEN
IF
(
ANY
(
ABS
(
aimag
(
exchange
))
.GT.
1
0.0
**
-10
)
)
THEN
IF
(
irank
==
0
)
WRITE
(
6
,
'(A)'
)
'exchangeCore: Warning! Unusually large imaginary component.'
WRITE
(
*
,
*
)
MAXVAL
(
ABS
(
aimag
(
exchange
)))
STOP
'exchangeCore: Unusually large imaginary component.'
...
...
hybrid/kp_perturbation.F90
View file @
567d369d
...
...
@@ -65,7 +65,7 @@
COMPLEX
::
cj
,
cdj
COMPLEX
::
denom
,
enum
COMPLEX
::
cdum
,
cdum1
,
cdum2
COMPLEX
,
PARAMETER
::
img
=
(
0
d0
,
1d
0
)
COMPLEX
,
PARAMETER
::
img
=
(
0
.0
,
1.
0
)
! - local arrays -
INTEGER
::
lmp_start
(
atoms
%
ntype
)
REAL
::
alo
(
atoms
%
nlod
,
atoms
%
ntype
),
blo
(
atoms
%
nlod
,
atoms
%
ntype
),&
...
...
@@ -647,7 +647,7 @@
INTEGER
::
point
(
-1
:
1
)
REAL
::
rfac
COMPLEX
::
cfac
COMPLEX
,
PARAMETER
::
img
=
(
0
d0
,
1d
0
)
COMPLEX
,
PARAMETER
::
img
=
(
0
.0
,
1.
0
)
rfac
=
sqrt
(
tpi_const
/
3
)
DO
j
=
-1
,
1
...
...
@@ -781,10 +781,10 @@
DO
iband1
=
bandi1
,
bandf1
DO
iband2
=
bandi2
,
bandf2
rdum
=
eig_irr
(
iband2
,
nk
)
-
eig_irr
(
iband1
,
nk
)
IF
(
abs
(
rdum
)
.gt.
1
d-6
)
THEN
!1d
-6
IF
(
abs
(
rdum
)
.gt.
1
0.0
**
-6
)
THEN
!10.0**
-6
dcprod
(
iband2
,
iband1
,:)
=
dcprod
(
iband2
,
iband1
,:)
/
rdum
ELSE
dcprod
(
iband2
,
iband1
,:)
=
0
d
0
dcprod
(
iband2
,
iband1
,:)
=
0
.
0
END
IF
END
DO
END
DO
...
...
@@ -846,7 +846,7 @@
INTEGER
::
lm_0
,
lm_1
,
lm0
,
lm1
,
lm2
,
lm3
,
n0
,
nn
,
n
,
l1
,
l2
,
m1
,
m2
,
ikpt1
INTEGER
::
irecl_cmt
,
irecl_z
,
m
COMPLEX
::
cdum
COMPLEX
::
img
=
(
0
d0
,
1d
0
)
COMPLEX
::
img
=
(
0
.0
,
1.
0
)
! - local arrays -
INTEGER
::
gpt
(
3
,
lapw
%
nv
(
jsp
))
...
...
@@ -886,12 +886,12 @@
DO
l
=
0
,
atoms
%
lmaxd
DO
M
=
-
l
,
l
lm
=
lm
+
1
fcoeff
(
lm
,
-1
)
=
-
sqrt
(
1
d
0
*
(
l
+
M
+1
)
*
(
l
+
M
+2
)
/
(
2
*
(
2
*
l
+1
)
*
(
2
*
l
+3
))
)
fcoeff
(
lm
,
0
)
=
sqrt
(
1
d
0
*
(
l
-
M
+1
)
*
(
l
+
M
+1
)
/
((
2
*
l
+1
)
*
(
2
*
l
+3
))
)
fcoeff
(
lm
,
1
)
=
-
sqrt
(
1
d
0
*
(
l
-
M
+1
)
*
(
l
-
M
+2
)
/
(
2
*
(
2
*
l
+1
)
*
(
2
*
l
+3
))
)
gcoeff
(
lm
,
-1
)
=
sqrt
(
1
d
0
*
(
l
-
M
)
*
(
l
-
M
-1
)
/
(
2
*
(
2
*
l
-1
)
*
(
2
*
l
+1
))
)
gcoeff
(
lm
,
0
)
=
sqrt
(
1
d
0
*
(
l
-
M
)
*
(
l
+
M
)
/
((
2
*
l
-1
)
*
(
2
*
l
+1
))
)
gcoeff
(
lm
,
1
)
=
sqrt
(
1
d
0
*
(
l
+
M
)
*
(
l
+
M
-1
)
/
(
2
*
(
2
*
l
-1
)
*
(
2
*
l
+1
))
)
fcoeff
(
lm
,
-1
)
=
-
sqrt
(
1
.
0
*
(
l
+
M
+1
)
*
(
l
+
M
+2
)
/
(
2
*
(
2
*
l
+1
)
*
(
2
*
l
+3
))
)
fcoeff
(
lm
,
0
)
=
sqrt
(
1
.
0
*
(
l
-
M
+1
)
*
(
l
+
M
+1
)
/
((
2
*
l
+1
)
*
(
2
*
l
+3
))
)
fcoeff
(
lm
,
1
)
=
-
sqrt
(
1
.
0
*
(
l
-
M
+1
)
*
(
l
-
M
+2
)
/
(
2
*
(
2
*
l
+1
)
*
(
2
*
l
+3
))
)
gcoeff
(
lm
,
-1
)
=
sqrt
(
1
.
0
*
(
l
-
M
)
*
(
l
-
M
-1
)
/
(
2
*
(
2
*
l
-1
)
*
(
2
*
l
+1
))
)
gcoeff
(
lm
,
0
)
=
sqrt
(
1
.
0
*
(
l
-
M
)
*
(
l
+
M
)
/
((
2
*
l
-1
)
*
(
2
*
l
+1
))
)
gcoeff
(
lm
,
1
)
=
sqrt
(
1
.
0
*
(
l
+
M
)
*
(
l
+
M
-1
)
/
(
2
*
(
2
*
l
-1
)
*
(
2
*
l
+1
))
)
END
DO
END
DO
...
...
@@ -1014,11 +1014,11 @@
! Transform to cartesian coordinates
hlp
=
0
hlp
(
1
,
1
)
=
1
d0
/
sqrt
(
2d
0
)
hlp
(
1
,
3
)
=
-1
d0
/
sqrt
(
2d
0
)
hlp
(
2
,
1
)
=
-
img
/
sqrt
(
2
d
0
)
hlp
(
2
,
3
)
=
-
img
/
sqrt
(
2
d
0
)
hlp
(
3
,
2
)
=
1
d
0
hlp
(
1
,
1
)
=
1
.0
/
sqrt
(
2.
0
)
hlp
(
1
,
3
)
=
-1
.0
/
sqrt
(
2.
0
)
hlp
(
2
,
1
)
=
-
img
/
sqrt
(
2
.
0
)
hlp
(
2
,
3
)
=
-
img
/
sqrt
(
2
.
0
)
hlp
(
3
,
2
)
=
1
.
0
DO
iband1
=
bandi1
,
bandf1
DO
iband2
=
bandi2
,
bandf2
momentum
(
iband2
,
iband1
,:)
=
-
img
*
matmul
(
momentum
(
iband2
,
iband1
,:),
transpose
(
hlp
))
...
...
hybrid/mixedbasis.F90
View file @
567d369d
...
...
@@ -648,7 +648,7 @@ SUBROUTINE mixedbasis(atoms,kpts,DIMENSION,input,cell,sym,xcpot,hybrid,enpara,mp
rdum
=
rdum
+
2
*
rdum1
**
2
END
IF
IF
(
rdum1
.GT.
1
d
-6
)
THEN
IF
(
rdum1
.GT.
1
0.0
**
-6
)
THEN
IF
(
mpi
%
irank
==
0
)
THEN
WRITE
(
6
,
'(A)'
)
'mixedbasis: Bad orthonormality of '
&
//
lchar
(
l
)//
'-product basis. Increase tolerance.'
...
...
hybrid/olap.F90
View file @
567d369d
...
...
@@ -22,7 +22,7 @@
! - local -
INTEGER
::
i
,
j
,
itype
,
icent
,
ineq
REAL
::
g
,
r
,
fgr
COMPLEX
,
PARAMETER
::
img
=
(
0
d0
,
1d
0
)
COMPLEX
,
PARAMETER
::
img
=
(
0
.0
,
1.
0
)
INTEGER
::
dg
(
3
)
...
...
@@ -90,7 +90,7 @@
! - local -
INTEGER
::
i
,
j
,
k
,
itype
,
icent
,
ineq
REAL
::
g
,
r
,
fgr
COMPLEX
,
PARAMETER
::
img
=
(
0
d0
,
1d
0
)
COMPLEX
,
PARAMETER
::
img
=
(
0
.0
,
1.
0
)
INTEGER
::
dg
(
3
)
if
(
l_real
)
THEN
...
...
@@ -263,8 +263,8 @@ else
wfolap_inv
=
wfolap_inv
+
dot_product
(
cpw1
,
matmul
(
olappw
,
cpw2
))
! CALL dgemv('N',ngpt1,ngpt2,1
d0,olappw,ngpt1,real(cpw2),1,0d
0,rarr1,1)
! CALL dgemv('N',ngpt1,ngpt2,1
d0,olappw,ngpt1,aimag(cpw2),1,0d
0,rarr2,1)
! CALL dgemv('N',ngpt1,ngpt2,1
.0,olappw,ngpt1,real(cpw2),1,0.
0,rarr1,1)
! CALL dgemv('N',ngpt1,ngpt2,1
.0,olappw,ngpt1,aimag(cpw2),1,0.
0,rarr2,1)
!
! rdum1 = dotprod(cpw1,rarr1)
! rdum2 = dotprod(cpw1,rarr2)
...
...
@@ -319,8 +319,8 @@ else
wfolap_noinv
=
wfolap_noinv
+
dot_product
(
cpw1
,
matmul
(
olappw
,
cpw2
))
! CALL dgemv('N',ngpt1,ngpt2,1
d0,olappw,ngpt1,real(cpw2),1,0d
0,rarr1,1)
! CALL dgemv('N',ngpt1,ngpt2,1
d0,olappw,ngpt1,aimag(cpw2),1,0d
0,rarr2,1)
! CALL dgemv('N',ngpt1,ngpt2,1
.0,olappw,ngpt1,real(cpw2),1,0.
0,rarr1,1)
! CALL dgemv('N',ngpt1,ngpt2,1
.0,olappw,ngpt1,aimag(cpw2),1,0.
0,rarr2,1)
!
! rdum1 = dotprod(cpw1,rarr1)
! rdum2 = dotprod(cpw1,rarr2)
...
...
hybrid/spmvec.F90
View file @
567d369d
...
...
@@ -48,7 +48,7 @@
REAL
::
gnorm
COMPLEX
,
PARAMETER
::
img
=
(
0
d0
,
1d
0
)
COMPLEX
,
PARAMETER
::
img
=
(
0
.0
,
1.
0
)
! - local arrays -
REAL
::
vecinhlp
(
hybrid
%
nbasm
(
ikpt
))
...
...
@@ -150,11 +150,11 @@
indx0
=
sum
(
(/
(
((
2
*
l
+1
)
*
atoms
%
neq
(
itype
),
l
=
0
,
hybrid
%
lcutm1
(
itype
)),
itype
=
1
,
atoms
%
ntype
)
/)
)
+
hybrid
%
ngptm
indx1
=
sum
(
(/
(
((
2
*
l
+1
)
*
atoms
%
neq
(
itype
),
l
=
0
,
hybrid
%
lcutm1
(
itype
)),
itype
=
1
,
atoms
%
ntype
)
/)
)
CALL
DGEMV
(
'N'
,
indx1
,
indx0
,
1
d
0
,
coulomb_mtir
,(
hybrid
%
maxlcutm1
+1
)
**
2
*
atoms
,&
&
vecinhlp
(
ibasm
+1
:),
1
,
0
d
0
,
vecout
(
ibasm
+1
:),
1
)
CALL
DGEMV
(
'N'
,
indx1
,
indx0
,
1
.
0
,
coulomb_mtir
,(
hybrid
%
maxlcutm1
+1
)
**
2
*
atoms
,&
&
vecinhlp
(
ibasm
+1
:),
1
,
0
.
0
,
vecout
(
ibasm
+1
:),
1
)
CALL
DGEMV
(
'T'
,
indx1
,
hybrid
,
1
d
0
,
coulomb_mtir
(:
indx1
,
indx1
+1
:),&
&
indx1
,
vecinhlp
(
ibasm
+1
:),
1
,
0
d
0
,
vecout
(
ibasm
+
indx1
+1
:),
1
)
CALL
DGEMV
(
'T'
,
indx1
,
hybrid
,
1
.
0
,
coulomb_mtir
(:
indx1
,
indx1
+1
:),&
&
indx1
,
vecinhlp
(
ibasm
+1
:),
1
,
0
.
0
,
vecout
(
ibasm
+
indx1
+1
:),
1
)
! vecout(ibasm+1:ibasm+indx1) = matmul( coulomb_mtir(:indx1,:indx0),vecinhlp(ibasm+1:ibasm+indx0) )
! vecout(ibasm+indx1+1:ibasm+indx0) = matmul( conjg(transpose(coulomb_mtir(:indx1,indx1+1:indx0))),
...
...
@@ -295,7 +295,7 @@
REAL
::
gnorm
COMPLEX
,
PARAMETER
::
img
=
(
0
d0
,
1d
0
)
COMPLEX
,
PARAMETER
::
img
=
(
0
.0
,
1.
0
)
! - local arrays -
REAL
::
vecr
(
hybrid
%
maxindxm1
-1
),
veci
(
hybrid
%
maxindxm1
-1
)
...
...
@@ -399,12 +399,12 @@
indx0
=
sum
(
(/
(
((
2
*
l
+1
)
*
atoms
%
neq
(
itype
),
l
=
0
,
hybrid
%
lcutm1
(
itype
)),
itype
=
1
,
atoms
%
ntype
)
/)
)
+
hybrid
%
ngptm
indx1
=
sum
(
(/
(
((
2
*
l
+1
)
*
atoms
%
neq
(
itype
),
l
=
0
,
hybrid
%
lcutm1
(
itype
)),
itype
=
1
,
atoms
%
ntype
)
/)
)
CALL
ZGEMV
(
'N'
,
indx1
,
indx0
,(
1
d0
,
0d
0
),
coulomb_mtir
,&
CALL
ZGEMV
(
'N'
,
indx1
,
indx0
,(
1
.0
,
0.
0
),
coulomb_mtir
,&
&
(
hybrid
%
maxlcutm1
+1
)
**
2
*
atoms
,
vecinhlp
(
ibasm
+1
:),&
&
1
,(
0
d0
,
0d
0
),
vecout
(
ibasm
+1
:),
1
)
&
1
,(
0
.0
,
0.
0
),
vecout
(
ibasm
+1
:),
1
)
CALL
ZGEMV
(
'C'
,
indx1
,
hybrid
,(
1
d0
,
0d
0
),
coulomb_mtir
(:
indx1
,
indx1
+1
:)&
&
,
indx1
,
vecinhlp
(
ibasm
+1
:),
1
,(
0
d0
,
0d
0
),&
CALL
ZGEMV
(
'C'
,
indx1
,
hybrid
,(
1
.0
,
0.
0
),
coulomb_mtir
(:
indx1
,
indx1
+1
:)&
&
,
indx1
,
vecinhlp
(
ibasm
+1
:),
1
,(
0
.0
,
0.
0
),&
&
vecout
(
ibasm
+
indx1
+1
:),
1
)
! vecout(ibasm+1:ibasm+indx1) = matmul( coulomb_mtir(:indx1,:indx0),vecinhlp(ibasm+1:ibasm+indx0) )
...
...
@@ -426,7 +426,7 @@
indx1
=
sum
(
(/
(
((
2
*
l
+1
)
*
atoms
%
neq
(
itype
),
l
=
0
,
hybrid
%
lcutm1
(
itype
)),&
&
itype
=
1
,
atoms
%
ntype
)
/)
)
+
hybrid
%
ngptm
(
ikpt
)
call
zhpmv
(
'U'
,
indx1
,(
1
d0
,
0d0
),
coulomb_mtir
,
vecinhlp
(
ibasm
+1
),
1
,(
0d0
,
0d
0
),
vecout
(
ibasm
+1
),
1
)
call
zhpmv
(
'U'
,
indx1
,(
1
.0
,
0.0
),
coulomb_mtir
,
vecinhlp
(
ibasm
+1
),
1
,(
0.0
,
0.
0
),
vecout
(
ibasm
+1
),
1
)
#endif
iatom
=
0
...
...
hybrid/subvxc.F90
View file @
567d369d
...
...
@@ -382,7 +382,7 @@ CONTAINS
END
DO
END
DO
rc
=
CMPLX
(
0
d0
,
1d
0
)
**
(
l2
-
l1
)
! adjusts to a/b/ccof-scaling
rc
=
CMPLX
(
0
.0
,
1.
0
)
**
(
l2
-
l1
)
! adjusts to a/b/ccof-scaling
! ic counts the entry in vxc
ic
=
icentry
...
...
@@ -427,7 +427,7 @@ CONTAINS
END
DO
END
DO
rc
=
CMPLX
(
0
d0
,
1d
0
)
**
(
lp
-
l1
)
! adjusts to a/b/ccof-scaling
rc
=
CMPLX
(
0
.0
,
1.
0
)
**
(
lp
-
l1
)
! adjusts to a/b/ccof-scaling
IF
(
hmat
%
l_real
)
THEN
vxc
(
ic
)
=
vxc
(
ic
)
+
invsfct
*
REAL
(
rr
*
rc
*
bascof_lo
(
pp
,
mp
,
ikvecp
,
ilop
,
iatom
)
*
&
...
...
@@ -466,7 +466,7 @@ CONTAINS
END
DO
END
DO
rc
=
CMPLX
(
0
d0
,
1d
0
)
**
(
lp
-
l1
)
! adjusts to a/b/ccof-scaling
rc
=
CMPLX
(
0
.0
,
1.
0
)
**
(
lp
-
l1
)
! adjusts to a/b/ccof-scaling
IF
(
hmat
%
l_real
)
THEN
vxc
(
ic
)
=
vxc
(
ic
)
+
invsfct
*
REAL
(
rr
*
rc
*
bascof_lo
(
pp
,
mp
,
ikvecp
,
ilop
,
iatom
)
*
&
...
...
hybrid/symm_hf.F90
View file @
567d369d
...
...
@@ -122,7 +122,7 @@ SUBROUTINE symm_hf(kpts,nk,sym,dimension,hybdat,eig_irr,atoms,hybrid,cell,&
REAL
::
tolerance
,
pi
COMPLEX
::
cdum
COMPLEX
,
PARAMETER
::
img
=
(
0
d0
,
1d
0
)
COMPLEX
,
PARAMETER
::
img
=
(
0
.0
,
1.
0
)
! - local arrays -
INTEGER
::
neqvkpt
(
kpts
%
nkptf
)
...
...
hybrid/symmetrizeh.F90
View file @
567d369d
...
...
@@ -41,7 +41,7 @@ SUBROUTINE symmetrizeh(atoms,bk,DIMENSION,jsp,lapw,sym,kveclo,cell,nsymop,psym,h
INTEGER
::
igpt
,
igpt1
,
igpt2
,
igpt3
INTEGER
::
invsfct
,
invsfct1
,
idum
INTEGER
::
l
,
l1
,
lm
,
j
,
ok
,
ratom
,
ratom1
,
nrgpt
COMPLEX
,
PARAMETER
::
img
=
(
0
d0
,
1d
0
)
COMPLEX
,
PARAMETER
::
img
=
(
0
.0
,
1.
0
)
COMPLEX
::
cdum
,
cdum2
! local arrays
...
...
@@ -107,7 +107,7 @@ SUBROUTINE symmetrizeh(atoms,bk,DIMENSION,jsp,lapw,sym,kveclo,cell,nsymop,psym,h
rtaual
=
MATMUL
(
rot
(:,:,
isym
),
atoms
%
taual
(:,
iatom
))
+
trans
(:,
isym
)
iatom1
=
0
DO
ieq1
=
1
,
atoms
%
neq
(
itype
)
IF
(
ALL
(
ABS
(
MODULO
(
rtaual
-
atoms
%
taual
(:,
iiatom
+
ieq1
)
+1
d-12
,
1d0
))
.LT.
1d-10
))
THEN
!The 1d
-12 is a dirty fix.
IF
(
ALL
(
ABS
(
MODULO
(
rtaual
-
atoms
%
taual
(:,
iiatom
+
ieq1
)
+1
0.0
**
-12
,
1.0
))
.LT.
10.0
**
-10
))
THEN
!The 10.0**
-12 is a dirty fix.
iatom1
=
iiatom
+
ieq1
END
IF
END
DO
...
...
@@ -552,10 +552,10 @@ SUBROUTINE symmetrizeh(atoms,bk,DIMENSION,jsp,lapw,sym,kveclo,cell,nsymop,psym,h
REAL
::
stheta
,
ctheta
,
sphi
,
cphi
,
r
,
rvec1
(
3
)
INTEGER
::
l
,
lm
COMPLEX
::
c
COMPLEX
,
PARAMETER
::
img
=
(
0
d0
,
1d
0
)
COMPLEX
,
PARAMETER
::
img
=
(
0
.0
,
1.
0
)
call
timestart
(
"symm harmonics"
)
Y
(
1
)
=
0.282094791773878
d0
Y
(
1
)
=
0.282094791773878
IF
(
ll
.EQ.
0
)
RETURN
stheta
=
0
...
...
@@ -563,16 +563,16 @@ SUBROUTINE symmetrizeh(atoms,bk,DIMENSION,jsp,lapw,sym,kveclo,cell,nsymop,psym,h
sphi
=
0
cphi
=
0
r
=
SQRT
(
SUM
(
rvec
**
2
))
IF
(
r
.GT.
1
d
-16
)
THEN
IF
(
r
.GT.
1
0.0
**
-16
)
THEN
rvec1
=
rvec
/
r
ctheta
=
rvec1
(
3
)