From c2221170b8e8553a8ea38fe4b47a5a8077401ec1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrg=20Hutter?= Date: Tue, 25 Sep 2018 07:02:03 +0000 Subject: [PATCH] KG: towards atomic integration grid svn-origin-rev: 18614 --- src/input_cp2k_dft.F | 7 + src/integration_grid_types.F | 179 ++++++++++++++++++ src/kg_correction.F | 6 + src/kg_environment.F | 109 ++++++++++- src/kg_environment_types.F | 25 ++- src/qs_energy_init.F | 31 +++- src/qs_environment.F | 16 +- src/qs_grid_atom.F | 347 ++++++++++++++++++++++++++++++++++- src/xc/xc_libxc.F | 12 +- 9 files changed, 697 insertions(+), 35 deletions(-) create mode 100644 src/integration_grid_types.F diff --git a/src/input_cp2k_dft.F b/src/input_cp2k_dft.F index 874316b7ba..474e4f3bca 100644 --- a/src/input_cp2k_dft.F +++ b/src/input_cp2k_dft.F @@ -5286,6 +5286,13 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, name="INTEGRATION_GRID", & + description="Grid [small,medium,large,huge]to be used for the TNADD integration.", & + usage="INTEGRATION_GRID MEDIUM", & + default_c_val="MEDIUM") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + CALL section_create(subsection, name="PRINT", & description="Print section", & n_keywords=0, n_subsections=1, repeats=.FALSE.) diff --git a/src/integration_grid_types.F b/src/integration_grid_types.F new file mode 100644 index 0000000000..e7c897faaa --- /dev/null +++ b/src/integration_grid_types.F @@ -0,0 +1,179 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright (C) 2000 - 2018 CP2K developers group ! +!--------------------------------------------------------------------------------------------------! +! ************************************************************************************************** +MODULE integration_grid_types + + USE kinds, ONLY: dp +#include "./base/base_uses.f90" + + IMPLICIT NONE + + PRIVATE + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'integration_grid_types' + + TYPE grid_batch_val_1d_type + INTEGER :: np1 + REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: val1d + END TYPE grid_batch_val_1d_type + + TYPE grid_batch_val_2d_type + INTEGER :: np1, np2 + REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: val2d + END TYPE grid_batch_val_2d_type + + TYPE gnlist_type + INTEGER, DIMENSION(:), ALLOCATABLE :: atom_list + REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: atom_pos + END TYPE gnlist_type + + TYPE grid_batch_info_type + INTEGER :: np + INTEGER :: ref_atom + INTEGER :: ibatch + TYPE(gnlist_type) :: gnlist + REAL(KIND=dp), DIMENSION(3) :: rcenter + REAL(KIND=dp) :: radius + REAL(dp), DIMENSION(:, :), ALLOCATABLE :: rco + REAL(dp), DIMENSION(:), ALLOCATABLE :: weight + REAL(dp), DIMENSION(:), ALLOCATABLE :: wref + REAL(dp), DIMENSION(:), ALLOCATABLE :: wsum + END TYPE grid_batch_info_type + + TYPE integration_grid_type + INTEGER :: nbatch + TYPE(grid_batch_info_type), DIMENSION(:), ALLOCATABLE :: grid_batch + END TYPE integration_grid_type + + TYPE integration_grid_value_type + INTEGER :: nbatch + TYPE(grid_batch_val_1d_type), DIMENSION(:), ALLOCATABLE :: grid_val_1d + TYPE(grid_batch_val_2d_type), DIMENSION(:), ALLOCATABLE :: grid_val_2d + END TYPE integration_grid_value_type + + PUBLIC :: integration_grid_type, allocate_intgrid, deallocate_intgrid + PUBLIC :: integration_grid_value_type, allocate_intgrid_val, deallocate_intgrid_val + +! ************************************************************************************************** + +CONTAINS + +! ************************************************************************************************** +!> \brief Initialize integration_grid_type +!> \param int_grid ... +!> \date 02.2018 +!> \param +!> \author JGH +!> \version 1.0 +! ************************************************************************************************** + SUBROUTINE allocate_intgrid(int_grid) + + TYPE(integration_grid_type), POINTER :: int_grid + + CHARACTER(len=*), PARAMETER :: routineN = 'allocate_intgrid', & + routineP = moduleN//':'//routineN + + IF (ASSOCIATED(int_grid)) CALL deallocate_intgrid(int_grid) + ALLOCATE (int_grid) + int_grid%nbatch = 0 + + END SUBROUTINE allocate_intgrid + +! ************************************************************************************************** +!> \brief Deallocate integration_grid_type +!> \param int_grid ... +!> \date 02.2018 +!> \param +!> \author JGH +!> \version 1.0 +! ************************************************************************************************** + SUBROUTINE deallocate_intgrid(int_grid) + TYPE(integration_grid_type), POINTER :: int_grid + + CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_intgrid', & + routineP = moduleN//':'//routineN + + INTEGER :: i + + IF (ASSOCIATED(int_grid)) THEN + IF (ALLOCATED(int_grid%grid_batch)) THEN + DO i = 1, int_grid%nbatch + IF (ALLOCATED(int_grid%grid_batch(i)%rco)) DEALLOCATE (int_grid%grid_batch(i)%rco) + IF (ALLOCATED(int_grid%grid_batch(i)%weight)) DEALLOCATE (int_grid%grid_batch(i)%weight) + IF (ALLOCATED(int_grid%grid_batch(i)%wref)) DEALLOCATE (int_grid%grid_batch(i)%wref) + IF (ALLOCATED(int_grid%grid_batch(i)%wsum)) DEALLOCATE (int_grid%grid_batch(i)%wsum) + ! + IF (ALLOCATED(int_grid%grid_batch(i)%gnlist%atom_list)) DEALLOCATE (int_grid%grid_batch(i)%gnlist%atom_list) + IF (ALLOCATED(int_grid%grid_batch(i)%gnlist%atom_pos)) DEALLOCATE (int_grid%grid_batch(i)%gnlist%atom_pos) + END DO + DEALLOCATE (int_grid%grid_batch) + END IF + DEALLOCATE (int_grid) + ELSE + CALL cp_abort(__LOCATION__, & + "The pointer int_grid is not associated and "// & + "cannot be deallocated") + END IF + END SUBROUTINE deallocate_intgrid + +! ************************************************************************************************** +!> \brief Initialize integration_grid_value_type +!> \param int_grid ... +!> \date 02.2018 +!> \param +!> \author JGH +!> \version 1.0 +! ************************************************************************************************** + SUBROUTINE allocate_intgrid_val(int_grid) + + TYPE(integration_grid_value_type), POINTER :: int_grid + + CHARACTER(len=*), PARAMETER :: routineN = 'allocate_intgrid_val', & + routineP = moduleN//':'//routineN + + IF (ASSOCIATED(int_grid)) CALL deallocate_intgrid_val(int_grid) + ALLOCATE (int_grid) + int_grid%nbatch = 0 + + END SUBROUTINE allocate_intgrid_val + +! ************************************************************************************************** +!> \brief Deallocate integration_grid_value_type +!> \param int_grid ... +!> \date 02.2018 +!> \param +!> \author JGH +!> \version 1.0 +! ************************************************************************************************** + SUBROUTINE deallocate_intgrid_val(int_grid) + TYPE(integration_grid_value_type), POINTER :: int_grid + + CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_intgrid_val', & + routineP = moduleN//':'//routineN + + INTEGER :: i + + IF (ASSOCIATED(int_grid)) THEN + IF (ALLOCATED(int_grid%grid_val_1d)) THEN + DO i = 1, int_grid%nbatch + IF (ALLOCATED(int_grid%grid_val_1d(i)%val1d)) DEALLOCATE (int_grid%grid_val_1d(i)%val1d) + END DO + DEALLOCATE (int_grid%grid_val_1d) + END IF + IF (ALLOCATED(int_grid%grid_val_2d)) THEN + DO i = 1, int_grid%nbatch + IF (ALLOCATED(int_grid%grid_val_2d(i)%val2d)) DEALLOCATE (int_grid%grid_val_2d(i)%val2d) + END DO + DEALLOCATE (int_grid%grid_val_2d) + END IF + DEALLOCATE (int_grid) + ELSE + CALL cp_abort(__LOCATION__, & + "The pointer int_grid is not associated and "// & + "cannot be deallocated") + END IF + END SUBROUTINE deallocate_intgrid_val + +END MODULE integration_grid_types diff --git a/src/kg_correction.F b/src/kg_correction.F index d40ed2f8bd..258a358d6e 100644 --- a/src/kg_correction.F +++ b/src/kg_correction.F @@ -486,6 +486,9 @@ CONTAINS ekin_mol = -ekin_imol xcvirial(1:3, 1:3) = 0.0_dp IF (use_virial) xcvirial(1:3, 1:3) = xcvirial(1:3, 1:3)+virial%pv_xc(1:3, 1:3) +!deb +! WRITE(6,*) " E KIN (full) ",-ekin_mol +!deb ! loop over all subsets ALLOCATE (atomlist(natom)) @@ -504,6 +507,9 @@ CONTAINS CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=kg_env%xc_section_kg, & vxc_rho=vxc_rho, vxc_tau=vxc_tau, exc=ekin_imol) ekin_mol = ekin_mol+ekin_imol +!deb +! WRITE(6,*) " E KIN (molecule) ",isub,ekin_imol +!deb DO ispin = 1, nspins vxc_rho(ispin)%pw%cr3d = -vxc_rho(ispin)%pw%cr3d*vxc_rho(ispin)%pw%pw_grid%dvol diff --git a/src/kg_environment.F b/src/kg_environment.F index 4a1948c748..a0184b9260 100644 --- a/src/kg_environment.F +++ b/src/kg_environment.F @@ -24,6 +24,9 @@ MODULE kg_environment USE cp_control_types, ONLY: dft_control_type USE cp_files, ONLY: close_file,& open_file + USE cp_log_handling, ONLY: cp_get_default_logger,& + cp_logger_get_default_io_unit,& + cp_logger_type USE cp_para_types, ONLY: cp_para_env_type USE distribution_1d_types, ONLY: distribution_1d_type USE distribution_2d_types, ONLY: distribution_2d_type @@ -39,6 +42,8 @@ MODULE kg_environment USE input_section_types, ONLY: section_vals_get_subs_vals,& section_vals_type,& section_vals_val_get + USE integration_grid_types, ONLY: allocate_intgrid,& + integration_grid_type USE kg_environment_types, ONLY: kg_environment_type USE kg_vertex_coloring_methods, ONLY: kg_vertex_coloring USE kinds, ONLY: dp,& @@ -58,6 +63,7 @@ MODULE kg_environment USE qs_dispersion_utils, ONLY: qs_dispersion_env_set USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type + USE qs_grid_atom, ONLY: initialize_atomic_grid USE qs_interactions, ONLY: init_interaction_radii_orb_basis USE qs_kind_types, ONLY: get_qs_kind,& qs_kind_type @@ -130,12 +136,17 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'init_kg_env', routineP = moduleN//':'//routineN - INTEGER :: handle, ikind, nkind - REAL(KIND=dp) :: eps_pgf_orb + CHARACTER(LEN=10) :: intgrid + INTEGER :: handle, i, iatom, ib, ikind, iunit, n, & + na, natom, nbatch, nkind, np, nr + INTEGER, ALLOCATABLE, DIMENSION(:, :) :: bid + REAL(KIND=dp) :: eps_pgf_orb, load, radb, rmax TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + TYPE(cp_logger_type), POINTER :: logger TYPE(cp_para_env_type), POINTER :: para_env TYPE(dft_control_type), POINTER :: dft_control TYPE(gto_basis_set_type), POINTER :: basis_set, harris_basis, lri_aux_basis + TYPE(integration_grid_type), POINTER :: ig_full, ig_mol TYPE(qs_dispersion_type), POINTER :: dispersion_env TYPE(qs_kind_type), POINTER :: qs_kind TYPE(section_vals_type), POINTER :: lri_section, nl_section, pp_section, & @@ -152,6 +163,9 @@ CONTAINS NULLIFY (kg_env%subset) NULLIFY (kg_env%tnadd_mat) NULLIFY (kg_env%lri_env) + NULLIFY (kg_env%int_grid_atom) + NULLIFY (kg_env%int_grid_molecules) + NULLIFY (kg_env%int_grid_full) NULLIFY (kg_env%lri_density) NULLIFY (kg_env%ec_env%sab_orb, kg_env%ec_env%sac_ppl, kg_env%ec_env%sap_ppnl) NULLIFY (kg_env%ec_env%matrix_ks, kg_env%ec_env%matrix_h, kg_env%ec_env%matrix_s) @@ -290,17 +304,106 @@ CONTAINS IF (kg_env%tnadd_method == kg_tnadd_embed_ri) THEN ! initialize the LRI environment - ! Check if LRI_AUX basis is available, auto-generate if needed + ! Check if LRI_AUX basis is available + rmax = 0.0_dp nkind = SIZE(qs_kind_set) DO ikind = 1, nkind qs_kind => qs_kind_set(ikind) NULLIFY (lri_aux_basis) CALL get_qs_kind(qs_kind, basis_set=lri_aux_basis, basis_type="LRI_AUX") CPASSERT(ASSOCIATED(lri_aux_basis)) + CALL get_gto_basis_set(gto_basis_set=lri_aux_basis, kind_radius=radb) + rmax = MAX(rmax, radb) END DO + rmax = 1.25_dp*rmax lri_section => section_vals_get_subs_vals(input, "DFT%KG_METHOD%LRIGPW") CALL lri_env_init(kg_env%lri_env, lri_section) CALL lri_env_basis("LRI", qs_env, kg_env%lri_env, qs_kind_set) + ! + ! integration grid + ! + CALL section_vals_val_get(input, "DFT%KG_METHOD%INTEGRATION_GRID", c_val=intgrid) + CALL uppercase(intgrid) + SELECT CASE (intgrid) + CASE ("SMALL") + nr = 50 + na = 38 + CASE ("MEDIUM") + nr = 100 + na = 110 + CASE ("LARGE") + nr = 200 + na = 302 + CASE ("HUGE") + nr = 400 + na = 590 + CASE DEFAULT + CPABORT("KG:INTEGRATION_GRID") + END SELECT + NULLIFY (logger) + logger => cp_get_default_logger() + iunit = cp_logger_get_default_io_unit(logger) + CALL initialize_atomic_grid(kg_env%int_grid_atom, nr, na, rmax, iunit=iunit) + ! load balancing + CALL get_qs_env(qs_env=qs_env, natom=natom, para_env=para_env) + np = para_env%num_pe + load = REAL(natom, KIND=dp)*kg_env%int_grid_atom%ntot/REAL(np, KIND=dp) + ! + CALL allocate_intgrid(kg_env%int_grid_full) + ig_full => kg_env%int_grid_full + CALL allocate_intgrid(kg_env%int_grid_molecules) + ig_mol => kg_env%int_grid_molecules + nbatch = (natom*kg_env%int_grid_atom%nbatch)/np + nbatch = NINT((nbatch+1)*1.2_dp) + ALLOCATE (bid(2, nbatch)) + nbatch = 0 + DO iatom = 1, natom + DO ib = 1, kg_env%int_grid_atom%nbatch + IF (para_env%mepos == MOD(iatom+ib, np)) THEN + nbatch = nbatch+1 + CPASSERT(nbatch <= SIZE(bid, 2)) + bid(1, nbatch) = iatom + bid(2, nbatch) = ib + END IF + END DO + END DO + ! + ig_full%nbatch = nbatch + ALLOCATE (ig_full%grid_batch(nbatch)) + ! + ig_mol%nbatch = nbatch + ALLOCATE (ig_mol%grid_batch(nbatch)) + ! + DO i = 1, nbatch + iatom = bid(1, i) + ib = bid(2, i) + ! + ig_full%grid_batch(i)%ref_atom = iatom + ig_full%grid_batch(i)%ibatch = ib + ig_full%grid_batch(i)%np = kg_env%int_grid_atom%batch(ib)%np + ig_full%grid_batch(i)%radius = kg_env%int_grid_atom%batch(ib)%rad + ig_full%grid_batch(i)%rcenter(1:3) = kg_env%int_grid_atom%batch(ib)%rcenter(1:3) + n = ig_full%grid_batch(i)%np + ALLOCATE (ig_full%grid_batch(i)%rco(3, n)) + ALLOCATE (ig_full%grid_batch(i)%weight(n)) + ALLOCATE (ig_full%grid_batch(i)%wref(n)) + ALLOCATE (ig_full%grid_batch(i)%wsum(n)) + ig_full%grid_batch(i)%weight(:) = kg_env%int_grid_atom%batch(ib)%weight(:) + ! + ig_mol%grid_batch(i)%ref_atom = iatom + ig_mol%grid_batch(i)%ibatch = ib + ig_mol%grid_batch(i)%np = kg_env%int_grid_atom%batch(ib)%np + ig_mol%grid_batch(i)%radius = kg_env%int_grid_atom%batch(ib)%rad + ig_mol%grid_batch(i)%rcenter(1:3) = kg_env%int_grid_atom%batch(ib)%rcenter(1:3) + n = ig_mol%grid_batch(i)%np + ALLOCATE (ig_mol%grid_batch(i)%rco(3, n)) + ALLOCATE (ig_mol%grid_batch(i)%weight(n)) + ALLOCATE (ig_mol%grid_batch(i)%wref(n)) + ALLOCATE (ig_mol%grid_batch(i)%wsum(n)) + ig_mol%grid_batch(i)%weight(:) = kg_env%int_grid_atom%batch(ib)%weight(:) + END DO + ! + DEALLOCATE (bid) END IF CALL timestop(handle) diff --git a/src/kg_environment_types.F b/src/kg_environment_types.F index 1f216fe623..27879a1c95 100644 --- a/src/kg_environment_types.F +++ b/src/kg_environment_types.F @@ -14,6 +14,8 @@ MODULE kg_environment_types USE cp_dbcsr_operations, ONLY: dbcsr_deallocate_matrix_set USE dbcsr_api, ONLY: dbcsr_p_type USE input_section_types, ONLY: section_vals_type + USE integration_grid_types, ONLY: deallocate_intgrid,& + integration_grid_type USE kinds, ONLY: dp USE lri_environment_types, ONLY: lri_density_release,& lri_density_type,& @@ -22,6 +24,8 @@ MODULE kg_environment_types USE molecule_types, ONLY: molecule_type USE qs_dispersion_types, ONLY: qs_dispersion_release,& qs_dispersion_type + USE qs_grid_atom, ONLY: atom_integration_grid_type,& + deallocate_atom_int_grid USE qs_neighbor_list_types, ONLY: deallocate_neighbor_list_set,& neighbor_list_set_p_type USE task_list_types, ONLY: deallocate_task_list,& @@ -100,20 +104,27 @@ MODULE kg_environment_types TYPE(section_vals_type), POINTER :: xc_section_kg INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_to_molecule TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set - TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: tnadd_mat + INTEGER :: tnadd_method TYPE(neighbor_list_set_p_type), & DIMENSION(:), POINTER :: sab_orb_full, sac_kin + ! INTEGER, DIMENSION(:), POINTER :: subset_of_mol TYPE(subset_type), DIMENSION(:), POINTER :: subset INTEGER :: nsubsets INTEGER :: maxdegree INTEGER :: coloring_method - INTEGER :: tnadd_method + ! + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: tnadd_mat + ! LOGICAL :: energy_correction TYPE(energy_correction_type) :: ec_env ! LRI TYPE(lri_environment_type), POINTER :: lri_env TYPE(lri_density_type), POINTER :: lri_density + ! atomic grid + TYPE(atom_integration_grid_type), POINTER :: int_grid_atom + TYPE(integration_grid_type), POINTER :: int_grid_molecules + TYPE(integration_grid_type), POINTER :: int_grid_full END TYPE kg_environment_type CONTAINS @@ -170,6 +181,16 @@ CONTAINS IF (ASSOCIATED(kg_env%lri_density)) THEN CALL lri_density_release(kg_env%lri_density) END IF + ! atom grids + IF (ASSOCIATED(kg_env%int_grid_atom)) THEN + CALL deallocate_atom_int_grid(kg_env%int_grid_atom) + END IF + IF (ASSOCIATED(kg_env%int_grid_molecules)) THEN + CALL deallocate_intgrid(kg_env%int_grid_molecules) + END IF + IF (ASSOCIATED(kg_env%int_grid_full)) THEN + CALL deallocate_intgrid(kg_env%int_grid_full) + END IF ! energy correction IF (kg_env%energy_correction) THEN diff --git a/src/qs_energy_init.F b/src/qs_energy_init.F index 9911b74da0..a5ea7633ec 100644 --- a/src/qs_energy_init.F +++ b/src/qs_energy_init.F @@ -185,7 +185,7 @@ CONTAINS END IF SELECT CASE (kg_env%tnadd_method) - CASE (kg_tnadd_embed, kg_tnadd_embed_ri) + CASE (kg_tnadd_embed) ! allocate the subset list IF (.NOT. ASSOCIATED(kg_env%subset_of_mol)) THEN ALLOCATE (kg_env%subset_of_mol(SIZE(molecule_set))) @@ -198,6 +198,27 @@ CONTAINS CALL kg_build_neighborlist(qs_env, sab_orb=kg_env%subset(isubset)%sab_orb, molecular=.TRUE., & subset_of_mol=kg_env%subset_of_mol, current_subset=isubset) END DO + CASE (kg_tnadd_embed_ri) +!deb should be deleted as soon as atomic grids work + ! allocate the subset list + IF (.NOT. ASSOCIATED(kg_env%subset_of_mol)) THEN + ALLOCATE (kg_env%subset_of_mol(SIZE(molecule_set))) + END IF + ! + CALL kg_build_subsets(kg_env, para_env) + ! + DO isubset = 1, kg_env%nsubsets + ! build the (new) molecular neighborlist of the current subset + CALL kg_build_neighborlist(qs_env, sab_orb=kg_env%subset(isubset)%sab_orb, molecular=.TRUE., & + subset_of_mol=kg_env%subset_of_mol, current_subset=isubset) + END DO +!deb + ! LRI neighborlist + NULLIFY (soo_list) + CALL kg_build_neighborlist(qs_env, sab_orb=soo_list, molecular=.TRUE.) + kg_env%lri_env%soo_list => soo_list + CALL calculate_lri_integrals(kg_env%lri_env, qs_env) + ! Atomic grids CASE (kg_tnadd_atomic) ! build the A-C list for the nonadditive kinetic energy potential CALL kg_build_neighborlist(qs_env, sac_kin=kg_env%sac_kin) @@ -207,14 +228,6 @@ CONTAINS CPABORT("KG:TNADD METHOD") END SELECT - ! LRI neighborlist - IF (kg_env%tnadd_method == kg_tnadd_embed_ri) THEN - NULLIFY (soo_list) - CALL kg_build_neighborlist(qs_env, sab_orb=soo_list, molecular=.TRUE.) - kg_env%lri_env%soo_list => soo_list - CALL calculate_lri_integrals(kg_env%lri_env, qs_env) - END IF - CALL timestop(handle) END SUBROUTINE qs_energies_init_kg diff --git a/src/qs_environment.F b/src/qs_environment.F index 88fca1cd2e..0947a6de3c 100644 --- a/src/qs_environment.F +++ b/src/qs_environment.F @@ -471,8 +471,8 @@ CONTAINS nelectron, nkind, output_unit, tnadd_method INTEGER, DIMENSION(2) :: n_mo, nelectron_spin LOGICAL :: all_potential_present, be_silent, & - do_kpoints, has_unit_metric, lribas, & - was_present + do_kpoints, e1terms, has_unit_metric, & + lribas, was_present REAL(dp) :: ewald_rcut, maxocc, verlet_skin TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(cp_logger_type), POINTER :: logger @@ -672,14 +672,19 @@ CONTAINS END IF lribas = .FALSE. - IF (dft_control%qs_control%method_id == do_method_lrigpw) lribas = .TRUE. + e1terms = .FALSE. + IF (dft_control%qs_control%method_id == do_method_lrigpw) THEN + lribas = .TRUE. + CALL get_qs_env(qs_env, lri_env=lri_env) + e1terms = lri_env%exact_1c_terms + END IF IF (dft_control%qs_control%do_kg) THEN CALL section_vals_val_get(dft_section, "KG_METHOD%TNADD_METHOD", i_val=tnadd_method) IF (tnadd_method == kg_tnadd_embed_ri) lribas = .TRUE. END IF IF (lribas) THEN ! Check if LRI_AUX basis is available, auto-generate if needed - CALL get_qs_env(qs_env, nkind=nkind, lri_env=lri_env) + CALL get_qs_env(qs_env, nkind=nkind) DO ikind = 1, nkind NULLIFY (lri_aux_basis) qs_kind => qs_kind_set(ikind) @@ -689,8 +694,7 @@ CONTAINS CALL cp_warn(__LOCATION__, "Automatic Generation of LRI_AUX basis. "// & "This is experimental code.") ! Generate a default basis - CALL create_lri_aux_basis_set(lri_aux_basis, qs_kind, & - dft_control%auto_basis_lri_aux, lri_env%exact_1c_terms) + CALL create_lri_aux_basis_set(lri_aux_basis, qs_kind, dft_control%auto_basis_lri_aux, e1terms) CALL add_basis_set_to_container(qs_kind%basis_sets, lri_aux_basis, "LRI_AUX") END IF END DO diff --git a/src/qs_grid_atom.F b/src/qs_grid_atom.F index e4cc6305fb..893007248a 100644 --- a/src/qs_grid_atom.F +++ b/src/qs_grid_atom.F @@ -9,7 +9,8 @@ MODULE qs_grid_atom do_gapw_gct,& do_gapw_log USE kinds, ONLY: dp - USE lebedev, ONLY: lebedev_grid + USE lebedev, ONLY: get_number_of_lebedev_grid,& + lebedev_grid USE mathconstants, ONLY: pi USE memory_utilities, ONLY: reallocate #include "./base/base_uses.f90" @@ -20,8 +21,25 @@ MODULE qs_grid_atom CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_grid_atom' - TYPE grid_atom_type + TYPE grid_batch_type + INTEGER :: np + REAL(KIND=dp), DIMENSION(3) :: rcenter + REAL(KIND=dp) :: rad + REAL(dp), DIMENSION(:, :), ALLOCATABLE :: rco + REAL(dp), DIMENSION(:), ALLOCATABLE :: weight + END TYPE grid_batch_type + TYPE atom_integration_grid_type + INTEGER :: nr, na + INTEGER :: np, ntot + INTEGER :: lebedev_grid + REAL(dp), DIMENSION(:), ALLOCATABLE :: rr + REAL(dp), DIMENSION(:), ALLOCATABLE :: wr, wa + INTEGER :: nbatch + TYPE(grid_batch_type), DIMENSION(:), ALLOCATABLE :: batch + END TYPE atom_integration_grid_type + + TYPE grid_atom_type INTEGER :: nr, ng_sphere REAL(dp), DIMENSION(:), POINTER :: rad, rad2, & wr, wa, & @@ -29,19 +47,19 @@ MODULE qs_grid_atom pol, cos_pol, sin_pol, usin_azi REAL(dp), DIMENSION(:, :), & POINTER :: rad2l, oorad2l, weight - END TYPE grid_atom_type - PUBLIC :: allocate_grid_atom, & - create_grid_atom, & - deallocate_grid_atom - + PUBLIC :: allocate_grid_atom, create_grid_atom, deallocate_grid_atom PUBLIC :: grid_atom_type + PUBLIC :: initialize_atomic_grid + PUBLIC :: atom_integration_grid_type, deallocate_atom_int_grid + +! ************************************************************************************************** CONTAINS ! ************************************************************************************************** -!> \brief Deallocate a Gaussian-type orbital (GTO) basis set data set. +!> \brief Initialize components of the grid_atom_type structure !> \param grid_atom ... !> \date 03.11.2000 !> \author MK @@ -257,6 +275,317 @@ CONTAINS END SUBROUTINE create_grid_atom ! ************************************************************************************************** +!> \brief Initialize atomic grid +!> \param int_grid ... +!> \param nr ... +!> \param na ... +!> \param rmax ... +!> \param quadrature ... +!> \param iunit ... +!> \date 02.2018 +!> \author JGH +!> \version 1.0 +! ************************************************************************************************** + SUBROUTINE initialize_atomic_grid(int_grid, nr, na, rmax, quadrature, iunit) + TYPE(atom_integration_grid_type), POINTER :: int_grid + INTEGER, INTENT(IN) :: nr, na + REAL(KIND=dp), INTENT(IN) :: rmax + INTEGER, INTENT(IN), OPTIONAL :: quadrature, iunit + + CHARACTER(len=*), PARAMETER :: routineN = 'initialize_atomic_grid', & + routineP = moduleN//':'//routineN + + INTEGER :: ia, ig, ir, ix, iy, iz, la, ll, my_quad, & + n1, n2, n3, nbatch, ng, no, np, ntot, & + nu, nx + INTEGER, ALLOCATABLE, DIMENSION(:) :: icell + REAL(KIND=dp) :: ag, dd, dmax, r1, r2, r3 + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: rad, rad2, wa, wc, wr + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rang, rco + REAL(KIND=dp), DIMENSION(10) :: dco + REAL(KIND=dp), DIMENSION(3) :: rm + TYPE(atom_integration_grid_type), POINTER :: igr + + ALLOCATE (igr) + + ! type of quadrature grid + IF (PRESENT(quadrature)) THEN + my_quad = quadrature + ELSE + my_quad = do_gapw_log + END IF + + ! radial grid + CPASSERT(nr > 1) + ALLOCATE (rad(nr), rad2(nr), wr(nr)) + CALL radial_grid(nr, rad, rad2, wr, my_quad) + ! + igr%nr = nr + ALLOCATE (igr%rr(nr)) + ALLOCATE (igr%wr(nr)) + ! store grid points always in ascending order + IF (rad(1) > rad(nr)) THEN + DO ir = nr, 1, -1 + igr%rr(nr-ir+1) = rad(ir) + igr%wr(nr-ir+1) = wr(ir) + END DO + ELSE + igr%rr(1:nr) = rad(1:nr) + igr%wr(1:nr) = wr(1:nr) + END IF + ! only include grid points smaller than rmax + np = 0 + DO ir = 1, nr + IF (igr%rr(ir) < rmax) THEN + np = np+1 + rad(np) = igr%rr(ir) + wr(np) = igr%wr(ir) + END IF + END DO + igr%np = np + ! + ! angular grid + CPASSERT(na > 1) + ll = get_number_of_lebedev_grid(n=na) + np = lebedev_grid(ll)%n + la = lebedev_grid(ll)%l + ALLOCATE (rang(3, np), wa(np)) + wa(1:na) = 4._dp*pi*lebedev_grid(ll)%w(1:np) + rang(1:3, 1:np) = lebedev_grid(ll)%r(1:3, 1:np) + igr%lebedev_grid = ll + ALLOCATE (igr%wa(np)) + igr%na = np + igr%wa(1:np) = wa(1:np) + ! + ! total grid points + ntot = igr%na*igr%np + igr%ntot = ntot + ALLOCATE (rco(3, ntot), wc(ntot)) + ig = 0 + DO ir = 1, igr%np + DO ia = 1, igr%na + ig = ig+1 + rco(1:3, ig) = rang(1:3, ia)*rad(ir) + wc(ig) = wa(ia)*wr(ir) + END DO + END DO + ! grid for batches, odd number of cells + ng = NINT((REAL(ntot, dp)/32._dp)**0.33333_dp) + ng = ng+MOD(ng+1, 2) + ! avarage number of points along radial grid + dco = 0.0_dp + ag = REAL(igr%np, dp)/ng + CPASSERT(SIZE(dco) >= (ng+1)/2) + DO ig = 1, ng, 2 + ir = MIN(NINT(ag)*ig, igr%np) + ia = (ig+1)/2 + dco(ia) = rad(ir) + END DO + ! batches + ALLOCATE (icell(ntot)) + icell = 0 + nx = (ng-1)/2 + DO ig = 1, ntot + ix = grid_coord(rco(1, ig), dco, nx+1)+nx + iy = grid_coord(rco(2, ig), dco, nx+1)+nx + iz = grid_coord(rco(3, ig), dco, nx+1)+nx + icell(ig) = iz*ng*ng+iy*ng+ix+1 + END DO + ! + igr%nbatch = ng*ng*ng + ALLOCATE (igr%batch(igr%nbatch)) + igr%batch(:)%np = 0 + DO ig = 1, ntot + ia = icell(ig) + igr%batch(ia)%np = igr%batch(ia)%np+1 + END DO + DO ig = 1, igr%nbatch + np = igr%batch(ig)%np + ALLOCATE (igr%batch(ig)%rco(3, np), igr%batch(ig)%weight(np)) + igr%batch(ig)%np = 0 + END DO + DO ig = 1, ntot + ia = icell(ig) + igr%batch(ia)%np = igr%batch(ia)%np+1 + np = igr%batch(ia)%np + igr%batch(ia)%rco(1:3, np) = rco(1:3, ig) + igr%batch(ia)%weight(np) = wc(ig) + END DO + ! + DEALLOCATE (rad, rad2, rang, wr, wa) + DEALLOCATE (rco, wc, icell) + ! + IF (ASSOCIATED(int_grid)) CALL deallocate_atom_int_grid(int_grid) + ALLOCATE (int_grid) + ALLOCATE (int_grid%rr(igr%nr), int_grid%wr(igr%nr), int_grid%wa(igr%na)) + int_grid%nr = igr%nr + int_grid%na = igr%na + int_grid%np = igr%np + int_grid%ntot = igr%ntot + int_grid%lebedev_grid = igr%lebedev_grid + int_grid%rr(:) = igr%rr(:) + int_grid%wr(:) = igr%wr(:) + int_grid%wa(:) = igr%wa(:) + ! + ! count batches + nbatch = 0 + DO ig = 1, igr%nbatch + IF (igr%batch(ig)%np == 0) THEN + ! empty batch + ELSE IF (igr%batch(ig)%np <= 48) THEN + ! single batch + nbatch = nbatch+1 + ELSE + ! multiple batches + nbatch = nbatch+NINT(igr%batch(ig)%np/32._dp) + END IF + END DO + int_grid%nbatch = nbatch + ALLOCATE (int_grid%batch(nbatch)) + ! fill batches + n1 = 0 + DO ig = 1, igr%nbatch + IF (igr%batch(ig)%np == 0) THEN + ! empty batch + ELSE IF (igr%batch(ig)%np <= 48) THEN + ! single batch + n1 = n1+1 + np = igr%batch(ig)%np + ALLOCATE (int_grid%batch(n1)%rco(3, np), int_grid%batch(n1)%weight(np)) + int_grid%batch(n1)%np = np + int_grid%batch(n1)%rco(1:3, 1:np) = igr%batch(ig)%rco(1:3, 1:np) + int_grid%batch(n1)%weight(1:np) = igr%batch(ig)%weight(1:np) + ELSE + ! multiple batches + n2 = NINT(igr%batch(ig)%np/32._dp) + n3 = igr%batch(ig)%np/n2 + DO ia = n1+1, n1+n2 + nu = (ia-n1-1)*n3+1 + no = nu+n3-1 + IF (ia == n1+n2) no = igr%batch(ig)%np + np = no-nu+1 + ALLOCATE (int_grid%batch(ia)%rco(3, np), int_grid%batch(ia)%weight(np)) + int_grid%batch(ia)%np = np + int_grid%batch(ia)%rco(1:3, 1:np) = igr%batch(ig)%rco(1:3, nu:no) + int_grid%batch(ia)%weight(1:np) = igr%batch(ig)%weight(nu:no) + END DO + n1 = n1+n2 + END IF + END DO + CPASSERT(nbatch == n1) + ! batch center and radius + DO ig = 1, int_grid%nbatch + np = int_grid%batch(ig)%np + IF (np > 0) THEN + rm(1) = SUM(int_grid%batch(ig)%rco(1, 1:np)) + rm(2) = SUM(int_grid%batch(ig)%rco(2, 1:np)) + rm(3) = SUM(int_grid%batch(ig)%rco(3, 1:np)) + rm(1:3) = rm(1:3)/REAL(np, KIND=dp) + END IF + int_grid%batch(ig)%rcenter(1:3) = rm(1:3) + dmax = 0.0_dp + DO ia = 1, np + dd = SUM((int_grid%batch(ig)%rco(1:3, ia)-rm(1:3))**2) + dmax = MAX(dd, dmax) + END DO + int_grid%batch(ig)%rad = SQRT(dmax) + END DO + ! + CALL deallocate_atom_int_grid(igr) + ! + IF (PRESENT(iunit)) THEN + IF (iunit > 0) THEN + WRITE (iunit, "(/,A)") " Atomic Integration Grid Information" + WRITE (iunit, "(A,T51,3I10)") " Number of grid points [radial,angular,total]", & + int_grid%np, int_grid%na, int_grid%ntot + WRITE (iunit, "(A,T71,I10)") " Lebedev grid number", int_grid%lebedev_grid + WRITE (iunit, "(A,T61,F20.5)") " Maximum of radial grid [Bohr]", & + int_grid%rr(int_grid%np) + nbatch = int_grid%nbatch + WRITE (iunit, "(A,T71,I10)") " Total number of gridpoint batches", nbatch + n1 = int_grid%ntot + n2 = 0 + n3 = NINT(REAL(int_grid%ntot, dp)/REAL(nbatch, dp)) + DO ig = 1, nbatch + n1 = MIN(n1, int_grid%batch(ig)%np) + n2 = MAX(n2, int_grid%batch(ig)%np) + END DO + WRITE (iunit, "(A,T51,3I10)") " Number of grid points/batch [min,max,ave]", n1, n2, n3 + r1 = 1000._dp + r2 = 0.0_dp + r3 = 0.0_dp + DO ig = 1, int_grid%nbatch + r1 = MIN(r1, int_grid%batch(ig)%rad) + r2 = MAX(r2, int_grid%batch(ig)%rad) + r3 = r3+int_grid%batch(ig)%rad + END DO + r3 = r3/REAL(ng*ng*ng, KIND=dp) + WRITE (iunit, "(A,T51,3F10.2)") " Batch radius (bohr) [min,max,ave]", r1, r2, r3 + END IF + END IF + + END SUBROUTINE initialize_atomic_grid + +! ************************************************************************************************** +!> \brief ... +!> \param x ... +!> \param dco ... +!> \param ng ... +!> \return ... +!> \retval igrid ... +! ************************************************************************************************** + FUNCTION grid_coord(x, dco, ng) RESULT(igrid) + REAL(KIND=dp), INTENT(IN) :: x + REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: dco + INTEGER, INTENT(IN) :: ng + INTEGER :: igrid + + INTEGER :: ig + REAL(KIND=dp) :: xval + + xval = ABS(x) + igrid = ng + DO ig = 1, ng + IF (xval <= dco(ig)) THEN + igrid = ig-1 + EXIT + END IF + END DO + IF (x < 0.0_dp) igrid = -igrid + CPASSERT(ABS(igrid) < ng) + END FUNCTION grid_coord + +! ************************************************************************************************** +!> \brief ... +!> \param int_grid ... +! ************************************************************************************************** + SUBROUTINE deallocate_atom_int_grid(int_grid) + TYPE(atom_integration_grid_type), POINTER :: int_grid + + CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_atom_int_grid', & + routineP = moduleN//':'//routineN + + INTEGER :: ib + + IF (ASSOCIATED(int_grid)) THEN + IF (ALLOCATED(int_grid%rr)) DEALLOCATE (int_grid%rr) + IF (ALLOCATED(int_grid%wr)) DEALLOCATE (int_grid%wr) + IF (ALLOCATED(int_grid%wa)) DEALLOCATE (int_grid%wa) + ! batch + IF (ALLOCATED(int_grid%batch)) THEN + DO ib = 1, SIZE(int_grid%batch) + IF (ALLOCATED(int_grid%batch(ib)%rco)) DEALLOCATE (int_grid%batch(ib)%rco) + IF (ALLOCATED(int_grid%batch(ib)%weight)) DEALLOCATE (int_grid%batch(ib)%weight) + END DO + DEALLOCATE (int_grid%batch) + END IF + ! + DEALLOCATE (int_grid) + NULLIFY (int_grid) + END IF + + END SUBROUTINE deallocate_atom_int_grid +! ************************************************************************************************** !> \brief Generate a radial grid with n points by a quadrature rule. !> \param n ... !> \param r ... @@ -275,7 +604,7 @@ CONTAINS SUBROUTINE radial_grid(n, r, r2, wr, radial_quadrature) INTEGER, INTENT(IN) :: n - REAL(dp), DIMENSION(:), POINTER :: r, r2, wr + REAL(dp), DIMENSION(:), INTENT(INOUT) :: r, r2, wr INTEGER, INTENT(IN) :: radial_quadrature CHARACTER(len=*), PARAMETER :: routineN = 'radial_grid', routineP = moduleN//':'//routineN diff --git a/src/xc/xc_libxc.F b/src/xc/xc_libxc.F index 9633dca4be..499816a60e 100644 --- a/src/xc/xc_libxc.F +++ b/src/xc/xc_libxc.F @@ -1336,7 +1336,7 @@ CONTAINS DO ii = 1, npoints IF ((rho(ii) > epsilon_rho) .AND. (tau(ii) > epsilon_tau)) THEN sigma = norm_drho(ii)**2 - my_tau = MAX(tau(ii), sigma(1)/(8.0_dp*rho(ii))) + my_tau(1) = MAX(tau(ii), sigma(1)/(8.0_dp*rho(ii))) CALL xc_f03_mgga_exc(xc_func, 1, rho(ii), sigma, & laplace_rho(ii), my_tau, exc) e_0(ii) = e_0(ii)+sc*exc(1)*rho(ii) @@ -1348,7 +1348,7 @@ CONTAINS DO ii = 1, npoints IF ((rho(ii) > epsilon_rho) .AND. (tau(ii) > epsilon_tau)) THEN sigma = norm_drho(ii)**2 - my_tau = MAX(tau(ii), sigma(1)/(8.0_dp*rho(ii))) + my_tau(1) = MAX(tau(ii), sigma(1)/(8.0_dp*rho(ii))) CALL xc_f03_mgga_vxc(xc_func, 1, rho(ii), sigma, & laplace_rho(ii), my_tau, vrho, vsigma, vlapl, vtau) e_rho(ii) = e_rho(ii)+sc*vrho(1) @@ -1362,8 +1362,8 @@ CONTAINS !$OMP DO DO ii = 1, npoints IF ((rho(ii) > epsilon_rho) .AND. (tau(ii) > epsilon_tau)) THEN - sigma = norm_drho(ii)**2 - my_tau = MAX(tau(ii), sigma(1)/(8.0_dp*rho(ii))) + sigma(1) = norm_drho(ii)**2 + my_tau(1) = MAX(tau(ii), sigma(1)/(8.0_dp*rho(ii))) IF (no_exc) THEN CALL xc_f03_mgga_vxc(xc_func, 1, rho(ii), sigma, & laplace_rho(ii), my_tau, vrho, vsigma, vlapl, vtau) @@ -1385,7 +1385,7 @@ CONTAINS DO ii = 1, npoints IF ((rho(ii) > epsilon_rho) .AND. (tau(ii) > epsilon_tau)) THEN sigma = norm_drho(ii)**2 - my_tau = MAX(tau(ii), sigma(1)/(8.0_dp*rho(ii))) + my_tau(1) = MAX(tau(ii), sigma(1)/(8.0_dp*rho(ii))) IF (no_exc) THEN CALL xc_f03_mgga_vxc(xc_func, 1, rho(ii), sigma, & laplace_rho(ii), my_tau, vrho, vsigma, vlapl, vtau) @@ -1419,7 +1419,7 @@ CONTAINS DO ii = 1, npoints IF ((rho(ii) > epsilon_rho) .AND. (tau(ii) > epsilon_tau)) THEN sigma = norm_drho(ii)**2 - my_tau = MAX(tau(ii), sigma(1)/(8.0_dp*rho(ii))) + my_tau(1) = MAX(tau(ii), sigma(1)/(8.0_dp*rho(ii))) IF (no_exc) THEN CALL xc_f03_mgga_vxc(xc_func, 1, rho(ii), sigma, & laplace_rho(ii), my_tau, vrho, vsigma, vlapl, vtau)