scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

fft_index_to_3d for RISM

This commit is contained in:
Satomichi Nishihara 2021-01-24 13:03:21 +09:00 committed by Minoru Otani
parent 0d7d749371
commit ff4397bffd
7 changed files with 118 additions and 471 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
Modules/corrdipole_laue.f90
View File

@ -21,6 +21,7 @@ SUBROUTINE corrdipole_laue(rismt, lextract, ierr)
!
USE cell_base, ONLY : alat
USE constants, ONLY : K_BOLTZMANN_RY
USE fft_types, ONLY : fft_index_to_3d
USE err_rism, ONLY : IERR_RISM_NULL, IERR_RISM_INCORRECT_DATA_TYPE
USE kinds, ONLY : DP
USE mp, ONLY : mp_sum
@ -182,40 +183,21 @@ CONTAINS
LOGICAL, INTENT(IN) :: modify_csr
!
INTEGER :: ir
INTEGER :: idx
INTEGER :: i1, i2, i3
INTEGER :: ii2, ii3
INTEGER :: iz
INTEGER :: isite
LOGICAL :: offrange
!
IF (rismt%nsite < 1) THEN
RETURN
END IF
!
! ... update R-space
ii2 = rismt%dfft%my_i0r2p
ii3 = rismt%dfft%my_i0r3p
!
!$omp parallel do default(shared) private(ir, idx, i1, i2, i3, iz, isite)
DO ir = 1, rismt%dfft%nr1x * rismt%dfft%my_nr2p * rismt%dfft%my_nr3p
!$omp parallel do default(shared) private(ir, i1, i2, i3, iz, isite, offrange)
DO ir = 1, rismt%dfft%nnr
!
idx = ir - 1
i3 = idx / (rismt%dfft%nr1x * rismt%dfft%my_nr2p)
idx = idx - (rismt%dfft%nr1x * rismt%dfft%my_nr2p) * i3
i3 = i3 + ii3
IF (i3 >= rismt%dfft%nr3) THEN
CYCLE
END IF
!
i2 = idx / rismt%dfft%nr1x
idx = idx - rismt%dfft%nr1x * i2
i2 = i2 + ii2
IF (i2 >= rismt%dfft%nr2) THEN
CYCLE
END IF
!
i1 = idx
IF (i1 >= rismt%dfft%nr1) THEN
CALL fft_index_to_3d(ir, rismt%dfft, i1, i2, i3, offrange)
IF (offrange) THEN
CYCLE
END IF
!

68
Modules/corrgxy0_laue.f90
View File

@ -17,10 +17,11 @@ SUBROUTINE corrgxy0_laue(rismt, lextract, ar, ag0, ierr)
! ... lextract: if .TRUE. extract Gxy=0 terms of correlations
! ... if .FALSE. sum A(r) + A(gxy=0,z)
!
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 err_rism, ONLY : IERR_RISM_NULL, IERR_RISM_INCORRECT_DATA_TYPE
USE fft_types, ONLY : fft_index_to_3d
USE kinds, ONLY : DP
USE mp, ONLY : mp_sum
USE rism, ONLY : rism_type, ITYPE_LAUERISM
!
IMPLICIT NONE
!
@ -69,12 +70,11 @@ CONTAINS
REAL(DP), INTENT(OUT) :: ag0(rismt%nrzl, 1:*)
!
INTEGER :: ir
INTEGER :: idx
INTEGER :: i1, i2, i3
INTEGER :: ii2, ii3
INTEGER :: iz
INTEGER :: iiz
INTEGER :: isite
LOGICAL :: offrange
REAL(DP), ALLOCATABLE :: bg0(:,:)
#if defined(_OPENMP)
REAL(DP), ALLOCATABLE :: bg1(:,:)
@ -84,37 +84,19 @@ CONTAINS
RETURN
END IF
!
ii2 = rismt%dfft%my_i0r2p
ii3 = rismt%dfft%my_i0r3p
!
ALLOCATE(bg0(rismt%dfft%nr3, rismt%nsite))
bg0 = 0.0_DP
!
!$omp parallel default(shared) private(ir, idx, i1, i2, i3, iz, isite, bg1)
!$omp parallel default(shared) private(ir, i1, i2, i3, iz, isite, offrange, bg1)
#if defined(_OPENMP)
ALLOCATE(bg1(rismt%dfft%nr3, rismt%nsite))
bg1 = 0.0_DP
#endif
!$omp do
DO ir = 1, rismt%dfft%nr1x * rismt%dfft%my_nr2p * rismt%dfft%my_nr3p
DO ir = 1, rismt%dfft%nnr
!
idx = ir - 1
i3 = idx / (rismt%dfft%nr1x * rismt%dfft%my_nr2p)
idx = idx - (rismt%dfft%nr1x * rismt%dfft%my_nr2p) * i3
i3 = i3 + ii3
IF (i3 >= rismt%dfft%nr3) THEN
CYCLE
END IF
!
i2 = idx / rismt%dfft%nr1x
idx = idx - rismt%dfft%nr1x * i2
i2 = i2 + ii2
IF (i2 >= rismt%dfft%nr2) THEN
CYCLE
END IF
!
i1 = idx
IF (i1 >= rismt%dfft%nr1) THEN
CALL fft_index_to_3d(ir, rismt%dfft, i1, i2, i3, offrange)
IF (offrange) THEN
CYCLE
END IF
!
@ -148,6 +130,7 @@ CONTAINS
bg0 = bg0 / DBLE(rismt%dfft%nr1 * rismt%dfft%nr2)
!
DO isite = 1, rismt%nsite
ag0(:, isite) = 0.0_DP
DO iz = rismt%lfft%izcell_start, rismt%lfft%izcell_end
iiz = iz - rismt%lfft%izcell_start + 1
ag0(iz, isite) = bg0(iiz, isite)
@ -164,39 +147,20 @@ CONTAINS
REAL(DP), INTENT(IN) :: ag0(rismt%nrzl, 1:*)
!
INTEGER :: ir
INTEGER :: idx
INTEGER :: i1, i2, i3
INTEGER :: ii2, ii3
INTEGER :: iz
INTEGER :: isite
LOGICAL :: offrange
!
IF (rismt%nsite < 1) THEN
RETURN
END IF
!
ii2 = rismt%dfft%my_i0r2p
ii3 = rismt%dfft%my_i0r3p
!
!$omp parallel do default(shared) private(ir, idx, i1, i2, i3, iz, isite)
DO ir = 1, rismt%dfft%nr1x * rismt%dfft%my_nr2p * rismt%dfft%my_nr3p
!$omp parallel do default(shared) private(ir, i1, i2, i3, iz, isite, offrange)
DO ir = 1, rismt%dfft%nnr
!
idx = ir - 1
i3 = idx / (rismt%dfft%nr1x * rismt%dfft%my_nr2p)
idx = idx - (rismt%dfft%nr1x * rismt%dfft%my_nr2p) * i3
i3 = i3 + ii3
IF (i3 >= rismt%dfft%nr3) THEN
CYCLE
END IF
!
i2 = idx / rismt%dfft%nr1x
idx = idx - rismt%dfft%nr1x * i2
i2 = i2 + ii2
IF (i2 >= rismt%dfft%nr2) THEN
CYCLE
END IF
!
i1 = idx
IF (i1 >= rismt%dfft%nr1) THEN
CALL fft_index_to_3d(ir, rismt%dfft, i1, i2, i3, offrange)
IF (offrange) THEN
CYCLE
END IF
!

