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

Atom code: ADMM analysis

Browse source 
svn-origin-rev: 17119
This commit is contained in:
Jürg Hutter 2016-07-18 08:03:40 +00:00
parent 60b02676b4
commit b39e976d1b
10 changed files with 691 additions and 77 deletions
Download patch file Download diff file Expand all files Collapse all files

15
src/atom_energy.F
View file

@ -5,6 +5,7 @@
! **************************************************************************************************
MODULE atom_energy
USE atom_admm_methods, ONLY: atom_admm
USE atom_electronic_structure, ONLY: calculate_atom
USE atom_fit, ONLY: atom_fit_density,&
atom_fit_kgpot
@ -95,9 +96,9 @@ CONTAINS
TYPE(atom_potential_type), POINTER :: ae_pot, p_pot
TYPE(atom_state), POINTER :: state
TYPE(cp_logger_type), POINTER :: logger
TYPE(section_vals_type), POINTER :: basis_section, external_vxc_section, method_section, &
opt_section, potential_section, powell_section, xc_section, zmp_restart_section, &
zmp_section
TYPE(section_vals_type), POINTER :: admm_section, basis_section, external_vxc_section, &
method_section, opt_section, potential_section, powell_section, xc_section, &
zmp_restart_section, zmp_section
CALL timeset(routineN, handle)
@ -449,6 +450,14 @@ CONTAINS
END IF
CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%ANALYZE_BASIS")
! Analyse ADMM approximation
iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%ADMM", extension=".log")
admm_section => section_vals_get_subs_vals(atom_section, "PRINT%ADMM")
IF (iw > 0) THEN
CALL atom_admm(atom_info, admm_section, iw)
END IF
CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%ADMM")
! clean up
IF (had_ae) THEN
CALL atom_int_release(ae_int)

2
src/atom_fit.F
View file

@ -846,7 +846,7 @@ CONTAINS
fopt = HUGE(0._dp)
ostate%state = 0
CALL powell_optimize(ostate%nvar, x, ostate)
DO
IF (ostate%state == 2) THEN

83
src/atom_operators.F
View file

@ -15,9 +15,9 @@ MODULE atom_operators
sg_coulomb, sg_erf, sg_erfc, sg_exchange, sg_gpot, sg_kinetic, sg_kinnuc, sg_nuclear, &
sg_overlap, sg_proj_ol, sto_kinetic, sto_nuclear, sto_overlap
USE atom_types, ONLY: &
atom_basis_gridrep, atom_basis_type, atom_integrals, atom_potential_type, atom_state, &
cgto_basis, ecp_pseudo, gth_pseudo, gto_basis, no_pseudo, num_basis, release_atom_basis, &
sto_basis, upf_pseudo
atom_basis_gridrep, atom_basis_type, atom_compare_grids, atom_integrals, &
atom_potential_type, atom_state, cgto_basis, ecp_pseudo, gth_pseudo, gto_basis, no_pseudo, &
num_basis, release_atom_basis, sto_basis, upf_pseudo
USE atom_utils, ONLY: &
atom_solve, contract2, contract2add, contract4, coulomb_potential_numeric, integrate_grid, &
numpot_matrix, slater_density, wigner_slater_functional
@ -43,7 +43,7 @@ MODULE atom_operators
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_operators'
PUBLIC :: atom_int_setup, atom_ppint_setup, atom_int_release, atom_ppint_release
PUBLIC :: atom_relint_setup, atom_relint_release
PUBLIC :: atom_relint_setup, atom_relint_release, atom_basis_projection_overlap
! **************************************************************************************************
@ -933,6 +933,81 @@ CONTAINS
END DO
END SUBROUTINE dkh_integrals
! **************************************************************************************************
!> \brief ...
!> \param ovlap ...
!> \param basis_a ...
!> \param basis_b ...
! **************************************************************************************************
SUBROUTINE atom_basis_projection_overlap(ovlap, basis_a, basis_b)
REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(OUT) :: ovlap
TYPE(atom_basis_type), INTENT(IN) :: basis_a, basis_b
INTEGER :: i, j, l, ma, mb, na, nb
LOGICAL :: ebas
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: omat
ebas = (basis_a%basis_type == basis_b%basis_type)
ovlap = 0.0_dp
IF (ebas) THEN
SELECT CASE (basis_a%basis_type)
CASE DEFAULT
CPABORT("")
CASE (GTO_BASIS)
DO l = 0, 3
na = basis_a%nbas(l)
nb = basis_b%nbas(l)
CALL sg_overlap(ovlap(1:na, 1:nb, l), l, basis_a%am(1:na, l), basis_b%am(1:nb, l))
END DO
CASE (CGTO_BASIS)
DO l = 0, 3
na = basis_a%nbas(l)
nb = basis_b%nbas(l)
ma = basis_a%nprim(l)
mb = basis_b%nprim(l)
ALLOCATE (omat(ma, mb))
CALL sg_overlap(omat(1:ma, 1:mb), l, basis_a%am(1:ma, l), basis_b%am(1:mb, l))
ovlap(1:na, 1:nb, l) = MATMUL(TRANSPOSE(basis_a%cm(1:ma, 1:na, l)), &
MATMUL(omat(1:ma, 1:mb), basis_b%cm(1:mb, 1:nb, l)))
DEALLOCATE (omat)
END DO
CASE (STO_BASIS)
DO l = 0, 3
na = basis_a%nbas(l)
nb = basis_b%nbas(l)
CALL sto_overlap(ovlap(1:na, 1:nb, l), basis_a%ns(1:na, l), basis_b%as(1:nb, l), &
basis_a%ns(1:na, l), basis_b%as(1:nb, l))
END DO
CASE (NUM_BASIS)
CPASSERT(atom_compare_grids(basis_a%grid, basis_b%grid))
DO l = 0, 3
na = basis_a%nbas(l)
nb = basis_b%nbas(l)
DO j = 1, nb
DO i = 1, na
ovlap(i, j, l) = integrate_grid(basis_a%bf(:, i, l), basis_b%bf(:, j, l), basis_a%grid)
END DO
END DO
END DO
END SELECT
ELSE
CPASSERT(atom_compare_grids(basis_a%grid, basis_b%grid))
DO l = 0, 3
na = basis_a%nbas(l)
nb = basis_b%nbas(l)
DO j = 1, nb
DO i = 1, na
ovlap(i, j, l) = integrate_grid(basis_a%bf(:, i, l), basis_b%bf(:, j, l), basis_a%grid)
END DO
END DO
END DO
END IF
END SUBROUTINE atom_basis_projection_overlap
! **************************************************************************************************
!> \brief ...
!> \param integrals ...

