m_pawxc/pawxc [ Functions ]

Start from the density or spin-density, and compute xc correlation
potential and energies inside a paw sphere.
USE THE DENSITY OVER A WHOLE SPHERICAL GRID (r,theta,phi)
Driver of XC functionals.

INPUTS

corexc(nrad)=core density on radial grid
ixc= choice of exchange-correlation scheme
lm_size=size of density array rhor (see below)
lmselect(lm_size)=select the non-zero LM-moments of input density rhor
nhat(nrad,lm_size,nspden)=compensation density
(total in 1st half and spin-up in 2nd half if nspden=2)
nkxc=second dimension of the kxc array. If /=0, the exchange-correlation kernel must be computed
non_magnetic_xc= true if usepawu==4
nrad=size of radial mesh for densities/potentials (might be different from pawrad%mesh_size)
nspden=number of spin-density components
option=0 compute both XC energies (direct+double-counting) and potential
1 compute only XC potential
2 compute only XC energies (direct+double-counting)
3 compute only XC energy by direct scheme
4 compute only XC energy by direct scheme for spherical part of the density
5 compute only XC potential for spherical part of the density
pawang <type(pawang_type)>=paw angular mesh and related data
pawrad <type(pawrad_type)>=paw radial mesh and related data
rhor(nrad,lm_size,nspden)=electron density in real space in electrons/bohr**3
(total in 1st half and spin-up in 2nd half if nspden=2)
usecore= 1 if core density has to be used in Exc/Vxc ; 0 otherwise
usexcnhat= 0 if compensation density does not have to be used
1 if compensation density has to be used in double counting energy term only
2 if compensation density (nhat) has to be used in Exc/Vxc and double counting energy term
xclevel= XC functional level
xc_denpos= lowest allowed density (usually for the computation of the XC functionals)

m_pawxc/pawxc_drivexc_libxc [ Functions ]

5602 subroutine pawxc_drivexc_libxc()
5603 5604 5605 !This section has been created automatically by the script Abilint (TD).
5606 !Do not modify the following lines by hand.
5607 #undef ABI_FUNC
5608 #define ABI_FUNC 'pawxc_drivexc_libxc'
5609 !End of the abilint section
5610 5611 implicit none
5612 5613 ! *************************************************************************
5614 5615 !Check the compatibility of input arguments
5616 if (libxc_functionals_ismgga()) then
5617 msg='MGGA is not yet coded in pawxc_drivexc_wrapper/LIBXC'
5618 MSG_ERROR(msg)
5619 end if
5620 if (ixc>=0) then
5621 msg='ixc argument should be negative!'
5622 MSG_BUG(msg)
5623 end if
5624 if (ixc/=libxc_functionals_ixc()) then
5625 msg='The value of ixc differs from the one used to initialize the functional!'
5626 MSG_BUG(msg)
5627 end if
5628 if ((order<1.and.order/=-2).or.order>4) then
5629 msg='The only allowed values for order are 1, 2, -2, or 3!'
5630 MSG_BUG(msg)
5631 end if
5632 if ((order**2>1).and.(.not.present(dvxc))) then
5633 msg='The value of order is not compatible with the presence of the array dvxc!'
5634 MSG_BUG(msg)
5635 end if
5636 if ((order==3).and.(.not.present(d2vxc))) then
5637 msg='The value of order is not compatible with the presence of the array d2vxc!'
5638 MSG_BUG(msg)
5639 end if
5640 if (libxc_functionals_isgga()) then
5641 if ((.not.present(grho2)).or.(.not.present(vxcgrho)).or.(nvxcgrho==0)) then
5642 write(msg,'(3a)') 'At least one of the functionals is a GGA,',ch10, &
5643 & 'but not all the necessary optional arguments are present.'
5644 MSG_BUG(msg)
5645 end if
5646 if (ngr2==0.or.nvxcgrho/=3) then
5647 msg='The values of nvxcgrho or ngr2 are not compatible with GGA!'
5648 MSG_BUG(msg)
5649 end if
5650 end if
5651 5652 !Call LibXC routines
5653 if (libxc_functionals_isgga()) then
5654 if (order**2<=1) then
5655 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho,exc,vxcrho,&
5656 & grho2=grho2,vxcgr=vxcgrho)
5657 else
5658 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho,exc,vxcrho,&
5659 & grho2=grho2,vxcgr=vxcgrho,dvxc=dvxc)
5660 end if
5661 else
5662 if (order**2<=1) then
5663 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho,exc,vxcrho)
5664 else if (order**2<=4) then
5665 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho,exc,vxcrho,&
5666 & dvxc=dvxc)
5667 else
5668 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho,exc,vxcrho,&
5669 & dvxc=dvxc,d2vxc=d2vxc)
5670 end if
5671 end if
5672 5673 end subroutine pawxc_drivexc_libxc