42
Modules/do_3drism.f90
View File

@ -25,6 +25,7 @@ SUBROUTINE do_3drism(rismt, maxiter, rmsconv, nbox, eta, title, ierr)
USE control_flags, ONLY : iverbosity, gamma_only
USE err_rism, ONLY : IERR_RISM_NULL, IERR_RISM_INCORRECT_DATA_TYPE, IERR_RISM_NOT_CONVERGED
USE fft_interfaces, ONLY : fwfft, invfft
USE fft_types, ONLY : fft_index_to_3d
USE io_global, ONLY : stdout
USE kinds, ONLY : DP
USE mdiis, ONLY : mdiis_type, allocate_mdiis, deallocate_mdiis, update_by_mdiis, reset_mdiis
@ -412,53 +413,24 @@ CONTAINS
SUBROUTINE clean_out_of_range()
IMPLICIT NONE
INTEGER :: ir
INTEGER :: idx
INTEGER :: i1, i2, i3
INTEGER :: ii2, ii3
LOGICAL :: offrange
!
IF (rismt%nsite < 1) THEN
RETURN
END IF
!
ii2 = rismt%dfft%my_i0r2p
ii3 = rismt%dfft%my_i0r3p
!
!$omp parallel do default(shared) private(ir, idx, i1, i2, i3)
DO ir = 1, rismt%dfft%nr1x * rismt%dfft%my_nr2p * rismt%dfft%my_nr3p
!$omp parallel do default(shared) private(ir, i1, i2, i3, offrange)
DO ir = 1, rismt%dfft%nnr
!
idx = ir - 1
i3 = idx / (rismt%dfft%nr1x * rismt%dfft%my_nr2p)
idx = idx - (rismt%dfft%nr1x * rismt%dfft%my_nr2p) * i3
i3 = i3 + ii3
IF (i3 >= rismt%dfft%nr3) THEN
CALL fft_index_to_3d(ir, rismt%dfft, i1, i2, i3, offrange)
!
IF (offrange) THEN
rismt%csr(ir, :) = 0.0_DP
rismt%hr (ir, :) = 0.0_DP
rismt%gr (ir, :) = 0.0_DP
csr (ir, :) = 0.0_DP
dcsr (ir, :) = 0.0_DP
CYCLE
END IF
!
i2 = idx / rismt%dfft%nr1x
idx = idx - rismt%dfft%nr1x * i2
i2 = i2 + ii2
IF (i2 >= rismt%dfft%nr2) THEN
rismt%csr(ir, :) = 0.0_DP
rismt%hr (ir, :) = 0.0_DP
rismt%gr (ir, :) = 0.0_DP
csr (ir, :) = 0.0_DP
dcsr (ir, :) = 0.0_DP
CYCLE
END IF
!
i1 = idx
IF (i1 >= rismt%dfft%nr1) THEN
rismt%csr(ir, :) = 0.0_DP
rismt%hr (ir, :) = 0.0_DP
rismt%gr (ir, :) = 0.0_DP
csr (ir, :) = 0.0_DP
dcsr (ir, :) = 0.0_DP
CYCLE
END IF
!
END DO

