Add mixed potential to MP2 (#1819)

* Add mixed coulomb-longranged potential to MP2

* Add tests for mixed cl potential in MP2

* Accelerate tests/QS/regtest-rs-dhft/H2O-B2P3LYP.inp
This commit is contained in:
Frederick Stein 2021-12-17 18:40:49 +01:00 committed by GitHub
parent 9347a4eb7d
commit 2ced3eabf8
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
9 changed files with 316 additions and 12 deletions

View file

@ -26,12 +26,13 @@ MODULE input_cp2k_mp2
USE cp_units, ONLY: cp_unit_to_cp2k
USE input_constants, ONLY: &
do_eri_gpw, do_eri_mme, do_eri_os, do_potential_coulomb, do_potential_id, &
do_potential_long, do_potential_short, do_potential_truncated, do_potential_tshpsc, &
eri_default, gaussian, gw_no_print_exx, gw_pade_approx, gw_print_exx, gw_read_exx, &
gw_skip_for_regtest, gw_two_pole_model, kp_weights_W_auto, kp_weights_W_tailored, &
kp_weights_W_uniform, mp2_method_direct, mp2_method_gpw, mp2_method_none, numerical, &
ri_default, ri_rpa_g0w0_crossing_bisection, ri_rpa_g0w0_crossing_newton, &
ri_rpa_g0w0_crossing_z_shot, wfc_mm_style_gemm, wfc_mm_style_syrk
do_potential_long, do_potential_mix_cl, do_potential_short, do_potential_truncated, &
do_potential_tshpsc, eri_default, gaussian, gw_no_print_exx, gw_pade_approx, gw_print_exx, &
gw_read_exx, gw_skip_for_regtest, gw_two_pole_model, kp_weights_W_auto, &
kp_weights_W_tailored, kp_weights_W_uniform, mp2_method_direct, mp2_method_gpw, &
mp2_method_none, numerical, ri_default, ri_rpa_g0w0_crossing_bisection, &
ri_rpa_g0w0_crossing_newton, ri_rpa_g0w0_crossing_z_shot, wfc_mm_style_gemm, &
wfc_mm_style_syrk
USE input_cp2k_hfx, ONLY: create_hfx_section
USE input_cp2k_kpoints, ONLY: create_kpoint_set_section
USE input_keyword_types, ONLY: keyword_create,&
@ -1371,12 +1372,13 @@ CONTAINS
description="Which interaction potential should be used "// &
"(Coulomb, TShPSC operator).", &
usage="POTENTIAL_TYPE TSHPSC", &
enum_c_vals=s2a("COULOMB", "TShPSC", "LONGRANGE", "SHORTRANGE", "TRUNCATED"), &
enum_c_vals=s2a("COULOMB", "TShPSC", "LONGRANGE", "SHORTRANGE", "TRUNCATED", "MIX_CL"), &
enum_i_vals=(/do_potential_coulomb, &
do_potential_TShPSC, &
do_potential_long, &
do_potential_short, &
do_potential_truncated/), &
do_potential_truncated, &
do_potential_mix_cl/), &
enum_desc=s2a("Coulomb potential: 1/r", &
"TShPSC:<ul><li>1/x - s/Rc for x &le; Rc</li>"// &
"<li>(1 - s)/Rc - (x - Rc)/Rc^2 + (x - Rc)^2/Rc^3 - "// &
@ -1386,7 +1388,8 @@ CONTAINS
"<li>0 for x &gt; n*Rc</li></ul>", &
"Longrange Coulomb potential: erf(w*r)/r)", &
"Shortrange Coulomb potential: erfc(w*r)/r", &
"Truncated Coulomb potential"), &
"Truncated Coulomb potential", &
"Mixed Coulomb/Longrange Coulomb potential"), &
default_i_val=do_potential_coulomb)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
@ -1422,6 +1425,26 @@ CONTAINS
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create( &
keyword, __LOCATION__, &
name="SCALE_COULOMB", &
description="Scaling factor of (truncated) Coulomb potential in mixed (truncated) Coulomb/Longrange potential."// &
"Only valid when mixed potential is requested.", &
usage="OMEGA 0.5", type_of_var=real_t, &
default_r_val=1.0_dp)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create( &
keyword, __LOCATION__, &
name="SCALE_LONGRANGE", &
description="Scaling factor of longrange Coulomb potential in mixed (truncated) Coulomb/Longrange potential."// &
"Only valid when mixed potential is requested.", &
usage="OMEGA 0.5", type_of_var=real_t, &
default_r_val=1.0_dp)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
END SUBROUTINE create_mp2_potential
! **************************************************************************************************

