Skip to content

Commit

Permalink
Refactor and extent WFC gradient and Linres code
Browse files Browse the repository at this point in the history 
Refactor matrix treatment in WFC and Linres code
Refactor some communication routines in WFC
Refactor spin treatment in WFC
Use Linres routines for MP2 gradients
Add XC and ADMM kernel to MP2 gradients
Remove unused variables and dead code
Add gradient and stress tensor for double-hybrid functionals
  • Loading branch information
@fstein93
Frederick Stein authored and fstein93 committed Sep 21, 2021
1 parent d0bd076 commit 549e906
Show file tree
Hide file tree
Showing 61 changed files with 6,527 additions and 9,577 deletions.
25 changes: 24 additions & 1 deletion src/common/cp_para_env.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@ MODULE cp_para_env
USE cp_para_types, ONLY: cp_para_cart_type,&
cp_para_env_type
USE message_passing, ONLY: mp_comm_free,&
mp_comm_split_direct,&
mp_environ
#include "../base/base_uses.f90"

Expand All @@ -25,7 +26,7 @@ MODULE cp_para_env
LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_para_env'

PUBLIC :: cp_para_env_retain, cp_para_env_release, cp_para_env_create
PUBLIC :: cp_para_env_retain, cp_para_env_release, cp_para_env_create, cp_para_env_split
PUBLIC :: cp_cart_create, cp_cart_release
!***
CONTAINS
Expand Down Expand Up @@ -82,6 +83,28 @@ SUBROUTINE cp_para_env_retain(para_env)
para_env%ref_count = para_env%ref_count + 1
END SUBROUTINE cp_para_env_retain

! **************************************************************************************************
!> \brief splits the communicator of the template para_env according to the provided color
!> \param para_env para_env with new sub communicator
!> \param para_env_template para_env to be split
!> \param color all processes with same color will belong to the same sub communicator
! **************************************************************************************************
SUBROUTINE cp_para_env_split(para_env, para_env_template, color)
TYPE(cp_para_env_type), POINTER :: para_env, para_env_template
INTEGER, INTENT(IN) :: color

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

INTEGER :: comm

CPASSERT(ASSOCIATED(para_env_template))
CALL mp_comm_split_direct(para_env_template%group, comm, color)
NULLIFY (para_env)
CALL cp_para_env_create(para_env, comm)

END SUBROUTINE cp_para_env_split

! **************************************************************************************************
!> \brief releases the para object (to be called when you don't want anymore
!> the shared copy of this object)
Expand Down
4 changes: 2 additions & 2 deletions src/efield_utils.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ MODULE efield_utils
custom_env,&
gaussian_env,&
ramp_env
USE kahan_sum, ONLY: accurate_sum
USE kahan_sum, ONLY: accurate_dot_product
USE kinds, ONLY: dp
USE mathconstants, ONLY: pi
USE message_passing, ONLY: mp_sum
Expand Down Expand Up @@ -107,7 +107,7 @@ SUBROUTINE efield_potential(qs_env, v_efield_rspace)
END DO
efield_ener = 0.0_dp
DO i = 1, dft_control%nspins
efield_ener = efield_ener + accurate_sum(v_efield_rspace%pw%cr3d*rho_r(i)%pw%cr3d)*dvol
efield_ener = efield_ener + accurate_dot_product(v_efield_rspace%pw%cr3d, rho_r(i)%pw%cr3d)*dvol
END DO
CALL mp_sum(efield_ener, para_env%group)
energy%efield = efield_ener
Expand Down
125 changes: 81 additions & 44 deletions src/energy_corrections.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -174,7 +174,7 @@ MODULE energy_corrections

LOGICAL, PARAMETER :: debug_forces = .TRUE.

PUBLIC :: energy_correction
PUBLIC :: energy_correction, calc_ks_response_potentials

CONTAINS

Expand Down Expand Up @@ -778,8 +778,7 @@ SUBROUTINE ec_build_ks_matrix_force(qs_env, ec_env)

