Skip to content

Commit

Permalink
Harris Method (cp2k#3665)
Browse files Browse the repository at this point in the history
  • Loading branch information
@juerghutter
juerghutter authored Sep 9, 2024
1 parent 17fb60b commit f2bc97d
Show file tree
Hide file tree
Showing 29 changed files with 1,841 additions and 144 deletions.
3 changes: 3 additions & 0 deletions src/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -246,6 +246,7 @@ list(
input_cp2k_colvar.F
input_cp2k_constraints.F
input_cp2k_dft.F
input_cp2k_harris.F
input_cp2k_print_dft.F
input_cp2k_scf.F
input_cp2k_as.F
Expand Down Expand Up @@ -575,6 +576,8 @@ list(
qs_gcp_utils.F
qs_grid_atom.F
qs_gspace_mixing.F
qs_harris_types.F
qs_harris_utils.F
qs_harmonics_atom.F
qs_hash_table_functions.F
qs_initial_guess.F
Expand Down
1 change: 0 additions & 1 deletion src/aobasis/aux_basis_set.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -220,7 +220,6 @@ SUBROUTINE create_aux_basis(aux_basis, bsname, nsets, lmin, lmax, nl, npgf, zet)
ALLOCATE (aux_basis%scon(maxco, nsgf))
ALLOCATE (aux_basis%norm_cgf(ncgf))
aux_basis%norm_type = 2
! CALL init_orb_basis_set(aux_basis)

END SUBROUTINE create_aux_basis

Expand Down
5 changes: 4 additions & 1 deletion src/aobasis/basis_set_container_types.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,8 @@ MODULE basis_set_container_types
p_lri_aux_basis = 116, &
aux_opt_basis = 117, &
min_basis = 118, &
tda_k_basis = 119
tda_k_basis = 119, &
rhoin_basis = 120
! **************************************************************************************************
TYPE basis_set_container_type
PRIVATE
Expand Down Expand Up @@ -133,6 +134,8 @@ FUNCTION get_basis_type(basis_set_type) RESULT(basis_type_nr)
basis_type_nr = ri_xas_basis
CASE ("AUX_OPT")
basis_type_nr = aux_opt_basis
CASE ("RHOIN")
basis_type_nr = rhoin_basis
CASE DEFAULT
basis_type_nr = unknown_basis
END SELECT
Expand Down
9 changes: 8 additions & 1 deletion src/aobasis/basis_set_types.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -279,13 +279,15 @@ END SUBROUTINE copy_gto_basis_set
!> \brief ...
!> \param basis_set ...
!> \param pbasis ...
!> \param lmax ...
! **************************************************************************************************
SUBROUTINE create_primitive_basis_set(basis_set, pbasis)
SUBROUTINE create_primitive_basis_set(basis_set, pbasis, lmax)

! Create a primitives only basis set

TYPE(gto_basis_set_type), INTENT(IN) :: basis_set
TYPE(gto_basis_set_type), POINTER :: pbasis
INTEGER, INTENT(IN), OPTIONAL :: lmax

INTEGER :: i, ico, ip, ipgf, iset, ishell, l, lm, &
lshell, m, maxco, mpgf, nc, ncgf, ns, &
Expand Down Expand Up @@ -332,6 +334,11 @@ SUBROUTINE create_primitive_basis_set(basis_set, pbasis)

CALL allocate_gto_basis_set(pbasis)

IF (PRESENT(lmax)) THEN
! if requested, reduce max l val
lm = MIN(lm, lmax)
END IF

pbasis%name = basis_set%name//"_primitive"
pbasis%kind_radius = basis_set%kind_radius
pbasis%short_kind_radius = basis_set%short_kind_radius
Expand Down
135 changes: 74 additions & 61 deletions src/atom_fit.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,7 @@ MODULE atom_fit
evolt
USE powell, ONLY: opt_state_type,&
powell_optimize
USE qs_grid_atom, ONLY: grid_atom_type
#include "./base/base_uses.f90"

IMPLICIT NONE
Expand All @@ -69,24 +70,27 @@ MODULE atom_fit
!> \param num_gto number of Gaussian basis functions
!> \param norder ...
!> \param iunit output file unit
!> \param agto Gaussian exponents
!> \param powell_section POWELL input section
!> \param results (output) array of num_gto+2 real numbers in the following format:
!> starting exponent, factor of geometrical series, expansion coefficients(1:num_gto)
! **************************************************************************************************
SUBROUTINE atom_fit_density(atom, num_gto, norder, iunit, powell_section, results)
SUBROUTINE atom_fit_density(atom, num_gto, norder, iunit, agto, powell_section, results)
TYPE(atom_type), POINTER :: atom
INTEGER, INTENT(IN) :: num_gto, norder, iunit
REAL(KIND=dp), DIMENSION(:), OPTIONAL :: agto
TYPE(section_vals_type), OPTIONAL, POINTER :: powell_section
REAL(KIND=dp), DIMENSION(:), OPTIONAL :: results

INTEGER :: n10
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: co
INTEGER :: i, n10
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: co, pe
REAL(KIND=dp), DIMENSION(2) :: x
TYPE(opgrid_type), POINTER :: density
TYPE(opt_state_type) :: ostate

ALLOCATE (co(num_gto))
ALLOCATE (co(num_gto), pe(num_gto))
co = 0._dp
pe = 0._dp
NULLIFY (density)
CALL create_opgrid(density, atom%basis%grid)
CALL atom_denmat(atom%orbitals%pmat, atom%orbitals%wfn, atom%basis%nbas, atom%state%occupation, &
Expand All @@ -98,52 +102,60 @@ SUBROUTINE atom_fit_density(atom, num_gto, norder, iunit, powell_section, result
density%op = density%op*atom%basis%grid%rad**norder
END IF

ostate%nf = 0
ostate%nvar = 2
x(1) = 0.10_dp !starting point of geometric series
x(2) = 2.00_dp !factor of series
IF (PRESENT(powell_section)) THEN
CALL section_vals_val_get(powell_section, "ACCURACY", r_val=ostate%rhoend)
CALL section_vals_val_get(powell_section, "STEP_SIZE", r_val=ostate%rhobeg)
CALL section_vals_val_get(powell_section, "MAX_FUN", i_val=ostate%maxfun)
IF (PRESENT(agto)) THEN
CPASSERT(num_gto <= SIZE(agto))
pe(1:num_gto) = agto(1:num_gto)
ELSE
ostate%rhoend = 1.e-8_dp
ostate%rhobeg = 5.e-2_dp
ostate%maxfun = 1000
END IF
ostate%iprint = 1
ostate%unit = iunit
ostate%nf = 0
ostate%nvar = 2
x(1) = 0.10_dp !starting point of geometric series
x(2) = 2.00_dp !factor of series
IF (PRESENT(powell_section)) THEN
CALL section_vals_val_get(powell_section, "ACCURACY", r_val=ostate%rhoend)
CALL section_vals_val_get(powell_section, "STEP_SIZE", r_val=ostate%rhobeg)
CALL section_vals_val_get(powell_section, "MAX_FUN", i_val=ostate%maxfun)
ELSE
ostate%rhoend = 1.e-8_dp
ostate%rhobeg = 5.e-2_dp
ostate%maxfun = 1000
END IF
ostate%iprint = 1
ostate%unit = iunit

ostate%state = 0
IF (iunit > 0) THEN
WRITE (iunit, '(/," POWELL| Start optimization procedure")')
END IF
n10 = ostate%maxfun/10
ostate%state = 0
IF (iunit > 0) THEN
WRITE (iunit, '(/," POWELL| Start optimization procedure")')
END IF
n10 = ostate%maxfun/10

DO
DO

IF (ostate%state == 2) THEN
CALL density_fit(density, atom, num_gto, x(1), x(2), co, ostate%f)
END IF
IF (ostate%state == 2) THEN
pe(1:num_gto) = [(x(1)*x(2)**i, i=1, num_gto)]
CALL density_fit(density%op, density%grid, num_gto, pe, co, ostate%f)
END IF

IF (ostate%state == -1) EXIT
IF (ostate%state == -1) EXIT

CALL powell_optimize(ostate%nvar, x, ostate)
CALL powell_optimize(ostate%nvar, x, ostate)

IF (MOD(ostate%nf, n10) == 0 .AND. iunit > 0) THEN
WRITE (iunit, '(" POWELL| Reached",i4,"% of maximal function calls")') &
INT(REAL(ostate%nf, dp)/REAL(ostate%maxfun, dp)*100._dp)
END IF
IF (MOD(ostate%nf, n10) == 0 .AND. iunit > 0) THEN
WRITE (iunit, '(" POWELL| Reached",i4,"% of maximal function calls")') &
INT(REAL(ostate%nf, dp)/REAL(ostate%maxfun, dp)*100._dp)
END IF

END DO
END DO

ostate%state = 8
CALL powell_optimize(ostate%nvar, x, ostate)
CALL density_fit(density, atom, num_gto, x(1), x(2), co, ostate%f)
ostate%state = 8
CALL powell_optimize(ostate%nvar, x, ostate)
pe(1:num_gto) = [(x(1)*x(2)**i, i=1, num_gto)]
END IF

CALL density_fit(density%op, density%grid, num_gto, pe, co, ostate%f)

CALL release_opgrid(density)

IF (iunit > 0) THEN
IF (iunit > 0 .AND. .NOT. PRESENT(agto)) THEN
WRITE (iunit, '(" POWELL| Number of function evaluations",T71,I10)') ostate%nf
WRITE (iunit, '(" POWELL| Final value of function",T61,G20.10)') ostate%fopt
WRITE (iunit, '(" Optimized representation of density using a Geometrical Gaussian basis")')
Expand All @@ -154,13 +166,18 @@ SUBROUTINE atom_fit_density(atom, num_gto, norder, iunit, powell_section, result
END IF

IF (PRESENT(results)) THEN
CPASSERT(SIZE(results) >= num_gto + 2)
results(1) = x(1)
results(2) = x(2)
results(3:2 + num_gto) = co(1:num_gto)
IF (PRESENT(agto)) THEN
CPASSERT(SIZE(results) >= num_gto)
results(1:num_gto) = co(1:num_gto)
ELSE
CPASSERT(SIZE(results) >= num_gto + 2)
results(1) = x(1)
results(2) = x(2)
results(3:2 + num_gto) = co(1:num_gto)
END IF
END IF

DEALLOCATE (co)
DEALLOCATE (co, pe)

END SUBROUTINE atom_fit_density
! **************************************************************************************************
Expand Down Expand Up @@ -1286,50 +1303,46 @@ END SUBROUTINE opt_nlcc_param
! **************************************************************************************************
!> \brief Low level routine to fit an atomic electron density.
!> \param density electron density on an atomic radial grid 'atom%basis%grid'
!> \param atom information about the atomic kind
!> (only 'atom%basis%grid' is accessed, why not to pass it instead?)
!> \param grid information about the atomic grid
!> \param n number of Gaussian basis functions
!> \param aval exponent of the first Gaussian basis function in the series
!> \param cval factor of geometrical series
!> \param pe exponents of Gaussian basis functions
!> \param co computed expansion coefficients
!> \param aerr error in fitted density \f[\sqrt{\sum_{i}^{nr} (density_fitted(i)-density(i))^2 }\f]
! **************************************************************************************************
SUBROUTINE density_fit(density, atom, n, aval, cval, co, aerr)
TYPE(opgrid_type), POINTER :: density
TYPE(atom_type), POINTER :: atom
SUBROUTINE density_fit(density, grid, n, pe, co, aerr)
REAL(dp), DIMENSION(:), INTENT(IN) :: density
TYPE(grid_atom_type) :: grid
INTEGER, INTENT(IN) :: n
REAL(dp), INTENT(IN) :: aval, cval
REAL(dp), DIMENSION(:), INTENT(IN) :: pe
REAL(dp), DIMENSION(:), INTENT(INOUT) :: co
REAL(dp), INTENT(OUT) :: aerr

INTEGER :: i, info, ip, iq, k, nr
INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv
REAL(dp) :: a, rk, zval
REAL(dp), ALLOCATABLE, DIMENSION(:) :: den, pe, uval
REAL(dp), ALLOCATABLE, DIMENSION(:) :: den, uval
REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: bf, smat, tval

ALLOCATE (pe(n))
nr = atom%basis%grid%nr
nr = grid%nr
ALLOCATE (bf(nr, n), den(nr))
bf = 0._dp

DO i = 1, n
pe(i) = aval*cval**i
a = pe(i)
DO k = 1, nr
rk = atom%basis%grid%rad(k)
rk = grid%rad(k)
bf(k, i) = EXP(-a*rk**2)
END DO
END DO

! total charge
zval = integrate_grid(density%op, atom%basis%grid)
zval = integrate_grid(density, grid)

! allocate vectors and matrices for overlaps
ALLOCATE (tval(n + 1, 1), uval(n), smat(n + 1, n + 1))
DO i = 1, n
uval(i) = (pi/pe(i))**1.5_dp
tval(i, 1) = integrate_grid(density%op, bf(:, i), atom%basis%grid)
tval(i, 1) = integrate_grid(density, bf(:, i), grid)
END DO
tval(n + 1, 1) = zval

Expand All @@ -1354,10 +1367,10 @@ SUBROUTINE density_fit(density, atom, n, aval, cval, co, aerr)
den(:) = den(:) + co(i)*bf(:, i)
END DO
den(:) = den(:)*fourpi
den(:) = (den(:) - density%op(:))**2
aerr = SQRT(integrate_grid(den, atom%basis%grid))
den(:) = (den(:) - density(:))**2
aerr = SQRT(integrate_grid(den, grid))

DEALLOCATE (pe, bf, den)
DEALLOCATE (bf, den)

DEALLOCATE (tval, uval, smat)

Expand Down
47 changes: 30 additions & 17 deletions src/atom_kind_orbitals.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -472,25 +472,27 @@ END SUBROUTINE calculate_atomic_orbitals
!> \param qs_kind ...
!> \param ngto ...
!> \param iunit ...
!> \param optbasis ... Default=T, if basis should be optimized, if not basis is given in input (density)
!> \param allelectron ...
!> \param confine ...
! **************************************************************************************************
SUBROUTINE calculate_atomic_density(density, atomic_kind, qs_kind, ngto, iunit, allelectron, confine)
SUBROUTINE calculate_atomic_density(density, atomic_kind, qs_kind, ngto, iunit, &
optbasis, allelectron, confine)
REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: density
TYPE(atomic_kind_type), POINTER :: atomic_kind
TYPE(qs_kind_type), POINTER :: qs_kind
INTEGER, INTENT(IN) :: ngto
INTEGER, INTENT(IN), OPTIONAL :: iunit
LOGICAL, INTENT(IN), OPTIONAL :: allelectron, confine
LOGICAL, INTENT(IN), OPTIONAL :: optbasis, allelectron, confine
INTEGER, PARAMETER :: num_gto = 40
INTEGER :: i, ii, k, l, ll, m, mb, mo, ngp, nn, nr, &
quadtype, relativistic, z
INTEGER :: i, ii, iw, k, l, ll, m, mb, mo, ngp, nn, &
nr, quadtype, relativistic, z
INTEGER, DIMENSION(0:lmat) :: starti
INTEGER, DIMENSION(0:lmat, 10) :: ncalc, ncore, nelem
INTEGER, DIMENSION(:), POINTER :: econf
LOGICAL :: ecp_semi_local
LOGICAL :: do_basopt, ecp_semi_local
REAL(KIND=dp) :: al, aval, cc, cval, ear, rk, xx, zeff
REAL(KIND=dp), DIMENSION(num_gto+2) :: results
TYPE(all_potential_type), POINTER :: all_potential
Expand All @@ -513,13 +515,24 @@ SUBROUTINE calculate_atomic_density(density, atomic_kind, qs_kind, ngto, iunit,
gth_potential=gth_potential, &
sgp_potential=sgp_potential)
IF (PRESENT(iunit)) THEN
iw = iunit
ELSE
iw = -1
END IF
IF (PRESENT(allelectron)) THEN
IF (allelectron) THEN
NULLIFY (gth_potential)
zeff = z
END IF
END IF
do_basopt = .TRUE.
IF (PRESENT(optbasis)) THEN
do_basopt = optbasis
END IF
CPASSERT(ngto <= num_gto)
IF (ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) THEN
Expand Down Expand Up @@ -715,21 +728,21 @@ SUBROUTINE calculate_atomic_density(density, atomic_kind, qs_kind, ngto, iunit,
CALL create_atom_orbs(orbitals, mb, mo)
CALL set_atom(atom, orbitals=orbitals)
IF (PRESENT(iunit)) THEN
CALL calculate_atom(atom, iunit)
CALL atom_fit_density(atom, ngto, 0, iunit, results=results)
CALL calculate_atom(atom, iw)
IF (do_basopt) THEN
CALL atom_fit_density(atom, ngto, 0, iw, results=results)
xx = results(1)
cc = results(2)
DO i = 1, ngto
density(i, 1) = xx*cc**i
density(i, 2) = results(2 + i)
END DO
ELSE
CALL calculate_atom(atom, -1)
CALL atom_fit_density(atom, ngto, 0, -1, results=results)
CALL atom_fit_density(atom, ngto, 0, iw, agto=density(:, 1), results=results)
density(1:ngto, 2) = results(1:ngto)
END IF
xx = results(1)
cc = results(2)
DO i = 1, ngto
density(i, 1) = xx*cc**i
density(i, 2) = results(2 + i)
END DO
! clean up
CALL atom_int_release(integrals)
CALL atom_ppint_release(integrals)
Expand Down
Loading

0 comments on commit f2bc97d

Please sign in to comment.