Commit c6f8513c authored by Daniel Wortmann's avatar Daniel Wortmann

Started some redesign of hsmt in spherical part. Generated new subroutine to...

Started some redesign of hsmt in spherical part. Generated new subroutine to calculate the overlap matrix
parent 5d1544b7
......@@ -7,6 +7,8 @@ eigen/hlomat.F90
eigen/hsint.F90
eigen/hsmt.F90
eigen/hsmt_extra.F90
eigen/hsmt_fjgj.F90
eigen/hsmt_overlap.F90
eigen/hsmt_hlptomat.F90
eigen/hsmt_nonsph.F90
eigen/hsmt_soc_init.F90
......
......@@ -69,6 +69,8 @@ CONTAINS
USE m_hsmt_sph
USE m_hsmt_extra
USE m_types
USE m_hsmt_fjgj
USE m_hsmt_overlap
IMPLICIT NONE
TYPE(t_mpi),INTENT(IN) :: mpi
TYPE(t_dimension),INTENT(IN) :: DIMENSION
......@@ -119,6 +121,11 @@ CONTAINS
TYPE(t_rsoc):: rsoc
#if 1==2
REAL, ALLOCATABLE :: fj_test(:,:,:,:),gj_test(:,:,:,:)
COMPLEX,ALLOCATABLE :: bb_test(:)
#endif
!
CALL timestart("hsmt init")
l_socfirst= noco%l_soc .AND. noco%l_noco .AND. (.NOT. noco%l_ss)
......@@ -182,11 +189,36 @@ CONTAINS
CALL timestop("hsmt init")
DO isp = jsp_start,jsp_end
#if 1==2
ALLOCATE(bb_test(size(bb)))
bb_test=0.0
bb_test=bb
fj=0.0
gj=0.0
#endif
CALL timestart("hsmt spherical")
CALL hsmt_sph(DIMENSION,atoms,SUB_COMM,n_size,n_rank,sphhar,isp,ab_dim,&
input,hlpmsize,noco,l_socfirst,cell,nintsp,lapw,enpara%el0,usdus,&
vr,gk,rsoc,isigma, aa,bb,fj,gj)
CALL timestop("hsmt spherical")
#if 1==2
ALLOCATE(fj_test(size(fj,1),0:size(fj,2)-1,size(fj,3),size(fj,4)))
ALLOCATE(gj_test(size(fj,1),0:size(fj,2)-1,size(fj,3),size(fj,4)))
fj_test=0.0
gj_test=0.0
CALL timestart("hsmt_fjgj")
CALL hsmt_fjgj(atoms,isp,noco,l_socfirst,cell,nintsp, lapw,usdus,fj_test,gj_test)
CALL timestop("hsmt_fjgj")
print *,"fj diff",maxval(abs(fj-fj_test))
print *,"gj diff",maxval(abs(gj-gj_test))
deallocate(fj_test,gj_test)
CALL timestart("hsmt_overlap")
CALL hsmt_overlap(input,atoms,n_size,n_rank,isp,l_socfirst,hlpmsize,ab_dim,&
noco,cell,nintsp, lapw,usdus,gk,fj,gj,bb_test)
CALL timestop("hsmt_overlap")
print *,"bb_diff",maxval(abs(bb-bb_test))
deallocate(bb_test)
#endif
IF (.NOT.input%secvar) THEN
CALL timestart("hsmt extra")
IF (ANY(atoms%nlo>0).OR.ANY(atoms%lda_u%l.GE.0)) &
......@@ -202,7 +234,10 @@ CONTAINS
CALL timestop("hsmt non-spherical")
ENDIF
ENDDO
#if 1==2
print *,"hsmt checkmode"
call judft_end("test",0)
#endif
RETURN
END SUBROUTINE hsmt
END MODULE m_hsmt
!--------------------------------------------------------------------------------
! 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_hsmt_fjgj
use m_juDFT
implicit none
CONTAINS
SUBROUTINE hsmt_fjgj(atoms,isp,noco,l_socfirst,cell,nintsp, lapw,usdus,fj,gj)
!Calculate the fj&gj array which contain the part of the A,B matching coeff. depending on the
!radial functions at the MT boundary as contained in usdus
USE m_constants, ONLY : fpi_const
USE m_sphbes
USE m_dsphbs
USE m_types
IMPLICIT NONE
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_lapw),INTENT(IN) :: lapw!lpaw%nv_tot is updated
TYPE(t_usdus),INTENT(IN) :: usdus
! ..
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: isp
INTEGER, INTENT (IN) :: nintsp
LOGICAL,INTENT(IN) :: l_socfirst
REAL,INTENT(OUT) :: fj(:,0:,:,:),gj(:,0:,:,:)
! ..
! .. Local Scalars ..
REAL con1,ff,gg,gs,ws
INTEGER k,l,lo,n,iintsp
! .. Local Arrays ..
REAL gb(0:atoms%lmaxd), fb(0:atoms%lmaxd)
LOGICAL apw(0:atoms%lmaxd)
! ..
con1 = fpi_const/sqrt(cell%omtil)
DO iintsp = 1,nintsp
!$OMP parallel do DEFAULT(SHARED)&
!$OMP PRIVATE(n,l,apw,lo,k,gs,fb,gb,ws,ff,gg)
DO n = 1,atoms%ntype
DO l = 0,atoms%lmax(n)
apw(l)=any(atoms%l_dulo(:atoms%nlo(n),n))
#ifdef CPP_APW
IF (atoms%lapw_l(n).GE.l) apw(l) = .false.
#endif
ENDDO
DO lo = 1,atoms%nlo(n)
IF (atoms%l_dulo(lo,n)) apw(atoms%llo(lo,n)) = .true.
ENDDO
DO k = 1,lapw%nv(iintsp)
gs = lapw%rk(k,iintsp)*atoms%rmt(n)
CALL sphbes(atoms%lmax(n),gs, fb)
CALL dsphbs(atoms%lmax(n),gs,fb, gb)
DO l = 0,atoms%lmax(n)
!---> set up wronskians for the matching conditions for each ntype
ws = con1/(usdus%uds(l,n,isp)*usdus%dus(l,n,isp)&
- usdus%us(l,n,isp)*usdus%duds(l,n,isp))
ff = fb(l)
gg = lapw%rk(k,iintsp)*gb(l)
IF ( apw(l) ) THEN
fj(k,l,n,iintsp) = 1.0*con1 * ff / usdus%us(l,n,isp)
gj(k,l,n,iintsp) = 0.0d0
ELSE
IF (noco%l_constr.or.l_socfirst) THEN
!---> in a constrained calculation fj and gj are needed
!---> both local spin directions (isp) at the same time
fj(k,l,n,isp) = ws * ( usdus%uds(l,n,isp)*gg - usdus%duds(l,n,isp)*ff )
gj(k,l,n,isp) = ws * ( usdus%dus(l,n,isp)*ff - usdus%us(l,n,isp)*gg )
ELSE
!---> in a spin-spiral calculation fj and gj are needed
!---> both interstitial spin directions at the same time
!---> In all other cases iintsp runs from 1 to 1.
fj(k,l,n,iintsp) = ws * ( usdus%uds(l,n,isp)*gg - usdus%duds(l,n,isp)*ff )
gj(k,l,n,iintsp) = ws * ( usdus%dus(l,n,isp)*ff - usdus%us(l,n,isp)*gg )
ENDIF
ENDIF
ENDDO
ENDDO ! k = 1, lapw%nv
ENDDO ! n = 1,atoms%ntype
!$OMP end parallel do
ENDDO ! iintsp = 1,nintsp
RETURN
END SUBROUTINE hsmt_fjgj
END MODULE m_hsmt_fjgj
!--------------------------------------------------------------------------------
! 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_hsmt_overlap
use m_juDFT
implicit none
CONTAINS
SUBROUTINE hsmt_overlap(input,atoms,n_size,n_rank,isp,l_socfirst,hlpmsize,ab_dim,&
noco,cell,nintsp, lapw,usdus,gk,fj,gj,bb)
!Calculate overlap matrix
#include"cpp_double.h"
USE m_constants, ONLY : fpi_const,tpi_const
USE m_types
USE m_hsmt_spinor
IMPLICIT NONE
TYPE(t_input),INTENT(IN) :: input
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_lapw),INTENT(INOUT) :: lapw!lpaw%nv_tot is updated
TYPE(t_usdus),INTENT(INOUT) :: usdus
! ..
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: isp
INTEGER, INTENT (IN) :: n_size,n_rank,ab_dim
INTEGER, INTENT (IN) :: hlpmsize,nintsp
LOGICAL, INTENT (IN) :: l_socfirst
! ..
! .. Array Arguments ..
REAL,INTENT(IN) :: gk(:,:,:)
REAL,INTENT(IN) :: fj(:,0:,:,:),gj(:,0:,:,:)
#ifdef CPP_INVERSION
REAL, INTENT (INOUT) :: bb(:)!(matsize)
#else
COMPLEX, INTENT (INOUT) :: bb(:)
#endif
! ..
! .. Local Scalars ..
REAL tnn(3),fct,ski(3),fjkiln,gjkiln
COMPLEX chi11,chi21,chi22
INTEGER ii,iii,ij,k,ki,kj,l,n,n0,n1,nn,kjmax, nsp,iintsp,jintsp
INTEGER nc ,kii
! ..
! .. Local Arrays ..
COMPLEX chi(2,2),bbwa(maxval(lapw%nv))
REAL qssbti(3),qssbtj(3)
REAL fleg1(0:atoms%lmaxd),fleg2(0:atoms%lmaxd),fl2p1(0:atoms%lmaxd)
REAL, ALLOCATABLE :: plegend(:,:)
REAL, ALLOCATABLE :: rph(:,:),cph(:,:)
COMPLEX, ALLOCATABLE :: bbhlp(:)
IF ( noco%l_noco .AND. (.NOT. noco%l_ss) ) THEN
!---> pk non-collinear
!---> initialize help array
ALLOCATE ( bbhlp(hlpmsize) )
bbhlp=cmplx(0.0,0.0)
ENDIF
! ..
DO l = 0,atoms%lmaxd
fleg1(l) = real(l+l+1)/real(l+1)
fleg2(l) = real(l)/real(l+1)
fl2p1(l) = real(l+l+1)/fpi_const
END DO
!
!$OMP PARALLEL DEFAULT(shared)&
!$OMP PRIVATE(kii,ki,nc,iii,kjmax,ski,kj,plegend,l,n1,n)&
!$OMP PRIVATE(n0,rph,cph,nn,tnn,fjkiln,gjkiln)&
!$OMP PRIVATE(fct,ij)&
!$OMP PRIVATE(chi,chi11,chi21,chi22,nsp)&
!$OMP PRIVATE(bbwa,ii)
ALLOCATE(rph(maxval(lapw%nv),ab_dim))
ALLOCATE(cph(maxval(lapw%nv),ab_dim))
ALLOCATE(plegend(maxval(lapw%nv),0:atoms%lmaxd))
plegend=0.0
DO iintsp = 1,nintsp
IF (iintsp.EQ.1) THEN
qssbti = - noco%qss/2
ELSE
qssbti = + noco%qss/2
ENDIF
DO jintsp = 1,iintsp
IF (jintsp.EQ.1) THEN
qssbtj = - noco%qss/2
ELSE
qssbtj = + noco%qss/2
ENDIF
nc = 0 ! maybe IF (iintsp.EQ
IF ( noco%l_noco .AND. (n_size.GT.1) ) THEN
!---> for EV-parallelization & noco ( see comments at top )
lapw%nv_tot = lapw%nv(1) + lapw%nv(2)
IF (noco%l_ss) CALL juDFT_error("ev-|| & spin-spiral !",calledby ="hssphn")
ELSE
lapw%nv_tot = lapw%nv(iintsp)
ENDIF
!
!$OMP DO SCHEDULE(DYNAMIC,1)
DO kii = n_rank, lapw%nv_tot-1, n_size
ki = mod(kii,lapw%nv(iintsp)) + 1
nc = nc + 1
!$ nc= 1+ (kii-n_rank)/n_size
iii = nc*(nc-1)/2 * n_size - (nc-1)*(n_size - n_rank - 1)
! ii = nc*(nc+1)/2 * n_size - nc *(n_size - n_rank - 1)
IF (noco%l_ss.OR.noco%l_constr.OR.l_socfirst) THEN
kjmax = lapw%nv(jintsp)
ELSE
kjmax = ki
ENDIF
ski = (/lapw%k1(ki,iintsp),lapw%k2(ki,iintsp),lapw%k3(ki,iintsp)/) + qssbti
!---> legendre polynomials
plegend(:,0)=1.0
DO kj = 1,kjmax
plegend(kj,1) = dot_product(gk(kj,:,jintsp),gk(ki,:,iintsp))
END DO
DO l = 1,maxval(atoms%lmax) - 1
plegend(:,l+1) = fleg1(l)*plegend(:,1)*plegend(:,l) - fleg2(l)*plegend(:,l-1)
END DO
!include factor fl2p1 already here
DO l=0,maxval(atoms%lmax)
plegend(:kjmax,l)=plegend(:kjmax,l)*fl2p1(l)
ENDDO
!---> loop over equivalent atoms
n1 = 0
DO n = 1,atoms%ntype
IF (noco%l_noco) THEN
!---> pk non-collinear
!---> set up the spinors of this atom within global
!---> spin-coordinateframe
call hsmt_spinor(isp,n, noco, input,chi, chi11, chi21, chi22)
ENDIF
!---> pk non-collinear
n0 = n1 + 1
n1 = n1 + atoms%neq(n)
rph(:,1) = 0.0
cph(:,1) = 0.0
DO nn = n0,n1
tnn = tpi_const*atoms%taual(:,nn)
!---> set up phase factors
DO kj = 1,kjmax
rph(kj,1) = rph(kj,1) +&
cos(dot_product(ski-(/lapw%k1(kj,jintsp),lapw%k2(kj,jintsp),lapw%k3(kj,jintsp)/)+qssbtj,tnn))
#ifndef CPP_INVERSION
!---> if the system does not posses inversion symmetry
!---> the complex part of the exponential is needed.
cph(kj,1) = cph(kj,1) +&
sin(dot_product((/lapw%k1(kj,jintsp),lapw%k2(kj,jintsp),lapw%k3(kj,jintsp)/)+qssbtj-ski,tnn))
#endif
END DO
END DO
!---> update overlap and l-diagonal hamiltonian matrix
DO l = 0,atoms%lmax(n)
IF (noco%l_constr.or.l_socfirst) THEN
fjkiln = fj(ki,l,n,isp)
gjkiln = gj(ki,l,n,isp)*usdus%ddn(l,n,isp)
ELSE
fjkiln = fj(ki,l,n,iintsp)
gjkiln = gj(ki,l,n,iintsp)*usdus%ddn(l,n,isp)
ENDIF
!
!
IF ( noco%l_noco .AND. (.NOT. noco%l_ss) ) THEN
!---> pk non-collinear
#ifndef CPP_INVERSION
IF (noco%l_constr.or.l_socfirst) THEN
DO kj = 1,ki
fct = plegend(kj,l)*&
( fjkiln*fj(kj,l,n,isp) + gjkiln*gj(kj,l,n,isp) )
bbwa(kj) = cmplx(rph(kj,1),cph(kj,1))*fct
ENDDO
ELSE
DO kj = 1,ki
fct = plegend(kj,l)*&
( fjkiln*fj(kj,l,n,jintsp) + gjkiln*gj(kj,l,n,jintsp) )
bbwa(kj) = cmplx(rph(kj,1),cph(kj,1))*fct
ENDDO
ENDIF
!+||
IF ( kii+1.LE.lapw%nv(1) ) THEN
!---> spin-up spin-up part
CALL CPP_BLAS_caxpy(ki,chi11,bbwa,1,bb(iii+1),1)
!---> spin-down spin-up part, upper triangle.
!---> the help array is used to allow vectorization on
!---> Cray PVP systems. it is mapped onto the hamiltonian
!---> matrix after the setup is complete.
DO kj = 1,ki - 1
ii = iii + kj
bbhlp(ii)=bbhlp(ii)+conjg(bbwa(kj))*chi21
END DO
ENDIF
IF ( (kii+1.GT.lapw%nv(1)).OR.(n_size.EQ.1) ) THEN
IF (n_size.EQ.1) THEN
ii = (lapw%nv(1)+atoms%nlotot+ki-1)*(lapw%nv(1)+atoms%nlotot+ki)/2
ELSE
ii = iii
ENDIF
!---> spin-down spin-up part, lower triangle
CALL CPP_BLAS_caxpy(ki,chi21,bbwa,1,bb(ii+1),1)
!---> spin-down spin-down part
ii = ii + lapw%nv(1)+atoms%nlotot
CALL CPP_BLAS_caxpy(ki,chi22,bbwa,1,bb(ii+1),1)
ENDIF
!-||
!---> when fj and gj are available for both local spins
!---> (isp), add the contribution from the constraint
!---> B-field.
ELSEIF ( noco%l_noco .AND. noco%l_ss ) THEN
IF ( iintsp.EQ.2 .AND. jintsp.EQ.1 ) THEN
kjmax = lapw%nv(1)
ELSE
kjmax = ki
ENDIF
DO kj = 1,kjmax
fct = plegend(kj,l)* ( fjkiln*fj(kj,l,n,jintsp) +&
gjkiln*gj(kj,l,n,jintsp) )
bbwa(kj) = cmplx(rph(kj,1),cph(kj,1))*fct
ENDDO
IF ( iintsp.EQ.1 .AND. jintsp.EQ.1 ) THEN
!---> spin-up spin-up part
ii = (ki-1)*(ki)/2
DO kj = 1,ki
bb(ii+kj) = bb(ii+kj) + bbwa(kj)*chi11
ENDDO
ELSEIF ( iintsp.EQ.2 .AND. jintsp.EQ.2 ) THEN
!---> spin-down spin-down part
ii = (lapw%nv(1)+atoms%nlotot+ki-1)*(lapw%nv(1)+atoms%nlotot+ki)/2 +&
lapw%nv(1)+atoms%nlotot
DO kj = 1,ki
bb(ii+kj) = bb(ii+kj) + bbwa(kj)*chi22
ENDDO
ELSE
!---> spin-down spin-up part
ii = (lapw%nv(1)+atoms%nlotot+ki-1)*(lapw%nv(1)+atoms%nlotot+ki)/2
DO kj = 1,lapw%nv(1)
bb(ii+kj) = bb(ii+kj) + bbwa(kj)*chi21
ENDDO
ENDIF
#endif
!---> pk non-collinear
ELSE
DO kj = 1,ki
fct = plegend(kj,l)*( fjkiln*fj(kj,l,n,jintsp) + gjkiln*gj(kj,l,n,jintsp) )
ij = iii + kj
#ifdef CPP_INVERSION
bb(ij) = bb(ij) + rph(kj,1) * fct
!-APW
#else
bb(ij) = bb(ij) + cmplx(rph(kj,1),cph(kj,1))*fct
#endif
END DO
ENDIF
!---> end loop over l
enddo
!---> end loop over atom types (ntype)
enddo
!---> end loop over ki
enddo
!$OMP END DO
!---> end loops over interstitial spins
ENDDO ! jintsp = 1,iintsp
ENDDO ! iintsp = 1,nintsp
deallocate(plegend)
deallocate(rph,cph)
!$OMP END PARALLEL
!---> pk non-collinear
IF ( noco%l_noco .AND. (.NOT. noco%l_ss) ) THEN
!---> map the hamiltonian help array onto the hamitonian matrix
IF (n_size.EQ.1) THEN
DO ki = 1,lapw%nv(1)
ii = (ki-1)*(ki)/2
DO kj = 1,ki-1
ij = (lapw%nv(1)+atoms%nlotot+kj-1)*(lapw%nv(1)+atoms%nlotot+kj)/2 + ki
bb(ij) = bb(ij) + bbhlp(ii+kj)
ENDDO
ENDDO
ELSE
#ifdef CPP_MPI
CALL mingeselle(SUB_COMM,n_size,n_rank,lapw%nv, bbhlp, bb)
#endif
ENDIF
ENDIF
RETURN
END SUBROUTINE hsmt_overlap
END MODULE m_hsmt_overlap
MODULE m_hsmt_sph
use m_juDFT
implicit none
CONTAINS
SUBROUTINE hsmt_sph(dimension,atoms,SUB_COMM,n_size,n_rank,sphhar,isp,ab_dim,&
input,hlpmsize,noco,l_socfirst,cell,nintsp, lapw,el,usdus,vr,gk,rsoc,isigma, aa,bb,fj,gj)
#include"cpp_double.h"
USE m_constants, ONLY : fpi_const,tpi_const
USE m_sphbes
USE m_dsphbs
USE m_ylm
USE m_hsmt_socinit,ONLY:t_rsoc
USE m_hsmt_spinor
USE m_radovlp
#ifdef CPP_MPI
USE m_mingeselle
#endif
USE m_types
IMPLICIT NONE
TYPE(t_dimension),INTENT(IN):: dimension
TYPE(t_input),INTENT(IN) :: input
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_lapw),INTENT(INOUT) :: lapw!lpaw%nv_tot is updated
TYPE(t_usdus),INTENT(INOUT) :: usdus
! ..
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: isp,ab_dim
INTEGER, INTENT (IN) :: SUB_COMM,n_size,n_rank
INTEGER, INTENT (IN) :: hlpmsize,nintsp
LOGICAL, INTENT (IN) :: l_socfirst
! ..
! .. Array Arguments ..
REAL, INTENT (IN) :: el(0:atoms%lmaxd,atoms%ntypd,dimension%jspd)
REAL, INTENT (IN) :: vr(atoms%jmtd,0:sphhar%nlhd,atoms%ntypd,dimension%jspd)
REAL,INTENT(IN) :: gk(:,:,:)
COMPLEX,INTENT(IN) :: isigma(2,2,3)
TYPE(t_rsoc),INTENT(IN) :: rsoc
#ifdef CPP_INVERSION
REAL, INTENT (INOUT) :: aa(:),bb(:)!(matsize)
#else
COMPLEX, INTENT (INOUT) :: aa(:),bb(:)
#endif
REAL,INTENT(OUT) :: fj(:,0:,:,:),gj(:,0:,:,:)
! ..
! .. Local Scalars ..
REAL con1,ff,gg,gs,ws,tnn(3), elall,fct,fjkiln,gjkiln,ddnln,ski(3)
REAL apw_lo1,apw_lo2,apw1,w1
COMPLEX chi11,chi21,chi22,capw1
INTEGER ii,iii,ij,k,ki,kj,l,lo,n,n0,n1,nn,kjmax, nsp,iintsp,jintsp
INTEGER nc ,kii
! ..
! .. Local Arrays ..
REAL fb(0:atoms%lmaxd),fleg1(0:atoms%lmaxd),fleg2(0:atoms%lmaxd),fl2p1(0:atoms%lmaxd)
REAL gb(0:atoms%lmaxd)
REAL, ALLOCATABLE :: plegend(:,:)
REAL, ALLOCATABLE :: rph(:),cph(:)
! ..
con1 = fpi_const/sqrt(cell%omtil)
DO l = 0,atoms%lmaxd
fleg1(l) = real(l+l+1)/real(l+1)
fleg2(l) = real(l)/real(l+1)
fl2p1(l) = real(l+l+1)/fpi_const
END DO
DO iintsp = 1,nintsp
DO n = 1,atoms%ntype
DO k = 1,lapw%nv(iintsp)
gs = lapw%rk(k,iintsp)*atoms%rmt(n)
CALL sphbes(atoms%lmax(n),gs, fb)
CALL dsphbs(atoms%lmax(n),gs,fb, gb)
DO l = 0,atoms%lmax(n)
!---> set up wronskians for the matching conditions for each ntype
ws = con1/(usdus%uds(l,n,isp)*usdus%dus(l,n,isp)&
- usdus%us(l,n,isp)*usdus%duds(l,n,isp))
ff = fb(l)
gg = lapw%rk(k,iintsp)*gb(l)
! IF ( apw(l) ) THEN
! fj(k,l,n,iintsp) = 1.0*con1 * ff / usdus%us(l,n,isp)
! gj(k,l,n,iintsp) = 0.0d0
! ELSE
!---> in a spin-spiral calculation fj and gj are needed
!---> both interstitial spin directions at the same time
!---> In all other cases iintsp runs from 1 to 1.
fj(k,l,n,iintsp) = ws * ( usdus%uds(l,n,isp)*gg - usdus%duds(l,n,isp)*ff )
gj(k,l,n,iintsp) = ws * ( usdus%dus(l,n,isp)*ff - usdus%us(l,n,isp)*gg )
! ENDIF
ENDDO
ENDDO ! k = 1, lapw%nv
ENDDO ! n = 1,atoms%ntype
ENDDO ! iintsp = 1,nintsp
!
ALLOCATE(rph(dimension%nvd))
ALLOCATE(cph(dimension%nvd))
ALLOCATE(plegend(dimension%nvd,0:atoms%lmaxd))
plegend=0.0
plegend(:,0)=1.0
DO iintsp = 1,nintsp
DO jintsp = 1,iintsp
nc = 0 ! maybe IF (iintsp.EQ
lapw%nv_tot = lapw%nv(iintsp)
DO kii = n_rank, lapw%nv_tot-1, n_size
ki = mod(kii,lapw%nv(iintsp)) + 1
nc = nc + 1
!$ nc= 1+ (kii-n_rank)/n_size
iii = nc*(nc-1)/2 * n_size - (nc-1)*(n_size - n_rank - 1)
! ii = nc*(nc+1)/2 * n_size - nc *(n_size - n_rank - 1)
kjmax = ki
ski = (/lapw%k1(ki,iintsp),lapw%k2(ki,iintsp),lapw%k3(ki,iintsp)/)
!---> legendre polynomials
DO kj = 1,kjmax