INTEGER :: handle, i, iounit, ispin, nao, natom, &
norb, nspins
INTEGER, DIMENSION(2, 3) :: bo
LOGICAL :: lsd, use_virial
LOGICAL :: use_virial
REAL(dp) :: dehartree, eexc, eovrl, focc, &
total_rhoout
REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ftot
Expand All @@ -791,25 +790,23 @@ SUBROUTINE ec_build_ks_matrix_force(qs_env, ec_env)
TYPE(cp_fm_type), POINTER :: mo_coeff
TYPE(cp_logger_type), POINTER :: logger
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s, scrm
TYPE(dbcsr_p_type) :: scrm
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
TYPE(dft_control_type), POINTER :: dft_control
TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_p_type) :: dv_hartree_rspace, rho_tot_gspace, &
v_hartree_gspace, v_hartree_rspace, &
vtot_rspace
TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_g, rho_r, rhoout_g, rhoout_r, tau_r, &
TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_g, rho_r, rhoout_g, rhoout_r, &
v_rspace, v_tau_rspace, v_xc
TYPE(pw_poisson_type), POINTER :: poisson_env
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(qs_rho_type), POINTER :: rho
TYPE(section_vals_type), POINTER :: input, xc_fun_section, xc_section
TYPE(section_vals_type), POINTER :: input, xc_section
TYPE(virial_type), POINTER :: virial
TYPE(xc_derivative_set_type), POINTER :: deriv_set
TYPE(xc_rho_cflags_type) :: needs
TYPE(xc_rho_set_type), POINTER :: rho1_set, rho_set

CALL timeset(routineN, handle)

Expand Down Expand Up @@ -852,7 +849,7 @@ SUBROUTINE ec_build_ks_matrix_force(qs_env, ec_env)
! calculate output density on grid
! rho_in(R): CALL qs_rho_get(rho, rho_r=rho_r)
! rho_in(G): CALL qs_rho_get(rho, rho_g=rho_g)
CALL qs_rho_get(rho, rho_r=rho_r, rho_g=rho_g, tau_r=tau_r)
CALL qs_rho_get(rho, rho_r=rho_r, rho_g=rho_g)
ALLOCATE (rhoout_r(nspins), rhoout_g(nspins))
DO ispin = 1, nspins
NULLIFY (rhoout_r(ispin)%pw, rhoout_g(ispin)%pw)
Expand Down Expand Up @@ -908,28 +905,7 @@ SUBROUTINE ec_build_ks_matrix_force(qs_env, ec_env)
END DO
CALL get_qs_env(qs_env, input=input)
xc_section => ec_env%xc_section
lsd = (nspins == 2)
NULLIFY (deriv_set, rho_set, rho1_set)
CALL xc_prep_2nd_deriv(deriv_set, rho_set, rho_r, auxbas_pw_pool, xc_section=xc_section)
bo = rhoout_r(1)%pw%pw_grid%bounds_local
CALL xc_rho_set_create(rho1_set, bo, &
rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))

xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.)

! calculate the arguments needed by the functionals
CALL xc_rho_set_update(rho1_set, rhoout_r, rhoout_g, tau_r, needs, &
section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
auxbas_pw_pool)
CALL xc_calc_2nd_deriv(v_xc, deriv_set, rho_set, &
rho1_set, auxbas_pw_pool, xc_section=xc_section)
CALL xc_dset_release(deriv_set)
CALL xc_rho_set_release(rho_set)
CALL xc_rho_set_release(rho1_set)
CALL calc_ks_response_potentials(v_xc, xc_section, qs_env, rhoout_r, rhoout_g)
!
CALL get_qs_env(qs_env=qs_env, rho=rho, matrix_s_kp=matrix_s)
CALL dbcsr_allocate_matrix_set(ec_env%matrix_hz, nspins)
Expand Down Expand Up @@ -993,11 +969,10 @@ SUBROUTINE ec_build_ks_matrix_force(qs_env, ec_env)
CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_tot_gspace%pw)

! initialize src matrix
NULLIFY (scrm)
CALL dbcsr_allocate_matrix_set(scrm, 1, 1)
ALLOCATE (scrm(1, 1)%matrix)
CALL dbcsr_create(scrm(1, 1)%matrix, template=ec_env%matrix_s(1, 1)%matrix)
CALL cp_dbcsr_alloc_block_from_nbl(scrm(1, 1)%matrix, ec_env%sab_orb)
NULLIFY (scrm%matrix)
ALLOCATE (scrm%matrix)
CALL dbcsr_create(scrm%matrix, template=ec_env%matrix_s(1, 1)%matrix)
CALL cp_dbcsr_alloc_block_from_nbl(scrm%matrix, ec_env%sab_orb)

