TABLE OF CONTENTS
- ABINIT/dfptff_bec
- ABINIT/dfptff_die
- ABINIT/dfptff_ebp
- ABINIT/dfptff_edie
- ABINIT/dfptff_gbefd
- ABINIT/dfptff_gradberry
- ABINIT/dfptff_initberry
- ABINIT/m_dfpt_fef
- ABINIT/qmatrix
ABINIT/dfptff_bec [ Functions ]
NAME
dfptff_bec
FUNCTION
calculate Born effective charge tensor in Eq.(33) in PRB 75, 115116(2007) [[cite:Wang2007]].
INPUTS
cg(2,mpw*nspinor*mband_mem*mkmem_rbz*nsppol) = planewave coefficients of wavefunctions cg1(2,mpw1*nspinor*mband_mem*mk1mem*nsppol) = pw coefficients of RF wavefunctions at k,q. dtefield = variables related to response Berry-phase calculation idirpert = the current coloumn of the dielectric permittivity tensor mband = maximum number of bands mband_mem = maximum number of bands on this cpu mkmem_rbz = maximum number of k-points in core memory mpw = maximum number of plane waves mpw1 = maximum number of plane waves for response wavefunctions nkpt = number of k points npwarr(nkpt) = number of planewaves in basis and boundary at this k point npwar1(nkpt) = number of planewaves in basis and boundary for response wfs nspinor = 1 for scalar wfs, 2 for spinor wfs nsppol = 1 for unpolarized, 2 for spin-polarized qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3) = inverse of the overlap matrix pwindall(max(mpw,mpw1)*mkmem_rbz,8,3) = array used to compute the overlap matrices pwindall(:,1,:) <- <u^(0)_i|u^(0)_i+1> pwindall(:,2,:) <- <u^(0)_i|u^(0)_i-1> pwindall(:,3,:) <- <u^(1)_i|u^(1)_i+1> pwindall(:,4,:) <- <u^(1)_i|u^(1)_i-1> pwindall(:,5,:) <- <u^(1)_i|u^(0)_i+n+1> pwindall(:,6,:) <- <u^(1)_i|u^(0)_i+n-1> pwindall(:,7,:) <- <u^(0)_i|u^(1)_i-n+1> pwindall(:,8,:) <- <u^(0)_i|u^(1)_i-n-1> rprimd(3,3)=dimensional primitive translations in real space (bohr)
OUTPUT
d2lo(1,1:3,natom+5,1:3,1:natom) = Born effective charge tensor
SOURCE
2624 subroutine dfptff_bec(cg,cg1,dtefield,natom,d2lo,idirpert,ipert,mband,mband_mem,mkmem_rbz,& 2625 & mpi_enreg,mpw,mpw1,mpert,nkpt,npwarr,npwar1,nsppol,nspinor,pwindall,qmat,rprimd) 2626 2627 !Arguments ---------------------------------------- 2628 !scalars 2629 integer,intent(in) :: idirpert,ipert,mband,mkmem_rbz,mpert,mpw,mpw1,natom,nkpt 2630 integer,intent(in) :: mband_mem 2631 integer,intent(in) :: nspinor,nsppol 2632 type(efield_type),intent(in) :: dtefield 2633 type(MPI_type),intent(in) :: mpi_enreg 2634 !arrays 2635 integer,intent(in) :: npwar1(nkpt),npwarr(nkpt) 2636 integer,intent(in) :: pwindall(max(mpw,mpw1)*mkmem_rbz,8,3) 2637 real(dp),intent(in) :: cg(2,mpw*mband_mem*mkmem_rbz*nspinor*nsppol) 2638 real(dp),intent(in) :: cg1(2,mpw1*mband_mem*mkmem_rbz*nspinor*nsppol) 2639 real(dp),intent(in) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3) 2640 real(dp),intent(in) :: rprimd(3,3) 2641 real(dp),intent(inout) :: d2lo(2,3,mpert,3,mpert) !vz_i 2642 2643 !Local variables ---------------------------------- 2644 !scalars 2645 integer :: ialpha,iband,icg,icg1,idir,ikpt,ikpt1,jband,mpw_tmp,npw_k1 2646 integer :: npw_k2,pwmax,pwmin 2647 integer :: me_band,ierr, iband_me, jband_me 2648 real(dp) :: doti,dotr,e0,fac 2649 !arrays 2650 integer,allocatable :: pwind_tmp(:) 2651 integer :: band_procs(mband) 2652 real(dp) :: edir(3) 2653 real(dp),allocatable :: s1mat(:,:,:),vect1(:,:),vect2(:,:) 2654 2655 ! ************************************************************************* 2656 2657 !calculate s1 matrices ----------------------------- 2658 mpw_tmp=max(mpw,mpw1) 2659 ABI_MALLOC(s1mat,(2,dtefield%mband_occ,dtefield%mband_occ)) 2660 ABI_MALLOC(vect1,(2,0:mpw_tmp)) 2661 ABI_MALLOC(vect2,(2,0:mpw_tmp)) 2662 ABI_MALLOC(pwind_tmp,(mpw_tmp)) 2663 vect1(:,0) = zero ; vect2(:,0) = zero 2664 2665 me_band = mpi_enreg%me_band 2666 2667 edir(:)=zero 2668 2669 !TODO MJV: where is the loop over sppol? proc_distrb_band_procs chooses isppol 1 below 2670 do ikpt=1,nkpt 2671 call proc_distrb_band(band_procs,mpi_enreg%proc_distrb,ikpt,1,mband,& 2672 & me_band,mpi_enreg%me_kpt,mpi_enreg%comm_band) 2673 2674 do idir=1,3 2675 2676 ! compute <u^(0)_{k_j}|u^(1)_{k_j+1,q}> matrix--- q=0 ---------------------------------------- 2677 2678 ! prepare to calculate overlap matrix 2679 ikpt1 = dtefield%ikpt_dk(ikpt,1,idir) 2680 icg = dtefield%cgindex(ikpt,1) 2681 icg1 = dtefield%cgindex(ikpt1,1+nsppol) 2682 npw_k1 = npwarr(ikpt) 2683 npw_k2 = npwar1(ikpt1) 2684 pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,7,idir) 2685 vect1(:,0) = zero ; vect2(:,0) = zero 2686 jband_me = 0 2687 s1mat = zero 2688 do jband = 1, dtefield%mband_occ 2689 if (band_procs(jband) == me_band) then 2690 jband_me = jband_me + 1 2691 vect2(:,1:npw_k2) = & 2692 & cg1(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor) 2693 end if 2694 call xmpi_bcast(vect2,band_procs(jband), mpi_enreg%comm_band,ierr) 2695 2696 ! now everyone has vect2 for present jband 2697 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 2698 2699 iband_me = 0 2700 do iband = 1, dtefield%mband_occ 2701 if (band_procs(iband) /= me_band) cycle 2702 iband_me = iband_me + 1 2703 2704 pwmin = (iband_me-1)*npw_k1*nspinor 2705 pwmax = pwmin + npw_k1*nspinor 2706 vect1(:,1:npw_k1) = & 2707 & cg(:,icg + 1 + pwmin:icg + pwmax) 2708 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 2709 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 2710 & vect1,vect2) 2711 s1mat(1,iband,jband) = dotr 2712 s1mat(2,iband,jband) = doti 2713 end do ! iband 2714 end do !jband 2715 ! for the moment s1mat is only filled for iband on current cpu 2716 2717 2718 ! compute <u^(1)_{k_j,q}|u^(0)_{k_j+1}> matrix-- q=0 ------------------------------------- 2719 2720 ! prepare to calculate overlap matrix 2721 ikpt1 = dtefield%ikpt_dk(ikpt,1,idir) 2722 icg = dtefield%cgindex(ikpt,1+nsppol) 2723 icg1 = dtefield%cgindex(ikpt1,1) 2724 npw_k1 = npwar1(ikpt) 2725 npw_k2 = npwarr(ikpt1) 2726 pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,5,idir) 2727 vect1(:,0) = zero ; vect2(:,0) = zero 2728 jband_me = 0 2729 do jband = 1, dtefield%mband_occ 2730 if (band_procs(jband) == me_band) then 2731 jband_me = jband_me + 1 2732 vect2(:,1:npw_k2) = & 2733 & cg(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor) 2734 end if 2735 call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr) 2736 2737 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 2738 iband_me = 0 2739 do iband = 1, dtefield%mband_occ 2740 if (band_procs(iband) /= me_band) cycle 2741 iband_me = iband_me + 1 2742 2743 pwmin = (iband_me-1)*npw_k1*nspinor 2744 pwmax = pwmin + npw_k1*nspinor 2745 vect1(:,1:npw_k1) = & 2746 & cg1(:,icg + 1 + pwmin:icg + pwmax) 2747 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 2748 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 2749 & vect1,vect2) 2750 2751 s1mat(1,iband,jband) = s1mat(1,iband,jband) + dotr 2752 s1mat(2,iband,jband) = s1mat(2,iband,jband) + doti 2753 end do ! iband 2754 end do !jband 2755 2756 ! recompose full s1mat on each proc 2757 call xmpi_sum(s1mat,mpi_enreg%comm_band,ierr) 2758 2759 ! sum over the whole------------------------------------------------------------ 2760 2761 e0=zero 2762 2763 do iband=1,dtefield%mband_occ 2764 do jband=1,dtefield%mband_occ 2765 e0 = e0 + (s1mat(1,iband,jband)*qmat(2,jband,iband,ikpt,1,idir)& 2766 & + s1mat(2,iband,jband)*qmat(1,jband,iband,ikpt,1,idir)) 2767 end do 2768 end do 2769 2770 do ialpha=1,3 2771 fac = rprimd(ialpha,idir)/& 2772 & (dble(dtefield%nstr(idir))*pi) 2773 2774 edir(ialpha)=edir(ialpha)+ e0*fac 2775 end do 2776 2777 end do ! idir 2778 end do ! ikpt 2779 2780 d2lo(1,1:3,natom+2,idirpert,ipert)=edir(:) 2781 2782 ABI_FREE(s1mat) 2783 ABI_FREE(vect1) 2784 ABI_FREE(vect2) 2785 ABI_FREE(pwind_tmp) 2786 2787 end subroutine dfptff_bec
ABINIT/dfptff_die [ Functions ]
NAME
dfptff_die
FUNCTION
calculate electric susceptibility tensor in Eq.(28) in PRB 75, 115116(2007) [[cite:Wang2007]].
INPUTS
cg(2,mpw*nspinor*mband*mkmem_rbz*nsppol) = planewave coefficients of wavefunctions cg1(2,mpw1*nspinor*mband*mk1mem*nsppol) = pw coefficients of RF wavefunctions at k,q. dtefield = variables related to response Berry-phase calculation idirpert = the current coloumn of the dielectric permittivity tensor mband = maximum number of bands mkmem_rbz = maximum number of k-points in core memory mpw = maximum number of plane waves mpw1 = maximum number of plane waves for response wavefunctions nkpt = number of k points npwarr(nkpt) = number of planewaves in basis and boundary at this k point npwar1(nkpt) = number of planewaves in basis and boundary for response wfs nspinor = 1 for scalar wfs, 2 for spinor wfs nsppol = 1 for unpolarized, 2 for spin-polarized qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3) = inverse of the overlap matrix pwindall(max(mpw,mpw1)*mkmem_rbz,8,3) = array used to compute the overlap matrices pwindall(:,1,:) <- <u^(0)_i|u^(0)_i+1> pwindall(:,2,:) <- <u^(0)_i|u^(0)_i-1> pwindall(:,3,:) <- <u^(1)_i|u^(1)_i+1> pwindall(:,4,:) <- <u^(1)_i|u^(1)_i-1> pwindall(:,5,:) <- <u^(1)_i|u^(0)_i+n+1> pwindall(:,6,:) <- <u^(1)_i|u^(0)_i+n-1> pwindall(:,7,:) <- <u^(0)_i|u^(1)_i-n+1> pwindall(:,8,:) <- <u^(0)_i|u^(1)_i-n-1> rprimd(3,3)=dimensional primitive translations in real space (bohr)
OUTPUT
diet = electric susceptibility tensor
SOURCE
2426 subroutine dfptff_die(cg,cg1,dtefield,d2lo,idirpert,ipert,mband,mband_mem,mkmem_rbz,& 2427 & mpi_enreg,mpw,mpw1,mpert,nkpt,npwarr,npwar1,nsppol,nspinor,pwindall,qmat,rprimd) 2428 2429 !Arguments ---------------------------------------- 2430 !scalars 2431 integer,intent(in) :: idirpert,ipert,mband,mkmem_rbz,mpert,mpw,mpw1,nkpt,nspinor 2432 integer,intent(in) :: mband_mem 2433 integer,intent(in) :: nsppol 2434 type(efield_type),intent(in) :: dtefield 2435 type(MPI_type),intent(in) :: mpi_enreg 2436 !arrays 2437 integer,intent(in) :: npwar1(nkpt),npwarr(nkpt) 2438 integer,intent(in) :: pwindall(max(mpw,mpw1)*mkmem_rbz,8,3) 2439 real(dp),intent(in) :: cg(2,mpw*mband_mem*mkmem_rbz*nspinor*nsppol) 2440 real(dp),intent(in) :: cg1(2,mpw1*mband_mem*mkmem_rbz*nspinor*nsppol) 2441 real(dp),intent(in) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3) 2442 real(dp),intent(in) :: rprimd(3,3) 2443 real(dp),intent(inout) :: d2lo(2,3,mpert,3,mpert) !vz_i 2444 2445 !Local variables ---------------------------------- 2446 !scalars 2447 integer :: ialpha,iband,icg,icg1,idir,ikpt,ikpt1,jband,mpw_tmp,npw_k1 2448 integer :: jband_me, iband_me, me_band, ierr 2449 integer :: npw_k2,pwmax,pwmin 2450 real(dp) :: doti,dotr,e0,fac 2451 !arrays 2452 integer :: band_procs(mband) 2453 integer,allocatable :: pwind_tmp(:) 2454 real(dp) :: edir(3) 2455 real(dp),allocatable :: s1mat(:,:,:),vect1(:,:),vect2(:,:) 2456 2457 ! ************************************************************************* 2458 2459 !calculate s1 matrices ----------------------------- 2460 mpw_tmp=max(mpw,mpw1) 2461 ABI_MALLOC(s1mat,(2,dtefield%mband_occ,dtefield%mband_occ)) 2462 ABI_MALLOC(vect1,(2,0:mpw_tmp)) 2463 ABI_MALLOC(vect2,(2,0:mpw_tmp)) 2464 ABI_MALLOC(pwind_tmp,(mpw_tmp)) 2465 vect1(:,0) = zero ; vect2(:,0) = zero 2466 2467 edir(:)=zero 2468 2469 me_band = mpi_enreg%me_band 2470 2471 !TODO no info on sppol here - I default to isppol 1 2472 do ikpt=1,nkpt 2473 call proc_distrb_band(band_procs,mpi_enreg%proc_distrb,ikpt,1,mband,& 2474 & me_band,mpi_enreg%me_kpt,mpi_enreg%comm_band) 2475 do idir=1,3 2476 ! compute <u^(0)_{k_j}|u^(1)_{k_j+1,q}> matrix--- q=0 ---------------------------------------- 2477 2478 ! prepare to calculate overlap matrix 2479 ikpt1 = dtefield%ikpt_dk(ikpt,1,idir) 2480 icg = dtefield%cgindex(ikpt,1) 2481 icg1 = dtefield%cgindex(ikpt1,1+nsppol) 2482 npw_k1 = npwarr(ikpt) 2483 npw_k2 = npwar1(ikpt1) 2484 pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,7,idir) 2485 2486 vect1(:,0) = zero ; vect2(:,0) = zero 2487 jband_me = 0 2488 do jband = 1, dtefield%mband_occ 2489 if (band_procs(jband) == me_band) then 2490 jband_me = jband_me + 1 2491 vect2(:,1:npw_k2) = & 2492 & cg1(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor) 2493 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 2494 end if 2495 call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr) 2496 2497 iband_me = 0 2498 do iband = 1, dtefield%mband_occ 2499 if (band_procs(iband) /= me_band) cycle 2500 iband_me = iband_me + 1 2501 pwmin = (iband_me-1)*npw_k1*nspinor 2502 pwmax = pwmin + npw_k1*nspinor 2503 vect1(:,1:npw_k1) = & 2504 & cg(:,icg + 1 + pwmin:icg + pwmax) 2505 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 2506 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 2507 & vect1,vect2) 2508 s1mat(1,iband,jband) = dotr 2509 s1mat(2,iband,jband) = doti 2510 end do ! iband 2511 end do !jband 2512 2513 ! compute <u^(1)_{k_j,q}|u^(0)_{k_j+1}> matrix-- q=0 ------------------------------------- 2514 2515 ! prepare to calculate overlap matrix 2516 ikpt1 = dtefield%ikpt_dk(ikpt,1,idir) 2517 icg = dtefield%cgindex(ikpt,1+nsppol) 2518 icg1 = dtefield%cgindex(ikpt1,1) 2519 npw_k1 = npwar1(ikpt) 2520 npw_k2 = npwarr(ikpt1) 2521 pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,5,idir) 2522 vect1(:,0) = zero ; vect2(:,0) = zero 2523 jband_me = 0 2524 do jband = 1, dtefield%mband_occ 2525 if (band_procs(jband) == me_band) then 2526 jband_me = jband_me + 1 2527 vect2(:,1:npw_k2) = & 2528 & cg(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor) 2529 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 2530 end if 2531 call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr) 2532 2533 iband_me = 0 2534 do iband = 1, dtefield%mband_occ 2535 if (band_procs(iband) /= me_band) cycle 2536 iband_me = iband_me + 1 2537 pwmin = (iband_me-1)*npw_k1*nspinor 2538 pwmax = pwmin + npw_k1*nspinor 2539 vect1(:,1:npw_k1) = & 2540 & cg1(:,icg + 1 + pwmin:icg + pwmax) 2541 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 2542 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 2543 & vect1,vect2) 2544 s1mat(1,iband,jband) = s1mat(1,iband,jband) + dotr 2545 s1mat(2,iband,jband) = s1mat(2,iband,jband) + doti 2546 end do ! iband 2547 end do !jband 2548 2549 ! recompose full s1mat on each proc, summing over iband 2550 call xmpi_sum(s1mat,mpi_enreg%comm_band,ierr) 2551 2552 ! sum over the whole------------------------------------------------------------ 2553 2554 e0=zero 2555 2556 do iband=1,dtefield%mband_occ 2557 do jband=1,dtefield%mband_occ 2558 e0 = e0 + (s1mat(1,iband,jband)*qmat(2,jband,iband,ikpt,1,idir)& 2559 & + s1mat(2,iband,jband)*qmat(1,jband,iband,ikpt,1,idir)) 2560 2561 end do 2562 end do 2563 2564 do ialpha=1,3 2565 fac = rprimd(ialpha,idir)/& 2566 & (dble(dtefield%nstr(idir))*pi) 2567 edir(ialpha)=edir(ialpha)+ e0*fac 2568 end do 2569 2570 end do !idir 2571 end do !ikpt 2572 2573 d2lo(1,1:3,ipert,idirpert,ipert)=edir(:) 2574 2575 ABI_FREE(s1mat) 2576 ABI_FREE(vect1) 2577 ABI_FREE(vect2) 2578 ABI_FREE(pwind_tmp) 2579 2580 end subroutine dfptff_die
ABINIT/dfptff_ebp [ Functions ]
NAME
dfptff_ebp
FUNCTION
calculation of the energy from the term \Omega E \cdot P
INPUTS
cg(2,mpw*nspinor*mband*mkmem_rbz*nsppol) = planewave coefficients of wavefunctions cg1(2,mpw1*nspinor*mband*mk1mem*nsppol) = pw coefficients of RF wavefunctions at k,q. dtefield = variables related to response Berry-phase calculation ikpt = the index of the current k point isppol = the index of the spin component mband = maximum number of bands mkmem_rbz = maximum number of k-points in core memory mpw = maximum number of plane waves mpw1 = maximum number of plane waves for response wavefunctions nkpt = number of k points npwarr(nkpt) = number of planewaves in basis and boundary at this k point npwar1(nkpt) = number of planewaves in basis and boundary for response wfs nspinor = 1 for scalar wfs, 2 for spinor wfs nsppol = 1 for unpolarized, 2 for spin-polarized qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3) = inverse of the overlap matrix pwindall(max(mpw,mpw1)*mkmem_rbz,8,3) = array used to compute the overlap matrices pwindall(:,1,:) <- <u^(0)_i|u^(0)_i+1> pwindall(:,2,:) <- <u^(0)_i|u^(0)_i-1> pwindall(:,3,:) <- <u^(1)_i|u^(1)_i+1> pwindall(:,4,:) <- <u^(1)_i|u^(1)_i-1> pwindall(:,5,:) <- <u^(1)_i|u^(0)_i+n+1> pwindall(:,6,:) <- <u^(1)_i|u^(0)_i+n-1> pwindall(:,7,:) <- <u^(0)_i|u^(1)_i-n+1> pwindall(:,8,:) <- <u^(0)_i|u^(1)_i-n-1>
OUTPUT
grad_berry(2,mpw1,dtefield%mband_occ) = the gradient of the Berry phase term
SOURCE
2104 subroutine dfptff_ebp(cg,cg1,dtefield,eberry,mband,mband_mem,mkmem_rbz,& 2105 & mpi_enreg,mpw,mpw1,nkpt,npwarr,npwar1,nsppol,nspinor,pwindall,qmat) 2106 2107 !Arguments ---------------------------------------- 2108 !scalars 2109 integer,intent(in) :: mband,mkmem_rbz,mpw,mpw1,nkpt,nspinor,nsppol 2110 integer,intent(in) :: mband_mem 2111 real(dp),intent(out) :: eberry 2112 type(efield_type),intent(in) :: dtefield 2113 type(MPI_type),intent(in) :: mpi_enreg 2114 !arrays 2115 integer,intent(in) :: npwar1(nkpt),npwarr(nkpt) 2116 integer,intent(in) :: pwindall(max(mpw,mpw1)*mkmem_rbz,8,3) 2117 real(dp),intent(in) :: cg(2,mpw*mband_mem*mkmem_rbz*nspinor*nsppol) 2118 real(dp),intent(in) :: cg1(2,mpw1*mband_mem*mkmem_rbz*nspinor*nsppol) 2119 real(dp),intent(in) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3) 2120 2121 !Local variables ---------------------------------- 2122 !scalars 2123 integer :: iband,icg,icg1,idir 2124 integer :: jband_me, iband_me, me_band, ierr 2125 integer :: ikpt,ikpt1,ikptn,ikptnm 2126 integer :: jband,kband,mpw_tmp,npw_k1,npw_k2,pwmax,pwmin 2127 real(dp) :: doti,dotr,e0,fac 2128 !arrays 2129 integer :: band_procs(mband) 2130 integer,allocatable :: pwind_tmp(:) 2131 real(dp) :: z1(2) 2132 real(dp),allocatable :: Amat(:,:,:),umat(:,:,:,:),vect1(:,:),vect2(:,:) 2133 2134 ! ************************************************************************* 2135 2136 !calculate 4 matrices ----------------------------- 2137 mpw_tmp=max(mpw,mpw1) 2138 ABI_MALLOC(umat,(2,dtefield%mband_occ,dtefield%mband_occ,4)) 2139 ABI_MALLOC(vect1,(2,0:mpw_tmp)) 2140 ABI_MALLOC(vect2,(2,0:mpw_tmp)) 2141 ABI_MALLOC(pwind_tmp,(mpw_tmp)) 2142 ABI_MALLOC(Amat,(2,dtefield%mband_occ,dtefield%mband_occ)) 2143 vect1(:,0) = zero ; vect2(:,0) = zero 2144 eberry=zero 2145 2146 me_band = mpi_enreg%me_band 2147 2148 !TODO no info on sppol here - I default to isppol 1 2149 do ikpt=1,nkpt 2150 call proc_distrb_band(band_procs,mpi_enreg%proc_distrb,ikpt,1,mband,& 2151 & me_band,mpi_enreg%me_kpt,mpi_enreg%comm_band) 2152 2153 2154 do idir=1,3 2155 2156 fac = dtefield%efield_dot(idir)/& 2157 & (dble(dtefield%nstr(idir))*four_pi) 2158 2159 ! compute <u^(1)_{k_j,q}|u^(1)_{k_j+1,q}> matrix--------- 2160 2161 ! prepare to calculate overlap matrix 2162 ikpt1 = dtefield%ikpt_dk(ikpt,1,idir) 2163 icg = dtefield%cgindex(ikpt,1+nsppol) 2164 icg1 = dtefield%cgindex(ikpt1,1+nsppol) 2165 npw_k1 = npwar1(ikpt) 2166 npw_k2 = npwar1(ikpt1) 2167 pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,3,idir) 2168 vect1(:,0) = zero ; vect2(:,0) = zero 2169 jband_me = 0 2170 do jband = 1, dtefield%mband_occ 2171 if (band_procs(jband) == me_band) then 2172 jband_me = jband_me + 1 2173 vect2(:,1:npw_k2) = cg1(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor) 2174 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 2175 end if 2176 call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr) 2177 2178 2179 iband_me = 0 2180 do iband = 1, dtefield%mband_occ 2181 if (band_procs(iband) /= me_band) cycle 2182 iband_me = iband_me + 1 2183 2184 pwmin = (iband_me-1)*npw_k1*nspinor 2185 pwmax = pwmin + npw_k1*nspinor 2186 vect1(:,1:npw_k1) = & 2187 & cg1(:,icg + 1 + pwmin:icg + pwmax) 2188 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 2189 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 2190 & vect1,vect2) 2191 umat(1,iband,jband,1) = dotr 2192 umat(2,iband,jband,1) = doti 2193 2194 end do ! iband 2195 end do !jband 2196 2197 ! compute <u^(0)_{k_j}|u^(1)_{k_j-n+1,q}> matrix---------------------------------------------------- 2198 2199 ! prepare to calculate overlap matrix 2200 ikpt1 = dtefield%ikpt_dk(ikpt,5,idir) 2201 icg = dtefield%cgindex(ikpt,1) 2202 icg1 = dtefield%cgindex(ikpt1,1+nsppol) 2203 npw_k1 = npwarr(ikpt) 2204 npw_k2 = npwar1(ikpt1) 2205 pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,7,idir) 2206 2207 vect1(:,0) = zero ; vect2(:,0) = zero 2208 jband_me = 0 2209 do jband = 1, dtefield%mband_occ 2210 if (band_procs(jband) == me_band) then 2211 jband_me = jband_me + 1 2212 vect2(:,1:npw_k2) = & 2213 & cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor) 2214 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 2215 end if 2216 call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr) 2217 2218 iband_me = 0 2219 do iband = 1, dtefield%mband_occ 2220 if (band_procs(iband) /= me_band) cycle 2221 iband_me = iband_me + 1 2222 2223 pwmin = (iband_me-1)*npw_k1*nspinor 2224 pwmax = pwmin + npw_k1*nspinor 2225 vect1(:,1:npw_k1) = & 2226 & cg(:,icg + 1 + pwmin:icg + pwmax) 2227 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 2228 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 2229 & vect1,vect2) 2230 umat(1,iband,jband,2) = dotr 2231 umat(2,iband,jband,2) = doti 2232 2233 end do ! iband 2234 end do !jband 2235 2236 ! compute <u^(1)_{k_j-n,q}|u^(0)_{k_j+1}> matrix---------------------------------------------------- 2237 2238 ! prepare to calculate overlap matrix 2239 ikptn = dtefield%ikpt_dk(ikpt,8,idir) 2240 ikpt1 = dtefield%ikpt_dk(ikpt,1,idir) 2241 icg = dtefield%cgindex(ikptn,1+nsppol) 2242 icg1 = dtefield%cgindex(ikpt1,1) 2243 npw_k1 = npwar1(ikptn) 2244 npw_k2 = npwarr(ikpt1) 2245 pwind_tmp(1:npw_k1) = pwindall((ikptn-1)*mpw_tmp+1:(ikptn-1)*mpw_tmp+npw_k1,5,idir) 2246 2247 vect1(:,0) = zero ; vect2(:,0) = zero 2248 jband_me = 0 2249 do jband = 1, dtefield%mband_occ 2250 if (band_procs(jband) == me_band) then 2251 jband_me = jband_me + 1 2252 vect2(:,1:npw_k2) = & 2253 & cg(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor) 2254 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 2255 end if 2256 call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr) 2257 2258 iband_me = 0 2259 do iband = 1, dtefield%mband_occ 2260 if (band_procs(iband) /= me_band) cycle 2261 iband_me = iband_me + 1 2262 2263 pwmin = (iband_me-1)*npw_k1*nspinor 2264 pwmax = pwmin + npw_k1*nspinor 2265 vect1(:,1:npw_k1) = & 2266 & cg1(:,icg + 1 + pwmin:icg + pwmax) 2267 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 2268 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 2269 & vect1,vect2) 2270 umat(1,iband,jband,3) = dotr 2271 umat(2,iband,jband,3) = doti 2272 2273 end do ! iband 2274 end do !jband 2275 2276 ! compute <u^(0)_{-k_j-n+1}|u^(1)_{-k_j,q}> matrix---------------------------------------------------- 2277 2278 ! prepare to calculate overlap matrix 2279 ikptn = dtefield%ikpt_dk(ikpt,5,idir) 2280 ikptnm = dtefield%ikpt_dk(ikptn,9,idir) 2281 ikpt1 = dtefield%ikpt_dk(ikpt,9,idir) 2282 icg = dtefield%cgindex(ikptnm,1) 2283 icg1 = dtefield%cgindex(ikpt1,1+nsppol) 2284 npw_k1 = npwarr(ikptnm) 2285 npw_k2 = npwar1(ikpt1) 2286 pwind_tmp(1:npw_k1) = pwindall((ikptnm-1)*mpw_tmp+1:(ikptnm-1)*mpw_tmp+npw_k1,7,idir) 2287 2288 vect1(:,0) = zero ; vect2(:,0) = zero 2289 jband_me = 0 2290 do jband = 1, dtefield%mband_occ 2291 if (band_procs(jband) == me_band) then 2292 jband_me = jband_me + 1 2293 vect2(:,1:npw_k2) = & 2294 & cg1(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor) 2295 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 2296 end if 2297 call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr) 2298 2299 iband_me = 0 2300 do iband = 1, dtefield%mband_occ 2301 if (band_procs(iband) /= me_band) cycle 2302 iband_me = iband_me + 1 2303 2304 pwmin = (iband_me-1)*npw_k1*nspinor 2305 pwmax = pwmin + npw_k1*nspinor 2306 vect1(:,1:npw_k1) = & 2307 & cg(:,icg + 1 + pwmin:icg + pwmax) 2308 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 2309 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 2310 & vect1,vect2) 2311 umat(1,iband,jband,4) = dotr 2312 umat(2,iband,jband,4) = doti 2313 2314 end do ! iband 2315 end do !jband 2316 2317 ! recompose full umat on each proc, summing over iband 2318 call xmpi_sum(umat,mpi_enreg%comm_band,ierr) 2319 2320 ! sum over the whole------------------------------------------------------------ 2321 2322 e0=zero 2323 do iband=1,dtefield%mband_occ 2324 do jband=1,dtefield%mband_occ 2325 e0 = e0 + 4_dp*(umat(1,iband,jband,1)*qmat(2,jband,iband,ikpt,1,idir)& 2326 & + umat(2,iband,jband,1)*qmat(1,jband,iband,ikpt,1,idir)) 2327 2328 end do 2329 end do 2330 2331 eberry = eberry - e0*fac 2332 2333 e0=zero 2334 2335 ikptn=dtefield%ikpt_dk(ikpt,8,idir) 2336 2337 Amat(:,:,:)=zero 2338 2339 ! calculate Amat 2340 do iband=1, dtefield%mband_occ 2341 do jband=1, dtefield%mband_occ 2342 do kband=1, dtefield%mband_occ 2343 Amat(1,iband,jband) = Amat(1,iband,jband) + (umat(1,iband,kband,3))*& 2344 & qmat(1,kband,jband,ikpt,1,idir)& 2345 & - (umat(2,iband,kband,3))*qmat(2,kband,jband,ikpt,1,idir) 2346 Amat(2,iband,jband) = Amat(2,iband,jband) + (umat(1,iband,kband,3))*& 2347 & qmat(2,kband,jband,ikpt,1,idir)& 2348 & + (umat(2,iband,kband,3))*qmat(1,kband,jband,ikpt,1,idir) 2349 end do 2350 end do 2351 end do 2352 2353 do iband=1, dtefield%mband_occ 2354 do jband=1, dtefield%mband_occ 2355 do kband=1, dtefield%mband_occ 2356 2357 z1(1) = (umat(1,jband,iband,4)+umat(1,iband,jband,2))*& 2358 & qmat(1,jband,kband,ikptn,1,idir)& 2359 & - (umat(2,jband,iband,4)+umat(2,iband,jband,2))*& 2360 & qmat(2,jband,kband,ikptn,1,idir) 2361 z1(2) = (umat(1,jband,iband,4)+umat(1,iband,jband,2))*& 2362 & qmat(2,jband,kband,ikptn,1,idir)& 2363 & + (umat(2,jband,iband,4)+umat(2,iband,jband,2))*& 2364 & qmat(1,jband,kband,ikptn,1,idir) 2365 2366 e0 = e0 - 4_dp*(z1(1)*Amat(2,kband,iband)+z1(2)*Amat(1,kband,iband)) 2367 2368 end do 2369 end do 2370 end do 2371 2372 eberry = eberry - e0*fac 2373 2374 end do !end idir 2375 end do !end ikpt 2376 2377 ABI_FREE(umat) 2378 ABI_FREE(vect1) 2379 ABI_FREE(vect2) 2380 ABI_FREE(pwind_tmp) 2381 ABI_FREE(Amat) 2382 2383 end subroutine dfptff_ebp
ABINIT/dfptff_edie [ Functions ]
NAME
dfptff_edie
FUNCTION
calculate the second order energy from the contribution of \Omega E \cdot P term, Eq.(6) in PRB 75, 115116(2007) [[cite:Wang2007]].
INPUTS
cg(2,mpw*nspinor*mband*mkmem_rbz*nsppol) = planewave coefficients of wavefunctions cg1(2,mpw1*nspinor*mband*mk1mem*nsppol) = pw coefficients of RF wavefunctions at k,q. dtefield = variables related to response Berry-phase calculation ikpt = the index of the current k point isppol = the index of the spin component mband = maximum number of bands mkmem_rbz = maximum number of k-points in core memory mpw = maximum number of plane waves mpw1 = maximum number of plane waves for response wavefunctions nkpt = number of k points npwarr(nkpt) = number of planewaves in basis and boundary at this k point npwar1(nkpt) = number of planewaves in basis and boundary for response wfs nspinor = 1 for scalar wfs, 2 for spinor wfs nsppol = 1 for unpolarized, 2 for spin-polarized qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3) = inverse of the overlap matrix pwindall(max(mpw,mpw1)*mkmem_rbz,8,3) = array used to compute the overlap matrices pwindall(:,1,:) <- <u^(0)_i|u^(0)_i+1> pwindall(:,2,:) <- <u^(0)_i|u^(0)_i-1> pwindall(:,3,:) <- <u^(1)_i|u^(1)_i+1> pwindall(:,4,:) <- <u^(1)_i|u^(1)_i-1> pwindall(:,5,:) <- <u^(1)_i|u^(0)_i+n+1> pwindall(:,6,:) <- <u^(1)_i|u^(0)_i+n-1> pwindall(:,7,:) <- <u^(0)_i|u^(1)_i-n+1> pwindall(:,8,:) <- <u^(0)_i|u^(1)_i-n-1>
OUTPUT
eberry = the energy of the Berry phase term
SOURCE
1713 subroutine dfptff_edie(cg,cg1,dtefield,eberry,idir_efield,mband,mband_mem,mkmem_rbz,& 1714 & mpi_enreg,mpw,mpw1,nkpt,npwarr,npwar1,nsppol,nspinor,pwindall,qmat,rprimd) 1715 1716 !Arguments ---------------------------------------- 1717 !scalars 1718 integer,intent(in) :: idir_efield,mband,mkmem_rbz,mpw,mpw1,nkpt,nspinor,nsppol 1719 integer,intent(in) :: mband_mem 1720 real(dp),intent(out) :: eberry 1721 type(efield_type),intent(in) :: dtefield 1722 type(MPI_type),intent(in) :: mpi_enreg 1723 !arrays 1724 integer,intent(in) :: npwar1(nkpt),npwarr(nkpt) 1725 integer,intent(in) :: pwindall(max(mpw,mpw1)*mkmem_rbz,8,3) 1726 real(dp),intent(in) :: cg(2,mpw*mband_mem*mkmem_rbz*nspinor*nsppol) 1727 real(dp),intent(in) :: cg1(2,mpw1*mband_mem*mkmem_rbz*nspinor*nsppol) 1728 real(dp),intent(in) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3) 1729 real(dp),intent(in) :: rprimd(3,3) 1730 1731 !Local variables ---------------------------------- 1732 !scalars 1733 integer :: iband,icg,icg1,idir 1734 integer :: jband_me, iband_me, me_band, ierr 1735 integer :: ikpt,ikpt1,ikptn,ikptnm 1736 integer :: jband,kband,mpw_tmp,npw_k1,npw_k2,pwmax,pwmin 1737 real(dp) :: doti,dotr,e0,fac 1738 !arrays 1739 integer,allocatable :: pwind_tmp(:) 1740 integer :: band_procs(mband) 1741 real(dp) :: z1(2) 1742 real(dp),allocatable :: Amat(:,:,:),umat(:,:,:,:),vect1(:,:),vect2(:,:) 1743 1744 ! ************************************************************************* 1745 1746 !calculate 4 matrices ----------------------------- 1747 mpw_tmp=max(mpw,mpw1) 1748 ABI_MALLOC(umat,(2,dtefield%mband_occ,dtefield%mband_occ,4)) 1749 ABI_MALLOC(vect1,(2,0:mpw_tmp)) 1750 ABI_MALLOC(vect2,(2,0:mpw_tmp)) 1751 ABI_MALLOC(pwind_tmp,(mpw_tmp)) 1752 ABI_MALLOC(Amat,(2,dtefield%mband_occ,dtefield%mband_occ)) 1753 vect1(:,0) = zero ; vect2(:,0) = zero 1754 eberry=zero 1755 1756 me_band = mpi_enreg%me_band 1757 1758 do ikpt=1,nkpt 1759 1760 call proc_distrb_band(band_procs,mpi_enreg%proc_distrb,ikpt,1,mband,& 1761 & me_band,mpi_enreg%me_kpt,mpi_enreg%comm_band) 1762 1763 do idir=1,3 1764 fac = dtefield%efield_dot(idir)/& 1765 & (dble(dtefield%nstr(idir))*four_pi) 1766 1767 ! compute <u^(1)_{k_j,q}|u^(1)_{k_j+1,q}> matrix---------------------------------------------------- 1768 1769 ! prepare to calculate overlap matrix 1770 ikpt1 = dtefield%ikpt_dk(ikpt,1,idir) 1771 icg = dtefield%cgindex(ikpt,1+nsppol) 1772 icg1 = dtefield%cgindex(ikpt1,1+nsppol) 1773 npw_k1 = npwar1(ikpt) 1774 npw_k2 = npwar1(ikpt1) 1775 pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,3,idir) 1776 vect1(:,0) = zero ; vect2(:,0) = zero 1777 jband_me = 0 1778 do jband = 1, dtefield%mband_occ 1779 if (band_procs(jband) == me_band) then 1780 jband_me = jband_me + 1 1781 vect2(:,1:npw_k2) = & 1782 & cg1(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor) 1783 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 1784 end if 1785 call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr) 1786 1787 iband_me = 0 1788 do iband = 1, dtefield%mband_occ 1789 if (band_procs(iband) /= me_band) cycle 1790 iband_me = iband_me + 1 1791 pwmin = (iband_me-1)*npw_k1*nspinor 1792 pwmax = pwmin + npw_k1*nspinor 1793 vect1(:,1:npw_k1) = & 1794 & cg1(:,icg + 1 + pwmin:icg + pwmax) 1795 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 1796 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 1797 & vect1,vect2) 1798 umat(1,iband,jband,1) = dotr 1799 umat(2,iband,jband,1) = doti 1800 end do ! iband 1801 end do !jband 1802 1803 ! compute <u^(0)_{k_j}|u^(1)_{k_j-n+1,q}> matrix---------------------------------------------------- 1804 1805 ! prepare to calculate overlap matrix 1806 ikpt1 = dtefield%ikpt_dk(ikpt,5,idir) 1807 icg = dtefield%cgindex(ikpt,1) 1808 icg1 = dtefield%cgindex(ikpt1,1+nsppol) 1809 npw_k1 = npwarr(ikpt) 1810 npw_k2 = npwar1(ikpt1) 1811 pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,7,idir) 1812 1813 vect1(:,0) = zero ; vect2(:,0) = zero 1814 jband_me = 0 1815 do jband = 1, dtefield%mband_occ 1816 if (band_procs(jband) == me_band) then 1817 jband_me = jband_me + 1 1818 vect2(:,1:npw_k2) = & 1819 & cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor) 1820 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 1821 end if 1822 call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr) 1823 1824 iband_me = 0 1825 do iband = 1, dtefield%mband_occ 1826 if (band_procs(iband) /= me_band) cycle 1827 iband_me = iband_me + 1 1828 pwmin = (iband_me-1)*npw_k1*nspinor 1829 pwmax = pwmin + npw_k1*nspinor 1830 vect1(:,1:npw_k1) = & 1831 & cg(:,icg + 1 + pwmin:icg + pwmax) 1832 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 1833 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 1834 & vect1,vect2) 1835 umat(1,iband,jband,2) = dotr 1836 umat(2,iband,jband,2) = doti 1837 end do ! iband 1838 end do !jband 1839 1840 ! compute <u^(1)_{k_j-n,q}|u^(0)_{k_j+1}> matrix---------------------------------------------------- 1841 1842 ! prepare to calculate overlap matrix 1843 ikptn = dtefield%ikpt_dk(ikpt,8,idir) 1844 ikpt1 = dtefield%ikpt_dk(ikpt,1,idir) 1845 icg = dtefield%cgindex(ikptn,1+nsppol) 1846 icg1 = dtefield%cgindex(ikpt1,1) 1847 npw_k1 = npwar1(ikptn) 1848 npw_k2 = npwarr(ikpt1) 1849 pwind_tmp(1:npw_k1) = pwindall((ikptn-1)*mpw_tmp+1:(ikptn-1)*mpw_tmp+npw_k1,5,idir) 1850 1851 vect1(:,0) = zero ; vect2(:,0) = zero 1852 jband_me = 0 1853 do jband = 1, dtefield%mband_occ 1854 if (band_procs(jband) == me_band) then 1855 jband_me = jband_me + 1 1856 vect2(:,1:npw_k2) = & 1857 & cg(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor) 1858 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 1859 end if 1860 call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr) 1861 1862 iband_me = 0 1863 do iband = 1, dtefield%mband_occ 1864 if (band_procs(iband) /= me_band) cycle 1865 iband_me = iband_me + 1 1866 pwmin = (iband_me-1)*npw_k1*nspinor 1867 pwmax = pwmin + npw_k1*nspinor 1868 vect1(:,1:npw_k1) = & 1869 & cg1(:,icg + 1 + pwmin:icg + pwmax) 1870 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 1871 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 1872 & vect1,vect2) 1873 umat(1,iband,jband,3) = dotr 1874 umat(2,iband,jband,3) = doti 1875 end do ! iband 1876 end do !jband 1877 1878 ! compute <u^(0)_{-k_j-n+1}|u^(1)_{-k_j,q}> matrix---------------------------------------------------- 1879 1880 ! prepare to calculate overlap matrix 1881 ikptn = dtefield%ikpt_dk(ikpt,5,idir) 1882 ikptnm = dtefield%ikpt_dk(ikptn,9,idir) 1883 ikpt1 = dtefield%ikpt_dk(ikpt,9,idir) 1884 icg = dtefield%cgindex(ikptnm,1) 1885 icg1 = dtefield%cgindex(ikpt1,1+nsppol) 1886 npw_k1 = npwarr(ikptnm) 1887 npw_k2 = npwar1(ikpt1) 1888 pwind_tmp(1:npw_k1) = pwindall((ikptnm-1)*mpw_tmp+1:(ikptnm-1)*mpw_tmp+npw_k1,7,idir) 1889 vect1(:,0) = zero ; vect2(:,0) = zero 1890 jband_me = 0 1891 do jband = 1, dtefield%mband_occ 1892 if (band_procs(jband) == me_band) then 1893 jband_me = jband_me + 1 1894 vect2(:,1:npw_k2) = & 1895 & cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor) 1896 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 1897 end if 1898 call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr) 1899 1900 iband_me = 0 1901 do iband = 1, dtefield%mband_occ 1902 if (band_procs(iband) /= me_band) cycle 1903 iband_me = iband_me + 1 1904 pwmin = (iband_me-1)*npw_k1*nspinor 1905 pwmax = pwmin + npw_k1*nspinor 1906 vect1(:,1:npw_k1) = & 1907 & cg(:,icg + 1 + pwmin:icg + pwmax) 1908 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 1909 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 1910 & vect1,vect2) 1911 umat(1,iband,jband,4) = dotr 1912 umat(2,iband,jband,4) = doti 1913 end do ! iband 1914 end do !jband 1915 1916 ! sum up over all iband and all procs 1917 call xmpi_sum(umat,mpi_enreg%comm_band,ierr) 1918 1919 ! sum over the whole------------------------------------------------------------ 1920 1921 e0=zero 1922 do iband=1,dtefield%mband_occ 1923 do jband=1,dtefield%mband_occ 1924 e0 = e0 + 2_dp*(umat(1,iband,jband,1)*qmat(2,jband,iband,ikpt,1,idir)& 1925 & + umat(2,iband,jband,1)*qmat(1,jband,iband,ikpt,1,idir)) 1926 end do 1927 end do 1928 eberry = eberry - e0*fac 1929 e0=zero 1930 ikptn=dtefield%ikpt_dk(ikpt,8,idir) 1931 Amat(:,:,:)=zero 1932 ! calculate Amat 1933 do iband=1, dtefield%mband_occ 1934 do jband=1, dtefield%mband_occ 1935 do kband=1, dtefield%mband_occ 1936 Amat(1,iband,jband) = Amat(1,iband,jband) + (umat(1,iband,kband,3))*qmat(1,kband,jband,ikpt,1,idir)& 1937 & - (umat(2,iband,kband,3))*qmat(2,kband,jband,ikpt,1,idir) 1938 Amat(2,iband,jband) = Amat(2,iband,jband) + (umat(1,iband,kband,3))*qmat(2,kband,jband,ikpt,1,idir)& 1939 & + (umat(2,iband,kband,3))*qmat(1,kband,jband,ikpt,1,idir) 1940 end do 1941 end do 1942 end do 1943 1944 do iband=1, dtefield%mband_occ 1945 do jband=1, dtefield%mband_occ 1946 do kband=1, dtefield%mband_occ 1947 z1(1) = (umat(1,jband,iband,4)+umat(1,iband,jband,2))*qmat(1,jband,kband,ikptn,1,idir)& 1948 & - (umat(2,jband,iband,4)+umat(2,iband,jband,2))*qmat(2,jband,kband,ikptn,1,idir) 1949 z1(2) = (umat(1,jband,iband,4)+umat(1,iband,jband,2))*qmat(2,jband,kband,ikptn,1,idir)& 1950 & + (umat(2,jband,iband,4)+umat(2,iband,jband,2))*qmat(1,jband,kband,ikptn,1,idir) 1951 1952 e0 = e0 - 4_dp*(z1(1)*Amat(2,kband,iband)+z1(2)*Amat(1,kband,iband)) 1953 end do 1954 end do 1955 end do 1956 1957 eberry = eberry - e0*fac 1958 1959 ! !---------------------------------last part--------------------------------------------- 1960 1961 fac = rprimd(idir_efield,idir)/& 1962 & (dble(dtefield%nstr(idir))*two_pi) 1963 1964 ! compute <u^(1)_{k_j-n,q}|u^(0)_{k_j+1}> matrix---------------------------------------------------- 1965 1966 ! prepare to calculate overlap matrix 1967 ikptn = dtefield%ikpt_dk(ikpt,8,idir) 1968 ikpt1 = dtefield%ikpt_dk(ikpt,1,idir) 1969 icg = dtefield%cgindex(ikptn,1+nsppol) 1970 icg1 = dtefield%cgindex(ikpt1,1) 1971 npw_k1 = npwar1(ikptn) 1972 npw_k2 = npwarr(ikpt1) 1973 pwind_tmp(1:npw_k1) = pwindall((ikptn-1)*mpw_tmp+1:(ikptn-1)*mpw_tmp+npw_k1,5,idir) 1974 1975 vect1(:,0) = zero ; vect2(:,0) = zero 1976 jband_me = 0 1977 do jband = 1, dtefield%mband_occ 1978 if (band_procs(jband) == me_band) then 1979 jband_me = jband_me + 1 1980 vect2(:,1:npw_k2) = & 1981 & cg(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor) 1982 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 1983 end if 1984 call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr) 1985 1986 iband_me = 0 1987 do iband = 1, dtefield%mband_occ 1988 if (band_procs(iband) /= me_band) cycle 1989 iband_me = iband_me + 1 1990 pwmin = (iband_me-1)*npw_k1*nspinor 1991 pwmax = pwmin + npw_k1*nspinor 1992 vect1(:,1:npw_k1) = & 1993 & cg1(:,icg + 1 + pwmin:icg + pwmax) 1994 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 1995 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 1996 & vect1,vect2) 1997 umat(1,iband,jband,1) = dotr 1998 umat(2,iband,jband,1) = doti 1999 end do ! iband 2000 end do !jband 2001 2002 ! compute <u^(0)_{k_j}|u^(1)_{k_j-n+1,q}> matrix---------------------------------------------------- 2003 2004 ! prepare to calculate overlap matrix 2005 ikpt1 = dtefield%ikpt_dk(ikpt,5,idir) 2006 icg = dtefield%cgindex(ikpt,1) 2007 icg1 = dtefield%cgindex(ikpt1,1+nsppol) 2008 npw_k1 = npwarr(ikpt) 2009 npw_k2 = npwar1(ikpt1) 2010 pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,7,idir) 2011 vect1(:,0) = zero ; vect2(:,0) = zero 2012 jband_me = 0 2013 do jband = 1, dtefield%mband_occ 2014 if (band_procs(jband) == me_band) then 2015 jband_me = jband_me + 1 2016 vect2(:,1:npw_k2) = & 2017 & cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor) 2018 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 2019 end if 2020 call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr) 2021 2022 iband_me = 0 2023 do iband = 1, dtefield%mband_occ 2024 if (band_procs(iband) /= me_band) cycle 2025 iband_me = iband_me + 1 2026 pwmin = (iband_me-1)*npw_k1*nspinor 2027 pwmax = pwmin + npw_k1*nspinor 2028 vect1(:,1:npw_k1) = & 2029 & cg(:,icg + 1 + pwmin:icg + pwmax) 2030 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 2031 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 2032 & vect1,vect2) 2033 umat(1,iband,jband,1) = umat(1,iband,jband,1) + dotr 2034 umat(2,iband,jband,1) = umat(2,iband,jband,1) + doti 2035 end do ! iband 2036 end do !jband 2037 2038 ! sum up over all iband and all procs 2039 call xmpi_sum(umat,mpi_enreg%comm_band,ierr) 2040 2041 e0=zero 2042 2043 do iband=1,dtefield%mband_occ 2044 do jband=1,dtefield%mband_occ 2045 e0 = e0 + (umat(1,iband,jband,1)*qmat(2,jband,iband,ikpt,1,idir)& 2046 & + umat(2,iband,jband,1)*qmat(1,jband,iband,ikpt,1,idir)) 2047 end do 2048 end do 2049 2050 eberry = eberry - e0*fac 2051 2052 end do !end idir 2053 end do !end ikpt 2054 2055 ABI_FREE(umat) 2056 ABI_FREE(vect1) 2057 ABI_FREE(vect2) 2058 ABI_FREE(pwind_tmp) 2059 ABI_FREE(Amat) 2060 2061 end subroutine dfptff_edie
ABINIT/dfptff_gbefd [ Functions ]
NAME
dfptff_gbefd
FUNCTION
calculate the gradient of the second order \Omega E \cdot P term, Eq.(23) in PRB 75, 115116(2007) [[cite:Wang2007]].
INPUTS
cg(2,mpw*nspinor*mband*mkmem_rbz*nsppol) = planewave coefficients of wavefunctions cg1(2,mpw1*nspinor*mband*mk1mem*nsppol) = pw coefficients of RF wavefunctions at k,q. dtefield = variables related to response Berry-phase calculation ikpt = the index of the current k point isppol = the index of the spin component mband = maximum number of bands mband_mem = maximum number of bands on this cpu mkmem_rbz = maximum number of k-points in core memory mpi_enreg = parallel distribution datastructure mpw = maximum number of plane waves mpw1 = maximum number of plane waves for response wavefunctions nkpt = number of k points npwarr(nkpt) = number of planewaves in basis and boundary at this k point npwar1(nkpt) = number of planewaves in basis and boundary for response wfs nspinor = 1 for scalar wfs, 2 for spinor wfs nsppol = 1 for unpolarized, 2 for spin-polarized qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3) = inverse of the overlap matrix pwindall(max(mpw,mpw1)*mkmem_rbz,8,3) = array used to compute the overlap matrices pwindall(:,1,:) <- <u^(0)_i|u^(0)_i+1> pwindall(:,2,:) <- <u^(0)_i|u^(0)_i-1> pwindall(:,3,:) <- <u^(1)_i|u^(1)_i+1> pwindall(:,4,:) <- <u^(1)_i|u^(1)_i-1> pwindall(:,5,:) <- <u^(1)_i|u^(0)_i+n+1> pwindall(:,6,:) <- <u^(1)_i|u^(0)_i+n-1> pwindall(:,7,:) <- <u^(0)_i|u^(1)_i-n+1> pwindall(:,8,:) <- <u^(0)_i|u^(1)_i-n-1>
OUTPUT
grad_berry = the gradient of the Berry phase term
SOURCE
1218 subroutine dfptff_gbefd(cg,cg1,dtefield,grad_berry,idir_efield,ikpt,isppol,& 1219 & mband,mband_mem,mpw,mpw1,mkmem_rbz,mk1mem,& 1220 & mpi_enreg,nkpt,npwarr,npwar1,nspinor,nsppol,qmat,pwindall,rprimd) 1221 1222 !Arguments ---------------------------------------- 1223 !scalars 1224 integer,intent(in) :: idir_efield,ikpt,isppol,mband,mk1mem,mkmem_rbz,mpw,mpw1,nkpt 1225 integer,intent(in) :: mband_mem 1226 integer,intent(in) :: nspinor,nsppol 1227 type(efield_type),intent(in) :: dtefield 1228 type(MPI_type),intent(in) :: mpi_enreg 1229 !arrays 1230 integer,intent(in) :: npwar1(nkpt),npwarr(nkpt) 1231 integer,intent(in) :: pwindall(max(mpw,mpw1)*mkmem_rbz,8,3) 1232 real(dp),intent(in) :: cg(2,mpw*nspinor*mband_mem*mkmem_rbz*nsppol) 1233 real(dp),intent(in) :: cg1(2,mpw1*nspinor*mband_mem*mk1mem*nsppol) 1234 real(dp),intent(in) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3) 1235 real(dp),intent(in) :: rprimd(3,3) 1236 !TODO MJV: this array is still npw*nband, and should be parallelized as well at some point... 1237 real(dp),intent(out) :: grad_berry(2,mpw1,dtefield%mband_occ) 1238 1239 !Local variables ------------------------- 1240 !scalars 1241 integer :: iband,icg,icg1,idir,ikpt1 1242 integer :: iband_me, jband_me, ierr 1243 integer :: ikptn,ikptnp1,ipw,jband,jpw,kband 1244 integer :: mpw_tmp,npw_k1,npw_k2,pwmax,pwmin 1245 real(dp) :: doti,dotr,fac,wfi,wfr 1246 !arrays 1247 integer,allocatable :: pwind_tmp(:) 1248 integer :: band_procs(mband) 1249 real(dp) :: z1(2),z2(2) 1250 real(dp),allocatable :: Amat(:,:,:),Bmat(:,:,:),s1mat(:,:,:),vect1(:,:) 1251 real(dp),allocatable :: vect2(:,:) 1252 1253 ! ************************************************************************* 1254 1255 mpw_tmp=max(mpw,mpw1) 1256 ABI_MALLOC(vect1,(2,0:mpw_tmp)) 1257 ABI_MALLOC(vect2,(2,0:mpw_tmp)) 1258 ABI_MALLOC(s1mat,(2,dtefield%mband_occ,dtefield%mband_occ)) 1259 ABI_MALLOC(pwind_tmp,(mpw_tmp)) 1260 ABI_MALLOC(Amat,(2,dtefield%mband_occ,dtefield%mband_occ)) 1261 ABI_MALLOC(Bmat,(2,dtefield%mband_occ,dtefield%mband_occ)) 1262 vect1(:,0) = zero ; vect2(:,0) = zero 1263 s1mat(:,:,:)=zero 1264 grad_berry(:,:,:) = zero 1265 1266 call proc_distrb_band(band_procs,mpi_enreg%proc_distrb,ikpt,isppol,mband,& 1267 & mpi_enreg%me_band,mpi_enreg%me_kpt,mpi_enreg%comm_band) 1268 1269 do idir=1,3 1270 fac = dtefield%efield_dot(idir)*dble(nkpt)/& 1271 & (dble(dtefield%nstr(idir))*four_pi) 1272 1273 ! prepare 1274 ikpt1 = dtefield%ikpt_dk(ikpt,1,idir) 1275 icg1 = dtefield%cgindex(ikpt1,isppol+nsppol) 1276 npw_k1 = npwar1(ikpt) 1277 npw_k2 = npwar1(ikpt1) 1278 pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,3,idir) 1279 1280 do ipw = 1, npw_k1 1281 jpw = pwind_tmp(ipw) 1282 if (jpw > 0) then 1283 iband_me = 0 1284 do iband = 1, dtefield%mband_occ 1285 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,isppol,mpi_enreg%me_band)) cycle 1286 iband_me = iband_me + 1 1287 wfr = cg1(1,icg1 + (iband_me - 1)*npw_k2*nspinor + jpw) 1288 wfi = cg1(2,icg1 + (iband_me - 1)*npw_k2*nspinor + jpw) 1289 1290 do jband = 1, dtefield%mband_occ 1291 grad_berry(1,ipw,jband) = & 1292 & grad_berry(1,ipw,jband) + & 1293 & fac*qmat(1,iband,jband,ikpt,1,idir)*wfr - fac*qmat(2,iband,jband,ikpt,1,idir)*wfi 1294 1295 grad_berry(2,ipw,jband) = & 1296 & grad_berry(2,ipw,jband) + & 1297 & fac*qmat(1,iband,jband,ikpt,1,idir)*wfi + fac*qmat(2,iband,jband,ikpt,1,idir)*wfr 1298 end do 1299 end do 1300 end if 1301 end do 1302 ! accumulate full sum over iband below at the end 1303 1304 ! compute <u^(0)_{k_j}|u^(1)_{k_j+1}> matrix---------------------------------------------------- 1305 1306 ! prepare to calculate overlap matrix 1307 ikptn = ikpt 1308 ikpt1 = dtefield%ikpt_dk(ikpt,1,idir) 1309 icg = dtefield%cgindex(ikptn,isppol) 1310 icg1 = dtefield%cgindex(ikpt1,isppol+nsppol) 1311 npw_k1 = npwarr(ikptn) 1312 npw_k2 = npwar1(ikpt1) 1313 pwind_tmp(1:npw_k1) = pwindall((ikptn-1)*mpw_tmp+1:(ikptn-1)*mpw_tmp+npw_k1,7,idir) 1314 1315 vect1(:,0) = zero ; vect2(:,0) = zero 1316 jband_me = 0 1317 do jband = 1, dtefield%mband_occ 1318 if (band_procs(jband) == mpi_enreg%me_band) then 1319 jband_me = jband_me + 1 1320 vect2(:,1:npw_k2) = & 1321 & cg1(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor) 1322 end if 1323 call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr) 1324 1325 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 1326 iband_me = 0 1327 do iband = 1, dtefield%mband_occ 1328 if (band_procs(iband) /= mpi_enreg%me_band) cycle 1329 iband_me = iband_me + 1 1330 pwmin = (iband_me-1)*npw_k1*nspinor 1331 pwmax = pwmin + npw_k1*nspinor 1332 vect1(:,1:npw_k1) = cg(:,icg + 1 + pwmin:icg + pwmax) 1333 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 1334 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,vect1,vect2) 1335 s1mat(1,iband,jband) = dotr 1336 s1mat(2,iband,jband) = doti 1337 end do ! iband 1338 end do !jband 1339 1340 ! compute <u^(1)_{k_j}|u^(0)_{k_j+1}> matrix----------------------------------------------------- 1341 1342 ! prepare to calculate overlap matrix 1343 ikpt1 = dtefield%ikpt_dk(ikpt,1,idir) 1344 icg = dtefield%cgindex(ikpt,isppol+nsppol) 1345 icg1 = dtefield%cgindex(ikpt1,isppol) 1346 npw_k1 = npwar1(ikpt) 1347 npw_k2 = npwarr(ikpt1) 1348 pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,5,idir) 1349 1350 vect1(:,0) = zero ; vect2(:,0) = zero 1351 jband_me = 0 1352 do jband = 1, dtefield%mband_occ 1353 if (band_procs(jband) == mpi_enreg%me_band) then 1354 jband_me = jband_me + 1 1355 vect2(:,1:npw_k2) = & 1356 & cg(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor) 1357 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 1358 end if 1359 call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr) 1360 1361 iband_me = 0 1362 do iband = 1, dtefield%mband_occ 1363 if (band_procs(iband) /= mpi_enreg%me_band) cycle 1364 iband_me = iband_me + 1 1365 pwmin = (iband_me-1)*npw_k1*nspinor 1366 pwmax = pwmin + npw_k1*nspinor 1367 vect1(:,1:npw_k1) = & 1368 & cg1(:,icg + 1 + pwmin:icg + pwmax) 1369 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 1370 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,vect1,vect2) 1371 s1mat(1,jband,iband) = s1mat(1,jband,iband) + dotr 1372 s1mat(2,jband,iband) = s1mat(2,jband,iband) + doti 1373 end do ! iband 1374 end do !jband 1375 ! accumulate all iband values on each proc 1376 call xmpi_sum(s1mat,mpi_enreg%comm_band,ierr) 1377 1378 Amat(:,:,:)=zero 1379 1380 ! calculate Amat 1381 do iband=1, dtefield%mband_occ 1382 do jband=1, dtefield%mband_occ 1383 do kband=1, dtefield%mband_occ 1384 Amat(1,iband,jband) = Amat(1,iband,jband) + s1mat(1,iband,kband)*qmat(1,kband,jband,ikpt,1,idir)& 1385 & - s1mat(2,iband,kband)*qmat(2,kband,jband,ikpt,1,idir) 1386 Amat(2,iband,jband) = Amat(2,iband,jband) + s1mat(1,iband,kband)*qmat(2,kband,jband,ikpt,1,idir)& 1387 & + s1mat(2,iband,kband)*qmat(1,kband,jband,ikpt,1,idir) 1388 end do 1389 end do 1390 end do 1391 1392 Bmat(:,:,:)=zero 1393 1394 ! calculate Bmat 1395 ikptn = dtefield%ikpt_dk(ikpt,7,idir) 1396 do iband=1, dtefield%mband_occ 1397 do jband=1, dtefield%mband_occ 1398 do kband=1, dtefield%mband_occ 1399 Bmat(1,jband,kband) = Bmat(1,jband,kband) + Amat(1,iband,kband)*qmat(1,jband,iband,ikptn,1,idir)& 1400 & - Amat(2,iband,kband)*qmat(2,jband,iband,ikptn,1,idir) 1401 Bmat(2,jband,kband) = Bmat(2,jband,kband) + Amat(1,iband,kband)*qmat(2,jband,iband,ikptn,1,idir)& 1402 & + Amat(2,iband,kband)*qmat(1,jband,iband,ikptn,1,idir) 1403 end do 1404 end do 1405 end do 1406 1407 ! calc. the second term of gradient------------------------------ 1408 1409 ! preparation 1410 1411 ikptnp1 = dtefield%ikpt_dk(ikpt,3,idir) 1412 icg = dtefield%cgindex(ikptnp1,isppol) 1413 npw_k1 = npwar1(ikpt) 1414 npw_k2 = npwarr(ikptnp1) 1415 pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,5,idir) 1416 1417 z1(:) = zero 1418 z2(:) = zero 1419 1420 do ipw = 1, npw_k1 1421 jpw = pwind_tmp(ipw) 1422 if (jpw > 0) then 1423 iband_me = 0 1424 do iband = 1, dtefield%mband_occ 1425 if (band_procs(iband) /= mpi_enreg%me_band) cycle 1426 iband_me = iband_me + 1 1427 wfr = cg(1,icg + (iband_me - 1)*npw_k2*nspinor + jpw) 1428 wfi = cg(2,icg + (iband_me - 1)*npw_k2*nspinor + jpw) 1429 1430 do jband=1, dtefield%mband_occ 1431 grad_berry(1,ipw,jband) = grad_berry(1,ipw,jband) - fac*(Bmat(1,iband,jband)*wfr - Bmat(2,iband,jband)*wfi) 1432 grad_berry(2,ipw,jband) = grad_berry(2,ipw,jband) - fac*(Bmat(1,iband,jband)*wfi + Bmat(2,iband,jband)*wfr) 1433 end do 1434 end do 1435 end if 1436 end do 1437 1438 ! Second part of gradient of Berry phase++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1439 1440 vect1(:,0) = zero ; vect2(:,0) = zero 1441 1442 ! prepare 1443 ikpt1 = dtefield%ikpt_dk(ikpt,2,idir) 1444 icg1 = dtefield%cgindex(ikpt1,isppol+nsppol) 1445 npw_k1 = npwar1(ikpt) 1446 npw_k2 = npwar1(ikpt1) 1447 pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,4,idir) 1448 1449 do ipw = 1, npw_k1 1450 jpw = pwind_tmp(ipw) 1451 if (jpw > 0) then 1452 iband_me = 0 1453 do iband = 1, dtefield%mband_occ 1454 if (band_procs(iband) /= mpi_enreg%me_band) cycle 1455 iband_me = iband_me + 1 1456 wfr = cg1(1,icg1 + (iband_me - 1)*npw_k2*nspinor + jpw) 1457 wfi = cg1(2,icg1 + (iband_me - 1)*npw_k2*nspinor + jpw) 1458 1459 do jband = 1, dtefield%mband_occ 1460 grad_berry(1,ipw,jband) = & 1461 & grad_berry(1,ipw,jband) - & 1462 & fac*qmat(1,iband,jband,ikpt,2,idir)*wfr + fac*qmat(2,iband,jband,ikpt,2,idir)*wfi 1463 1464 grad_berry(2,ipw,jband) = & 1465 & grad_berry(2,ipw,jband) - & 1466 & fac*qmat(1,iband,jband,ikpt,2,idir)*wfi - fac*qmat(2,iband,jband,ikpt,2,idir)*wfr 1467 end do 1468 end do 1469 end if 1470 end do 1471 1472 ! compute <u^(0)_{k_j}|u^(1)_{k_j-1}> matrix---------------------------------------------------- 1473 1474 ! prepare to calculate overlap matrix 1475 ikptn = ikpt 1476 ikpt1 = dtefield%ikpt_dk(ikpt,2,idir) 1477 icg = dtefield%cgindex(ikptn,isppol) 1478 icg1 = dtefield%cgindex(ikpt1,isppol+nsppol) 1479 npw_k1 = npwarr(ikptn) 1480 npw_k2 = npwar1(ikpt1) 1481 pwind_tmp(1:npw_k1) =pwindall((ikptn-1)*mpw_tmp+1:(ikptn-1)*mpw_tmp+npw_k1,8,idir) 1482 1483 vect1(:,0) = zero ; vect2(:,0) = zero 1484 jband_me = 0 1485 do jband = 1, dtefield%mband_occ 1486 if (band_procs(jband) == mpi_enreg%me_band) then 1487 jband_me = jband_me + 1 1488 vect2(:,1:npw_k2) = & 1489 & cg1(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor) 1490 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 1491 end if 1492 call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr) 1493 1494 1495 iband_me = 0 1496 do iband = 1, dtefield%mband_occ 1497 if (band_procs(iband) /= mpi_enreg%me_band) cycle 1498 iband_me = iband_me + 1 1499 pwmin = (iband_me-1)*npw_k1*nspinor 1500 pwmax = pwmin + npw_k1*nspinor 1501 vect1(:,1:npw_k1) = & 1502 & cg(:,icg + 1 + pwmin:icg + pwmax) 1503 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 1504 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 1505 & vect1,vect2) 1506 s1mat(1,iband,jband) = dotr 1507 s1mat(2,iband,jband) = doti 1508 end do ! iband 1509 end do !jband 1510 1511 ! compute <u^(1)_{k_j}|u^(0)_{k_j-1}> matrix----------------------------------------------------- 1512 1513 ! prepare to calculate overlap matrix 1514 ikpt1 = dtefield%ikpt_dk(ikpt,2,idir) 1515 icg = dtefield%cgindex(ikpt,isppol+nsppol) 1516 icg1 = dtefield%cgindex(ikpt1,isppol) 1517 npw_k1 = npwarr(ikpt) 1518 npw_k2 = npwar1(ikpt1) 1519 pwind_tmp(1:npw_k1) =pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,6,idir) 1520 1521 vect1(:,0) = zero ; vect2(:,0) = zero 1522 jband_me = 0 1523 do jband = 1, dtefield%mband_occ 1524 if (band_procs(jband) == mpi_enreg%me_band) then 1525 jband_me = jband_me + 1 1526 vect2(:,1:npw_k2) = & 1527 & cg(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor) 1528 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 1529 end if 1530 call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr) 1531 1532 do iband = 1, dtefield%mband_occ 1533 pwmin = (iband-1)*npw_k1*nspinor 1534 pwmax = pwmin + npw_k1*nspinor 1535 vect1(:,1:npw_k1) = & 1536 & cg1(:,icg + 1 + pwmin:icg + pwmax) 1537 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 1538 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 1539 & vect1,vect2) 1540 s1mat(1,jband,iband) = s1mat(1,jband,iband) + dotr 1541 s1mat(2,jband,iband) = s1mat(2,jband,iband) + doti 1542 end do ! iband 1543 end do !jband 1544 1545 Amat(:,:,:)=zero 1546 1547 ! calculate Amat 1548 do iband=1, dtefield%mband_occ 1549 do jband=1, dtefield%mband_occ 1550 do kband=1, dtefield%mband_occ 1551 Amat(1,iband,jband) = Amat(1,iband,jband) + s1mat(1,iband,kband)*qmat(1,kband,jband,ikpt,2,idir)& 1552 & - s1mat(2,iband,kband)*qmat(2,kband,jband,ikpt,2,idir) 1553 Amat(2,iband,jband) = Amat(2,iband,jband) + s1mat(1,iband,kband)*qmat(2,kband,jband,ikpt,2,idir)& 1554 & + s1mat(2,iband,kband)*qmat(1,kband,jband,ikpt,2,idir) 1555 end do 1556 end do 1557 end do 1558 1559 Bmat(:,:,:)=zero 1560 1561 ! calculate Bmat 1562 ikptn = dtefield%ikpt_dk(ikpt,7,idir) 1563 do iband=1, dtefield%mband_occ 1564 do jband=1, dtefield%mband_occ 1565 do kband=1, dtefield%mband_occ 1566 Bmat(1,jband,kband) = Bmat(1,jband,kband) + Amat(1,iband,kband)*qmat(1,jband,iband,ikptn,2,idir)& 1567 & - Amat(2,iband,kband)*qmat(2,jband,iband,ikptn,2,idir) 1568 Bmat(2,jband,kband) = Bmat(2,jband,kband) + Amat(1,iband,kband)*qmat(2,jband,iband,ikptn,2,idir)& 1569 + Amat(2,iband,kband)*qmat(1,jband,iband,ikptn,2,idir) 1570 end do 1571 end do 1572 end do 1573 1574 ! calc. the second term of gradient------------------------------ 1575 1576 ! preparation 1577 1578 ikptnp1 = dtefield%ikpt_dk(ikpt,4,idir) 1579 icg = dtefield%cgindex(ikptnp1,isppol) 1580 npw_k1 = npwar1(ikpt) 1581 npw_k2 = npwarr(ikptnp1) 1582 pwind_tmp(1:npw_k1) =pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,6,idir) 1583 z1(:) = zero 1584 z2(:) = zero 1585 do ipw = 1, npw_k1 1586 jpw = pwind_tmp(ipw) 1587 if (jpw > 0) then 1588 do iband = 1, dtefield%mband_occ 1589 wfr = cg(1,icg + (iband - 1)*npw_k2*nspinor + jpw) 1590 wfi = cg(2,icg + (iband - 1)*npw_k2*nspinor + jpw) 1591 do jband=1, dtefield%mband_occ 1592 grad_berry(1,ipw,jband) = grad_berry(1,ipw,jband) + fac*(Bmat(1,iband,jband)*wfr - Bmat(2,iband,jband)*wfi) 1593 grad_berry(2,ipw,jband) = grad_berry(2,ipw,jband) + fac*(Bmat(1,iband,jband)*wfi + Bmat(2,iband,jband)*wfr) 1594 end do 1595 end do 1596 end if 1597 end do 1598 1599 end do !idir 1600 1601 !!----------------------------------------third part of gradient------------------------------------------------------ 1602 do idir=1,3 1603 fac = rprimd(idir_efield,idir)*dble(nkpt)/& 1604 & (dble(dtefield%nstr(idir))*four_pi) 1605 1606 ! prepare 1607 ikpt1 = dtefield%ikpt_dk(ikpt,1,idir) 1608 icg1 = dtefield%cgindex(ikpt1,isppol) 1609 npw_k1 = npwar1(ikpt) 1610 npw_k2 = npwarr(ikpt1) 1611 pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,5,idir) 1612 do ipw = 1, npw_k1 1613 jpw = pwind_tmp(ipw) 1614 if (jpw > 0) then 1615 do iband = 1, dtefield%mband_occ 1616 wfr = cg(1,icg1 + (iband - 1)*npw_k2*nspinor + jpw) 1617 wfi = cg(2,icg1 + (iband - 1)*npw_k2*nspinor + jpw) 1618 do jband = 1, dtefield%mband_occ 1619 grad_berry(1,ipw,jband) = & 1620 & grad_berry(1,ipw,jband) + & 1621 & fac*qmat(1,iband,jband,ikpt,1,idir)*wfr - fac*qmat(2,iband,jband,ikpt,1,idir)*wfi 1622 grad_berry(2,ipw,jband) = & 1623 & grad_berry(2,ipw,jband) + & 1624 & fac*qmat(1,iband,jband,ikpt,1,idir)*wfi + fac*qmat(2,iband,jband,ikpt,1,idir)*wfr 1625 end do 1626 end do 1627 end if 1628 end do 1629 1630 ! prepare 1631 ikpt1 = dtefield%ikpt_dk(ikpt,2,idir) 1632 icg1 = dtefield%cgindex(ikpt1,isppol) 1633 npw_k1 = npwar1(ikpt) 1634 npw_k2 = npwarr(ikpt1) 1635 pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,6,idir) 1636 1637 do ipw = 1, npw_k1 1638 jpw = pwind_tmp(ipw) 1639 if (jpw > 0) then 1640 do iband = 1, dtefield%mband_occ 1641 wfr = cg(1,icg1 + (iband - 1)*npw_k2*nspinor + jpw) 1642 wfi = cg(2,icg1 + (iband - 1)*npw_k2*nspinor + jpw) 1643 do jband = 1, dtefield%mband_occ 1644 grad_berry(1,ipw,jband) = & 1645 & grad_berry(1,ipw,jband) - & 1646 & fac*qmat(1,iband,jband,ikpt,2,idir)*wfr + fac*qmat(2,iband,jband,ikpt,2,idir)*wfi 1647 1648 grad_berry(2,ipw,jband) = & 1649 & grad_berry(2,ipw,jband) - & 1650 & fac*qmat(1,iband,jband,ikpt,2,idir)*wfi - fac*qmat(2,iband,jband,ikpt,2,idir)*wfr 1651 1652 end do 1653 end do 1654 end if 1655 end do 1656 1657 end do !idir 1658 1659 ! accumulate full sum over iband in each jband entry 1660 call xmpi_sum(grad_berry,mpi_enreg%comm_band,ierr) 1661 1662 ABI_FREE(vect1) 1663 ABI_FREE(vect2) 1664 ABI_FREE(s1mat) 1665 ABI_FREE(Amat) 1666 ABI_FREE(Bmat) 1667 ABI_FREE(pwind_tmp) 1668 1669 end subroutine dfptff_gbefd
ABINIT/dfptff_gradberry [ Functions ]
NAME
dfptff_gradberry
FUNCTION
Calculation of the gradient of Berry-phase term in finite electric field.
COPYRIGHT
Copyright (C) 2004-2024 ABINIT group (XW). 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 .
INPUTS
cg(2,mpw*nspinor*mband*mkmem_rbz*nsppol) = planewave coefficients of wavefunctions cg1(2,mpw1*nspinor*mband*mk1mem*nsppol) = pw coefficients of RF wavefunctions at k,q. dtefield = variables related to finite electric field calculation ikpt = the index of the current k point isppol=1 for unpolarized, 2 for spin-polarized mband = maximum number of bands mkmem_rbz = maximum number of k-points in core memory mpw = maximum number of plane waves mpw1 = maximum number of plane waves for response wavefunctions nkpt = number of k points npwarr(nkpt) = number of planewaves in basis and boundary at this k point npwar1(nkpt) = number of planewaves in basis and boundary for response wfs nspinor = 1 for scalar wfs, 2 for spinor wfs nsppol = 1 for unpolarized, 2 for spin-polarized qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3) = inverse of the overlap matrix pwindall(max(mpw,mpw1)*mkmem_rbz,8,3) = array used to compute the overlap matrices pwindall(:,1,:) <- <u^(0)_i|u^(0)_i+1> pwindall(:,2,:) <- <u^(0)_i|u^(0)_i-1> pwindall(:,3,:) <- <u^(1)_i|u^(1)_i+1> pwindall(:,4,:) <- <u^(1)_i|u^(1)_i-1> pwindall(:,5,:) <- <u^(1)_i|u^(0)_i+n+1> pwindall(:,6,:) <- <u^(1)_i|u^(0)_i+n-1> pwindall(:,7,:) <- <u^(0)_i|u^(1)_i-n+1> pwindall(:,8,:) <- <u^(0)_i|u^(1)_i-n-1>
OUTPUT
grad_berry(2,mpw1,dtefield%mband_occ) = the gradient of the Berry phase term
SOURCE
735 subroutine dfptff_gradberry(cg,cg1,dtefield,grad_berry,ikpt,isppol,& 736 & mband,mband_mem,mpw,mpw1,mkmem_rbz,mk1mem,& 737 & mpi_enreg,nkpt,& 738 & npwarr,npwar1,nspinor,nsppol,qmat,pwindall) 739 740 !Arguments ---------------------------------------- 741 !scalars 742 integer,intent(in) :: ikpt,isppol,mband,mk1mem,mkmem_rbz,mpw,mpw1,nkpt,nspinor 743 integer,intent(in) :: mband_mem 744 integer,intent(in) :: nsppol 745 type(efield_type),intent(in) :: dtefield 746 type(MPI_type),intent(in) :: mpi_enreg 747 !arrays 748 integer,intent(in) :: npwar1(nkpt),npwarr(nkpt) 749 integer,intent(in) :: pwindall(max(mpw,mpw1)*mkmem_rbz,8,3) 750 real(dp),intent(in) :: cg(2,mpw*nspinor*mband_mem*mkmem_rbz*nsppol) 751 real(dp),intent(in) :: cg1(2,mpw1*nspinor*mband_mem*mk1mem*nsppol) 752 real(dp),intent(in) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3) 753 !TODO MJV: grad_berry should also be dimensioned with mband_mem 754 real(dp),intent(out) :: grad_berry(2,mpw1,dtefield%mband_occ) 755 756 !Local variables ------------------------- 757 !scalars 758 integer :: iband,icg,icg1,idir,ikpt1 759 integer :: iband_me, jband_me, me_band, ierr 760 integer :: ikpt1m,ikptn,ikptnm,ikptnp1,ipw,jband,jpw,kband 761 integer :: mpw_tmp,npw_k1,npw_k2,pwmax,pwmin 762 real(dp) :: doti,dotr,fac,wfi,wfr 763 !arrays 764 integer :: band_procs(mband) 765 integer,allocatable :: pwind_tmp(:) 766 real(dp) :: z1(2),z2(2) 767 real(dp),allocatable :: Amat(:,:,:),Bmat(:,:,:),s1mat(:,:,:),vect1(:,:) 768 real(dp),allocatable :: vect2(:,:) 769 770 ! ************************************************************************* 771 772 mpw_tmp=max(mpw,mpw1) 773 ABI_MALLOC(vect1,(2,0:mpw_tmp)) 774 ABI_MALLOC(vect2,(2,0:mpw_tmp)) 775 ABI_MALLOC(s1mat,(2,dtefield%mband_occ,dtefield%mband_occ)) 776 ABI_MALLOC(pwind_tmp,(mpw_tmp)) 777 ABI_MALLOC(Amat,(2,dtefield%mband_occ,dtefield%mband_occ)) 778 ABI_MALLOC(Bmat,(2,dtefield%mband_occ,dtefield%mband_occ)) 779 vect1(:,0) = zero ; vect2(:,0) = zero 780 s1mat(:,:,:)=zero 781 grad_berry(:,:,:) = zero 782 783 me_band = mpi_enreg%me_band 784 call proc_distrb_band(band_procs,mpi_enreg%proc_distrb,ikpt,isppol,mband,me_band,mpi_enreg%me_kpt,mpi_enreg%comm_band) 785 786 do idir=1,3 787 fac = dtefield%efield_dot(idir)*dble(nkpt)/& 788 & (dble(dtefield%nstr(idir))*four_pi) 789 790 ! prepare 791 ikpt1 = dtefield%ikpt_dk(ikpt,1,idir) 792 icg1 = dtefield%cgindex(ikpt1,isppol+nsppol) 793 npw_k1 = npwar1(ikpt) 794 npw_k2 = npwar1(ikpt1) 795 pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,3,idir) 796 797 !TODO: this algorithm is very inefficient: the looping is in the wrong order. 798 ! should be possible to store a temp vector and make it into a BLAS call... 799 ! basically boils down to an internal sum over iband 800 ! grad_berry(ipw,jband) = fac * wf(ipw,:) * qmat(:,jband) 801 do ipw = 1, npw_k1 802 jpw = pwind_tmp(ipw) 803 804 if (jpw > 0) then 805 iband_me = 0 806 do iband = 1, dtefield%mband_occ 807 if (band_procs(iband) /= me_band) cycle 808 iband_me = iband_me + 1 809 wfr = cg1(1,icg1 + (iband_me - 1)*npw_k2*nspinor + jpw) 810 wfi = cg1(2,icg1 + (iband_me - 1)*npw_k2*nspinor + jpw) 811 812 do jband = 1, dtefield%mband_occ 813 814 grad_berry(1,ipw,jband) = & 815 & grad_berry(1,ipw,jband) + & 816 & fac*qmat(1,iband,jband,ikpt,1,idir)*wfr - fac*qmat(2,iband,jband,ikpt,1,idir)*wfi 817 818 grad_berry(2,ipw,jband) = & 819 & grad_berry(2,ipw,jband) + & 820 & fac*qmat(1,iband,jband,ikpt,1,idir)*wfi + fac*qmat(2,iband,jband,ikpt,1,idir)*wfr 821 822 end do 823 end do 824 end if 825 826 end do 827 828 ! compute <u^(0)_{k_j+n}|u^(1)_{k_j+1,q}> matrix---------------------------------------------------- 829 830 ! prepare to calculate overlap matrix 831 ikptn = dtefield%ikpt_dk(ikpt,7,idir) 832 ikpt1 = dtefield%ikpt_dk(ikpt,1,idir) 833 icg = dtefield%cgindex(ikptn,isppol) 834 icg1 = dtefield%cgindex(ikpt1,isppol+nsppol) 835 npw_k1 = npwarr(ikptn) 836 npw_k2 = npwar1(ikpt1) 837 pwind_tmp(1:npw_k1) = pwindall((ikptn-1)*mpw_tmp+1:(ikptn-1)*mpw_tmp+npw_k1,7,idir) 838 839 840 vect1(:,0) = zero ; vect2(:,0) = zero 841 do jband = 1, dtefield%mband_occ 842 vect2(:,1:npw_k2) = & 843 & cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor) 844 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 845 846 do iband = 1, dtefield%mband_occ 847 848 pwmin = (iband-1)*npw_k1*nspinor 849 pwmax = pwmin + npw_k1*nspinor 850 vect1(:,1:npw_k1) = & 851 & cg(:,icg + 1 + pwmin:icg + pwmax) 852 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 853 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 854 & vect1,vect2) 855 s1mat(1,iband,jband) = dotr 856 s1mat(2,iband,jband) = doti 857 858 end do ! iband 859 end do !jband 860 861 ! compute <u^(0)_{-k_j+1}|u^(1)_{-k_j+n},q> matrix-------------------- 862 863 ! prepare to calculate overlap matrix 864 ikptn = dtefield%ikpt_dk(ikpt,7,idir) 865 ikptnm= dtefield%ikpt_dk(ikptn,9,idir) 866 ikpt1 = dtefield%ikpt_dk(ikpt,1,idir) 867 ikpt1m= dtefield%ikpt_dk(ikpt1,9,idir) 868 869 icg = dtefield%cgindex(ikpt1m,isppol) 870 icg1 = dtefield%cgindex(ikptnm,isppol+nsppol) 871 npw_k1 = npwarr(ikpt1m) 872 npw_k2 = npwar1(ikptnm) 873 pwind_tmp(1:npw_k1) = pwindall((ikpt1m-1)*mpw_tmp+1:(ikpt1m-1)*mpw_tmp+npw_k1,7,idir) 874 875 vect1(:,0) = zero ; vect2(:,0) = zero 876 do jband = 1, dtefield%mband_occ 877 vect2(:,1:npw_k2) = & 878 & cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor) 879 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 880 881 do iband = 1, dtefield%mband_occ 882 883 pwmin = (iband-1)*npw_k1*nspinor 884 pwmax = pwmin + npw_k1*nspinor 885 vect1(:,1:npw_k1) = & 886 & cg(:,icg + 1 + pwmin:icg + pwmax) 887 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 888 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 889 & vect1,vect2) 890 891 s1mat(1,jband,iband) = s1mat(1,jband,iband) + dotr 892 s1mat(2,jband,iband) = s1mat(2,jband,iband) + doti 893 894 end do ! iband 895 end do !jband 896 897 Amat(:,:,:)=zero 898 899 ! calculate Amat 900 do iband=1, dtefield%mband_occ 901 do jband=1, dtefield%mband_occ 902 do kband=1, dtefield%mband_occ 903 Amat(1,iband,jband) = Amat(1,iband,jband) + s1mat(1,iband,kband)*& 904 & qmat(1,kband,jband,ikpt,1,idir)& 905 & - s1mat(2,iband,kband)*qmat(2,kband,jband,ikpt,1,idir) 906 Amat(2,iband,jband) = Amat(2,iband,jband) + s1mat(1,iband,kband)*& 907 & qmat(2,kband,jband,ikpt,1,idir)& 908 & + s1mat(2,iband,kband)*qmat(1,kband,jband,ikpt,1,idir) 909 end do 910 end do 911 end do 912 913 Bmat(:,:,:)=zero 914 915 ! calculate Bmat 916 ikptn = dtefield%ikpt_dk(ikpt,7,idir) 917 do iband=1, dtefield%mband_occ 918 do jband=1, dtefield%mband_occ 919 do kband=1, dtefield%mband_occ 920 Bmat(1,jband,kband) = Bmat(1,jband,kband) + Amat(1,iband,kband)*& 921 & qmat(1,jband,iband,ikptn,1,idir)& 922 & - Amat(2,iband,kband)*qmat(2,jband,iband,ikptn,1,idir) 923 Bmat(2,jband,kband) = Bmat(2,jband,kband) + Amat(1,iband,kband)*& 924 & qmat(2,jband,iband,ikptn,1,idir)& 925 & + Amat(2,iband,kband)*qmat(1,jband,iband,ikptn,1,idir) 926 end do 927 end do 928 end do 929 930 ! calc. the second term of gradient------------------------------ 931 932 ! preparation 933 934 ikptnp1 = dtefield%ikpt_dk(ikpt,3,idir) 935 icg = dtefield%cgindex(ikptnp1,isppol) 936 npw_k1 = npwar1(ikpt) 937 npw_k2 = npwarr(ikptnp1) 938 pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,5,idir) 939 940 z1(:) = zero 941 z2(:) = zero 942 943 do ipw = 1, npw_k1 944 945 jpw = pwind_tmp(ipw) 946 947 if (jpw > 0) then 948 949 iband_me = 0 950 do iband = 1, dtefield%mband_occ 951 if (band_procs(iband) /= me_band) cycle 952 iband_me = iband_me + 1 953 wfr = cg(1,icg + (iband_me - 1)*npw_k2*nspinor + jpw) 954 wfi = cg(2,icg + (iband_me - 1)*npw_k2*nspinor + jpw) 955 956 do jband=1, dtefield%mband_occ 957 958 grad_berry(1,ipw,jband) = grad_berry(1,ipw,jband) & 959 & - fac*(Bmat(1,iband,jband)*wfr - Bmat(2,iband,jband)*wfi) 960 grad_berry(2,ipw,jband) = grad_berry(2,ipw,jband) & 961 & - fac*(Bmat(1,iband,jband)*wfi + Bmat(2,iband,jband)*wfr) 962 963 end do 964 end do 965 end if 966 end do 967 968 ! Second part of gradient of Berry phase++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 969 970 vect1(:,0) = zero ; vect2(:,0) = zero 971 972 ! prepare 973 ikpt1 = dtefield%ikpt_dk(ikpt,2,idir) 974 icg1 = dtefield%cgindex(ikpt1,isppol+nsppol) 975 npw_k1 = npwar1(ikpt) 976 npw_k2 = npwar1(ikpt1) 977 pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,4,idir) 978 979 ! write(std_out,*)'dfpt_cgwf:pwind_tmp',pwind_tmp 980 ! stop 981 982 do ipw = 1, npw_k1 983 984 jpw = pwind_tmp(ipw) 985 986 if (jpw > 0) then 987 iband_me = 0 988 do iband = 1, dtefield%mband_occ 989 if (band_procs(iband) /= me_band) cycle 990 iband_me = iband_me + 1 991 wfr = cg1(1,icg1 + (iband_me - 1)*npw_k2*nspinor + jpw) 992 wfi = cg1(2,icg1 + (iband_me - 1)*npw_k2*nspinor + jpw) 993 994 do jband = 1, dtefield%mband_occ 995 996 grad_berry(1,ipw,jband) = & 997 & grad_berry(1,ipw,jband) - & 998 & fac*qmat(1,iband,jband,ikpt,2,idir)*wfr + fac*qmat(2,iband,jband,ikpt,2,idir)*wfi 999 1000 grad_berry(2,ipw,jband) = & 1001 & grad_berry(2,ipw,jband) - & 1002 & fac*qmat(1,iband,jband,ikpt,2,idir)*wfi - fac*qmat(2,iband,jband,ikpt,2,idir)*wfr 1003 1004 end do 1005 end do 1006 end if 1007 end do 1008 1009 1010 ! compute <u^(0)_{k_j+n}|u^(1)_{k_j-1,q}> matrix---------------------------------------------------- 1011 1012 ! prepare to calculate overlap matrix 1013 ikptn = dtefield%ikpt_dk(ikpt,7,idir) 1014 ikpt1 = dtefield%ikpt_dk(ikpt,2,idir) 1015 icg = dtefield%cgindex(ikptn,isppol) 1016 icg1 = dtefield%cgindex(ikpt1,isppol+nsppol) 1017 npw_k1 = npwarr(ikptn) 1018 npw_k2 = npwar1(ikpt1) 1019 pwind_tmp(1:npw_k1) =pwindall((ikptn-1)*mpw_tmp+1:(ikptn-1)*mpw_tmp+npw_k1,8,idir) 1020 1021 vect1(:,0) = zero ; vect2(:,0) = zero 1022 jband_me = 0 1023 do jband = 1, dtefield%mband_occ 1024 if (band_procs(jband) == me_band) then 1025 jband_me = jband_me + 1 1026 vect2(:,1:npw_k2) = & 1027 & cg1(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor) 1028 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 1029 end if 1030 call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr) 1031 1032 iband_me = 0 1033 do iband = 1, dtefield%mband_occ 1034 if (band_procs(iband) /= mpi_enreg%me_band) cycle 1035 iband_me = iband_me + 1 1036 pwmin = (iband-1)*npw_k1*nspinor 1037 pwmax = pwmin + npw_k1*nspinor 1038 vect1(:,1:npw_k1) = & 1039 & cg(:,icg + 1 + pwmin:icg + pwmax) 1040 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 1041 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 1042 & vect1,vect2) 1043 s1mat(1,iband,jband) = dotr 1044 s1mat(2,iband,jband) = doti 1045 1046 end do ! iband 1047 end do !jband 1048 1049 ! compute <u^(0)_{-k_j-1}|u^(1)_{-k_j+n,q}> matrix----------------------------------------------------- 1050 1051 ! prepare to calculate overlap matrix 1052 ikptn = dtefield%ikpt_dk(ikpt,7,idir) 1053 ikptnm= dtefield%ikpt_dk(ikptn,9,idir) 1054 ikpt1 = dtefield%ikpt_dk(ikpt,2,idir) 1055 ikpt1m= dtefield%ikpt_dk(ikpt1,9,idir) 1056 icg = dtefield%cgindex(ikpt1m,isppol) 1057 icg1 = dtefield%cgindex(ikptnm,isppol+nsppol) 1058 npw_k1 = npwarr(ikpt1m) 1059 npw_k2 = npwar1(ikptnm) 1060 pwind_tmp(1:npw_k1) =pwindall((ikpt1m-1)*mpw_tmp+1:(ikpt1m-1)*mpw_tmp+npw_k1,8,idir) 1061 1062 1063 1064 vect1(:,0) = zero ; vect2(:,0) = zero 1065 jband_me = 0 1066 do jband = 1, dtefield%mband_occ 1067 if (band_procs(jband) == me_band) then 1068 jband_me = jband_me + 1 1069 vect2(:,1:npw_k2) = & 1070 & cg1(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor) 1071 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 1072 end if 1073 call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr) 1074 1075 iband_me = 0 1076 do iband = 1, dtefield%mband_occ 1077 if (band_procs(iband) /= mpi_enreg%me_band) cycle 1078 iband_me = iband_me + 1 1079 pwmin = (iband_me-1)*npw_k1*nspinor 1080 pwmax = pwmin + npw_k1*nspinor 1081 vect1(:,1:npw_k1) = & 1082 & cg(:,icg + 1 + pwmin:icg + pwmax) 1083 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 1084 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 1085 & vect1,vect2) 1086 s1mat(1,jband,iband) = s1mat(1,jband,iband) + dotr 1087 s1mat(2,jband,iband) = s1mat(2,jband,iband) + doti 1088 1089 end do ! iband 1090 end do !jband 1091 ! accumulate all iband values on each proc 1092 call xmpi_sum(s1mat,mpi_enreg%comm_band,ierr) 1093 1094 Amat(:,:,:)=zero 1095 1096 ! calculate Amat: s1mat is complete so sum over all bands here 1097 ! TODO: could reduce one loop over bands 1098 ! TODO: replace this with a BLAS call: A = s1*q in complex numbers 1099 do iband=1, dtefield%mband_occ 1100 do jband=1, dtefield%mband_occ 1101 do kband=1, dtefield%mband_occ 1102 Amat(1,iband,jband) = Amat(1,iband,jband) + s1mat(1,iband,kband)*& 1103 & qmat(1,kband,jband,ikpt,2,idir)& 1104 & - s1mat(2,iband,kband)*qmat(2,kband,jband,ikpt,2,idir) 1105 Amat(2,iband,jband) = Amat(2,iband,jband) + s1mat(1,iband,kband)*& 1106 & qmat(2,kband,jband,ikpt,2,idir)& 1107 & + s1mat(2,iband,kband)*qmat(1,kband,jband,ikpt,2,idir) 1108 end do 1109 end do 1110 end do 1111 1112 Bmat(:,:,:)=zero 1113 1114 ! calculate Bmat: as above, all matrices A B q are band-complete at this stage 1115 ! TODO: could reduce one loop over bands 1116 ! TODO: replace this with a BLAS call: B = q*A in complex numbers 1117 ikptn = dtefield%ikpt_dk(ikpt,7,idir) 1118 do iband=1, dtefield%mband_occ 1119 do jband=1, dtefield%mband_occ 1120 do kband=1, dtefield%mband_occ 1121 Bmat(1,jband,kband) = Bmat(1,jband,kband) + Amat(1,iband,kband)*& 1122 & qmat(1,jband,iband,ikptn,2,idir)& 1123 & - Amat(2,iband,kband)*qmat(2,jband,iband,ikptn,2,idir) 1124 Bmat(2,jband,kband) = Bmat(2,jband,kband) + Amat(1,iband,kband)*& 1125 & qmat(2,jband,iband,ikptn,2,idir)& 1126 + Amat(2,iband,kband)*qmat(1,jband,iband,ikptn,2,idir) 1127 end do 1128 end do 1129 end do 1130 1131 ! calc. the second term of gradient------------------------------ 1132 1133 ! preparation 1134 1135 ikptnp1 = dtefield%ikpt_dk(ikpt,4,idir) 1136 icg = dtefield%cgindex(ikptnp1,isppol) 1137 npw_k1 = npwar1(ikpt) 1138 npw_k2 = npwarr(ikptnp1) 1139 pwind_tmp(1:npw_k1) =pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,6,idir) 1140 1141 z1(:) = zero 1142 z2(:) = zero 1143 do ipw = 1, npw_k1 1144 jpw = pwind_tmp(ipw) 1145 if (jpw > 0) then 1146 iband_me = 0 1147 do iband = 1, dtefield%mband_occ 1148 if (band_procs(iband) /= me_band) cycle 1149 iband_me = iband_me + 1 1150 wfr = cg(1,icg + (iband_me - 1)*npw_k2*nspinor + jpw) 1151 wfi = cg(2,icg + (iband_me - 1)*npw_k2*nspinor + jpw) 1152 1153 do jband=1, dtefield%mband_occ 1154 grad_berry(1,ipw,jband) = grad_berry(1,ipw,jband) + fac*(Bmat(1,iband,jband)*wfr - Bmat(2,iband,jband)*wfi) 1155 grad_berry(2,ipw,jband) = grad_berry(2,ipw,jband) + fac*(Bmat(1,iband,jband)*wfi + Bmat(2,iband,jband)*wfr) 1156 end do 1157 end do 1158 end if 1159 end do 1160 1161 end do !idir 1162 1163 call xmpi_sum(grad_berry,mpi_enreg%comm_band,ierr) ! sum over iband for all previous loops 1164 1165 ABI_FREE(vect1) 1166 ABI_FREE(vect2) 1167 ABI_FREE(s1mat) 1168 ABI_FREE(Amat) 1169 ABI_FREE(Bmat) 1170 ABI_FREE(pwind_tmp) 1171 1172 end subroutine dfptff_gradberry
ABINIT/dfptff_initberry [ Functions ]
NAME
dfptff_initberry
FUNCTION
Initialization of response calculations in finite electric field.
INPUTS
dtset <type(dataset_type)> = all input variables in this dataset gmet(3,3) = reciprocal space metric tensor in bohr**-2 kg(3,mpw*mkmem_rbz) = reduced (integer) coordinates of G vecs in basis sphere kg1(3,mpw1*mkmem_rbz) = reduced (integer) coordinates of G vecs for response wfs mband = maximum number of bands mkmem_rbz = maximum number of k-points in core memory mpw = maximum number of plane waves mpw1 = maximum number of plane waves for response wavefunctions nkpt = number of k points npwarr(nkpt) = number of planewaves in basis and boundary at this k point npwar1(nkpt) = number of planewaves in basis and boundary for response wfs nsppol = 1 for unpolarized, 2 for spin-polarized occ(mband*nkpt*nsppol) = occup number for each band at each k point rprimd(3,3) = dimensional primitive vectors
OUTPUT
dtefield=variables related to response Berry-phase calculation pwindall(max(mpw,mpw1)*mkmem_rbz,8,3) = array used to compute the overlap matrices pwindall(:,1,:) <- <u^(0)_i|u^(0)_i+1> pwindall(:,2,:) <- <u^(0)_i|u^(0)_i-1> pwindall(:,3,:) <- <u^(1)_i|u^(1)_i+1> pwindall(:,4,:) <- <u^(1)_i|u^(1)_i-1> pwindall(:,5,:) <- <u^(1)_i|u^(0)_i+n+1> pwindall(:,6,:) <- <u^(1)_i|u^(0)_i+n-1> pwindall(:,7,:) <- <u^(0)_i|u^(1)_i-n+1> pwindall(:,8,:) <- <u^(0)_i|u^(1)_i-n-1>
NOTES
this duplicates in part initberry for the init of the dtefield - should be made into a common constructor in m_dtefield or somethin
SOURCE
95 subroutine dfptff_initberry(dtefield,dtset,gmet,kg,kg1,mband,mkmem_rbz,mpi_enreg,& 96 & mpw,mpw1,nkpt,npwarr,npwar1,nsppol,occ,pwindall,rprimd) 97 98 !Arguments ---------------------------------------- 99 !scalars 100 integer,intent(in) :: mband,mkmem_rbz,mpw,mpw1,nkpt,nsppol 101 type(MPI_type),intent(inout) :: mpi_enreg 102 type(dataset_type),intent(in) :: dtset 103 type(efield_type),intent(inout) :: dtefield !vz_i needs efield2 104 !arrays 105 integer,intent(in) :: kg(3,mpw*mkmem_rbz),kg1(3,mpw1*mkmem_rbz),npwar1(nkpt) 106 integer,intent(in) :: npwarr(nkpt) 107 integer,intent(out) :: pwindall(max(mpw,mpw1)*mkmem_rbz,8,3) 108 real(dp),intent(in) :: gmet(3,3),occ(mband*nkpt*nsppol),rprimd(3,3) 109 110 !Local variables ---------------------------------- 111 !scalars 112 integer :: flag,iband,icg,idir,ifor,ikg,ikg1,ikpt,ikpt1,ikpt2,ikstr 113 integer :: index,ipw,isppol,istr,iunmark,jpw,nband_k,mband_occ_k,nkstr,npw_k 114 integer :: npw_k1,orig 115 real(dp) :: ecut_eff,occ_val 116 character(len=500) :: message 117 !arrays 118 integer :: dg(3) 119 integer,allocatable :: kg_tmp(:,:),kpt_mark(:),npwarr_tmp(:),npwtot(:) 120 real(dp) :: diffk(3),dk(3) 121 real(dp),allocatable :: kpt1(:,:) 122 123 ! ************************************************************************* 124 125 !----------------------------------------------------------- 126 !---------------- initilize dtefield //--------------------- 127 !----------------------------------------------------------- 128 129 !Initialization of efield_type variables 130 dtefield%efield_dot(:) = zero 131 dtefield%dkvecs(:,:) = zero 132 dtefield%maxnstr = 0 ; dtefield%maxnkstr = 0 133 dtefield%nstr(:) = 0 ; dtefield%nkstr(:) = 0 134 ABI_MALLOC(dtefield%ikpt_dk,(nkpt,9,3)) 135 ABI_MALLOC(dtefield%cgindex,(nkpt,nsppol*2)) 136 ABI_MALLOC(dtefield%kgindex,(nkpt)) 137 dtefield%ikpt_dk(:,:,:) = 0 138 dtefield%cgindex(:,:) = 0 139 dtefield%mband_occ = 0 140 ABI_MALLOC(dtefield%nband_occ,(nsppol)) 141 dtefield%nband_occ = 0 142 pwindall(:,:,:) = 0 143 144 !Compute the number of occupied bands and check that-------- 145 !it is the same for each k-point---------------------------- 146 147 occ_val = two/(dtset%nsppol*one) 148 149 index = 0 150 do isppol = 1, nsppol 151 do ikpt = 1, nkpt 152 153 mband_occ_k = 0 154 nband_k = dtset%nband(ikpt) 155 156 do iband = 1, nband_k 157 index = index + 1 158 if (abs(occ(index) - occ_val) < tol8) mband_occ_k = mband_occ_k + 1 159 end do 160 161 if (ikpt > 1) then 162 if (dtefield%nband_occ(isppol) /= mband_occ_k) then 163 message = ' The number of valence bands is not the same for every k-point for present spin' 164 ABI_ERROR(message) 165 end if 166 else 167 dtefield%mband_occ = max(dtefield%mband_occ,mband_occ_k) 168 dtefield%nband_occ(isppol) = mband_occ_k 169 end if 170 171 end do 172 end do 173 174 !DEBUG 175 !write(std_out,*)'dfptff_initberry:nkpt',nkpt 176 !ENDDEBUG 177 178 !Compute the location of zero-order wavefunction -------------- 179 icg = 0 180 do isppol = 1, nsppol 181 do ikpt = 1, nkpt 182 183 dtefield%cgindex(ikpt,isppol) = icg 184 nband_k = dtset%nband(ikpt) 185 npw_k = npwarr(ikpt) 186 icg = icg + dtset%nspinor*npw_k*& 187 & proc_distrb_nband(mpi_enreg%proc_distrb,ikpt,nband_k,1,mpi_enreg%me_band) 188 189 end do 190 end do 191 192 !Compute the location of kg -------------- 193 ikg = 0 194 do ikpt = 1, nkpt 195 196 dtefield%kgindex(ikpt) = ikg 197 npw_k = npwarr(ikpt) 198 ikg = ikg + npw_k 199 200 end do 201 202 !Compute the location of first-order wavefunction ------------------- 203 icg = 0 204 do isppol = 1, nsppol 205 do ikpt = 1, nkpt 206 207 dtefield%cgindex(ikpt,isppol+nsppol) = icg 208 nband_k = dtset%nband(ikpt) 209 npw_k = npwar1(ikpt) 210 icg = icg + dtset%nspinor*npw_k* & 211 & proc_distrb_nband(mpi_enreg%proc_distrb,ikpt,nband_k,1,mpi_enreg%me_band) 212 213 end do 214 end do 215 216 !Compute the reciprocal lattice coordinates of the electric field----- 217 218 dtefield%efield_dot(1) = dot_product(dtset%efield(:),rprimd(:,1)) 219 dtefield%efield_dot(2) = dot_product(dtset%efield(:),rprimd(:,2)) 220 dtefield%efield_dot(3) = dot_product(dtset%efield(:),rprimd(:,3)) 221 222 write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,& 223 & ' initberry: Reciprocal lattice coordinates of the electric field',ch10,& 224 & ' efield_dot(1:3) = ',dtefield%efield_dot(1:3),ch10 225 call wrtout(std_out,message,'COLL') 226 227 !find the related k points to every k point in full BZ----------------- 228 !TODO: import hash tables for k-points, to make lookup fast and check for missing 229 ! matches 230 231 !loop over three reciprocal directions 232 do idir = 1, 3 233 234 if (dtset%rfdir(idir) == 1) then 235 236 ! Compute dk(:), the vector between a k-point and its nearest 237 ! neighbour along the direction idir 238 239 dk(:) = zero 240 dk(idir) = 1000_dp ! initialize with a big number 241 do ikpt = 2, nkpt 242 diffk(:) = abs(dtset%kptns(:,ikpt) - dtset%kptns(:,1)) 243 if ((diffk(1) < dk(1)+tol8).and.(diffk(2) < dk(2)+tol8).and.& 244 & (diffk(3) < dk(3)+tol8)) dk(:) = diffk(:) 245 end do 246 dtefield%dkvecs(:,idir) = dk(:) 247 248 ! For each k point, find k_prim such that k_prim= k + dk mod(G) 249 ! where G is a vector of the reciprocal lattice 250 251 do ikpt = 1, nkpt 252 253 ! First: k + dk 254 do ikpt1 = 1, nkpt 255 diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) - dk(:)) 256 if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then 257 dtefield%ikpt_dk(ikpt,1,idir) = ikpt1 258 exit 259 end if 260 end do 261 262 ! Second: k - dk 263 do ikpt1 = 1, nkpt 264 diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) + dk(:)) 265 if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then 266 dtefield%ikpt_dk(ikpt,2,idir) = ikpt1 267 exit 268 end if 269 end do 270 271 ! new 272 ! 3rd: k + (n+1)dk 273 274 do ikpt1 = 1, nkpt 275 diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) - dk(:) - dtset%qptn(:)) 276 if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then 277 dtefield%ikpt_dk(ikpt,3,idir) = ikpt1 278 exit 279 end if 280 end do 281 282 ! 6th: k - (n+1)dk 283 284 do ikpt1 = 1, nkpt 285 diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) + dk(:) + dtset%qptn(:)) 286 if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then 287 dtefield%ikpt_dk(ikpt,6,idir) = ikpt1 288 exit 289 end if 290 end do 291 292 ! 4th: k + (n-1)dk 293 294 do ikpt1 = 1, nkpt 295 diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) + dk(:) - dtset%qptn(:)) 296 if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then 297 dtefield%ikpt_dk(ikpt,4,idir) = ikpt1 298 exit 299 end if 300 end do 301 302 ! 5th: k - (n-1)dk 303 304 do ikpt1 = 1, nkpt 305 diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) - dk(:) + dtset%qptn(:)) 306 if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then 307 dtefield%ikpt_dk(ikpt,5,idir) = ikpt1 308 exit 309 end if 310 end do 311 312 ! 7th: k+n dk 313 314 do ikpt1 = 1, nkpt 315 diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) - dtset%qptn(:)) 316 if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then 317 dtefield%ikpt_dk(ikpt,7,idir) = ikpt1 318 exit 319 end if 320 end do 321 322 ! 8th: k-n dk 323 324 do ikpt1 = 1, nkpt 325 diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) + dtset%qptn(:)) 326 if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then 327 dtefield%ikpt_dk(ikpt,8,idir) = ikpt1 328 exit 329 end if 330 end do 331 332 ! 9th: -k 333 334 do ikpt1 = 1, nkpt 335 diffk(:) = abs(dtset%kptns(:,ikpt1) + dtset%kptns(:,ikpt)) 336 if(sum(diffk(:)) < 3*tol8) then 337 dtefield%ikpt_dk(ikpt,9,idir) = ikpt1 338 exit 339 end if 340 end do 341 342 end do ! ikpt 343 344 345 346 ! Find the string length, starting from k point 1 347 ! (all strings must have the same number of points) 348 349 nkstr = 1 350 ikpt1 = 1 351 do ikpt = 1, nkpt 352 ikpt1 = dtefield%ikpt_dk(ikpt1,1,idir) 353 if (ikpt1 == 1) exit 354 nkstr = nkstr + 1 355 end do 356 357 ! Check that the string length is a divisor of nkpt 358 if(mod(nkpt,nkstr) /= 0) then 359 write(message,'(a,i0,a,i0)')' The string length = ',nkstr,', is not a divisor of nkpt =',nkpt 360 ABI_BUG(message) 361 end if 362 dtefield%nkstr(idir) = nkstr 363 dtefield%nstr(idir) = nkpt/nkstr 364 365 end if ! dtset%rfdir(idir) == 1 366 367 write(message,'(a,i1,a,i3,a,i3)')& 368 & ' dfptff_initberry: for direction ',idir,', nkstr = ',dtefield%nkstr(idir),& 369 & ', nstr = ',dtefield%nstr(idir) 370 call wrtout(std_out,message,'COLL') 371 372 end do ! close loop over idir 373 374 dtefield%maxnstr = maxval(dtefield%nstr(:)) 375 dtefield%maxnkstr = maxval(dtefield%nkstr(:)) 376 ABI_MALLOC(dtefield%idxkstr,(dtefield%maxnkstr,dtefield%maxnstr,3)) 377 dtefield%idxkstr(:,:,:) = 0 378 379 380 !Build the different strings------------------------------------------ 381 382 ABI_MALLOC(kpt_mark,(nkpt)) 383 do idir = 1, 3 384 385 if (dtset%rfdir(idir) == 1) then 386 387 iunmark = 1 388 kpt_mark(:) = 0 389 do istr = 1, dtefield%nstr(idir) 390 391 do while(kpt_mark(iunmark) /= 0) 392 iunmark = iunmark + 1 393 end do 394 dtefield%idxkstr(1,istr,idir) = iunmark 395 kpt_mark(iunmark) = 1 396 do ikstr = 2, dtefield%nkstr(idir) 397 ikpt1 = dtefield%idxkstr(ikstr-1,istr,idir) 398 ikpt2 = dtefield%ikpt_dk(ikpt1,1,idir) 399 dtefield%idxkstr(ikstr,istr,idir) = ikpt2 400 kpt_mark(ikpt2) = 1 401 end do 402 403 end do ! istr 404 405 end if ! rfdir(idir) == 1 406 407 end do ! close loop over idir 408 409 ABI_FREE(kpt_mark) 410 411 412 !Build the array pwindall that is needed to compute the different overlap matrices 413 !at k +- dk 414 415 ABI_MALLOC(kg_tmp,(3,max(mpw,mpw1)*mkmem_rbz)) 416 ABI_MALLOC(kpt1,(3,nkpt)) 417 ABI_MALLOC(npwarr_tmp,(nkpt)) 418 ABI_MALLOC(npwtot,(nkpt)) 419 ecut_eff = dtset%ecut*(dtset%dilatmx)**2 420 421 do idir = 1, 3 422 423 if (dtset%rfdir(idir) == 1) then 424 425 dk(:) = dtefield%dkvecs(:,idir) 426 427 do ifor = 1, 2 428 429 if (ifor == 2) dk(:) = -1_dp*dk(:) 430 431 ! Build the array kpt1 = kptns + dk 432 ! all k-poins of kpt1 must be in the same BZ as those of kptns 433 kpt1(:,:) = zero 434 do ikpt = 1, nkpt 435 ikpt1 = dtefield%ikpt_dk(ikpt,ifor,idir) 436 kpt1(:,ikpt) = dtset%kptns(:,ikpt1) 437 end do ! close loop over ikpt 438 439 ! Set up the basis sphere of plane waves at kpt1 440 kg_tmp(:,:) = 0 441 call kpgio(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kg_tmp,& 442 & kpt1,mkmem_rbz,dtset%nband,nkpt,'PERS',mpi_enreg,mpw,& 443 & npwarr_tmp,npwtot,dtset%nsppol) 444 445 ikg = 0 ; ikg1 = 0 446 do ikpt = 1, nkpt 447 448 nband_k = dtset%nband(ikpt) 449 ikpt1 = dtefield%ikpt_dk(ikpt,ifor,idir) 450 451 452 dg(:) = nint(dtset%kptns(:,ikpt) + dk(:) + tol10 - & 453 & dtset%kptns(:,ikpt1)) 454 455 flag = 0; orig = 1 456 if (dg(1)*dg(1) + dg(2)*dg(2) + dg(3)*dg(3) > 0) flag = 1 457 458 ikpt1 = dtefield%ikpt_dk(ikpt,ifor,idir) 459 npw_k = npwarr(ikpt) 460 npw_k1 = npwarr_tmp(ikpt) 461 462 do ipw = 1, npw_k 463 do jpw = orig, npw_k1 464 if ((kg(1,ikg + ipw) == kg_tmp(1,ikg1 + jpw) - dg(1)).and. & 465 (kg(2,ikg + ipw) == kg_tmp(2,ikg1 + jpw) - dg(2)).and. & 466 (kg(3,ikg + ipw) == kg_tmp(3,ikg1 + jpw) - dg(3))) then 467 pwindall((ikpt-1)*max(mpw,mpw1) + ipw,ifor,idir) = jpw 468 if (flag == 0) orig = jpw 469 exit 470 end if 471 end do 472 end do 473 474 ikg = ikg + npw_k 475 ikg1 = ikg1 + npw_k1 476 477 end do ! close loop over ikpt 478 479 end do ! close loop over ifor 480 481 end if ! rfdir(idir) == 1 482 483 end do ! close loop over idir 484 485 !------------------------------------------------------------------------------ 486 !<u_q|u_q> 487 !at k +- dk 488 489 do idir = 1, 3 490 491 if (dtset%rfdir(idir) == 1) then 492 493 dk(:) = dtefield%dkvecs(:,idir) 494 495 do ifor = 1, 2 496 497 if (ifor == 2) dk(:) = -1_dp*dk(:) 498 499 ! Build the array kpt1 = kptns + qptn + dk 500 ! all k-poins of kpt1 must be in the same BZ as those of kptns 501 kpt1(:,:) = zero 502 do ikpt = 1, nkpt 503 ikpt1 = dtefield%ikpt_dk(ikpt,ifor,idir) 504 kpt1(:,ikpt) = dtset%kptns(:,ikpt1)+dtset%qptn(:) 505 end do ! close loop over ikpt 506 507 ! Set UP THE BASIS SPHERE OF PLANE waves at kpt1 508 kg_tmp(:,:) = 0 509 call kpgio(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kg_tmp,& 510 & kpt1,mkmem_rbz,dtset%nband,nkpt,'PERS',mpi_enreg,mpw1,& 511 & npwarr_tmp,npwtot,dtset%nsppol) 512 513 514 ikg = 0 ; ikg1 = 0 515 do ikpt = 1, nkpt 516 517 nband_k = dtset%nband(ikpt) 518 ikpt1 = dtefield%ikpt_dk(ikpt,ifor,idir) 519 520 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,1,mpi_enreg%me)) cycle 521 522 dg(:) = nint(dtset%kptns(:,ikpt) + dk(:) + tol10 - & 523 & dtset%kptns(:,ikpt1)) 524 525 flag = 0; orig = 1 526 if (dg(1)*dg(1) + dg(2)*dg(2) + dg(3)*dg(3) > 0) flag = 1 527 528 ikpt1 = dtefield%ikpt_dk(ikpt,ifor,idir) 529 npw_k = npwar1(ikpt) 530 npw_k1 = npwarr_tmp(ikpt) 531 532 do ipw = 1, npw_k 533 do jpw = orig, npw_k1 534 if ((kg1(1,ikg + ipw) == kg_tmp(1,ikg1 + jpw) - dg(1)).and. & 535 (kg1(2,ikg + ipw) == kg_tmp(2,ikg1 + jpw) - dg(2)).and. & 536 (kg1(3,ikg + ipw) == kg_tmp(3,ikg1 + jpw) - dg(3))) then 537 pwindall((ikpt-1)*max(mpw,mpw1) +ipw,ifor+2,idir) = jpw 538 if (flag == 0) orig = jpw 539 exit 540 end if 541 end do 542 end do 543 544 ikg = ikg + npw_k 545 ikg1 = ikg1 + npw_k1 546 547 end do ! close loop over ikpt 548 end do ! close loop over ifor 549 end if ! rfdir(idir) == 1 550 end do ! close loop over idir 551 552 553 !--------------------------------------------------------------------------- 554 555 do idir = 1, 3 556 557 if (dtset%rfdir(idir) == 1) then 558 559 dk(:) = dtset%qptn(:) + dtefield%dkvecs(:,idir) 560 561 do ifor = 1, 2 562 563 if (ifor == 2) dk(:) = dtset%qptn(:) - dtefield%dkvecs(:,idir) 564 565 ! Build the array kpt1 = kptns + qptn + dk 566 ! all k-poins of kpt1 must be in the same BZ as those of kptns 567 kpt1(:,:) = zero 568 do ikpt = 1, nkpt 569 ikpt1 = dtefield%ikpt_dk(ikpt,ifor+2,idir) 570 kpt1(:,ikpt) = dtset%kptns(:,ikpt1) 571 end do ! close loop over ikpt 572 573 ! Set UP THE BASIS SPHERE OF PLANE waves at kpt1 574 kg_tmp(:,:) = 0 575 call kpgio(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kg_tmp,& 576 & kpt1,mkmem_rbz,dtset%nband,nkpt,'PERS',mpi_enreg,mpw,& 577 & npwarr_tmp,npwtot,dtset%nsppol) 578 579 ikg = 0 ; ikg1 = 0 580 do ikpt = 1, nkpt 581 582 nband_k = dtset%nband(ikpt) 583 ikpt1 = dtefield%ikpt_dk(ikpt,ifor+2,idir) 584 585 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,1,mpi_enreg%me)) cycle 586 587 dg(:) = nint(dtset%kptns(:,ikpt) + dk(:) + tol10 - & 588 & dtset%kptns(:,ikpt1)) 589 590 flag = 0; orig = 1 591 if (dg(1)*dg(1) + dg(2)*dg(2) + dg(3)*dg(3) > 0) flag = 1 592 593 ikpt1 = dtefield%ikpt_dk(ikpt,ifor+2,idir) 594 npw_k = npwar1(ikpt) 595 npw_k1 = npwarr_tmp(ikpt) 596 597 do ipw = 1, npw_k 598 do jpw = orig, npw_k1 599 if ((kg1(1,ikg + ipw) == kg_tmp(1,ikg1 + jpw) - dg(1)).and. & 600 (kg1(2,ikg + ipw) == kg_tmp(2,ikg1 + jpw) - dg(2)).and. & 601 (kg1(3,ikg + ipw) == kg_tmp(3,ikg1 + jpw) - dg(3))) then 602 pwindall((ikpt-1)*max(mpw,mpw1)+ipw,ifor+4,idir) = jpw 603 if (flag == 0) orig = jpw 604 exit 605 end if 606 end do 607 end do 608 ikg = ikg + npw_k 609 ikg1 = ikg1 + npw_k1 610 611 end do ! close loop over ikpt 612 end do ! close loop over ifor 613 end if ! rfdir(idir) == 1 614 end do ! close loop over idir=============================================================== 615 616 617 !Build the array pwind3 that is needed to compute the overlap matrices=============================== 618 619 do idir = 1, 3 620 621 if (dtset%rfdir(idir) == 1) then 622 623 dk(:) = - dtset%qptn(:) + dtefield%dkvecs(:,idir) 624 625 do ifor = 1, 2 626 627 if (ifor == 2) dk(:) = - dtset%qptn(:) - dtefield%dkvecs(:,idir) 628 629 ! Build the array kpt1 = kptns + qptn + dk 630 ! all k-poins of kpt1 must be in the same BZ as those of kptns 631 kpt1(:,:) = zero 632 do ikpt = 1, nkpt 633 ikpt1 = dtefield%ikpt_dk(ikpt,ifor+4,idir) 634 kpt1(:,ikpt) = dtset%kptns(:,ikpt1)+ dtset%qptn(:) 635 end do ! close loop over ikpt 636 637 ! Set UP THE BASIS SPHERE OF PLANE waves at kpt1 638 kg_tmp(:,:) = 0 639 call kpgio(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kg_tmp,& 640 & kpt1,mkmem_rbz,dtset%nband,nkpt,'PERS',mpi_enreg,mpw1,& 641 & npwarr_tmp,npwtot,dtset%nsppol) 642 643 ikg = 0 ; ikg1 = 0 644 do ikpt = 1, nkpt 645 646 nband_k = dtset%nband(ikpt) 647 ikpt1 = dtefield%ikpt_dk(ikpt,ifor+4,idir) 648 649 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,1,mpi_enreg%me)) cycle 650 651 dg(:) = nint(dtset%kptns(:,ikpt) + dk(:) + tol10 - & 652 & dtset%kptns(:,ikpt1)) 653 654 flag = 0; orig = 1 655 if (dg(1)*dg(1) + dg(2)*dg(2) + dg(3)*dg(3) > 0) flag = 1 656 657 ikpt1 = dtefield%ikpt_dk(ikpt,ifor+4,idir) 658 npw_k = npwarr(ikpt) 659 npw_k1 = npwarr_tmp(ikpt) 660 661 do ipw = 1, npw_k 662 do jpw = orig, npw_k1 663 if ((kg(1,ikg + ipw) == kg_tmp(1,ikg1 + jpw) - dg(1)).and. & 664 (kg(2,ikg + ipw) == kg_tmp(2,ikg1 + jpw) - dg(2)).and. & 665 (kg(3,ikg + ipw) == kg_tmp(3,ikg1 + jpw) - dg(3))) then 666 pwindall((ikpt-1)*max(mpw,mpw1) + ipw,ifor+6,idir) = jpw 667 if (flag == 0) orig = jpw 668 exit 669 end if 670 end do 671 end do 672 ikg = ikg + npw_k 673 ikg1 = ikg1 + npw_k1 674 675 end do ! close loop over ikpt 676 end do ! close loop over ifor 677 end if ! rfdir(idir) == 1 678 end do ! close loop over idir==================================================================== 679 680 ABI_FREE(kg_tmp) 681 ABI_FREE(kpt1) 682 ABI_FREE(npwarr_tmp) 683 ABI_FREE(npwtot) 684 685 end subroutine dfptff_initberry
ABINIT/m_dfpt_fef [ Modules ]
NAME
m_dfpt_fef
FUNCTION
Response calculations in finite electric field.
COPYRIGHT
Copyright (C) 2004-2024 ABINIT group (XW). This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_dfpt_fef 23 24 use defs_basis 25 use m_abicore 26 use m_errors 27 use m_efield 28 use m_dtset 29 30 use defs_abitypes, only : MPI_type 31 use m_kg, only : kpgio 32 use m_cgtools, only : overlap_g 33 use m_xmpi 34 use m_mpinfo 35 36 implicit none 37 38 private
ABINIT/qmatrix [ Functions ]
NAME
qmatrix
FUNCTION
calculation of the inverse of the overlap matrix
INPUTS
cg(2,mpw*nspinor*mband_mem*mkmem_rbz*nsppol) = planewave coefficients of wavefunctions RF wavefunctions at k,q. dtefield = variables related to response Berry-phase calculation ikpt = the index of the current k point mband = maximum number of bands mband_mem = maximum number of bands on this cpu mkmem_rbz = maximum number of k-points in core memory mpw = maximum number of plane waves mpw1 = maximum number of plane waves for response wavefunctions nkpt = number of k points npwarr(nkpt) = number of planewaves in basis and boundary at this k point npwar1(nkpt) = number of planewaves in basis and boundary for response wfs nspinor = 1 for scalar wfs, 2 for spinor wfs nsppol = 1 for unpolarized, 2 for spin-polarized pwindall(max(mpw,mpw1)*mkmem_rbz,8,3) = array used to compute the overlap matrices
OUTPUT
qmat(2,dtefield%nband_occ,dtefield%nband_occ,nkpt,2,3) = inverse of the overlap matrix
SOURCE
2819 subroutine qmatrix(cg,dtefield,qmat,mpi_enreg,mpw,mpw1,mkmem_rbz,mband,mband_mem,npwarr,nkpt,nspinor,nsppol,pwindall) 2820 2821 use m_hide_lapack, only : dzgedi, dzgefa 2822 2823 !Arguments ---------------------------------------- 2824 !scalars 2825 integer,intent(in) :: mband,mkmem_rbz,mpw,mpw1,nkpt,nspinor,nsppol 2826 integer,intent(in) :: mband_mem 2827 type(efield_type),intent(in) :: dtefield 2828 type(MPI_type),intent(in) :: mpi_enreg 2829 !arrays 2830 integer,intent(in) :: npwarr(nkpt),pwindall(max(mpw,mpw1)*mkmem_rbz,8,3) 2831 real(dp),intent(in) :: cg(2,mpw*nspinor*mband_mem*mkmem_rbz*nsppol) 2832 real(dp),intent(out) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3) 2833 2834 !Local variables ------------------------- 2835 !scalars 2836 integer :: iband,icg,icg1,idir,ifor,ikpt,ikpt2,info,jband,job 2837 integer :: iband_me, jband_me, me_band, ierr 2838 integer :: npw_k1,npw_k2,pwmax,pwmin 2839 integer :: isppol 2840 real(dp) :: doti,dotr 2841 !arrays 2842 integer,allocatable :: ipvt(:),pwind_k(:) 2843 real(dp) :: det(2,2) 2844 integer :: band_procs(mband) 2845 real(dp),allocatable :: sinv(:,:,:),smat_k(:,:,:),vect1(:,:),vect2(:,:) 2846 real(dp),allocatable :: zgwork(:,:) 2847 2848 ! ************************************************************************* 2849 2850 ABI_MALLOC(ipvt,(dtefield%mband_occ)) 2851 ABI_MALLOC(sinv,(2,dtefield%mband_occ,dtefield%mband_occ)) 2852 ABI_MALLOC(zgwork,(2,dtefield%mband_occ)) 2853 ABI_MALLOC(vect1,(2,0:mpw)) 2854 ABI_MALLOC(vect2,(2,0:mpw)) 2855 ABI_MALLOC(smat_k,(2,dtefield%mband_occ,dtefield%mband_occ)) 2856 ABI_MALLOC(pwind_k,(max(mpw,mpw1))) 2857 vect1(:,0) = zero ; vect2(:,0) = zero 2858 2859 job = 11 2860 2861 me_band = mpi_enreg%me_band 2862 2863 !loop over k points 2864 do isppol = 1, nsppol 2865 do ikpt = 1, nkpt 2866 call proc_distrb_band(band_procs,mpi_enreg%proc_distrb,ikpt,isppol,mband,me_band,mpi_enreg%me_kpt,mpi_enreg%comm_band) 2867 2868 npw_k1 = npwarr(ikpt) 2869 icg = dtefield%cgindex(ikpt,1) 2870 do idir = 1, 3 2871 do ifor = 1, 2 2872 2873 ikpt2 = dtefield%ikpt_dk(ikpt,ifor,idir) 2874 npw_k2 = npwarr(ikpt2) 2875 icg1 = dtefield%cgindex(ikpt2,1) 2876 pwind_k(1:npw_k1) = pwindall((ikpt-1)*max(mpw,mpw1)+1:(ikpt-1)*max(mpw,mpw1)+npw_k1,ifor,idir) 2877 2878 smat_k = zero 2879 jband_me = 0 2880 do jband = 1, dtefield%nband_occ(isppol) 2881 if (band_procs(jband) == me_band) then 2882 jband_me = jband_me + 1 2883 vect2(:,1:npw_k2) = cg(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor) 2884 end if 2885 call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr) 2886 2887 if (npw_k2 < mpw) vect2(:,npw_k2+1:mpw) = zero 2888 2889 iband_me = 0 2890 do iband = 1, dtefield%nband_occ(isppol) 2891 if (band_procs(iband) /= me_band) cycle 2892 iband_me = iband_me + 1 2893 2894 pwmin = (iband_me-1)*npw_k1*nspinor 2895 pwmax = pwmin + npw_k1*nspinor 2896 vect1(:,1:npw_k1) = cg(:,icg + 1 + pwmin:icg + pwmax) 2897 if (npw_k1 < mpw) vect1(:,npw_k1+1:mpw) = zero 2898 call overlap_g(doti,dotr,mpw,npw_k1,npw_k2,nspinor,pwind_k,vect1,vect2) 2899 smat_k(1,iband,jband) = dotr 2900 smat_k(2,iband,jband) = doti 2901 end do ! iband 2902 end do !jband 2903 2904 ! recompose full s1mat on each proc 2905 call xmpi_sum(smat_k,mpi_enreg%comm_band,ierr) 2906 2907 sinv(:,:,:) = smat_k(:,:,:) 2908 2909 call dzgefa(sinv,dtefield%mband_occ,dtefield%nband_occ(isppol),ipvt,info) 2910 call dzgedi(sinv,dtefield%mband_occ,dtefield%nband_occ(isppol),ipvt,det,zgwork,job) 2911 2912 qmat(:,:,:,ikpt,ifor,idir) = sinv(:,:,:) 2913 end do 2914 end do 2915 end do !end loop over k 2916 end do 2917 2918 ABI_FREE(ipvt) 2919 ABI_FREE(sinv) 2920 ABI_FREE(zgwork) 2921 ABI_FREE(vect1) 2922 ABI_FREE(vect2) 2923 ABI_FREE(smat_k) 2924 ABI_FREE(pwind_k) 2925 2926 end subroutine qmatrix