quantum-espresso/Modules/guess_3drism.f90

212 lines
5.3 KiB
Fortran

!
! Copyright (C) 2015-2016 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 guess_3drism(rismt, ierr)
!---------------------------------------------------------------------------
!
! ... create initial guess of 3D-RISM or Laue-RISM.
!
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
USE solvmol, ONLY : get_nuniq_in_solVs, solVs, &
& iuniq_to_isite, isite_to_isolV, isite_to_iatom
!
IMPLICIT NONE
!
TYPE(rism_type), INTENT(INOUT) :: rismt
INTEGER, INTENT(OUT) :: ierr
!
INTEGER :: nq
INTEGER :: iq
INTEGER :: iiq
INTEGER :: iv
INTEGER :: isolV
INTEGER :: iatom
INTEGER :: ir
INTEGER :: i1, i2, i3
LOGICAL :: offrange
LOGICAL :: laue
REAL(DP) :: beta
REAL(DP) :: qv
REAL(DP) :: vlj
REAL(DP) :: cs0
REAL(DP) :: csmax
REAL(DP) :: erf0
!
REAL(DP), PARAMETER :: VLJ_MAX = eps4 ! Ry
REAL(DP), PARAMETER :: CS_SCALE = 0.1_DP
!
! ... number of sites in solvents
nq = get_nuniq_in_solVs()
!
! ... check data type
IF (rismt%itype /= ITYPE_3DRISM .AND. 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%nr < rismt%dfft%nnr) THEN
ierr = IERR_RISM_INCORRECT_DATA_TYPE
RETURN
END IF
!
! ... if no data, return as normally done
IF (rismt%nsite < 1) THEN
GOTO 1
END IF
!
! ... Laue-RISM or not
laue = .FALSE.
IF (rismt%itype == ITYPE_LAUERISM) THEN
laue = .TRUE.
END IF
!
! ... beta = 1 / (kB * T)
beta = 1.0_DP / K_BOLTZMANN_RY / rismt%temp
!
! ... create guess for each solvent's site
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)
isolV = isite_to_isolV(iv)
iatom = isite_to_iatom(iv)
qv = solVs(isolV)%charge(iatom)
!
! ... 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
!
CALL fft_index_to_3d(ir, rismt%dfft, i1, i2, i3, offrange)
IF (offrange) THEN
CYCLE
END IF
!
vlj = rismt%uljr(ir, iiq)
IF (laue) THEN ! add LJ-wall
vlj = vlj + rismt%uwr(ir, iiq)
END IF
!
IF (vlj >= VLJ_MAX) THEN
cs0 = beta * qv * rismt%vlr(ir)
csmax = MAX(csmax, ABS(cs0))
rismt%csr(ir, iiq) = cs0
END IF
!
END DO
!
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
!
CALL fft_index_to_3d(ir, rismt%dfft, i1, i2, i3, offrange)
IF (offrange) THEN
CYCLE
END IF
!
IF (csmax > 0.0_DP) THEN
cs0 = rismt%csr(ir, iiq)
erf0 = erf(ABS(cs0) / (CS_SCALE * csmax))
rismt%csr(ir, iiq) = cs0 * erf0 * erf0
END IF
!
END DO
!
END DO
!
! ... correction for Laue-RISM
IF (laue) THEN
CALL correct_edge()
END IF
!
! ... set dipole part for Laue-RISM
IF (laue) THEN
IF (rismt%nsite > 0) THEN
rismt%cda = 0.0_DP
END IF
END IF
!
! ... set Gxy=0 term for Laue-RISM
IF (laue) THEN
IF (rismt%nrzl * rismt%nsite > 0) THEN
rismt%csg0 = 0.0_DP
END IF
!
CALL corrgxy0_laue(rismt, .TRUE., rismt%csr, rismt%csg0, ierr)
IF (ierr /= IERR_RISM_NULL) THEN
RETURN
END IF
END IF
!
! ... normally done
1 CONTINUE
ierr = IERR_RISM_NULL
!
CONTAINS
!
SUBROUTINE correct_edge()
IMPLICIT NONE
!
INTEGER :: iz
REAL(DP) :: z
REAL(DP) :: z0
!
REAL(DP), PARAMETER :: Z_SCALE = 5.0_DP ! bohr
!
IF (rismt%nsite < 1) THEN
RETURN
END IF
!
z0 = 0.5_DP * at(3, 3)
!
DO ir = 1, rismt%dfft%nr1x * rismt%dfft%my_nr2p * rismt%dfft%my_nr3p
!
CALL fft_index_to_3d(ir, rismt%dfft, i1, i2, i3, offrange)
IF (offrange) THEN
CYCLE
END IF
!
IF (i3 < (rismt%dfft%nr3 - (rismt%dfft%nr3 / 2))) THEN
iz = i3 + (rismt%dfft%nr3 / 2)
ELSE
iz = i3 - rismt%dfft%nr3 + (rismt%dfft%nr3 / 2)
END IF
iz = iz + rismt%lfft%izcell_start
!
z = rismt%lfft%zleft + rismt%lfft%zoffset + rismt%lfft%zstep * DBLE(iz - 1)
!
IF (rismt%lfft%xright) THEN
erf0 = erf(alat * (z0 - z) / Z_SCALE)
rismt%csr(ir, :) = rismt%csr(ir, :) * (erf0 * erf0)
END IF
!
IF (rismt%lfft%xleft) THEN
erf0 = erf(alat * (z + z0) / Z_SCALE)
rismt%csr(ir, :) = rismt%csr(ir, :) * (erf0 * erf0)
END IF
!
END DO
!
END SUBROUTINE correct_edge
!
END SUBROUTINE guess_3drism