178
Modules/do_lauerism.f90
View File

@ -27,6 +27,7 @@ SUBROUTINE do_lauerism(rismt, maxiter, rmsconv, nbox, eta, charge, lboth, iref,
!
USE check_stop, ONLY : check_stop_now, stopped_by_user
USE control_flags, ONLY : iverbosity
USE fft_types, ONLY : fft_index_to_3d
USE err_rism, ONLY : IERR_RISM_NULL, IERR_RISM_INCORRECT_DATA_TYPE, IERR_RISM_NOT_CONVERGED
USE io_global, ONLY : stdout
USE kinds, ONLY : DP
@ -511,54 +512,20 @@ CONTAINS
INTEGER, INTENT(OUT) :: ngrid
!
INTEGER :: ir
INTEGER :: idx
INTEGER :: i1, i2, i3
INTEGER :: ii2, ii3
INTEGER :: iz
INTEGER :: mgrid
LOGICAL :: offrange
!
! ... for R-space
mgrid = 0
!
ii2 = rismt%dfft%my_i0r2p
ii3 = rismt%dfft%my_i0r3p
!
!$omp parallel do default(shared) private(ir, idx, i1, i2, i3, iz) reduction(+:mgrid)
DO ir = 1, rismt%dfft%nr1x * rismt%dfft%my_nr2p * rismt%dfft%my_nr3p
!$omp parallel do default(shared) private(ir, i1, i2, i3, iz, offrange) reduction(+:mgrid)
DO ir = 1, rismt%dfft%nnr
!
idx = ir - 1
i3 = idx / (rismt%dfft%nr1x * rismt%dfft%my_nr2p)
idx = idx - (rismt%dfft%nr1x * rismt%dfft%my_nr2p) * i3
i3 = i3 + ii3
IF (i3 >= rismt%dfft%nr3) THEN
IF (rismt%nsite > 0) THEN
rismt%csr (ir, :) = 0.0_DP
rismt%csdr(ir, :) = 0.0_DP
rismt%hr (ir, :) = 0.0_DP
rismt%gr (ir, :) = 0.0_DP
cst (ir, :) = 0.0_DP
dcst (ir, :) = 0.0_DP
END IF
CYCLE
END IF
CALL fft_index_to_3d(ir, rismt%dfft, i1, i2, i3, offrange)
!
i2 = idx / rismt%dfft%nr1x
idx = idx - rismt%dfft%nr1x * i2
i2 = i2 + ii2
IF (i2 >= rismt%dfft%nr2) THEN
IF (rismt%nsite > 0) THEN
rismt%csr (ir, :) = 0.0_DP
rismt%csdr(ir, :) = 0.0_DP
rismt%hr (ir, :) = 0.0_DP
rismt%gr (ir, :) = 0.0_DP
cst (ir, :) = 0.0_DP
dcst (ir, :) = 0.0_DP
END IF
CYCLE
END IF
!
i1 = idx
IF (i1 >= rismt%dfft%nr1) THEN
IF (offrange) THEN
IF (rismt%nsite > 0) THEN
rismt%csr (ir, :) = 0.0_DP
rismt%csdr(ir, :) = 0.0_DP
@ -638,40 +605,22 @@ CONTAINS
!
SUBROUTINE barrier_gr()
IMPLICIT NONE
INTEGER :: ir
INTEGER :: idx
INTEGER :: i1, i2, i3
INTEGER :: ii2, ii3
INTEGER :: iz
INTEGER :: ir
INTEGER :: i1, i2, i3
INTEGER :: iz
LOGICAL :: offrange
!
IF (rismt%nsite < 1) THEN
RETURN
END IF
!
! ... for R-space
ii2 = rismt%dfft%my_i0r2p
ii3 = rismt%dfft%my_i0r3p
!
!$omp parallel do default(shared) private(ir, idx, i1, i2, i3, iz)
DO ir = 1, rismt%dfft%nr1x * rismt%dfft%my_nr2p * rismt%dfft%my_nr3p
!$omp parallel do default(shared) private(ir, i1, i2, i3, iz, offrange)
DO ir = 1, rismt%dfft%nnr
!
idx = ir - 1
i3 = idx / (rismt%dfft%nr1x * rismt%dfft%my_nr2p)
idx = idx - (rismt%dfft%nr1x * rismt%dfft%my_nr2p) * i3
i3 = i3 + ii3
IF (i3 >= rismt%dfft%nr3) THEN
CYCLE
END IF
!
i2 = idx / rismt%dfft%nr1x
idx = idx - rismt%dfft%nr1x * i2
i2 = i2 + ii2
IF (i2 >= rismt%dfft%nr2) THEN
CYCLE
END IF
!
i1 = idx
IF (i1 >= rismt%dfft%nr1) THEN
CALL fft_index_to_3d(ir, rismt%dfft, i1, i2, i3, offrange)
IF (offrange) THEN
CYCLE
END IF
!
@ -706,39 +655,20 @@ CONTAINS
!
SUBROUTINE correct_csr()
IMPLICIT NONE
INTEGER :: ir
INTEGER :: idx
INTEGER :: i1, i2, i3
INTEGER :: ii2, ii3
INTEGER :: iz
INTEGER :: ir
INTEGER :: i1, i2, i3
INTEGER :: iz
LOGICAL :: offrange
!
IF (rismt%nsite < 1) THEN
RETURN
END IF
!
ii2 = rismt%dfft%my_i0r2p
ii3 = rismt%dfft%my_i0r3p
!
!$omp parallel do default(shared) private(ir, idx, i1, i2, i3, iz)
DO ir = 1, rismt%dfft%nr1x * rismt%dfft%my_nr2p * rismt%dfft%my_nr3p
!$omp parallel do default(shared) private(ir, i1, i2, i3, iz, offrange)
DO ir = 1, rismt%dfft%nnr
!
idx = ir - 1
i3 = idx / (rismt%dfft%nr1x * rismt%dfft%my_nr2p)
idx = idx - (rismt%dfft%nr1x * rismt%dfft%my_nr2p) * i3
i3 = i3 + ii3
IF (i3 >= rismt%dfft%nr3) THEN
CYCLE
END IF
!
i2 = idx / rismt%dfft%nr1x
idx = idx - rismt%dfft%nr1x * i2
i2 = i2 + ii2
IF (i2 >= rismt%dfft%nr2) THEN
CYCLE
END IF
!
i1 = idx
IF (i1 >= rismt%dfft%nr1) THEN
CALL fft_index_to_3d(ir, rismt%dfft, i1, i2, i3, offrange)
IF (offrange) THEN
CYCLE
END IF
!
@ -762,39 +692,20 @@ CONTAINS
!
SUBROUTINE correct_hr()
IMPLICIT NONE
INTEGER :: ir
INTEGER :: idx
INTEGER :: i1, i2, i3
INTEGER :: ii2, ii3
INTEGER :: iz
INTEGER :: ir
INTEGER :: i1, i2, i3
INTEGER :: iz
LOGICAL :: offrange
!
IF (rismt%nsite < 1) THEN
RETURN
END IF
!
ii2 = rismt%dfft%my_i0r2p
ii3 = rismt%dfft%my_i0r3p
!
!$omp parallel do default(shared) private(ir, idx, i1, i2, i3, iz)
DO ir = 1, rismt%dfft%nr1x * rismt%dfft%my_nr2p * rismt%dfft%my_nr3p
!$omp parallel do default(shared) private(ir, i1, i2, i3, iz, offrange)
DO ir = 1, rismt%dfft%nnr
!
idx = ir - 1
i3 = idx / (rismt%dfft%nr1x * rismt%dfft%my_nr2p)
idx = idx - (rismt%dfft%nr1x * rismt%dfft%my_nr2p) * i3
i3 = i3 + ii3
IF (i3 >= rismt%dfft%nr3) THEN
CYCLE
END IF
!
i2 = idx / rismt%dfft%nr1x
idx = idx - rismt%dfft%nr1x * i2
i2 = i2 + ii2
IF (i2 >= rismt%dfft%nr2) THEN
CYCLE
END IF
!
i1 = idx
IF (i1 >= rismt%dfft%nr1) THEN
CALL fft_index_to_3d(ir, rismt%dfft, i1, i2, i3, offrange)
IF (offrange) THEN
CYCLE
END IF
!
@ -824,12 +735,11 @@ CONTAINS
INTEGER :: isolV
INTEGER :: natom
INTEGER :: ir
INTEGER :: idx
INTEGER :: i1, i2, i3
INTEGER :: ii2, ii3
INTEGER :: iz
INTEGER :: itsta, itend
INTEGER :: izsta, izend
LOGICAL :: offrange
REAL(DP) :: rhov1
REAL(DP) :: rhov2
REAL(DP) :: rhovt1
@ -856,9 +766,6 @@ CONTAINS
CALL errore('do_lauerism', 'rhovt2 is not positive', 1)
END IF
!
ii2 = rismt%dfft%my_i0r2p
ii3 = rismt%dfft%my_i0r3p
!
DO iq = rismt%mp_site%isite_start, rismt%mp_site%isite_end
iiq = iq - rismt%mp_site%isite_start + 1
iv = iuniq_to_isite(1, iq)
@ -872,26 +779,11 @@ CONTAINS
cst( :, iiq) = 0.0_DP
dcst(:, iiq) = 0.0_DP
!
!$omp parallel do default(shared) private(ir, idx, i1, i2, i3, iz)
DO ir = 1, rismt%dfft%nr1x * rismt%dfft%my_nr2p * rismt%dfft%my_nr3p
!$omp parallel do default(shared) private(ir, i1, i2, i3, iz, offrange)
DO ir = 1, rismt%dfft%nnr
!
idx = ir - 1
i3 = idx / (rismt%dfft%nr1x * rismt%dfft%my_nr2p)
idx = idx - (rismt%dfft%nr1x * rismt%dfft%my_nr2p) * i3
i3 = i3 + ii3
IF (i3 >= rismt%dfft%nr3) THEN
CYCLE
END IF
!
i2 = idx / rismt%dfft%nr1x
idx = idx - rismt%dfft%nr1x * i2
i2 = i2 + ii2
IF (i2 >= rismt%dfft%nr2) THEN
CYCLE
END IF
!
i1 = idx
IF (i1 >= rismt%dfft%nr1) THEN
CALL fft_index_to_3d(ir, rismt%dfft, i1, i2, i3, offrange)
IF (offrange) THEN
CYCLE
END IF
!

73
Modules/guess_3drism.f90
View File

@ -15,6 +15,7 @@ SUBROUTINE guess_3drism(rismt, ierr)
USE cell_base, ONLY : at, alat
USE constants, ONLY : eps4, K_BOLTZMANN_RY
USE err_rism, ONLY : IERR_RISM_NULL, IERR_RISM_INCORRECT_DATA_TYPE
USE fft_types, ONLY : fft_index_to_3d
USE kinds, ONLY : DP
USE mp, ONLY : mp_max
USE rism, ONLY : rism_type, ITYPE_3DRISM, ITYPE_LAUERISM
@ -33,9 +34,8 @@ SUBROUTINE guess_3drism(rismt, ierr)
INTEGER :: isolV
INTEGER :: iatom
INTEGER :: ir
INTEGER :: idx
INTEGER :: i1, i2, i3
INTEGER :: ii2, ii3
LOGICAL :: offrange
LOGICAL :: laue
REAL(DP) :: beta
REAL(DP) :: qv
@ -90,33 +90,15 @@ SUBROUTINE guess_3drism(rismt, ierr)
iatom = isite_to_iatom(iv)
qv = solVs(isolV)%charge(iatom)
!
ii2 = rismt%dfft%my_i0r2p
ii3 = rismt%dfft%my_i0r3p
!
! ... set csr to be zero
csmax = 0.0_DP
rismt%csr(:, iiq) = 0.0_DP
!
! ... set csr initially
DO ir = 1, rismt%dfft%nr1x * rismt%dfft%my_nr2p * rismt%dfft%my_nr3p
DO ir = 1, rismt%dfft%nnr
!
idx = ir - 1
i3 = idx / (rismt%dfft%nr1x * rismt%dfft%my_nr2p)
idx = idx - (rismt%dfft%nr1x * rismt%dfft%my_nr2p) * i3
i3 = i3 + ii3
IF (i3 >= rismt%dfft%nr3) THEN
CYCLE
END IF
!
i2 = idx / rismt%dfft%nr1x
idx = idx - rismt%dfft%nr1x * i2
i2 = i2 + ii2
IF (i2 >= rismt%dfft%nr2) THEN
CYCLE
END IF
!
i1 = idx
IF (i1 >= rismt%dfft%nr1) THEN
CALL fft_index_to_3d(ir, rismt%dfft, i1, i2, i3, offrange)
IF (offrange) THEN
CYCLE
END IF
!
@ -136,25 +118,10 @@ SUBROUTINE guess_3drism(rismt, ierr)
CALL mp_max(csmax, rismt%mp_site%intra_sitg_comm)
!
! ... correct csr to be smooth
DO ir = 1, rismt%dfft%nr1x * rismt%dfft%my_nr2p * rismt%dfft%my_nr3p
DO ir = 1, rismt%dfft%nnr
!
idx = ir - 1
i3 = idx / (rismt%dfft%nr1x * rismt%dfft%my_nr2p)
idx = idx - (rismt%dfft%nr1x * rismt%dfft%my_nr2p) * i3
i3 = i3 + ii3
IF (i3 >= rismt%dfft%nr3) THEN
CYCLE
END IF
!
i2 = idx / rismt%dfft%nr1x
idx = idx - rismt%dfft%nr1x * i2
i2 = i2 + ii2
IF (i2 >= rismt%dfft%nr2) THEN
CYCLE
END IF
!
i1 = idx
IF (i1 >= rismt%dfft%nr1) THEN
CALL fft_index_to_3d(ir, rismt%dfft, i1, i2, i3, offrange)
IF (offrange) THEN
CYCLE
END IF
!
@ -213,28 +180,10 @@ CONTAINS
!
z0 = 0.5_DP * at(3, 3)
!
ii2 = rismt%dfft%my_i0r2p
ii3 = rismt%dfft%my_i0r3p
!
DO ir = 1, rismt%dfft%nr1x * rismt%dfft%my_nr2p * rismt%dfft%my_nr3p
DO ir = 1, rismt%dfft%nnr
!
idx = ir - 1
i3 = idx / (rismt%dfft%nr1x * rismt%dfft%my_nr2p)
idx = idx - (rismt%dfft%nr1x * rismt%dfft%my_nr2p) * i3
i3 = i3 + ii3
IF (i3 >= rismt%dfft%nr3) THEN
CYCLE
END IF
!
i2 = idx / rismt%dfft%nr1x
idx = idx - rismt%dfft%nr1x * i2
i2 = i2 + ii2
IF (i2 >= rismt%dfft%nr2) THEN
CYCLE
END IF
!
i1 = idx
IF (i1 >= rismt%dfft%nr1) THEN
CALL fft_index_to_3d(ir, rismt%dfft, i1, i2, i3, offrange)
IF (offrange) THEN
CYCLE
END IF
!

153
Modules/lj_solute.f90
View File

@ -233,6 +233,7 @@ SUBROUTINE lj_setup_solU_vlj_x(iq, rismt, rsmax, laue)
! ... for a solvent's site
!
USE cell_base, ONLY : alat, at
USE fft_types, ONLY : fft_index_to_3d
USE kinds, ONLY : DP
USE rism, ONLY : rism_type
USE solute, ONLY : solU_nat, solU_tau, solU_ljeps, solU_ljsig, isup_to_iuni
@ -251,13 +252,11 @@ SUBROUTINE lj_setup_solU_vlj_x(iq, rismt, rsmax, laue)
INTEGER :: iv
INTEGER :: isolV
INTEGER :: iatom
INTEGER :: ir
INTEGER :: idx
INTEGER :: ir, nr
INTEGER :: i1, i2, i3
INTEGER :: ii2, ii3
INTEGER :: n1, n2, n3
INTEGER :: nx1, nx2, nx3
INTEGER :: ia, iia
LOGICAL :: offrange
REAL(DP) :: tau_r(3)
REAL(DP) :: r1, r2, r3
REAL(DP) :: r3offset
@ -271,14 +270,10 @@ SUBROUTINE lj_setup_solU_vlj_x(iq, rismt, rsmax, laue)
REAL(DP) :: vlj
!
! ... FFT box
n1 = rismt%dfft%nr1
n2 = rismt%dfft%nr2
n3 = rismt%dfft%nr3
nx1 = rismt%dfft%nr1x
nx2 = rismt%dfft%my_nr2p
nx3 = rismt%dfft%my_nr3p
ii2 = rismt%dfft%my_i0r2p
ii3 = rismt%dfft%my_i0r3p
n1 = rismt%dfft%nr1
n2 = rismt%dfft%nr2
n3 = rismt%dfft%nr3
nr = rismt%dfft%nnr
!
! ... solvent properties
iiq = iq - rismt%mp_site%isite_start + 1
@ -302,30 +297,13 @@ SUBROUTINE lj_setup_solU_vlj_x(iq, rismt, rsmax, laue)
END IF
!
! ... calculate potential on each FFT grid
!$omp parallel do default(shared) private(ir, idx, i1, i2, i3, r1, r2, r3, tau_r, vlj, &
!$omp parallel do default(shared) private(ir, i1, i2, i3, offrange, r1, r2, r3, tau_r, vlj, &
!$omp ia, iia, su, suv, rmax, rmin, xuv, yuv, zuv, ruv2, eu, euv, sr2, sr6, sr12)
DO ir = 1, nx1 * nx2 * nx3
DO ir = 1, nr
!
! ... create coordinate of a FFT grid
idx = ir - 1
i3 = idx / (nx1 * nx2)
idx = idx - (nx1 * nx2) * i3
i3 = i3 + ii3
IF (i3 >= n3) THEN
rismt%uljr(ir, iiq) = 0.0_DP
CYCLE
END IF
!
i2 = idx / nx1
idx = idx - nx1 * i2
i2 = i2 + ii2
IF (i2 >= n2) THEN
rismt%uljr(ir, iiq) = 0.0_DP
CYCLE
END IF
!
i1 = idx
IF (i1 >= n1) THEN
CALL fft_index_to_3d(ir, rismt%dfft, i1, i2, i3, offrange)
IF (offrange) THEN
rismt%uljr(ir, iiq) = 0.0_DP
CYCLE
END IF
@ -446,6 +424,7 @@ SUBROUTINE lj_setup_wall_x(iq, rismt, rsmax)
!
USE cell_base, ONLY : alat, at
USE constants, ONLY : tpi
USE fft_types, ONLY : fft_index_to_3d
USE kinds, ONLY : DP
USE rism, ONLY : rism_type
USE solute, ONLY : iwall, wall_tau, wall_rho, wall_ljeps, wall_ljsig, &
@ -464,12 +443,10 @@ SUBROUTINE lj_setup_wall_x(iq, rismt, rsmax)
INTEGER :: iv
INTEGER :: isolV
INTEGER :: iatom
INTEGER :: ir
INTEGER :: idx
INTEGER :: ir, nr
INTEGER :: i1, i2, i3
INTEGER :: ii2, ii3
INTEGER :: n1, n2, n3
INTEGER :: nx1, nx2, nx3
LOGICAL :: offrange
REAL(DP) :: tau_z
REAL(DP) :: r3
REAL(DP) :: r3offset
@ -502,11 +479,7 @@ SUBROUTINE lj_setup_wall_x(iq, rismt, rsmax)
n1 = rismt%dfft%nr1
n2 = rismt%dfft%nr2
n3 = rismt%dfft%nr3
nx1 = rismt%dfft%nr1x
nx2 = rismt%dfft%my_nr2p
nx3 = rismt%dfft%my_nr3p
ii2 = rismt%dfft%my_i0r2p
ii3 = rismt%dfft%my_i0r3p
nr = rismt%dfft%nnr
!
! ... solvent properties
iiq = iq - rismt%mp_site%isite_start + 1
@ -539,30 +512,13 @@ SUBROUTINE lj_setup_wall_x(iq, rismt, rsmax)
#endif
!
! ... calculate potential on each FFT grid
!$omp parallel do default(shared) private(ir, idx, i1, i2, i3, r3, tau_z, &
!$omp parallel do default(shared) private(ir, i1, i2, i3, offrange, r3, tau_z, &
!$omp vw, zuv, sr, sr2, sr3, sr6, sr9)
DO ir = 1, nx1 * nx2 * nx3
DO ir = 1, nr
!
! ... create coordinate of a FFT grid
idx = ir - 1
i3 = idx / (nx1 * nx2)
idx = idx - (nx1 * nx2) * i3
i3 = i3 + ii3
IF (i3 >= n3) THEN
rismt%uwr(ir, iiq) = 0.0_DP
CYCLE
END IF
!
i2 = idx / nx1
idx = idx - nx1 * i2
i2 = i2 + ii2
IF (i2 >= n2) THEN
rismt%uwr(ir, iiq) = 0.0_DP
CYCLE
END IF
!
i1 = idx
IF (i1 >= n1) THEN
CALL fft_index_to_3d(ir, rismt%dfft, i1, i2, i3, offrange)
IF (offrange) THEN
rismt%uwr(ir, iiq) = 0.0_DP
CYCLE
END IF
@ -771,6 +727,7 @@ SUBROUTINE lj_get_force_x(iq, rismt, force, rsmax, laue)
!
USE cell_base, ONLY : alat, at, omega
USE constants, ONLY : eps32
USE fft_types, ONLY : fft_index_to_3d
#if defined(_OPENMP)
USE ions_base, ONLY : nat
#endif
@ -795,14 +752,13 @@ SUBROUTINE lj_get_force_x(iq, rismt, force, rsmax, laue)
INTEGER :: nv
INTEGER :: isolV
INTEGER :: iatom
INTEGER :: ir
INTEGER :: idx
INTEGER :: ir, nr
INTEGER :: i1, i2, i3
INTEGER :: ii2, ii3
INTEGER :: n1, n2, n3
INTEGER :: nx1, nx2, nx3
INTEGER :: iz
INTEGER :: ia, iia
LOGICAL :: offrange
REAL(DP) :: fac
REAL(DP) :: weight
REAL(DP) :: rho_right
@ -826,11 +782,7 @@ SUBROUTINE lj_get_force_x(iq, rismt, force, rsmax, laue)
n1 = rismt%dfft%nr1
n2 = rismt%dfft%nr2
n3 = rismt%dfft%nr3
nx1 = rismt%dfft%nr1x
nx2 = rismt%dfft%my_nr2p
nx3 = rismt%dfft%my_nr3p
ii2 = rismt%dfft%my_i0r2p
ii3 = rismt%dfft%my_i0r3p
nr = rismt%dfft%nnr
!
weight = omega / DBLE(n1 * n2 * n3)
!
@ -859,7 +811,7 @@ SUBROUTINE lj_get_force_x(iq, rismt, force, rsmax, laue)
END IF
!
! ... calculate potential on each FFT grid
!$omp parallel default(shared) private(ir, idx, i1, i2, i3, r1, r2, r3, tau_r, rhog, &
!$omp parallel default(shared) private(ir, i1, i2, i3, offrange, r1, r2, r3, tau_r, rhog, &
!$omp ia, iia, su, suv, rmax, rmin, xuv, yuv, zuv, ruv2, eu, euv, sr2, sr6, sr12, &
!$omp fac, fromp)
#if defined(_OPENMP)
@ -867,26 +819,11 @@ SUBROUTINE lj_get_force_x(iq, rismt, force, rsmax, laue)
fromp = 0.0_DP
#endif
!$omp do
DO ir = 1, nx1 * nx2 * nx3
DO ir = 1, nr
!
! ... create coordinate of a FFT grid
idx = ir - 1
i3 = idx / (nx1 * nx2)
idx = idx - (nx1 * nx2) * i3
i3 = i3 + ii3
IF (i3 >= n3) THEN
CYCLE
END IF
!
i2 = idx / nx1
idx = idx - nx1 * i2
i2 = i2 + ii2
IF (i2 >= n2) THEN
CYCLE
END IF
!
i1 = idx
IF (i1 >= n1) THEN
CALL fft_index_to_3d(ir, rismt%dfft, i1, i2, i3, offrange)
IF (offrange) THEN
CYCLE
END IF
!
@ -1045,6 +982,7 @@ SUBROUTINE lj_get_stress_x(iq, rismt, sigma, rsmax, laue)
!
USE cell_base, ONLY : alat, at, omega
USE constants, ONLY : eps32
USE fft_types, ONLY : fft_index_to_3d
USE kinds, ONLY : DP
USE rism, ONLY : rism_type
USE solute, ONLY : solU_nat, solU_tau, solU_ljeps, solU_ljsig, isup_to_iuni
@ -1066,14 +1004,12 @@ SUBROUTINE lj_get_stress_x(iq, rismt, sigma, rsmax, laue)
INTEGER :: nv
INTEGER :: isolV
INTEGER :: iatom
INTEGER :: ir
INTEGER :: idx
INTEGER :: ir, nr
INTEGER :: i1, i2, i3
INTEGER :: ii2, ii3
INTEGER :: n1, n2, n3
INTEGER :: nx1, nx2, nx3
INTEGER :: iz
INTEGER :: ia, iia
LOGICAL :: offrange
REAL(DP) :: fac
REAL(DP) :: weight
REAL(DP) :: rho_right
@ -1097,11 +1033,7 @@ SUBROUTINE lj_get_stress_x(iq, rismt, sigma, rsmax, laue)
n1 = rismt%dfft%nr1
n2 = rismt%dfft%nr2
n3 = rismt%dfft%nr3
nx1 = rismt%dfft%nr1x
nx2 = rismt%dfft%my_nr2p
nx3 = rismt%dfft%my_nr3p
ii2 = rismt%dfft%my_i0r2p
ii3 = rismt%dfft%my_i0r3p
nr = rismt%dfft%nnr
!
weight = omega / DBLE(n1 * n2 * n3)
!
@ -1130,33 +1062,18 @@ SUBROUTINE lj_get_stress_x(iq, rismt, sigma, rsmax, laue)
END IF
!
! ... calculate potential on each FFT grid
!$omp parallel default(shared) private(ir, idx, i1, i2, i3, r1, r2, r3, tau_r, rhog, &
!$omp parallel default(shared) private(ir, i1, i2, i3, offrange, r1, r2, r3, tau_r, rhog, &
!$omp ia, iia, su, suv, rmax, rmin, xuv, yuv, zuv, ruv2, eu, euv, sr2, sr6, sr12, &
!$omp fac, sgomp)
#if defined(_OPENMP)
sgomp = 0.0_DP
#endif
!$omp do
DO ir = 1, nx1 * nx2 * nx3
DO ir = 1, nr
!
! ... create coordinate of a FFT grid
idx = ir - 1
i3 = idx / (nx1 * nx2)
idx = idx - (nx1 * nx2) * i3
i3 = i3 + ii3
IF (i3 >= n3) THEN
CYCLE
END IF
!
i2 = idx / rismt%dfft%nr1x
idx = idx - rismt%dfft%nr1x * i2
i2 = i2 + ii2
IF (i2 >= n2) THEN
CYCLE
END IF
!
i1 = idx
IF (i1 >= n1) THEN
CALL fft_index_to_3d(ir, rismt%dfft, i1, i2, i3, offrange)
IF (offrange) THEN
CYCLE
END IF
!

45
Modules/solvavg.f90
View File

@ -15,7 +15,7 @@ MODULE solvavg
!
USE constants, ONLY : BOHR_RADIUS_ANGS
USE cell_base, ONLY : at, alat
USE fft_types, ONLY : fft_type_descriptor
USE fft_types, ONLY : fft_type_descriptor, fft_index_to_3d
USE io_global, ONLY : ionode
USE kinds, ONLY : DP
USE lauefft, ONLY : lauefft_type
@ -238,15 +238,13 @@ CONTAINS
REAL(DP), INTENT(IN) :: zuv(:) !dimension(nnr)
!
LOGICAL :: laue
INTEGER :: ir
INTEGER :: idx
INTEGER :: ir, nr
INTEGER :: i1, i2, i3
INTEGER :: ii2, ii3
INTEGER :: n1, n2, n3
INTEGER :: nx1, nx2, nx3
INTEGER :: nrz
INTEGER :: irz
INTEGER :: irz0
LOGICAL :: offrange
REAL(DP) :: area_xy
REAL(DP), ALLOCATABLE :: ztmp(:)
!
@ -265,18 +263,14 @@ CONTAINS
n1 = lfft%dfft%nr1
n2 = lfft%dfft%nr2
n3 = lfft%dfft%nr3
nx1 = lfft%dfft%nr1x
nx2 = lfft%dfft%my_nr2p
nx3 = lfft%dfft%my_nr3p
nr = lfft%dfft%nnr
nrz = lfft%nrz
irz0 = lfft%izcell_start
ELSE
n1 = dfft%nr1
n2 = dfft%nr2
n3 = dfft%nr3
nx1 = dfft%nr1x
nx2 = dfft%my_nr2p
nx3 = dfft%my_nr3p
nr = dfft%nnr
nrz = dfft%nr3
irz0 = 1
END IF
@ -287,33 +281,10 @@ CONTAINS
ztmp = 0.0_DP
!
! ... calculate planar average
IF (laue) THEN
ii2 = lfft%dfft%my_i0r2p
ii3 = lfft%dfft%my_i0r3p
ELSE
ii2 = dfft%my_i0r2p
ii3 = dfft%my_i0r3p
END IF
!
DO ir = 1, nx1 * nx2 * nx3
DO ir = 1, nr
!
idx = ir - 1
i3 = idx / (nx1 * nx2)
idx = idx - (nx1 * nx2) * i3
i3 = i3 + ii3
IF (i3 >= n3) THEN
CYCLE
END IF
!
i2 = idx / nx1
idx = idx - nx1 * i2
i2 = i2 + ii2
IF (i2 >= n2) THEN
CYCLE
END IF
!
i1 = idx
IF (i1 >= n1) THEN
CALL fft_index_to_3d(ir, rismt%dfft, i1, i2, i3, offrange)
IF (offrange) THEN
CYCLE
END IF
!