12
src/atom_output.F
View file

@ -705,25 +705,25 @@ CONTAINS
END IF
WRITE (iw, '(A2,A)') gthpot%symbol, ADJUSTL(TRIM(gthpot%pname))
WRITE (iw, '(4I5)') gthpot%econf
WRITE (iw, '(F18.12,I8,5F18.12)') gthpot%rc, gthpot%ncl, (gthpot%cl(i), i=1, gthpot%ncl)
WRITE (iw, '(F20.14,I8,5F20.14)') gthpot%rc, gthpot%ncl, (gthpot%cl(i), i=1, gthpot%ncl)
IF (gthpot%lpotextended) THEN
WRITE (iw, '(A,I5)') " LPOT", gthpot%nexp_lpot
DO i = 1, gthpot%nexp_lpot
WRITE (iw, '(F18.12,I8,5F18.12)') gthpot%alpha_lpot(i), gthpot%nct_lpot(i), &
WRITE (iw, '(F20.14,I8,5F20.14)') gthpot%alpha_lpot(i), gthpot%nct_lpot(i), &
(gthpot%cval_lpot(j, i), j=1, gthpot%nct_lpot(i))
END DO
END IF
IF (gthpot%lsdpot) THEN
WRITE (iw, '(A,I5)') " LSD ", gthpot%nexp_lsd
DO i = 1, gthpot%nexp_lsd
WRITE (iw, '(F18.12,I8,5F18.12)') gthpot%alpha_lsd(i), gthpot%nct_lsd(i), &
WRITE (iw, '(F20.14,I8,5F20.14)') gthpot%alpha_lsd(i), gthpot%nct_lsd(i), &
(gthpot%cval_lsd(j, i), j=1, gthpot%nct_lsd(i))
END DO
END IF
IF (gthpot%nlcc) THEN
WRITE (iw, '(A,I5)') " NLCC ", gthpot%nexp_nlcc
DO i = 1, gthpot%nexp_nlcc
WRITE (iw, '(F18.12,I8,5F18.12)') gthpot%alpha_nlcc(i), gthpot%nct_nlcc(i), &
WRITE (iw, '(F20.14,I8,5F20.14)') gthpot%alpha_nlcc(i), gthpot%nct_nlcc(i), &
(gthpot%cval_nlcc(j, i)*4.0_dp*pi, j=1, gthpot%nct_nlcc(i))
END DO
END IF
@ -736,9 +736,9 @@ CONTAINS
END DO
WRITE (iw, '(I8)') n
DO i = 0, n-1
WRITE (iw, '(F18.12,I8,5F18.12)') gthpot%rcnl(i), gthpot%nl(i), (gthpot%hnl(1, k, i), k=1, gthpot%nl(i))
WRITE (iw, '(F20.14,I8,5F20.14)') gthpot%rcnl(i), gthpot%nl(i), (gthpot%hnl(1, k, i), k=1, gthpot%nl(i))
DO j = 2, gthpot%nl(i)
WRITE (iw, '(T20,5F18.12)') (gthpot%hnl(j, k, i), k=j, gthpot%nl(i))
WRITE (iw, '(T20,5F20.14)') (gthpot%hnl(j, k, i), k=j, gthpot%nl(i))
END DO
END DO
IF (.NOT. PRESENT(iunit)) CALL close_file(unit_number=iw)

145
src/atom_types.F
View file

