mirror of https://gitlab.com/QEF/q-e.git
681 lines
22 KiB
Fortran
681 lines
22 KiB
Fortran
!
|
|
! Copyright (C) 2016 National Institute of Advanced Industrial Science and Technology (AIST)
|
|
! [ This code is written by Satomichi Nishihara. ]
|
|
!
|
|
! This file is distributed under the terms of the
|
|
! GNU General Public License. See the file `License'
|
|
! in the root directory of the present distribution,
|
|
! or http://www.gnu.org/copyleft/gpl.txt .
|
|
!
|
|
!---------------------------------------------------------------------------
|
|
SUBROUTINE eqn_lauelong(rismt, lboth, ierr)
|
|
!---------------------------------------------------------------------------
|
|
!
|
|
! ... solve long-range part of Laue-RISM equation, which is defined as
|
|
! ...
|
|
! ... /+inf
|
|
! ... hl1(gxy,z) = -beta * | dz' ul2(gxy,z') * x21(gxy,z'-z)
|
|
! ... /-inf
|
|
! ...
|
|
!
|
|
USE constants, ONLY : K_BOLTZMANN_RY, tpi
|
|
USE cell_base, ONLY : alat
|
|
USE err_rism, ONLY : IERR_RISM_NULL, IERR_RISM_INCORRECT_DATA_TYPE
|
|
USE kinds, ONLY : DP
|
|
USE mp, ONLY : mp_sum
|
|
USE rism, ONLY : rism_type, ITYPE_LAUERISM
|
|
USE solvmol, ONLY : get_nuniq_in_solVs, iuniq_to_isite, isite_to_isolV, isite_to_iatom, solVs
|
|
!
|
|
IMPLICIT NONE
|
|
!
|
|
TYPE(rism_type), INTENT(INOUT) :: rismt
|
|
LOGICAL, INTENT(IN) :: lboth ! both-hands calculation, or not
|
|
INTEGER, INTENT(OUT) :: ierr
|
|
!
|
|
INTEGER :: nq
|
|
INTEGER :: iq1
|
|
INTEGER :: iiq1
|
|
INTEGER :: igxy_sta
|
|
INTEGER :: igxy_end
|
|
INTEGER :: igxy_len
|
|
INTEGER :: izright_sta
|
|
INTEGER :: izright_end
|
|
INTEGER :: izleft_sta
|
|
INTEGER :: izleft_end
|
|
REAL(DP), ALLOCATABLE :: xz(:)
|
|
REAL(DP), ALLOCATABLE :: xgz(:)
|
|
REAL(DP), ALLOCATABLE :: xgzt(:,:)
|
|
REAL(DP), ALLOCATABLE :: yz(:)
|
|
REAL(DP), ALLOCATABLE :: ygz(:)
|
|
REAL(DP), ALLOCATABLE :: ygzt(:,:)
|
|
#if defined (__DEBUG_RISM)
|
|
!
|
|
CALL start_clock('3DRISM_eqnL')
|
|
#endif
|
|
!
|
|
! ... number of sites in solvents
|
|
nq = get_nuniq_in_solVs()
|
|
!
|
|
! ... check data type
|
|
IF (rismt%itype /= ITYPE_LAUERISM) THEN
|
|
ierr = IERR_RISM_INCORRECT_DATA_TYPE
|
|
RETURN
|
|
END IF
|
|
!
|
|
IF (rismt%mp_site%nsite < nq) THEN
|
|
ierr = IERR_RISM_INCORRECT_DATA_TYPE
|
|
RETURN
|
|
END IF
|
|
!
|
|
IF (rismt%ngxy < rismt%lfft%ngxy) THEN
|
|
ierr = IERR_RISM_INCORRECT_DATA_TYPE
|
|
RETURN
|
|
END IF
|
|
!
|
|
IF (rismt%ngs < rismt%lfft%nglxy) THEN
|
|
ierr = IERR_RISM_INCORRECT_DATA_TYPE
|
|
RETURN
|
|
END IF
|
|
!
|
|
IF (rismt%nrzl < rismt%lfft%nrz) THEN
|
|
ierr = IERR_RISM_INCORRECT_DATA_TYPE
|
|
RETURN
|
|
END IF
|
|
!
|
|
! ... allocate working memory
|
|
IF (rismt%nrzl > 0) THEN
|
|
ALLOCATE(xz(rismt%nrzl))
|
|
END IF
|
|
IF (rismt%nrzl * rismt%ngs > 0) THEN
|
|
ALLOCATE(xgz(rismt%nrzl * rismt%ngs))
|
|
END IF
|
|
IF (rismt%nrzl * rismt%ngs * rismt%nsite > 0) THEN
|
|
ALLOCATE(xgzt(rismt%nrzl * rismt%ngs, rismt%nsite))
|
|
END IF
|
|
!
|
|
IF (rismt%nrzl > 0) THEN
|
|
ALLOCATE(yz(rismt%nrzl))
|
|
END IF
|
|
IF (lboth) THEN
|
|
IF (rismt%nrzl * rismt%ngs > 0) THEN
|
|
ALLOCATE(ygz(rismt%nrzl * rismt%ngs))
|
|
END IF
|
|
IF (rismt%nrzl * rismt%ngs * rismt%nsite > 0) THEN
|
|
ALLOCATE(ygzt(rismt%nrzl * rismt%ngs, rismt%nsite))
|
|
END IF
|
|
END IF
|
|
!
|
|
! ----
|
|
! ... reduce x21 : x1 = > q2 * x21
|
|
! ----
|
|
! 2
|
|
DO iq1 = 1, nq
|
|
! ... properties of site1
|
|
IF (rismt%mp_site%isite_start <= iq1 .AND. iq1 <= rismt%mp_site%isite_end) THEN
|
|
iiq1 = iq1 - rismt%mp_site%isite_start + 1
|
|
ELSE
|
|
iiq1 = 0
|
|
END IF
|
|
!
|
|
! x21 -> x1 (the right-hand side)
|
|
CALL reduce_xgs(rismt%xgs(:, :, iq1), rismt%nrzl * rismt%ngs)
|
|
IF (iiq1 > 0) THEN
|
|
IF (rismt%nrzl * rismt%ngs > 0) THEN
|
|
xgzt(:, iiq1) = xgz(:)
|
|
END IF
|
|
END IF
|
|
!
|
|
IF (.NOT. lboth) THEN
|
|
CYCLE
|
|
END IF
|
|
!
|
|
! x21 -> x1 (the left-hand side)
|
|
CALL reduce_xgs(rismt%ygs(:, :, iq1), rismt%nrzl * rismt%ngs)
|
|
IF (iiq1 > 0) THEN
|
|
IF (rismt%nrzl * rismt%ngs > 0) THEN
|
|
ygzt(:, iiq1) = xgz(:)
|
|
END IF
|
|
END IF
|
|
!
|
|
END DO
|
|
!
|
|
! ...
|
|
! ... when Gxy /= 0
|
|
! ...
|
|
! ... set region of Gxy
|
|
igxy_sta = rismt%lfft%gxystart
|
|
igxy_end = rismt%lfft%ngxy
|
|
igxy_len = rismt%lfft%ngxy
|
|
!
|
|
! ... set integral regions as index of long Z-stick (i.e. expanded cell)
|
|
izright_sta = rismt%lfft%izright_start
|
|
izright_end = rismt%lfft%nrz
|
|
izleft_sta = 1
|
|
izleft_end = rismt%lfft%izleft_end
|
|
!
|
|
! ... integrate long-range part of Laue-RISM eq.
|
|
CALL int_lauelong()
|
|
!
|
|
! ...
|
|
! ... when Gxy = 0
|
|
! ...
|
|
! ... set region of Gxy
|
|
igxy_sta = 1
|
|
igxy_end = rismt%lfft%gxystart - 1
|
|
igxy_len = igxy_end - igxy_sta + 1
|
|
!
|
|
! ... set integral regions as index of long Z-stick (i.e. expanded cell)
|
|
izright_sta = rismt%lfft%izright_start0
|
|
izright_end = rismt%lfft%nrz
|
|
izleft_sta = 1
|
|
izleft_end = rismt%lfft%izleft_end0
|
|
!
|
|
! ... integrate long-range part of Laue-RISM eq.
|
|
CALL int_lauelong()
|
|
!
|
|
! ... deallocate working memory
|
|
IF (rismt%nrzl > 0) THEN
|
|
DEALLOCATE(xz)
|
|
END IF
|
|
IF (rismt%nrzl * rismt%ngs > 0) THEN
|
|
DEALLOCATE(xgz)
|
|
END IF
|
|
IF (rismt%nrzl * rismt%ngs * rismt%nsite > 0) THEN
|
|
DEALLOCATE(xgzt)
|
|
END IF
|
|
!
|
|
IF (rismt%nrzl > 0) THEN
|
|
DEALLOCATE(yz)
|
|
END IF
|
|
IF (lboth) THEN
|
|
IF (rismt%nrzl * rismt%ngs > 0) THEN
|
|
DEALLOCATE(ygz)
|
|
END IF
|
|
IF (rismt%nrzl * rismt%ngs * rismt%nsite > 0) THEN
|
|
DEALLOCATE(ygzt)
|
|
END IF
|
|
END IF
|
|
!
|
|
! ... normally done
|
|
ierr = IERR_RISM_NULL
|
|
#if defined (__DEBUG_RISM)
|
|
!
|
|
CALL stop_clock('3DRISM_eqnL')
|
|
#endif
|
|
!
|
|
CONTAINS
|
|
!
|
|
SUBROUTINE reduce_xgs(xgs, ndim)
|
|
IMPLICIT NONE
|
|
REAL(DP), INTENT(IN) :: xgs(1:ndim, 1:*)
|
|
INTEGER, INTENT(IN) :: ndim
|
|
!
|
|
INTEGER :: iq2
|
|
INTEGER :: iiq2
|
|
INTEGER :: iv2
|
|
INTEGER :: isolV2
|
|
INTEGER :: iatom2
|
|
REAL(DP) :: qv2
|
|
!
|
|
IF (rismt%nrzl * rismt%ngs > 0) THEN
|
|
xgz = 0.0_DP
|
|
END IF
|
|
!
|
|
DO iq2 = rismt%mp_site%isite_start, rismt%mp_site%isite_end
|
|
! ... properties of site2
|
|
iiq2 = iq2 - rismt%mp_site%isite_start + 1
|
|
iv2 = iuniq_to_isite(1, iq2)
|
|
isolV2 = isite_to_isolV(iv2)
|
|
iatom2 = isite_to_iatom(iv2)
|
|
qv2 = solVs(isolV2)%charge(iatom2)
|
|
!
|
|
! ... calculate x21 -> x1
|
|
IF (rismt%nrzl * rismt%ngs > 0) THEN
|
|
xgz(:) = xgz(:) + qv2 * xgs(1:(rismt%nrzl * rismt%ngs), iiq2)
|
|
END IF
|
|
!
|
|
END DO
|
|
!
|
|
IF (rismt%nrzl * rismt%ngs > 0) THEN
|
|
CALL mp_sum(xgz, rismt%mp_site%inter_sitg_comm)
|
|
END IF
|
|
!
|
|
END SUBROUTINE reduce_xgs
|
|
!
|
|
SUBROUTINE int_lauelong()
|
|
IMPLICIT NONE
|
|
INTEGER :: iq1
|
|
INTEGER :: iiq1
|
|
INTEGER :: igxy
|
|
INTEGER :: jgxy
|
|
INTEGER :: iglxy
|
|
INTEGER :: iglxy_old
|
|
INTEGER :: jglxy
|
|
INTEGER :: irz
|
|
INTEGER :: iz1, iz2
|
|
INTEGER :: nzint
|
|
INTEGER :: izint1
|
|
INTEGER :: izint2
|
|
INTEGER :: izdelt
|
|
INTEGER :: nzright
|
|
INTEGER :: nzleft
|
|
REAL(DP) :: beta
|
|
REAL(DP) :: z
|
|
REAL(DP) :: gxy
|
|
REAL(DP) :: rstep
|
|
COMPLEX(DP) :: zstep
|
|
COMPLEX(DP), ALLOCATABLE :: x21(:,:)
|
|
COMPLEX(DP), ALLOCATABLE :: x21out(:,:)
|
|
COMPLEX(DP), ALLOCATABLE :: vl2(:)
|
|
COMPLEX(DP), ALLOCATABLE :: hl1(:,:)
|
|
!
|
|
COMPLEX(DP), PARAMETER :: C_ZERO = CMPLX(0.0_DP, 0.0_DP, kind=DP)
|
|
COMPLEX(DP), PARAMETER :: C_ONE = CMPLX(1.0_DP, 0.0_DP, kind=DP)
|
|
!
|
|
EXTERNAL :: zgemv
|
|
!
|
|
! ... beta = 1 / (kB * T)
|
|
beta = 1.0_DP / K_BOLTZMANN_RY / rismt%temp
|
|
!
|
|
! ... set dz (in a.u.)
|
|
rstep = alat * rismt%lfft%zstep
|
|
zstep = CMPLX(rstep, 0.0_DP, kind=DP)
|
|
!
|
|
! ... count integral points along Z
|
|
nzright = MAX(izright_end - izright_sta + 1, 0)
|
|
nzleft = MAX(izleft_end - izleft_sta + 1, 0)
|
|
nzint = nzright + nzleft
|
|
!
|
|
! ... allocate working memory
|
|
IF (nzint > 0) THEN
|
|
ALLOCATE(x21(nzint, nzint))
|
|
END IF
|
|
IF (rismt%lfft%nrz * nzint > 0) THEN
|
|
ALLOCATE(x21out(rismt%lfft%nrz, nzint))
|
|
END IF
|
|
IF (rismt%lfft%nrz > 0) THEN
|
|
ALLOCATE(vl2(rismt%lfft%nrz))
|
|
END IF
|
|
IF (nzint * igxy_len > 0) THEN
|
|
ALLOCATE(hl1(nzint, igxy_len))
|
|
END IF
|
|
!
|
|
! ... calculate integral equation:
|
|
!
|
|
! /+inf
|
|
! hl1(gxy,z1) = -beta * | dz2 vl(gxy,z2) * x1(gxy,z2-z1)
|
|
! /-inf
|
|
!
|
|
DO iq1 = rismt%mp_site%isite_start, rismt%mp_site%isite_end
|
|
iiq1 = iq1 - rismt%mp_site%isite_start + 1
|
|
xgz(:) = xgzt(:, iiq1)
|
|
IF (lboth) THEN
|
|
ygz(:) = ygzt(:, iiq1)
|
|
END IF
|
|
!
|
|
! ... in case of zleft <= z <= zright
|
|
!
|
|
! /zright
|
|
! hl1(gxy,z1) = -beta * | dz2 vl(gxy,z2) * x1(gxy,z2-z1)
|
|
! /zleft
|
|
!
|
|
! ... solve Laue-RISM equation for each igxy
|
|
iglxy_old = -1
|
|
!
|
|
DO igxy = igxy_sta, igxy_end
|
|
jgxy = rismt%nrzl * (igxy - 1)
|
|
iglxy = rismt%lfft%igtonglxy(igxy)
|
|
jglxy = rismt%nrzl * (iglxy - 1)
|
|
gxy = rismt%lfft%gnxy(igxy)
|
|
!
|
|
! ... x(z2-z1)
|
|
IF (iglxy /= iglxy_old) THEN
|
|
iglxy_old = iglxy
|
|
!
|
|
IF (.NOT. lboth) THEN
|
|
! ... susceptibility matrix of one-hand calculation, which is symmetric
|
|
xz(1:rismt%nrzl) = xgz((1 + jglxy):(rismt%nrzl + jglxy))
|
|
!
|
|
!$omp parallel do default(shared) private(iz1, izint1, iz2, izint2, izdelt)
|
|
DO iz1 = izleft_sta, izleft_end
|
|
izint1 = iz1 - izleft_sta + 1
|
|
DO iz2 = izleft_sta, iz1
|
|
izint2 = iz2 - izleft_sta + 1
|
|
izdelt = iz1 - iz2 + 1
|
|
x21(izint2, izint1) = CMPLX(xz(izdelt), 0.0_DP, kind=DP)
|
|
END DO
|
|
END DO
|
|
!$omp end parallel do
|
|
!
|
|
!$omp parallel do default(shared) private(iz1, izint1, iz2, izint2, izdelt)
|
|
DO iz1 = izright_sta, izright_end
|
|
izint1 = nzleft + iz1 - izright_sta + 1
|
|
DO iz2 = izright_sta, iz1
|
|
izint2 = nzleft + iz2 - izright_sta + 1
|
|
izdelt = iz1 - iz2 + 1
|
|
x21(izint2, izint1) = CMPLX(xz(izdelt), 0.0_DP, kind=DP)
|
|
END DO
|
|
END DO
|
|
!$omp end parallel do
|
|
!
|
|
!$omp parallel do default(shared) private(iz1, izint1, iz2, izint2, izdelt)
|
|
DO iz1 = izright_sta, izright_end
|
|
izint1 = nzleft + iz1 - izright_sta + 1
|
|
DO iz2 = izleft_sta, izleft_end
|
|
izint2 = iz2 - izleft_sta + 1
|
|
izdelt = iz1 - iz2 + 1
|
|
x21(izint2, izint1) = CMPLX(xz(izdelt), 0.0_DP, kind=DP)
|
|
END DO
|
|
END DO
|
|
!$omp end parallel do
|
|
!
|
|
!$omp parallel do default(shared) private(izint1, izint2)
|
|
DO izint1 = 1, nzint
|
|
DO izint2 = 1, (izint1 - 1)
|
|
x21(izint1, izint2) = x21(izint2, izint1)
|
|
END DO
|
|
END DO
|
|
!$omp end parallel do
|
|
!
|
|
ELSE !IF (lboth) THEN
|
|
! ... susceptibility matrix of both-hands calculation
|
|
xz(1:rismt%nrzl) = xgz((1 + jglxy):(rismt%nrzl + jglxy))
|
|
yz(1:rismt%nrzl) = ygz((1 + jglxy):(rismt%nrzl + jglxy))
|
|
!
|
|
!$omp parallel do default(shared) private(iz1, izint1, iz2, izint2, izdelt)
|
|
DO iz1 = izleft_sta, izleft_end
|
|
izint1 = iz1 - izleft_sta + 1
|
|
DO iz2 = izleft_sta, izleft_end
|
|
izint2 = iz2 - izleft_sta + 1
|
|
izdelt = ABS(iz1 - iz2) + 1
|
|
x21(izint2, izint1) = CMPLX(yz(izdelt), 0.0_DP, kind=DP)
|
|
END DO
|
|
END DO
|
|
!$omp end parallel do
|
|
!
|
|
!$omp parallel do default(shared) private(iz1, izint1, iz2, izint2, izdelt)
|
|
DO iz1 = izright_sta, izright_end
|
|
izint1 = nzleft + iz1 - izright_sta + 1
|
|
DO iz2 = izright_sta, izright_end
|
|
izint2 = nzleft + iz2 - izright_sta + 1
|
|
izdelt = ABS(iz1 - iz2) + 1
|
|
x21(izint2, izint1) = CMPLX(xz(izdelt), 0.0_DP, kind=DP)
|
|
END DO
|
|
END DO
|
|
!$omp end parallel do
|
|
!
|
|
!$omp parallel do default(shared) private(iz1, izint1, iz2, izint2, izdelt)
|
|
DO iz1 = izleft_sta, izleft_end
|
|
izint1 = iz1 - izleft_sta + 1
|
|
DO iz2 = izright_sta, izright_end
|
|
izint2 = nzleft + iz2 - izright_sta + 1
|
|
izdelt = iz2 - iz1 + 1
|
|
x21(izint2, izint1) = CMPLX(yz(izdelt), 0.0_DP, kind=DP)
|
|
END DO
|
|
END DO
|
|
!$omp end parallel do
|
|
!
|
|
!$omp parallel do default(shared) private(iz1, izint1, iz2, izint2, izdelt)
|
|
DO iz1 = izright_sta, izright_end
|
|
izint1 = nzleft + iz1 - izright_sta + 1
|
|
DO iz2 = izleft_sta, izleft_end
|
|
izint2 = iz2 - izleft_sta + 1
|
|
izdelt = iz1 - iz2 + 1
|
|
x21(izint2, izint1) = CMPLX(xz(izdelt), 0.0_DP, kind=DP)
|
|
END DO
|
|
END DO
|
|
!$omp end parallel do
|
|
!
|
|
END IF
|
|
!
|
|
END IF
|
|
!
|
|
! ... vl(z2)
|
|
!$omp parallel do default(shared) private(iz2, izint2)
|
|
DO iz2 = izleft_sta, izleft_end
|
|
izint2 = iz2 - izleft_sta + 1
|
|
vl2(izint2) = rismt%vlgz(iz2 + jgxy)
|
|
END DO
|
|
!$omp end parallel do
|
|
!
|
|
!$omp parallel do default(shared) private(iz2, izint2)
|
|
DO iz2 = izright_sta, izright_end
|
|
izint2 = nzleft + iz2 - izright_sta + 1
|
|
vl2(izint2) = rismt%vlgz(iz2 + jgxy)
|
|
END DO
|
|
!$omp end parallel do
|
|
!
|
|
! ... hl(z1)
|
|
IF (nzint > 0) THEN
|
|
CALL zgemv('T', nzint, nzint, zstep, x21, nzint, vl2, 1, C_ZERO, hl1(1, igxy), 1)
|
|
END IF
|
|
!
|
|
END DO
|
|
!
|
|
! ... in case of z < zleft
|
|
!
|
|
! /zleft
|
|
! hl1(gxy,z1) = -beta * | dz2 vl(gxy,z2) * x1(gxy,z2-z1)
|
|
! /-inf
|
|
!
|
|
! ... solve Laue-RISM equation for each igxy
|
|
IF (rismt%lfft%xleft) THEN
|
|
iglxy_old = -1
|
|
!
|
|
DO igxy = igxy_sta, igxy_end
|
|
IF (.NOT. rismt%do_vleft(igxy)) THEN
|
|
CYCLE
|
|
END IF
|
|
!
|
|
jgxy = rismt%nrzl * (igxy - 1)
|
|
iglxy = rismt%lfft%igtonglxy(igxy)
|
|
jglxy = rismt%nrzl * (iglxy - 1)
|
|
gxy = rismt%lfft%gnxy(igxy)
|
|
!
|
|
! ... x(z2-z1)
|
|
IF (iglxy /= iglxy_old) THEN
|
|
iglxy_old = iglxy
|
|
x21out = C_ZERO
|
|
!
|
|
xz(1:rismt%nrzl) = xgz((1 + jglxy):(rismt%nrzl + jglxy))
|
|
IF (.NOT. lboth) THEN
|
|
yz(1:rismt%nrzl) = xgz((1 + jglxy):(rismt%nrzl + jglxy))
|
|
ELSE
|
|
yz(1:rismt%nrzl) = ygz((1 + jglxy):(rismt%nrzl + jglxy))
|
|
END IF
|
|
!
|
|
!$omp parallel do default(shared) private(iz1, izint1, iz2, izint2, izdelt)
|
|
DO iz1 = izleft_sta, izleft_end
|
|
izint1 = iz1 - izleft_sta + 1
|
|
DO iz2 = (iz1 - rismt%lfft%nrz + 1), 0
|
|
izint2 = iz2 + rismt%lfft%nrz
|
|
izdelt = iz1 - iz2 + 1
|
|
x21out(izint2, izint1) = CMPLX(yz(izdelt), 0.0_DP, kind=DP)
|
|
END DO
|
|
END DO
|
|
!$omp end parallel do
|
|
!
|
|
!$omp parallel do default(shared) private(iz1, izint1, iz2, izint2, izdelt)
|
|
DO iz1 = izright_sta, izright_end
|
|
izint1 = nzleft + iz1 - izright_sta + 1
|
|
DO iz2 = (iz1 - rismt%lfft%nrz + 1), 0
|
|
izint2 = iz2 + rismt%lfft%nrz
|
|
izdelt = iz1 - iz2 + 1
|
|
x21out(izint2, izint1) = CMPLX(xz(izdelt), 0.0_DP, kind=DP)
|
|
END DO
|
|
END DO
|
|
!$omp end parallel do
|
|
!
|
|
END IF
|
|
!
|
|
! ... vl(z2)
|
|
IF (rismt%lfft%gxystart > 1 .AND. igxy == 1) THEN
|
|
! ... Gxy = 0
|
|
!$omp parallel do default(shared) private(iz2, izint2, z)
|
|
DO iz2 = (-rismt%lfft%nrz + 1), 0
|
|
izint2 = iz2 + rismt%lfft%nrz
|
|
z = rismt%lfft%zleft + rismt%lfft%zoffset + rismt%lfft%zstep * DBLE(iz2 - 1)
|
|
vl2(izint2) = CMPLX(DBLE(rismt%vleft(igxy)) * z &
|
|
& + AIMAG(rismt%vleft(igxy)), 0.0_DP, kind=DP)
|
|
END DO
|
|
!$omp end parallel do
|
|
!
|
|
ELSE
|
|
! ... Gxy /= 0
|
|
!$omp parallel do default(shared) private(iz2, izint2, z)
|
|
DO iz2 = (-rismt%lfft%nrz + 1), 0
|
|
izint2 = iz2 + rismt%lfft%nrz
|
|
z = rismt%lfft%zleft + rismt%lfft%zoffset + rismt%lfft%zstep * DBLE(iz2 - 1)
|
|
vl2(izint2) = EXP( tpi * gxy * (z - rismt%lfft%zleft)) * rismt%vleft(igxy)
|
|
END DO
|
|
!$omp end parallel do
|
|
!
|
|
END IF
|
|
!
|
|
! ... hl(z1)
|
|
IF (nzint > 0) THEN
|
|
CALL zgemv('T', rismt%lfft%nrz, nzint, zstep, &
|
|
& x21out, rismt%lfft%nrz, vl2, 1, C_ONE, hl1(1, igxy), 1)
|
|
END IF
|
|
!
|
|
END DO
|
|
END IF
|
|
!
|
|
! ... in case of zright < z
|
|
!
|
|
! /+inf
|
|
! hl1(gxy,z1) = -beta * | dz2 vl(gxy,z2) * x1(gxy,z2-z1)
|
|
! /zright
|
|
!
|
|
! ... solve Laue-RISM equation for each igxy
|
|
IF (rismt%lfft%xright) THEN
|
|
iglxy_old = -1
|
|
!
|
|
DO igxy = igxy_sta, igxy_end
|
|
IF (.NOT. rismt%do_vright(igxy)) THEN
|
|
CYCLE
|
|
END IF
|
|
!
|
|
jgxy = rismt%nrzl * (igxy - 1)
|
|
iglxy = rismt%lfft%igtonglxy(igxy)
|
|
jglxy = rismt%nrzl * (iglxy - 1)
|
|
gxy = rismt%lfft%gnxy(igxy)
|
|
!
|
|
! ... x(z2-z1)
|
|
IF (iglxy /= iglxy_old) THEN
|
|
iglxy_old = iglxy
|
|
x21out = C_ZERO
|
|
!
|
|
xz(1:rismt%nrzl) = xgz((1 + jglxy):(rismt%nrzl + jglxy))
|
|
IF (.NOT. lboth) THEN
|
|
yz(1:rismt%nrzl) = xgz((1 + jglxy):(rismt%nrzl + jglxy))
|
|
ELSE
|
|
yz(1:rismt%nrzl) = ygz((1 + jglxy):(rismt%nrzl + jglxy))
|
|
END IF
|
|
!
|
|
!$omp parallel do default(shared) private(iz1, izint1, iz2, izint2, izdelt)
|
|
DO iz1 = izleft_sta, izleft_end
|
|
izint1 = iz1 - izleft_sta + 1
|
|
DO iz2 = (rismt%lfft%nrz + 1), (iz1 + rismt%lfft%nrz - 1)
|
|
izint2 = iz2 - rismt%lfft%nrz
|
|
izdelt = iz2 - iz1 + 1
|
|
x21out(izint2, izint1) = CMPLX(yz(izdelt), 0.0_DP, kind=DP)
|
|
END DO
|
|
END DO
|
|
!$omp end parallel do
|
|
!
|
|
!$omp parallel do default(shared) private(iz1, izint1, iz2, izint2, izdelt)
|
|
DO iz1 = izright_sta, izright_end
|
|
izint1 = nzleft + iz1 - izright_sta + 1
|
|
DO iz2 = (rismt%lfft%nrz + 1), (iz1 + rismt%lfft%nrz - 1)
|
|
izint2 = iz2 - rismt%lfft%nrz
|
|
izdelt = iz2 - iz1 + 1
|
|
x21out(izint2, izint1) = CMPLX(xz(izdelt), 0.0_DP, kind=DP)
|
|
END DO
|
|
END DO
|
|
!$omp end parallel do
|
|
!
|
|
END IF
|
|
!
|
|
! ... vl(z2)
|
|
IF (rismt%lfft%gxystart > 1 .AND. igxy == 1) THEN
|
|
! ... Gxy = 0
|
|
!$omp parallel do default(shared) private(iz2, izint2, z)
|
|
DO iz2 = (rismt%lfft%nrz + 1), (2 * rismt%lfft%nrz)
|
|
izint2 = iz2 - rismt%lfft%nrz
|
|
z = rismt%lfft%zleft + rismt%lfft%zoffset + rismt%lfft%zstep * DBLE(iz2 - 1)
|
|
vl2(izint2) = CMPLX(DBLE(rismt%vright(igxy)) * z &
|
|
& + AIMAG(rismt%vright(igxy)), 0.0_DP, kind=DP)
|
|
END DO
|
|
!$omp end parallel do
|
|
!
|
|
ELSE
|
|
! ... Gxy /= 0
|
|
!$omp parallel do default(shared) private(iz2, izint2, z)
|
|
DO iz2 = (rismt%lfft%nrz + 1), (2 * rismt%lfft%nrz)
|
|
izint2 = iz2 - rismt%lfft%nrz
|
|
z = rismt%lfft%zleft + rismt%lfft%zoffset + rismt%lfft%zstep * DBLE(iz2 - 1)
|
|
vl2(izint2) = EXP(-tpi * gxy * (z - rismt%lfft%zright)) * rismt%vright(igxy)
|
|
END DO
|
|
!$omp end parallel do
|
|
!
|
|
END IF
|
|
!
|
|
! ... hl(z1)
|
|
IF (nzint > 0) THEN
|
|
CALL zgemv('T', rismt%lfft%nrz, nzint, zstep, &
|
|
& x21out, rismt%lfft%nrz, vl2, 1, C_ONE, hl1(1, igxy), 1)
|
|
END IF
|
|
!
|
|
END DO
|
|
END IF
|
|
!
|
|
! ... copy hl1 -> hlgz
|
|
DO igxy = igxy_sta, igxy_end
|
|
jgxy = rismt%nrzl * (igxy - 1)
|
|
!
|
|
!$omp parallel do default(shared) private(iz1)
|
|
DO iz1 = 1, rismt%nrzl
|
|
rismt%hlgz(iz1 + jgxy, iiq1) = C_ZERO
|
|
END DO
|
|
!$omp end parallel do
|
|
END DO
|
|
!
|
|
DO igxy = igxy_sta, igxy_end
|
|
jgxy = rismt%nrzl * (igxy - 1)
|
|
!
|
|
!$omp parallel do default(shared) private(iz1, izint1)
|
|
DO iz1 = izleft_sta, izleft_end
|
|
izint1 = iz1 - izleft_sta + 1
|
|
rismt%hlgz(iz1 + jgxy, iiq1) = -beta * hl1(izint1, igxy)
|
|
END DO
|
|
!$omp end parallel do
|
|
!
|
|
!$omp parallel do default(shared) private(iz1, izint1)
|
|
DO iz1 = izright_sta, izright_end
|
|
izint1 = nzleft + iz1 - izright_sta + 1
|
|
rismt%hlgz(iz1 + jgxy, iiq1) = -beta * hl1(izint1, igxy)
|
|
END DO
|
|
!$omp end parallel do
|
|
END DO
|
|
!
|
|
END DO
|
|
!
|
|
! ... deallocate working memory
|
|
IF (nzint > 0) THEN
|
|
DEALLOCATE(x21)
|
|
END IF
|
|
IF (rismt%lfft%nrz * nzint > 0) THEN
|
|
DEALLOCATE(x21out)
|
|
END IF
|
|
IF (rismt%lfft%nrz > 0) THEN
|
|
DEALLOCATE(vl2)
|
|
END IF
|
|
IF (nzint * igxy_len > 0) THEN
|
|
DEALLOCATE(hl1)
|
|
END IF
|
|
!
|
|
END SUBROUTINE int_lauelong
|
|
!
|
|
END SUBROUTINE eqn_lauelong
|