Skip to content

Commit

Permalink
Unify initialisation of GAPW densities in both ground-state and time-…
Browse files Browse the repository at this point in the history 
…dependent DFT
  • Loading branch information
@schulkov
schulkov committed Jan 23, 2020
1 parent 77949a1 commit 5068987
Show file tree
Hide file tree
Showing 9 changed files with 67 additions and 121 deletions.
28 changes: 8 additions & 20 deletions src/pw_env_methods.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -163,7 +163,7 @@ SUBROUTINE pw_env_rebuild(pw_env, qs_env, external_para_env)
set_vdw_pool, should_output, &
smooth_required, spherical, uf_grid, &
use_ref_cell
REAL(KIND=dp) :: cutilev, max_rpgf0_s, rel_cutoff, zet0
REAL(KIND=dp) :: cutilev, rel_cutoff
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: radius
REAL(KIND=dp), DIMENSION(:), POINTER :: cutoff
TYPE(cell_type), POINTER :: cell, cell_ref, my_cell
Expand All @@ -180,7 +180,6 @@ SUBROUTINE pw_env_rebuild(pw_env, qs_env, external_para_env)
POINTER :: rs_descs
TYPE(realspace_grid_input_type) :: input_settings
TYPE(realspace_grid_p_type), DIMENSION(:), POINTER :: rs_grids
TYPE(rho0_mpole_type), POINTER :: rho0_mpole
TYPE(section_vals_type), POINTER :: efg_section, input, linres_section, &
poisson_section, print_section, &
rs_grid_section, xc_section
Expand All @@ -199,7 +198,7 @@ SUBROUTINE pw_env_rebuild(pw_env, qs_env, external_para_env)
!
!
NULLIFY (cutoff, cell, pw_grid, old_pw_grid, dft_control, qs_kind_set, &
pw_pools, rho0_mpole, rs_descs, para_env, cell_ref, vdw_ref_grid, &
pw_pools, rs_descs, para_env, cell_ref, vdw_ref_grid, &
mt_super_ref_grid, input, poisson_section, xc_super_ref_grid, &
dct_pw_grid, vdw_grid, super_ref_grid, my_cell, rs_grids, dispersion_env)

Expand Down Expand Up @@ -235,12 +234,6 @@ SUBROUTINE pw_env_rebuild(pw_env, qs_env, external_para_env)
ELSE
ALLOCATE (pw_env%gridlevel_info)
END IF
IF (dft_control%qs_control%gapw) THEN
CALL get_qs_env(qs_env=qs_env, rho0_mpole=rho0_mpole)
CPASSERT(ASSOCIATED(rho0_mpole))
CALL get_rho0_mpole(rho0_mpole=rho0_mpole, &
zet0_h=zet0, max_rpgf0_s=max_rpgf0_s)
END IF