@ -65,15 +65,15 @@ MODULE atom_types
INTEGER :: basis_type
INTEGER, DIMENSION(0:3) :: nbas
INTEGER, DIMENSION(0:3) :: nprim
REAL(KIND=dp), DIMENSION(:, :), POINTER :: am !GTO exponents
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: cm !Contraction coeffs
REAL(KIND=dp), DIMENSION(:, :), POINTER :: as !STO exponents
INTEGER, DIMENSION(:, :), POINTER :: ns !STO n-quantum numbers
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: bf !num. bsf
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dbf !derivatives (num)
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: ddbf !2nd derivatives (num)
REAL(KIND=dp), DIMENSION(:, :), POINTER :: am !GTO exponents
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: cm !Contraction coeffs
REAL(KIND=dp), DIMENSION(:, :), POINTER :: as !STO exponents
INTEGER, DIMENSION(:, :), POINTER :: ns !STO n-quantum numbers
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: bf !num. bsf
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dbf !derivatives (num)
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: ddbf !2nd derivatives (num)
REAL(KIND=dp) :: eps_eig
TYPE(grid_atom_type), POINTER :: grid
TYPE(grid_atom_type), POINTER :: grid
LOGICAL :: geometrical
REAL(KIND=dp) :: aval, cval
INTEGER, DIMENSION(0:3) :: start
@ -82,8 +82,8 @@ MODULE atom_types
!> Provides all information about a pseudopotential
! **************************************************************************************************
TYPE atom_gthpot_type
CHARACTER(LEN=2) :: symbol
CHARACTER(LEN=default_string_length) :: pname
CHARACTER(LEN=2) :: symbol
CHARACTER(LEN=default_string_length) :: pname
INTEGER, DIMENSION(0:3) :: econf
REAL(dp) :: zion
REAL(dp) :: rc
@ -91,31 +91,31 @@ MODULE atom_types
REAL(dp), DIMENSION(5) :: cl
INTEGER, DIMENSION(0:3) :: nl
REAL(dp), DIMENSION(0:3) :: rcnl
REAL(dp), DIMENSION(4, 4, 0:3) :: hnl
REAL(dp), DIMENSION(4, 4, 0:3) :: hnl
! type extensions
! NLCC
LOGICAL :: nlcc
INTEGER :: nexp_nlcc
REAL(KIND=dp), DIMENSION(10) :: alpha_nlcc
INTEGER, DIMENSION(10) :: nct_nlcc
REAL(KIND=dp), DIMENSION(4, 10) :: cval_nlcc
LOGICAL :: nlcc
INTEGER :: nexp_nlcc
REAL(KIND=dp), DIMENSION(10) :: alpha_nlcc
INTEGER, DIMENSION(10) :: nct_nlcc
REAL(KIND=dp), DIMENSION(4, 10) :: cval_nlcc
! LSD potential
LOGICAL :: lsdpot
INTEGER :: nexp_lsd
REAL(KIND=dp), DIMENSION(10) :: alpha_lsd
INTEGER, DIMENSION(10) :: nct_lsd
REAL(KIND=dp), DIMENSION(4, 10) :: cval_lsd
LOGICAL :: lsdpot
INTEGER :: nexp_lsd
REAL(KIND=dp), DIMENSION(10) :: alpha_lsd
INTEGER, DIMENSION(10) :: nct_lsd
REAL(KIND=dp), DIMENSION(4, 10) :: cval_lsd
! extended local potential
LOGICAL :: lpotextended
INTEGER :: nexp_lpot
REAL(KIND=dp), DIMENSION(10) :: alpha_lpot
INTEGER, DIMENSION(10) :: nct_lpot
REAL(KIND=dp), DIMENSION(4, 10) :: cval_lpot
LOGICAL :: lpotextended
INTEGER :: nexp_lpot
REAL(KIND=dp), DIMENSION(10) :: alpha_lpot
INTEGER, DIMENSION(10) :: nct_lpot
REAL(KIND=dp), DIMENSION(4, 10) :: cval_lpot
END TYPE atom_gthpot_type
TYPE atom_ecppot_type
CHARACTER(LEN=2) :: symbol
CHARACTER(LEN=default_string_length) :: pname
CHARACTER(LEN=2) :: symbol
CHARACTER(LEN=default_string_length) :: pname
INTEGER, DIMENSION(0:3) :: econf
REAL(dp) :: zion
INTEGER :: lmax
@ -124,9 +124,9 @@ MODULE atom_types
REAL(dp), DIMENSION(1:10) :: aloc ! coefficient
REAL(dp), DIMENSION(1:10) :: bloc ! exponent
INTEGER, DIMENSION(0:10) :: npot ! # terms
INTEGER, DIMENSION(1:10, 0:10) :: nrpot ! r**(n-2)
REAL(dp), DIMENSION(1:10, 0:10) :: apot ! coefficient
REAL(dp), DIMENSION(1:10, 0:10) :: bpot ! exponent
INTEGER, DIMENSION(1:10, 0:10) :: nrpot ! r**(n-2)
REAL(dp), DIMENSION(1:10, 0:10) :: apot ! coefficient
REAL(dp), DIMENSION(1:10, 0:10) :: bpot ! exponent
END TYPE atom_ecppot_type
TYPE atom_potential_type
@ -144,21 +144,21 @@ MODULE atom_types
!> Provides all information on states and occupation
! **************************************************************************************************
TYPE atom_state
REAL(KIND=dp), DIMENSION(0:3, 10) :: occ
REAL(KIND=dp), DIMENSION(0:3, 10) :: core
REAL(KIND=dp), DIMENSION(0:3, 10) :: occupation
REAL(KIND=dp), DIMENSION(0:3, 10) :: occ
REAL(KIND=dp), DIMENSION(0:3, 10) :: core
REAL(KIND=dp), DIMENSION(0:3, 10) :: occupation
INTEGER :: maxl_occ
INTEGER, DIMENSION(0:3) :: maxn_occ
INTEGER, DIMENSION(0:3) :: maxn_occ
INTEGER :: maxl_calc
INTEGER, DIMENSION(0:3) :: maxn_calc
INTEGER, DIMENSION(0:3) :: maxn_calc
INTEGER :: multiplicity
REAL(KIND=dp), DIMENSION(0:3, 10) :: occa, occb
REAL(KIND=dp), DIMENSION(0:3, 10) :: occa, occb
END TYPE atom_state
!> Holds atomic integrals
! **************************************************************************************************
TYPE eri
REAL(KIND=dp), DIMENSION(:, :), POINTER :: int
REAL(KIND=dp), DIMENSION(:, :), POINTER :: int
END TYPE eri
TYPE atom_integrals
@ -168,43 +168,43 @@ MODULE atom_types
LOGICAL :: eri_exchange
LOGICAL :: all_nu
INTEGER, DIMENSION(0:3) :: n, nne
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: ovlp, kin, core, clsd
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: utrans, uptrans
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: hnl
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: conf
TYPE(eri), DIMENSION(20) :: ceri
TYPE(eri), DIMENSION(20) :: eeri
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: ovlp, kin, core, clsd
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: utrans, uptrans
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: hnl
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: conf
TYPE(eri), DIMENSION(20) :: ceri
TYPE(eri), DIMENSION(20) :: eeri
INTEGER :: dkhstat = 0
INTEGER :: zorastat = 0
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: tzora
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: hdkh
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: tzora
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: hdkh
END TYPE atom_integrals
!> Holds atomic orbitals and energies
! **************************************************************************************************
TYPE atom_orbitals
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: wfn, wfna, wfnb
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: pmat, pmata, pmatb
REAL(KIND=dp), DIMENSION(:, :), POINTER :: ener, enera, enerb
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: refene, refchg, refnod
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: wrefene, wrefchg, wrefnod
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: crefene, crefchg, crefnod
REAL(KIND=dp), DIMENSION(:, :), POINTER :: wpsir0
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: rcmax
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: wfn, wfna, wfnb
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: pmat, pmata, pmatb
REAL(KIND=dp), DIMENSION(:, :), POINTER :: ener, enera, enerb
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: refene, refchg, refnod
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: wrefene, wrefchg, wrefnod
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: crefene, crefchg, crefnod
REAL(KIND=dp), DIMENSION(:, :), POINTER :: wpsir0
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: rcmax
END TYPE atom_orbitals
!> Operator matrices
! **************************************************************************************************
TYPE opmat_type
INTEGER, DIMENSION(0:3) :: n
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: op
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: op
END TYPE opmat_type
!> Operator grids
! **************************************************************************************************
TYPE opgrid_type
REAL(KIND=dp), DIMENSION(:), POINTER :: op
TYPE(grid_atom_type), POINTER :: grid
REAL(KIND=dp), DIMENSION(:), POINTER :: op
TYPE(grid_atom_type), POINTER :: grid
END TYPE opgrid_type
!> All energies
@ -273,8 +273,9 @@ MODULE atom_types
PUBLIC :: atom_p_type, atom_type, atom_basis_type, atom_state, atom_integrals
PUBLIC :: atom_orbitals, eri, atom_potential_type, atom_gthpot_type
PUBLIC :: atom_optimization_type
PUBLIC :: atom_compare_grids
PUBLIC :: create_atom_type, release_atom_type, set_atom
PUBLIC :: create_atom_orbs
PUBLIC :: create_atom_orbs, release_atom_orbs
PUBLIC :: init_atom_basis, atom_basis_gridrep, release_atom_basis
PUBLIC :: init_atom_potential, release_atom_potential
PUBLIC :: read_atom_opt_section
@ -2735,5 +2736,33 @@ CONTAINS
END SUBROUTINE read_ecp_potential
! **************************************************************************************************
!> \brief ...
!> \param grid1 ...
!> \param grid2 ...
!> \retval is_equal ...
! **************************************************************************************************
FUNCTION atom_compare_grids(grid1, grid2) RESULT(is_equal)
TYPE(grid_atom_type) :: grid1, grid2
LOGICAL :: is_equal
INTEGER :: i
REAL(KIND=dp) :: dr, dw
is_equal = .TRUE.
IF (grid1%nr == grid2%nr) THEN
DO i = 1, grid2%nr
dr = ABS(grid1%rad(i)-grid2%rad(i))
dw = ABS(grid1%wr(i)-grid2%wr(i))
IF (dr+dw > 1.0e-12_dp) THEN
is_equal = .FALSE.
EXIT
END IF
END DO
ELSE
is_equal = .FALSE.
END IF
END FUNCTION atom_compare_grids
! **************************************************************************************************
END MODULE atom_types