! v_rspace and v_tau_rspace are generated from the auxbas pool
NULLIFY (v_rspace, v_tau_rspace)
Expand All @@ -1013,18 +988,14 @@ SUBROUTINE ec_build_ks_matrix_force(qs_env, ec_env)
ENDDO
END IF

CALL qs_rho_get(rho, rho_r=rho_r)
IF (ASSOCIATED(v_tau_rspace)) THEN
CALL qs_rho_get(rho, tau_r=tau_r)
END IF
IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
DO ispin = 1, nspins
! Add v_hartree + v_xc = v_rspace
CALL pw_scale(v_rspace(ispin)%pw, v_rspace(ispin)%pw%pw_grid%dvol)
CALL pw_axpy(v_hartree_rspace%pw, v_rspace(ispin)%pw)
! integrate over potential <a|V|b>
CALL integrate_v_rspace(v_rspace=v_rspace(ispin), &
hmat=scrm(1, 1), &
hmat=scrm, &
pmat=ec_env%matrix_p(ispin, 1), &
qs_env=qs_env, &
calculate_forces=.TRUE., &
Expand All @@ -1034,7 +1005,7 @@ SUBROUTINE ec_build_ks_matrix_force(qs_env, ec_env)
IF (ASSOCIATED(v_tau_rspace)) THEN
! integrate over Tau-potential <nabla.a|V|nabla.b>
CALL pw_scale(v_tau_rspace(ispin)%pw, v_tau_rspace(ispin)%pw%pw_grid%dvol)
CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), hmat=scrm(1, 1), &
CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), hmat=scrm, &
pmat=ec_env%matrix_p(ispin, 1), &
qs_env=qs_env, calculate_forces=.TRUE., compute_tau=.TRUE., &
basis_type="HARRIS", &
Expand All @@ -1049,7 +1020,8 @@ SUBROUTINE ec_build_ks_matrix_force(qs_env, ec_env)
END IF

! delete scr matrix
CALL dbcsr_deallocate_matrix_set(scrm)
CALL dbcsr_release(scrm%matrix)
DEALLOCATE (scrm%matrix)

! return pw grids
CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_rspace%pw)
Expand Down Expand Up @@ -1085,6 +1057,71 @@ SUBROUTINE ec_build_ks_matrix_force(qs_env, ec_env)

END SUBROUTINE ec_build_ks_matrix_force

! **************************************************************************************************
!> \brief ...
!> \param v_xc ...
!> \param xc_section ...
!> \param qs_env ...
!> \param rho1_r ...
!> \param rho1_g ...
! **************************************************************************************************
SUBROUTINE calc_ks_response_potentials(v_xc, xc_section, qs_env, rho1_r, rho1_g, compute_virial, virial)
TYPE(pw_p_type), DIMENSION(:), INTENT(IN), POINTER :: v_xc
TYPE(section_vals_type), INTENT(IN), POINTER :: xc_section
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(pw_p_type), DIMENSION(:), INTENT(IN), POINTER :: rho1_r, rho1_g
LOGICAL, INTENT(IN), OPTIONAL :: compute_virial
REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), OPTIONAL :: virial

CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_ks_response_potentials'

INTEGER :: handle, nspins
INTEGER, DIMENSION(2, 3) :: bo
LOGICAL :: lsd
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_r, tau_r
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(qs_rho_type), POINTER :: rho
TYPE(section_vals_type), POINTER :: xc_fun_section
TYPE(xc_derivative_set_type), POINTER :: deriv_set
TYPE(xc_rho_cflags_type) :: needs
TYPE(xc_rho_set_type), POINTER :: rho1_set, rho_set

CALL timeset(routineN, handle)

CALL get_qs_env(qs_env, ks_env=ks_env, rho=rho, pw_env=pw_env)
CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
CALL qs_rho_get(rho, rho_r=rho_r, tau_r=tau_r)

nspins = SIZE(v_xc, 1)
lsd = (nspins == 2)
NULLIFY (deriv_set, rho_set, rho1_set)
CALL xc_prep_2nd_deriv(deriv_set, rho_set, rho_r, auxbas_pw_pool, xc_section=xc_section)
bo = rho1_r(1)%pw%pw_grid%bounds_local
CALL xc_rho_set_create(rho1_set, bo, &
rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))

xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.)

