Skip to content
Permalink
Browse files

XAS_TDP| 1) Refactoring speed of light into common/physcon.F 2) Clean-up

  • Loading branch information...
abussy authored and pseewald committed Sep 26, 2019
1 parent 2eb5016 commit e01c96e24dd4e0ea618b6c8c3df6df4cdcb79845
@@ -31,7 +31,7 @@ MODULE atom_operators
sqrt2
USE mathlib, ONLY: jacobi
USE periodic_table, ONLY: ptable
USE physcon, ONLY: a_fine
USE physcon, ONLY: c_light_au
USE qs_grid_atom, ONLY: grid_atom_type
#include "./base/base_uses.f90"

@@ -813,16 +813,14 @@ SUBROUTINE atom_relint_setup(integrals, basis, reltyp, zcore, alpha)
routineP = moduleN//':'//routineN

INTEGER :: dkhorder, handle, i, k1, k2, l, m, n, nl
REAL(dp) :: ascal, cspeed
REAL(dp) :: ascal
REAL(dp), ALLOCATABLE, DIMENSION(:) :: cpot
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: modpot
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ener, sps
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: hmat, pvp, sp, tp, vp, wfn

CALL timeset(routineN, handle)

cspeed = 1._dp/a_fine

SELECT CASE (reltyp)
CASE DEFAULT
CPABORT("")
@@ -856,7 +854,7 @@ SUBROUTINE atom_relint_setup(integrals, basis, reltyp, zcore, alpha)
ALLOCATE (modpot(1:m), cpot(1:m))
CALL calculate_model_potential(modpot, basis%grid, zcore)
! Zora potential
cpot(1:m) = modpot(1:m)/(4._dp*cspeed*cspeed-2._dp*modpot(1:m))
cpot(1:m) = modpot(1:m)/(4._dp*c_light_au*c_light_au-2._dp*modpot(1:m))
cpot(1:m) = cpot(1:m)/basis%grid%rad2(1:m)
CALL numpot_matrix(integrals%tzora, cpot, basis, 0)
DO l = 0, lmat
@@ -875,7 +873,7 @@ SUBROUTINE atom_relint_setup(integrals, basis, reltyp, zcore, alpha)
! eigenvalues and eigenvectors
CALL atom_solve(hmat, integrals%utrans, wfn, ener, basis%nbas, integrals%nne, lmat)
! relativistic kinetic energy
cpot(1:m) = cspeed*cspeed/(2._dp*cspeed*cspeed-modpot(1:m))**2
cpot(1:m) = c_light_au*c_light_au/(2._dp*c_light_au*c_light_au-modpot(1:m))**2
cpot(1:m) = cpot(1:m)/basis%grid%rad2(1:m)
pvp = 0.0_dp
CALL numpot_matrix(pvp, cpot, basis, 0)
@@ -14,6 +14,7 @@
!> bohr : Conversion factor [Angstrom] -> [Bohr]
!> boltzmann : Boltzmann constant [J/K]
!> c_light : Speed of light in vacuum [m/s]
!> c_light_au : Speed of light in vacuum [a.u.]
!> e_charge : Elementary charge [C]
!> e_mass : Electron mass [kg]
!> e_gfactor : Electron g factor [ ]
@@ -71,7 +72,7 @@ MODULE physcon

PRIVATE

PUBLIC :: a_bohr, a_fine, a_mass, angstrom, atm, bar, bohr, boltzmann, &
PUBLIC :: a_bohr, a_fine, a_mass, angstrom, atm, bar, bohr, boltzmann, c_light_au, &
debye, e_charge, e_gfactor, e_mass, evolt, femtoseconds, h_bar, &
hertz, joule, kcalmol, kelvin, kjmol, massunit, mu_perm, n_avogadro, newton, &
p_mass, pascal, picoseconds, seconds, vibfac, wavenumbers
@@ -83,6 +84,8 @@ MODULE physcon
! Exact constants
! Speed of light in vacuum [m/s]
REAL(KIND=dp), PARAMETER :: c_light = 299792458.0_dp
! Speed of light in vacuum, in atomic units (=1/a_fine)
REAL(KIND=dp), PARAMETER :: c_light_au = 137.035999679_dp