417
src/atom_xc.F
View file

@ -9,7 +9,8 @@
! **************************************************************************************************
MODULE atom_xc
USE atom_types, ONLY: atom_type,&
USE atom_types, ONLY: atom_basis_type,&
atom_type,&
gth_pseudo,&
opmat_type,&
upf_pseudo
@ -27,7 +28,8 @@ MODULE atom_xc
section_vals_type,&
section_vals_val_get
USE kinds, ONLY: dp
USE mathconstants, ONLY: fourpi
USE mathconstants, ONLY: fourpi,&
pi
USE xc_atom, ONLY: xc_rho_set_atom_update
USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
xc_dset_create,&
@ -53,6 +55,7 @@ MODULE atom_xc
! ZMP added public subroutines calculate_atom_zmp, calculate_atom_ext_vxc
PUBLIC :: calculate_atom_vxc_lda, calculate_atom_vxc_lsd, &
calculate_atom_zmp, calculate_atom_ext_vxc
PUBLIC :: atom_vxc_lda, atom_vxc_lsd, atom_dpot_lda
! **************************************************************************************************
@ -482,9 +485,9 @@ CONTAINS
END IF
IF (needs%laplace_rho_spin) THEN
ALLOCATE (lap(nr, 2))
CALL atom_density(lap(:, 1), atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, &
CALL atom_density(lap(:, 1), atom%orbitals%pmata, atom%basis, atom%state%maxl_occ, &
typ="LAP", rr=atom%basis%grid%rad2)
CALL atom_density(lap(:, 2), atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, &
CALL atom_density(lap(:, 2), atom%orbitals%pmatb, atom%basis, atom%state%maxl_occ, &
typ="LAP", rr=atom%basis%grid%rad2)
IF (nlcc) THEN
xfun = 0.0_dp
@ -607,6 +610,412 @@ CONTAINS
END SUBROUTINE calculate_atom_vxc_lsd
! **************************************************************************************************
!> \brief ...
!> \param basis ...
!> \param pmat ...
!> \param maxl ...
!> \param xc_section ...
!> \param fexc ...
!> \param xcmat ...
! **************************************************************************************************
SUBROUTINE atom_vxc_lda(basis, pmat, maxl, xc_section, fexc, xcmat)
TYPE(atom_basis_type), POINTER :: basis
REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: pmat
INTEGER, INTENT(IN) :: maxl
TYPE(section_vals_type), POINTER :: xc_section
REAL(KIND=dp), INTENT(OUT) :: fexc
REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: xcmat
CHARACTER(LEN=*), PARAMETER :: routineN = 'atom_vxc_lda', routineP = moduleN//':'//routineN
INTEGER :: deriv_order, handle, l, myfun, n1, n2, &
n3, nr, nspins
INTEGER, DIMENSION(2, 3) :: bounds
LOGICAL :: lsd
REAL(KIND=dp) :: density_cut, gradient_cut, tau_cut
REAL(KIND=dp), DIMENSION(:, :), POINTER :: drho, lap, rho, tau
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: taumat, xcpot
TYPE(section_vals_type), POINTER :: xc_fun_section
TYPE(xc_derivative_set_type), POINTER :: deriv_set
TYPE(xc_derivative_type), POINTER :: deriv
TYPE(xc_rho_cflags_type) :: needs
TYPE(xc_rho_set_type), POINTER :: rho_set
! -------------------------------------------------------------------------
CALL timeset(routineN, handle)
CPASSERT(ASSOCIATED(xc_section))
NULLIFY (rho_set)
xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
IF (myfun == xc_none) THEN
fexc = 0._dp
ELSE
CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", r_val=density_cut)
CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", r_val=gradient_cut)
CALL section_vals_val_get(xc_section, "TAU_CUTOFF", r_val=tau_cut)
lsd = .FALSE.
nspins = 1
needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, add_basic_components=.FALSE.)
! Prepare the structures needed to calculate and store the xc derivatives
! Array dimension: here anly one dimensional arrays are used,
! i.e. only the first column of deriv_data is read.
! The other to dimensions are set to size equal 1
nr = basis%grid%nr
bounds(1:2, 1:3) = 1
bounds(2, 1) = nr
! create a place where to put the derivatives
NULLIFY (deriv_set)
CALL xc_dset_create(deriv_set, local_bounds=bounds)
! create the place where to store the argument for the functionals
CALL xc_rho_set_create(rho_set, bounds, rho_cutoff=density_cut, &
drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
! allocate the required 3d arrays where to store rho and drho
CALL xc_rho_set_atom_update(rho_set, needs, nspins, bounds)
NULLIFY (rho, drho, tau)
IF (needs%rho) THEN
ALLOCATE (rho(nr, 1))
CALL atom_density(rho(:, 1), pmat, basis, maxl, typ="RHO")
END IF
IF (needs%norm_drho) THEN
ALLOCATE (drho(nr, 1))
CALL atom_density(drho(:, 1), pmat, basis, maxl, typ="DER")
END IF
IF (needs%tau) THEN
ALLOCATE (tau(nr, 1))
CALL atom_density(tau(:, 1), pmat, basis, maxl, typ="KIN", rr=basis%grid%rad2)
END IF
IF (needs%laplace_rho) THEN
ALLOCATE (lap(nr, 1))
CALL atom_density(lap(:, 1), pmat, basis, maxl, typ="LAP", rr=basis%grid%rad2)
END IF
CALL fill_rho_set(rho_set, nspins, needs, rho, drho, tau, lap, nr)
CALL xc_dset_zero_all(deriv_set)
deriv_order = 1
CALL xc_functionals_eval(xc_fun_section, lsd=lsd, rho_set=rho_set, deriv_set=deriv_set, &
deriv_order=deriv_order)
! Integration to get the matrix elements and energy
deriv => xc_dset_get_derivative(deriv_set, "", allocate_deriv=.FALSE.)
CALL xc_derivative_get(deriv, deriv_data=xcpot)
fexc = fourpi*integrate_grid(xcpot(:, 1, 1), basis%grid)
IF (needs%rho) THEN
deriv => xc_dset_get_derivative(deriv_set, "(rho)", allocate_deriv=.FALSE.)
CALL xc_derivative_get(deriv, deriv_data=xcpot)
CALL numpot_matrix(xcmat, xcpot(:, 1, 1), basis, 0)
DEALLOCATE (rho)
END IF
IF (needs%norm_drho) THEN
deriv => xc_dset_get_derivative(deriv_set, "(norm_drho)", allocate_deriv=.FALSE.)
CALL xc_derivative_get(deriv, deriv_data=xcpot)
CALL numpot_matrix(xcmat, xcpot(:, 1, 1), basis, 1)
DEALLOCATE (drho)
END IF
IF (needs%tau) THEN
deriv => xc_dset_get_derivative(deriv_set, "(tau)", allocate_deriv=.FALSE.)
CALL xc_derivative_get(deriv, deriv_data=xcpot)
n1 = SIZE(xcmat, 1)
n2 = SIZE(xcmat, 2)
n3 = SIZE(xcmat, 3)
ALLOCATE (taumat(n1, n2, 0:n3-1))
taumat = 0._dp
xcpot(:, 1, 1) = 0.5_dp*xcpot(:, 1, 1)
CALL numpot_matrix(xcmat, xcpot(:, 1, 1), basis, 2)
xcpot(:, 1, 1) = xcpot(:, 1, 1)/basis%grid%rad2(:)
CALL numpot_matrix(taumat, xcpot(:, 1, 1), basis, 0)
DO l = 0, 3
xcmat(:, :, l) = xcmat(:, :, l)+REAL(l*(l+1), dp)*taumat(:, :, l)
END DO
DEALLOCATE (tau)
DEALLOCATE (taumat)
END IF
! Release the xc structure used to store the xc derivatives
CALL xc_dset_release(deriv_set)
CALL xc_rho_set_release(rho_set)
END IF !xc_none
CALL timestop(handle)
END SUBROUTINE atom_vxc_lda
! **************************************************************************************************
!> \brief ...
!> \param basis ...
!> \param pmata ...
!> \param pmatb ...
!> \param maxl ...
!> \param xc_section ...
!> \param fexc ...
!> \param xcmata ...
!> \param xcmatb ...
! **************************************************************************************************
SUBROUTINE atom_vxc_lsd(basis, pmata, pmatb, maxl, xc_section, fexc, xcmata, xcmatb)
TYPE(atom_basis_type), POINTER :: basis
REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: pmata, pmatb
INTEGER, INTENT(IN) :: maxl
TYPE(section_vals_type), POINTER :: xc_section
REAL(KIND=dp), INTENT(OUT) :: fexc
REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: xcmata, xcmatb
CHARACTER(LEN=*), PARAMETER :: routineN = 'atom_vxc_lsd', routineP = moduleN//':'//routineN
INTEGER :: deriv_order, handle, l, myfun, n1, n2, &
n3, nr, nspins
INTEGER, DIMENSION(2, 3) :: bounds
LOGICAL :: lsd
REAL(KIND=dp) :: density_cut, gradient_cut, tau_cut
REAL(KIND=dp), DIMENSION(:, :), POINTER :: drho, lap, rho, tau
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: taumat, xcpot
TYPE(section_vals_type), POINTER :: xc_fun_section
TYPE(xc_derivative_set_type), POINTER :: deriv_set
TYPE(xc_derivative_type), POINTER :: deriv
TYPE(xc_rho_cflags_type) :: needs
TYPE(xc_rho_set_type), POINTER :: rho_set
! -------------------------------------------------------------------------
CALL timeset(routineN, handle)
CPASSERT(ASSOCIATED(xc_section))
NULLIFY (rho_set)
xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
IF (myfun == xc_none) THEN
fexc = 0._dp
ELSE
CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", r_val=density_cut)
CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", r_val=gradient_cut)
CALL section_vals_val_get(xc_section, "TAU_CUTOFF", r_val=tau_cut)
lsd = .TRUE.
nspins = 2
needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, add_basic_components=.FALSE.)
! Prepare the structures needed to calculate and store the xc derivatives
! Array dimension: here anly one dimensional arrays are used,
! i.e. only the first column of deriv_data is read.
! The other to dimensions are set to size equal 1
nr = basis%grid%nr
bounds(1:2, 1:3) = 1
bounds(2, 1) = nr
! create a place where to put the derivatives
NULLIFY (deriv_set)
CALL xc_dset_create(deriv_set, local_bounds=bounds)
! create the place where to store the argument for the functionals
CALL xc_rho_set_create(rho_set, bounds, rho_cutoff=density_cut, &
drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
! allocate the required 3d arrays where to store rho and drho
CALL xc_rho_set_atom_update(rho_set, needs, nspins, bounds)
NULLIFY (rho, drho, tau)
IF (needs%rho_spin) THEN
ALLOCATE (rho(nr, 2))
CALL atom_density(rho(:, 1), pmata, basis, maxl, typ="RHO")
CALL atom_density(rho(:, 2), pmatb, basis, maxl, typ="RHO")
END IF
IF (needs%norm_drho_spin) THEN
ALLOCATE (drho(nr, 2))
CALL atom_density(drho(:, 1), pmata, basis, maxl, typ="DER")
CALL atom_density(drho(:, 2), pmatb, basis, maxl, typ="DER")
END IF
IF (needs%tau_spin) THEN
ALLOCATE (tau(nr, 2))
CALL atom_density(tau(:, 1), pmata, basis, maxl, typ="KIN", rr=basis%grid%rad2)
CALL atom_density(tau(:, 2), pmatb, basis, maxl, typ="KIN", rr=basis%grid%rad2)
END IF
IF (needs%laplace_rho_spin) THEN
ALLOCATE (lap(nr, 2))
CALL atom_density(lap(:, 1), pmata, basis, maxl, typ="LAP", rr=basis%grid%rad2)
CALL atom_density(lap(:, 2), pmatb, basis, maxl, typ="LAP", rr=basis%grid%rad2)
END IF
CALL fill_rho_set(rho_set, nspins, needs, rho, drho, tau, lap, nr)
CALL xc_dset_zero_all(deriv_set)
deriv_order = 1
CALL xc_functionals_eval(xc_fun_section, lsd=lsd, rho_set=rho_set, deriv_set=deriv_set, &
deriv_order=deriv_order)
! Integration to get the matrix elements and energy
deriv => xc_dset_get_derivative(deriv_set, "", allocate_deriv=.FALSE.)
CALL xc_derivative_get(deriv, deriv_data=xcpot)
fexc = fourpi*integrate_grid(xcpot(:, 1, 1), basis%grid)
IF (needs%rho_spin) THEN
deriv => xc_dset_get_derivative(deriv_set, "(rhoa)", allocate_deriv=.FALSE.)
CALL xc_derivative_get(deriv, deriv_data=xcpot)
CALL numpot_matrix(xcmata, xcpot(:, 1, 1), basis, 0)
deriv => xc_dset_get_derivative(deriv_set, "(rhob)", allocate_deriv=.FALSE.)
CALL xc_derivative_get(deriv, deriv_data=xcpot)
CALL numpot_matrix(xcmatb, xcpot(:, 1, 1), basis, 0)
DEALLOCATE (rho)
END IF
IF (needs%norm_drho_spin) THEN
! drhoa
NULLIFY (deriv)
deriv => xc_dset_get_derivative(deriv_set, "(norm_drhoa)", allocate_deriv=.FALSE.)
CALL xc_derivative_get(deriv, deriv_data=xcpot)
CALL numpot_matrix(xcmata, xcpot(:, 1, 1), basis, 1)
! drhob
NULLIFY (deriv)
deriv => xc_dset_get_derivative(deriv_set, "(norm_drhob)", allocate_deriv=.FALSE.)
CALL xc_derivative_get(deriv, deriv_data=xcpot)
CALL numpot_matrix(xcmatb, xcpot(:, 1, 1), basis, 1)
! Cross Terms
NULLIFY (deriv)
deriv => xc_dset_get_derivative(deriv_set, "(norm_drho)")
IF (ASSOCIATED(deriv)) THEN
CALL xc_derivative_get(deriv, deriv_data=xcpot)
CALL numpot_matrix(xcmata, xcpot(:, 1, 1), basis, 1)
CALL numpot_matrix(xcmatb, xcpot(:, 1, 1), basis, 1)
END IF
DEALLOCATE (drho)
END IF
IF (needs%tau_spin) THEN
n1 = SIZE(xcmata, 1)
n2 = SIZE(xcmata, 2)
n3 = SIZE(xcmata, 3)
ALLOCATE (taumat(n1, n2, 0:n3-1))
deriv => xc_dset_get_derivative(deriv_set, "(tau_a)", allocate_deriv=.FALSE.)
CALL xc_derivative_get(deriv, deriv_data=xcpot)
taumat = 0._dp
xcpot(:, 1, 1) = 0.5_dp*xcpot(:, 1, 1)
CALL numpot_matrix(xcmata, xcpot(:, 1, 1), basis, 2)
xcpot(:, 1, 1) = xcpot(:, 1, 1)/basis%grid%rad2(:)
CALL numpot_matrix(taumat, xcpot(:, 1, 1), basis, 0)
DO l = 0, 3
xcmata(:, :, l) = xcmata(:, :, l)+REAL(l*(l+1), dp)*taumat(:, :, l)
END DO
deriv => xc_dset_get_derivative(deriv_set, "(tau_b)", allocate_deriv=.FALSE.)
CALL xc_derivative_get(deriv, deriv_data=xcpot)
taumat = 0._dp
xcpot(:, 1, 1) = 0.5_dp*xcpot(:, 1, 1)
CALL numpot_matrix(xcmatb, xcpot(:, 1, 1), basis, 2)
xcpot(:, 1, 1) = xcpot(:, 1, 1)/basis%grid%rad2(:)
CALL numpot_matrix(taumat, xcpot(:, 1, 1), basis, 0)
DO l = 0, 3
xcmatb(:, :, l) = xcmatb(:, :, l)+REAL(l*(l+1), dp)*taumat(:, :, l)
END DO
DEALLOCATE (tau)
DEALLOCATE (taumat)
END IF
! Release the xc structure used to store the xc derivatives
CALL xc_dset_release(deriv_set)
CALL xc_rho_set_release(rho_set)
END IF !xc_none
CALL timestop(handle)
END SUBROUTINE atom_vxc_lsd
! **************************************************************************************************
!> \brief ...
!> \param basis0 ...
!> \param pmat0 ...
!> \param basis1 ...
!> \param pmat1 ...
!> \param maxl ...
!> \param functional ...
!> \param dfexc ...
!> \param linxpar ...
! **************************************************************************************************
SUBROUTINE atom_dpot_lda(basis0, pmat0, basis1, pmat1, maxl, functional, dfexc, linxpar)
TYPE(atom_basis_type), POINTER :: basis0
REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: pmat0
TYPE(atom_basis_type), POINTER :: basis1
REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: pmat1
INTEGER, INTENT(IN) :: maxl
CHARACTER(LEN=*), INTENT(IN) :: functional
REAL(KIND=dp), INTENT(OUT) :: dfexc
REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: linxpar
CHARACTER(LEN=*), PARAMETER :: routineN = 'atom_dpot_lda', routineP = moduleN//':'//routineN
REAL(KIND=dp), PARAMETER :: alx = 0.20_dp, blx = -0.06_dp, &
fs = -0.75_dp*(3._dp/pi)**(1._dp/3._dp)
INTEGER :: handle, ir, nr
REAL(KIND=dp) :: a, b, fx
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: delta, drho0, drho1, pot0, pot1, rho0, &
rho1
! -------------------------------------------------------------------------
CALL timeset(routineN, handle)
nr = basis0%grid%nr
ALLOCATE (rho0(nr), drho0(nr))
CALL atom_density(rho0, pmat0, basis0, maxl, typ="RHO")
CALL atom_density(drho0, pmat0, basis0, maxl, typ="DER")
!
ALLOCATE (rho1(nr), drho1(nr))
CALL atom_density(rho1, pmat1, basis1, maxl, typ="RHO")
CALL atom_density(drho1, pmat1, basis1, maxl, typ="DER")
!
ALLOCATE (delta(nr))
fx = 4.0_dp/3.0_dp
delta(1:nr) = fs*(rho0(1:nr)**fx-rho1(1:nr)**fx)
SELECT CASE (TRIM (functional))
CASE ("LINX")
IF (PRESENT(linxpar)) THEN
a = linxpar(1)
b = linxpar(2)
ELSE
a = alx
b = blx
END IF
ALLOCATE (pot0(nr), pot1(nr))
DO ir = 1, nr
IF (rho0(ir) > 1.e-12_dp) THEN
pot0(ir) = 0.5_dp*drho0(ir)/(3._dp*pi*pi*rho0(ir)**fx)
ELSE
pot0(ir) = 0._dp
END IF
END DO
pot1(1:nr) = 1._dp+(a*pot0(1:nr)**2)/(1._dp+b*pot0(1:nr)**2)
pot1(1:nr) = pot1(1:nr)*delta(1:nr)
dfexc = fourpi*integrate_grid(pot1(1:nr), basis0%grid)
DEALLOCATE (pot0, pot1)
CASE DEFAULT
CPABORT("Unknown functional in atom_dpot_lda")
END SELECT
DEALLOCATE (rho0, rho1, drho0, drho1, delta)
CALL timestop(handle)
END SUBROUTINE atom_dpot_lda
! **************************************************************************************************
!> \brief ...
!> \param rho_set ...

28
src/input_cp2k_atom.F
View file

@ -196,7 +196,7 @@ CONTAINS
routineP = moduleN//':'//routineN
TYPE(keyword_type), POINTER :: keyword
TYPE(section_type), POINTER :: print_key
TYPE(section_type), POINTER :: print_key, subsection
CPASSERT(.NOT. ASSOCIATED(section))
CALL section_create(section, name="print", &
@ -332,6 +332,32 @@ CONTAINS
print_level=debug_print_level, filename="__STD_OUT__")
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)
! ADMM Analysis
CALL cp_print_key_section_create(print_key, "ADMM", &
description="Analysis of ADMM approximation to exact exchange", &
print_level=high_print_level, filename="__STD_OUT__")
NULLIFY (subsection)
CALL section_create(subsection, name="ADMM_BASIS", &
description="Section of basis set information for ADMM calculations.", &
n_keywords=0, n_subsections=0, repeats=.FALSE.)
CALL atom_basis_section(subsection)
CALL section_add_subsection(print_key, subsection)
CALL section_release(subsection)
! CALL keyword_create(keyword, name="OVERLAP_CONDITION_NUMBER", &
! description="Condition number of the basis set overlap matrix calculated for a cubic crystal", &
! usage="OVERLAP_CONDITION_NUMBER <logical>", type_of_var=logical_t, default_l_val=.FALSE.)
! CALL section_add_keyword(print_key, keyword)
! CALL keyword_release(keyword)
! CALL keyword_create(keyword, name="COMPLETENESS", &
! description="Calculate a completeness estimate for the basis set.", &
! usage="COMPLETENESS <logical>", type_of_var=logical_t, default_l_val=.FALSE.)
! CALL section_add_keyword(print_key, keyword)
! CALL keyword_release(keyword)
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)
END SUBROUTINE create_atom_print_section
! **************************************************************************************************

