Skip to content
Permalink
Browse files

Enabled ADMM with GAPW

  • Loading branch information
abussy authored and pseewald committed Dec 13, 2019
1 parent 8a41e87 commit 20168dcf5cca43901e8533cd0d19ccf4c61e7435
Showing with 1,961 additions and 315 deletions.
  1. +8 −3 src/admm_dm_methods.F
  2. +310 −25 src/admm_methods.F
  3. +66 −2 src/admm_types.F
  4. +4 −1 src/aobasis/basis_set_container_types.F
  5. +123 −8 src/hfx_admm_utils.F
  6. +3 −1 src/mp2.F
  7. +2 −2 src/mp2_gpw.F
  8. +1 −1 src/mp2_integrals.F
  9. +3 −5 src/mp2_ri_2c.F
  10. +1 −1 src/optbas_fenv_manipulation.F
  11. +8 −10 src/qs_collocate_density.F
  12. +2 −0 src/qs_energy_types.F
  13. +6 −2 src/qs_environment.F
  14. +47 −42 src/qs_gapw_densities.F
  15. +6 −10 src/qs_integrate_potential_product.F
  16. +3 −8 src/qs_interactions.F
  17. +17 −8 src/qs_kind_types.F
  18. +2 −1 src/qs_kpp1_env_methods.F
  19. +54 −42 src/qs_ks_atom.F
  20. +27 −16 src/qs_ks_methods.F
  21. +25 −7 src/qs_ks_utils.F
  22. +3 −3 src/qs_linres_epr_nablavks.F
  23. +6 −3 src/qs_linres_methods.F
  24. +4 −1 src/qs_p_env_methods.F
  25. +8 −11 src/qs_rho0_methods.F
  26. +75 −60 src/qs_rho_atom_methods.F
  27. +6 −4 src/qs_scf_output.F
  28. +48 −31 src/qs_vxc_atom.F
  29. +11 −3 src/rtp_admm_methods.F
  30. +245 −1 src/subsys/external_potential_types.F
  31. +1 −1 src/xc_adiabatic_utils.F
  32. +2 −2 src/xc_pot_saop.F
  33. +84 −0 tests/QS/regtest-admm-gapw/CH3-BP-CAUCHY.inp
  34. +84 −0 tests/QS/regtest-admm-gapw/CH3-BP-MO_DIAG.inp
  35. +84 −0 tests/QS/regtest-admm-gapw/CH3-BP-MO_NO_DIAG.inp
  36. +84 −0 tests/QS/regtest-admm-gapw/CH4-BP-CAUCHY_SUBSPACE.inp
  37. +79 −0 tests/QS/regtest-admm-gapw/H2+-BLOCKED-PURIFY-FULL.inp
  38. +77 −0 tests/QS/regtest-admm-gapw/H2-geo-ALL.inp
  39. +69 −0 tests/QS/regtest-admm-gapw/H2O-ADMMP-OPTX.inp
  40. +91 −0 tests/QS/regtest-admm-gapw/MD-1.inp
  41. +92 −0 tests/QS/regtest-admm-gapw/MD-2-ALL.inp
  42. +72 −0 tests/QS/regtest-admm-gapw/O2-triplet-ADMMS.inp
  43. +17 −0 tests/QS/regtest-admm-gapw/TEST_FILES
  44. 0 tests/QS/regtest-admm-gapw/TEST_FILES_RESET
  45. +1 −0 tests/TEST_DIRS
@@ -31,6 +31,7 @@ MODULE admm_dm_methods
USE qs_rho_types, ONLY: qs_rho_get,&
qs_rho_set,&
qs_rho_type
USE task_list_types, ONLY: task_list_type
#include "./base/base_uses.f90"

IMPLICIT NONE
@@ -247,13 +248,16 @@ SUBROUTINE update_rho_aux(ks_env)
TYPE(dft_control_type), POINTER :: dft_control
TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_g_aux, rho_r_aux
TYPE(qs_rho_type), POINTER :: rho_aux
TYPE(task_list_type), POINTER :: task_list_aux_fit

NULLIFY (dft_control, admm_dm, rho_aux, rho_ao_aux, rho_r_aux, rho_g_aux, tot_rho_r_aux)
NULLIFY (dft_control, admm_dm, rho_aux, rho_ao_aux, rho_r_aux, rho_g_aux, tot_rho_r_aux, &
task_list_aux_fit)

