Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Gaussian-product geometry and radial seeds shared by four-center evaluation and screening.
10 : ! **************************************************************************************************
11 : MODULE four_center_integral_preparation
12 : USE four_center_integral_types, ONLY: eri4_fundamentals_size,&
13 : eri4_pgf_image_type,&
14 : eri4_pgf_pair_type,&
15 : eri4_pgf_product_type
16 : USE input_constants, ONLY: do_potential_id
17 : USE integral_library_types, ONLY: coulomb_operator_type
18 : USE kinds, ONLY: dp
19 : USE mathconstants, ONLY: pi
20 : USE r12_fundamentals, ONLY: compute_r12_fundamentals
21 : #include "./base/base_uses.f90"
22 :
23 : IMPLICIT NONE
24 : PRIVATE
25 : PUBLIC :: prepare_eri4_pair, prepare_eri4_product, prepare_eri4_primitive
26 :
27 : CONTAINS
28 :
29 : ! **************************************************************************************************
30 : !> \brief Prepare one Gaussian pair, retaining the caller's image topology and screening metadata.
31 : !> \param ra first center
32 : !> \param rb second center in the selected image
33 : !> \param zeta first Gaussian exponent
34 : !> \param zetb second Gaussian exponent
35 : !> \param rab2 squared center distance, reusable across primitive pairs
36 : !> \param pair exponent data; image-list ownership is unchanged
37 : !> \param image product center and logarithmic Gaussian attenuation for this image
38 : ! **************************************************************************************************
39 249169194 : SUBROUTINE prepare_eri4_pair(ra, rb, zeta, zetb, rab2, pair, image)
40 : REAL(dp), INTENT(IN) :: ra(3), rb(3), zeta, zetb, rab2
41 : TYPE(eri4_pgf_pair_type), INTENT(INOUT) :: pair
42 : TYPE(eri4_pgf_image_type), INTENT(INOUT) :: image
43 :
44 249169194 : pair%zeta = zeta
45 249169194 : pair%zetb = zetb
46 249169194 : pair%zetapzetb = zeta + zetb
47 249169194 : pair%ZetaInv = 1.0_dp/pair%zetapzetb
48 996676776 : image%ra = ra
49 996676776 : image%rb = rb
50 249169194 : image%rab2 = rab2
51 249169194 : image%S1234 = -zeta*zetb*pair%ZetaInv*rab2
52 996676776 : image%P = (zeta*ra + zetb*rb)*pair%ZetaInv
53 249169194 : END SUBROUTINE prepare_eri4_pair
54 :
55 : ! **************************************************************************************************
56 : !> \brief Prepare the image-independent quartet data once per combination of Gaussian pairs.
57 : !> \param pair_a bra Gaussian pair
58 : !> \param pair_b ket Gaussian pair
59 : !> \param m_max maximum radial order, including any derivative increment
60 : !> \param product reusable template; centers and radial seeds are filled for each selected image
61 : ! **************************************************************************************************
62 186930594 : SUBROUTINE prepare_eri4_product(pair_a, pair_b, m_max, product)
63 : TYPE(eri4_pgf_pair_type), INTENT(IN) :: pair_a, pair_b
64 : INTEGER, INTENT(IN) :: m_max
65 : TYPE(eri4_pgf_product_type), INTENT(INOUT) :: product
66 :
67 186930594 : CPASSERT(m_max < eri4_fundamentals_size)
68 :
69 186930594 : product%prepared_max_order = m_max
70 934652970 : product%exponents = [pair_a%zeta, pair_a%zetb, pair_b%zeta, pair_b%zetb]
71 186930594 : product%ZetaInv = pair_a%ZetaInv
72 186930594 : product%EtaInv = pair_b%ZetaInv
73 186930594 : product%ZetapEtaInv = pair_a%zetapzetb + pair_b%zetapzetb
74 186930594 : product%ZetapEtaInv = 1.0_dp/product%ZetapEtaInv
75 186930594 : product%Rho = pair_a%zetapzetb*pair_b%zetapzetb*product%ZetapEtaInv
76 186930594 : END SUBROUTINE prepare_eri4_product
77 :
78 : ! **************************************************************************************************
79 : !> \brief Complete a prepared quartet for one selected image, including normalized radial seeds.
80 : !> \param product image-independent quartet template
81 : !> \param image_a bra pair geometry
82 : !> \param image_b ket pair geometry after the caller's image translation
83 : !> \param attenuation exponential of the two pair logarithmic attenuations, reusable across translations
84 : !> \param pq2 squared product-center distance, already available from cutoff screening
85 : !> \param potential_parameter radial operator parameters
86 : !> \param params complete primitive in original center order, ready for backend permutation
87 : ! **************************************************************************************************
88 402358382 : SUBROUTINE prepare_eri4_primitive(product, image_a, image_b, attenuation, pq2, potential_parameter, params)
89 : TYPE(eri4_pgf_product_type), INTENT(IN) :: product
90 : TYPE(eri4_pgf_image_type), INTENT(IN) :: image_a, image_b
91 : REAL(dp), INTENT(IN) :: attenuation, pq2
92 : TYPE(coulomb_operator_type), INTENT(IN) :: potential_parameter
93 : TYPE(eri4_pgf_product_type), INTENT(INOUT) :: params
94 :
95 : INTEGER :: i, m_max
96 : REAL(dp) :: den, fundamentals_work(eri4_fundamentals_size), G(3), normalization, num, &
97 : prefactor, ssss, T
98 :
99 402358382 : m_max = product%prepared_max_order
100 402358382 : params%prepared_max_order = m_max
101 2011791910 : params%exponents = product%exponents
102 402358382 : params%ZetaInv = product%ZetaInv
103 402358382 : params%EtaInv = product%EtaInv
104 402358382 : params%ZetapEtaInv = product%ZetapEtaInv
105 402358382 : params%Rho = product%Rho
106 1609433528 : params%centers(:, 1) = image_a%ra
107 1609433528 : params%centers(:, 2) = image_a%rb
108 1609433528 : params%centers(:, 3) = image_b%ra
109 1609433528 : params%centers(:, 4) = image_b%rb
110 1609433528 : params%P = image_a%P
111 1609433528 : params%Q = image_b%P
112 :
113 : ASSOCIATE (A => image_a%ra, B => image_a%rb, C => image_b%ra, D => image_b%rb, &
114 : zeta_a => product%exponents(1), zeta_b => product%exponents(2), &
115 : zeta_c => product%exponents(3), zeta_d => product%exponents(4))
116 1609433528 : params%W = ((zeta_a + zeta_b)*params%P + (zeta_c + zeta_d)*params%Q)*params%ZetapEtaInv
117 402358382 : T = params%Rho*pq2
118 804716764 : IF (potential_parameter%potential_type == do_potential_id) THEN
119 294732 : num = zeta_a*zeta_b
120 294732 : den = zeta_a + zeta_b
121 1178928 : ssss = -num/den*SUM((A - B)**2)
122 294732 : num = den*zeta_c
123 294732 : den = den + zeta_c
124 1178928 : ssss = ssss - num/den*SUM((params%P - C)**2)
125 1178928 : G = (zeta_a*A + zeta_b*B + zeta_c*C)/den
126 294732 : num = den*zeta_d
127 294732 : den = den + zeta_d
128 1178928 : ssss = ssss - num/den*SUM((G - D)**2)
129 6484104 : params%Fm(:) = EXP(ssss)
130 : ! The ID seed already contains the four-Gaussian attenuation.
131 294732 : normalization = 1.0_dp
132 294732 : IF (attenuation > EPSILON(0.0_dp)) normalization = 1.0_dp/attenuation
133 : ELSE
134 : CALL compute_r12_fundamentals(params%Fm, normalization, m_max, potential_parameter%potential_type, &
135 : params%Rho, T, fundamentals_work, &
136 : omega=potential_parameter%omega, &
137 : cutoff_radius=potential_parameter%cutoff_radius, &
138 : scale_coulomb=potential_parameter%scale_coulomb, &
139 : scale_longrange=potential_parameter%scale_longrange, &
140 402063650 : scale_gaussian=potential_parameter%scale_gaussian)
141 : END IF
142 : END ASSOCIATE
143 :
144 402358382 : prefactor = normalization*attenuation*SQRT((Pi*params%ZetapEtaInv)**3)
145 1467983387 : DO i = 1, m_max + 1
146 1467983387 : params%Fm(i) = params%Fm(i)*prefactor
147 : END DO
148 402358382 : END SUBROUTINE prepare_eri4_primitive
149 :
150 : END MODULE four_center_integral_preparation
|