IF (ASSOCIATED(pw_env%cube_info)) THEN
DO igrid_level = 1, SIZE(pw_env%cube_info)
Expand Down Expand Up @@ -336,7 +329,6 @@ SUBROUTINE pw_env_rebuild(pw_env, qs_env, external_para_env)
! methods that require smoothing or nearest neighbor have to use a plane distributed setup
! find the xc properties (FIXME this could miss other xc sections that operate on the grid ...)
CALL get_qs_env(qs_env=qs_env, input=input)
xc_section => section_vals_get_subs_vals(input, "DFT%XC")
xc_deriv_method_id = section_get_ival(xc_section, "XC_GRID%XC_DERIV")
xc_smooth_method_id = section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO")
Expand Down Expand Up @@ -746,6 +738,7 @@ END SUBROUTINE pw_env_rebuild
!> \param qs_env ...
!> \par History
!> 10.2010 refactored [Joost VandeVondele]
!> 01.2020 igrid_zet0_s initialisation code is reused in rho0_s_grid_create() [Sergey Chulkov]
!> \author Joost VandeVondele
! **************************************************************************************************
SUBROUTINE compute_max_radius(radius, pw_env, qs_env)
Expand All @@ -766,7 +759,7 @@ SUBROUTINE compute_max_radius(radius, pw_env, qs_env)
INTEGER, DIMENSION(:), POINTER :: npgfa, npgfb, nshella, nshellb
INTEGER, DIMENSION(:, :), POINTER :: lshella, lshellb
REAL(KIND=dp) :: alpha, core_charge, eps_gvg, eps_rho, &
max_rpgf0_s, maxradius, zet0, zetp
max_rpgf0_s, maxradius, zet0_h, zetp
REAL(KIND=dp), DIMENSION(:, :), POINTER :: zeta, zetb
TYPE(dft_control_type), POINTER :: dft_control
TYPE(gto_basis_set_type), POINTER :: orb_basis_set
Expand All @@ -790,20 +783,15 @@ SUBROUTINE compute_max_radius(radius, pw_env, qs_env)
IF (dft_control%qs_control%gapw) THEN
CALL get_qs_env(qs_env=qs_env, rho0_mpole=rho0_mpole)
CPASSERT(ASSOCIATED(rho0_mpole))
CALL get_rho0_mpole(rho0_mpole=rho0_mpole, zet0_h=zet0, max_rpgf0_s=max_rpgf0_s)
CALL get_rho0_mpole(rho0_mpole=rho0_mpole, max_rpgf0_s=max_rpgf0_s, zet0_h=zet0_h)
igrid_zet0_s = gaussian_gridlevel(pw_env%gridlevel_info, 2.0_dp*zet0_h)
rho0_mpole%igrid_zet0_s = igrid_zet0_s
END IF
ngrid_level = SIZE(radius)
nkind = SIZE(qs_kind_set)
! Find the grid level suitable for rho0_soft
IF (dft_control%qs_control%gapw) THEN
igrid_zet0_s = gaussian_gridlevel(pw_env%gridlevel_info, 2.0_dp*zet0)
rho0_mpole%igrid_zet0_s = igrid_zet0_s
ELSE
igrid_zet0_s = 0
END IF
! try to predict the maximum radius of the gaussians to be mapped on the grid
! up to now, it is not yet very good
radius = 0.0_dp
Expand Down
14 changes: 9 additions & 5 deletions src/qs_environment.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -129,6 +129,7 @@ MODULE qs_environment
write_qs_particle_coordinates,&
write_structure_data
USE particle_types, ONLY: particle_type
USE pw_env_types, ONLY: pw_env_type
USE qmmm_types_low, ONLY: qmmm_env_qm_type
USE qs_dftb_parameters, ONLY: qs_dftb_param_init
USE qs_dftb_types, ONLY: qs_dftb_pairpot_type
Expand Down Expand Up @@ -161,6 +162,7 @@ MODULE qs_environment
qs_ks_env_type,&
qs_ks_release,&
set_ks_env
USE qs_local_rho_types, ONLY: local_rho_type
USE qs_mo_types, ONLY: allocate_mo_set,&
mo_set_p_type
USE qs_rho0_ggrid, ONLY: rho0_s_grid_create
Expand Down Expand Up @@ -519,12 +521,14 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
TYPE(gto_basis_set_type), POINTER :: aux_fit_basis, lri_aux_basis, &
ri_hfx_basis, ri_xas_basis, &
tmp_basis_set
TYPE(local_rho_type), POINTER :: local_rho_set
TYPE(lri_environment_type), POINTER :: lri_env
TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos, mos_aux_fit
TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
TYPE(nddo_mpole_type), POINTER :: se_nddo_mpole
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(pw_env_type), POINTER :: pw_env
TYPE(qs_control_type), POINTER :: qs_control
TYPE(qs_dftb_pairpot_type), DIMENSION(:, :), &
POINTER :: dftb_potential
Expand Down Expand Up @@ -989,9 +993,9 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
CALL init_rho_atom(rho_atom_set, qs_kind_set, qs_env, gapw_control)
CALL set_qs_env(qs_env=qs_env, rho_atom_set=rho_atom_set)
IF (dft_control%qs_control%method_id /= do_method_gapw_xc) THEN
CALL get_qs_env(qs_env=qs_env, natom=natom)
CALL get_qs_env(qs_env=qs_env, local_rho_set=local_rho_set, natom=natom)
! *** Allocate and initialize the compensation density rho0 ***
CALL init_rho0(qs_env, gapw_control)
CALL init_rho0(local_rho_set, qs_env, gapw_control, .FALSE.)
! *** Allocate and Initialize the local coulomb term ***
CALL init_coulomb_local(qs_env%hartree_local, natom)
END IF
Expand Down Expand Up @@ -1491,10 +1495,10 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
! Sets up pw_env, qs_charges, mpools ...
CALL qs_env_setup(qs_env)
! Allocate and Initialie rho0 soft on the global grid
! Allocate and initialise rho0 soft on the global grid
IF (dft_control%qs_control%method_id == do_method_gapw) THEN
CALL get_qs_env(qs_env=qs_env, rho0_mpole=rho0_mpole)
CALL rho0_s_grid_create(qs_env, rho0_mpole)
CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, rho0_mpole=rho0_mpole)
CALL rho0_s_grid_create(pw_env, rho0_mpole)
END IF
IF (output_unit > 0) CALL m_flush(output_unit)
Expand Down
33 changes: 6 additions & 27 deletions src/qs_environment_methods.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,6 +57,7 @@ MODULE qs_environment_methods
set_ks_env
USE qs_matrix_pools, ONLY: mpools_rebuild_fm_pools
USE qs_outer_scf, ONLY: outer_loop_variables_count
USE qs_rho0_ggrid, ONLY: rho0_s_grid_create
USE qs_rho0_types, ONLY: rho0_mpole_type
USE scf_control_types, ONLY: scf_control_type
#include "./base/base_uses.f90"
Expand Down Expand Up @@ -215,9 +216,9 @@ SUBROUTINE qs_env_rebuild_pw_env(qs_env)
TYPE(ewald_environment_type), POINTER :: ewald_env
TYPE(ewald_pw_type), POINTER :: ewald_pw
TYPE(pw_env_type), POINTER :: new_pw_env
TYPE(pw_p_type), POINTER :: embed_pot, external_vxc, rho0_gs, &
rho0_rs, rho_core, rho_nlcc, &
rho_nlcc_g, spin_embed_pot, vee, vppl
TYPE(pw_p_type), POINTER :: embed_pot, external_vxc, rho_core, &
rho_nlcc, rho_nlcc_g, spin_embed_pot, &
vee, vppl
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
TYPE(pw_type), POINTER :: v_hartree_rspace
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
Expand All @@ -227,7 +228,7 @@ SUBROUTINE qs_env_rebuild_pw_env(qs_env)
CALL timeset(routineN, handle)
! rebuild pw_env
NULLIFY (dft_control, cell, ks_env, v_hartree_rspace, auxbas_pw_pool)
NULLIFY (rho0_mpole, rho0_gs, rho0_rs)
NULLIFY (rho0_mpole)
NULLIFY (ewald_env, ewald_pw, new_pw_env, external_vxc, rho_core, rho_nlcc, rho_nlcc_g, vee, vppl, &
embed_pot, spin_embed_pot)