CALL get_ks_env(ks_env, &
admm_dm=admm_dm, &
dft_control=dft_control, &
rho_aux_fit=rho_aux)
rho_aux_fit=rho_aux, &
task_list_aux_fit=task_list_aux_fit)

CALL qs_rho_get(rho_aux, &
rho_ao=rho_ao_aux, &
@@ -268,7 +272,8 @@ SUBROUTINE update_rho_aux(ks_env)
rho_gspace=rho_g_aux(ispin), &
total_rho=tot_rho_r_aux(ispin), &
soft_valid=.FALSE., &
basis_type="AUX_FIT")
basis_type="AUX_FIT", &
task_list_external=task_list_aux_fit)
END DO

CALL qs_rho_set(rho_aux, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)

Large diffs are not rendered by default.

@@ -7,6 +7,7 @@
!> \brief Types and set/get functions for auxiliary denisty matrix methods
!> \par History
!> 05.2008 created [Manuel Guidon]
!> 12.2019 Made GAPW compatiblae [Augustin Bussy]
!> \author Manuel Guidon
! **************************************************************************************************
MODULE admm_types
@@ -27,13 +28,21 @@ MODULE admm_types
USE input_section_types, ONLY: section_vals_release,&
section_vals_type
USE kinds, ONLY: dp
USE qs_kind_types, ONLY: deallocate_qs_kind_set,&
qs_kind_type
USE qs_local_rho_types, ONLY: local_rho_set_release,&
local_rho_type
USE qs_mo_types, ONLY: get_mo_set,&
mo_set_p_type
USE qs_oce_types, ONLY: deallocate_oce_set,&
oce_matrix_type
USE task_list_types, ONLY: deallocate_task_list,&
task_list_type
#include "./base/base_uses.f90"

IMPLICIT NONE
PRIVATE
PUBLIC :: admm_env_create, admm_env_release, admm_type
PUBLIC :: admm_env_create, admm_env_release, admm_type, admm_gapw_type

CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'admm_types'

@@ -45,6 +54,20 @@ MODULE admm_types
TYPE(eigvals_type), POINTER :: eigvals
END TYPE

! **************************************************************************************************
!> \brief A subtype of the admm_env that contains the extra data needed for an ADMM GAPW calculation
!> \param admm_kind_set gets its own qs_kind set to store all relevant basis/grid/etc info
!> \param local_rho_set caontains soft and hard AUX_FIT atomoc densities
!> \param task_list the task list used for all soft density pw operations
!> \param oce stores the precomputed oce integrals
! **************************************************************************************************
TYPE admm_gapw_type
TYPE(qs_kind_type), DIMENSION(:), POINTER :: admm_kind_set => Null()
TYPE(local_rho_type), POINTER :: local_rho_set => Null()
TYPE(task_list_type), POINTER :: task_list => Null()
TYPE(oce_matrix_type), POINTER :: oce => Null()
END TYPE admm_gapw_type

! **************************************************************************************************
!> \brief stores some data used in wavefunction fitting
!> \param S overlap matrix for auxiliary fit basis set
@@ -67,6 +90,8 @@ MODULE admm_types
!> \param nmo number of molecular orbitals per spin
!> \param eigvals_lamda eigenvalues of lambda matrix
!> \param gsi contains ratio N_dens_m/N_aux_dens_m
!> \param admm_gapw_env the type containing ADMM GAPW specific data
!> \param do_gapw an internal logical switch for GAPW
!> \par History
!> 05.2008 created [Manuel Guidon]
!> \author Manuel Guidon
@@ -129,7 +154,9 @@ MODULE admm_types
REAL(KIND=dp), DIMENSION(3) :: aux_x_param
LOGICAL :: block_dm = .FALSE.
LOGICAL :: block_fit = .FALSE.
INTEGER, DIMENSION(:, :), POINTER :: block_map => Null()
INTEGER, DIMENSION(:, :), POINTER :: block_map => Null()
TYPE(admm_gapw_type), POINTER :: admm_gapw_env
LOGICAL :: do_gapw = .FALSE.
END TYPE

