From 2ced3eabf80b1f31124662a38694cecfb826a4c8 Mon Sep 17 00:00:00 2001 From: Frederick Stein <43850145+fstein93@users.noreply.github.com> Date: Fri, 17 Dec 2021 18:40:49 +0100 Subject: [PATCH] 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 --- src/input_cp2k_mp2.F | 41 +++++-- src/libint_2c_3c.F | 2 + src/mp2_eri_gpw.F | 22 +++- src/mp2_setup.F | 4 + src/pw/pw_methods.F | 59 +++++++++++ .../regtest-mp2-lr-stress/H2O_mp2_mix_cl.inp | 100 ++++++++++++++++++ tests/QS/regtest-mp2-lr-stress/TEST_FILES | 1 + tests/QS/regtest-rs-dhft/H2O-B2P3LYP.inp | 98 +++++++++++++++++ tests/QS/regtest-rs-dhft/TEST_FILES | 1 + 9 files changed, 316 insertions(+), 12 deletions(-) create mode 100644 tests/QS/regtest-mp2-lr-stress/H2O_mp2_mix_cl.inp create mode 100644 tests/QS/regtest-rs-dhft/H2O-B2P3LYP.inp diff --git a/src/input_cp2k_mp2.F b/src/input_cp2k_mp2.F index 4b0479fe45..100710ed62 100644 --- a/src/input_cp2k_mp2.F +++ b/src/input_cp2k_mp2.F @@ -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:", & "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 ! ************************************************************************************************** diff --git a/src/libint_2c_3c.F b/src/libint_2c_3c.F index 05f5b10c5d..ffd83af9fc 100644 --- a/src/libint_2c_3c.F +++ b/src/libint_2c_3c.F @@ -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 diff --git a/src/mp2_eri_gpw.F b/src/mp2_eri_gpw.F index 446a41b385..30c9849218 100644 --- a/src/mp2_eri_gpw.F +++ b/src/mp2_eri_gpw.F @@ -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 diff --git a/src/mp2_setup.F b/src/mp2_setup.F index 2005e7bc44..59108cc728 100644 --- a/src/mp2_setup.F +++ b/src/mp2_setup.F @@ -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) diff --git a/src/pw/pw_methods.F b/src/pw/pw_methods.F index d4f9af7576..b309afd31c 100644 --- a/src/pw/pw_methods.F +++ b/src/pw/pw_methods.F @@ -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 diff --git a/tests/QS/regtest-mp2-lr-stress/H2O_mp2_mix_cl.inp b/tests/QS/regtest-mp2-lr-stress/H2O_mp2_mix_cl.inp new file mode 100644 index 0000000000..99997fd094 --- /dev/null +++ b/tests/QS/regtest-mp2-lr-stress/H2O_mp2_mix_cl.inp @@ -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 diff --git a/tests/QS/regtest-mp2-lr-stress/TEST_FILES b/tests/QS/regtest-mp2-lr-stress/TEST_FILES index 3fbe5d7ea3..a3c3538a27 100644 --- a/tests/QS/regtest-mp2-lr-stress/TEST_FILES +++ b/tests/QS/regtest-mp2-lr-stress/TEST_FILES @@ -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 diff --git a/tests/QS/regtest-rs-dhft/H2O-B2P3LYP.inp b/tests/QS/regtest-rs-dhft/H2O-B2P3LYP.inp new file mode 100644 index 0000000000..b1b7f2e05d --- /dev/null +++ b/tests/QS/regtest-rs-dhft/H2O-B2P3LYP.inp @@ -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 diff --git a/tests/QS/regtest-rs-dhft/TEST_FILES b/tests/QS/regtest-rs-dhft/TEST_FILES index 9fb6661be4..0675c407d8 100644 --- a/tests/QS/regtest-rs-dhft/TEST_FILES +++ b/tests/QS/regtest-rs-dhft/TEST_FILES @@ -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