Expand Down Expand Up @@ -263,29 +264,7 @@ SUBROUTINE qs_env_rebuild_pw_env(qs_env)
CALL set_ks_env(ks_env, rho_core=rho_core)
END IF
CALL get_qs_env(qs_env=qs_env, rho0_mpole=rho0_mpole)
IF (ASSOCIATED(rho0_mpole)) THEN
rho0_rs => rho0_mpole%rho0_s_rs
rho0_gs => rho0_mpole%rho0_s_gs
IF (ASSOCIATED(rho0_rs)) THEN
CALL pw_release(rho0_rs%pw)
DEALLOCATE (rho0_rs)
END IF
ALLOCATE (rho0_rs)
CALL pw_env_get(new_pw_env, auxbas_pw_pool=auxbas_pw_pool)
CALL pw_pool_create_pw(auxbas_pw_pool, rho0_rs%pw, &
use_data=REALDATA3D, in_space=REALSPACE)
rho0_mpole%rho0_s_rs => rho0_rs

IF (ASSOCIATED(rho0_gs)) THEN
CALL pw_release(rho0_gs%pw)
DEALLOCATE (rho0_gs)
END IF
ALLOCATE (rho0_gs)
CALL pw_env_get(new_pw_env, auxbas_pw_pool=auxbas_pw_pool)
CALL pw_pool_create_pw(auxbas_pw_pool, rho0_gs%pw, &
use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE)
rho0_mpole%rho0_s_gs => rho0_gs
END IF
CALL rho0_s_grid_create(new_pw_env, rho0_mpole)
ELSE IF (dft_control%qs_control%semi_empirical) THEN
IF (dft_control%qs_control%se_control%do_ewald .OR. &
dft_control%qs_control%se_control%do_ewald_gks) THEN
Expand Down
8 changes: 6 additions & 2 deletions src/qs_environment_types.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -245,7 +245,7 @@ MODULE qs_environment_types
TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos, mos_aux_fit
TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: mo_derivs_aux_fit
TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: mo_loc_history
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mo_derivs
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mo_derivs
TYPE(scf_control_type), POINTER :: scf_control
TYPE(rel_control_type), POINTER :: rel_control
! ZMP adding variables
Expand Down Expand Up @@ -424,6 +424,7 @@ MODULE qs_environment_types
!> \param subsys ...
!> \param cp_subsys ...
!> \param oce ...
!> \param local_rho_set ...
!> \param rho_atom_set ...
!> \param task_list ...
!> \param task_list_aux_fit ...
Expand Down Expand Up @@ -509,7 +510,7 @@ SUBROUTINE get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, ce
mpools, mpools_aux_fit, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, &
vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, id_nr, local_particles, &
local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, &
molecule_set, subsys, cp_subsys, oce, rho_atom_set, &
molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, &
task_list, &
task_list_aux_fit, &
task_list_soft, &
Expand Down Expand Up @@ -591,6 +592,7 @@ SUBROUTINE get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, ce
TYPE(qs_subsys_type), OPTIONAL, POINTER :: subsys
TYPE(cp_subsys_type), OPTIONAL, POINTER :: cp_subsys
TYPE(oce_matrix_type), OPTIONAL, POINTER :: oce
TYPE(local_rho_type), OPTIONAL, POINTER :: local_rho_set
TYPE(rho_atom_type), DIMENSION(:), OPTIONAL, &
POINTER :: rho_atom_set
TYPE(task_list_type), OPTIONAL, POINTER :: task_list, task_list_aux_fit, &
Expand Down Expand Up @@ -729,6 +731,8 @@ SUBROUTINE get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, ce
! Resp charges
IF (PRESENT(rhs)) rhs => qs_env%rhs