CONTAINS
@@ -322,6 +349,10 @@ SUBROUTINE admm_env_create(admm_env, admm_control, mos, mos_aux_fit, para_env, n
END DO
END DO
ENDIF

NULLIFY (admm_env%admm_gapw_env)
admm_env%do_gapw = .FALSE.

END SUBROUTINE admm_env_create

! **************************************************************************************************
@@ -418,9 +449,42 @@ SUBROUTINE admm_env_release(admm_env)
IF (ASSOCIATED(admm_env%xc_section_aux)) &
CALL section_vals_release(admm_env%xc_section_aux)

IF (ASSOCIATED(admm_env%admm_gapw_env)) CALL admm_gapw_env_release(admm_env%admm_gapw_env)

DEALLOCATE (admm_env)

END SUBROUTINE admm_env_release

! **************************************************************************************************
!> \brief Release the ADMM GAPW stuff
!> \param admm_gapw_env ...
! **************************************************************************************************
SUBROUTINE admm_gapw_env_release(admm_gapw_env)

TYPE(admm_gapw_type), POINTER :: admm_gapw_env

CHARACTER(LEN=*), PARAMETER :: routineN = 'admm_gapw_env_release', &
routineP = moduleN//':'//routineN

IF (ASSOCIATED(admm_gapw_env%admm_kind_set)) THEN
CALL deallocate_qs_kind_set(admm_gapw_env%admm_kind_set)
END IF

IF (ASSOCIATED(admm_gapw_env%local_rho_set)) THEN
CALL local_rho_set_release(admm_gapw_env%local_rho_set)
END IF

IF (ASSOCIATED(admm_gapw_env%task_list)) THEN
CALL deallocate_task_list(admm_gapw_env%task_list)
END IF

IF (ASSOCIATED(admm_gapw_env%oce)) THEN
CALL deallocate_oce_set(admm_gapw_env%oce)
END IF

DEALLOCATE (admm_gapw_env)

END SUBROUTINE admm_gapw_env_release

END MODULE admm_types

@@ -37,7 +37,8 @@ MODULE basis_set_container_types
aux_gw_basis = 110, &
ri_hxc_basis = 111, &
ri_k_basis = 112, &
ri_xas_basis = 113
ri_xas_basis = 113, &
aux_fit_soft_basis = 114
! **************************************************************************************************
TYPE basis_set_container_type
PRIVATE
@@ -104,6 +105,8 @@ FUNCTION get_basis_type(basis_set_type) RESULT(basis_type_nr)
basis_type_nr = lri_aux_basis
CASE ("AUX_FIT")
basis_type_nr = aux_fit_basis
CASE ("AUX_FIT_SOFT")
basis_type_nr = aux_fit_soft_basis
CASE ("SOFT")
basis_type_nr = soft_basis
CASE ("HARD")
@@ -9,15 +9,20 @@
!>
!> \par History
!> refactoring 03-2011 [MI]
!> Made GAPW compatible 12.2019 (A. Bussy)
!> \author MI
! **************************************************************************************************
MODULE hfx_admm_utils
USE admm_dm_types, ONLY: admm_dm_create,&
admm_dm_type
USE admm_methods, ONLY: scale_dm
USE admm_types, ONLY: admm_env_create,&
admm_gapw_type,&
admm_type
USE atomic_kind_types, ONLY: atomic_kind_type
USE basis_set_container_types, ONLY: add_basis_set_to_container
USE basis_set_types, ONLY: copy_gto_basis_set,&
gto_basis_set_type
USE cell_types, ONLY: cell_type
USE cp_control_types, ONLY: dft_control_type
USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
@@ -33,6 +38,7 @@ MODULE hfx_admm_utils
dbcsr_p_type,&
dbcsr_set,&
dbcsr_type
USE external_potential_types, ONLY: copy_potential
USE hfx_derivatives, ONLY: derivatives_four_center
USE hfx_energy_potential, ONLY: integrate_four_center
USE hfx_types, ONLY: hfx_type
@@ -71,11 +77,17 @@ MODULE hfx_admm_utils
qs_environment_type,&
set_qs_env
USE qs_force_types, ONLY: qs_force_type
USE qs_kind_types, ONLY: qs_kind_type
USE qs_interactions, ONLY: init_interaction_radii
USE qs_kind_types, ONLY: get_qs_kind,&
init_gapw_basis_set,&
init_gapw_nlcc,&
qs_kind_type
USE qs_ks_types, ONLY: qs_ks_env_type,&
set_ks_env
USE qs_local_rho_types, ONLY: local_rho_set_create
USE qs_mo_types, ONLY: get_mo_set,&
mo_set_p_type
USE qs_rho_atom_methods, ONLY: init_rho_atom
USE qs_rho_types, ONLY: qs_rho_get,&
qs_rho_type
USE rt_propagation_types, ONLY: rt_prop_type
@@ -157,12 +169,6 @@ SUBROUTINE hfx_admm_init(qs_env, do_mp2_hfx)
IF (.NOT. do_hfx) &
CPABORT("Wavefunction fitting requested without Hartree-Fock.")

! ** Method runs with GAPW only if no DFT exchange correction
IF (dft_control%qs_control%gapw .AND. &
dft_control%admm_control%aux_exch_func .NE. do_admm_aux_exch_func_none) THEN
CPABORT("ADMM only works with GAPW without DFT exchange correction")
END IF

CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
IF (n_rep_hf > 1) &
CPABORT("ADMM can handle only one HF section.")
@@ -177,8 +183,21 @@ SUBROUTINE hfx_admm_init(qs_env, do_mp2_hfx)
xc_section => section_vals_get_subs_vals(input, "DFT%XC")
CALL create_admm_xc_section(qs_env=qs_env, xc_section=xc_section, &
admm_env=admm_env)

! Initialize the GAPW data types
IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) &
CALL init_admm_gapw(qs_env)