! calculate the arguments needed by the functionals
CALL xc_rho_set_update(rho1_set, rho1_r, rho1_g, tau_r, needs, &
section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
auxbas_pw_pool)
CALL xc_calc_2nd_deriv(v_xc, deriv_set, rho_set, &
rho1_set, auxbas_pw_pool, xc_section=xc_section, compute_virial=compute_virial, virial_xc=virial)
CALL xc_dset_release(deriv_set)
CALL xc_rho_set_release(rho_set)
CALL xc_rho_set_release(rho1_set)

CALL timestop(handle)

END SUBROUTINE calc_ks_response_potentials

! **************************************************************************************************
!> \brief Solve KS equation for a given matrix
!> \param qs_env ...
Expand Down
4 changes: 2 additions & 2 deletions src/fm/cp_fm_basic_linalg.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -741,8 +741,8 @@ SUBROUTINE cp_fm_trace_a0b0t0(matrix_a, matrix_b, trace)
trace = accurate_sum(a(1:nrow_local, 1:ncol_local)* &
REAL(b_sp(1:nrow_local, 1:ncol_local), dp))
ELSE
trace = accurate_sum(a(1:nrow_local, 1:ncol_local)* &
b(1:nrow_local, 1:ncol_local))
trace = accurate_dot_product(a(1:nrow_local, 1:ncol_local), &
b(1:nrow_local, 1:ncol_local))
ENDIF

CALL mp_sum(trace, group)
Expand Down
3 changes: 2 additions & 1 deletion src/fm/cp_fm_types.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,8 @@ MODULE cp_fm_types
cp_fm_read_unformatted, & ! reads a full matrix from an open unit
cp_fm_setup, & ! allows to set flags for fms
cp_fm_get_mm_type, &
cp_fm_write_info
cp_fm_write_info, &
cp_fm_to_fm_submat_general ! copy matrix across different contexts

PUBLIC :: cp_fm_indxg2p, &
cp_fm_indxg2l, &
Expand Down
28 changes: 7 additions & 21 deletions src/group_dist_types.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@
! **************************************************************************************************
MODULE group_dist_types
USE cp_para_types, ONLY: cp_para_env_type
USE message_passing, ONLY: mp_sum
USE message_passing, ONLY: mp_allgather
USE util, ONLY: get_limit
#include "./base/base_uses.f90"

Expand Down Expand Up @@ -144,19 +144,12 @@ SUBROUTINE create_group_dist_d1_i3(this, starts, ends, sizes, para_env)
CALL timeset(routineN, handle)

ALLOCATE (this%starts(0:para_env%num_pe - 1))
this%starts = 0
ALLOCATE (this%ends(0:para_env%num_pe - 1))
this%ends = 0
ALLOCATE (this%sizes(0:para_env%num_pe - 1))
this%sizes = 0

this%starts(para_env%mepos) = starts
this%ends(para_env%mepos) = ends
this%sizes(para_env%mepos) = sizes

CALL mp_sum(this%starts, para_env%group)
CALL mp_sum(this%ends, para_env%group)
CALL mp_sum(this%sizes, para_env%group)
CALL mp_allgather(starts, this%starts, para_env%group)
CALL mp_allgather(ends, this%ends, para_env%group)
CALL mp_allgather(sizes, this%sizes, para_env%group)

CALL timestop(handle)

Expand All @@ -180,19 +173,12 @@ SUBROUTINE create_group_dist_d1_gd(this, group_dist_ext, para_env)
CALL timeset(routineN, handle)

ALLOCATE (this%starts(0:para_env%num_pe - 1))
this%starts = 0
ALLOCATE (this%ends(0:para_env%num_pe - 1))
this%ends = 0
ALLOCATE (this%sizes(0:para_env%num_pe - 1))
this%sizes = 0

this%starts(para_env%mepos) = group_dist_ext%starts
this%ends(para_env%mepos) = group_dist_ext%ends
this%sizes(para_env%mepos) = group_dist_ext%sizes

CALL mp_sum(this%starts, para_env%group)
CALL mp_sum(this%ends, para_env%group)
CALL mp_sum(this%sizes, para_env%group)
CALL mp_allgather(group_dist_ext%starts, this%starts, para_env%group)
CALL mp_allgather(group_dist_ext%ends, this%ends, para_env%group)
CALL mp_allgather(group_dist_ext%sizes, this%sizes, para_env%group)

CALL timestop(handle)

Expand Down

0 comments on commit 549e906

Please sign in to comment.