! Magnetic constant or permeability of vacuum [N/A**2]
REAL(KIND=dp), PARAMETER :: mu_perm = 4.0_dp*pi*1.0E-7_dp
@@ -25,7 +25,7 @@ MODULE dkh_main
cp_fm_type
USE cp_gemm_interface, ONLY: cp_gemm
USE kinds, ONLY: dp
USE physcon, ONLY: a_fine
USE physcon, ONLY: c_light_au
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
#include "./base/base_uses.f90"
@@ -87,10 +87,10 @@ MODULE dkh_main
!> ev2 (:,:) DKH-even2 matrix in position space
!> ove (:,:) scratch
!> aux (:,:) scratch
!> velit velocity of light 137 a.u.
!> prea prefactor, 1/137^2
!> con2 prefactor, 2/137^2
!> con prefactor, 137^2
!> c_light_au velocity of light 137 a.u.
!> prea prefactor, 1/137^2
!> con2 prefactor, 2/137^2
!> con prefactor, 137^2
!> \author
!> Jens Thar, Barbara Kirchner: Uni Bonn (09/2006)
!> Markus Reiher: ETH Zurich (09/2006)
@@ -105,7 +105,6 @@ SUBROUTINE DKH_full_transformation(qs_env, matrix_s, matrix_v, matrix_t, matrix_
routineP = moduleN//':'//routineN

INTEGER :: handle, i
REAL(KIND=dp) :: velit
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: aa, e, ev0t, rr, tt
TYPE(cp_blacs_env_type), POINTER :: blacs_env
TYPE(cp_fm_struct_type), POINTER :: matrix_full
@@ -114,12 +113,6 @@ SUBROUTINE DKH_full_transformation(qs_env, matrix_s, matrix_v, matrix_t, matrix_

CALL timeset(routineN, handle)
NULLIFY (blacs_env)
!-----------------------------------------------------------------------
! Define velocity of light
!-----------------------------------------------------------------------

! velit = 137.0359895_dp ! 1/a_fine physcon
velit = 1._dp/a_fine

!-----------------------------------------------------------------------
! Construct the matrix structure
@@ -182,7 +175,7 @@ SUBROUTINE DKH_full_transformation(qs_env, matrix_s, matrix_v, matrix_t, matrix_
! Calculate kinetic part of Hamiltonian in T-basis
!-----------------------------------------------------------------------

CALL kintegral(n, ev0t, tt, e, velit)
CALL kintegral(n, ev0t, tt, e)

!-----------------------------------------------------------------------
! Calculate reverse transformation matrix revt
@@ -205,8 +198,8 @@ SUBROUTINE DKH_full_transformation(qs_env, matrix_s, matrix_v, matrix_t, matrix_
!-----------------------------------------------------------------------

DO i = 1, n
aa(i) = SQRT((velit*velit+e(i))/(2.0_dp*e(i)))
rr(i) = SQRT(velit*velit)/(velit*velit+e(i))
aa(i) = SQRT((c_light_au**2+e(i))/(2.0_dp*e(i)))
rr(i) = SQRT(c_light_au**2)/(c_light_au**2+e(i))
END DO

!-----------------------------------------------------------------------
@@ -329,20 +322,18 @@ END SUBROUTINE DKH_full_transformation
!> \param ev0t ...
!> \param tt ...
!> \param e ...
!> \param velit ...
! **************************************************************************************************
SUBROUTINE kintegral(n, ev0t, tt, e, velit)
SUBROUTINE kintegral(n, ev0t, tt, e)
INTEGER, INTENT(IN) :: n
REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: ev0t
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: tt
REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: e
REAL(KIND=dp), INTENT(IN) :: velit

INTEGER :: i
REAL(KIND=dp) :: con, con2, prea, ratio, tv1, tv2, tv3, &
tv4

prea = 1/(velit*velit)
prea = 1/(c_light_au**2)
con2 = prea+prea
con = 1.0_dp/prea

@@ -355,7 +346,7 @@ SUBROUTINE kintegral(n, ev0t, tt, e, velit)
! cancellation, otherwise calculate SQRT directly

ev0t(i) = tt(i)
ratio = tt(i)/velit
ratio = tt(i)/c_light_au
IF (ratio .LE. 0.02_dp) THEN
tv1 = tt(i)
tv2 = -tv1*tt(i)*prea*0.5_dp
@@ -1329,10 +1320,10 @@ SUBROUTINE DKH_atom_transformation(s, v, h, pVp, n, dkh_order)
! ev2 (:,:) DKH-even2 matrix in position space *
! ove (:,:) scratch *
! aux (:,:) scratch *
! velit velocity of light 137 a.u. *
! prea prefactor, 1/137^2 *
! con2 prefactor, 2/137^2 *
! con prefactor, 137^2 *
! c_light_au velocity of light 137 a.u. *
! prea prefactor, 1/137^2 *
! con2 prefactor, 2/137^2 *
! con prefactor, 137^2 *
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
@@ -1344,7 +1335,6 @@ SUBROUTINE DKH_atom_transformation(s, v, h, pVp, n, dkh_order)
routineP = moduleN//':'//routineN

INTEGER :: i, j, k
REAL(KIND=dp) :: velit
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: aa, e, ev0t, rr, tt
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: aux, eig, ev1, ev1t, ev2, ev2t, ev3, &
ev3t, ev4, ev4t, ove, pev1tp, pVpt, &
@@ -1353,12 +1343,6 @@ SUBROUTINE DKH_atom_transformation(s, v, h, pVp, n, dkh_order)
IF (dkh_order .LT. 0) RETURN

!CAW pp: p^2-values (in momentum-space), stored as matrix!!
!-----------------------------------------------------------------------
! Define velocity of light
!-----------------------------------------------------------------------
! velit = 137.0359895_dp

velit = 1._dp/a_fine

!-----------------------------------------------------------------------
! Allocate some matrices
@@ -1402,7 +1386,7 @@ SUBROUTINE DKH_atom_transformation(s, v, h, pVp, n, dkh_order)
! Calculate kinetic part of Hamiltonian in T-basis
!-----------------------------------------------------------------------

CALL kintegral_a(n, ev0t, tt, e, velit)
CALL kintegral_a(n, ev0t, tt, e)

!-----------------------------------------------------------------------
! Calculate reverse transformation matrix revt
@@ -1430,8 +1414,8 @@ SUBROUTINE DKH_atom_transformation(s, v, h, pVp, n, dkh_order)
!-----------------------------------------------------------------------

DO i = 1, n
aa(i) = SQRT((velit*velit+e(i))/(2.0_dp*e(i)))
rr(i) = SQRT(velit*velit)/(velit*velit+e(i))
aa(i) = SQRT((c_light_au**2+e(i))/(2.0_dp*e(i)))
rr(i) = SQRT(c_light_au**2)/(c_light_au**2+e(i))
END DO

!-----------------------------------------------------------------------
@@ -1546,17 +1530,15 @@ END SUBROUTINE dkh_atom_transformation
!> \param ev0t ...
!> \param tt ...
!> \param e ...
!> \param velit ...
! **************************************************************************************************
SUBROUTINE kintegral_a(n, ev0t, tt, e, velit)
SUBROUTINE kintegral_a(n, ev0t, tt, e)

!-----------------------------------------------------------------------

INTEGER, INTENT(IN) :: n
REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: ev0t
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: tt
REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: e
REAL(KIND=dp), INTENT(IN) :: velit

INTEGER :: i
REAL(KIND=dp) :: con, con2, prea, ratio, tv1, tv2, tv3, &
@@ -1571,15 +1553,15 @@ SUBROUTINE kintegral_a(n, ev0t, tt, e, velit)

! Calculate some constants

prea = 1/(velit*velit)
prea = 1/(c_light_au**2)
con2 = prea+prea
con = 1.0_dp/prea

! If T is sufficiently small, use series expansion to avoid
! cancellation, otherwise calculate SQRT directly

ev0t(i) = tt(i)
ratio = tt(i)/velit
ratio = tt(i)/c_light_au
IF (ratio .LE. 0.02_dp) THEN
tv1 = tt(i)
tv2 = -tv1*tt(i)*prea/2.0_dp
@@ -11,7 +11,7 @@
!> \par History
!> created [06.2015]
!> 05.2019: Added truncated coulomb operator (A. Bussy)
!> \author Doro/hea Golze
!> \author Dorothea Golze
! **************************************************************************************************
MODULE generic_os_integrals
USE ai_contraction_sphi, ONLY: ab_contract,&
@@ -69,9 +69,13 @@ MODULE hfx_libint_wrapper
#endif
END TYPE

!This is a safeguard/hack because 3-center ERIs are bugged on libint 2.5
#if(__LIBINT)
!Fortran 3-center integrals are bugged in libint 2.5. More specifically, the dimensions of the
!following array are mixed up. This should be fixed from libint 2.6 on. Then, one will be able
!to simply USE libint2_build_3eri, as currently commented above
TYPE(C_FUNPTR), DIMENSION(0:LIBINT2_MAX_AM,0:LIBINT2_MAX_AM,0:LIBINT2_MAX_AM_3eri), &
BIND(C) :: libint2_build_3eri
#endif

CONTAINS

@@ -1,4 +1,3 @@
!TEMP COMMENT
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 - 2019 CP2K developers group !
@@ -6195,7 +6195,7 @@ SUBROUTINE create_print_embed_pot_cube(section)
description="The stride (X,Y,Z) used to write the cube file "// &
"(larger values result in smaller cube files). You can provide 3 numbers (for X,Y,Z) or"// &
" 1 number valid for all components.", &
usage="STRIDE 2 2 2", n_var=-1, default_i_vals=(/2, 2, 2/), type_of_var=integer_t)
usage="STRIDE 1 1 1", n_var=-1, default_i_vals=(/1, 1, 1/), type_of_var=integer_t)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)

@@ -6249,7 +6249,7 @@ SUBROUTINE create_print_simple_grid(section)
description="The stride (X,Y,Z) used to write the cube file "// &
"(larger values result in smaller cube files). You can provide 3 numbers (for X,Y,Z) or"// &
" 1 number valid for all components.", &
usage="STRIDE 1 1 1", n_var=-1, default_i_vals=(/2, 2, 2/), type_of_var=integer_t)
usage="STRIDE 1 1 1", n_var=-1, default_i_vals=(/1, 1, 1/), type_of_var=integer_t)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)

