TABLE OF CONTENTS
- ABINIT/m_cgtools
- m_blas/cgpaw_cholesky
- m_cgtools/cg_addtorho
- m_cgtools/cg_box2gsph
- m_cgtools/cg_dznrm2
- m_cgtools/cg_envlop
- m_cgtools/cg_filter
- m_cgtools/cg_from_reim
- m_cgtools/cg_fromcplx
- m_cgtools/cg_getspin
- m_cgtools/cg_gsph2box
- m_cgtools/cg_normev
- m_cgtools/cg_precon
- m_cgtools/cg_precon_block
- m_cgtools/cg_real_zdotc
- m_cgtools/cg_setaug_zero
- m_cgtools/cg_setval
- m_cgtools/cg_to_reim
- m_cgtools/cg_tocplx
- m_cgtools/cg_vlocpsi
- m_cgtools/cg_zaxpby
- m_cgtools/cg_zaxpy
- m_cgtools/cg_zcopy
- m_cgtools/cg_zdotc
- m_cgtools/cg_zdotu
- m_cgtools/cg_zgemm
- m_cgtools/cg_zgemv
- m_cgtools/cg_zscal
- m_cgtools/cgnc_cholesky
- m_cgtools/cgnc_gramschmidt
- m_cgtools/cgnc_gsortho
- m_cgtools/cgnc_normalize
- m_cgtools/cgpaw_gramschmidt
- m_cgtools/cgpaw_gsortho
- m_cgtools/cgpaw_normalize
- m_cgtools/cz_zprecon_block
- m_cgtools/dotprod_g
- m_cgtools/dotprod_v
- m_cgtools/dotprod_vn
- m_cgtools/fxphas_seq
- m_cgtools/matrixelmt_g
- m_cgtools/mean_fftr
- m_cgtools/overlap_g
- m_cgtools/projbd
- m_cgtools/set_istwfk
- m_cgtools/sqnorm_g
- m_cgtools/sqnorm_v
ABINIT/m_cgtools [ Modules ]
NAME
m_cgtools
FUNCTION
This module defines wrappers for BLAS routines. The arguments are stored using the "cg" convention, namely real array of shape cg(2,...)
COPYRIGHT
Copyright (C) 1992-2018 ABINIT group (MG, MT, XG, DCA, GZ, FB, MVer) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
NOTES
1) The convention about names of interfaced routine is: cg_<name>,
where <name> is equal to the name of the standard BLAS routine
2) Blas routines are called without an explicit interface on purpose since
a) The compiler should pass the base address of the array to the F77 BLAS
b) Any compiler would complain about type mismatch (REAL,COMPLEX)
if an explicit interface is given.
m_blas/cgpaw_cholesky [ Functions ]
[ Top ] [ m_blas ] [ Functions ]
NAME
cgpaw_cholesky
FUNCTION
Cholesky orthonormalization of the vectors stored in cgblock. (version for PAW wavefunctions).
INPUTS
npws=Size of each vector (usually npw*nspinor) nband=Number of band in cgblock and gsc istwfk=Storage mode for the wavefunctions. 1 for standard full mode me_g0=1 if this node has G=0. comm_pw=MPI communicator for the planewave group. Set to xmpi_comm_self for sequential mode.
SIDE EFFECTS
cgblock(2*npws*nband)
input: Input set of vectors |C>, S|C>
output: Orthonormalized set such as <C|S|C> = 1
gsc(2*npws*nband)
destroyed in output.
PARENTS
pw_orthon
CHILDREN
SOURCE
2836 subroutine cgpaw_cholesky(npws,nband,cgblock,gsc,istwfk,me_g0,comm_pw) 2837 2838 2839 !This section has been created automatically by the script Abilint (TD). 2840 !Do not modify the following lines by hand. 2841 #undef ABI_FUNC 2842 #define ABI_FUNC 'cgpaw_cholesky' 2843 !End of the abilint section 2844 2845 implicit none 2846 2847 !Arguments ------------------------------------ 2848 !scalars 2849 integer,intent(in) :: npws,nband,istwfk 2850 integer,intent(in) :: me_g0,comm_pw 2851 !arrays 2852 real(dp),intent(inout) :: cgblock(2*npws*nband) 2853 real(dp),intent(inout) :: gsc(2*npws*nband) 2854 2855 !Local variables ------------------------------ 2856 !scalars 2857 integer :: ierr,b1,b2 2858 character(len=500) :: msg 2859 !arrays 2860 real(dp),allocatable :: rovlp(:,:) 2861 real(dp) :: rcg0(nband),rg0sc(nband) 2862 complex(dpc),allocatable :: cf_ovlp(:,:) 2863 2864 ! ************************************************************************* 2865 2866 ! 1) Calculate O_ij = <phi_i|S|phi_j> 2867 ABI_MALLOC(cf_ovlp,(nband,nband)) 2868 2869 call ABI_ZGEMM("C","N",nband,nband,npws,cone,cgblock,npws,gsc,npws,czero,cf_ovlp,nband) 2870 2871 if (istwfk==1) then 2872 ! 2873 ! Sum the overlap if PW are distributed. 2874 if (comm_pw /= xmpi_comm_self) then 2875 call xmpi_sum(cf_ovlp,comm_pw,ierr) 2876 end if 2877 ! 2878 ! 2) Cholesky factorization: O = U^H U with U upper triangle matrix. 2879 call ZPOTRF('U',nband,cf_ovlp,nband,ierr) 2880 2881 if (ierr/=0) then 2882 write(msg,'(a,i0)')' ZPOTRF returned info= ',ierr 2883 MSG_ERROR(msg) 2884 end if 2885 2886 else 2887 ! overlap is real. Note that nspinor is always 1 in this case. 2888 ABI_MALLOC(rovlp,(nband,nband)) 2889 rovlp = two * REAL(cf_ovlp) 2890 2891 if (istwfk==2 .and. me_g0==1) then 2892 ! Extract the real part at G=0 and subtract its contribution to the overlap. 2893 call dcopy(nband,cgblock,2*npws,rcg0,1) 2894 call dcopy(nband,gsc,2*npws,rg0sc,1) 2895 do b2=1,nband 2896 do b1=1,b2 2897 rovlp(b1,b2) = rovlp(b1,b2) - rcg0(b1)*rg0sc(b2) 2898 end do 2899 end do 2900 end if 2901 ! 2902 ! Sum the overlap if PW are distributed. 2903 if (comm_pw /= xmpi_comm_self) then 2904 call xmpi_sum(rovlp,comm_pw,ierr) 2905 end if 2906 ! 2907 ! 2) Cholesky factorization: O = U^H U with U upper triangle matrix. 2908 call DPOTRF('U',nband,rovlp,nband,ierr) 2909 2910 if (ierr/=0) then 2911 write(msg,'(a,i0)')' DPOTRF returned info= ',ierr 2912 MSG_ERROR(msg) 2913 end if 2914 2915 cf_ovlp = DCMPLX(rovlp) 2916 ABI_FREE(rovlp) 2917 end if 2918 ! 2919 ! 3) Solve X U = cgblock. 2920 call ZTRSM('Right','Upper','Normal','Normal',npws,nband,cone,cf_ovlp,nband,cgblock,npws) 2921 2922 ! 4) Solve Y U = gsc. On exit <cgblock|gsc> = 1 2923 call ZTRSM('Right','Upper','Normal','Normal',npws,nband,cone,cf_ovlp,nband,gsc,npws) 2924 2925 ABI_FREE(cf_ovlp) 2926 2927 end subroutine cgpaw_cholesky
m_cgtools/cg_addtorho [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_addtorho
FUNCTION
Add |ur|**2 to the ground-states density rho.
rho = rho + weight_r * Re[ur]**2 + weight_i * Im[ur]**2
INPUTS
nx,ny,nz=physical dimension of the FFT box. ldx,ldy,ldz=leading dimensions of the arrays. ndat=number of contributions to accumulate. weight_r=weight used for the accumulation of the density in real space weight_i=weight used for the accumulation of the density in real space ur(2,ldx,ldy,ldz*ndat)=wavefunctions in real space
SIDE EFFECTS
rho(ldx,ldy,ldz) = contains the input density at input,
modified in input with the contribution gived by ur.
PARENTS
fourwf,m_dfti,m_fftw3
CHILDREN
SOURCE
2466 subroutine cg_addtorho(nx,ny,nz,ldx,ldy,ldz,ndat,weight_r,weight_i,ur,rho) 2467 2468 2469 !This section has been created automatically by the script Abilint (TD). 2470 !Do not modify the following lines by hand. 2471 #undef ABI_FUNC 2472 #define ABI_FUNC 'cg_addtorho' 2473 !End of the abilint section 2474 2475 implicit none 2476 2477 !Arguments ------------------------------------ 2478 !scalars 2479 integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,ndat 2480 real(dp),intent(in) :: weight_i,weight_r 2481 !arrays 2482 real(dp),intent(in) :: ur(2,ldx,ldy,ldz*ndat) 2483 real(dp),intent(inout) :: rho(ldx,ldy,ldz) 2484 2485 !Local variables------------------------------- 2486 !scalars 2487 integer :: ix,iy,iz,idat,izdat 2488 2489 ! ************************************************************************* 2490 2491 if (ndat==1) then 2492 !$OMP PARALLEL DO 2493 do iz=1,nz 2494 do iy=1,ny 2495 do ix=1,nx 2496 rho(ix,iy,iz) = rho(ix,iy,iz) + weight_r * ur(1,ix,iy,iz)**2 & 2497 & + weight_i * ur(2,ix,iy,iz)**2 2498 end do 2499 end do 2500 end do 2501 2502 else 2503 ! It would be nice to use $OMP PARALLEL DO PRIVATE(izdat) REDUCTION(+:rho) 2504 ! but it's risky as the private rho is allocated on the stack of the thread. 2505 !$OMP PARALLEL PRIVATE(izdat) 2506 do idat=1,ndat 2507 !$OMP DO 2508 do iz=1,nz 2509 izdat = iz + (idat-1)*ldz 2510 do iy=1,ny 2511 do ix=1,nx 2512 rho(ix,iy,iz) = rho(ix,iy,iz) + weight_r * ur(1,ix,iy,izdat)**2 & 2513 & + weight_i * ur(2,ix,iy,izdat)**2 2514 end do 2515 end do 2516 end do 2517 !$OMP END DO NOWAIT 2518 end do 2519 !$OMP END PARALLEL 2520 end if 2521 2522 end subroutine cg_addtorho
m_cgtools/cg_box2gsph [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_box2gsph
FUNCTION
INPUTS
nx,ny,nz=physical dimension of the FFT box. ldx,ldy,ldz=Logical dimensions of the arrays. ndat=number of data in iarrbox npw_k=Number of planewaves in the G-sphere. kg_k(3,npw_k)=Reduced coordinates of the G-vectoes. iarrbox(2,ldx,ldy,ldz*ndat)=Input arrays on the FFT box. [rscal] = Scaling factor
OUTPUT
oarrsph(2,npw_k*ndat)=Data defined on the G-sphere.
PARENTS
fourwf,m_dfti,m_fftw3
CHILDREN
SOURCE
2352 subroutine cg_box2gsph(nx,ny,nz,ldx,ldy,ldz,ndat,npw_k,kg_k,iarrbox,oarrsph,rscal) 2353 2354 2355 !This section has been created automatically by the script Abilint (TD). 2356 !Do not modify the following lines by hand. 2357 #undef ABI_FUNC 2358 #define ABI_FUNC 'cg_box2gsph' 2359 !End of the abilint section 2360 2361 implicit none 2362 2363 !Arguments ------------------------------------ 2364 !scalars 2365 integer,intent(in) :: npw_k,nx,ny,nz,ldx,ldy,ldz,ndat 2366 real(dp),optional,intent(in) :: rscal 2367 !arrays 2368 integer,intent(in) :: kg_k(3,npw_k) 2369 real(dp),intent(in) :: iarrbox(2,ldx*ldy*ldz*ndat) 2370 real(dp),intent(out) :: oarrsph(2,npw_k*ndat) 2371 2372 !Local variables------------------------------- 2373 !scalars 2374 integer :: ig,ix,iy,iz,idat,sph_pad,box_pad,ifft 2375 2376 ! ************************************************************************* 2377 2378 if (.not. PRESENT(rscal)) then 2379 ! 2380 if (ndat==1) then 2381 !$OMP PARALLEL DO PRIVATE(ix,iy,iz,ifft) 2382 do ig=1,npw_k 2383 ix=kg_k(1,ig); if (ix<0) ix=ix+nx; ix=ix+1 2384 iy=kg_k(2,ig); if (iy<0) iy=iy+ny; iy=iy+1 2385 iz=kg_k(3,ig); if (iz<0) iz=iz+nz; iz=iz+1 2386 ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy 2387 oarrsph(1,ig) = iarrbox(1,ifft) 2388 oarrsph(2,ig) = iarrbox(2,ifft) 2389 end do 2390 else 2391 !$OMP PARALLEL DO PRIVATE(sph_pad,box_pad,ix,iy,iz,ifft) 2392 do idat=1,ndat 2393 sph_pad = (idat-1)*npw_k 2394 box_pad = (idat-1)*ldx*ldy*ldz 2395 do ig=1,npw_k 2396 ix=kg_k(1,ig); if (ix<0) ix=ix+nx; ix=ix+1 2397 iy=kg_k(2,ig); if (iy<0) iy=iy+ny; iy=iy+1 2398 iz=kg_k(3,ig); if (iz<0) iz=iz+nz; iz=iz+1 2399 ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy 2400 oarrsph(1,ig+sph_pad) = iarrbox(1,ifft+box_pad) 2401 oarrsph(2,ig+sph_pad) = iarrbox(2,ifft+box_pad) 2402 end do 2403 end do 2404 end if 2405 ! 2406 else 2407 if (ndat==1) then 2408 !$OMP PARALLEL DO PRIVATE(ix,iy,iz,ifft) 2409 do ig=1,npw_k 2410 ix=kg_k(1,ig); if (ix<0) ix=ix+nx; ix=ix+1 2411 iy=kg_k(2,ig); if (iy<0) iy=iy+ny; iy=iy+1 2412 iz=kg_k(3,ig); if (iz<0) iz=iz+nz; iz=iz+1 2413 ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy 2414 oarrsph(1,ig) = iarrbox(1,ifft) * rscal 2415 oarrsph(2,ig) = iarrbox(2,ifft) * rscal 2416 end do 2417 else 2418 !$OMP PARALLEL DO PRIVATE(sph_pad,box_pad,ix,iy,iz,ifft) 2419 do idat=1,ndat 2420 sph_pad = (idat-1)*npw_k 2421 box_pad = (idat-1)*ldx*ldy*ldz 2422 do ig=1,npw_k 2423 ix=kg_k(1,ig); if (ix<0) ix=ix+nx; ix=ix+1 2424 iy=kg_k(2,ig); if (iy<0) iy=iy+ny; iy=iy+1 2425 iz=kg_k(3,ig); if (iz<0) iz=iz+nz; iz=iz+1 2426 ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy 2427 oarrsph(1,ig+sph_pad) = iarrbox(1,ifft+box_pad) * rscal 2428 oarrsph(2,ig+sph_pad) = iarrbox(2,ifft+box_pad) * rscal 2429 end do 2430 end do 2431 end if 2432 end if 2433 2434 end subroutine cg_box2gsph
m_cgtools/cg_dznrm2 [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_dznrm2
FUNCTION
returns the euclidean norm of a vector via the function name, so that DZNRM2 := sqrt( x**H*x )
INPUTS
n = Specifies the number of elements in vector x. x(2*x) = Input array.
OUTPUT
OUTPUT
PARENTS
SOURCE
627 function cg_dznrm2(n, x) result(res) 628 629 630 !This section has been created automatically by the script Abilint (TD). 631 !Do not modify the following lines by hand. 632 #undef ABI_FUNC 633 #define ABI_FUNC 'cg_dznrm2' 634 !End of the abilint section 635 636 implicit none 637 638 !Arguments ------------------------------------ 639 !scalars 640 integer,intent(in) :: n 641 real(dp) :: res 642 !arrays 643 real(dp),intent(in) :: x(2*n) 644 real(dp),external :: dznrm2 645 646 ! ************************************************************************* 647 648 res = dznrm2(n, x, 1) 649 650 end function cg_dznrm2
m_cgtools/cg_envlop [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_envlop
FUNCTION
Multiply random number values in cg by envelope function to lower initial kinetic energy. Envelope $\left( 1-\left( G/G_{\max }\right) ^2\right) ^{power}$ for |G|<= Gmax. Near G=0, little scaling, and goes to zero flatly near Gmax.Loop over perturbations
INPUTS
cg(2,npw*nband)=initial random number wavefunctions ecut=kinetic energy cutoff in Ha gmet(3,3)=reciprocal space metric (bohr^-2) icgmod=shift to be given to the location of data in cg kg(3,npw)=reduced coordinates of G vectors in basis sphere kpoint(3)=reduced coordinates of k point mcg=maximum second dimension of cg (at least npw*nband*nspinor) nband=number of bands being considered npw=number of planewaves in basis sphere nspinor=number of spinorial components of the wavefunctions
OUTPUT
cg(2,mcg)=revised values (not orthonormalized)
PARENTS
wfconv
CHILDREN
SOURCE
3677 subroutine cg_envlop(cg,ecut,gmet,icgmod,kg,kpoint,mcg,nband,npw,nspinor) 3678 3679 3680 !This section has been created automatically by the script Abilint (TD). 3681 !Do not modify the following lines by hand. 3682 #undef ABI_FUNC 3683 #define ABI_FUNC 'cg_envlop' 3684 !End of the abilint section 3685 3686 implicit none 3687 3688 !Arguments ------------------------------------ 3689 !scalars 3690 integer,intent(in) :: icgmod,mcg,nband,npw,nspinor 3691 real(dp),intent(in) :: ecut 3692 !arrays 3693 integer,intent(in) :: kg(3,npw) 3694 real(dp),intent(in) :: gmet(3,3),kpoint(3) 3695 real(dp),intent(inout) :: cg(2,mcg) 3696 3697 !Local variables------------------------------- 3698 !scalars 3699 integer,parameter :: re=1,im=2,power=12 3700 integer :: i1,i2,i3,ig,igs,ispinor,nn,spad 3701 real(dp) :: cutoff,gs,kpgsqc 3702 !character(len=500) :: msg 3703 !arrays 3704 real(dp),allocatable :: cut_pws(:) 3705 3706 ! ************************************************************************* 3707 3708 !$(k+G)^2$ cutoff from $(1/2)(2 Pi (k+G))^2 = ecut$ 3709 kpgsqc=ecut/(2.0_dp*pi**2) 3710 cutoff = kpgsqc 3711 3712 ABI_MALLOC(cut_pws,(npw)) 3713 3714 !Run through G vectors in basis 3715 !$OMP PARALLEL DO PRIVATE(gs) 3716 do ig=1,npw 3717 i1=kg(1,ig) ; i2=kg(2,ig) ; i3=kg(3,ig) 3718 !(k+G)^2 evaluated using metric and kpoint 3719 gs = gmet(1,1)*(kpoint(1)+dble(i1))**2+& 3720 & gmet(2,2)*(kpoint(2)+dble(i2))**2+& 3721 & gmet(3,3)*(kpoint(3)+dble(i3))**2+& 3722 & 2.0_dp*(gmet(2,1)*(kpoint(2)+dble(i2))*(kpoint(1)+dble(i1))+& 3723 & gmet(3,2)*(kpoint(3)+dble(i3))*(kpoint(2)+dble(i2))+& 3724 & gmet(1,3)*(kpoint(1)+dble(i1))*(kpoint(3)+dble(i3))) 3725 if (gs>cutoff) then 3726 cut_pws(ig) = zero 3727 else 3728 cut_pws(ig) = (1.0_dp-(gs/cutoff))**power 3729 end if 3730 end do 3731 3732 !Run through bands (real and imaginary components) 3733 !$OMP PARALLEL DO PRIVATE(spad,igs) 3734 do nn=1,nband 3735 spad = (nn-1)*npw*nspinor+icgmod 3736 do ispinor=1,nspinor 3737 do ig=1,npw 3738 igs=ig+(ispinor-1)*npw 3739 cg(1,igs+spad) = cg(1,igs+spad) * cut_pws(ig) 3740 cg(2,igs+spad) = cg(2,igs+spad) * cut_pws(ig) 3741 end do 3742 end do 3743 end do 3744 3745 ABI_FREE(cut_pws) 3746 3747 end subroutine cg_envlop
m_cgtools/cg_filter [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_filter
FUNCTION
Set all the elements of x to zero where mask is .TRUE.
INPUTS
n=Specifies the number of elements in vectors x and y. mask(n)=Logical array.
SIDE EFFECTS
x(n)=See description.
PARENTS
SOURCE
311 pure subroutine cg_filter(n, x, mask) 312 313 314 !This section has been created automatically by the script Abilint (TD). 315 !Do not modify the following lines by hand. 316 #undef ABI_FUNC 317 #define ABI_FUNC 'cg_filter' 318 !End of the abilint section 319 320 implicit none 321 322 !Arguments ------------------------------------ 323 !scalars 324 integer,intent(in) :: n 325 !arrays 326 real(dp),intent(inout) :: x(2,n) 327 logical,intent(in) :: mask(n) 328 329 ! ************************************************************************* 330 331 where (mask) 332 x(1,:) = zero 333 x(2,:) = zero 334 end where 335 336 end subroutine cg_filter
m_cgtools/cg_from_reim [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_from_reim
FUNCTION
INPUTS
PARENTS
CHILDREN
SOURCE
472 subroutine cg_from_reim(npw,ndat,reim,factor,cg) 473 474 475 !This section has been created automatically by the script Abilint (TD). 476 !Do not modify the following lines by hand. 477 #undef ABI_FUNC 478 #define ABI_FUNC 'cg_from_reim' 479 !End of the abilint section 480 481 implicit none 482 483 !Arguments ------------------------------------ 484 !scalars 485 integer,intent(in) :: npw,ndat 486 real(dp),intent(in) :: factor 487 !arrays 488 real(dp),intent(in) :: reim(npw*ndat*2) 489 real(dp),intent(out) :: cg(2*npw*ndat) 490 491 ! ************************************************************************* 492 493 ! UnPack real and imaginary part and multiply by scale factor if /= one. 494 ! Could use blocking but oh well 495 call dcopy(npw*ndat, reim(1), 1, cg(1), 2) 496 call dcopy(npw*ndat, reim(npw*ndat+1), 1, cg(2), 2) 497 498 if (factor /= one) call dscal(2*npw*ndat,factor, cg(1), 1) 499 500 end subroutine cg_from_reim
m_cgtools/cg_fromcplx [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_fromcplx
FUNCTION
Convert a complex array to a real array with (real,imag) part
INPUTS
n = Specifies the number of elements in icplx and ocg. icplx(n)=Input complex array.
OUTPUT
ocg(2*n)=Output array with real and imaginary part.
PARENTS
CHILDREN
SOURCE
257 subroutine cg_fromcplx(n,icplx,ocg) 258 259 260 !This section has been created automatically by the script Abilint (TD). 261 !Do not modify the following lines by hand. 262 #undef ABI_FUNC 263 #define ABI_FUNC 'cg_fromcplx' 264 !End of the abilint section 265 266 implicit none 267 268 !Arguments ------------------------------------ 269 !scalars 270 integer,intent(in) :: n 271 !arrays 272 real(dp),intent(out) :: ocg(2*n) 273 complex(dpc),intent(in) :: icplx(n) 274 275 !Local variables ------------------------------ 276 !scalars 277 integer :: ii,idx 278 279 ! ************************************************************************* 280 281 !$OMP PARALLEL DO PRIVATE(ii,idx) 282 do ii=1,n 283 idx = 2*ii-1 284 ocg(idx ) = DBLE (icplx(ii)) 285 ocg(idx+1) = AIMAG(icplx(ii)) 286 end do 287 288 end subroutine cg_fromcplx
m_cgtools/cg_getspin [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_getspin
FUNCTION
Sandwich a single wave function on the Pauli matrices
INPUTS
npw_k = number of plane waves cgcband = coefficients of spinorial wave function
OUTPUT
spin = 3-vector of spin components for this state cgcmat = outer spin product of spinorial wf with itself
PARENTS
m_cut3d,partial_dos_fractions
CHILDREN
SOURCE
2099 subroutine cg_getspin(cgcband, npw_k, spin, cgcmat) 2100 2101 2102 !This section has been created automatically by the script Abilint (TD). 2103 !Do not modify the following lines by hand. 2104 #undef ABI_FUNC 2105 #define ABI_FUNC 'cg_getspin' 2106 !End of the abilint section 2107 2108 implicit none 2109 2110 !Arguments ------------------------------------ 2111 !scalars 2112 integer, intent(in) :: npw_k 2113 real(dp), intent(in) :: cgcband(2,2*npw_k) 2114 complex(dpc), intent(out),optional :: cgcmat(2,2) 2115 real(dp), intent(out) :: spin(3) 2116 2117 !Local variables------------------------------- 2118 !scalars 2119 complex(dpc),parameter :: pauli_0(2,2) = reshape([cone,czero,czero,cone], [2,2]) 2120 complex(dpc),parameter :: pauli_x(2,2) = reshape([czero,cone,cone,czero], [2,2]) 2121 complex(dpc),parameter :: pauli_y(2,2) = reshape([czero,j_dpc,-j_dpc,czero], [2,2]) 2122 complex(dpc),parameter :: pauli_z(2,2) = reshape([cone,czero,czero,-cone], [2,2]) 2123 complex(dpc) :: cspin(0:3), cgcmat_(2,2) 2124 ! *********************************************************************** 2125 2126 ! cgcmat_ = cgcband * cgcband^T* i.e. 2x2 matrix of spin components (dpcomplex) 2127 cgcmat_ = czero 2128 call zgemm('n','c',2,2,npw_k,cone,cgcband,2,cgcband,2,czero,cgcmat_,2) 2129 2130 ! spin(*) = sum_{si sj pi} cgcband(si,pi)^* pauli_*(si,sj) cgcband(sj,pi) 2131 cspin(0) = cgcmat_(1,1)*pauli_0(1,1) + cgcmat_(2,1)*pauli_0(2,1) & 2132 & + cgcmat_(1,2)*pauli_0(1,2) + cgcmat_(2,2)*pauli_0(2,2) 2133 cspin(1) = cgcmat_(1,1)*pauli_x(1,1) + cgcmat_(2,1)*pauli_x(2,1) & 2134 & + cgcmat_(1,2)*pauli_x(1,2) + cgcmat_(2,2)*pauli_x(2,2) 2135 cspin(2) = cgcmat_(1,1)*pauli_y(1,1) + cgcmat_(2,1)*pauli_y(2,1) & 2136 & + cgcmat_(1,2)*pauli_y(1,2) + cgcmat_(2,2)*pauli_y(2,2) 2137 cspin(3) = cgcmat_(1,1)*pauli_z(1,1) + cgcmat_(2,1)*pauli_z(2,1) & 2138 & + cgcmat_(1,2)*pauli_z(1,2) + cgcmat_(2,2)*pauli_z(2,2) 2139 !write(std_out,*) 'cgmat: ', cgcmat_ 2140 !write(std_out,*) 'real(spin): ', real(cspin) 2141 !write(std_out,*) 'aimag(spin): ', aimag(cspin) 2142 2143 spin = real(cspin(1:3)) 2144 if (present(cgcmat)) cgcmat = cgcmat_ 2145 2146 end subroutine cg_getspin
m_cgtools/cg_gsph2box [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_gsph2box
FUNCTION
Array iarrsph is defined in sphere with npw_k points. Insert iarrsph inside box of nx*ny*nz points to define array oarrbox for fft box. rest of oarrbox is filled with 0 s.
INPUTS
iarrsph(2,npw_k*ndat)= contains values for npw_k G vectors in basis sphere ndat=number of FFT to perform. npw_k=number of G vectors in basis at this k point oarrbox(2,ldx*ldy*ldz*ndat) = fft box nx,ny,nz=physical dimension of the box (oarrbox) ldx,ldy,ldz=memory dimension of oarrbox kg_k(3,npw_k)=integer coordinates of G vectors in basis sphere istwf_k=option parameter that describes the storage of wfs
OUTPUT
oarrbox(ldx*ldy*ldz*ndat)
NOTES
If istwf_k differs from 1, then special storage modes must be taken into account, for symmetric wavefunctions coming from k=(0 0 0) or other special k points.
PARENTS
CHILDREN
SOURCE
2183 subroutine cg_gsph2box(nx,ny,nz,ldx,ldy,ldz,ndat,npw_k,istwf_k,kg_k,iarrsph,oarrbox) 2184 2185 2186 !This section has been created automatically by the script Abilint (TD). 2187 !Do not modify the following lines by hand. 2188 #undef ABI_FUNC 2189 #define ABI_FUNC 'cg_gsph2box' 2190 !End of the abilint section 2191 2192 implicit none 2193 2194 !Arguments ------------------------------------ 2195 !scalars 2196 integer,intent(in) :: istwf_k,nx,ny,nz,ldx,ldy,ldz,ndat,npw_k 2197 !arrays 2198 integer,intent(in) :: kg_k(3,npw_k) 2199 real(dp),intent(in) :: iarrsph(2,npw_k*ndat) 2200 real(dp),intent(out) :: oarrbox(2,ldx*ldy*ldz*ndat) 2201 2202 !Local variables------------------------------- 2203 !scalars 2204 integer,parameter :: me_g0=1 2205 integer :: ix,ixinv,iy,iyinv,iz,izinv,dat,ipw,npwmin,pad_box,pad_sph,ifft,ifft_inv,ldxyz 2206 character(len=500) :: msg 2207 !arrays 2208 integer,allocatable :: ixinver(:),iyinver(:),izinver(:) 2209 2210 ! ************************************************************************* 2211 2212 !In the case of special k-points, invariant under time-reversal, 2213 !but not Gamma, initialize the inverse coordinates 2214 !Remember indeed that 2215 !u_k(G) = u_{k+G0}(G-G0); u_{-k}(G) = u_k(G)^* 2216 !and therefore: 2217 !u_{G0/2}(G) = u_{G0/2}(-G-G0)^*. 2218 if (istwf_k>=2) then 2219 ABI_MALLOC(ixinver,(nx)) 2220 ABI_MALLOC(iyinver,(ny)) 2221 ABI_MALLOC(izinver,(nz)) 2222 if ( ANY(istwf_k==(/2,4,6,8/)) ) then 2223 ixinver(1)=1 2224 do ix=2,nx 2225 ixinver(ix)=nx+2-ix 2226 end do 2227 else 2228 do ix=1,nx 2229 ixinver(ix)=nx+1-ix 2230 end do 2231 end if 2232 if (istwf_k>=2 .and. istwf_k<=5) then 2233 iyinver(1)=1 2234 do iy=2,ny 2235 iyinver(iy)=ny+2-iy 2236 end do 2237 else 2238 do iy=1,ny 2239 iyinver(iy)=ny+1-iy 2240 end do 2241 end if 2242 if ( ANY(istwf_k==(/2,3,6,7/)) ) then 2243 izinver(1)=1 2244 do iz=2,nz 2245 izinver(iz)=nz+2-iz 2246 end do 2247 else 2248 do iz=1,nz 2249 izinver(iz)=nz+1-iz 2250 end do 2251 end if 2252 end if 2253 2254 ldxyz = ldx*ldy*ldz 2255 2256 if (istwf_k==1) then 2257 2258 !$OMP PARALLEL DO PRIVATE(pad_sph,pad_box,ix,iy,iz,ifft) 2259 do dat=1,ndat 2260 pad_sph = (dat-1)*npw_k 2261 pad_box = (dat-1)*ldxyz 2262 oarrbox(:,1+pad_box:ldxyz+pad_box) = zero ! zero the sub-array 2263 do ipw=1,npw_k 2264 ix=kg_k(1,ipw); if (ix<0) ix=ix+nx; ix=ix+1 2265 iy=kg_k(2,ipw); if (iy<0) iy=iy+ny; iy=iy+1 2266 iz=kg_k(3,ipw); if (iz<0) iz=iz+nz; iz=iz+1 2267 ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy 2268 #if defined __INTEL_COMPILER && defined HAVE_OPENMP 2269 if (ifft==0) then 2270 MSG_ERROR("prevent ifort+OMP from miscompiling this section on cronos") 2271 end if 2272 #endif 2273 oarrbox(1,ifft+pad_box) = iarrsph(1,ipw+pad_sph) 2274 oarrbox(2,ifft+pad_box) = iarrsph(2,ipw+pad_sph) 2275 end do 2276 end do 2277 2278 else if (istwf_k>=2) then 2279 ! 2280 npwmin=1 2281 if(istwf_k==2 .and. me_g0==1) then ! If gamma point, then oarrbox must be completed 2282 do dat=1,ndat 2283 pad_sph = (dat-1)*npw_k 2284 pad_box = (dat-1)*ldxyz 2285 oarrbox(1,1+pad_box) = iarrsph(1,1+pad_sph) 2286 oarrbox(2,1+pad_box) = zero 2287 end do 2288 npwmin=2 2289 end if 2290 2291 !$OMP PARALLEL DO PRIVATE(pad_sph,pad_box,ix,iy,iz,ixinv,iyinv,izinv,ifft,ifft_inv) 2292 do dat=1,ndat 2293 pad_sph = (dat-1)*npw_k 2294 pad_box = (dat-1)*ldxyz 2295 oarrbox(:,npwmin+pad_box:ldxyz+pad_box) = zero 2296 do ipw=npwmin,npw_k 2297 ix=kg_k(1,ipw); if(ix<0)ix=ix+nx; ix=ix+1 2298 iy=kg_k(2,ipw); if(iy<0)iy=iy+ny; iy=iy+1 2299 iz=kg_k(3,ipw); if(iz<0)iz=iz+nz; iz=iz+1 2300 ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy 2301 ! Construct the coordinates of -k-G 2302 ixinv=ixinver(ix); iyinv=iyinver(iy); izinv=izinver(iz) 2303 ifft_inv = ixinv + (iyinv-1)*ldx + (izinv-1)*ldx*ldy 2304 2305 oarrbox(:,ifft +pad_box) = iarrsph(:,ipw+pad_sph) 2306 oarrbox(1,ifft_inv+pad_box) = iarrsph(1,ipw+pad_sph) 2307 oarrbox(2,ifft_inv+pad_box) = -iarrsph(2,ipw+pad_sph) 2308 end do 2309 end do 2310 ! 2311 else 2312 write(msg,'(a,i0)')"Wrong istwfk ",istwf_k 2313 MSG_ERROR(msg) 2314 end if 2315 2316 if (istwf_k>=2) then 2317 ABI_FREE(ixinver) 2318 ABI_FREE(iyinver) 2319 ABI_FREE(izinver) 2320 end if 2321 2322 end subroutine cg_gsph2box
m_cgtools/cg_normev [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_normev
FUNCTION
Normalize a set of nband eigenvectors of complex length npw (real length 2*npw) and set phases to make cg(i,i) real and positive. Near convergence, cg(i,j) approaches delta(i,j).
INPUTS
cg(2*npw,nband)=unnormalized eigenvectors npw=dimension of cg as shown nband=number of eigenvectors and complex length thereof.
OUTPUT
cg(2*npw,nband)=nband normalized eigenvectors
PARENTS
subdiago
CHILDREN
SOURCE
3776 subroutine cg_normev(cg,npw,nband) 3777 3778 3779 !This section has been created automatically by the script Abilint (TD). 3780 !Do not modify the following lines by hand. 3781 #undef ABI_FUNC 3782 #define ABI_FUNC 'cg_normev' 3783 !End of the abilint section 3784 3785 implicit none 3786 3787 !Arguments ------------------------------------ 3788 !scalars 3789 integer,intent(in) :: npw,nband 3790 !arrays 3791 real(dp),intent(inout) :: cg(2*npw,nband) 3792 3793 !Local variables------------------------------- 3794 !scalars 3795 integer :: ii,jj 3796 real(dp) :: den,evim,evre,phim,phre,xnorm 3797 character(len=500) :: message 3798 3799 ! ************************************************************************* 3800 3801 !Loop over vectors 3802 do ii=1,nband 3803 ! find norm 3804 xnorm=0.0d0 3805 do jj=1,2*npw 3806 xnorm=xnorm+cg(jj,ii)**2 3807 end do 3808 3809 if((xnorm-one)**2>tol6)then 3810 write(message,'(6a,i6,a,es16.6,3a)' )ch10,& 3811 & 'normev: ',ch10,& 3812 & 'Starting xnorm should be close to one (tol is tol6).',ch10,& 3813 & 'However, for state number',ii,', xnorm=',xnorm,ch10,& 3814 & 'It might be that your LAPACK library has not been correctly installed.' 3815 MSG_BUG(message) 3816 end if 3817 3818 xnorm=1.0d0/sqrt(xnorm) 3819 ! Set up phase 3820 phre=cg(2*ii-1,ii) 3821 phim=cg(2*ii,ii) 3822 if (phim/=0.0d0) then 3823 den=1.0d0/sqrt(phre**2+phim**2) 3824 phre=phre*xnorm*den 3825 phim=phim*xnorm*den 3826 else 3827 ! give xnorm the same sign as phre (negate if negative) 3828 phre=sign(xnorm,phre) 3829 phim=0.0d0 3830 end if 3831 ! normalize with phase change 3832 do jj=1,2*npw,2 3833 evre=cg(jj,ii) 3834 evim=cg(jj+1,ii) 3835 cg(jj,ii)=phre*evre+phim*evim 3836 cg(jj+1,ii)=phre*evim-phim*evre 3837 end do 3838 end do 3839 3840 end subroutine cg_normev
m_cgtools/cg_precon [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_precon
FUNCTION
precondition $<G|(H-e_{n,k})|C_{n,k}>$
INPUTS
$cg(2,npw)=<G|C_{n,k}>$.
$eval=current band eigenvalue=<C_{n,k}|H|C_{n,k}>$.
istwf_k=option parameter that describes the storage of wfs
kinpw(npw)=(modified) kinetic energy for each plane wave (Hartree)
nspinor=number of spinorial components of the wavefunctions
$vect(2,npw)=<G|H|C_{n,k}>$.
npw=number of planewaves at this k point.
optekin= 1 if the kinetic energy used in preconditionning is modified
according to Kresse, Furthmuller, PRB 54, 11169 (1996)
0 otherwise
mg_g0=1 if the node treats G0.
comm=MPI communicator
OUTPUT
pcon(npw)=preconditionning matrix
vect(2,npw*nspinor)=<G|(H-eval)|C_{n,k}>*(polynomial ratio)
PARENTS
cgwf,dfpt_cgwf
CHILDREN
SOURCE
3878 subroutine cg_precon(cg,eval,istwf_k,kinpw,npw,nspinor,me_g0,optekin,pcon,vect,comm) 3879 3880 3881 !This section has been created automatically by the script Abilint (TD). 3882 !Do not modify the following lines by hand. 3883 #undef ABI_FUNC 3884 #define ABI_FUNC 'cg_precon' 3885 use interfaces_18_timing 3886 !End of the abilint section 3887 3888 implicit none 3889 3890 !Arguments ------------------------------------ 3891 !scalars 3892 integer,intent(in) :: istwf_k,npw,nspinor,optekin,me_g0,comm 3893 real(dp),intent(in) :: eval 3894 !arrays 3895 real(dp),intent(in) :: cg(2,npw*nspinor),kinpw(npw) 3896 real(dp),intent(inout) :: pcon(npw),vect(2,npw*nspinor) 3897 3898 !Local variables------------------------------- 3899 !scalars 3900 integer :: ierr,ig,igs,ipw1,ispinor 3901 real(dp) :: ek0,ek0_inv,fac,poly,xx 3902 character(len=500) :: message 3903 !arrays 3904 real(dp) :: tsec(2) 3905 3906 ! ************************************************************************* 3907 3908 !Compute mean kinetic energy of band 3909 if(istwf_k==1)then 3910 ek0=zero 3911 do ispinor=1,nspinor 3912 igs=(ispinor-1)*npw 3913 ! !$OMP PARALLEL DO PRIVATE(ig) REDUCTION(+:ek0) SHARED(cg,igs,kinpw,npw) 3914 do ig=1+igs,npw+igs 3915 if(kinpw(ig-igs)<huge(0.0_dp)*1.d-11)then 3916 ek0=ek0+kinpw(ig-igs)*(cg(1,ig)**2+cg(2,ig)**2) 3917 end if 3918 end do 3919 end do 3920 3921 else if (istwf_k>=2)then 3922 if (istwf_k==2 .and. me_g0 == 1)then 3923 ek0=zero ; ipw1=2 3924 if(kinpw(1)<huge(0.0_dp)*1.d-11)ek0=0.5_dp*kinpw(1)*cg(1,1)**2 3925 else 3926 ek0=zero ; ipw1=1 3927 end if 3928 do ispinor=1,nspinor 3929 igs=(ispinor-1)*npw 3930 ! !$OMP PARALLEL DO PRIVATE(ig) REDUCTION(+:ek0) SHARED(cg,ipw1,kinpw,npw) 3931 do ig=ipw1+igs,npw+igs 3932 if(kinpw(ig)<huge(0.0_dp)*1.d-11)then 3933 ek0=ek0+kinpw(ig)*(cg(1,ig)**2+cg(2,ig)**2) 3934 end if 3935 end do 3936 end do 3937 ek0=2.0_dp*ek0 3938 end if 3939 3940 call timab(48,1,tsec) 3941 call xmpi_sum(ek0,comm,ierr) 3942 call timab(48,2,tsec) 3943 3944 if(ek0<1.0d-10)then 3945 write(message,'(3a)')& 3946 & 'The mean kinetic energy of a wavefunction vanishes.',ch10,& 3947 & 'It is reset to 0.1Ha.' 3948 MSG_WARNING(message) 3949 ek0=0.1_dp 3950 end if 3951 3952 if (optekin==1) then 3953 ek0_inv=2.0_dp/(3._dp*ek0) 3954 else 3955 ek0_inv=1.0_dp/ek0 3956 end if 3957 3958 !Carry out preconditioning 3959 do ispinor=1,nspinor 3960 igs=(ispinor-1)*npw 3961 !$OMP PARALLEL DO PRIVATE(fac,ig,poly,xx) SHARED(cg,ek0_inv,eval,kinpw,igs,npw,vect,pcon) 3962 do ig=1+igs,npw+igs 3963 if(kinpw(ig-igs)<huge(0.0_dp)*1.d-11)then 3964 xx=kinpw(ig-igs)*ek0_inv 3965 ! Teter polynomial ratio 3966 poly=27._dp+xx*(18._dp+xx*(12._dp+xx*8._dp)) 3967 fac=poly/(poly+16._dp*xx**4) 3968 if (optekin==1) fac=two*fac 3969 pcon(ig-igs)=fac 3970 vect(1,ig)=( vect(1,ig)-eval*cg(1,ig) )*fac 3971 vect(2,ig)=( vect(2,ig)-eval*cg(2,ig) )*fac 3972 else 3973 pcon(ig-igs)=zero 3974 vect(1,ig)=zero 3975 vect(2,ig)=zero 3976 end if 3977 end do 3978 end do 3979 3980 end subroutine cg_precon
m_cgtools/cg_precon_block [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_precon_block
FUNCTION
precondition $<G|(H-e_{n,k})|C_{n,k}>$ for a block of band (band-FFT parallelisation) in the case of real WFs (istwfk/=1)
INPUTS
blocksize= size of blocks of bands
$cg(vectsize,blocksize)=<G|C_{n,k}> for a block of bands$.
$eval(blocksize,blocksize)=current block of bands eigenvalues=<C_{n,k}|H|C_{n,k}>$.
$ghc(vectsize,blocksize)=<G|H|C_{n,k}> for a block of bands$.
iterationnumber=number of iterative minimizations in LOBPCG
kinpw(npw)=(modified) kinetic energy for each plane wave (Hartree)
nspinor=number of spinorial components of the wavefunctions (on current proc)
$vect(vectsize,blocksize)=<G|H|C_{n,k}> for a block of bands$.
npw=number of planewaves at this k point.
optekin= 1 if the kinetic energy used in preconditionning is modified
according to Kresse, Furthmuller, PRB 54, 11169 (1996)
0 otherwise
optpcon= 0 the TPA preconditionning matrix does not depend on band
1 the TPA preconditionning matrix (not modified)
2 the TPA preconditionning matrix is independant of iterationnumber
vectsize= size of vectors
mg_g0=1 if this node has Gamma, 0 otherwise.
OUTPUT
vect(2,npw)=<g|(h-eval)|c_{n,k}>*(polynomial ratio)
SIDE EFFECTS
pcon(npw,blocksize)=preconditionning matrix
input if optpcon=0,2 and iterationnumber/=1
output if optpcon=0,2 and iterationnumber==1
PARENTS
m_lobpcg
CHILDREN
SOURCE
4026 subroutine cg_precon_block(cg,eval,blocksize,iterationnumber,kinpw,& 4027 & npw,nspinor,me_g0,optekin,optpcon,pcon,ghc,vect,vectsize,comm) 4028 4029 4030 !This section has been created automatically by the script Abilint (TD). 4031 !Do not modify the following lines by hand. 4032 #undef ABI_FUNC 4033 #define ABI_FUNC 'cg_precon_block' 4034 use interfaces_18_timing 4035 !End of the abilint section 4036 4037 implicit none 4038 4039 !Arguments ------------------------------------ 4040 !scalars 4041 integer,intent(in) :: blocksize,iterationnumber,npw,nspinor,optekin 4042 integer,intent(in) :: optpcon,vectsize,me_g0,comm 4043 !arrays 4044 real(dp),intent(in) :: cg(vectsize,blocksize),eval(blocksize,blocksize) 4045 real(dp),intent(in) :: ghc(vectsize,blocksize),kinpw(npw) 4046 real(dp),intent(inout) :: pcon(npw,blocksize),vect(vectsize,blocksize) 4047 4048 !Local variables------------------------------- 4049 !scalars 4050 integer :: iblocksize,ierr,ig,igs,ipw1,ispinor 4051 real(dp) :: fac,poly,xx 4052 character(len=500) :: message 4053 !arrays 4054 real(dp) :: tsec(2) 4055 real(dp),allocatable :: ek0(:),ek0_inv(:) 4056 4057 ! ************************************************************************* 4058 4059 call timab(536,1,tsec) 4060 4061 !In this case, the Teter, Allan and Payne preconditioner is approximated: 4062 !the factor xx=Ekin(G) and no more Ekin(G)/Ekin(iband) 4063 if (optpcon==0) then 4064 do ispinor=1,nspinor 4065 igs=(ispinor-1)*npw 4066 if (me_g0 == 1) then 4067 do ig=1+igs,1+igs !g=0 4068 if (iterationnumber==1) then 4069 if(kinpw(ig-igs)<huge(0.0_dp)*1.d-11)then 4070 xx=kinpw(ig-igs) 4071 ! teter polynomial ratio 4072 poly=27._dp+xx*(18._dp+xx*(12._dp+xx*8._dp)) 4073 fac=poly/(poly+16._dp*xx**4) 4074 if (optekin==1) fac=two*fac 4075 pcon(ig-igs,1)=fac 4076 do iblocksize=1,blocksize 4077 vect(ig,iblocksize)=(ghc(ig,iblocksize)-& 4078 & eval(iblocksize,iblocksize)*cg(ig,iblocksize))*pcon(ig-igs,1) 4079 end do 4080 else 4081 pcon(ig-igs,1)=zero 4082 vect(ig,:)=0.0_dp 4083 end if 4084 else 4085 do iblocksize=1,blocksize 4086 vect(ig,iblocksize)=(ghc(ig,iblocksize)-& 4087 & eval(iblocksize,iblocksize)*cg(ig,iblocksize))*pcon(ig-igs,1) 4088 end do 4089 end if 4090 end do 4091 do ig=2+igs,npw+igs 4092 if (iterationnumber==1) then 4093 if(kinpw(ig-igs)<huge(0.0_dp)*1.d-11)then 4094 xx=kinpw(ig-igs) 4095 ! teter polynomial ratio 4096 poly=27._dp+xx*(18._dp+xx*(12._dp+xx*8._dp)) 4097 fac=poly/(poly+16._dp*xx**4) 4098 if (optekin==1) fac=two*fac 4099 pcon(ig-igs,1)=fac 4100 do iblocksize=1,blocksize 4101 vect(ig,iblocksize)=(ghc(ig,iblocksize)-& 4102 & eval(iblocksize,iblocksize)*cg(ig,iblocksize))*pcon(ig-igs,1) 4103 vect(ig+npw-1,iblocksize)=(ghc(ig+npw-1,iblocksize)-& 4104 & eval(iblocksize,iblocksize)*cg(ig+npw-1,iblocksize))*pcon(ig-igs,1) 4105 end do 4106 else 4107 pcon(ig-igs,1)=zero 4108 vect(ig,:)=zero 4109 vect(ig+npw-1,:)=zero 4110 end if 4111 else 4112 do iblocksize=1,blocksize 4113 vect(ig,iblocksize)=(ghc(ig,iblocksize)-& 4114 & eval(iblocksize,iblocksize)*cg(ig,iblocksize))*pcon(ig-igs,1) 4115 vect(ig+npw-1,iblocksize)=(ghc(ig+npw-1,iblocksize)-& 4116 & eval(iblocksize,iblocksize)*cg(ig+npw-1,iblocksize))*pcon(ig-igs,1) 4117 end do 4118 end if 4119 end do 4120 else 4121 do ig=1+igs,npw+igs 4122 if (iterationnumber==1) then 4123 if(kinpw(ig-igs)<huge(0.0_dp)*1.d-11)then 4124 xx=kinpw(ig-igs) 4125 ! teter polynomial ratio 4126 poly=27._dp+xx*(18._dp+xx*(12._dp+xx*8._dp)) 4127 fac=poly/(poly+16._dp*xx**4) 4128 if (optekin==1) fac=two*fac 4129 pcon(ig-igs,1)=fac 4130 do iblocksize=1,blocksize 4131 vect(ig,iblocksize)=(ghc(ig,iblocksize)-& 4132 & eval(iblocksize,iblocksize)*cg(ig,iblocksize))*pcon(ig-igs,1) 4133 vect(ig+npw,iblocksize)=(ghc(ig+npw,iblocksize)-& 4134 & eval(iblocksize,iblocksize)*cg(ig+npw,iblocksize))*pcon(ig-igs,1) 4135 end do 4136 else 4137 pcon(ig-igs,:)=zero 4138 vect(ig,:)=zero 4139 vect(ig+npw,:)=zero 4140 end if 4141 else 4142 do iblocksize=1,blocksize 4143 vect(ig,iblocksize)=(ghc(ig,iblocksize)-& 4144 & eval(iblocksize,iblocksize)*cg(ig,iblocksize))*pcon(ig-igs,1) 4145 vect(ig+npw,iblocksize)=(ghc(ig+npw,iblocksize)-& 4146 & eval(iblocksize,iblocksize)*cg(ig+npw,iblocksize))*pcon(ig-igs,1) 4147 end do 4148 end if 4149 end do 4150 end if 4151 end do 4152 4153 else if (optpcon>0) then 4154 ! Compute mean kinetic energy of all bands 4155 ABI_ALLOCATE(ek0,(blocksize)) 4156 ABI_ALLOCATE(ek0_inv,(blocksize)) 4157 if (iterationnumber==1.or.optpcon==1) then 4158 do iblocksize=1,blocksize 4159 if (me_g0 == 1)then 4160 ek0(iblocksize)=0.0_dp ; ipw1=2 4161 if(kinpw(1)<huge(0.0_dp)*1.d-11)ek0(iblocksize)=0.5_dp*kinpw(1)*cg(1,iblocksize)**2 4162 do ig=ipw1,npw 4163 if(kinpw(ig)<huge(0.0_dp)*1.d-11)then 4164 ek0(iblocksize)=ek0(iblocksize)+& 4165 & kinpw(ig)*(cg(ig,iblocksize)**2+cg(ig+npw-1,iblocksize)**2) 4166 end if 4167 end do 4168 else 4169 ek0(iblocksize)=0.0_dp ; ipw1=1 4170 do ig=ipw1,npw 4171 if(kinpw(ig)<huge(0.0_dp)*1.d-11)then 4172 ek0(iblocksize)=ek0(iblocksize)+& 4173 & kinpw(ig)*(cg(ig,iblocksize)**2+cg(ig+npw,iblocksize)**2) 4174 end if 4175 end do 4176 end if 4177 end do 4178 4179 call xmpi_sum(ek0,comm,ierr) 4180 4181 do iblocksize=1,blocksize 4182 if(ek0(iblocksize)<1.0d-10)then 4183 write(message, '(4a)' )ch10,& 4184 & 'cg_precon_block: the mean kinetic energy of a wavefunction vanishes.',ch10,& 4185 & 'it is reset to 0.1ha.' 4186 MSG_WARNING(message) 4187 ek0(iblocksize)=0.1_dp 4188 end if 4189 end do 4190 if (optekin==1) then 4191 ek0_inv(:)=2.0_dp/(3._dp*ek0(:)) 4192 else 4193 ek0_inv(:)=1.0_dp/ek0(:) 4194 end if 4195 end if !iterationnumber==1.or.optpcon==1 4196 4197 ! Carry out preconditioning 4198 do iblocksize=1,blocksize 4199 do ispinor=1,nspinor 4200 igs=(ispinor-1)*npw 4201 if (me_g0 == 1) then 4202 do ig=1+igs,1+igs !g=0 4203 if (iterationnumber==1.or.optpcon==1) then 4204 if(kinpw(ig-igs)<huge(0.0_dp)*1.d-11)then 4205 xx=kinpw(ig-igs)*ek0_inv(iblocksize) 4206 ! teter polynomial ratio 4207 poly=27._dp+xx*(18._dp+xx*(12._dp+xx*8._dp)) 4208 fac=poly/(poly+16._dp*xx**4) 4209 if (optekin==1) fac=two*fac 4210 pcon(ig-igs,iblocksize)=fac 4211 vect(ig,iblocksize)=(ghc(ig,iblocksize)-& 4212 & eval(iblocksize,iblocksize)*cg(ig,iblocksize))*fac 4213 else 4214 pcon(ig-igs,iblocksize)=zero 4215 vect(ig,iblocksize)=0.0_dp 4216 end if 4217 else 4218 vect(ig,iblocksize)=(ghc(ig,iblocksize)-& 4219 & eval(iblocksize,iblocksize)*cg(ig,iblocksize))*pcon(ig-igs,iblocksize) 4220 end if 4221 end do 4222 do ig=2+igs,npw+igs 4223 if (iterationnumber==1.or.optpcon==1) then 4224 if(kinpw(ig-igs)<huge(0.0_dp)*1.d-11)then 4225 xx=kinpw(ig-igs)*ek0_inv(iblocksize) 4226 ! teter polynomial ratio 4227 poly=27._dp+xx*(18._dp+xx*(12._dp+xx*8._dp)) 4228 fac=poly/(poly+16._dp*xx**4) 4229 if (optekin==1) fac=two*fac 4230 pcon(ig-igs,iblocksize)=fac 4231 vect(ig,iblocksize)=(ghc(ig,iblocksize)-& 4232 & eval(iblocksize,iblocksize)*cg(ig,iblocksize))*fac 4233 vect(ig+npw-1,iblocksize)=(ghc(ig+npw-1,iblocksize)-& 4234 & eval(iblocksize,iblocksize)*cg(ig+npw-1,iblocksize))*fac 4235 else 4236 pcon(ig-igs,iblocksize)=zero 4237 vect(ig,iblocksize)=zero 4238 vect(ig+npw-1,iblocksize)=zero 4239 end if 4240 else 4241 vect(ig,iblocksize)=(ghc(ig,iblocksize)-& 4242 & eval(iblocksize,iblocksize)*cg(ig,iblocksize))*pcon(ig-igs,iblocksize) 4243 vect(ig+npw-1,iblocksize)=(ghc(ig+npw-1,iblocksize)-& 4244 & eval(iblocksize,iblocksize)*cg(ig+npw-1,iblocksize))*pcon(ig-igs,iblocksize) 4245 end if 4246 end do 4247 else 4248 do ig=1+igs,npw+igs 4249 if (iterationnumber==1.or.optpcon==1) then 4250 if(kinpw(ig-igs)<huge(0.0_dp)*1.d-11)then 4251 xx=kinpw(ig-igs)*ek0_inv(iblocksize) 4252 ! teter polynomial ratio 4253 poly=27._dp+xx*(18._dp+xx*(12._dp+xx*8._dp)) 4254 fac=poly/(poly+16._dp*xx**4) 4255 if (optekin==1) fac=two*fac 4256 pcon(ig-igs,iblocksize)=fac 4257 vect(ig,iblocksize)=(ghc(ig,iblocksize)-& 4258 & eval(iblocksize,iblocksize)*cg(ig,iblocksize))*fac 4259 vect(ig+npw,iblocksize)=(ghc(ig+npw,iblocksize)-& 4260 & eval(iblocksize,iblocksize)*cg(ig+npw,iblocksize))*fac 4261 else 4262 pcon(ig-igs,iblocksize)=zero 4263 vect(ig,iblocksize)=zero 4264 vect(ig+npw,iblocksize)=zero 4265 end if 4266 else 4267 vect(ig,iblocksize)=(ghc(ig,iblocksize)-& 4268 & eval(iblocksize,iblocksize)*cg(ig,iblocksize))*pcon(ig-igs,iblocksize) 4269 vect(ig+npw,iblocksize)=(ghc(ig+npw,iblocksize)-& 4270 & eval(iblocksize,iblocksize)*cg(ig+npw,iblocksize))*pcon(ig-igs,iblocksize) 4271 end if 4272 end do 4273 end if 4274 end do 4275 end do 4276 ABI_DEALLOCATE(ek0) 4277 ABI_DEALLOCATE(ek0_inv) 4278 end if !optpcon 4279 4280 call timab(536,2,tsec) 4281 4282 end subroutine cg_precon_block
m_cgtools/cg_real_zdotc [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_real_zdotc
FUNCTION
Perform a vector-vector operation defined as res = REAL (\Sigma (conjg(x)*y)) where x and y are n-element vectors.
INPUTS
n = Specifies the number of elements in vector x and y x,y = Input arrays.
OUTPUT
res=Real part of the scalar product.
PARENTS
SOURCE
738 function cg_real_zdotc(n,x,y) result(res) 739 740 741 !This section has been created automatically by the script Abilint (TD). 742 !Do not modify the following lines by hand. 743 #undef ABI_FUNC 744 #define ABI_FUNC 'cg_real_zdotc' 745 !End of the abilint section 746 747 implicit none 748 749 !Arguments ------------------------------------ 750 !scalars 751 integer,intent(in) :: n 752 !arrays 753 real(dp),intent(in) :: x(2,n) 754 real(dp),intent(in) :: y(2,n) 755 real(dp) :: res 756 757 !Local variables------------------------------- 758 real(dp),external :: ddot 759 760 ! ************************************************************************* 761 762 res = ddot(2*n,x,1,y,1) 763 764 end function cg_real_zdotc
m_cgtools/cg_setaug_zero [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_setaug_zero
FUNCTION
Set to zero all elements of the array that are not in the FFT box.
INPUTS
nx,ny,nz=physical dimensions of the FFT box ldx,ldy,ldx=memory dimension of arr ndat=number of FFTs SIDE EFFECT arr(2,ldx,ldy,ldz*ndat)= all entries in the augmented region are set to zero
PARENTS
SOURCE
360 pure subroutine cg_setaug_zero(cplex,nx,ny,nz,ldx,ldy,ldz,ndat,arr) 361 362 363 !This section has been created automatically by the script Abilint (TD). 364 !Do not modify the following lines by hand. 365 #undef ABI_FUNC 366 #define ABI_FUNC 'cg_setaug_zero' 367 !End of the abilint section 368 369 implicit none 370 371 !Arguments ------------------------------------ 372 !scalars 373 integer,intent(in) :: cplex,nx,ny,nz,ldx,ldy,ldz,ndat 374 !arrays 375 real(dp),intent(inout) :: arr(cplex,ldx,ldy,ldz*ndat) 376 377 !Local variables------------------------------- 378 integer :: iy,iz,dat,padat 379 380 ! ************************************************************************* 381 382 if (nx /= ldx) then 383 do iz=1,ldz*ndat 384 do iy=1,ldy 385 arr(:,nx+1:ldx,iy,iz) = zero 386 end do 387 end do 388 end if 389 390 if (ny /= ldy) then 391 do iz=1,ldz*ndat 392 arr(:,:,ny+1:ldy,iz) = zero 393 end do 394 end if 395 396 if (nz /= ldz) then 397 do dat=1,ndat 398 padat = ldz*(dat-1) 399 do iz=nz+1,ldz 400 arr(:,:,:,iz+padat) = zero 401 end do 402 end do 403 end if 404 405 end subroutine cg_setaug_zero
m_cgtools/cg_setval [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_setval
FUNCTION
Set cg=alpha.
INPUTS
OUTPUT
PARENTS
CHILDREN
SOURCE
142 subroutine cg_setval(n,cg,alpha) 143 144 145 !This section has been created automatically by the script Abilint (TD). 146 !Do not modify the following lines by hand. 147 #undef ABI_FUNC 148 #define ABI_FUNC 'cg_setval' 149 !End of the abilint section 150 151 implicit none 152 153 !Arguments ------------------------------------ 154 !scalars 155 integer,intent(in) :: n 156 !arrays 157 real(dp),optional,intent(in) :: alpha(2) 158 real(dp),intent(inout) :: cg(2,n) 159 160 ! ************************************************************************* 161 162 if (PRESENT(alpha)) then 163 !$OMP PARALLEL 164 !$OMP WORKSHARE 165 cg(1,:)=alpha(1) 166 cg(2,:)=alpha(2) 167 !$OMP END WORKSHARE 168 !$OMP END PARALLEL 169 else 170 !$OMP PARALLEL 171 !$OMP WORKSHARE 172 cg(:,:)=zero 173 !$OMP END WORKSHARE 174 !$OMP END PARALLEL 175 end if 176 177 end subroutine cg_setval
m_cgtools/cg_to_reim [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_to_reim
FUNCTION
INPUTS
PARENTS
CHILDREN
SOURCE
424 subroutine cg_to_reim(npw,ndat,cg,factor,reim) 425 426 427 !This section has been created automatically by the script Abilint (TD). 428 !Do not modify the following lines by hand. 429 #undef ABI_FUNC 430 #define ABI_FUNC 'cg_to_reim' 431 !End of the abilint section 432 433 implicit none 434 435 !Arguments ------------------------------------ 436 !scalars 437 integer,intent(in) :: npw,ndat 438 real(dp),intent(in) :: factor 439 !arrays 440 real(dp),intent(in) :: cg(2*npw*ndat) 441 real(dp),intent(out) :: reim(npw*ndat*2) 442 443 ! ************************************************************************* 444 445 ! Pack real and imaginary part of the wavefunctions. 446 ! and multiply by scale factor if factor /= one. Could block but oh well 447 call dcopy(npw*ndat, cg(1), 2, reim(1), 1) 448 if (factor /= one) call dscal(npw*ndat,factor,reim(1),1) 449 450 call dcopy(npw*ndat, cg(2), 2, reim(npw*ndat+1), 1) 451 if (factor /= one) call dscal(npw*ndat,factor,reim(npw*ndat+1),1) 452 453 end subroutine cg_to_reim
m_cgtools/cg_tocplx [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_tocplx
FUNCTION
Convert a real array with (real,imag) part to complex.
INPUTS
n = Specifies the number of elements in cg and ocplx cg(2*n)=Input array with real and imaginary part.
OUTPUT
ocplx(n)=Output complex array.
PARENTS
CHILDREN
SOURCE
202 subroutine cg_tocplx(n, cg, ocplx) 203 204 205 !This section has been created automatically by the script Abilint (TD). 206 !Do not modify the following lines by hand. 207 #undef ABI_FUNC 208 #define ABI_FUNC 'cg_tocplx' 209 !End of the abilint section 210 211 implicit none 212 213 !Arguments ------------------------------------ 214 !scalars 215 integer,intent(in) :: n 216 !arrays 217 real(dp),intent(in) :: cg(2*n) 218 complex(dpc),intent(out) :: ocplx(n) 219 220 !Local variables ------------------------------ 221 !scalars 222 integer :: ii,idx 223 224 ! ************************************************************************* 225 226 !$OMP PARALLEL DO PRIVATE(ii,idx) 227 do ii=1,n 228 idx = 2*ii-1 229 ocplx(ii) = DCMPLX(cg(idx),cg(idx+1)) 230 end do 231 232 end subroutine cg_tocplx
m_cgtools/cg_vlocpsi [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_vlocpsi
FUNCTION
Apply the local part of the potentatil to the wavefunction in real space.
INPUTS
nx,ny,nz=physical dimension of the FFT box. ldx,ldy,ldz=leading dimensions of the arrays. ndat=number of wavefunctions. cplex= 1 if vloc is real, 2 for complex vloc(cplex*ldx,ldy,ldz)=Local potential on the FFT box.
SIDE EFFECTS
ur(2,ldx,ldy,ldz*ndat)=
Input = wavefunctions in real space.
Output= vloc |ur>
PARENTS
m_dfti,m_fftw3
CHILDREN
SOURCE
2553 subroutine cg_vlocpsi(nx,ny,nz,ldx,ldy,ldz,ndat,cplex,vloc,ur) 2554 2555 2556 !This section has been created automatically by the script Abilint (TD). 2557 !Do not modify the following lines by hand. 2558 #undef ABI_FUNC 2559 #define ABI_FUNC 'cg_vlocpsi' 2560 !End of the abilint section 2561 2562 implicit none 2563 2564 !Arguments ------------------------------------ 2565 !scalars 2566 integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,ndat,cplex 2567 !arrays 2568 real(dp),intent(in) :: vloc(cplex*ldx,ldy,ldz) 2569 real(dp),intent(inout) :: ur(2,ldx,ldy,ldz*ndat) 2570 2571 !Local variables------------------------------- 2572 !scalars 2573 integer :: idat,ix,iy,iz,padat 2574 real(dp) :: fim,fre 2575 2576 ! ************************************************************************* 2577 2578 if (cplex==1) then 2579 ! 2580 if (ndat==1) then 2581 !$OMP PARALLEL DO 2582 do iz=1,nz 2583 do iy=1,ny 2584 do ix=1,nx 2585 ur(1,ix,iy,iz) = vloc(ix,iy,iz) * ur(1,ix,iy,iz) 2586 ur(2,ix,iy,iz) = vloc(ix,iy,iz) * ur(2,ix,iy,iz) 2587 end do 2588 end do 2589 end do 2590 ! 2591 else 2592 ! 2593 !$OMP PARALLEL DO PRIVATE(padat) 2594 do idat=1,ndat 2595 padat = ldz*(idat-1) 2596 do iz=1,nz 2597 do iy=1,ny 2598 do ix=1,nx 2599 ur(1,ix,iy,iz+padat) = vloc(ix,iy,iz) * ur(1,ix,iy,iz+padat) 2600 ur(2,ix,iy,iz+padat) = vloc(ix,iy,iz) * ur(2,ix,iy,iz+padat) 2601 end do 2602 end do 2603 end do 2604 end do 2605 ! 2606 end if 2607 ! 2608 else if (cplex==2)then 2609 ! 2610 if (ndat==1) then 2611 !$OMP PARALLEL DO PRIVATE(fre,fim) 2612 do iz=1,nz 2613 do iy=1,ny 2614 do ix=1,nx 2615 fre = ur(1,ix,iy,iz) 2616 fim = ur(2,ix,iy,iz) 2617 ur(1,ix,iy,iz) = vloc(2*ix-1,iy,iz)*fre - vloc(2*ix,iy,iz)*fim 2618 ur(2,ix,iy,iz) = vloc(2*ix-1,iy,iz)*fim + vloc(2*ix,iy,iz)*fre 2619 end do 2620 end do 2621 end do 2622 else 2623 !$OMP PARALLEL DO PRIVATE(padat,fre,fim) 2624 do idat=1,ndat 2625 padat = ldz*(idat-1) 2626 do iz=1,nz 2627 do iy=1,ny 2628 do ix=1,nx 2629 fre = ur(1,ix,iy,iz+padat) 2630 fim = ur(2,ix,iy,iz+padat) 2631 ur(1,ix,iy,iz+padat) = vloc(2*ix-1,iy,iz)*fre - vloc(2*ix,iy,iz)*fim 2632 ur(2,ix,iy,iz+padat) = vloc(2*ix-1,iy,iz)*fim + vloc(2*ix,iy,iz)*fre 2633 end do 2634 end do 2635 end do 2636 end do 2637 end if 2638 ! 2639 else 2640 ur = huge(one) 2641 !MSG_BUG("Wrong cplex") 2642 end if 2643 2644 end subroutine cg_vlocpsi
m_cgtools/cg_zaxpby [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_zaxpby
FUNCTION
Scales two vectors, adds them to one another and stores result in the vector. y := a*x + b*y
INPUTS
n = the number of elements in vectors x and y. a = Specifies the scalar a. x = Array. b = Specifies the scalar b. y = Array
OUTPUT
y Contains the updated vector y.
PARENTS
CHILDREN
SOURCE
916 subroutine cg_zaxpby(n,a,x,b,y) 917 918 919 !This section has been created automatically by the script Abilint (TD). 920 !Do not modify the following lines by hand. 921 #undef ABI_FUNC 922 #define ABI_FUNC 'cg_zaxpby' 923 !End of the abilint section 924 925 implicit none 926 927 !Arguments ------------------------------------ 928 !scalars 929 integer,intent(in) :: n 930 real(dp),intent(in) :: a(2),b(2) 931 !arrays 932 real(dp),intent(in) :: x(2*n) 933 real(dp),intent(inout) :: y(2*n) 934 935 ! ************************************************************************* 936 937 #ifdef HAVE_LINALG_AXPBY 938 call zaxpby(n, a, x, 1, b, y, 1) 939 #else 940 call zscal(n, b, y, 1) 941 call zaxpy(n, a, x, 1, y,1) 942 #endif 943 944 end subroutine cg_zaxpby
m_cgtools/cg_zaxpy [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_zaxpy
FUNCTION
Computes y = alpha*x + y
INPUTS
n = Specifies the number of elements in vectors x and y. alpha = Specifies the scalar alpha. x = Array
SIDE EFFECTS
y = Array. In output, y contains the updated vector.
PARENTS
cgwf,dfpt_cgwf,lapackprof,rf2_init
CHILDREN
SOURCE
857 subroutine cg_zaxpy(n,alpha,x,y) 858 859 860 !This section has been created automatically by the script Abilint (TD). 861 !Do not modify the following lines by hand. 862 #undef ABI_FUNC 863 #define ABI_FUNC 'cg_zaxpy' 864 !End of the abilint section 865 866 implicit none 867 868 !Arguments ------------------------------------ 869 !scalars 870 integer,intent(in) :: n 871 real(dp),intent(in) :: alpha(2) 872 !arrays 873 real(dp),intent(in) :: x(2*n) 874 real(dp),intent(inout) :: y(2*n) 875 876 !local variables 877 ! integer :: ii 878 879 ! ************************************************************************* 880 881 if (alpha(2) == zero) then 882 call daxpy(2*n,alpha(1),x,1,y,1) 883 else 884 call zaxpy(n,alpha,x,1,y,1) 885 end if 886 887 end subroutine cg_zaxpy
m_cgtools/cg_zcopy [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_zcopy
FUNCTION
Perform y = x, where x and y are vectors.
INPUTS
n = Specifies the number of elements in vectors x and y. x = Input Array
OUTPUT
y = In output, y contains a copy of the values of x.
PARENTS
cgwf,corrmetalwf1,dfpt_cgwf,dfpt_mkrho,dfpt_vtowfk,lapackprof,m_iowf
CHILDREN
SOURCE
526 subroutine cg_zcopy(n, x, y) 527 528 529 !This section has been created automatically by the script Abilint (TD). 530 !Do not modify the following lines by hand. 531 #undef ABI_FUNC 532 #define ABI_FUNC 'cg_zcopy' 533 !End of the abilint section 534 535 implicit none 536 537 !Arguments ------------------------------------ 538 !scalars 539 integer,intent(in) :: n 540 !arrays 541 real(dp),intent(in) :: x(2*n) 542 real(dp),intent(out) :: y(2*n) 543 544 ! ************************************************************************* 545 546 call zcopy(n,x,1,y,1) 547 548 end subroutine cg_zcopy
m_cgtools/cg_zdotc [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_zdotc
FUNCTION
Perform a vector-vector operation defined as res = \Sigma (conjg(x)*y) where x and y are n-element vectors.
INPUTS
n = Specifies the number of elements in vector x and y x,y = Input arrays.
OUTPUT
res(2)=Real and Imaginary part of the scalar product.
PARENTS
SOURCE
672 function cg_zdotc(n,x,y) result(res) 673 674 675 !This section has been created automatically by the script Abilint (TD). 676 !Do not modify the following lines by hand. 677 #undef ABI_FUNC 678 #define ABI_FUNC 'cg_zdotc' 679 !End of the abilint section 680 681 implicit none 682 683 !Arguments ------------------------------------ 684 !scalars 685 integer,intent(in) :: n 686 !arrays 687 real(dp),intent(in) :: x(2,n) 688 real(dp),intent(in) :: y(2,n) 689 real(dp) :: res(2) 690 691 !Local variables------------------------------- 692 #ifdef HAVE_LINALG_ZDOTC_BUG 693 integer :: ii 694 #endif 695 complex(dpc) :: cres 696 complex(dpc),external :: zdotc 697 698 ! ************************************************************************* 699 700 #ifdef HAVE_LINALG_ZDOTC_BUG 701 ! Workaround for veclib on MacOSx 702 res = zero 703 !$OMP PARALLEL DO PRIVATE(ii) REDUCTION(+:res) 704 do ii=1,n 705 res(1) = res(1) + x(1,ii)*y(1,ii) + x(2,ii)*y(2,ii) 706 res(2) = res(2) + x(1,ii)*y(2,ii) - x(2,ii)*y(1,ii) 707 end do 708 709 #else 710 cres = zdotc(n, x, 1, y, 1) 711 res(1) = REAL(cres) 712 res(2) = AIMAG(cres) 713 #endif 714 715 end function cg_zdotc
m_cgtools/cg_zdotu [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_zdotu
FUNCTION
Perform a vector-vector operation defined as res = \Sigma (x*y) where x and y are n-element vectors. Note that x is unconjugated.
INPUTS
n = Specifies the number of elements in vector x and y x,y = Input arrays.
OUTPUT
res(2)=Real and Imaginary part of the scalar product.
PARENTS
SOURCE
788 function cg_zdotu(n, x, y) result(res) 789 790 791 !This section has been created automatically by the script Abilint (TD). 792 !Do not modify the following lines by hand. 793 #undef ABI_FUNC 794 #define ABI_FUNC 'cg_zdotu' 795 !End of the abilint section 796 797 implicit none 798 799 !Arguments ------------------------------------ 800 !scalars 801 integer,intent(in) :: n 802 !arrays 803 real(dp),intent(in) :: x(2,n) 804 real(dp),intent(in) :: y(2,n) 805 real(dp) :: res(2) 806 807 !Local variables------------------------------- 808 #ifdef HAVE_LINALG_ZDOTU_BUG 809 integer :: ii 810 #endif 811 complex(dpc) :: cres 812 complex(dpc),external :: zdotu 813 814 ! ************************************************************************* 815 816 #ifdef HAVE_LINALG_ZDOTU_BUG 817 ! Workaround for veclib on MacOSx 818 res = zero 819 !$OMP PARALLEL DO PRIVATE(ii) REDUCTION(+:res) 820 do ii=1,n 821 res(1) = res(1) + x(1,ii)*y(1,ii) - x(2,ii)*y(2,ii) 822 res(2) = res(2) + x(1,ii)*y(2,ii) + x(2,ii)*y(1,ii) 823 end do 824 #else 825 cres = zdotu(n, x, 1, y, 1) 826 res(1) = REAL(cres) 827 res(2) = AIMAG(cres) 828 #endif 829 830 end function cg_zdotu
m_cgtools/cg_zgemm [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_zgemm
FUNCTION
The ?gemm routines perform a matrix-matrix operation with general matrices. The operation is defined as C := alpha*op(A)*op(B) + beta*C, where: op(x) is one of op(x) = x, or op(x) = x', or op(x) = conjg(x'), alpha and beta are scalars, A, B and C are matrices: op(A) is an m-by-k matrix, op(B) is a k-by-n matrix, C is an m-by-n matrix.
INPUTS
OUTPUT
PARENTS
lapackprof,m_cgtools
CHILDREN
SOURCE
1057 subroutine cg_zgemm(transa,transb,npws,ncola,ncolb,cg_a,cg_b,cg_c,alpha,beta) 1058 1059 1060 !This section has been created automatically by the script Abilint (TD). 1061 !Do not modify the following lines by hand. 1062 #undef ABI_FUNC 1063 #define ABI_FUNC 'cg_zgemm' 1064 !End of the abilint section 1065 1066 implicit none 1067 1068 !Arguments ------------------------------------ 1069 !scalars 1070 integer,intent(in) :: npws,ncola,ncolb 1071 real(dp),optional,intent(in) :: alpha(2),beta(2) 1072 character(len=1),intent(in) :: transa,transb 1073 !arrays 1074 real(dp),intent(in) :: cg_a(2,npws*ncola) 1075 real(dp),intent(in) :: cg_b(2,npws*ncolb) 1076 real(dp),intent(inout) :: cg_c(2,*) 1077 1078 !Local variables------------------------------- 1079 !scalars 1080 integer :: mm,nn,kk,lda,ldb,ldc 1081 real(dp) :: my_alpha(2),my_beta(2) 1082 1083 ! ************************************************************************* 1084 1085 lda = npws 1086 ldb = npws 1087 1088 mm = npws 1089 nn = ncolb 1090 kk = ncola 1091 1092 if (toupper(transa) /= 'N') then 1093 mm = ncola 1094 kk = npws 1095 end if 1096 if (toupper(transb) /= 'N') nn = npws 1097 1098 ldc = mm 1099 1100 my_alpha = cg_cone; if (PRESENT(alpha)) my_alpha = alpha 1101 my_beta = cg_czero; if (PRESENT(beta)) my_beta = beta 1102 1103 call ZGEMM(transa,transb,mm,nn,kk,my_alpha,cg_a,lda,cg_b,ldb,my_beta,cg_c,ldc) 1104 1105 end subroutine cg_zgemm
m_cgtools/cg_zgemv [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_zgemv
FUNCTION
The ?gemv routines perform a matrix-vector operation defined as y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, or y := alpha*conjg(A')*x + beta*y, where: alpha and beta are scalars, x and y are vectors, A is an m-by-n matrix.
INPUTS
OUTPUT
PARENTS
lapackprof,m_cgtools
CHILDREN
SOURCE
975 subroutine cg_zgemv(trans,nrows,ncols,cgmat,vec,matvec,alpha,beta) 976 977 978 !This section has been created automatically by the script Abilint (TD). 979 !Do not modify the following lines by hand. 980 #undef ABI_FUNC 981 #define ABI_FUNC 'cg_zgemv' 982 !End of the abilint section 983 984 implicit none 985 986 !Arguments ------------------------------------ 987 !scalars 988 integer,intent(in) :: nrows,ncols 989 real(dp),optional,intent(in) :: alpha(2),beta(2) 990 character(len=1),intent(in) :: trans 991 !arrays 992 real(dp),intent(in) :: cgmat(2,nrows*ncols) 993 real(dp),intent(in) :: vec(2,*) 994 real(dp),intent(inout) :: matvec(2,*) 995 996 !Local variables------------------------------- 997 !scalars 998 integer :: mm,nn,kk,lda,ldb,ldc 999 real(dp) :: my_alpha(2),my_beta(2) 1000 1001 ! ************************************************************************* 1002 1003 lda = nrows 1004 mm = nrows 1005 nn = 1 1006 kk = ncols 1007 1008 if (toupper(trans) /= 'N') then 1009 mm = ncols 1010 kk = nrows 1011 end if 1012 1013 ldb = kk 1014 ldc = mm 1015 1016 my_alpha = cg_cone; if (PRESENT(alpha)) my_alpha = alpha 1017 my_beta = cg_czero; if (PRESENT(beta)) my_beta = beta 1018 1019 call ZGEMM(trans,"N",mm,nn,kk,my_alpha,cgmat,lda,vec,ldb,my_beta,matvec,ldc) 1020 ! ZGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC) 1021 1022 !call ZGEMV(trans,mm,nn,my_alpha,cgmat,lda,vec,1,my_beta,matvec,1) 1023 1024 end subroutine cg_zgemv
m_cgtools/cg_zscal [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_zscal
FUNCTION
Perform x = a*x
INPUTS
n = Specifies the number of elements in vector x. a(2)= The scalar a. If a(2) is zero, x = a*x is computed via zdscal
OUTPUT
x = Updated vector.
OUTPUT
PARENTS
cgwf,m_cgtools
CHILDREN
SOURCE
576 subroutine cg_zscal(n, a, x) 577 578 579 !This section has been created automatically by the script Abilint (TD). 580 !Do not modify the following lines by hand. 581 #undef ABI_FUNC 582 #define ABI_FUNC 'cg_zscal' 583 !End of the abilint section 584 585 implicit none 586 587 !Arguments ------------------------------------ 588 !scalars 589 integer,intent(in) :: n 590 real(dp),intent(in) :: a(2) 591 !arrays 592 real(dp),intent(inout) :: x(2*n) 593 594 ! ************************************************************************* 595 596 if (a(2) == zero) then 597 call zdscal(n, a, x, 1) 598 else 599 call zscal(n, a, x, 1) 600 end if 601 602 end subroutine cg_zscal
m_cgtools/cgnc_cholesky [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cgnc_cholesky
FUNCTION
Cholesky orthonormalization of the vectors stored in cgblock (version optimized for NC wavefunctions).
INPUTS
npws=Size of each vector (usually npw*nspinor) nband=Number of band in cgblock istwfk=Storage mode for the wavefunctions. 1 for standard full mode me_g0=1 if this node has G=0. comm_pw=MPI communicator for the planewave group. Set to xmpi_comm_self for sequential mode.
SIDE EFFECTS
cgblock(2*npws*nband)
input: Input set of vectors.
output: Orthonormalized set.
PARENTS
pw_orthon
CHILDREN
SOURCE
2675 subroutine cgnc_cholesky(npws,nband,cgblock,istwfk,me_g0,comm_pw,use_gemm) 2676 2677 2678 !This section has been created automatically by the script Abilint (TD). 2679 !Do not modify the following lines by hand. 2680 #undef ABI_FUNC 2681 #define ABI_FUNC 'cgnc_cholesky' 2682 use interfaces_14_hidewrite 2683 !End of the abilint section 2684 2685 implicit none 2686 2687 !Arguments ------------------------------------ 2688 !scalars 2689 integer,intent(in) :: npws,nband,istwfk 2690 integer,intent(in) :: comm_pw,me_g0 2691 logical,optional,intent(in) :: use_gemm 2692 !arrays 2693 real(dp),intent(inout) :: cgblock(2*npws*nband) 2694 2695 !Local variables ------------------------------ 2696 !scalars 2697 integer :: ierr,b1,b2 2698 #ifdef DEBUG_MODE 2699 integer :: ptr 2700 #endif 2701 logical :: my_usegemm 2702 character(len=500) :: msg 2703 !arrays 2704 real(dp) :: rcg0(nband) 2705 real(dp),allocatable :: rovlp(:,:) 2706 complex(dpc),allocatable :: cf_ovlp(:,:) 2707 2708 ! ************************************************************************* 2709 2710 #ifdef DEBUG_MODE 2711 if (istwfk==2) then 2712 ierr = 0 2713 do b1=1,nband 2714 ptr = 2 + 2*(b1-1)*npws 2715 if (ABS(cgblock(ptr)) > zero) then 2716 ierr = ierr + 1 2717 write(msg,'(a,i0,es13.6)')" Input b1, Im u(g=0) should be zero ",b1,cgblock(ptr) 2718 call wrtout(std_out,msg,"COLL") 2719 !cgblock(ptr) = zero 2720 end if 2721 end do 2722 ABI_CHECK(ierr==0,"Non zero imag part") 2723 end if 2724 #endif 2725 2726 my_usegemm=.FALSE.; if (PRESENT(use_gemm)) my_usegemm = use_gemm 2727 2728 ABI_MALLOC(cf_ovlp,(nband,nband)) 2729 ! 2730 ! 1) Calculate O_ij = <phi_i|phi_j> 2731 if (my_usegemm) then 2732 call ABI_ZGEMM("Conjugate","Normal",nband,nband,npws,cone,cgblock,npws,cgblock,npws,czero,cf_ovlp,nband) 2733 else 2734 call ZHERK("U","C",nband,npws,one,cgblock,npws,zero,cf_ovlp,nband) 2735 end if 2736 2737 if (istwfk==1) then 2738 ! 2739 ! Sum the overlap if PW are distributed. 2740 if (comm_pw /= xmpi_comm_self) then 2741 call xmpi_sum(cf_ovlp,comm_pw,ierr) 2742 end if 2743 ! 2744 ! 2) Cholesky factorization: O = U^H U with U upper triangle matrix. 2745 call ZPOTRF('U',nband,cf_ovlp,nband,ierr) 2746 2747 if (ierr/=0) then 2748 write(msg,'(a,i0)')' ZPOTRF returned info = ',ierr 2749 MSG_ERROR(msg) 2750 end if 2751 2752 else 2753 ! overlap is real. Note that nspinor is always 1 in this case. 2754 ABI_MALLOC(rovlp,(nband,nband)) 2755 rovlp = two * REAL(cf_ovlp) 2756 2757 if (istwfk==2 .and. me_g0==1) then 2758 ! Extract the real part at G=0 and subtract its contribution to the overlap. 2759 call dcopy(nband,cgblock,2*npws,rcg0,1) 2760 do b2=1,nband 2761 do b1=1,b2 2762 rovlp(b1,b2) = rovlp(b1,b2) - rcg0(b1)*rcg0(b2) 2763 end do 2764 end do 2765 end if 2766 2767 ! Sum the overlap if PW are distributed. 2768 if (comm_pw /= xmpi_comm_self) then 2769 call xmpi_sum(rovlp,comm_pw,ierr) 2770 end if 2771 ! 2772 ! 2) Cholesky factorization: O = U^H U with U upper triangle matrix. 2773 call DPOTRF('U',nband,rovlp,nband,ierr) 2774 2775 if (ierr/=0) then 2776 write(msg,'(a,i0)')' DPOTRF returned info = ',ierr 2777 MSG_ERROR(msg) 2778 end if 2779 2780 cf_ovlp = DCMPLX(rovlp) 2781 ABI_FREE(rovlp) 2782 end if 2783 ! 2784 ! 3) Solve X U = cgblock. On exit cgblock is orthonormalized. 2785 call ZTRSM('Right','Upper','Normal','Normal',npws,nband,cone,cf_ovlp,nband,cgblock,npws) 2786 2787 #ifdef DEBUG_MODE 2788 if (istwfk==2) then 2789 ierr = 0 2790 do b1=1,nband 2791 ptr = 2 + 2*(b1-1)*npws 2792 if (ABS(cgblock(ptr)) > zero) then 2793 ierr = ierr + 1 2794 write(msg,'(a,i0,es13.6)')" Output b1, Im u(g=0) should be zero ",b1,cgblock(ptr) 2795 end if 2796 end do 2797 ABI_CHECK(ierr==0,"Non zero imag part") 2798 end if 2799 #endif 2800 2801 ABI_FREE(cf_ovlp) 2802 2803 end subroutine cgnc_cholesky
m_cgtools/cgnc_gramschmidt [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cgnc_grortho
FUNCTION
Gram-Schmidt orthonormalization of the vectors stored in cgblock
INPUTS
npws=Size of each vector (usually npw*nspinor) nband=Number of band in cgblock istwfk=Storage mode for the wavefunctions. 1 for standard full mode me_g0=1 if this node has G=0. comm_pw=MPI communicator for the planewave group. Set to xmpi_comm_self for sequential mode.
SIDE EFFECTS
cgblock(2*npws*nband)
input: Input set of vectors.
output: Orthonormalized set.
PARENTS
SOURCE
3148 subroutine cgnc_gramschmidt(npws,nband,cgblock,istwfk,me_g0,comm_pw) 3149 3150 3151 !This section has been created automatically by the script Abilint (TD). 3152 !Do not modify the following lines by hand. 3153 #undef ABI_FUNC 3154 #define ABI_FUNC 'cgnc_gramschmidt' 3155 !End of the abilint section 3156 3157 implicit none 3158 3159 !Arguments ------------------------------------ 3160 !scalars 3161 integer,intent(in) :: npws,nband,istwfk 3162 integer,intent(in) :: comm_pw,me_g0 3163 !arrays 3164 real(dp),intent(inout) :: cgblock(2*npws*nband) 3165 3166 !Local variables ------------------------------ 3167 !scalars 3168 integer :: b1,nb2,opt 3169 logical :: normalize 3170 3171 ! ************************************************************************* 3172 3173 ! Normalize the first vector. 3174 call cgnc_normalize(npws,1,cgblock(1),istwfk,me_g0,comm_pw) 3175 if (nband == 1) RETURN 3176 3177 ! Orthogonaluze b1 wrt to the bands in [1,b1-1]. 3178 normalize = .TRUE. 3179 do b1=2,nband 3180 opt = 1 + 2*npws*(b1-1) 3181 nb2=b1-1 3182 call cgnc_gsortho(npws,nb2,cgblock(1),1,cgblock(opt),istwfk,normalize,me_g0,comm_pw) 3183 end do 3184 3185 end subroutine cgnc_gramschmidt
m_cgtools/cgnc_gsortho [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cgnc_gsortho
FUNCTION
INPUTS
npws=Size of each vector (usually npw*nspinor) nband1=Number of vectors in icg1 nband1=Number of vectors in cg2 comm_pw=MPI communicator.
SIDE EFFECTS
cg2(2*npws*nband2)
icg1(2*npws*nband1)
input: Input set of vectors.
output: Orthonormalized set.
PARENTS
m_cgtools
CHILDREN
SOURCE
3051 subroutine cgnc_gsortho(npws,nband1,icg1,nband2,iocg2,istwfk,normalize,me_g0,comm_pw) 3052 3053 3054 !This section has been created automatically by the script Abilint (TD). 3055 !Do not modify the following lines by hand. 3056 #undef ABI_FUNC 3057 #define ABI_FUNC 'cgnc_gsortho' 3058 !End of the abilint section 3059 3060 implicit none 3061 3062 !Arguments ------------------------------------ 3063 !scalars 3064 integer,intent(in) :: npws,nband1,nband2,istwfk,me_g0 3065 integer,optional,intent(in) :: comm_pw 3066 logical,intent(in) :: normalize 3067 !arrays 3068 real(dp),intent(in) :: icg1(2*npws*nband1) 3069 real(dp),intent(inout) :: iocg2(2*npws*nband2) 3070 3071 !Local variables ------------------------------ 3072 !scalars 3073 integer :: ierr,b1,b2 3074 !arrays 3075 real(dp) :: r_icg1(nband1),r_iocg2(nband2) 3076 real(dp),allocatable :: proj(:,:,:) 3077 3078 ! ************************************************************************* 3079 3080 ABI_MALLOC(proj,(2,nband1,nband2)) 3081 !proj = zero 3082 3083 ! 1) Calculate <cg1|cg2> 3084 call cg_zgemm("C","N",npws,nband1,nband2,icg1,iocg2,proj) 3085 3086 if (istwfk>1) then 3087 ! nspinor is always 1 in this case. 3088 ! Account for the missing G and set the imaginary part to zero since wavefunctions are real. 3089 proj(1,:,:) = two * proj(1,:,:) 3090 proj(2,:,:) = zero 3091 ! 3092 if (istwfk==2 .and. me_g0==1) then 3093 ! Extract the real part at G=0 and subtract its contribution. 3094 call dcopy(nband1,icg1, 2*npws,r_icg1, 1) 3095 call dcopy(nband2,iocg2,2*npws,r_iocg2,1) 3096 do b2=1,nband2 3097 do b1=1,nband1 3098 proj(1,b1,b2) = proj(1,b1,b2) - r_icg1(b1) * r_iocg2(b2) 3099 end do 3100 end do 3101 end if 3102 ! 3103 end if 3104 ! 3105 ! This is for the MPI version 3106 if (comm_pw /= xmpi_comm_self) then 3107 call xmpi_sum(proj,comm_pw,ierr) 3108 end if 3109 3110 ! 2) cg2 = cg2 - <cg1|cg2> cg1 3111 call cg_zgemm("N","N",npws,nband1,nband2,icg1,proj,iocg2,alpha=-cg_cone,beta=cg_cone) 3112 3113 ABI_FREE(proj) 3114 3115 ! 3) Normalize iocg2 if required. 3116 if (normalize) then 3117 call cgnc_normalize(npws,nband2,iocg2,istwfk,me_g0,comm_pw) 3118 end if 3119 3120 end subroutine cgnc_gsortho
m_cgtools/cgnc_normalize [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cgnc_normalize
FUNCTION
INPUTS
npws=Size of each vector (usually npw*nspinor) nband=Number of vectors in icg1
SIDE EFFECTS
PARENTS
m_cgtools
CHILDREN
SOURCE
2951 subroutine cgnc_normalize(npws,nband,cg,istwfk,me_g0,comm_pw) 2952 2953 2954 !This section has been created automatically by the script Abilint (TD). 2955 !Do not modify the following lines by hand. 2956 #undef ABI_FUNC 2957 #define ABI_FUNC 'cgnc_normalize' 2958 !End of the abilint section 2959 2960 implicit none 2961 2962 !Arguments ------------------------------------ 2963 !scalars 2964 integer,intent(in) :: npws,nband,istwfk,me_g0,comm_pw 2965 !arrays 2966 real(dp),intent(inout) :: cg(2*npws*nband) 2967 2968 !Local variables ------------------------------ 2969 !scalars 2970 integer :: ptr,ierr,band 2971 character(len=500) :: msg 2972 !arrays 2973 real(dp) :: norm(nband),alpha(2) 2974 2975 ! ************************************************************************* 2976 2977 !$OMP PARALLEL DO PRIVATE(ptr) IF (nband > 1) 2978 do band=1,nband 2979 ptr = 1 + 2*npws*(band-1) 2980 norm(band) = cg_dznrm2(npws, cg(ptr)) 2981 norm(band) = norm(band) ** 2 2982 !norm(band) = cg_real_zdotc(npws, cg(ptr), cg(ptr)) 2983 end do 2984 2985 if (istwfk>1) then 2986 norm = two * norm 2987 if (istwfk==2 .and. me_g0==1) then 2988 !$OMP PARALLEL DO PRIVATE(ptr) IF (nband >1) 2989 do band=1,nband 2990 ptr = 1 + 2*npws*(band-1) 2991 norm(band) = norm(band) - cg(ptr)**2 2992 end do 2993 end if 2994 end if 2995 2996 if (comm_pw /= xmpi_comm_self) then 2997 call xmpi_sum(norm,comm_pw,ierr) 2998 end if 2999 3000 ierr = 0 3001 do band=1,nband 3002 if (norm(band) > zero) then 3003 norm(band) = SQRT(norm(band)) 3004 else 3005 ierr = ierr + 1 3006 end if 3007 end do 3008 3009 if (ierr/=0) then 3010 write(msg,'(a,i0,a)')" Found ",ierr," vectors with norm <= zero!" 3011 MSG_ERROR(msg) 3012 end if 3013 3014 !$OMP PARALLEL DO PRIVATE(ptr,alpha) IF (nband > 1) 3015 do band=1,nband 3016 ptr = 1 + 2*npws*(band-1) 3017 alpha = (/one/norm(band), zero/) 3018 call cg_zscal(npws,alpha,cg(ptr)) 3019 end do 3020 3021 end subroutine cgnc_normalize
m_cgtools/cgpaw_gramschmidt [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cgpaw_grortho
FUNCTION
Gram-Schmidt orthonormalization of the vectors stored in cg
INPUTS
npws=Size of each vector (usually npw*nspinor) nband=Number of bands in cg istwfk=Storage mode for the wavefunctions. 1 for standard full mode me_g0=1 if this node has G=0. comm_pw=MPI communicator for the planewave group. Set to xmpi_comm_self for sequential mode.
SIDE EFFECTS
cg(2*npws*nband), gsc(2*npws*nband)
input: Input set of vectors.
output: Orthonormalized set.
PARENTS
SOURCE
3424 subroutine cgpaw_gramschmidt(npws,nband,cg,gsc,istwfk,me_g0,comm_pw) 3425 3426 3427 !This section has been created automatically by the script Abilint (TD). 3428 !Do not modify the following lines by hand. 3429 #undef ABI_FUNC 3430 #define ABI_FUNC 'cgpaw_gramschmidt' 3431 !End of the abilint section 3432 3433 implicit none 3434 3435 !Arguments ------------------------------------ 3436 !scalars 3437 integer,intent(in) :: npws,nband,istwfk,comm_pw,me_g0 3438 !arrays 3439 real(dp),intent(inout) :: cg(2*npws*nband),gsc(2*npws*nband) 3440 3441 !Local variables ------------------------------ 3442 !scalars 3443 integer :: b1,nb2,opt 3444 logical :: normalize 3445 3446 ! ************************************************************************* 3447 3448 ! Normalize the first vector. 3449 call cgpaw_normalize(npws,1,cg(1),gsc(1),istwfk,me_g0,comm_pw) 3450 if (nband == 1) RETURN 3451 3452 ! Orthogonalize b1 wrt to the bands in [1,b1-1]. 3453 normalize = .TRUE. 3454 do b1=2,nband 3455 opt = 1 + 2*npws*(b1-1) 3456 nb2=b1-1 3457 call cgpaw_gsortho(npws,nb2,cg(1),gsc(1),1,cg(opt),gsc(opt),istwfk,normalize,me_g0,comm_pw) 3458 end do 3459 3460 end subroutine cgpaw_gramschmidt
m_cgtools/cgpaw_gsortho [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cgpaw_gsortho
FUNCTION
This routine uses the Gram-Schmidt method to orthogonalize a set of PAW wavefunctions. with respect to an input block of states.
INPUTS
npws=Size of each vector (usually npw*nspinor) nband1=Number of vectors in the input block icg1 icg1(2*npws*nband1)=Input block of vectors. igsc1(2*npws*nband1)= S|C> for C in icg1. nband2=Number of vectors to orthogonalize normalize=True if output wavefunction must be normalized. istwfk=Storage mode for the wavefunctions. 1 for standard full mode me_g0=1 if this node has G=0. comm_pw=MPI communicator for the planewave group. Set to xmpi_comm_self for sequential mode.
SIDE EFFECTS
iocg2(2*npws*nband2), iogsc2(2*npws*nband1)
input: set of |C> and S|C> wher |C> is the set of states to orthogonalize
output: Orthonormalized set.
PARENTS
m_cgtools
CHILDREN
SOURCE
3325 subroutine cgpaw_gsortho(npws,nband1,icg1,igsc1,nband2,iocg2,iogsc2,istwfk,normalize,me_g0,comm_pw) 3326 3327 3328 !This section has been created automatically by the script Abilint (TD). 3329 !Do not modify the following lines by hand. 3330 #undef ABI_FUNC 3331 #define ABI_FUNC 'cgpaw_gsortho' 3332 !End of the abilint section 3333 3334 implicit none 3335 3336 !Arguments ------------------------------------ 3337 !scalars 3338 integer,intent(in) :: npws,nband1,nband2,istwfk,me_g0 3339 integer,optional,intent(in) :: comm_pw 3340 logical,intent(in) :: normalize 3341 !arrays 3342 real(dp),intent(in) :: icg1(2*npws*nband1),igsc1(2*npws*nband1) 3343 real(dp),intent(inout) :: iocg2(2*npws*nband2),iogsc2(2*npws*nband2) 3344 3345 !Local variables ------------------------------ 3346 !scalars 3347 integer :: ierr,b1,b2 3348 !arrays 3349 real(dp) :: r_icg1(nband1),r_iocg2(nband2) 3350 real(dp),allocatable :: proj(:,:,:) 3351 3352 ! ************************************************************************* 3353 3354 ABI_MALLOC(proj,(2,nband1,nband2)) 3355 3356 ! 1) Calculate <cg1|cg2> 3357 call cg_zgemm("C","N",npws,nband1,nband2,igsc1,iocg2,proj) 3358 3359 if (istwfk>1) then 3360 ! nspinor is always 1 in this case. 3361 ! Account for the missing G and set the imaginary part to zero since wavefunctions are real. 3362 proj(1,:,:) = two * proj(1,:,:) 3363 proj(2,:,:) = zero 3364 ! 3365 if (istwfk==2 .and. me_g0==1) then 3366 ! Extract the real part at G=0 and subtract its contribution. 3367 call dcopy(nband1,igsc1,2*npws,r_icg1, 1) 3368 call dcopy(nband2,iocg2,2*npws,r_iocg2,1) 3369 do b2=1,nband2 3370 do b1=1,nband1 3371 proj(1,b1,b2) = proj(1,b1,b2) - r_icg1(b1) * r_iocg2(b2) 3372 end do 3373 end do 3374 end if 3375 ! 3376 end if 3377 ! 3378 ! This is for the MPI version 3379 if (comm_pw /= xmpi_comm_self) then 3380 call xmpi_sum(proj,comm_pw,ierr) 3381 end if 3382 3383 ! 2) 3384 ! cg2 = cg2 - <Scg1|cg2> cg1 3385 ! S cg2 = S cg2 - <Scg1|cg2> S cg1 3386 call cg_zgemm("N","N",npws,nband1,nband2,icg1,proj,iocg2,alpha=-cg_cone,beta=cg_cone) 3387 call cg_zgemm("N","N",npws,nband1,nband2,igsc1,proj,iogsc2,alpha=-cg_cone,beta=cg_cone) 3388 3389 ABI_FREE(proj) 3390 3391 ! 3) Normalize iocg2 and iogsc2 if required. 3392 if (normalize) then 3393 call cgpaw_normalize(npws,nband2,iocg2,iogsc2,istwfk,me_g0,comm_pw) 3394 end if 3395 3396 end subroutine cgpaw_gsortho
m_cgtools/cgpaw_normalize [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cgpaw_normalize
FUNCTION
Normalize a set of PAW pseudo wavefunctions.
INPUTS
npws=Size of each vector (usually npw*nspinor) nband=Number of band in cgblock and gsc istwfk=Storage mode for the wavefunctions. 1 for standard full mode me_g0=1 if this node has G=0. comm_pw=MPI communicator for the planewave group. Set to xmpi_comm_self for sequential mode.
SIDE EFFECTS
cg(2*npws*nband)
input: Input set of vectors |C>
output: Normalized set such as <C|S|C> = 1
gsc(2*npws*nband)
input: Input set of vectors S|C>
output: New S|C> compute with the new |C>
PARENTS
m_cgtools
CHILDREN
SOURCE
3219 subroutine cgpaw_normalize(npws,nband,cg,gsc,istwfk,me_g0,comm_pw) 3220 3221 3222 !This section has been created automatically by the script Abilint (TD). 3223 !Do not modify the following lines by hand. 3224 #undef ABI_FUNC 3225 #define ABI_FUNC 'cgpaw_normalize' 3226 !End of the abilint section 3227 3228 implicit none 3229 3230 !Arguments ------------------------------------ 3231 !scalars 3232 integer,intent(in) :: npws,nband,istwfk,me_g0,comm_pw 3233 !arrays 3234 real(dp),intent(inout) :: cg(2*npws*nband),gsc(2*npws*nband) 3235 3236 !Local variables ------------------------------ 3237 !scalars 3238 integer :: ptr,ierr,band 3239 character(len=500) :: msg 3240 !arrays 3241 real(dp) :: norm(nband),alpha(2) 3242 3243 ! ************************************************************************* 3244 3245 !$OMP PARALLEL DO PRIVATE(ptr) IF (nband > 1) 3246 do band=1,nband 3247 ptr = 1 + 2*npws*(band-1) 3248 norm(band) = cg_real_zdotc(npws,gsc(ptr),cg(ptr)) 3249 end do 3250 3251 if (istwfk>1) then 3252 norm = two * norm 3253 if (istwfk==2 .and. me_g0==1) then 3254 !$OMP PARALLEL DO PRIVATE(ptr) IF (nband > 1) 3255 do band=1,nband 3256 ptr = 1 + 2*npws*(band-1) 3257 norm(band) = norm(band) - gsc(ptr) * cg(ptr) 3258 end do 3259 end if 3260 end if 3261 3262 if (comm_pw /= xmpi_comm_self) then 3263 call xmpi_sum(norm,comm_pw,ierr) 3264 end if 3265 3266 ierr = 0 3267 do band=1,nband 3268 if (norm(band) > zero) then 3269 norm(band) = SQRT(norm(band)) 3270 else 3271 ierr = ierr + 1 3272 end if 3273 end do 3274 3275 if (ierr/=0) then 3276 write(msg,'(a,i0,a)')" Found ",ierr," vectors with norm <= zero!" 3277 MSG_ERROR(msg) 3278 end if 3279 3280 ! Scale |C> and S|C>. 3281 !$OMP PARALLEL DO PRIVATE(ptr,alpha) IF (nband > 1) 3282 do band=1,nband 3283 ptr = 1 + 2*npws*(band-1) 3284 alpha = (/one/norm(band), zero/) 3285 call cg_zscal(npws,alpha,cg(ptr)) 3286 call cg_zscal(npws,alpha,gsc(ptr)) 3287 end do 3288 3289 end subroutine cgpaw_normalize
m_cgtools/cz_zprecon_block [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
cg_zprecon_block
FUNCTION
precondition $<G|(H-e_{n,k})|C_{n,k}>$ for a block of band (band-FFT parallelisation)
INPUTS
blocksize= size of blocks of bands
$cg(vectsize,blocksize)=<G|C_{n,k}> for a block of bands$.
$eval(blocksize,blocksize)=current block of bands eigenvalues=<C_{n,k}|H|C_{n,k}>$.
$ghc(vectsize,blocksize)=<G|H|C_{n,k}> for a block of bands$.
iterationnumber=number of iterative minimizations in LOBPCG
kinpw(npw)=(modified) kinetic energy for each plane wave (Hartree)
nspinor=number of spinorial components of the wavefunctions (on current proc)
$vect(vectsize,blocksize)=<G|H|C_{n,k}> for a block of bands$.
npw=number of planewaves at this k point.
optekin= 1 if the kinetic energy used in preconditionning is modified
according to Kresse, Furthmuller, PRB 54, 11169 (1996)
0 otherwise
optpcon= 0 the TPA preconditionning matrix does not depend on band
1 the TPA preconditionning matrix (not modified)
2 the TPA preconditionning matrix is independant of iterationnumber
vectsize= size of vectors
comm=MPI communicator.
OUTPUT
vect(2,npw)=<g|(h-eval)|c_{n,k}>*(polynomial ratio)
SIDE EFFECTS
pcon(npw,blocksize)=preconditionning matrix
input if optpcon=0,2 and iterationnumber/=1
output if optpcon=0,2 and iterationnumber==1
PARENTS
m_lobpcg
CHILDREN
SOURCE
4327 subroutine cg_zprecon_block(cg,eval,blocksize,iterationnumber,kinpw,& 4328 & npw,nspinor,optekin,optpcon,pcon,ghc,vect,vectsize,comm) 4329 4330 4331 !This section has been created automatically by the script Abilint (TD). 4332 !Do not modify the following lines by hand. 4333 #undef ABI_FUNC 4334 #define ABI_FUNC 'cg_zprecon_block' 4335 use interfaces_18_timing 4336 !End of the abilint section 4337 4338 implicit none 4339 4340 !Arguments ------------------------------------ 4341 !scalars 4342 integer,intent(in) :: blocksize,iterationnumber,npw,nspinor,optekin 4343 integer,intent(in) :: optpcon,vectsize,comm 4344 !arrays 4345 real(dp),intent(in) :: kinpw(npw) 4346 real(dp),intent(inout) :: pcon(npw,blocksize) 4347 complex(dpc),intent(in) :: cg(vectsize,blocksize),eval(blocksize,blocksize) 4348 complex(dpc),intent(in) :: ghc(vectsize,blocksize) 4349 complex(dpc),intent(inout) :: vect(vectsize,blocksize) 4350 4351 !Local variables------------------------------- 4352 !scalars 4353 integer :: iblocksize,ierr,ig,igs,ispinor 4354 real(dp) :: fac,poly,xx 4355 character(len=500) :: message 4356 !arrays 4357 real(dp) :: tsec(2) 4358 real(dp),allocatable :: ek0(:),ek0_inv(:) 4359 4360 ! ************************************************************************* 4361 4362 call timab(536,1,tsec) 4363 4364 !In this case, the Teter, Allan and Payne preconditioner is approximated: 4365 !the factor xx=Ekin(G) and no more Ekin(G)/Ekin(iband) 4366 if (optpcon==0) then 4367 do ispinor=1,nspinor 4368 igs=(ispinor-1)*npw 4369 do ig=1+igs,npw+igs 4370 if (iterationnumber==1) then 4371 if(kinpw(ig-igs)<huge(0.0_dp)*1.d-11)then 4372 xx=kinpw(ig-igs) 4373 ! teter polynomial ratio 4374 poly=27._dp+xx*(18._dp+xx*(12._dp+xx*8._dp)) 4375 fac=poly/(poly+16._dp*xx**4) 4376 if (optekin==1) fac=two*fac 4377 pcon(ig-igs,1)=fac 4378 do iblocksize=1,blocksize 4379 vect(ig,iblocksize)=(ghc(ig,iblocksize)-eval(iblocksize,iblocksize)*cg(ig,iblocksize))*pcon(ig-igs,1) 4380 end do 4381 else 4382 pcon(ig-igs,1)=zero 4383 vect(ig,:)=dcmplx(0.0_dp,0.0_dp) 4384 end if 4385 else 4386 do iblocksize=1,blocksize 4387 vect(ig,iblocksize)=(ghc(ig,iblocksize)-eval(iblocksize,iblocksize)*cg(ig,iblocksize))*pcon(ig-igs,1) 4388 end do 4389 end if 4390 end do 4391 end do 4392 4393 else if (optpcon>0) then 4394 ! Compute mean kinetic energy of all bands 4395 ABI_ALLOCATE(ek0,(blocksize)) 4396 ABI_ALLOCATE(ek0_inv,(blocksize)) 4397 if (iterationnumber==1.or.optpcon==1) then 4398 do iblocksize=1,blocksize 4399 ek0(iblocksize)=0.0_dp 4400 do ispinor=1,nspinor 4401 igs=(ispinor-1)*npw 4402 do ig=1+igs,npw+igs 4403 if(kinpw(ig-igs)<huge(0.0_dp)*1.d-11)then 4404 ek0(iblocksize)=ek0(iblocksize)+kinpw(ig-igs)*& 4405 & (real(cg(ig,iblocksize))**2+aimag(cg(ig,iblocksize))**2) 4406 end if 4407 end do 4408 end do 4409 end do 4410 4411 call xmpi_sum(ek0,comm,ierr) 4412 4413 do iblocksize=1,blocksize 4414 if(ek0(iblocksize)<1.0d-10)then 4415 message = 'the mean kinetic energy of a wavefunction vanishes. it is reset to 0.1ha.' 4416 MSG_WARNING(message) 4417 ek0(iblocksize)=0.1_dp 4418 end if 4419 end do 4420 if (optekin==1) then 4421 ek0_inv(:)=2.0_dp/(3._dp*ek0(:)) 4422 else 4423 ek0_inv(:)=1.0_dp/ek0(:) 4424 end if 4425 end if !iterationnumber==1.or.optpcon==1 4426 4427 ! Carry out preconditioning 4428 do iblocksize=1,blocksize 4429 do ispinor=1,nspinor 4430 igs=(ispinor-1)*npw 4431 do ig=1+igs,npw+igs 4432 if (iterationnumber==1.or.optpcon==1) then 4433 if(kinpw(ig-igs)<huge(0.0_dp)*1.d-11)then 4434 xx=kinpw(ig-igs)*ek0_inv(iblocksize) 4435 ! teter polynomial ratio 4436 poly=27._dp+xx*(18._dp+xx*(12._dp+xx*8._dp)) 4437 fac=poly/(poly+16._dp*xx**4) 4438 if (optekin==1) fac=two*fac 4439 pcon(ig-igs,iblocksize)=fac 4440 vect(ig,iblocksize)=(ghc(ig,iblocksize)-& 4441 & eval(iblocksize,iblocksize)*cg(ig,iblocksize))*pcon(ig-igs,iblocksize) 4442 else 4443 pcon(ig-igs,iblocksize)=zero 4444 vect(ig,iblocksize)=dcmplx(0.0_dp,0.0_dp) 4445 end if 4446 else 4447 vect(ig,iblocksize)=(ghc(ig,iblocksize)-& 4448 & eval(iblocksize,iblocksize)*cg(ig,iblocksize))*pcon(ig-igs,iblocksize) 4449 end if 4450 end do 4451 end do 4452 end do 4453 ABI_DEALLOCATE(ek0) 4454 ABI_DEALLOCATE(ek0_inv) 4455 end if !optpcon 4456 4457 call timab(536,2,tsec) 4458 4459 end subroutine cg_zprecon_block
m_cgtools/dotprod_g [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
dotprod_g
FUNCTION
Compute scalar product <vec1|vect2> of complex vectors vect1 and vect2 (can be the same) Take into account the storage mode of the vectors (istwf_k) If option=1, compute only real part, if option=2 compute also imaginary part. If the number of calls to the dot product scales quadratically with the volume of the system, it is preferable not to call the present routine, but but to write a specially optimized routine, that will avoid many branches related to the existence of 'option' and 'istwf_k'.
INPUTS
istwf_k=option parameter that describes the storage of wfs
vect1(2,npw)=first vector (one should take its complex conjugate)
vect2(2,npw)=second vector
npw= (effective) number of planewaves at this k point (including spinorial level)
option= 1 if only real part to be computed,
2 if both real and imaginary.
3 if in case istwf_k==1 must compute real and imaginary parts,
but if istwf_k >1 must compute only real part
me_g0=1 if this processor treats G=0, 0 otherwise
comm=MPI communicator used to reduce the results.
OUTPUT
$doti=\Im ( <vect1|vect2> )$ , output only if option=2 and eventually option=3. $dotr=\Re ( <vect1|vect2> )$
PARENTS
cgwf,chebfi,corrmetalwf1,d2frnl,dfpt_cgwf,dfpt_nsteltwf,dfpt_nstpaw dfpt_nstwf,dfpt_vtowfk,dfpt_wfkfermi,dfptnl_resp,dotprod_set_cgcprj dotprodm_sumdiag_cgcprj,eig2stern,extrapwf,fock2ACE,fock_ACE_getghc fock_getghc,m_efmas,m_gkk,m_phgamma,m_phpi,m_rf2,m_sigmaph,mkresi nonlop_gpu,nonlop_test,rf2_init
CHILDREN
SOURCE
1304 subroutine dotprod_g(dotr,doti,istwf_k,npw,option,vect1,vect2,me_g0,comm) 1305 1306 1307 !This section has been created automatically by the script Abilint (TD). 1308 !Do not modify the following lines by hand. 1309 #undef ABI_FUNC 1310 #define ABI_FUNC 'dotprod_g' 1311 !End of the abilint section 1312 1313 implicit none 1314 1315 !Arguments ------------------------------------ 1316 !scalars 1317 integer,intent(in) :: istwf_k,npw,option,me_g0,comm 1318 real(dp),intent(out) :: doti,dotr 1319 !arrays 1320 real(dp),intent(in) :: vect1(2,npw),vect2(2,npw) 1321 1322 !Local variables------------------------------- 1323 !scalars 1324 integer :: ierr 1325 !arrays 1326 real(dp) :: dotarr(2) 1327 1328 ! ************************************************************************* 1329 1330 if (istwf_k==1) then ! General k-point 1331 1332 if(option==1)then 1333 dotr = cg_real_zdotc(npw,vect1,vect2) 1334 else 1335 dotarr = cg_zdotc(npw,vect1,vect2) 1336 dotr = dotarr(1) 1337 doti = dotarr(2) 1338 end if 1339 1340 else if (istwf_k==2 .and. me_g0==1) then ! Gamma k-point and I have G=0 1341 1342 dotr=half*vect1(1,1)*vect2(1,1) 1343 dotr = dotr + cg_real_zdotc(npw-1,vect1(1,2),vect2(1,2)) 1344 dotr = two*dotr 1345 if (option==2) doti=zero 1346 1347 else ! Other TR k-points 1348 1349 dotr = cg_real_zdotc(npw,vect1,vect2) 1350 dotr=two*dotr 1351 if (option==2) doti=zero 1352 end if 1353 1354 !Reduction in case of parallelism 1355 if (xmpi_comm_size(comm)>1) then 1356 if (option==1.or.istwf_k/=1) then 1357 call xmpi_sum(dotr,comm,ierr) 1358 else 1359 dotarr(1)=dotr ; dotarr(2)=doti 1360 call xmpi_sum(dotarr,comm,ierr) 1361 dotr=dotarr(1) ; doti=dotarr(2) 1362 end if 1363 end if 1364 1365 end subroutine dotprod_g
m_cgtools/dotprod_v [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
dotprod_v
FUNCTION
Compute dot product of two potentials (integral over FFT grid), to obtain a square residual-like quantity (so the sum of product of values is NOT divided by the number of FFT points, and NOT multiplied by the primitive cell volume). Take into account the spin components of the potentials (nspden), and sum over them.
INPUTS
cplex=if 1, real space functions on FFT grid are REAL, if 2, COMPLEX nfft= (effective) number of FFT grid points (for this processor) nspden=number of spin-density components opt_storage: 0, if potentials are stored as V^up-up, V^dn-dn, Re[V^up-dn], Im[V^up-dn] 1, if potentials are stored as V, B_x, B_y, Bz (B=magn. field) pot1(cplex*nfft,nspden)=first real space potential on FFT grid pot2(cplex*nfft,nspden)=second real space potential on FFT grid comm=MPI communicator in which results will be reduced.
OUTPUT
dotr= value of the dot product
PARENTS
m_epjdos
CHILDREN
SOURCE
1567 subroutine dotprod_v(cplex,dotr,nfft,nspden,opt_storage,pot1,pot2,comm) 1568 1569 1570 !This section has been created automatically by the script Abilint (TD). 1571 !Do not modify the following lines by hand. 1572 #undef ABI_FUNC 1573 #define ABI_FUNC 'dotprod_v' 1574 !End of the abilint section 1575 1576 implicit none 1577 1578 !Arguments ------------------------------------ 1579 !scalars 1580 integer,intent(in) :: cplex,nfft,nspden,opt_storage,comm 1581 real(dp),intent(out) :: dotr 1582 !arrays 1583 real(dp),intent(in) :: pot1(cplex*nfft,nspden),pot2(cplex*nfft,nspden) 1584 1585 !Local variables------------------------------- 1586 !scalars 1587 integer :: ierr,ifft,ispden 1588 real(dp) :: ar 1589 !arrays 1590 1591 ! ************************************************************************* 1592 1593 !Real or complex inputs are coded 1594 1595 dotr=zero 1596 !$OMP PARALLEL DO COLLAPSE(2) REDUCTION(+:dotr) 1597 do ispden=1,min(nspden,2) 1598 do ifft=1,cplex*nfft 1599 dotr =dotr + pot1(ifft,ispden)*pot2(ifft,ispden) 1600 end do 1601 end do 1602 1603 if (nspden==4) then 1604 ar=zero 1605 !$OMP PARALLEL DO COLLAPSE(2) REDUCTION(+:ar) 1606 do ispden=3,4 1607 do ifft=1,cplex*nfft 1608 ar = ar + pot1(ifft,ispden)*pot2(ifft,ispden) 1609 end do 1610 end do 1611 1612 if (opt_storage==0) then 1613 if (cplex==1) then 1614 dotr = dotr+two*ar 1615 else 1616 dotr = dotr+ar 1617 end if 1618 else 1619 dotr = half*(dotr+ar) 1620 end if 1621 end if 1622 1623 !MPIWF reduction (addition) on dotr is needed here 1624 if (xmpi_comm_size(comm)>1) then 1625 call xmpi_sum(dotr,comm,ierr) 1626 end if 1627 1628 end subroutine dotprod_v
m_cgtools/dotprod_vn [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
dotprod_vn
FUNCTION
Compute dot product of potential and density (integral over FFT grid), to obtain an energy-like quantity (so the usual dotproduct is divided by the number of FFT points, and multiplied by the primitive cell volume). Take into account the spin components of the density and potentials (nspden), and sum correctly over them. Note that the storage of densities and potentials is different : for potential, one stores the matrix components, while for the density, one stores the trace, and then, either the up component (if nspden=2), or the magnetization vector (if nspden=4).
INPUTS
cplex=if 1, real space functions on FFT grid are REAL, if 2, COMPLEX dens(cplex*nfft,nspden)=real space density on FFT grid mpi_enreg=informations about MPI parallelization nfft= (effective) number of FFT grid points (for this processor) nfftot= total number of FFT grid points nspden=number of spin-density components option= if 1, only the real part is computed if 2, both real and imaginary parts are computed (not yet coded) pot(cplex*nfft,nspden)=real space potential on FFT grid (will be complex conjugated if cplex=2 and option=2) ucvol=unit cell volume (Bohr**3)
OUTPUT
doti= imaginary part of the dot product, output only if option=2 (and cplex=2). dotr= real part
NOTES
Concerning storage when nspden=4:
cplex=1:
V is stored as : V^11, V^22, Re[V^12], Im[V^12] (complex, hermitian)
N is stored as : n, m_x, m_y, m_z (real)
cplex=2:
V is stored as : V^11, V^22, V^12, i.V^21 (complex)
N is stored as : n, m_x, m_y, mz (complex)
PARENTS
dfpt_dyxc1,dfpt_eltfrxc,dfpt_nselt,dfpt_nstdy,dfpt_nstpaw,dfpt_rhotov dfptnl_loop,energy,newfermie1,odamix,prcrskerker2,rhotov,rhotoxc,setvtr
CHILDREN
xmpi_sum
SOURCE
1680 subroutine dotprod_vn(cplex,dens,dotr,doti,nfft,nfftot,nspden,option,pot,ucvol, & 1681 mpi_comm_sphgrid) ! Optional 1682 1683 1684 !This section has been created automatically by the script Abilint (TD). 1685 !Do not modify the following lines by hand. 1686 #undef ABI_FUNC 1687 #define ABI_FUNC 'dotprod_vn' 1688 !End of the abilint section 1689 1690 implicit none 1691 1692 !Arguments ------------------------------------ 1693 !scalars 1694 integer,intent(in) :: cplex,nfft,nfftot,nspden,option 1695 integer,intent(in),optional :: mpi_comm_sphgrid 1696 real(dp),intent(in) :: ucvol 1697 real(dp),intent(out) :: doti,dotr 1698 !arrays 1699 real(dp),intent(in) :: dens(cplex*nfft,nspden),pot(cplex*nfft,nspden) 1700 1701 !Local variables------------------------------- 1702 !scalars 1703 integer :: ierr,ifft,jfft 1704 real(dp) :: dim11,dim12,dim21,dim22,dim_dn,dim_up,dre11,dre12,dre21,dre22 1705 real(dp) :: dre_dn,dre_up,factor,nproc_sphgrid,pim11,pim12,pim21,pim22,pim_dn,pim_up,pre11 1706 real(dp) :: pre12,pre21,pre22,pre_dn,pre_up 1707 real(dp) :: bx_re,bx_im,by_re,by_im,bz_re,bz_im,v0_re,v0_im 1708 !arrays 1709 real(dp) :: buffer2(2) 1710 1711 ! ************************************************************************* 1712 1713 !Real or complex inputs are coded 1714 DBG_CHECK(ANY(cplex==(/1,2/)),"Wrong cplex") 1715 DBG_CHECK(ANY(nspden==(/1,2,4/)),"Wrong nspden") 1716 1717 !Real or complex output are coded 1718 DBG_CHECK(ANY(option==(/1,2/)),"Wrong option") 1719 1720 dotr=zero; doti=zero 1721 1722 if(nspden==1)then 1723 1724 if(option==1 .or. cplex==1 )then 1725 !$OMP PARALLEL DO REDUCTION(+:dotr) 1726 do ifft=1,cplex*nfft 1727 dotr=dotr + pot(ifft,1)*dens(ifft,1) 1728 end do 1729 ! dotr = ddot(cplex*nfft,pot,1,dens,1) 1730 1731 else ! option==2 and cplex==2 : one builds the imaginary part, from complex den/pot 1732 1733 !$OMP PARALLEL DO PRIVATE(jfft) REDUCTION(+:dotr,doti) 1734 do ifft=1,nfft 1735 jfft=2*ifft 1736 dotr=dotr + pot(jfft-1,1)*dens(jfft-1,1) + pot(jfft,1)*dens(jfft ,1) 1737 doti=doti + pot(jfft-1,1)*dens(jfft ,1) - pot(jfft,1)*dens(jfft-1,1) 1738 end do 1739 1740 end if 1741 1742 else if(nspden==2)then 1743 1744 if(option==1 .or. cplex==1 )then 1745 !$OMP PARALLEL DO REDUCTION(+:dotr) 1746 do ifft=1,cplex*nfft 1747 dotr=dotr + pot(ifft,1)* dens(ifft,2) & ! This is the spin up contribution 1748 & + pot(ifft,2)*(dens(ifft,1)-dens(ifft,2)) ! This is the spin down contribution 1749 end do 1750 1751 else ! option==2 and cplex==2 : one builds the imaginary part, from complex den/pot 1752 1753 !$OMP PARALLEL DO DEFAULT(PRIVATE) SHARED(nfft,dens,pot) REDUCTION(+:dotr,doti) 1754 do ifft=1,nfft 1755 1756 jfft=2*ifft 1757 dre_up=dens(jfft-1,2) 1758 dim_up=dens(jfft ,2) 1759 dre_dn=dens(jfft-1,1)-dre_up 1760 dim_dn=dens(jfft ,1)-dim_up 1761 pre_up=pot(jfft-1,1) 1762 pim_up=pot(jfft ,1) 1763 pre_dn=pot(jfft-1,2) 1764 pim_dn=pot(jfft ,2) 1765 1766 dotr=dotr + pre_up * dre_up & 1767 & + pim_up * dim_up & 1768 & + pre_dn * dre_dn & 1769 & + pim_dn * dim_dn 1770 doti=doti + pre_up * dim_up & 1771 & - pim_up * dre_up & 1772 & + pre_dn * dim_dn & 1773 & - pim_dn * dre_dn 1774 1775 end do 1776 end if 1777 1778 else if(nspden==4)then 1779 ! \rho{\alpha,\beta} V^{\alpha,\beta} = 1780 ! rho*(V^{11}+V^{22})/2$ 1781 ! + m_x Re(V^{12})- m_y Im{V^{12}}+ m_z(V^{11}-V^{22})/2 1782 if (cplex==1) then 1783 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(nfft,dens,pot) REDUCTION(+:dotr) 1784 do ifft=1,nfft 1785 dotr=dotr + & 1786 & (pot(ifft,1) + pot(ifft,2))*half*dens(ifft,1) & ! This is the density contrib 1787 & + pot(ifft,3) *dens(ifft,2) & ! This is the m_x contrib 1788 & - pot(ifft,4) *dens(ifft,3) & ! This is the m_y contrib 1789 & +(pot(ifft,1) - pot(ifft,2))*half*dens(ifft,4) ! This is the m_z contrib 1790 end do 1791 else ! cplex=2 1792 ! Note concerning storage when cplex=2: 1793 ! V is stored as : v^11, v^22, V^12, i.V^21 (each are complex) 1794 ! N is stored as : n, m_x, m_y, mZ (each are complex) 1795 if (option==1) then 1796 !$OMP PARALLEL DO DEFAULT(PRIVATE) SHARED(nfft,dens,pot) REDUCTION(+:dotr) 1797 do ifft=1,nfft 1798 jfft=2*ifft 1799 dre11=half*(dens(jfft-1,1)+dens(jfft-1,4)) 1800 dim11=half*(dens(jfft ,1)+dens(jfft-1,4)) 1801 dre22=half*(dens(jfft-1,1)-dens(jfft-1,4)) 1802 dim22=half*(dens(jfft ,1)-dens(jfft-1,4)) 1803 dre12=half*(dens(jfft-1,2)+dens(jfft ,3)) 1804 dim12=half*(dens(jfft ,2)-dens(jfft-1,3)) 1805 dre21=half*(dens(jfft-1,2)-dens(jfft ,3)) 1806 dim21=half*(dens(jfft ,2)+dens(jfft-1,3)) 1807 pre11= pot(jfft-1,1) 1808 pim11= pot(jfft ,1) 1809 pre22= pot(jfft-1,2) 1810 pim22= pot(jfft ,2) 1811 pre12= pot(jfft-1,3) 1812 pim12= pot(jfft ,3) 1813 pre21= pot(jfft ,4) 1814 pim21=-pot(jfft-1,4) 1815 1816 v0_re=half*(pre11+pre22) 1817 v0_im=half*(pim11+pim22) 1818 bx_re=half*(pre12+pre21) 1819 bx_im=half*(pim12+pim21) 1820 by_re=half*(-pim12+pim21) 1821 by_im=half*(pre12-pre21) 1822 bz_re=half*(pre11-pre22) 1823 bz_im=half*(pim11-pim22) 1824 dotr=dotr+v0_re * dens(jfft-1,1)& 1825 & + v0_im * dens(jfft ,1) & 1826 & + bx_re * dens(jfft-1,2) & 1827 & + bx_im * dens(jfft ,2) & 1828 & + by_re * dens(jfft-1,3) & 1829 & + by_im * dens(jfft ,3) & 1830 & + bz_re * dens(jfft-1,4) & 1831 & + bz_im * dens(jfft ,4) 1832 ! dotr=dotr + pre11 * dre11 & 1833 !& + pim11 * dim11 & 1834 !& + pre22 * dre22 & 1835 !& + pim22 * dim22 & 1836 !& + pre12 * dre12 & 1837 !& + pim12 * dim12 & 1838 !& + pre21 * dre21 & 1839 !& + pim21 * dim21 1840 end do 1841 else ! option=2 1842 !$OMP PARALLEL DO DEFAULT(PRIVATE) SHARED(nfft,dens,pot) REDUCTION(+:dotr,doti) 1843 do ifft=1,nfft 1844 jfft=2*ifft 1845 dre11=half*(dens(jfft-1,1)+dens(jfft-1,4)) 1846 dim11=half*(dens(jfft ,1)+dens(jfft-1,4)) 1847 dre22=half*(dens(jfft-1,1)-dens(jfft-1,4)) 1848 dim22=half*(dens(jfft ,1)-dens(jfft-1,4)) 1849 dre12=half*(dens(jfft-1,2)+dens(jfft ,3)) 1850 dim12=half*(dens(jfft ,2)-dens(jfft-1,3)) 1851 dre21=half*(dens(jfft-1,2)-dens(jfft ,3)) 1852 dim21=half*(dens(jfft ,2)+dens(jfft-1,3)) 1853 pre11= pot(jfft-1,1) 1854 pim11= pot(jfft ,1) 1855 pre22= pot(jfft-1,2) 1856 pim22= pot(jfft ,2) 1857 pre12= pot(jfft-1,3) 1858 pim12= pot(jfft ,3) 1859 pre21= pot(jfft ,4) 1860 pim21=-pot(jfft-1,4) 1861 dotr=dotr + pre11 * dre11 & 1862 & + pim11 * dim11 & 1863 & + pre22 * dre22 & 1864 & + pim22 * dim22 & 1865 & + pre12 * dre12 & 1866 & + pim12 * dim12 & 1867 & + pre21 * dre21 & 1868 & + pim21 * dim21 1869 doti=doti + pre11 * dim11 & 1870 & - pim11 * dre11 & 1871 & + pre22 * dim22 & 1872 & - pim22 * dre22 & 1873 & + pre12 * dim12 & 1874 & - pim12 * dre12 & 1875 & + pre21 * dim21 & 1876 & - pim21 * dre21 1877 end do 1878 end if ! option 1879 end if ! cplex 1880 end if ! nspden 1881 1882 factor=ucvol/dble(nfftot) 1883 dotr=factor*dotr 1884 doti=factor*doti 1885 1886 !MPIWF reduction (addition) on dotr, doti is needed here 1887 if(present(mpi_comm_sphgrid)) then 1888 nproc_sphgrid=xmpi_comm_size(mpi_comm_sphgrid) 1889 if(nproc_sphgrid>1) then 1890 buffer2(1)=dotr 1891 buffer2(2)=doti 1892 call xmpi_sum(buffer2,mpi_comm_sphgrid,ierr) 1893 dotr=buffer2(1) 1894 doti=buffer2(2) 1895 end if 1896 end if 1897 1898 end subroutine dotprod_vn
m_cgtools/fxphas_seq [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
fxphas_seq
FUNCTION
Fix phase of all bands. Keep normalization but maximize real part (minimize imag part). Also fix the sign of real part by setting the first non-zero element to be positive. This version has been stripped of all the mpi_enreg junk by MJV
INPUTS
cg(2,mcg)= contains the wavefunction |c> coefficients.
gsc(2,mgsc)= if useoverlap==1, contains the S|c> coefficients,
where S is an overlap matrix.
icg=shift to be applied on the location of data in the array cg
igsc=shift to be applied on the location of data in the array gsc
istwfk=input option parameter that describes the storage of wfs
(set to 1 if usual complex vectors)
mcg=size of second dimension of cg
mgsc=size of second dimension of gsc
nband_k=number of bands
npw_k=number of planewaves
useoverlap=describe the overlap of wavefunctions:
0: no overlap (S=Identi0,ty_matrix)
1: wavefunctions are overlapping
OUTPUT
cg(2,mcg)=same array with altered phase. gsc(2,mgsc)= same array with altered phase.
NOTES
When the sign of the real part was fixed (modif v3.1.3g.6), the test Tv3#5 , dataset 5, behaved differently than previously. This should be cleared up.
PARENTS
m_dynmat,rayleigh_ritz
CHILDREN
SOURCE
4506 subroutine fxphas_seq(cg,gsc,icg,igsc,istwfk,mcg,mgsc,nband_k,npw_k,useoverlap) 4507 4508 4509 !This section has been created automatically by the script Abilint (TD). 4510 !Do not modify the following lines by hand. 4511 #undef ABI_FUNC 4512 #define ABI_FUNC 'fxphas_seq' 4513 !End of the abilint section 4514 4515 implicit none 4516 4517 !Arguments ------------------------------------ 4518 !scalars 4519 integer,intent(in) :: icg,igsc,istwfk,mcg,mgsc,nband_k,npw_k,useoverlap 4520 !arrays 4521 real(dp),intent(inout) :: cg(2,mcg),gsc(2,mgsc*useoverlap) 4522 4523 !Local variables------------------------------- 4524 !scalars 4525 integer :: iband,ii,indx 4526 real(dp) :: cim,cre,gscim,gscre,quotient,root1,root2,saa,sab,sbb,theta 4527 real(dp) :: thppi,xx,yy 4528 character(len=500) :: message 4529 !arrays 4530 real(dp),allocatable :: cimb(:),creb(:),saab(:),sabb(:),sbbb(:) !,sarr(:,:) 4531 4532 ! ************************************************************************* 4533 4534 !The general case, where a complex phase indeterminacy is present 4535 if(istwfk==1)then 4536 4537 ABI_ALLOCATE(cimb,(nband_k)) 4538 ABI_ALLOCATE(creb,(nband_k)) 4539 ABI_ALLOCATE(saab,(nband_k)) 4540 ABI_ALLOCATE(sabb,(nband_k)) 4541 ABI_ALLOCATE(sbbb,(nband_k)) 4542 cimb(:)=zero ; creb(:)=zero 4543 4544 ! Loop over bands 4545 ! TODO: MG store saa arrays in sarr(3,nband_k) to reduce false sharing. 4546 do iband=1,nband_k 4547 indx=icg+(iband-1)*npw_k 4548 4549 ! Compute several sums over Re, Im parts of c 4550 saa=0.0_dp ; sbb=0.0_dp ; sab=0.0_dp 4551 do ii=1+indx,npw_k+indx 4552 saa=saa+cg(1,ii)*cg(1,ii) 4553 sbb=sbb+cg(2,ii)*cg(2,ii) 4554 sab=sab+cg(1,ii)*cg(2,ii) 4555 end do 4556 saab(iband)=saa 4557 sbbb(iband)=sbb 4558 sabb(iband)=sab 4559 end do ! iband 4560 4561 4562 do iband=1,nband_k 4563 4564 indx=icg+(iband-1)*npw_k 4565 4566 saa=saab(iband) 4567 sbb=sbbb(iband) 4568 sab=sabb(iband) 4569 4570 ! Get phase angle theta 4571 if (sbb+saa>tol8)then 4572 if(abs(sbb-saa)>tol8*(sbb+saa) .or. 2*abs(sab)>tol8*(sbb+saa))then 4573 if (abs(sbb-saa)>tol8*abs(sab)) then 4574 quotient=sab/(sbb-saa) 4575 theta=0.5_dp*atan(2.0_dp*quotient) 4576 else 4577 ! Taylor expansion of the atan in terms of inverse of its argument. Correct up to 1/x2, included. 4578 theta=0.25_dp*(pi-(sbb-saa)/sab) 4579 end if 4580 ! Check roots to get theta for max Re part 4581 root1=cos(theta)**2*saa+sin(theta)**2*sbb-2.0_dp*cos(theta)*sin(theta)*sab 4582 thppi=theta+0.5_dp*pi 4583 root2=cos(thppi)**2*saa+sin(thppi)**2*sbb-2.0_dp*cos(thppi)*sin(thppi)*sab 4584 if (root2>root1) theta=thppi 4585 else 4586 ! The real part vector and the imaginary part vector are orthogonal, and of same norm. Strong indeterminacy. 4587 ! Will determine the first non-zero coefficient, and fix its phase 4588 do ii=1+indx,npw_k+indx 4589 cre=cg(1,ii) 4590 cim=cg(2,ii) 4591 if(cre**2+cim**2>tol8**2*(saa+sbb))then 4592 if(cre**2>tol8**2**cim**2)then 4593 theta=atan(cim/cre) 4594 else 4595 ! Taylor expansion of the atan in terms of inverse of its argument. Correct up to 1/x2, included. 4596 theta=pi/2-cre/cim 4597 end if 4598 exit 4599 end if 4600 end do 4601 end if 4602 else 4603 write(message,'(a,i0,5a)')& 4604 & 'The eigenvector with band ',iband,' has zero norm.',ch10,& 4605 & 'This usually happens when the number of bands (nband) is comparable to the number of planewaves (mpw)',ch10,& 4606 & 'Action: Check the parameters of the calculation. If nband ~ mpw, then decrease nband or, alternatively, increase ecut' 4607 MSG_ERROR(message) 4608 end if 4609 4610 xx=cos(theta) 4611 yy=sin(theta) 4612 4613 ! Here, set the first non-zero element to be positive 4614 do ii=1+indx,npw_k+indx 4615 cre=cg(1,ii) 4616 cim=cg(2,ii) 4617 cre=xx*cre-yy*cim 4618 if(abs(cre)>tol8)exit 4619 end do 4620 if(cre<zero)then 4621 xx=-xx ; yy=-yy 4622 end if 4623 4624 creb(iband)=xx 4625 cimb(iband)=yy 4626 4627 end do 4628 4629 do iband=1,nband_k 4630 4631 indx=icg+(iband-1)*npw_k 4632 4633 xx=creb(iband) 4634 yy=cimb(iband) 4635 do ii=1+indx,npw_k+indx 4636 cre=cg(1,ii) 4637 cim=cg(2,ii) 4638 cg(1,ii)=xx*cre-yy*cim 4639 cg(2,ii)=xx*cim+yy*cre 4640 end do 4641 4642 ! Alter phase of array S|cg> 4643 if (useoverlap==1) then 4644 indx=igsc+(iband-1)*npw_k 4645 do ii=1+indx,npw_k+indx 4646 gscre=gsc(1,ii) 4647 gscim=gsc(2,ii) 4648 gsc(1,ii)=xx*gscre-yy*gscim 4649 gsc(2,ii)=xx*gscim+yy*gscre 4650 end do 4651 end if 4652 4653 end do ! iband 4654 4655 ABI_DEALLOCATE(cimb) 4656 ABI_DEALLOCATE(creb) 4657 ABI_DEALLOCATE(saab) 4658 ABI_DEALLOCATE(sabb) 4659 ABI_DEALLOCATE(sbbb) 4660 4661 ! ==================================================================== 4662 4663 ! Storages that take into account the time-reversal symmetry : the freedom is only a sign freedom 4664 else ! if istwfk/=1 4665 4666 ABI_ALLOCATE(creb,(nband_k)) 4667 creb(:)=zero 4668 ! Loop over bands 4669 do iband=1,nband_k 4670 4671 indx=icg+(iband-1)*npw_k 4672 4673 ! Here, set the first non-zero real element to be positive 4674 do ii=1+indx,npw_k+indx 4675 cre=cg(1,ii) 4676 if(abs(cre)>tol8)exit 4677 end do 4678 creb(iband)=cre 4679 4680 end do ! iband 4681 4682 do iband=1,nband_k 4683 4684 cre=creb(iband) 4685 if(cre<zero)then 4686 indx=icg+(iband-1)*npw_k 4687 do ii=1+indx,npw_k+indx 4688 cg(1,ii)=-cg(1,ii) 4689 cg(2,ii)=-cg(2,ii) 4690 end do 4691 if(useoverlap==1)then 4692 indx=igsc+(iband-1)*npw_k 4693 do ii=1+indx,npw_k+indx 4694 gsc(1,ii)=-gsc(1,ii) 4695 gsc(2,ii)=-gsc(2,ii) 4696 end do 4697 end if 4698 end if 4699 4700 end do ! iband 4701 4702 ABI_DEALLOCATE(creb) 4703 4704 end if ! istwfk 4705 4706 end subroutine fxphas_seq
m_cgtools/matrixelmt_g [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
matrixelmt_g
FUNCTION
Compute a matrix element of two wavefunctions, in reciprocal space, for an operator that is diagonal in reciprocal space: <wf1|op|wf2> For the time being, only spin-independent operators are treated.
INPUTS
diag(npw)=diagonal operator (real, spin-independent !) istwf_k=storage mode of the vectors needimag=0 if the imaginary part is not needed ; 1 if the imaginary part is needed npw=number of planewaves of the first vector nspinor=number of spinor components vect1(2,npw*nspinor)=first vector vect2(2,npw*nspinor)=second vector comm_fft=MPI communicator for the FFT me_g0=1 if this processors treats the G=0 component.
OUTPUT
ai=imaginary part of the matrix element ar=real part of the matrix element
PARENTS
dfpt_vtowfk
CHILDREN
SOURCE
1401 subroutine matrixelmt_g(ai,ar,diag,istwf_k,needimag,npw,nspinor,vect1,vect2,me_g0,comm_fft) 1402 1403 1404 !This section has been created automatically by the script Abilint (TD). 1405 !Do not modify the following lines by hand. 1406 #undef ABI_FUNC 1407 #define ABI_FUNC 'matrixelmt_g' 1408 !End of the abilint section 1409 1410 implicit none 1411 1412 !Arguments ------------------------------------ 1413 !scalars 1414 integer,intent(in) :: istwf_k,needimag,npw,nspinor,me_g0,comm_fft 1415 real(dp),intent(out) :: ai,ar 1416 !arrays 1417 real(dp),intent(in) :: diag(npw),vect1(2,npw*nspinor),vect2(2,npw*nspinor) 1418 1419 !Local variables------------------------------- 1420 !scalars 1421 integer :: i1,ierr,ipw 1422 character(len=500) :: message 1423 !arrays 1424 real(dp) :: buffer2(2) 1425 !real(dp),allocatable :: re_prod(:), im_prod(:) 1426 1427 ! ************************************************************************* 1428 1429 if (nspinor==2 .and. istwf_k/=1) then 1430 write(message,'(a,a,a,i6,a,i6)')& 1431 & 'When istwf_k/=1, nspinor must be 1,',ch10,& 1432 & 'however, nspinor=',nspinor,', and istwf_k=',istwf_k 1433 MSG_BUG(message) 1434 end if 1435 1436 #if 0 1437 !TODO 1438 ABI_MALLOC(re_prod,(npw*nspinor)) 1439 do ipw=1,npw*nspinor 1440 re_prod(ipw) = vect1(1,ipw)*vect2(1,ipw) + vect1(2,ipw)*vect2(2,ipw) 1441 end do 1442 1443 if (needimag == 1) then 1444 ABI_MALLOC(im_prod,(npw*nspinor)) 1445 do ipw=1,npw*nspinor 1446 im_prod(ipw) = vect1(1,ipw)*vect2(2,ipw) - vect1(2,ipw)*vect2(1,ipw) 1447 end do 1448 end if 1449 #endif 1450 1451 ar=zero 1452 if(needimag==1)ai=zero 1453 1454 !Normal storage mode 1455 if(istwf_k==1)then 1456 1457 ! Need only real part 1458 if(needimag==0)then 1459 1460 do ipw=1,npw 1461 ar=ar+diag(ipw)*(vect1(1,ipw)*vect2(1,ipw)+vect1(2,ipw)*vect2(2,ipw)) 1462 end do 1463 1464 if(nspinor==2)then 1465 do ipw=1+npw,2*npw 1466 ar=ar+diag(ipw-npw)*(vect1(1,ipw)*vect2(1,ipw)+vect1(2,ipw)*vect2(2,ipw)) 1467 end do 1468 end if 1469 1470 else ! Need also the imaginary part 1471 1472 do ipw=1,npw 1473 ar=ar+diag(ipw)*(vect1(1,ipw)*vect2(1,ipw)+vect1(2,ipw)*vect2(2,ipw)) 1474 ai=ai+diag(ipw)*(vect1(1,ipw)*vect2(2,ipw)-vect1(2,ipw)*vect2(1,ipw)) 1475 end do 1476 1477 if(nspinor==2)then 1478 do ipw=1+npw,2*npw 1479 ar=ar+diag(ipw-npw)*(vect1(1,ipw)*vect2(1,ipw)+vect1(2,ipw)*vect2(2,ipw)) 1480 ai=ai+diag(ipw-npw)*(vect1(1,ipw)*vect2(2,ipw)-vect1(2,ipw)*vect2(1,ipw)) 1481 end do 1482 end if 1483 1484 end if ! needimag 1485 1486 else if(istwf_k>=2)then 1487 1488 ! XG030513 : MPIWF need to know which proc has G=0 1489 1490 i1=1 1491 if(istwf_k==2 .and. me_g0==1)then 1492 ar=half*diag(1)*vect1(1,1)*vect2(1,1) ; i1=2 1493 end if 1494 1495 ! Need only real part 1496 if(needimag==0)then 1497 1498 do ipw=i1,npw 1499 ar=ar+diag(ipw)*(vect1(1,ipw)*vect2(1,ipw)+vect1(2,ipw)*vect2(2,ipw)) 1500 end do 1501 ar=two*ar 1502 1503 else ! Need also the imaginary part 1504 1505 do ipw=i1,npw 1506 ar=ar+diag(ipw)*(vect1(1,ipw)*vect2(1,ipw)+vect1(2,ipw)*vect2(2,ipw)) 1507 ai=ai+diag(ipw)*(vect1(1,ipw)*vect2(2,ipw)-vect1(2,ipw)*vect2(1,ipw)) 1508 end do 1509 ar=two*ar ; ai=two*ai 1510 1511 end if 1512 1513 end if ! istwf_k 1514 1515 #if 0 1516 ABI_FREE(re_prod) 1517 if (needimag == 1) then 1518 ABI_FREE(im_prod) 1519 end if 1520 #endif 1521 1522 !MPIWF need to make reduction on ar and ai . 1523 if (xmpi_comm_size(comm_fft)>1) then 1524 buffer2(1)=ai 1525 buffer2(2)=ar 1526 call xmpi_sum(buffer2,comm_fft ,ierr) 1527 ai=buffer2(1) 1528 ar=buffer2(2) 1529 end if 1530 1531 end subroutine matrixelmt_g
m_cgtools/mean_fftr [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
mean_fftr
FUNCTION
Compute the mean of an arraysp(nfft,nspden), over the FFT grid, for each component nspden,
and return it in meansp(nspden).
Take into account the spread of the array due to parallelism: the actual number of fft
points is nfftot, but the number of points on this proc is nfft only.
So : for ispden from 1 to nspden
meansp(ispden) = sum(ifft=1,nfftot) arraysp(ifft,ispden) / nfftot
INPUTS
arraysp(nfft,nspden)=the array whose average has to be computed nfft=number of FFT points stored by one proc nfftot=total number of FFT points nspden=number of spin-density components
OUTPUT
meansp(nspden)=mean value for each nspden component
PARENTS
fresid,newvtr,pawmknhat,prcref,prcref_PMA,psolver_rhohxc,rhohxcpositron rhotov,rhotoxc
CHILDREN
SOURCE
2028 subroutine mean_fftr(arraysp,meansp,nfft,nfftot,nspden,mpi_comm_sphgrid) 2029 2030 2031 !This section has been created automatically by the script Abilint (TD). 2032 !Do not modify the following lines by hand. 2033 #undef ABI_FUNC 2034 #define ABI_FUNC 'mean_fftr' 2035 !End of the abilint section 2036 2037 implicit none 2038 2039 !Arguments ------------------------------------ 2040 !scalars 2041 integer,intent(in) :: nfft,nfftot,nspden 2042 integer,intent(in),optional:: mpi_comm_sphgrid 2043 !arrays 2044 real(dp),intent(in) :: arraysp(nfft,nspden) 2045 real(dp),intent(out) :: meansp(nspden) 2046 2047 !Local variables------------------------------- 2048 !scalars 2049 integer :: ierr,ifft,ispden,nproc_sphgrid 2050 real(dp) :: invnfftot,tmean 2051 2052 ! ************************************************************************* 2053 2054 invnfftot=one/(dble(nfftot)) 2055 2056 do ispden=1,nspden 2057 tmean=zero 2058 !$OMP PARALLEL DO REDUCTION(+:tmean) 2059 do ifft=1,nfft 2060 tmean=tmean+arraysp(ifft,ispden) 2061 end do 2062 meansp(ispden)=tmean*invnfftot 2063 end do 2064 2065 !XG030514 : MPIWF The values of meansp(ispden) should 2066 !now be summed accross processors in the same WF group, and spread on all procs. 2067 if(present(mpi_comm_sphgrid)) then 2068 nproc_sphgrid=xmpi_comm_size(mpi_comm_sphgrid) 2069 if(nproc_sphgrid>1) then 2070 call xmpi_sum(meansp,nspden,mpi_comm_sphgrid,ierr) 2071 end if 2072 end if 2073 2074 end subroutine mean_fftr
m_cgtools/overlap_g [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
overlap_g
FUNCTION
Compute the scalar product between WF at two different k-points
< u_{n,k1} | u_{n,k2}>
INPUTS
mpw = maximum dimensioned size of npw npw_k1 = number of plane waves at k1 npw_k2 = number of plane waves at k2 nspinor = 1 for scalar, 2 for spinor wavefunctions pwind_k = array required to compute the scalar product (see initberry.f) vect1 = wavefunction at k1: | u_{n,k1} > vect2 = wavefunction at k1: | u_{n,k2} >
OUTPUT
doti = imaginary part of the scalarproduct dotr = real part of the scalarproduct
NOTES
In case a G-vector of the basis sphere of plane waves at k1 does not belong to the basis sphere of plane waves at k2, pwind = 0. Therefore, the dimensions of vect1 & vect2 are (1:2,0:mpw) and the element (1:2,0) MUST be set to zero. The current implementation if not compatible with TR-symmetry (i.e. istwfk/=1) !
PARENTS
dfptff_bec,dfptff_die,dfptff_ebp,dfptff_edie,dfptff_gbefd dfptff_gradberry,qmatrix,smatrix
CHILDREN
SOURCE
4746 subroutine overlap_g(doti,dotr,mpw,npw_k1,npw_k2,nspinor,pwind_k,vect1,vect2) 4747 4748 4749 !This section has been created automatically by the script Abilint (TD). 4750 !Do not modify the following lines by hand. 4751 #undef ABI_FUNC 4752 #define ABI_FUNC 'overlap_g' 4753 !End of the abilint section 4754 4755 implicit none 4756 4757 !Arguments ------------------------------------ 4758 !scalars 4759 integer,intent(in) :: mpw,npw_k1,npw_k2,nspinor 4760 real(dp),intent(out) :: doti,dotr 4761 !arrays 4762 integer,intent(in) :: pwind_k(mpw) 4763 real(dp),intent(in) :: vect1(1:2,0:mpw*nspinor),vect2(1:2,0:mpw*nspinor) 4764 4765 !Local variables------------------------------- 4766 !scalars 4767 integer :: ipw,ispinor,jpw,spnshft1,spnshft2 4768 character(len=500) :: message 4769 4770 ! ************************************************************************* 4771 4772 !Check if vect1(:,0) = 0 and vect2(:,0) = 0 4773 if ((abs(vect1(1,0)) > tol12).or.(abs(vect1(2,0)) > tol12).or. & 4774 & (abs(vect2(1,0)) > tol12).or.(abs(vect2(2,0)) > tol12)) then 4775 message = ' vect1(:,0) and/or vect2(:,0) are not equal to zero' 4776 MSG_BUG(message) 4777 end if 4778 4779 !Compute the scalar product 4780 dotr = zero; doti = zero 4781 do ispinor = 1, nspinor 4782 spnshft1 = (ispinor-1)*npw_k1 4783 spnshft2 = (ispinor-1)*npw_k2 4784 !$OMP PARALLEL DO PRIVATE(jpw) REDUCTION(+:doti,dotr) 4785 do ipw = 1, npw_k1 4786 jpw = pwind_k(ipw) 4787 dotr = dotr + vect1(1,spnshft1+ipw)*vect2(1,spnshft2+jpw) + vect1(2,spnshft1+ipw)*vect2(2,spnshft2+jpw) 4788 doti = doti + vect1(1,spnshft1+ipw)*vect2(2,spnshft2+jpw) - vect1(2,spnshft1+ipw)*vect2(1,spnshft2+jpw) 4789 end do 4790 end do 4791 4792 end subroutine overlap_g
m_cgtools/projbd [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
projbd
FUNCTION
Project out vector "direc" onto the bands contained in "cg".
if useoverlap==0
New direc=direc-$sum_{j/=i} { <cg_{j}|direc>.|cg_{j}> }$
if useoverlap==1 (use of overlap matrix S)
New direc=direc-$sum_{j/=i} { <cg_{j}|S|direc>.|cg_{j}> }$
(index i can be set to -1 to sum over all bands)
INPUTS
cg(2,mcg)=wavefunction coefficients for ALL bands
iband0=which particular band we are interested in ("i" in the above formula)
Can be set to -1 to sum over all bands...
icg=shift to be given to the location of the data in cg
iscg=shift to be given to the location of the data in cg
istwf_k=option parameter that describes the storage of wfs
mcg=maximum size of second dimension of cg
mscg=maximum size of second dimension of scg
nband=number of bands
npw=number of planewaves
nspinor=number of spinorial components (on current proc)
scg(2,mscg*useoverlap)=<G|S|band> for ALL bands,
where S is an overlap matrix
scprod_io=0 if scprod array has to be computed; 1 if it is input (already in memory)
tim_projbd=timing code of the calling subroutine(can be set to 0 if not attributed)
useoverlap=describe the overlap of wavefunctions:
0: no overlap (S=Identity_matrix)
1: wavefunctions are overlapping
me_g0=1 if this processors treats G=0, 0 otherwise.
comm=MPI communicator (used if G vectors are distributed.
SIDE EFFECTS
direc(2,npw)= input: vector to be orthogonalised with respect to cg (and S)
output: vector that has been orthogonalized wrt cg (and S)
scprod(2,nband)=scalar_product
if useoverlap==0: scalar_product_i=$<cg_{j}|direc_{i}>$
if useoverlap==1: scalar_product_i=$<cg_{j}|S|direc_{i}>$
if scprod_io=0, scprod is output
if scprod_io=1, scprod is input
NOTES
1) MPIWF Might have to be recoded for efficient paralellism
2) The new version employs BLAS2 routine so that the OMP parallelism is delegated to BLAS library.
3) Note for PAW: ref.= PRB 73, 235101 (2006), equations (71) and (72):
in normal use, projbd applies P_c projector
if cg and scg are inverted, projbd applies P_c+ projector
4) cg_zgemv wraps ZGEMM whose implementation is more efficient, especially in the threaded case.
PARENTS
cgwf,dfpt_cgwf,dfpt_nstpaw,getdc1,lapackprof
CHILDREN
SOURCE
3528 subroutine projbd(cg,direc,iband0,icg,iscg,istwf_k,mcg,mscg,nband,& 3529 & npw,nspinor,scg,scprod,scprod_io,tim_projbd,useoverlap,me_g0,comm) 3530 3531 3532 !This section has been created automatically by the script Abilint (TD). 3533 !Do not modify the following lines by hand. 3534 #undef ABI_FUNC 3535 #define ABI_FUNC 'projbd' 3536 use interfaces_18_timing 3537 !End of the abilint section 3538 3539 implicit none 3540 3541 !Arguments ------------------------------------ 3542 !scalars 3543 integer,intent(in) :: iband0,icg,iscg,istwf_k,mcg,mscg,nband,npw,nspinor 3544 integer,intent(in) :: scprod_io,tim_projbd,useoverlap,me_g0,comm 3545 !arrays 3546 real(dp),intent(in) :: cg(2,mcg),scg(2,mscg*useoverlap) 3547 real(dp),intent(inout) :: direc(2,npw*nspinor) 3548 real(dp),intent(inout) :: scprod(2,nband) 3549 3550 !Local variables------------------------------- 3551 !scalars 3552 integer :: nbandm,npw_sp,ierr 3553 !arrays 3554 real(dp) :: tsec(2),bkp_scprod(2),bkp_dirg0(2) 3555 3556 ! ************************************************************************* 3557 3558 call timab(210+tim_projbd,1,tsec) 3559 3560 npw_sp=npw*nspinor 3561 3562 nbandm=nband 3563 3564 if (istwf_k==1) then 3565 3566 if (scprod_io==0) then 3567 if (useoverlap==1) then 3568 call cg_zgemv("C",npw_sp,nbandm,scg(1,iscg+1),direc,scprod) 3569 else 3570 call cg_zgemv("C",npw_sp,nbandm,cg(1,icg+1), direc,scprod) 3571 end if 3572 call xmpi_sum(scprod,comm,ierr) 3573 end if 3574 3575 if (iband0>0.and.iband0<=nbandm) then 3576 bkp_scprod = scprod(:,iband0) 3577 scprod(:,iband0) = zero 3578 end if 3579 3580 call cg_zgemv("N",npw_sp,nbandm,cg(1,icg+1),scprod,direc,alpha=-cg_cone,beta=cg_cone) 3581 3582 if (iband0>0.and.iband0<=nbandm) scprod(:,iband0) = bkp_scprod ! Restore previous value as scprod is output. 3583 3584 else if (istwf_k>=2) then 3585 ! 3586 ! u_{G0/2}(G) = u_{G0/2}(-G-G0)^*; k = G0/2 3587 ! hence: 3588 ! sum_G f*(G) g(G) = 2 REAL sum_G^{IZONE} w(G) f*(G)g(G) 3589 ! where 3590 ! w(G) = 1 except for k=0 and G=0 where w(G=0) = 1/2. 3591 ! 3592 if (scprod_io==0) then 3593 3594 if (useoverlap==1) then 3595 3596 if (istwf_k==2 .and. me_g0==1) then 3597 bkp_dirg0 = direc(:,1) 3598 direc(1,1) = half * direc(1,1) 3599 direc(2,1) = zero 3600 end if 3601 3602 call cg_zgemv("C",npw_sp,nbandm,scg(1,iscg+1),direc,scprod) 3603 scprod = two * scprod 3604 scprod(2,:) = zero 3605 3606 if(istwf_k==2 .and. me_g0==1) direc(:,1) = bkp_dirg0 3607 3608 else 3609 3610 if (istwf_k==2 .and. me_g0==1) then 3611 bkp_dirg0 = direc(:,1) 3612 direc(1,1) = half * direc(1,1) 3613 direc(2,1) = zero 3614 end if 3615 3616 call cg_zgemv("C",npw_sp,nbandm,cg(1,icg+1),direc,scprod) 3617 scprod = two * scprod 3618 scprod(2,:) = zero 3619 3620 if (istwf_k==2 .and. me_g0==1) direc(:,1) = bkp_dirg0 3621 end if ! useoverlap 3622 3623 call xmpi_sum(scprod,comm,ierr) 3624 end if 3625 3626 if (iband0>0.and.iband0<=nbandm) then 3627 bkp_scprod = scprod(:,iband0) 3628 scprod(:,iband0) = zero 3629 end if 3630 3631 call cg_zgemv("N",npw_sp,nbandm,cg(1,icg+1),scprod,direc,alpha=-cg_cone,beta=cg_cone) 3632 3633 if (iband0>0.and.iband0<=nbandm) scprod(:,iband0) = bkp_scprod ! Restore previous value as scprod is output. 3634 3635 end if ! Test on istwf_k 3636 3637 call timab(210+tim_projbd,2,tsec) 3638 3639 end subroutine projbd
m_cgtools/set_istwfk [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
set_istwfk
FUNCTION
Returns the value of istwfk associated to the input k-point.
INPUTS
kpoint(3)=The k-point in reduced coordinates.
OUTPUT
istwfk= Integer flag internally used in the code to define the storage mode of the wavefunctions. It also define the algorithm used to apply an operator in reciprocal space as well as the FFT algorithm used to go from G- to r-space and vice versa. 1 => time-reversal cannot be used 2 => use time-reversal at the Gamma point. 3 => use time-reversal symmetry for k=(1/2, 0 , 0 ) 4 => use time-reversal symmetry for k=( 0 , 0 ,1/2) 5 => use time-reversal symmetry for k=(1/2, 0 ,1/2) 6 => use time-reversal symmetry for k=( 0 ,1/2, 0 ) 7 => use time-reversal symmetry for k=(1/2,1/2, 0 ) 8 => use time-reversal symmetry for k=( 0 ,1/2,1/2) 9 => use time-reversal symmetry for k=(1/2,1/2,1/2) Useful relations: u_k(G) = u_{k+G0}(G-G0); u_{-k}(G) = u_k(G)^* and therefore: u_{G0/2}(G) = u_{G0/2}(-G-G0)^*.
PARENTS
CHILDREN
SOURCE
1145 integer pure function set_istwfk(kpoint) result(istwfk) 1146 1147 1148 !This section has been created automatically by the script Abilint (TD). 1149 !Do not modify the following lines by hand. 1150 #undef ABI_FUNC 1151 #define ABI_FUNC 'set_istwfk' 1152 !End of the abilint section 1153 1154 implicit none 1155 1156 !Arguments ------------------------------------ 1157 real(dp),intent(in) :: kpoint(3) 1158 1159 !Local variables------------------------------- 1160 !scalars 1161 integer :: bit0,ii 1162 !arrays 1163 integer :: bit(3) 1164 1165 ! ************************************************************************* 1166 1167 bit0=1 1168 1169 do ii=1,3 1170 if (DABS(kpoint(ii))<tol10) then 1171 bit(ii)=0 1172 else if (DABS(kpoint(ii)-half)<tol10 ) then 1173 bit(ii)=1 1174 else 1175 bit0=0 1176 end if 1177 end do 1178 1179 if (bit0==0) then 1180 istwfk=1 1181 else 1182 istwfk=2+bit(1)+4*bit(2)+2*bit(3) ! Note the inversion between bit(2) and bit(3) 1183 end if 1184 1185 end function set_istwfk
m_cgtools/sqnorm_g [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
sqnorm_g
FUNCTION
Compute the square of the norm of one complex vector vecti, in reciprocal space Take into account the storage mode of the vector (istwf_k)
INPUTS
istwf_k=option parameter that describes the storage of wfs npwsp= (effective) number of planewaves at this k point. vect(2,npwsp)=the vector in reciprocal space (npw*nspinor, usually) me_g0=1 if this processors treats G=0, 0 otherwise. comm=MPI communicator for MPI sum.
OUTPUT
dotr= <vect|vect>
PARENTS
cgwf,dfpt_cgwf,dfpt_vtowfk,m_epjdos,mkresi,rf2_init
CHILDREN
SOURCE
1213 subroutine sqnorm_g(dotr,istwf_k,npwsp,vect,me_g0,comm) 1214 1215 1216 !This section has been created automatically by the script Abilint (TD). 1217 !Do not modify the following lines by hand. 1218 #undef ABI_FUNC 1219 #define ABI_FUNC 'sqnorm_g' 1220 !End of the abilint section 1221 1222 implicit none 1223 1224 !Arguments ------------------------------------ 1225 !scalars 1226 integer,intent(in) :: istwf_k,npwsp,me_g0,comm 1227 real(dp),intent(out) :: dotr 1228 !arrays 1229 real(dp),intent(in) :: vect(2,npwsp) 1230 1231 !Local variables------------------------------- 1232 !scalars 1233 integer :: ierr 1234 1235 ! ************************************************************************* 1236 1237 if (istwf_k==1) then ! General k-point 1238 !dotr = cg_real_zdotc(npwsp,vect,vect) 1239 dotr = cg_dznrm2(npwsp, vect) 1240 dotr = dotr * dotr 1241 1242 else 1243 if (istwf_k==2 .and. me_g0==1) then 1244 ! Gamma k-point and I have G=0 1245 dotr=half*vect(1,1)**2 1246 dotr = dotr + cg_real_zdotc(npwsp-1,vect(1,2),vect(1,2)) 1247 else 1248 ! Other TR k-points 1249 dotr = cg_real_zdotc(npwsp,vect,vect) 1250 end if 1251 dotr=two*dotr 1252 end if 1253 1254 if (xmpi_comm_size(comm)>1) then 1255 call xmpi_sum(dotr,comm,ierr) 1256 end if 1257 1258 end subroutine sqnorm_g
m_cgtools/sqnorm_v [ Functions ]
[ Top ] [ m_cgtools ] [ Functions ]
NAME
sqnorm_v
FUNCTION
Compute square of the norm of a potential (integral over FFT grid), to obtain a square residual-like quantity (so the sum of product of values is NOT divided by the number of FFT points, and NOT multiplied by the primitive cell volume). Take into account the spin components of the potentials (nspden), and sum over them.
INPUTS
cplex=if 1, real space function on FFT grid is REAL, if 2, COMPLEX nfft= (effective) number of FFT grid points (for this processor) nspden=number of spin-density components opt_storage: 0, if potential is stored as V^up-up, V^dn-dn, Re[V^up-dn], Im[V^up-dn] 1, if potential is stored as V, B_x, B_y, Bz (B=magn. field) pot(cplex*nfft,nspden)=real space potential on FFT grid
OUTPUT
norm2= value of the square of the norm
PARENTS
dfpt_rhotov,dfpt_vtorho,rhotov,vtorho
CHILDREN
SOURCE
1932 subroutine sqnorm_v(cplex,nfft,norm2,nspden,opt_storage,pot,mpi_comm_sphgrid) 1933 1934 1935 !This section has been created automatically by the script Abilint (TD). 1936 !Do not modify the following lines by hand. 1937 #undef ABI_FUNC 1938 #define ABI_FUNC 'sqnorm_v' 1939 !End of the abilint section 1940 1941 implicit none 1942 1943 !Arguments ------------------------------------ 1944 !scalars 1945 integer,intent(in) :: cplex,nfft,nspden,opt_storage 1946 integer,intent(in),optional :: mpi_comm_sphgrid 1947 real(dp),intent(out) :: norm2 1948 !arrays 1949 real(dp),intent(in) :: pot(cplex*nfft,nspden) 1950 1951 !Local variables------------------------------- 1952 !scalars 1953 integer :: ierr,ifft,ispden,nproc_sphgrid 1954 real(dp) :: ar 1955 1956 ! ************************************************************************* 1957 1958 !Real or complex inputs are coded 1959 1960 norm2=zero 1961 do ispden=1,min(nspden,2) 1962 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ispden,nfft,pot) REDUCTION(+:norm2) 1963 do ifft=1,cplex*nfft 1964 norm2=norm2 + pot(ifft,ispden)**2 1965 end do 1966 end do 1967 if (nspden==4) then 1968 ar=zero 1969 do ispden=3,4 1970 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ispden,nfft,pot) REDUCTION(+:ar) 1971 do ifft=1,cplex*nfft 1972 ar=ar + pot(ifft,ispden)**2 1973 end do 1974 end do 1975 if (opt_storage==0) then 1976 if (cplex==1) then 1977 norm2=norm2+two*ar 1978 else 1979 norm2=norm2+ar 1980 end if 1981 else 1982 norm2=half*(norm2+ar) 1983 end if 1984 end if 1985 1986 !MPIWF reduction (addition) on norm2 is needed here 1987 if(present(mpi_comm_sphgrid)) then 1988 nproc_sphgrid=xmpi_comm_size(mpi_comm_sphgrid) 1989 if(nproc_sphgrid>1)then 1990 call xmpi_sum(norm2,mpi_comm_sphgrid,ierr) 1991 end if 1992 end if 1993 1994 end subroutine sqnorm_v