Commit 9306e171 authored by Daniel Wortmann's avatar Daniel Wortmann

Modified sorting routine to use secondary key and used this in several places...

Modified sorting routine to use secondary key and used this in several places to simplify sorting of LAPW, stars...
parent d1782d3c
......@@ -52,7 +52,7 @@ inpgen/atom_sym.f inpgen/generator.f inpgen/read_record.f inpgen/soc_or_ssdw.f i
inpgen/bravais_symm.f inpgen/set_atom_core.f inpgen/spg_gen.f global/triang.f
inpgen/lapw_input.f inpgen/struct_input.f inpgen/write_struct.f
io/calculator.f global/ss_sym.f global/soc_sym.f math/inv3.f io/rw_symfile.f
global/sort.f kpoints/kptgen_hybrid.f kpoints/od_kptsgen.f kpoints/bravais.f kpoints/divi.f kpoints/brzone.f
kpoints/kptgen_hybrid.f kpoints/od_kptsgen.f kpoints/bravais.f kpoints/divi.f kpoints/brzone.f
kpoints/kptmop.f kpoints/kpttet.f init/bandstr1.F kpoints/ordstar.f kpoints/fulstar.f kpoints/kprep.f
kpoints/tetcon.f kpoints/kvecon.f init/boxdim.f global/radsra.f math/intgr.F global/differ.f math/inwint.f
math/outint.f xc-pot/gaunt.f math/grule.f
......@@ -60,7 +60,7 @@ kpoints/tetcon.f kpoints/kvecon.f init/boxdim.f global/radsra.f math/intgr.F glo
set(inpgen_F90 ${inpgen_F90} global/constants.f90 io/xsf_io.f90
eigen/orthoglo.F90 juDFT/usage_data.F90 math/ylm4.F90
global/chkmt.f90 inpgen/inpgen.f90 inpgen/set_inp.f90 inpgen/inpgen_help.f90 io/rw_inp.f90 juDFT/juDFT.F90 global/find_enpara.f90
global/sort.f90 global/chkmt.f90 inpgen/inpgen.f90 inpgen/set_inp.f90 inpgen/inpgen_help.f90 io/rw_inp.f90 juDFT/juDFT.F90 global/find_enpara.f90
inpgen/closure.f90 inpgen/inpgen_arguments.F90
juDFT/info.F90 juDFT/stop.F90 juDFT/args.F90 juDFT/time.F90 juDFT/init.F90 juDFT/sysinfo.F90 io/w_inpXML.f90 kpoints/julia.f90 global/utility.F90
init/compile_descr.F90 kpoints/kpoints.f90 io/xmlOutput.F90 kpoints/brzone2.f90 cdn/slab_dim.f90 cdn/slabgeom.f90 dos/nstm3.f90 cdn/int_21.f90
......
......@@ -61,7 +61,7 @@ CONTAINS
END DO
gvacl(n2) = SQRT(REAL(gvac(1)**2+gvac(2)**2))
ENDDO k_loop
CALL sort(n2,gvacl,gindex)
CALL sort(gindex,gvacl)
DO j = 1,n2
! gvac1d, gvac2d are now ordered by increasing length
gvac1d(j)=gvac1(gindex(j))
......
......@@ -55,7 +55,7 @@ CONTAINS
COMPLEX, ALLOCATABLE :: cph(:,:)
ALLOCATE(cph(MAXVAL(lapw%nv),2))
DO i=MIN(jintsp,iintsp),MAX(jintsp,iintsp)
CALL lapw%phase_factors(i,atoms%taual(:,na),noco%qss,cph(:,i))
CALL lapw%phase_factors(i,atoms%taual(:,na),-1*noco%qss,cph(:,i))
ENDDO
IF ((atoms%invsat(na).EQ.0) .OR. (atoms%invsat(na).EQ.1)) THEN
......
......@@ -179,7 +179,7 @@ CONTAINS
END DO
END DO
CALL sort(n,e,index)
CALL sort(index,e)
! Check if no deep eigenvalue is found
IF (e_min-MINVAL(e(1:n))>1.0) THEN
......
......@@ -7,7 +7,6 @@ global/radsra.f
global/radsrd.F
global/radsrdn.f
global/soc_sym.f
global/sort.f
global/ss_sym.f
global/starf.f
global/triang.f
......@@ -23,6 +22,7 @@ global/constants.f90
global/matrix_copy.F90
global/checkdop.F90
global/checkdopall.f90
global/sort.f90
global/genMTBasis.f90
global/chkmt.f90
global/convn.f90
......
......@@ -274,7 +274,7 @@
CALL juDFT_error("Too many atoms at same position.",calledby ="chkmt")
END IF
numNearestNeighbors(n) = MIN(maxCubeAtoms,iNeighborAtom)
CALL sort(iNeighborAtom,sqrDistances,distIndexList)
CALL sort(distIndexList(:iNeighborAtom),sqrDistances(:iNeighborAtom))
DO i = 1, numNearestNeighbors(n)
nearestNeighbors(i,n) = neighborAtoms(distIndexList(i))
nearestNeighborDists(i,n) = SQRT(sqrDistances(distIndexList(i)))
......@@ -293,7 +293,7 @@
! Sort distances and set MT radii for the atoms
CALL sort(atoms%ntype,nearestAtomDists,sortedDistList)
CALL sort(sortedDistList,nearestAtomDists)
rmt1 = -1.0
minRmts = -1.0
DO i = 1, atoms%ntype
......
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_sort
CONTAINS
SUBROUTINE sort(
> n,sk,
< js)
c********************************************************************
c heapsort routine
c input: n = number of objects
c sk = array of objects to be sorted
c output: js(i) = index of i'th smallest object m.w.
c modified to avoid numerical uncertainties for equal-lentgh
c vectors and to make sorting machine independent aug. 90, r.p.
c********************************************************************
IMPLICIT NONE
c
C .. Scalar Arguments ..
INTEGER, INTENT (IN) :: n
C ..
C .. Array Arguments ..
REAL, INTENT (IN) :: sk(n)
INTEGER, INTENT (OUT) :: js(n)
C ..
C .. Local Scalars ..
REAL eps,q
INTEGER i,ind,ir,j,l
C ..
C .. Data statements ..
DATA eps/1.e-10/
C ..
c
IF (n == 0) RETURN ! Nothing to do
IF (n == 1) THEN ! Not much to do
js(1) = 1
RETURN
END IF
DO i = 1,n
js(i) = i
ENDDO
c
l = n/2 + 1
ir = n
20 CONTINUE
IF (l.GT.1) THEN
l = l - 1
ind = js(l)
q = sk(ind)
ELSE
ind = js(ir)
q = sk(ind)
js(ir) = js(1)
ir = ir - 1
IF (ir.EQ.1) THEN
js(1) = ind
RETURN
END IF
END IF
i = l
j = l + l
30 IF (j.LE.ir) THEN
IF (j.LT.ir) THEN
c if(sk(js(j)).lt.sk(js(j+1))) j=j+1
IF ((sk(js(j+1))-sk(js(j))).GT.eps) j = j + 1
END IF
c if(q.lt.sk(js(j))) then
IF ((sk(js(j))-q).GT.eps) THEN
js(i) = js(j)
i = j
j = j + j
ELSE
j = ir + 1
END IF
GO TO 30
END IF
js(i) = ind
GO TO 20
END SUBROUTINE
END
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_sort
USE m_judft
CONTAINS
SUBROUTINE sort(ind,lv,lv1)
!********************************************************************
! heapsort routine
! input: lv = array of objects to be sorted
! lv1 = second array to use as secondary sort key
! output: ind(i) = index of i'th smallest object
!********************************************************************
IMPLICIT NONE
!
REAL, INTENT (IN) :: lv(:)
INTEGER, INTENT (OUT) :: ind(:)
REAL,INTENT(IN),OPTIONAL :: lv1(:)
! ..
! .. Local Scalars ..
REAL eps,q,q1
INTEGER i,idx,ir,j,l,n
REAL,ALLOCATABLE :: llv(:)
! ..
! .. Data statements ..
DATA eps/1.e-10/
! ..
!
n=SIZE(ind)
IF (n>SIZE(lv)) CALL judft_error("BUG: incosistent dimensions")
ALLOCATE(llv(n))
IF (PRESENT(lv1)) THEN
IF (n>SIZE(lv1)) CALL judft_error("BUG: incosistent dimensions")
llv=lv1
ELSE
llv=(/(1.*i,i=1,n)/)
END IF
IF (n == 0) RETURN ! Nothing to do
IF (n == 1) THEN ! Not much to do
ind(1) = 1
RETURN
END IF
DO i = 1,n
ind(i) = i
ENDDO
!
l = n/2 + 1
ir = n
DO
IF (l.GT.1) THEN
l = l - 1
idx = ind(l)
q = lv(idx)
q1= llv(idx)
ELSE
idx = ind(ir)
q = lv(idx)
q1= llv(idx)
ind(ir) = ind(1)
ir = ir - 1
IF (ir.EQ.1) THEN
ind(1) = idx
RETURN
END IF
END IF
i = l
j = l + l
DO WHILE(j.LE.ir)
IF (j.LT.ir) THEN
! if(lv(ind(j)).lt.lv(ind(j+1))) j=j+1
IF (((lv(ind(j+1))-lv(ind(j))).GE.eps).OR. &!Standard comparison
((ABS((lv(ind(j+1))-lv(ind(j))))<eps).AND.&!Same length, check second key
((llv(ind(j+1))-llv(ind(j))).GE.eps))) &
j=j+1
END IF
! if(q.lt.lv(ind(j))) then
IF ((lv(ind(j))-q).GE.eps.OR.&
(ABS((lv(ind(j))-q))<eps.AND.(llv(ind(j))-q1).GE.eps))THEN
ind(i) = ind(j)
i = j
j = j + j
ELSE
j = ir + 1
END IF
enddo
ind(i) = idx
ENDDO
END SUBROUTINE sort
END MODULE m_sort
......@@ -148,8 +148,10 @@ CONTAINS
ENDIF
!---> sort for increasing length sk2
CALL sort(stars%ng2,stars%sk2,index)
DO k = 1,stars%ng2
gsk3(k) = (stars%mx1+stars%kv2(1,k)) + (stars%mx2+stars%kv2(2,k))*(2*stars%mx1+1)
ENDDO
CALL sort(INDEX(:stars%ng2),stars%sk2(:stars%ng2),gsk3(:stars%ng2))
DO k = 1,stars%ng2
kv3rev(k,1) = stars%kv2(1,INDEX(k))
kv3rev(k,2) = stars%kv2(2,INDEX(k))
......@@ -163,33 +165,6 @@ CONTAINS
stars%phi2(k) = phi3(k)
ENDDO
!--roa ....
!......sort stars of equal length 2D .....
i=1
gla=0.
gsk3(1)=0.0
eps=1.e-10
DO k = 2,stars%ng2
IF (stars%sk2(k)-gla.GE.eps) i=i+1
gla=stars%sk2(k)
gmi = (stars%mx1+stars%kv2(1,k)) + (stars%mx2+stars%kv2(2,k))*(2*stars%mx1+1)
gsk3(k) = i * ((2*stars%mx1+1)*(2*stars%mx2+1)+9)+gmi
ENDDO
CALL sort(stars%ng2,gsk3,index2)
DO k = 1,stars%ng2
kv3rev(k,1) = stars%kv2(1,index2(k))
kv3rev(k,2) = stars%kv2(2,index2(k))
gsk3(k) = stars%sk2(index2(k))
phi3(k) = stars%phi2(index2(k))
ENDDO
DO k = 1,stars%ng2
stars%kv2(1,k) = kv3rev(k,1)
stars%kv2(2,k) = kv3rev(k,2)
stars%sk2(k) = gsk3(k)
stars%phi2(k) = phi3(k)
! if (index2(k).ne.k) write(*,*) ' ic2: ',k,index2(k)
ENDDO
!--roa ....
WRITE (6,'(/'' nq2='',i4/'' k,kv2(1,2), sk2, phi2''&
& /(3i4,f10.5,f10.5))')&
......@@ -252,8 +227,13 @@ CONTAINS
ENDDO
!---> sort for increasing length sk3
CALL sort(stars%ng3,stars%sk3,index)
! secondary key for equal length stars
DO k = 1,stars%ng3
gsk3(k) = (stars%mx1+stars%kv3(1,k)) +&
& (stars%mx2+stars%kv3(2,k))*(2*stars%mx1+1) +&
& (stars%mx3+stars%kv3(3,k))*(2*stars%mx1+1)*(2*stars%mx2+1)
ENDDO
CALL sort(index,stars%sk3,gsk3)
ALLOCATE (ig2p(stars%ng3))
......@@ -271,41 +251,7 @@ CONTAINS
stars%sk3(k) = gsk3(k)
stars%ig2(k) = ig2p(k)
ENDDO
!
!--roa ....
!......sort stars of equal length 3D .....
i=1
gla=0.
gsk3(1)=0.
eps=1.e-10
DO k = 2,stars%ng3
IF (stars%sk3(k)-gla.GE.eps) i=i+1
gla = stars%sk3(k)
gmi = (stars%mx1+stars%kv3(1,k)) +&
& (stars%mx2+stars%kv3(2,k))*(2*stars%mx1+1) +&
& (stars%mx3+stars%kv3(3,k))*(2*stars%mx1+1)*(2*stars%mx2+1)
gsk3(k) = i * (9.+(2*stars%mx1+1)*(2*stars%mx2+1)*(2*stars%mx3+1)) + gmi
ENDDO
CALL sort(stars%ng3,gsk3,index3)
DO k = 1,stars%ng3
kv3rev(k,1) = stars%kv3(1,index3(k))
kv3rev(k,2) = stars%kv3(2,index3(k))
kv3rev(k,3) = stars%kv3(3,index3(k))
gsk3(k) = stars%sk3(index3(k))
ig2p(k) = stars%ig2(index3(k))
ENDDO
DO k = 1,stars%ng3
stars%kv3(1,k) = kv3rev(k,1)
stars%kv3(2,k) = kv3rev(k,2)
stars%kv3(3,k) = kv3rev(k,3)
stars%sk3(k) = gsk3(k)
stars%ig2(k) = ig2p(k)
! if (index3(k).ne.k) write(*,*) ' ic: ',k,index3(k)
ENDDO
DEALLOCATE (ig2p)
!--roa ....
!
!---> determine true gmax and change old gmax to new gmax
!
......@@ -682,8 +628,13 @@ CONTAINS
ENDIF
!---> sort for increasing length sk3
!
CALL sort(stars%ng3,stars%sk3,index)
! secondary key for equal length stars
DO k = 1,stars%ng3
gsk3(k) = (stars%mx1+stars%kv3(1,k)) +&
& (stars%mx2+stars%kv3(2,k))*(2*stars%mx1+1) +&
& (stars%mx3+stars%kv3(3,k))*(2*stars%mx1+1)*(2*stars%mx2+1)
ENDDO
CALL sort(index,stars%sk3,gsk3)
DO k = 1,stars%ng3
kv3rev(k,1) = stars%kv3(1,INDEX(k))
kv3rev(k,2) = stars%kv3(2,INDEX(k))
......@@ -697,36 +648,6 @@ CONTAINS
stars%sk3(k) = gsk3(k)
ENDDO
!
!--roa ....
!......sort stars of equal length 3D .....
i=1
gla=0.
gsk3(1)=0.
eps=1.e-10
DO k = 2,stars%ng3
IF (stars%sk3(k)-gla.GE.eps) i=i+1
gla = stars%sk3(k)
gmi = (stars%mx1+stars%kv3(1,k)) +&
& (stars%mx2+stars%kv3(2,k))*(2*stars%mx1+1) +&
& (stars%mx3+stars%kv3(3,k))*(2*stars%mx1+1)*(2*stars%mx2+1)
gsk3(k) = i * (9.+(2*stars%mx1+1)*(2*stars%mx2+1)*(2*stars%mx3+1)) + gmi
ENDDO
CALL sort(stars%ng3,gsk3,index3)
DO k = 1,stars%ng3
kv3rev(k,1) = stars%kv3(1,index3(k))
kv3rev(k,2) = stars%kv3(2,index3(k))
kv3rev(k,3) = stars%kv3(3,index3(k))
gsk3(k) = stars%sk3(index3(k))
ENDDO
DO k = 1,stars%ng3
stars%kv3(1,k) = kv3rev(k,1)
stars%kv3(2,k) = kv3rev(k,2)
stars%kv3(3,k) = kv3rev(k,3)
stars%sk3(k) = gsk3(k)
! if (index3(k).ne.k) write(*,*) ' ic: ',k,index3(k)
ENDDO
!--roa ....
!
!---> determine true gmax and change old gmax to new gmax
!
WRITE (6,8060) stars%gmax, stars%sk3(stars%ng3)
......
......@@ -1481,7 +1481,7 @@ input%preconditioning_param = evaluateFirstOnly(xmlGetAttributeValue('/fleurInpu
DO iLLO = 1, speciesNLO(iSpecies)
speciesLLOReal(iLLO) = speciesLLO(iLLO)
END DO
CALL sort(speciesNLO(iSpecies),speciesLLOReal,loOrderList)
CALL sort(loOrderList(:speciesNLO(iSpecies)),speciesLLOReal(:speciesNLO(iSpecies)))
DEALLOCATE(speciesLLOReal)
! apply species parameters to atom groups
......
......@@ -150,7 +150,7 @@ CONTAINS
! ..
! .. Local Arrays ..
REAL :: s(3),sq(3)
REAL,ALLOCATABLE :: rk(:),rkq(:)
REAL,ALLOCATABLE :: rk(:),rkq(:),rkqq(:)
INTEGER,ALLOCATABLE :: gvec(:,:),index3(:)
......@@ -169,7 +169,7 @@ CONTAINS
CALL lapw%alloc(cell,input,noco)
ALLOCATE(gvec(3,SIZE(lapw%gvec,2)))
ALLOCATE(rk(SIZE(lapw%gvec,2)),rkq(SIZE(lapw%gvec,2)))
ALLOCATE(rk(SIZE(lapw%gvec,2)),rkq(SIZE(lapw%gvec,2)),rkqq(SIZE(lapw%gvec,2)))
ALLOCATE(index3(SIZE(lapw%gvec,2)))
......@@ -210,45 +210,19 @@ CONTAINS
ENDDO
lapw%nv(ispin) = n
!Sort according to k+g
! CALL sort(lapw%nv(ispin),rkq,index3)
! DO n=1,lapw%nv(ispin)
! lapw%gvec(:,n,ispin) = gvec(:,index3(n))
! lapw%rk(n,ispin) = rk(index3(n))
! ENDDO
!
!---> sort by shell-metzner
!
! (for spin-spirals & LO's we have to sort according to the k+G's (rkq), not
! the k+G+q's (rk). Otherwise we might couple an LO to k+G1+q and k+G2-q !)
! gb01
lapw%gvec(:,:,ispin)=gvec
lapw%rk(:,ispin)=rk
m = lapw%nv(ispin)
80 m = m/2
IF (m.LE.0) GO TO 130
k = n - m
j = 1
90 i = j
100 l = i + m
IF (rkq(i).GT.rkq(l)) GO TO 120
110 j = j + 1
IF (j.GT.k) GO TO 80
GO TO 90
120 t = rkq(i)
rkq(i) = rkq(l)
rkq(l) = t
t = lapw%rk(i,ispin)
lapw%rk(i,ispin) = lapw%rk(l,ispin)
lapw%rk(l,ispin) = t
itt = lapw%gvec(:,i,ispin)
lapw%gvec(:,i,ispin) = lapw%gvec(:,l,ispin)
lapw%gvec(:,l,ispin) = itt
i = i - m
IF (i.LT.1) GO TO 110
GO TO 100
130 CONTINUE
!+gu
!Sort according to k+g, first construct secondary sort key
DO k = 1,lapw%nv(ispin)
rkqq(k) = (mk1+gvec(1,k)) + (mk2+gvec(2,k))*(2*mk1+1) +&
(mk3+gvec(3,k))*(2*mk1+1)*(2*mk2+1)
ENDDO
CALL sort(index3(:lapw%nv(ispin)),rkq,rkqq)
PRINT *,"L:",SIZE(index3),lapw%nv(ispin)
PRINT *,"L:",MINVAL(index3),MAXVAL(index3),SUM(index3)
DO n=1,lapw%nv(ispin)
lapw%gvec(:,n,ispin) = gvec(:,index3(n))
lapw%rk(n,ispin) = rk(index3(n))
PRINT *,"RK:",n,lapw%rk(n,ispin)
ENDDO
!---> determine pairs of K-vectors, where K_z = K'_-z to use
!---> z-reflection
IF (l_zref) THEN
......@@ -266,75 +240,54 @@ CONTAINS
nred=n
IF (PRESENT(mpi)) THEN
IF (mpi%n_size.GT.1) THEN
!
!---> order K's in sequence K_1,...K_n | K_0,... | K_-1,....K_-n
!
n_inner = lapw%nv(ispin) - nred
IF (MOD(nred,mpi%n_size).EQ.0) THEN
n_bound = nred
ELSE
n_bound = (1+INT( nred/mpi%n_size ))*mpi%n_size
ENDIF
IF (lapw%nv(ispin) - nred + n_bound.GT.SIZE(lapw%gvec,2)) THEN
CALL juDFT_error("BUG:z-ref & ev || : dimension too small!" ,calledby ="types_lapw")
ENDIF
i = 1
j = 1
DO n = 1, nred
IF (lapw%matind(n,1).EQ.lapw%matind(n,2)) THEN
index3(lapw%matind(n,1)) = n_inner + i
i = i + 1
!
!---> order K's in sequence K_1,...K_n | K_0,... | K_-1,....K_-n
!
n_inner = lapw%nv(ispin) - nred
IF (MOD(nred,mpi%n_size).EQ.0) THEN
n_bound = nred
ELSE
index3(lapw%matind(n,1)) = j
index3(lapw%matind(n,2)) = j + n_bound
j = j + 1
n_bound = (1+INT( nred/mpi%n_size ))*mpi%n_size
ENDIF
IF (lapw%nv(ispin) - nred + n_bound.GT.SIZE(lapw%gvec,2)) THEN
CALL juDFT_error("BUG:z-ref & ev || : dimension too small!" ,calledby ="types_lapw")
ENDIF
ENDDO
!---> resort the rk,k1,k2,k3 and lapw%matind arrays:
DO n = 1, lapw%nv(ispin)
rk(n) = lapw%rk(n,ispin)
gvec(:,n) = lapw%gvec(:,n,ispin)
ENDDO
DO n = lapw%nv(ispin), 1, -1
lapw%rk(index3(n),ispin) = rk(n)
lapw%gvec(:,index3(n),ispin) = gvec(:,n)
ENDDO
DO n = nred + 1, n_bound
lapw%rk(n,ispin) = lapw%rk(lapw%nv(ispin),ispin)
lapw%gvec(:,n,ispin) = lapw%gvec(:,lapw%nv(ispin),ispin)
ENDDO
lapw%nv(ispin) = lapw%nv(ispin) - nred + n_bound
ENDIF
ENDIF
ENDIF
IF (noco%l_ss) THEN ! sort additionally like in strgn1... gb
i = 1
gla = 0.
rk(1) = 0.0
eps=1.e-10
DO k = 1,lapw%nv(ispin)
IF (rkq(k)-gla.GE.eps) i=i+1
gla = rkq(k)
gmi = (mk1+lapw%gvec(1,k,ispin)) + (mk2+lapw%gvec(2,k,ispin))*(2*mk1+1) +&
(mk3+lapw%gvec(3,k,ispin))*(2*mk1+1)*(2*mk2+1)
rk(k) = i * (9.+(2*mk1+1)*(2*mk2+1)*(2*mk3+1)) + gmi
ENDDO
CALL sort(lapw%nv(ispin),rk,index3)
DO k = 1,lapw%nv(ispin)
gvec(:,k) = lapw%gvec(:,index3(k),ispin)
rk(k) = lapw%rk(index3(k),ispin)
ENDDO
lapw%gvec(:,:lapw%nv(ispin),ispin) = gvec(:,:lapw%nv(ispin))
lapw%rk(:lapw%nv(ispin),ispin) = rk(:lapw%nv(ispin))
i = 1
j = 1
DO n = 1, nred
IF (lapw%matind(n,1).EQ.lapw%matind(n,2)) THEN
index3(lapw%matind(n,1)) = n_inner + i
i = i + 1
ELSE
index3(lapw%matind(n,1)) = j
index3(lapw%matind(n,2)) = j + n_bound
j = j + 1
ENDIF
ENDDO
!---> resort the rk,k1,k2,k3 and lapw%matind arrays:
DO n = 1, lapw%nv(ispin)
rk(n) = lapw%rk(n,ispin)
gvec(:,n) = lapw%gvec(:,n,ispin)
ENDDO
DO n = lapw%nv(ispin), 1, -1
lapw%rk(index3(n),ispin) = rk(n)
lapw%gvec(:,index3(n),ispin) = gvec(:,n)
ENDDO
DO n = nred + 1, n_bound
lapw%rk(n,ispin) = lapw%rk(lapw%nv(ispin),ispin)
lapw%gvec(:,n,ispin) = lapw%gvec(:,lapw%nv(ispin),ispin)
ENDDO
lapw%nv(ispin) = lapw%nv(ispin) - nred + n_bound
ENDIF
ENDIF
ENDIF
!-gu
DO k=1,lapw%nv(ispin)
lapw%vk(:,k,ispin)=lapw%bkpt+lapw%gvec(:,k,ispin)+(ispin-1.5)*noco%qss
lapw%gk(:,k,ispin)=MATMUL(TRANSPOSE(cell%bmat),lapw%vk(:,k,ispin))/MAX (lapw%rk(k,ispin),1.0e-30)
ENDDO
IF (.NOT.noco%l_ss.AND.input%jspins==2) THEN
!Second spin is the same
lapw%nv(2)=lapw%nv(1)
......@@ -372,41 +325,41 @@ CONTAINS
CONTAINS
SUBROUTINE priv_lo_basis_setup(lapw,atoms,sym,noco,cell)
USE m_types_setup
IMPLICIT NONE
TYPE(t_lapw),INTENT(INOUT):: lapw
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_noco),INTENT(IN) :: noco
SUBROUTINE priv_lo_basis_setup(lapw,atoms,sym,noco,cell)
USE m_types_setup
IMPLICIT NONE
TYPE(t_lapw),INTENT(INOUT):: lapw
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_noco),INTENT(IN) :: noco
INTEGER:: n,na,nn,np,lo,nkvec_sv,nkvec(atoms%nlod,2),iindex
IF (.NOT.ALLOCATED(lapw%kvec)) THEN
ALLOCATE(lapw%kvec(2*(2*atoms%llod+1),atoms%nlod,atoms%nat))
ALLOCATE(lapw%nkvec(atoms%nlod,atoms%nat))
ALLOCATE(lapw%index_lo(atoms%nlod,atoms%nat))
ENDIF
iindex=0
na=0
nkvec_sv=0
DO n=1,atoms%ntype
DO nn=1,atoms%neq(n)
na=na+1
if (atoms%invsat(na)>1) cycle
!np = MERGE(oneD%ods%ngopr(na),sym%invtab(atoms%ngopr(na)),oneD%odi%d1)
np=sym%invtab(atoms%ngopr(na))
CALL priv_vec_for_lo(atoms,sym,na,n,np,noco,lapw,cell)
DO lo = 1,atoms%nlo(n)
lapw%index_lo(lo,na)=iindex
iindex=iindex+lapw%nkvec(lo,na)
ENDDO
ENDDO
ENDDO
END SUBROUTINE priv_lo_basis_setup
INTEGER:: n,na,nn,np,lo,nkvec_sv,nkvec(atoms%nlod,2),iindex
IF (.NOT.ALLOCATED(lapw%kvec)) THEN
ALLOCATE(lapw%kvec(2*(2*atoms%llod+1),atoms%nlod,atoms%nat))
ALLOCATE(lapw%nkvec(atoms%nlod,atoms%nat))
ALLOCATE(lapw%index_lo(atoms%nlod,atoms%nat))
ENDIF
iindex=0
na=0
nkvec_sv=0
DO n=1,atoms%ntype
DO nn=1,atoms%neq(n)
na=na+1
if (atoms%invsat(na)>1) cycle
!np = MERGE(oneD%ods%ngopr(na),sym%invtab(atoms%ngopr(na)),oneD%odi%d1)
np=sym%invtab(atoms%ngopr(na))
CALL priv_vec_for_lo(atoms,sym,na,n,np,noco,lapw,cell)
DO lo = 1,atoms%nlo(n)
lapw%index_lo(lo,na)=iindex
iindex=iindex+lapw%nkvec(lo,na)
ENDDO
ENDDO
ENDDO
END SUBROUTINE priv_lo_basis_setup