@@ -8206,13 +8206,13 @@ SUBROUTINE create_xas_tdp_section(section)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, name="OPERATOR", variants=(/"OP"/), &
description="The type of operator used for exact exchange. The stendard "// &
description="The type of operator used for exact exchange. The standard "// &
"Coulomb operator cannot be used in periodic systems.", &
usage="OPERATOR {string}", &
repeats=.FALSE., &
default_i_val=do_potential_coulomb, &
enum_c_vals=s2a("COULOMB", "TRUNCATED", "SHORTRANGE"), &
enum_desc=s2a("Stendard Coulomb operator: 1/r", &
enum_desc=s2a("Standard Coulomb operator: 1/r", &
"Truncated Coulomb operator: 1/r if r < R_c, 0 otherwise ", &
"Short range: erfc(omega*r)/r"), &
enum_i_vals=(/do_potential_coulomb, do_potential_truncated, &
@@ -124,7 +124,7 @@ SUBROUTINE eri_3center(int_abc, la_min, la_max, npgfa, zeta, ra, &
zetk = zetc(kpgf)
c_start = (kpgf-1)*ncoset(lc_max)

!start with all the (c|ba) integrals (stendard order) and keep to lb >= la
!start with all the (c|ba) integrals (standard order) and keep to lb >= la
CALL set_params(lib, ra, rb, rc, la_max, lb_max, lc_max, zeti, zetj, zetk, op, &
params_out=params, r_cutoff=my_r_cutoff, omega=my_omega)

@@ -188,9 +188,9 @@ SUBROUTINE calculate_o3c_integrals(o3c, calculate_forces, matrix_p)
tvec(:, :) = 0.0_dp

rjk(1:3) = rik(1:3)-rij(1:3)
dij = SQRT(SUM(rij**2))
dik = SQRT(SUM(rik**2))
djk = SQRT(SUM(rjk**2))
dij = NORM2(rij)
dik = NORM2(rik)
djk = NORM2(rjk)

DO iset = 1, nseta
DO jset = 1, nsetb
@@ -67,7 +67,7 @@ MODULE xas_tdp_atom
USE orbital_transformation_matrices, ONLY: orbtramat
USE particle_methods, ONLY: get_particle_set
USE particle_types, ONLY: particle_type
USE physcon, ONLY: a_fine
USE physcon, ONLY: c_light_au
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_grid_atom, ONLY: allocate_grid_atom,&
@@ -3258,7 +3258,7 @@ SUBROUTINE integrate_soc_atoms(matrix_soc, xas_atom_env, qs_env)

INTEGER :: blk, handle, i, iat, ikind, ir, jat, &
maxso, na, nkind, nr, nset, nsgf
REAL(dp) :: cspeed, zeff
REAL(dp) :: zeff
REAL(dp), ALLOCATABLE, DIMENSION(:) :: Vr
REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: V
REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: intso
@@ -3281,7 +3281,6 @@ SUBROUTINE integrate_soc_atoms(matrix_soc, xas_atom_env, qs_env)
CALL get_qs_env(qs_env, nkind=nkind, qs_kind_set=qs_kind_set, matrix_s=matrix_s, &
particle_set=particle_set)
ALLOCATE (int_soc(nkind))
cspeed = 1.0_dp/a_fine

! Loop over the kinds to compute the integrals
DO ikind = 1, nkind
@@ -3302,7 +3301,7 @@ SUBROUTINE integrate_soc_atoms(matrix_soc, xas_atom_env, qs_env)
ALLOCATE (V(na, nr))
V = 0.0_dp
DO ir = 1, nr
CALL daxpy(na, Vr(ir)/(4.0_dp*cspeed**2-2.0_dp*Vr(ir)), grid%weight(1:na, ir), 1, &
CALL daxpy(na, Vr(ir)/(4.0_dp*c_light_au**2-2.0_dp*Vr(ir)), grid%weight(1:na, ir), 1, &
V(1:na, ir), 1)
END DO
DEALLOCATE (Vr)

0 comments on commit e01c96e

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