comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 22:25:32 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

KG: towards atomic integration grid

Browse source 
svn-origin-rev: 18614
This commit is contained in:
Jürg Hutter 2018-09-25 07:02:03 +00:00
parent 96e300a19e
commit c2221170b8
9 changed files with 697 additions and 35 deletions
Download patch file Download diff file Expand all files Collapse all files

7
src/input_cp2k_dft.F
View file

@ -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.)

179
src/integration_grid_types.F Normal file
View file

@ -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

6
src/kg_correction.F
View file

@ -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

109
src/kg_environment.F
View file

@ -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)

25
src/kg_environment_types.F
View file

@ -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

31
src/qs_energy_init.F
View file

@ -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

16
src/qs_environment.F
View file

@ -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

347
src/qs_grid_atom.F
View file

@ -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

12
src/xc/xc_libxc.F
View file

@ -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)