View file

@ -68,6 +68,8 @@ MODULE libint_2c_3c
REAL(dp) :: omega !SR: erfc(w*r)/r
REAL(dp) :: cutoff_radius !TC cutoff/effective SR range
CHARACTER(default_path_length) :: filename
REAL(dp) :: scale_coulomb ! Only, for WFC methods
REAL(dp) :: scale_longrange ! Only, for WFC methods
END TYPE
CONTAINS

View file

@ -25,6 +25,7 @@ MODULE mp2_eri_gpw
USE input_constants, ONLY: do_potential_coulomb,&
do_potential_id,&
do_potential_long,&
do_potential_mix_cl,&
do_potential_short,&
do_potential_truncated
USE kinds, ONLY: dp
@ -39,8 +40,8 @@ MODULE mp2_eri_gpw
pw_env_release,&
pw_env_type
USE pw_methods, ONLY: &
pw_compl_gauss_damp, pw_copy, pw_derive, pw_gauss_damp, pw_integral_ab, pw_scale, &
pw_transfer, pw_truncated, pw_zero
pw_compl_gauss_damp, pw_copy, pw_derive, pw_gauss_damp, pw_gauss_damp_mix, pw_integral_ab, &
pw_scale, pw_transfer, pw_truncated, pw_zero
USE pw_poisson_methods, ONLY: pw_poisson_solve
USE pw_poisson_types, ONLY: pw_poisson_type
USE pw_pool_types, ONLY: pw_pool_create_pw,&
@ -832,7 +833,8 @@ CONTAINS
TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
INTEGER :: i
REAL(KIND=dp) :: omega_2, rchalf, tmp
REAL(KIND=dp) :: omega_2, rchalf, scale_coul, scale_long, &
tmp
IF (.NOT. (pw%in_space == RECIPROCALSPACE .AND. pw%in_use == COMPLEXDATA1D)) &
CPABORT("pw in wrong space or wrong data type")
@ -857,6 +859,17 @@ CONTAINS
pw%cc(i) = pw%cc(i)*(0.5_dp*tmp - tmp**2/12.0_dp)
END IF
END DO
!$OMP END PARALLEL DO
CASE (do_potential_mix_cl)
omega_2 = 1.0_dp/(2.0_dp*potential_parameter%omega)**2
scale_coul = potential_parameter%scale_coulomb
scale_long = potential_parameter%scale_longrange
!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,tmp) &
!$OMP SHARED(pw,omega_2,scale_long,scale_coul)
DO i = 1, SIZE(pw%cc)
tmp = omega_2*pw%pw_grid%gsq(i)
pw%cc(i) = pw%cc(i)*(1.0_dp + scale_long*tmp*EXP(-tmp)/(scale_coul + scale_long*EXP(-tmp)))
END DO
!$OMP END PARALLEL DO
CASE (do_potential_truncated)
rchalf = 0.5_dp*potential_parameter%cutoff_radius
@ -1263,6 +1276,9 @@ CONTAINS
IF (my_potential_type == do_potential_long) CALL pw_gauss_damp(pot_g%pw, potential_parameter%omega)
IF (my_potential_type == do_potential_short) CALL pw_compl_gauss_damp(pot_g%pw, potential_parameter%omega)
IF (my_potential_type == do_potential_truncated) CALL pw_truncated(pot_g%pw, potential_parameter%cutoff_radius)
IF (my_potential_type == do_potential_mix_cl) CALL pw_gauss_damp_mix(pot_g%pw, potential_parameter%omega, &
potential_parameter%scale_coulomb, &
potential_parameter%scale_longrange)
IF (my_transfer) CALL pw_transfer(pot_g%pw, pot_r%pw)
ELSE
! If we use an overlap metric, make sure to use the correct potential=density on output

View file

@ -326,6 +326,10 @@ CONTAINS
c_val=mp2_env%potential_parameter%filename)
CALL section_vals_val_get(mp2_section, "INTEGRALS%INTERACTION_POTENTIAL%OMEGA", &
r_val=mp2_env%potential_parameter%omega)
CALL section_vals_val_get(mp2_section, "INTEGRALS%INTERACTION_POTENTIAL%SCALE_COULOMB", &
r_val=mp2_env%potential_parameter%scale_coulomb)
CALL section_vals_val_get(mp2_section, "INTEGRALS%INTERACTION_POTENTIAL%SCALE_LONGRANGE", &
r_val=mp2_env%potential_parameter%scale_longrange)
NULLIFY (mp2_env%eri_mme_param)
ALLOCATE (mp2_env%eri_mme_param)