ENDIF

IF (admm_env%do_gapw .AND. dft_control%do_admm_dm) THEN
CPABORT("GAPW ADMM not implemented for MCWEENY or NONE_DM purification.")
END IF

IF (admm_env%do_gapw .AND. qs_env%run_rtp) THEN
CPABORT("GAPW ADMM not implemented for real-time propagation.")
END IF

IF (dft_control%do_admm_dm .AND. .NOT. ASSOCIATED(admm_dm)) THEN
CALL admm_dm_create(admm_dm, dft_control%admm_control, nspins=nspins, natoms=natoms)
CALL set_ks_env(ks_env, admm_dm=admm_dm)
@@ -188,6 +207,95 @@ SUBROUTINE hfx_admm_init(qs_env, do_mp2_hfx)

END SUBROUTINE hfx_admm_init

! **************************************************************************************************
!> \brief Sets up the admm_gapw env
!> \param qs_env ...
! **************************************************************************************************
SUBROUTINE init_admm_gapw(qs_env)

TYPE(qs_environment_type), POINTER :: qs_env

CHARACTER(LEN=*), PARAMETER :: routineN = 'init_admm_gapw', routineP = moduleN//':'//routineN

INTEGER :: ikind, nkind
TYPE(admm_gapw_type), POINTER :: admm_gapw_env
TYPE(admm_type), POINTER :: admm_env
TYPE(dft_control_type), POINTER :: dft_control
TYPE(gto_basis_set_type), POINTER :: aux_fit_basis, aux_fit_soft_basis, &
orb_basis, soft_basis
TYPE(qs_kind_type), DIMENSION(:), POINTER :: admm_kind_set, qs_kind_set

NULLIFY (qs_kind_set, admm_kind_set, aux_fit_basis, orb_basis, dft_control, soft_basis, &
aux_fit_soft_basis)

CALL get_qs_env(qs_env, admm_env=admm_env, qs_kind_set=qs_kind_set, dft_control=dft_control)

admm_env%do_gapw = .TRUE.
ALLOCATE (admm_env%admm_gapw_env)
admm_gapw_env => admm_env%admm_gapw_env
NULLIFY (admm_gapw_env%local_rho_set)
NULLIFY (admm_gapw_env%admm_kind_set)
NULLIFY (admm_gapw_env%task_list)

!Create a new kind set for the ADMM stuff (paw_proj soft AUX_FIT basis, etc)
nkind = SIZE(qs_kind_set)
ALLOCATE (admm_gapw_env%admm_kind_set(nkind))
admm_kind_set => admm_gapw_env%admm_kind_set

!In this new kind set, we want the AUX_FIT basis to be known as ORB, such that GAPW routines work
DO ikind = 1, nkind
!copying over simple data of interest from qs_kind_set
admm_kind_set(ikind)%name = qs_kind_set(ikind)%name
admm_kind_set(ikind)%element_symbol = qs_kind_set(ikind)%element_symbol
admm_kind_set(ikind)%natom = qs_kind_set(ikind)%natom
admm_kind_set(ikind)%hard_radius = qs_kind_set(ikind)%hard_radius
admm_kind_set(ikind)%max_rad_local = qs_kind_set(ikind)%max_rad_local
admm_kind_set(ikind)%gpw_type_forced = qs_kind_set(ikind)%gpw_type_forced
admm_kind_set(ikind)%ngrid_rad = qs_kind_set(ikind)%ngrid_rad
admm_kind_set(ikind)%ngrid_ang = qs_kind_set(ikind)%ngrid_ang

