Line data    Source code

       1             : module m_gamma_double_gpt_loop
       2             :    use m_types 
       3             :    use m_juDFT
       4             :    use m_constants
       5             :    use m_glob_tofrom_loc
       6             : contains
       7          12 :    subroutine gamma_double_gpt_loop(fi, fmpi, hybdat, mpdata, sphbesmoment, gmat, ngptm1, pgptm1, pqnrm, coul) 
       8             :       implicit none
       9             :       type(t_fleurinput), intent(in)    :: fi
      10             :       TYPE(t_mpi), INTENT(IN)           :: fmpi
      11             :       type(t_hybdat), intent(in)        :: hybdat
      12             :       TYPE(t_mpdata), intent(in)        :: mpdata
      13             :       real, intent(in)                  :: sphbesmoment(0:, :, :), gmat(:,:)
      14             :       integer, intent(in)               :: ngptm1(:), pgptm1(:,:), pqnrm(:,:)
      15             :       complex, intent(inout)            :: coul(:,:)
      16             : 
      17             :       real        :: rdum
      18             :       integer     :: igpt0, igpt2, igptp2, iqnrm1, iqnrm2, igpt1, igptp1
      19             :       integer     :: ix, iy, ix_loc, iatm1, iatm2, itype1, itype2, pe_ix
      20             :       REAL        :: qnorm, rdum1
      21             :       complex     :: cdum
      22             :       
      23          12 :       real, allocatable :: qs(:,:), qnorms(:)
      24             : 
      25          12 :       call timestart("double gpt loop")
      26             : 
      27          12 :       call setup_q_and_qnorm(fi, mpdata, qs, qnorms)
      28             : 
      29          12 :       rdum = (fpi_const)**(1.5)/fi%cell%vol**2*gmat(1, 1)
      30             : 
      31             :       !$OMP PARALLEL DO default(none) schedule(dynamic) &
      32             :       !$OMP shared(fmpi, pgptm1, ngptm1, hybdat, fi, pqnrm, mpdata) &
      33             :       !$OMP shared(coul, rdum, sphbesmoment, qs, qnorms)&
      34             :       !$OMP private(igpt1, igpt2, iy, iqnrm1, igptp1, rdum1, iatm1, iqnrm2) &
      35          12 :       !$OMP private(itype1, iatm2, itype2, cdum, ix, pe_ix, ix_loc, igptp2)
      36             :       DO igpt0 = 1, ngptm1(1)
      37             :          igpt2 = pgptm1(igpt0, 1)
      38             :          if(igpt2 /= 1) then
      39             :             ix = hybdat%nbasp + igpt2
      40             :             pe_ix = mod((ix-1), fmpi%n_size)
      41             :             ix_loc = ((ix-1)/fmpi%n_size) +1
      42             :             if(pe_ix == fmpi%n_rank) then
      43             :                iqnrm2 = pqnrm(igpt2, 1)
      44             :                igptp2 = mpdata%gptm_ptr(igpt2, 1)
      45             :                DO igpt1 = 2, igpt2
      46             :                   iy = hybdat%nbasp + igpt1
      47             :                   iqnrm1 = pqnrm(igpt1, 1)
      48             :                   igptp1 = mpdata%gptm_ptr(igpt1, 1)
      49             :                   rdum1 = dot_PRODUCT(qs(:,igptp1), qs(:,igptp2))/(qnorms(igptp1)*qnorms(igptp2))
      50             :                   do iatm1 = 1,fi%atoms%nat
      51             :                      itype1 = fi%atoms%itype(iatm1)
      52             :                      do iatm2 = 1,fi%atoms%nat 
      53             :                         itype2 = fi%atoms%itype(iatm2)
      54             :                         cdum = EXP(CMPLX(0.0, 1.0)*tpi_const* &
      55             :                                  (-dot_PRODUCT(mpdata%g(:, igptp1), fi%atoms%taual(:, iatm1)) &
      56             :                                     + dot_PRODUCT(mpdata%g(:, igptp2), fi%atoms%taual(:, iatm2))))
      57             :                         coul(iy, ix_loc) = coul(iy, ix_loc) + rdum*cdum*( &
      58             :                                           -sphbesmoment(1, itype1, iqnrm1) &
      59             :                                           *sphbesmoment(1, itype2, iqnrm2)*rdum1/3 &
      60             :                                           - sphbesmoment(0, itype1, iqnrm1) &
      61             :                                           *sphbesmoment(2, itype2, iqnrm2)/6 &
      62             :                                           - sphbesmoment(2, itype1, iqnrm1) &
      63             :                                           *sphbesmoment(0, itype2, iqnrm2)/6 &
      64             :                                           + sphbesmoment(0, itype1, iqnrm1) &
      65             :                                           *sphbesmoment(1, itype2, iqnrm2)/qnorms(igptp2)/2 &
      66             :                                           + sphbesmoment(1, itype1, iqnrm1) &
      67             :                                           *sphbesmoment(0, itype2, iqnrm2)/qnorms(igptp1)/2)
      68             :                      END DO
      69             :                   END DO
      70             :                END DO
      71             :             endif !pe_ix
      72             :          endif
      73             :       END DO
      74             :       !$OMP end parallel do
      75          12 :       call timestop("double gpt loop")
      76          12 :    end subroutine gamma_double_gpt_loop 
      77             : 
      78          12 :    subroutine setup_q_and_qnorm(fi, mpdata, qs, qnorms)
      79             :       implicit none 
      80             :       type(t_fleurinput), intent(in)    :: fi
      81             :       TYPE(t_mpdata), intent(in)        :: mpdata
      82             :       real, intent(inout), allocatable  :: qs(:,:), qnorms(:)
      83             : 
      84          12 :       real, allocatable :: g(:,:), tmp(:,:)
      85             :       integer :: i, num_gs
      86             : 
      87          12 :       num_gs = size(mpdata%g,2)
      88        9508 :       allocate(qs(3,num_gs), source=0.0)
      89             : 
      90        9484 :       g = real(mpdata%g)
      91          12 :       call dgemm("T", "N", 3,num_gs,3, 1.0,   fi%cell%bmat, 3,   g, 3, 0.0,  qs, 3)
      92          12 :       qnorms = norm2(qs, dim=1)
      93          12 :    end subroutine setup_q_and_qnorm
      94             : end module m_gamma_double_gpt_loop