View file

@ -66,6 +66,7 @@ MODULE pw_methods
PUBLIC :: pw_zero, pw_structure_factor, pw_smoothing
PUBLIC :: pw_copy, pw_axpy, pw_transfer, pw_scale
PUBLIC :: pw_gauss_damp, pw_compl_gauss_damp, pw_derive, pw_dr2, pw_write, pw_multiply
PUBLIC :: pw_gauss_damp_mix
PUBLIC :: pw_integral_ab, pw_integral_a2b
PUBLIC :: pw_dr2_gg, pw_integrate_function
PUBLIC :: pw_set, pw_truncated
@ -483,6 +484,64 @@ CONTAINS
END SUBROUTINE pw_compl_gauss_damp
! **************************************************************************************************
!> \brief Multiply all data points with a Gaussian damping factor and mixes it with the original function
!> Needed for mixed longrange/Coulomb potential
!> V(\vec r)=(a+b*erf(omega*r))/r
!> V(\vec g)=\frac{4*\pi}{g**2}*(a+b*exp(-g**2/omega**2))
!> \param pw ...
!> \param omega ...
!> \param scale_coul ...
!> \param scale_long ...
!> \par History
!> Frederick Stein (16-Dec-2021) created
!> \author Frederick Stein (16-Dec-2021)
!> \note
!> Performs a Gaussian damping
!> PW has to be in RECIPROCALSPACE and data in use is COMPLEXDATA1D
! **************************************************************************************************
SUBROUTINE pw_gauss_damp_mix(pw, omega, scale_coul, scale_long)
TYPE(pw_type), INTENT(INOUT) :: pw
REAL(KIND=dp), INTENT(IN) :: omega, scale_coul, scale_long
CHARACTER(len=*), PARAMETER :: routineN = 'pw_gauss_damp_mix'
INTEGER :: cnt, handle, n_exp
REAL(KIND=dp) :: flop, omega_2
CALL timeset(routineN, handle)
CPASSERT(pw%ref_count > 0)
CPASSERT(omega >= 0)
flop = 0.0_dp
n_exp = 0
omega_2 = omega*omega
omega_2 = 0.25_dp/omega_2
IF (pw%in_space == RECIPROCALSPACE .AND. &
pw%in_use == COMPLEXDATA1D) THEN
cnt = SIZE(pw%cc)
!$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(cnt, pw, omega_2, scale_coul, scale_long)
pw%cc(:) = pw%cc(:)*(scale_coul + scale_long*EXP(-pw%pw_grid%gsq(:)*omega_2))
!$OMP END PARALLEL WORKSHARE
flop = flop + 4*cnt
n_exp = n_exp + cnt
ELSE
CPABORT("No suitable data field")
END IF
flop = flop*1.e-6_dp
CALL timestop(handle)
END SUBROUTINE pw_gauss_damp_mix
! **************************************************************************************************
!> \brief Multiply all data points with a complementary cosine
!> Needed for truncated Coulomb potential

View file

@ -0,0 +1,100 @@
&GLOBAL
PROJECT H2O_mp2_lr
PRINT_LEVEL LOW
RUN_TYPE ENERGY_FORCE
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
STRESS_TENSOR ANALYTICAL
&DFT
BASIS_SET_FILE_NAME HFX_BASIS
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF 100
REL_CUTOFF 30
&END MGRID
&QS
METHOD GPW
EPS_DEFAULT 1.0E-10
&END QS
&SCF
SCF_GUESS ATOMIC
EPS_SCF 1.0E-4
MAX_SCF 100
&END SCF
&XC
&XC_FUNCTIONAL NONE
&END XC_FUNCTIONAL
&HF
FRACTION 1.0
&INTERACTION_POTENTIAL
POTENTIAL_TYPE TRUNCATED
CUTOFF_RADIUS 1.99
T_C_G_DATA t_c_g.dat
&END
&SCREENING
EPS_SCHWARZ 1.0E-6
EPS_SCHWARZ_FORCES 1.0E-6
SCREEN_ON_INITIAL_P FALSE
&END SCREENING
&END HF
&WF_CORRELATION
&RI_MP2
BLOCK_SIZE 1
EPS_CANONICAL 0.0001
FREE_HFX_BUFFER .TRUE.
&CPHF
EPS_CONV 1.0E-4
MAX_ITER 10
&END
&END RI_MP2
&INTEGRALS
&INTERACTION_POTENTIAL
POTENTIAL_TYPE MIX_CL
OMEGA 0.2
SCALE_COULOMB 1.0
SCALE_LONGRANGE 1.5
&END INTERACTION_POTENTIAL
&WFC_GPW
CUTOFF 50
REL_CUTOFF 10
EPS_GRID 1.0E-6
EPS_FILTER 1.0E-6
&END WFC_GPW
&END INTEGRALS
MEMORY 200.
NUMBER_PROC 1
SCALE_S 1.0
SCALE_T 1.0
&END
&END XC
&END DFT
&PRINT
&STRESS_TENSOR
&END
&END
&SUBSYS
&CELL
ABC [angstrom] 5.0 5.0 5.0
&END CELL
&KIND H
BASIS_SET DZVP-GTH
BASIS_SET RI_AUX RI_DZVP-GTH
POTENTIAL GTH-PBE-q1
&END KIND
&KIND O
BASIS_SET DZVP-GTH
BASIS_SET RI_AUX RI_DZVP-GTH
POTENTIAL GTH-PBE-q6
&END KIND
&COORD
O 0.000000 0.000000 -0.211000
H 0.000000 -0.844000 0.495000
H 0.000000 0.744000 0.495000
&END
&TOPOLOGY
&CENTER_COORDINATES
&END
&END TOPOLOGY
&END SUBSYS
&END FORCE_EVAL

