Commit 02e8940e authored by Daniel Wortmann's avatar Daniel Wortmann

Finished merge with hsmt_simple branch. At least some tests run already.

Cleanup and bugfixes are required.
parent da2884ab
......@@ -336,8 +336,6 @@ CONTAINS
IF(l_cs.AND.jspin.EQ.1) CALL corespec_gaunt()
! calculation of core spectra (EELS) initializations -end-
ALLOCATE ( kveclo(atoms%nlotot) )
IF (mpi%irank==0) THEN
WRITE (6,FMT=8000) jspin
......@@ -570,7 +568,7 @@ CONTAINS
ikpt,jspin,zmat%nbasfcn,noco%l_ss,noco%l_noco,&
noccbd,n_start,n_end,&
ello,evdu,epar,&
lapw,wk,nbands,eig,zMat)
wk,nbands,eig,zMat)
#ifdef CPP_MPI
! Sinchronizes the RMA operations
if (l_evp) CALL MPI_BARRIER(mpi%mpi_comm,ie)
......
......@@ -19,7 +19,7 @@ CONTAINS
!>@author D. Wortmann
SUBROUTINE eigen(mpi,stars,sphhar,atoms,obsolete,xcpot,&
sym,kpts,DIMENSION, vacuum, input, cell, enpara_in,banddos, noco,jij, oneD,hybrid,&
it,eig_id,results,v,vx)
it,eig_id,results,inden,v,vx)
USE m_constants, ONLY : pi_const,sfp_const
USE m_types
USE m_lodpot
......@@ -61,6 +61,7 @@ CONTAINS
TYPE(t_kpts),INTENT(IN) :: kpts
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(INOUT) :: atoms!in u_setup n_u might be modified
TYPE(t_potden),INTENT(IN) :: inden
TYPE(t_potden),INTENT(INOUT) :: v,vx
#ifdef CPP_MPI
INCLUDE 'mpif.h'
......@@ -143,7 +144,7 @@ CONTAINS
!---> loop over spins
!---> set up k-point independent t(l'm',lm) matrices
!
CALL mt_setup(atoms,sym,sphhar,input,noco,enpara,v,mpi,results,DIMENSION,td,ud)
CALL mt_setup(atoms,sym,sphhar,input,noco,enpara,inden,v,mpi,results,DIMENSION,td,ud)
DO jsp = 1,MERGE(1,input%jspins,noco%l_noco)
k_loop:DO nk = mpi%n_start,kpts%nkpt,mpi%n_stride
......
......@@ -130,98 +130,7 @@ CONTAINS
COMPLEX,ALLOCATABLE:: ab(:,:),ab1(:,:),ab_select1(:,:),ab_select(:,:)
real :: rchi
<<<<<<< HEAD
!
!---> update hamiltonian and overlap matrices
nc = 0
IF ( noco%l_noco .AND. (n_size>1) ) THEN
lapw%nv_tot = lapw%nv(1) + lapw%nv(2)
ELSE
lapw%nv_tot = lapw%nv(iintsp)
ENDIF
kii=n_rank
DO WHILE(kii<lapw%nv_tot)
!DO kii = n_rank, nv_tot-1, n_size
ki = MOD(kii,lapw%nv(iintsp)) + 1
bsize=MIN(SIZE(aa_block,1),(lapw%nv(iintsp)-ki)/n_size+1) !Either use maximal blocksize or number of rows left to calculate
IF (bsize<1) EXIT !nothing more to do here
bsize2=bsize*n_size
bsize2=min(bsize2,lapw%nv(iintsp)-ki+1)
!aa_block(:bsize,:ki+bsize2-1)=matmul(a(ki:ki+bsize2-1:n_size,0:lmp,iintsp),conjg(transpose(ax(:ki+bsize2-1,0:lmp))))+ &
! matmul(b(ki:ki+bsize2-1:n_size,0:lmp,iintsp),conjg(transpose(bx(:ki+bsize2-1,0:lmp))))
IF (n_size==1) THEN !Make this a special case to avoid copy-in of a array
call zgemm("N","C",bsize,ki+bsize2-1,lmp+1,one,a(ki,0,iintsp),SIZE(a,1),ax(1,0),SIZE(ax,1),zero,aa_block,SIZE(aa_block,1))
call zgemm("N","C",bsize,ki+bsize2-1,lmp+1,one,b(ki,0,iintsp),SIZE(a,1),bx(1,0),SIZE(ax,1),one ,aa_block,SIZE(aa_block,1))
ELSE
CALL zgemm("N","C",bsize,ki+bsize2-1,lmp+1,one,a(ki:ki+bsize2-1:n_size,0:lmp,iintsp),SIZE(a(ki:ki+bsize2-1:n_size,0:lmp,iintsp),1),ax(1,0),SIZE(ax,1),zero,aa_block,SIZE(aa_block,1))
CALL zgemm("N","C",bsize,ki+bsize2-1,lmp+1,one,b(ki:ki+bsize2-1:n_size,0:lmp,iintsp),SIZE(a(ki:ki+bsize2-1:n_size,0:lmp,iintsp),1),bx(1,0),SIZE(ax,1),one,aa_block,SIZE(aa_block,1))
ENDIF
DO kb=1,bsize
IF ( noco%l_noco .AND. (.NOT. noco%l_ss) ) THEN
nc = 1+kii/n_size
ii = nc*(nc-1)/2*n_size-(nc-1)*(n_size-n_rank-1)
IF ( (n_size==1).OR.(kii+1<=lapw%nv(1)) ) THEN !
aahlp(ii+1:ii+ki) = aahlp(ii+1:ii+ki)+MATMUL(CONJG(ax(:ki,:lmp)),a(ki,:lmp,iintsp))+MATMUL(CONJG(bx(:ki,:lmp)),b(ki,:lmp,iintsp))
ELSE ! components for <2||2> block unused
aa_tmphlp(:ki) = MATMUL(CONJG(ax(:ki,:lmp)),a(ki,:lmp,iintsp))+MATMUL(CONJG(bx(:ki,:lmp)),b(ki,:lmp,iintsp))
!---> spin-down spin-down part
ij = ii + lapw%nv(1)
aa_c(ij+1:ij+ki)=aa_c(ij+1:ij+ki)+chi22*aa_tmphlp(:ki)
!---> spin-down spin-up part, lower triangle
ij = ii
aa_c(ij+1:ij+ki)=aa_c(ij+1:ij+ki)+chi21*aa_tmphlp(:ki)
ENDIF
!-||
ELSEIF ( noco%l_noco .AND. noco%l_ss ) THEN
IF ( iintsp==1 .AND. jintsp==1 ) THEN
!---> spin-up spin-up part
kjmax = ki
chihlp = chi11
ii = (ki-1)*(ki)/2
ELSEIF ( iintsp==2 .AND. jintsp==2 ) THEN
!---> spin-down spin-down part
kjmax = ki
chihlp = chi22
ii = (lapw%nv(1)+atoms%nlotot+ki-1)*(lapw%nv(1)+atoms%nlotot+ki)/2+&
lapw%nv(1)+atoms%nlotot
ELSE
!---> spin-down spin-up part
kjmax = lapw%nv(1)
chihlp = chi21
ii = (lapw%nv(1)+atoms%nlotot+ki-1)*(lapw%nv(1)+atoms%nlotot+ki)/2
ENDIF
aa_c(ii+1:ii+kjmax) = aa_c(ii+1:ii+kjmax) + chihlp*&
(MATMUL(CONJG(ax(:kjmax,:lmp)),a(ki,:lmp,iintsp))+MATMUL(CONJG(bx(:kjmax,:lmp)),b(ki,:lmp,iintsp)))
ELSE
nc = 1+kii/n_size
ii = nc*(nc-1)/2*n_size- (nc-1)*(n_size-n_rank-1)
if (l_real) THEN
aa_r(ii+1:ii+ki) = aa_r(ii+1:ii+ki) + aa_block(kb,:ki)
ELSE
aa_c(ii+1:ii+ki) = aa_c(ii+1:ii+ki) + aa_block(kb,:ki)
endif
!print*,ii,ki,kb
! IF (.not.apw(l)) THEN
!aa(ii+1:ii+ki) = aa(ii+1:ii+ki) + b(ki,lmp,iintsp)*bx(:ki)
! ENDIF
ENDIF
ki=ki+n_size
kii=kii+n_size
ENDDO
!---> end loop over ki
END DO
!---> end loops over interstitial spin
ENDDO
ENDDO
ENDIF ! atoms%invsat(na) = 0 or 1
!---> end loop over equivalent atoms
END DO
IF ( noco%l_noco .AND. (.NOT. noco%l_ss) ) CALL hsmt_hlptomat(atoms%nlotot,lapw%nv,sub_comm,chi11,chi21,chi22,aahlp,aa_c)
!---> end loop over atom types (ntype)
ENDDO ntyploop
=======
ALLOCATE(ab(MAXVAL(lapw%nv),2*atoms%lnonsph(n)*(atoms%lnonsph(n)+2)+2),ab1(lapw%nv(iintsp),2*atoms%lnonsph(n)*(atoms%lnonsph(n)+2)+2),ab_select(lapw%num_local_cols(jintsp),2*atoms%lnonsph(n)*(atoms%lnonsph(n)+2)+2))
>>>>>>> hsmt_simple
IF (iintsp.NE.jintsp) ALLOCATE(ab_select1(lapw%num_local_cols(jintsp),2*atoms%lnonsph(n)*(atoms%lnonsph(n)+2)+2))
......
......@@ -48,140 +48,12 @@ CONTAINS
CALL timestart("spherical setup")
<<<<<<< HEAD
! for Spin-orbit...
REAL, ALLOCATABLE :: dplegend(:,:)
REAL, ALLOCATABLE :: cross_k(:,:)
INTEGER :: j1,j2
COMPLEX :: isigma_x(2,2),isigma_y(2,2),isigma_z(2,2)
COMPLEX :: chi11so(2,2),chi21so(2,2),chi22so(2,2),angso(DIMENSION%nvd,2,2)
REAL :: tmp,tmp1,tmp2,tmp3
IF ( noco%l_noco .AND. (.NOT. noco%l_ss) ) ALLOCATE ( aahlp(hlpmsize),bbhlp(hlpmsize) )
! ..
con1 = fpi_const/SQRT(cell%omtil)
=======
>>>>>>> hsmt_simple
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
fl2p1bt(l) = fl2p1(l)*0.5
END DO
<<<<<<< HEAD
!---> calculate the overlap of the radial basis functions with different
!---> spin directions.
IF (noco%l_constr) THEN
ALLOCATE(uun21(0:atoms%lmaxd,atoms%ntype),udn21(0:atoms%lmaxd,atoms%ntype),&
dun21(0:atoms%lmaxd,atoms%ntype),ddn21(0:atoms%lmaxd,atoms%ntype) )
CALL rad_ovlp(atoms,usdus,input,vr,el(0:,:,:), uun21,udn21,dun21,ddn21)
ENDIF
!---> loop over each atom type
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) = .FALSE.
DO lo = 1,atoms%nlo(n)
IF (atoms%l_dulo(lo,n)) apw(l) = .TRUE.
ENDDO
IF ((input%l_useapw).AND.(atoms%lapw_l(n).GE.l)) apw(l) = .FALSE.
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
!
IF ( noco%l_noco .AND. (.NOT. noco%l_ss) ) THEN
!---> pk non-collinear
!---> initialize hamiltonian help array
aahlp=CMPLX(0.0,0.0)
bbhlp=CMPLX(0.0,0.0)
ENDIF
!$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(w1,apw_lo1,apw_lo2,ddnln,elall,fct,ij,apw1)&
!$OMP PRIVATE(cross_k,dplegend,chi,chi11,chi21,chi22,nsp,chj)&
!$OMP PRIVATE(isigma_x,isigma_y,isigma_z,j1,j2,chi11so,chi21so,chi22so)&
!$OMP PRIVATE(tmp1,tmp2,tmp3,tmp)&
!$OMP PRIVATE(aawa,bbwa,capw1,ii) IF (.not.l_socfirst)
!$ IF (l_socfirst) WRITE(*,*) "WARNING: first variation SOC does not work with OPENMP in hsmt_sph"
!$ IF (l_socfirst) WRITE(*,*) " switching off openmp parallelization"
ALLOCATE(rph(DIMENSION%nvd,ab_dim))
ALLOCATE(cph(DIMENSION%nvd,ab_dim))
ALLOCATE(plegend(DIMENSION%nvd,0:atoms%lmaxd))
IF (l_socfirst)THEN
ALLOCATE ( dplegend(DIMENSION%nvd,0:atoms%lmaxd),cross_k(DIMENSION%nvd,3))
dplegend(:,0)=0.e0
dplegend(:,1)=1.e0
ENDIF
plegend=0.0
plegend(:,0)=1.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-parallelization & spin-spiral !",calledby ="hsmt_sph")
ELSE
lapw%nv_tot = lapw%nv(iintsp)
ENDIF
=======
!$OMP PARALLEL DEFAULT(SHARED)&
!$OMP PRIVATE(kii,ki,ski,kj,plegend,l)&
!$OMP PRIVATE(cph,nn,tnn,fjkiln,gjkiln)&
......@@ -219,7 +91,6 @@ CONTAINS
DO l = 0,atoms%lmax(n)
fjkiln = fj(ki,l,iintsp)
gjkiln = gj(ki,l,iintsp)
>>>>>>> hsmt_simple
!
IF (input%l_useapw) THEN
w1 = 0.5 * ( usdus%uds(l,n,isp)*usdus%dus(l,n,isp) + &
......@@ -264,331 +135,7 @@ CONTAINS
hmat%data_c(kj,kii)=hmat%data_c(kj,kii) + capw1
ENDIF
END DO
<<<<<<< HEAD
!---> 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,l_socfirst,&
isigma,isigma_x,isigma_y,isigma_z,chi11so,chi21so,chi22so,angso,chj,cross_k)
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
!!$OMP SIMD PRIVATE(tmp1,tmp2,tmp3,tmp)
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))
tmp1 = (ski(1)-lapw%k1(kj,jintsp)+qssbtj(1)) * tnn(1)
tmp2 = (ski(2)-lapw%k2(kj,jintsp)+qssbtj(2)) * tnn(2)
tmp3 = (ski(3)-lapw%k3(kj,jintsp)+qssbtj(3)) * tnn(3)
tmp = cos(tmp1 + tmp2 + tmp3)
rph(kj,1) = rph(kj,1) + tmp
END DO
IF (.NOT.sym%invs) THEN
!---> if the system does not posses inversion symmetry
!---> the complex part of the exponential is needed.
!!$OMP SIMD PRIVATE(tmp1,tmp2,tmp3,tmp)
DO kj = 1,kjmax
! cph(kj,1) = cph(kj,1) +&
! SIN(DOT_PRODUCT((/lapw%k1(kj,jintsp),lapw%k2(kj,jintsp),lapw%k3(kj,jintsp)/)+qssbtj-ski,tnn))
tmp1 = (lapw%k1(kj,jintsp)+qssbtj(1)-ski(1)) * tnn(1)
tmp2 = (lapw%k2(kj,jintsp)+qssbtj(2)-ski(2)) * tnn(2)
tmp3 = (lapw%k3(kj,jintsp)+qssbtj(3)-ski(3)) * tnn(3)
tmp = sin(tmp1 + tmp2 + tmp3)
cph(kj,1) = cph(kj,1) + tmp
END DO
ENDIF
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)
ELSE
fjkiln = fj(ki,l,n,iintsp)
gjkiln = gj(ki,l,n,iintsp)
ENDIF
!
w1 = 0.5 * ( usdus%uds(l,n,isp)*usdus%dus(l,n,isp) + &
usdus%us(l,n,isp)*usdus%duds(l,n,isp) )
apw_lo1 = fl2p1(l) * 0.5 * atoms%rmt(n)**2 * ( gjkiln * w1 +&
fjkiln * usdus%us(l,n,isp) * usdus%dus(l,n,isp) )
apw_lo2 = fl2p1(l) * 0.5 * atoms%rmt(n)**2 * ( fjkiln * w1 +&
gjkiln * usdus%uds(l,n,isp) * usdus%duds(l,n,isp) )
!
ddnln = usdus%ddn(l,n,isp)
elall = el(l,n,isp)
IF ( noco%l_noco .AND. (.NOT. noco%l_ss) ) THEN
!---> pk non-collinear
IF (noco%l_constr.OR.l_socfirst) THEN
DO kj = 1,ki
fct = plegend(kj,l)*fl2p1(l)*&
( fjkiln*fj(kj,l,n,isp) + gjkiln*gj(kj,l,n,isp)*ddnln )
bbwa(kj) = CMPLX(rph(kj,1),cph(kj,1))*fct
aawa(kj) = CMPLX(rph(kj,1),cph(kj,1)) * ( &
fct*elall + plegend(kj,l)*fl2p1bt(l)*&
( fjkiln*gj(kj,l,n,isp) + gjkiln*fj(kj,l,n,isp) ) )
IF (input%l_useapw) THEN
capw1 = CMPLX(rph(kj,1),cph(kj,1))*plegend(kj,l)&
* ( apw_lo1 * fj(kj,l,n,isp) + apw_lo2 * gj(kj,l,n,isp) )
aawa(kj) = aawa(kj) + capw1
ENDIF
ENDDO
ELSE
DO kj = 1,ki
fct = plegend(kj,l)*fl2p1(l)*&
( fjkiln*fj(kj,l,n,jintsp) + gjkiln*gj(kj,l,n,jintsp)*ddnln )
bbwa(kj) = CMPLX(rph(kj,1),cph(kj,1))*fct
aawa(kj) = CMPLX(rph(kj,1),cph(kj,1)) * (&
fct*elall + plegend(kj,l)*fl2p1bt(l)*&
( fjkiln*gj(kj,l,n,jintsp) + gjkiln*fj(kj,l,n,jintsp) ) )
IF (input%l_useapw) THEN
capw1 = CMPLX(rph(kj,1),cph(kj,1))*plegend(kj,l)&
* ( apw_lo1 * fj(kj,l,n,jintsp) + apw_lo2 * gj(kj,l,n,jintsp) )
aawa(kj) = aawa(kj) + capw1
ENDIF
ENDDO
ENDIF
!+||
IF ( kii+1.LE.lapw%nv(1) ) THEN
!---> spin-up spin-up part
CALL CPP_BLAS_caxpy(ki,chi11,bbwa,1,hamOvlp%b_c(iii+1),1)
CALL CPP_BLAS_caxpy(ki,chi11,aawa,1,hamOvlp%a_c(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
aahlp(ii)=aahlp(ii)+CONJG(aawa(kj))*chi21
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,hamOvlp%b_c(ii+1),1)
CALL CPP_BLAS_caxpy(ki,chi21,aawa,1,hamOvlp%a_c(ii+1),1)
!---> spin-down spin-down part
ii = ii + lapw%nv(1)+atoms%nlotot
CALL CPP_BLAS_caxpy(ki,chi22,bbwa,1,hamOvlp%b_c(ii+1),1)
CALL CPP_BLAS_caxpy(ki,chi22,aawa,1,hamOvlp%a_c(ii+1),1)
ENDIF
!-||
!---> when fj and gj are available for both local spins
!---> (isp), add the contribution from the constraint
!---> B-field.
IF (noco%l_constr .AND. (isp .EQ. 2)) THEN
DO nsp = 1,input%jspins
IF (nsp.EQ.1) THEN
DO kj = 1,lapw%nv(1)
aawa(kj) = (-0.5)*CMPLX(noco%b_con(1,n),noco%b_con(2,n))*&
CMPLX(rph(kj,1),cph(kj,1))*&
plegend(kj,l)*fl2p1(l)*(&
+ fj(ki,l,n,2)*fj(kj,l,n,1)*uun21(l,n)&
+ fj(ki,l,n,2)*gj(kj,l,n,1)*udn21(l,n)&
+ gj(ki,l,n,2)*fj(kj,l,n,1)*dun21(l,n)&
+ gj(ki,l,n,2)*gj(kj,l,n,1)*ddn21(l,n))
ENDDO
ELSE
DO kj = 1,lapw%nv(1)
aawa(kj) = (-0.5)*CMPLX(noco%b_con(1,n),-noco%b_con(2,n))*&
CMPLX(rph(kj,1),cph(kj,1))*&
plegend(kj,l)*fl2p1(l)*(&
+ fj(ki,l,n,1)*fj(kj,l,n,2)*uun21(l,n)&
+ fj(ki,l,n,1)*gj(kj,l,n,2)*dun21(l,n)&
+ gj(ki,l,n,1)*fj(kj,l,n,2)*udn21(l,n)&
+ gj(ki,l,n,1)*gj(kj,l,n,2)*ddn21(l,n))
ENDDO
ENDIF
!---> spin-up spin-up part
ii = (ki-1)*(ki)/2
DO kj = 1,ki
hamOvlp%a_c(ii+kj) = hamOvlp%a_c(ii+kj) + aawa(kj)*chj(nsp,1,1,n)
ENDDO
!---> 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
hamOvlp%a_c(ii+kj) = hamOvlp%a_c(ii+kj) + aawa(kj)*chj(nsp,2,2,n)
ENDDO
!---> 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)
hamOvlp%a_c(ii+kj) = hamOvlp%a_c(ii+kj) + aawa(kj)*chj(nsp,2,1,n)
ENDDO
ENDDO
ENDIF
! Add spin-orbit coupling
IF (isp.EQ.2) THEN
IF ((l.GT.0).AND.l_socfirst) THEN
DO j1 = 1,input%jspins
DO j2 = 1,input%jspins
DO kj = 1,lapw%nv(1)
aawa(kj)=CMPLX(rph(kj,1),cph(kj,1))*(&
+ fj(ki,l,n,j1)*fj(kj,l,n,j2)*rsoc%rsopp(n,l,j1,j2)&
+ fj(ki,l,n,j1)*gj(kj,l,n,j2)*rsoc%rsopdp(n,l,j1,j2)&
+ gj(ki,l,n,j1)*fj(kj,l,n,j2)*rsoc%rsoppd(n,l,j1,j2)&
+ gj(ki,l,n,j1)*gj(kj,l,n,j2)*rsoc%rsopdpd(n,l,j1,j2))&
*dplegend(kj,l)*fl2p1(l)*angso(kj,j1,j2)
ENDDO
IF (n_size.EQ.1) THEN
!---> spin-up spin-up part
ii = (ki-1)*(ki)/2
DO kj = 1,ki
hamOvlp%a_c(ii+kj) = hamOvlp%a_c(ii+kj) + aawa(kj)*chi11so(j1,j2)
ENDDO
!---> 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
hamOvlp%a_c(ii+kj) = hamOvlp%a_c(ii+kj) + aawa(kj)*chi22so(j1,j2)
ENDDO
!---> 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)
hamOvlp%a_c(ii+kj) = hamOvlp%a_c(ii+kj) + aawa(kj)*chi21so(j1,j2)
ENDDO
ELSE ! eigenvalue parallelization
IF ( kii+1.LE.lapw%nv(1) ) THEN
!---> spin-up spin-up part
CALL CPP_BLAS_caxpy(ki,chi11so(j1,j2),aawa,1, hamOvlp%a_c(iii+1),1)
!---> spin-down spin-up part, upper triangle.
DO kj = 1,ki - 1
ii = iii + kj
aahlp(ii) = aahlp(ii) + CONJG(aawa(kj))*chi21so(j2,j1)
END DO
ENDIF
IF (kii+1.GT.lapw%nv(1)) THEN
ii = iii
!---> spin-down spin-up part, lower triangle
CALL CPP_BLAS_caxpy(ki,chi21so(j1,j2),aawa,1, hamOvlp%a_c(ii+1),1)
!---> spin-down spin-down part
ii = ii + lapw%nv(1)+atoms%nlotot
CALL CPP_BLAS_caxpy(ki,chi22so(j1,j2),aawa,1, hamOvlp%a_c(ii+1),1)
ENDIF
ENDIF ! eigenvalue par.
ENDDO ! j2
ENDDO ! j1
ENDIF ! ( l > 0 ) & socfirst
ENDIF ! (isp = 2)
! End spin-orbit
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)*fl2p1(l)* ( fjkiln*fj(kj,l,n,jintsp) +&
gjkiln*gj(kj,l,n,jintsp)*ddnln )
bbwa(kj) = CMPLX(rph(kj,1),cph(kj,1))*fct
aawa(kj) = CMPLX(rph(kj,1),cph(kj,1)) * ( &
fct*elall + plegend(kj,l)*fl2p1bt(l)*&
( fjkiln*gj(kj,l,n,jintsp) + gjkiln*fj(kj,l,n,jintsp) ) )
ENDDO
IF ( iintsp.EQ.1 .AND. jintsp.EQ.1 ) THEN
!---> spin-up spin-up part
ii = (ki-1)*(ki)/2
DO kj = 1,ki
hamOvlp%b_c(ii+kj) = hamOvlp%b_c(ii+kj) + bbwa(kj)*chi11
hamOvlp%a_c(ii+kj) = hamOvlp%a_c(ii+kj) + aawa(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
hamOvlp%b_c(ii+kj) = hamOvlp%b_c(ii+kj) + bbwa(kj)*chi22
hamOvlp%a_c(ii+kj) = hamOvlp%a_c(ii+kj) + aawa(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)
hamOvlp%b_c(ii+kj) = hamOvlp%b_c(ii+kj) + bbwa(kj)*chi21
hamOvlp%a_c(ii+kj) = hamOvlp%a_c(ii+kj) + aawa(kj)*chi21
ENDDO
ENDIF
!---> pk non-collinear
ELSE
IF (l_real) THEN
DO kj = 1,ki
fct = plegend(kj,l)*fl2p1(l)*&
( fjkiln*fj(kj,l,n,jintsp) + gjkiln*gj(kj,l,n,jintsp)*ddnln )
ij = iii + kj
hamOvlp%b_r(ij) = hamOvlp%b_r(ij) + rph(kj,1) * fct
hamOvlp%a_r(ij) = hamOvlp%a_r(ij) + rph(kj,1) * ( fct * elall + plegend(kj,l) * fl2p1bt(l) *&
( fjkiln*gj(kj,l,n,jintsp) + gjkiln*fj(kj,l,n,jintsp) ) )
!+APW
IF (input%l_useapw) THEN
apw1 = rph(kj,1) * plegend(kj,l) * &
( apw_lo1 * fj(kj,l,n,iintsp) + apw_lo2 * gj(kj,l,n,iintsp) )
hamOvlp%a_r(ij) = hamOvlp%a_r(ij) + apw1
ENDIF
!-APW
ENDDO
ELSE
DO kj = 1,ki
fct = plegend(kj,l)*fl2p1(l)*&
( fjkiln*fj(kj,l,n,jintsp) + gjkiln*gj(kj,l,n,jintsp)*ddnln )
ij = iii + kj
hamOvlp%b_c(ij) = hamOvlp%b_c(ij) + CMPLX(rph(kj,1),cph(kj,1))*fct
hamOvlp%a_c(ij) = hamOvlp%a_c(ij) + CMPLX(rph(kj,1),cph(kj,1)) * (fct*elall + plegend(kj,l)*fl2p1bt(l) *&
( fjkiln*gj(kj,l,n,jintsp) + gjkiln*fj(kj,l,n,jintsp) ) )
IF (input%l_useapw) THEN