m_pawxc/pawxc_mkdenpos_wrapper [ Functions ]

Make a ground-state density positive everywhere :
when the density (or spin-density) is smaller than xc_denpos,
set it to the value of xc_denpos

INPUTS

nfft=(effective) number of FFT grid points (for this processor)
nspden=number of spin-density components (max. 2)
option=0 if density rhonow is stored as (up,dn)
1 if density rhonow is stored as (up+dn,up)
Active only when nspden=2
xc_denpos= lowest allowed density (usually for the computation of the XC functionals)

OUTPUT

SIDE EFFECTS

Input/output
iwarn=At input: iwarn=0 a warning will be printed when rho is negative
iwarn>0 no warning will be printed out
At output: iwarn is increased by 1
rhonow(nfft,nspden)=electron (spin)-density in real space,
either on the unshifted grid (if ishift==0,
then equal to rhor),or on the shifted grid

m_pawxc/pawxc_rotate_mag [ Functions ]

Project (rotate) a non-collinear density (stored as density+magn.)
on a magnetization and give a collinear density (stored as [up,dn] or [up+dn,up]).

INPUTS

rho_in(vectsize,4)=input non-collinear density and magnetization
mag(vectsize,3)=magnetization used for projection
vectsize=size of vector fields
[rho_out_format]= 1=rho_out is stored as [up,dn]
2=rho_out is stored as [up+dn,up]
Default=1

OUTPUT

rho_out(vectsize,2)=output (projected, collinear) density
[mag_norm_out(vectsize)]= --optional-- norm of mag(:) at each point of the grid

ndvxc size of the array dvxc(npts,ndvxc) for allocation
ngr2 size of the array grho2_updn(npts,ngr2) for allocation
nd2vxc size of the array d2vxc(npts,nd2vxc) for allocation
nvxcdgr size of the array dvxcdgr(npts,nvxcdgr) for allocation

m_pawxc/pawxc_xcmult_wrapper [ Functions ]

In the case of GGA, multiply the different gradient of spin-density
by the derivative of the XC functional with respect
to the norm of the gradient, then divide it by the
norm of the gradient

INPUTS