View file

@ -1,4 +1,5 @@
H2O_mp2_lr.inp 31 2e-05 -9.32626734433E+00
CH_mp2_lr.inp 31 2e-04 -5.76778050488E+00
H2O_mp2_lr_Coulomb_metric.inp 31 2e-05 -9.32473681187E+00
H2O_mp2_mix_cl.inp 31 2e-05 -8.27068864541E+00
#EOF

View file

@ -0,0 +1,98 @@
@SET MY_OMEGA 0.5
&GLOBAL
PROJECT H2O-srLDAlrMP2
PRINT_LEVEL MEDIUM
RUN_TYPE ENERGY
&TIMINGS
THRESHOLD 0.01
&END
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME HFX_BASIS
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF 100
REL_CUTOFF 20
&END MGRID
&POISSON
PERIODIC XYZ
POISSON_SOLVER WAVELET
&END POISSON
&QS
METHOD GPW
EPS_DEFAULT 1.0E-15
EPS_PGF_ORB 1.0E-30
&END QS
&SCF
SCF_GUESS RESTART
EPS_SCF 1.0E-7
MAX_SCF 100
! ADDED_MOS 15000 15000
&END SCF
&XC
&XC_FUNCTIONAL
&GGA_X_B88
SCALE 0.47
&END
&GGA_C_LYP
SCALE 0.73
&END
&END XC_FUNCTIONAL
&HF
FRACTION 0.53
&INTERACTION_POTENTIAL
POTENTIAL_TYPE TRUNCATED
CUTOFF_RADIUS 1.99
&END
&SCREENING
EPS_SCHWARZ 1.0E-6
SCREEN_ON_INITIAL_P FALSE
&END SCREENING
&END HF
&WF_CORRELATION
&RI_MP2
&END RI_MP2
&INTEGRALS
&INTERACTION_POTENTIAL
POTENTIAL_TYPE MIX_CL
OMEGA 0.2
SCALE_COULOMB 1.0
SCALE_LONGRANGE 1.5
&END INTERACTION_POTENTIAL
&WFC_GPW
CUTOFF 50
REL_CUTOFF 10
&END WFC_GPW
&END INTEGRALS
MEMORY 200.
NUMBER_PROC 1
SCALE_S 0.27
SCALE_T 0.27
&END
&END XC
&END DFT
&SUBSYS
&CELL
ABC [angstrom] 5.000 5.000 5.000
PERIODIC XYZ
&END CELL
&KIND H
BASIS_SET DZVP-GTH
BASIS_SET RI_AUX RI_DZVP-GTH
POTENTIAL GTH-HF-q1
&END KIND
&KIND O
BASIS_SET DZVP-GTH
BASIS_SET RI_AUX RI_DZVP-GTH
POTENTIAL GTH-HF-q6
&END KIND
&TOPOLOGY
COORD_FILE_NAME H2O_gas.xyz
COORD_FILE_FORMAT xyz
&CENTER_COORDINATES
&END
&END TOPOLOGY
&END SUBSYS
&END FORCE_EVAL

View file

@ -1,3 +1,4 @@
CH3-rsLDAlrMP2.inp 11 1e-8 -7.724767496356234
H2O-rsLDAlrMP2.inp 11 1e-8 -16.899255585036251
H2O-B2P3LYP.inp 11 1e-8 -17.030982292980237
#EOF