!copying potentials of interest from qs_kind_set
IF (ASSOCIATED(qs_kind_set(ikind)%all_potential)) THEN
CALL copy_potential(qs_kind_set(ikind)%all_potential, admm_kind_set(ikind)%all_potential)
END IF
IF (ASSOCIATED(qs_kind_set(ikind)%gth_potential)) THEN
CALL copy_potential(qs_kind_set(ikind)%gth_potential, admm_kind_set(ikind)%gth_potential)
END IF
IF (ASSOCIATED(qs_kind_set(ikind)%sgp_potential)) THEN
CALL copy_potential(qs_kind_set(ikind)%sgp_potential, admm_kind_set(ikind)%sgp_potential)
END IF

NULLIFY (orb_basis)
CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_fit_basis, basis_type="AUX_FIT")
CALL copy_gto_basis_set(aux_fit_basis, orb_basis)
CALL add_basis_set_to_container(admm_kind_set(ikind)%basis_sets, orb_basis, "ORB")
END DO

!Create the corresponding soft basis set (and projectors)
CALL init_gapw_basis_set(admm_kind_set, dft_control%qs_control, qs_env%input, &
modify_qs_control=.FALSE.)

!Make sure the basis and the projectors are well initialized
CALL init_interaction_radii(dft_control%qs_control, admm_kind_set)

!We also init the atomic grids and harmonics
CALL local_rho_set_create(admm_gapw_env%local_rho_set)
CALL init_rho_atom(admm_gapw_env%local_rho_set%rho_atom_set, admm_kind_set, qs_env, &
dft_control%qs_control%gapw_control)

!Make sure that any NLCC potential is well initialized
CALL init_gapw_nlcc(admm_kind_set)

!Need to have access to the soft AUX_FIT basis from the qs_env => add it to the qs_kinds
DO ikind = 1, nkind
NULLIFY (aux_fit_soft_basis)
CALL get_qs_kind(admm_kind_set(ikind), soft_basis_set=soft_basis)
CALL copy_gto_basis_set(soft_basis, aux_fit_soft_basis)
CALL add_basis_set_to_container(qs_kind_set(ikind)%basis_sets, aux_fit_soft_basis, "AUX_FIT_SOFT")
END DO

END SUBROUTINE init_admm_gapw

! **************************************************************************************************
!> \brief Add the hfx contributions to the Hamiltonian
!>
@@ -276,7 +384,6 @@ SUBROUTINE hfx_ks_matrix(qs_env, matrix_ks, rho, energy, calculate_forces, &
IF (dft_control%do_admm) THEN
DO ispin = 1, nspins
CALL dbcsr_set(matrix_ks_aux_fit(ispin)%matrix, 0.0_dp)
CALL dbcsr_set(matrix_ks_aux_fit_hfx(ispin)%matrix, 0.0_dp)
END DO
END IF
DO ispin = 1, nspins
@@ -416,6 +523,14 @@ SUBROUTINE hfx_ks_matrix(qs_env, matrix_ks, rho, energy, calculate_forces, &
hf_energy, just_energy, calculate_forces, use_virial)
END IF ! do_adiabatic_rescaling

!update the hfx aux_fit matrixIF (dft_control%do_admm) THEN
IF (dft_control%do_admm) THEN
DO ispin = 1, nspins
CALL dbcsr_add(matrix_ks_aux_fit_hfx(ispin)%matrix, matrix_ks_aux_fit(ispin)%matrix, &
0.0_dp, 1.0_dp)
END DO
END IF

CALL timestop(handle)

END SUBROUTINE hfx_ks_matrix
@@ -1200,11 +1200,13 @@ SUBROUTINE calculate_exx(qs_env, unit_nr, do_gw, do_admm, E_ex_from_GW, t3)
END DO
! Remove the Exchange-correlation energy contributions from the total energy
energy%total = energy%total - (energy%exc + energy%exc1 + energy%ex + energy%exc_aux_fit)
energy%total = energy%total - (energy%exc + energy%exc1 + energy%ex + &
energy%exc_aux_fit + energy%exc1_aux_fit)
energy%exc = 0.0_dp
energy%exc1 = 0.0_dp
energy%exc_aux_fit = 0.0_dp
energy%exc1_aux_fit = 0.0_dp
energy%ex = 0.0_dp
! take the exact exchange energy from GW or calculate it

0 comments on commit 20168dc

Please sign in to comment.
You can’t perform that action at this time.