41
tests/ATOM/regtest-2/C_ADMM.inp Normal file
View file

@ -0,0 +1,41 @@
&GLOBAL
PROGRAM_NAME ATOM
&END GLOBAL
&ATOM
ELEMENT C
RUN_TYPE ENERGY
ELECTRON_CONFIGURATION [He] 2s2 2p2
MAX_ANGULAR_MOMENTUM 1
&METHOD
METHOD_TYPE KOHN-SHAM
RELATIVISTIC DKH(3)
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END METHOD
&OPTIMIZATION
EPS_SCF 1.e-8
&END
&PRINT
&ADMM
&ADMM_BASIS
BASIS_SET_FILE_NAME EMSL_BASIS_SETS
BASIS_TYPE CONTRACTED_GTO
BASIS_SET 3-21Gx
&END ADMM_BASIS
&END ADMM
&END PRINT
&AE_BASIS
BASIS_SET_FILE_NAME EMSL_BASIS_SETS
BASIS_TYPE CONTRACTED_GTO
BASIS_SET aug-cc-pVTZ
&END AE_BASIS
&END ATOM

24
tests/ATOM/regtest-2/DUMMY_BASIS
View file

@ -49,3 +49,27 @@ H QZV3P-GTH
0.8380000000 0.0000000000 1.0000000000 0.0000000000
0.2920000000 0.0000000000 0.0000000000 1.0000000000
#
Ru DUMMY
3
1 0 0 6 4
3.7326050700 -8.760467E-03 -5.122607E-03 -2.691963E-01 -1.532898E+00
1.8341903900 8.361481E-02 -3.252435E-02 1.269982E+00 5.342929E+00
0.8090639000 -2.452690E-01 1.302864E-01 -3.881830E+00 -7.876596E+00
0.3451510100 -2.404606E-01 8.596278E-01 6.248922E+00 5.472671E+00
0.1383665500 6.068951E-01 -2.840153E+00 -4.365986E+00 -2.354735E+00
0.0496701000 6.658400E-01 1.990853E+00 1.370297E+00 5.877098E-01
1 1 1 6 3
3.7326050700 -6.473343E-03 -6.201480E-03 -8.815307E-02
1.8341903900 3.344182E-02 1.171603E-02 3.762492E-01
0.8090639000 -5.292065E-02 -7.176380E-02 -1.155100E+00
0.3451510100 -4.009707E-02 4.493970E-01 2.915465E+00
0.1383665500 2.768114E-01 -1.868171E+00 -2.260303E+00
0.0496701000 8.071936E-01 1.360595E+00 7.696126E-01
1 2 2 6 4
3.7326050700 1.077335E-01 -3.143227E-02 2.571880E-02 6.782125E-02
1.8341903900 -3.465083E-01 8.716166E-02 -6.348120E-02 -3.118865E-01
0.8090639000 -3.933684E-01 1.145400E-01 -4.306530E-01 1.776388E+00
0.3451510100 -3.459824E-01 3.799333E-01 -1.307229E-01 -2.621994E+00
0.1383665500 -1.984013E-01 -4.769718E-01 1.391356E+00 1.735134E+00
0.0496701000 -3.527514E-02 -7.160104E-01 -1.077805E+00 -7.103047E-01

1
tests/ATOM/regtest-2/TEST_FILES
View file

@ -45,4 +45,5 @@ slater_10.inp 35 1e-12
#
Hg.inp 35 1e-12 -19602.261685089714
C.inp 35 1e-12 -37.800000362746
C_ADMM.inp 35 1e-12 -37.758721798835
#EOF