depsxc(nfft,nspgrad)=derivative of Exc with respect to the (spin-)density,
or to the norm of the gradient of the (spin-)density,
further divided by the norm of the gradient of the (spin-)density
The different components of depsxc will be
for nspden=1, depsxc(:,1)=d(rho.exc)/d(rho)
and if ngrad=2, depsxc(:,2)=1/2*1/|grad rho_up|*d(rho.exc)/d(|grad rho_up|)
+ 1/|grad rho|*d(rho.exc)/d(|grad rho|)
(do not forget : |grad rho| /= |grad rho_up| + |grad rho_down|
for nspden=2, depsxc(:,1)=d(rho.exc)/d(rho_up)
depsxc(:,2)=d(rho.exc)/d(rho_down)
and if ngrad=2, depsxc(:,3)=1/|grad rho_up|*d(rho.exc)/d(|grad rho_up|)
depsxc(:,4)=1/|grad rho_down|*d(rho.exc)/d(|grad rho_down|)
depsxc(:,5)=1/|grad rho|*d(rho.exc)/d(|grad rho|)
nfft=(effective) number of FFT grid points (for this processor)
ngrad = must be 2
nspden=number of spin-density components
nspgrad=number of spin-density and spin-density-gradient components

OUTPUT

(see side effects)

SIDE EFFECTS

rhonow(nfft,nspden,ngrad*ngrad)=
at input :
electron (spin)-density in real space and its gradient,
either on the unshifted grid (if ishift==0,
then equal to rhor), or on the shifted grid
rhonow(:,:,1)=electron density in electrons/bohr**3
rhonow(:,:,2:4)=gradient of electron density in el./bohr**4
at output :
rhonow(:,:,2:4) has been multiplied by the proper factor,
described above.

m_pawxc/pawxcm [ Functions ]

Start from the density or spin-density, and compute xc correlation
potential and energies inside a paw sphere.
LDA+GGA - USE A DEVELOPMENT OF THE DENSITY OVER (L,M) MOMENTS
Driver of XC functionals.

INPUTS

corexc(nrad)=core density on radial grid
exexch= choice of <<<local>>> exact exchange. Active if exexch=3 (only for PBE)
ixc= choice of exchange-correlation scheme
lm_size=size of density array rhor (see below)
lmselect(lm_size)=select the non-zero LM-moments of input density rhor
nhat(nrad,lm_size,nspden)=compensation density
(total in 1st half and spin-up in 2nd half if nspden=2)
nkxc=second dimension of the kxc array. If /=0, the exchange-correlation kernel must be computed
non_magnetic_xc= true if usepawu==4
nrad=size of radial mesh for densities/potentials (might be different from pawrad%mesh_size)
nspden=number of spin-density components
option=0 compute both XC energies (direct+double-counting) and potential (and Kernel)
1 compute only XC potential (and Kernel)
2 compute only XC energies (direct+double-counting)
3 compute only XC energy by direct scheme
4 compute only XC energy by direct scheme for spherical part of the density
5 compute only XC potential (and Kernel) for spherical part of the density
pawang <type(pawang_type)>=paw angular mesh and related data
pawrad <type(pawrad_type)>=paw radial mesh and related data
pawxcdev=order of Vxc development
rhor(nrad,lm_size,nspden)=electron density in real space in electrons/bohr**3
(total in 1st half and spin-up in 2nd half if nspden=2)
usecore= 1 if core density has to be used in Exc/Vxc ; 0 otherwise
usexcnhat= 0 if compensation density does not have to be used
1 if compensation density has to be used in double counting energy term only
2 if compensation density (nhat) has to be used in Exc/Vxc and double counting energy term
xclevel= XC functional level
xc_denpos= lowest allowed density (usually for the computation of the XC functionals)

m_pawxc/pawxcmpositron [ Functions ]

Compute electron-positron correlation potential and energies inside a PAW sphere
LDA+GGA - USE A DEVELOPMENT OF THE DENSITY OVER (L,M) MOMENTS
Driver of XC functionals.

INPUTS

calctype=type of electron-positron calculation:
calctype=1 : positron in electronic density
calctype=2 : electrons in positronic density
corexc(nrad)=electron core density on radial grid
ixcpositron=choice of electron-positron XC scheme
lm_size=size of density array rhor (see below)
lmselect (lm_size)=select the non-zero LM-moments of input density rhor (see below)
lmselect_ep(lm_size)=select the non-zero LM-moments of input density rhor_ep (see below)
nhat (nrad,lm_size,nspden)=compensation density corresponding to rhor
nhat_ep(nrad,lm_size,nspden)=compensation density corresponding to rhor_ep
nrad=size of radial mesh for densities/potentials (might be different from pawrad%mesh_size)
nspden=number of spin-density components
option=0 compute both XC energies (direct+double-counting) and potential
1 compute only XC potential
2 compute only XC energies (direct+double-counting)
3 compute only XC energy by direct scheme
4 compute only XC energy by direct scheme for spherical part of the density
pawang <type(pawang_type)>=paw angular mesh and related data
pawrad <type(pawrad_type)>=paw radial mesh and related data
pawxcdev=order of Vxc development
posdensity0_limit=True if we are in the zero positron density limit
rhor(nrad,lm_size,nspden)=electron (or positron) density in real space
(total in 1st half and spin-up in 2nd half if nspden=2)
Contents depends on calctype value:
calctype=1: rhor is the positronic density
calctype=2: rhor is the electronic density
rhor_ep(nrad,lm_size,nspden)=electron (or positron) density in real space
(total in 1st half and spin-up in 2nd half if nspden=2)
Contents depends on calctype value:
calctype=1: rhor_ep is the electronic density
calctype=2: rhor_ep is the positronic density
usecore= 1 if core density has to be used in Exc/Vxc for the electronic density ; 0 otherwise
usexcnhat= 0 if compensation density does not have to be used
1 if compensation density has to be used in double counting energy term only
2 if compensation density (nhat) has to be used in Exc/Vxc and double counting energy term
xc_denpos= lowest allowed density (usually for the computation of the XC functionals)

m_pawxc/pawxcpositron [ Functions ]

Compute electron-positron correlation potential and energies inside a PAW sphere
LDA ONLY - USE THE DENSITY OVER A WHOLE SPHERICAL GRID (r,theta,phi)
Driver of XC functionals.

INPUTS

calctype=type of electronpositron calculation:
calctype=1 : positron in electronic density
calctype=2 : electrons in positronic density
corexc(nrad)=electron core density on radial grid
ixcpositron=choice of electron-positron XC scheme
lm_size=size of density array rhor (see below)
lmselect (lm_size)=select the non-zero LM-moments of input density rhor (see below)
lmselect_ep(lm_size)=select the non-zero LM-moments of input density rhor_ep (see below)
nhat (nrad,lm_size,nspden)=compensation density corresponding to rhor
nhat_ep(nrad,lm_size,nspden)=compensation density corresponding to rhor_ep
nrad=size of radial mesh for densities/potentials (might be different from pawrad%mesh_size)
nspden=number of spin-density components
option=0 compute both XC energies (direct+double-counting) and potential
1 compute only XC potential
2 compute only XC energies (direct+double-counting)
3 compute only XC energy by direct scheme
4 compute only XC energy by direct scheme for spherical part of the density
pawang <type(pawang_type)>=paw angular mesh and related data
pawrad <type(pawrad_type)>=paw radial mesh and related data
posdensity0_limit=True if we are in the zero positron density limit
rhor(nrad,lm_size,nspden)=electron (or positron) density in real space
(total in 1st half and spin-up in 2nd half if nspden=2)
Contents depends on calctype value:
calctype=1: rhor is the positronic density
calctype=2: rhor is the electronic density
rhor_ep(nrad,lm_size,nspden)=electron (or positron) density in real space
(total in 1st half and spin-up in 2nd half if nspden=2)
Contents depends on calctype value:
calctype=1: rhor_ep is the electronic density
calctype=2: rhor_ep is the positronic density
usecore= 1 if core density has to be used in Exc/Vxc for the electronic density ; 0 otherwise
usexcnhat= 0 if compensation density does not have to be used
1 if compensation density has to be used in double counting energy term only
2 if compensation density (nhat) has to be used in Exc/Vxc and double counting energy term
xc_denpos= lowest allowed density (usually for the computation of the XC functionals)

m_pawxc/pawxcsum [ Functions ]

Compute useful sums of moments of densities needed to compute on-site contributions to XC energy and potential
First order sums:
Sum1(1)=Sum_L{Rho1_L(r)**2}
Sum1(2)=Sum_L{Rho1_L(r)*Rho2_L(r)}
Sum1(3)=Sum_L{Rho2_L(r)**2}
With L>0
Second order sums:
Sum2(L,1)=Sum_L1_L2{Rho1_L1(r)*Rho1_L1(r)*Gaunt_(L,L1,L2)}
Sum2(L,2)=Sum_L1_L2{Rho1_L1(r)*Rho2_L2(r)*Gaunt_(L,L1,L2)}
Sum2(L,3)=Sum_L1_L2{Rho2_L2(r)*Rho2_L2(r)*Gaunt_(L,L1,L2)}
With L1>0, L2>0

INPUTS

cplex1=if 1, density Rho1 is REAL, if 2, COMPLEX
cplex2=if 1, density Rho2 is REAL, if 2, COMPLEX
cplexsum=if 1, output sums (Sum1 and Sum2) are REAL, if 2, COMPLEX
lmselect1(lm_size)=select the non-zero LM-moments of input density Rho1
lmselect2(lm_size)=select the non-zero LM-moments of input density Rho2
lm_size=number of moments of the density
nrad=number of radial points
nsums=number of sums to compute:
if nsums=1, computes only
Sum1(1)=Sum_L{Rho1_L(r)*Rho2_L(r)}
Sum2(L,1)=Sum_L1_L2{Rho1_L1(r)*Rho2_L2(r)*Gaunt_(L,L1,L2)}
if nsums=3, computes all sums (Sum1(1:3), Sum2(1:3)
option= 1: compute first order sums
2: compute first and second order sums
pawang <type(pawang_type)>=paw angular mesh and related data
rho1(cplex1*nrad,lm_size)=moments of first density on each radial point
rho2(cplex2*nrad,lm_size)=moments of 2nd density on each radial point

OUTPUT

sum1(cplexsum*nrad,nsums)=first order sums
=== if option>=2
sum2(cplexsum*nrad,lm_size,nsums)=second order sums

3371 subroutine pawxcsum(cplex1,cplex2,cplexsum,lmselect1,lmselect2,lm_size,nrad,nsums,&
3372 & option,pawang,rho1,rho2,sum1,sum2)
3373 3374 3375 !This section has been created automatically by the script Abilint (TD).
3376 !Do not modify the following lines by hand.
3377 #undef ABI_FUNC
3378 #define ABI_FUNC 'pawxcsum'
3379 !End of the abilint section
3380 3381 implicit none
3382 3383 !Arguments ------------------------------------
3384 !scalars
3385 integer,intent(in) :: cplex1,cplex2,cplexsum,lm_size,nrad,nsums,option
3386 !arrays
3387 logical,intent(in) :: lmselect1(lm_size),lmselect2(lm_size)
3388 real(dp),intent(in) :: rho1(cplex1*nrad,lm_size),rho2(cplex2*nrad,lm_size)
3389 real(dp),intent(out) :: sum1(cplexsum*nrad,nsums),sum2(cplexsum*nrad,lm_size,nsums*(option/2))
3390 type(pawang_type),intent(in) :: pawang
3391 3392 !Local variables-------------------------------
3393 !scalars
3394 integer :: ilm,ilm1,ilm2,ir,i1r,i2r,i3r,isel
3395 real(dp) :: fact,ro1i,ro1r,ro2i,ro2r
3396 character(len=500) :: msg
3397 !arrays
3398 3399 !************************************************************************
3400 3401 if(nsums/=1.and.nsums/=3) then
3402 msg='nsums must be 1 or 3!'
3403 MSG_BUG(msg)
3404 end if
3405 if(pawang%gnt_option==0) then
3406 msg='pawang%gnt_option=0!'
3407 MSG_BUG(msg)
3408 end if
3409 3410 if (option>=1) then
3411 3412 ! SUM1(r)= Sum_L{Rho1_L(r)*Rho2_L(r)} (L>0)
3413 ! --------------------------------------------------
3414 sum1=zero
3415 3416 ! ===== All input/output densities are REAL ====
3417 if (cplex1==1.and.cplex2==1.and.cplexsum==1) then
3418 ! One sum to compute
3419 if (nsums==1) then
3420 do ilm=2,lm_size
3421 if (lmselect1(ilm).and.lmselect2(ilm)) then
3422 sum1(:,1)=sum1(:,1)+rho1(:,ilm)*rho2(:,ilm)
3423 end if
3424 end do
3425 ! Three sums to compute
3426 else
3427 do ilm=2,lm_size
3428 if (lmselect1(ilm)) then
3429 sum1(:,1)=sum1(:,1)+rho1(:,ilm)**2
3430 if (lmselect2(ilm)) sum1(:,2)=sum1(:,2)+rho1(:,ilm)*rho2(:,ilm)
3431 end if
3432 if (lmselect2(ilm)) sum1(:,3)=sum1(:,3)+rho2(:,ilm)**2
3433 end do
3434 end if
3435 3436 ! ===== At least one of Rho1 and Rho2 is COMPLEX ====
3437 else
3438 ! One sum to compute
3439 if (nsums==1) then
3440 do ilm=2,lm_size
3441 if (lmselect1(ilm).and.lmselect2(ilm)) then
3442 do ir=1,nrad
3443 i1r=cplex1*(ir-1)+1;i2r=cplex2*(ir-1)+1;i3r=cplexsum*(ir-1)+1
3444 ro1r=rho1(i1r,ilm);ro1i=zero;if (cplex1==2) ro1i=rho1(i1r+1,ilm)
3445 ro2r=rho2(i2r,ilm);ro2i=zero;if (cplex2==2) ro2i=rho2(i2r+1,ilm)
3446 sum1(i3r,1)=sum1(i3r,1)+ro1r*ro2r-ro1i*ro2i
3447 if (cplexsum==2) sum1(i3r+1,1)=sum1(i3r+1,1)+ro1r*ro2i+ro1i*ro2r
3448 end do
3449 end if
3450 end do
3451 ! Three sums to compute
3452 else
3453 do ilm=2,lm_size
3454 do ir=1,nrad
3455 i1r=cplex1*(ir-1)+1;i2r=cplex2*(ir-1)+1;i3r=cplexsum*(ir-1)+1
3456 ro1r=rho1(i1r,ilm);ro1i=zero;if (cplex1==2) ro1i=rho1(i1r+1,ilm)
3457 ro2r=rho2(i2r,ilm);ro2i=zero;if (cplex2==2) ro2i=rho2(i2r+1,ilm)
3458 if (lmselect1(ilm)) then
3459 sum1(i3r,1)=sum1(i3r,1)+ro1r**2-ro1i**2
3460 if (lmselect2(ilm)) sum1(i3r,2)=sum1(i3r,2)+ro1r*ro2r-ro1i*ro2i
3461 end if
3462 if (lmselect2(ilm)) sum1(i3r,3)=sum1(i3r,3)+ro2r**2-ro2i**2
3463 if (cplexsum==2) then
3464 if (lmselect1(ilm)) then
3465 sum1(i3r+1,1)=sum1(i3r+1,1)+two*ro1r*ro1i
3466 if (lmselect2(ilm)) sum1(i3r+1,2)=sum1(i3r+1,2)+ro1r*ro2i+ro1i*ro2r
3467 end if
3468 if (lmselect2(ilm)) sum1(i3r+1,3)=sum1(i3r+1,3)+two*ro2r*ro2i
3469 end if
3470 end do
3471 end do
3472 end if ! nsums
3473 end if ! cplex
3474 3475 end if !option
3476 3477 if (option>=2) then
3478 3479 ! SUM2(r,L)= Sum_L1_L2{Rho1_L1(r)*Rho2_L2(r)*Gaunt_(L,L1,L2)} (L1>0, L2>0)
3480 ! --------------------------------------------------
3481 sum2=zero
3482 ! ===== All input/output densities are REAL ====
3483 if (cplex1==1.and.cplex2==1.and.cplexsum==1) then
3484 ! One sum to compute
3485 if (nsums==1) then
3486 do ilm=1,lm_size
3487 do ilm1=2,lm_size
3488 if (lmselect1(ilm1)) then
3489 do ilm2=2,ilm1
3490 if (lmselect2(ilm2)) then
3491 isel=pawang%gntselect(ilm,ilm2+ilm1*(ilm1-1)/2)
3492 if (isel>0) then
3493 fact=pawang%realgnt(isel);if (ilm1/=ilm2) fact=two*fact
3494 sum2(:,ilm,1)=sum2(:,ilm,1)+fact*rho1(:,ilm1)*rho2(:,ilm2)
3495 end if
3496 end if
3497 end do
3498 end if
3499 end do
3500 end do
3501 ! Three sums to compute
3502 else
3503 do ilm=1,lm_size
3504 do ilm1=2,lm_size
3505 if (lmselect1(ilm1)) then
3506 do ilm2=2,ilm1
3507 if (lmselect1(ilm2)) then
3508 isel=pawang%gntselect(ilm,ilm2+ilm1*(ilm1-1)/2)
3509 if (isel>0) then
3510 fact=pawang%realgnt(isel);if (ilm1/=ilm2) fact=two*fact
3511 sum2(:,ilm,1)=sum2(:,ilm,1)+fact*rho1(:,ilm1)*rho1(:,ilm2)
3512 end if
3513 end if
3514 end do
3515 end if
3516 end do
3517 do ilm1=2,lm_size
3518 if (lmselect2(ilm1)) then
3519 do ilm2=2,ilm1
3520 if (lmselect2(ilm2)) then
3521 isel=pawang%gntselect(ilm,ilm2+ilm1*(ilm1-1)/2)
3522 if (isel>0) then
3523 fact=pawang%realgnt(isel);if (ilm1/=ilm2) fact=two*fact
3524 sum2(:,ilm,3)=sum2(:,ilm,3)+fact*rho2(:,ilm1)*rho2(:,ilm2)
3525 end if
3526 end if
3527 end do
3528 end if
3529 end do
3530 do ilm1=2,lm_size
3531 if (lmselect1(ilm1)) then
3532 do ilm2=2,ilm1
3533 if (lmselect2(ilm2)) then
3534 isel=pawang%gntselect(ilm,ilm2+ilm1*(ilm1-1)/2)
3535 if (isel>0) then
3536 fact=pawang%realgnt(isel)
3537 sum2(:,ilm,2)=sum2(:,ilm,2)+fact*rho1(:,ilm1)*rho2(:,ilm2)
3538 end if
3539 end if
3540 end do
3541 if (ilm1<lm_size) then
3542 do ilm2=ilm1+1,lm_size
3543 if (lmselect2(ilm2)) then
3544 isel=pawang%gntselect(ilm,ilm1+ilm2*(ilm2-1)/2)
3545 if (isel>0) then
3546 fact=pawang%realgnt(isel)
3547 sum2(:,ilm,2)=sum2(:,ilm,2)+fact*rho1(:,ilm1)*rho2(:,ilm2)
3548 end if
3549 end if
3550 end do
3551 end if
3552 end if
3553 end do
3554 end do
3555 end if ! nsums
3556 3557 ! ===== At least one of Rho1 and Rho2 is COMPLEX ====
3558 else
3559 ! One sum to compute
3560 if (nsums==1) then
3561 do ilm=1,lm_size
3562 do ilm1=2,lm_size
3563 if (lmselect1(ilm1)) then
3564 do ilm2=2,ilm1
3565 if (lmselect2(ilm2)) then
3566 isel=pawang%gntselect(ilm,ilm2+ilm1*(ilm1-1)/2)
3567 if (isel>0) then
3568 fact=pawang%realgnt(isel);if (ilm1/=ilm2) fact=two*fact
3569 do ir=1,nrad
3570 i1r=cplex1*(ir-1)+1;i2r=cplex2*(ir-1)+1;i3r=cplexsum*(ir-1)+1
3571 ro1r=rho1(i1r,ilm1);ro1i=zero;if (cplex1==2) ro1i=rho1(i1r+1,ilm1)
3572 ro2r=rho2(i2r,ilm2);ro2i=zero;if (cplex2==2) ro2i=rho2(i2r+1,ilm2)
3573 sum2(i3r,ilm,1)=sum2(i3r,ilm,1)+fact*(ro1r*ro2r-ro1i*ro2i)
3574 if (cplexsum==2) sum2(i3r+1,ilm,1)=sum2(i3r+1,ilm,1)+fact*(ro1r*ro2i+ro1i*ro2r)
3575 end do
3576 end if
3577 end if
3578 end do
3579 end if
3580 end do
3581 end do
3582 ! Three sums to compute
3583 else
3584 do ilm=2,lm_size
3585 do ir=1,nrad
3586 i1r=cplex1*(ir-1)+1;i2r=cplex2*(ir-1)+1;i3r=cplexsum*(ir-1)+1
3587 do ilm1=2,lm_size
3588 if (lmselect1(ilm1)) then
3589 ro1r=rho1(i1r,ilm1);ro1i=zero;if (cplex1==2) ro1i=rho1(i1r+1,ilm1)
3590 do ilm2=2,ilm1
3591 if (lmselect1(ilm2)) then
3592 isel=pawang%gntselect(ilm,ilm2+ilm1*(ilm1-1)/2)
3593 if (isel>0) then
3594 fact=pawang%realgnt(isel);if (ilm1/=ilm2) fact=two*fact
3595 ro2r=rho1(i1r,ilm2);ro2i=zero;if (cplex1==2) ro2i=rho1(i1r+1,ilm2)
3596 sum2(i3r,ilm,1)=sum2(i3r,ilm,1)+fact*(ro1r*ro2r-ro1i*ro2i)
3597 if (cplexsum==2) sum2(i3r+1,ilm,1)=sum2(i3r+1,ilm,1)+fact*(ro1r*ro2i+ro1i*ro2r)
3598 end if
3599 end if
3600 end do
3601 end if
3602 end do
3603 do ilm1=2,lm_size
3604 if (lmselect2(ilm1)) then
3605 ro1r=rho2(i2r,ilm1);ro1i=zero;if (cplex2==2) ro1i=rho2(i2r+1,ilm1)
3606 do ilm2=2,ilm1
3607 if (lmselect2(ilm2)) then
3608 isel=pawang%gntselect(ilm,ilm2+ilm1*(ilm1-1)/2)
3609 if (isel>0) then
3610 fact=pawang%realgnt(isel);if (ilm1/=ilm2) fact=two*fact
3611 ro2r=rho2(i2r,ilm2);ro2i=zero;if (cplex2==2) ro2i=rho2(i2r+1,ilm2)
3612 sum2(i3r,ilm,3)=sum2(i3r,ilm,3)+fact*(ro1r*ro2r-ro1i*ro2i)
3613 if (cplexsum==2) sum2(i3r+1,ilm,3)=sum2(i3r+1,ilm,3)+fact*(ro1r*ro2i+ro1i*ro2r)
3614 end if
3615 end if
3616 end do
3617 end if
3618 end do
3619 do ilm1=2,lm_size
3620 if (lmselect1(ilm1)) then
3621 ro1r=rho1(i1r,ilm1);ro1i=zero;if (cplex1==2) ro1i=rho1(i1r+1,ilm1)
3622 do ilm2=2,ilm1
3623 if (lmselect2(ilm2)) then
3624 isel=pawang%gntselect(ilm,ilm2+ilm1*(ilm1-1)/2)
3625 if (isel>0) then
3626 fact=pawang%realgnt(isel)
3627 ro2r=rho2(i2r,ilm2);ro2i=zero;if (cplex2==2) ro2i=rho2(i2r+1,ilm2)
3628 sum2(i3r,ilm,2)=sum2(i3r,ilm,2)+fact*(ro1r*ro2r-ro1i*ro2i)
3629 if (cplexsum==2) sum2(i3r+1,ilm,2)=sum2(i3r+1,ilm,2)+fact*(ro1r*ro2i+ro1i*ro2r)
3630 end if
3631 end if
3632 end do
3633 if (ilm1<lm_size) then
3634 do ilm2=ilm1+1,lm_size
3635 if (lmselect2(ilm2)) then
3636 isel=pawang%gntselect(ilm,ilm1+ilm2*(ilm2-1)/2)
3637 if (isel>0) then
3638 fact=pawang%realgnt(isel)
3639 ro2r=rho2(i2r,ilm2);ro2i=zero;if (cplex2==2) ro2i=rho2(i2r+1,ilm2)
3640 sum2(i3r,ilm,2)=sum2(i3r,ilm,2)+fact*(ro1r*ro2r-ro1i*ro2i)
3641 if (cplexsum==2) sum2(i3r+1,ilm,2)=sum2(i3r+1,ilm,2)+fact*(ro1r*ro2i+ro1i*ro2r)
3642 end if
3643 end if
3644 end do
3645 end if
3646 end if
3647 end do
3648 end do
3649 end do
3650 end if ! nsums
3651 3652 end if ! cplex
3653 3654 end if !option
3655 3656 end subroutine pawxcsum