IF (PRESENT(local_rho_set)) &
local_rho_set => qs_env%local_rho_set
IF (PRESENT(rho_atom_set)) &
CALL get_local_rho(qs_env%local_rho_set, rho_atom_set=rho_atom_set)
IF (PRESENT(rho0_atom_set)) &
Expand Down
9 changes: 6 additions & 3 deletions src/qs_p_env_methods.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,6 +61,7 @@ MODULE qs_p_env_methods
section_vals_type
USE kinds, ONLY: dp
USE preconditioner_types, ONLY: init_preconditioner
USE pw_env_types, ONLY: pw_env_type
USE qs_energy_types, ONLY: qs_energy_type
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
Expand All @@ -79,6 +80,7 @@ MODULE qs_p_env_methods
USE qs_mo_types, ONLY: get_mo_set,&
mo_set_p_type
USE qs_p_env_types, ONLY: qs_p_env_type
USE qs_rho0_ggrid, ONLY: rho0_s_grid_create
USE qs_rho0_methods, ONLY: init_rho0
USE qs_rho_atom_methods, ONLY: allocate_rho_atom_internals
USE qs_rho_methods, ONLY: qs_rho_rebuild,&
Expand Down Expand Up @@ -137,6 +139,7 @@ SUBROUTINE p_env_create(p_env, qs_env, kpp1_env, p1_option, &
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
TYPE(dft_control_type), POINTER :: dft_control
TYPE(pw_env_type), POINTER :: pw_env
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set

! code
Expand Down Expand Up @@ -258,10 +261,10 @@ SUBROUTINE p_env_create(p_env, qs_env, kpp1_env, p1_option, &
!----------------------!
IF (dft_control%qs_control%gapw) THEN
CALL local_rho_set_create(p_env%local_rho_set)
CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
CALL get_qs_env(qs_env, pw_env=pw_env, qs_kind_set=qs_kind_set)
CALL allocate_rho_atom_internals(qs_env, p_env%local_rho_set%rho_atom_set, qs_kind_set)
CALL init_rho0(qs_env, dft_control%qs_control%gapw_control, &
.TRUE., p_env%local_rho_set)
CALL init_rho0(p_env%local_rho_set, qs_env, dft_control%qs_control%gapw_control, .TRUE.)
CALL rho0_s_grid_create(pw_env, p_env%local_rho_set%rho0_mpole)
CALL hartree_local_create(p_env%hartree_local)
CALL get_qs_env(qs_env=qs_env, natom=natom)
CALL init_coulomb_local(p_env%hartree_local, natom)
Expand Down
28 changes: 12 additions & 16 deletions src/qs_rho0_ggrid.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ MODULE qs_rho0_ggrid
USE cp_control_types, ONLY: dft_control_type
USE cp_para_types, ONLY: cp_para_env_type
USE cube_utils, ONLY: cube_info_type
USE gaussian_gridlevels, ONLY: gaussian_gridlevel
USE input_constants, ONLY: tddfpt_singlet
USE kinds, ONLY: dp,&
int_8
Expand Down Expand Up @@ -279,33 +280,27 @@ END SUBROUTINE put_rho0_on_grid

! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param pw_env ...
!> \param rho0_mpole ...
!> \param tddft ...
! **************************************************************************************************
SUBROUTINE rho0_s_grid_create(qs_env, rho0_mpole, tddft)
SUBROUTINE rho0_s_grid_create(pw_env, rho0_mpole)

TYPE(qs_environment_type), POINTER :: qs_env
TYPE(pw_env_type), POINTER :: pw_env
TYPE(rho0_mpole_type), POINTER :: rho0_mpole
LOGICAL, INTENT(IN), OPTIONAL :: tddft

CHARACTER(len=*), PARAMETER :: routineN = 'rho0_s_grid_create', &
routineP = moduleN//':'//routineN

LOGICAL :: my_tddft
TYPE(pw_env_type), POINTER :: new_pw_env
INTEGER :: handle
TYPE(pw_p_type), POINTER :: rho_gs, rho_rs
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool

my_tddft = .FALSE.
IF (PRESENT(tddft)) my_tddft = tddft
CALL timeset(routineN, handle)

NULLIFY (new_pw_env)
CALL get_qs_env(qs_env, pw_env=new_pw_env)
CPASSERT(ASSOCIATED(new_pw_env))
CPASSERT(ASSOCIATED(pw_env))

NULLIFY (auxbas_pw_pool)
CALL pw_env_get(new_pw_env, auxbas_pw_pool=auxbas_pw_pool)
CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
CPASSERT(ASSOCIATED(auxbas_pw_pool))

! reallocate rho0 on the global grid in real and reciprocal space
Expand Down Expand Up @@ -335,9 +330,10 @@ SUBROUTINE rho0_s_grid_create(qs_env, rho0_mpole, tddft)
rho_gs%pw%in_space = RECIPROCALSPACE
rho0_mpole%rho0_s_gs => rho_gs

IF (my_tddft) THEN
rho0_mpole%igrid_zet0_s = qs_env%local_rho_set%rho0_mpole%igrid_zet0_s
END IF
! Find the grid level suitable for rho0_soft
rho0_mpole%igrid_zet0_s = gaussian_gridlevel(pw_env%gridlevel_info, 2.0_dp*rho0_mpole%zet0_h)

CALL timestop(handle)

END SUBROUTINE rho0_s_grid_create

Expand Down

0 comments on commit 5068987

Please sign in to comment.