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 Calculates 2-center integrals for different r12 operators comparing the Solid harmonic
10 : !> Gaussian integral scheme to the Obara-Saika (OS) scheme
11 : !> \author Dorothea Golze [05.2016]
12 : ! **************************************************************************************************
13 : MODULE shg_integrals_test
14 :
15 : USE basis_set_types, ONLY: allocate_gto_basis_set,&
16 : deallocate_gto_basis_set,&
17 : gto_basis_set_type,&
18 : init_orb_basis_set,&
19 : read_gto_basis_set
20 : USE constants_operator, ONLY: operator_coulomb,&
21 : operator_gauss,&
22 : operator_verf,&
23 : operator_verfc,&
24 : operator_vgauss
25 : USE cp_log_handling, ONLY: cp_to_string
26 : USE generic_os_integrals, ONLY: int_operators_r12_ab_os,&
27 : int_overlap_ab_os,&
28 : int_overlap_aba_os,&
29 : int_overlap_abb_os,&
30 : int_ra2m_ab_os
31 : USE generic_shg_integrals, ONLY: int_operators_r12_ab_shg,&
32 : int_overlap_ab_shg,&
33 : int_overlap_aba_shg,&
34 : int_overlap_abb_shg,&
35 : int_ra2m_ab_shg
36 : USE generic_shg_integrals_init, ONLY: contraction_matrix_shg,&
37 : contraction_matrix_shg_mix,&
38 : contraction_matrix_shg_rx2m,&
39 : get_clebsch_gordon_coefficients
40 : USE input_cp2k_subsys, ONLY: create_basis_section
41 : USE input_keyword_types, ONLY: keyword_create,&
42 : keyword_release,&
43 : keyword_type
44 : USE input_section_types, ONLY: &
45 : section_add_keyword, section_add_subsection, section_create, section_release, &
46 : section_type, section_vals_get, section_vals_get_subs_vals, section_vals_type, &
47 : section_vals_val_get
48 : USE input_val_types, ONLY: real_t
49 : USE kinds, ONLY: default_string_length,&
50 : dp
51 : USE orbital_pointers, ONLY: init_orbital_pointers
52 : USE orbital_transformation_matrices, ONLY: init_spherical_harmonics
53 : #include "./base/base_uses.f90"
54 :
55 : IMPLICIT NONE
56 :
57 : PRIVATE
58 :
59 : ! **************************************************************************************************
60 :
61 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'shg_integrals_test'
62 :
63 : PUBLIC :: create_shg_integrals_test_section, shg_integrals_perf_acc_test
64 :
65 : CONTAINS
66 :
67 : ! **************************************************************************************************
68 : !> \brief Create input section for unit testing
69 : !> \param section ...
70 : ! **************************************************************************************************
71 1386 : SUBROUTINE create_shg_integrals_test_section(section)
72 : TYPE(section_type), INTENT(INOUT), POINTER :: section
73 :
74 : TYPE(keyword_type), POINTER :: keyword
75 : TYPE(section_type), POINTER :: subsection
76 :
77 1386 : NULLIFY (keyword, subsection)
78 :
79 1386 : CPASSERT(.NOT. ASSOCIATED(section))
80 : CALL section_create(section, __LOCATION__, name="SHG_INTEGRALS_TEST", &
81 : description="Parameters for testing the SHG 2-center integrals for "// &
82 : "different r12 operators. Test w.r.t. performance and accurarcy.", &
83 1386 : n_keywords=4, n_subsections=1)
84 :
85 1386 : CALL create_basis_section(subsection)
86 1386 : CALL section_add_subsection(section, subsection)
87 1386 : CALL section_release(subsection)
88 :
89 : CALL keyword_create(keyword, __LOCATION__, &
90 : name="_SECTION_PARAMETERS_", &
91 : description="Controls the activation the SHG integral test. ", &
92 : default_l_val=.FALSE., &
93 1386 : lone_keyword_l_val=.TRUE.)
94 1386 : CALL section_add_keyword(section, keyword)
95 1386 : CALL keyword_release(keyword)
96 :
97 : CALL keyword_create(keyword, __LOCATION__, name="ABC", &
98 : description="Specify the lengths of the cell vectors A, B, and C. ", &
99 : usage="ABC 10.000 10.000 10.000", unit_str="angstrom", &
100 1386 : n_var=3, type_of_var=real_t)
101 1386 : CALL section_add_keyword(section, keyword)
102 1386 : CALL keyword_release(keyword)
103 :
104 : CALL keyword_create(keyword, __LOCATION__, name="NAB_MIN", &
105 : description="Minimum number of atomic distances to consider. ", &
106 1386 : default_i_val=8)
107 1386 : CALL section_add_keyword(section, keyword)
108 1386 : CALL keyword_release(keyword)
109 :
110 : CALL keyword_create(keyword, __LOCATION__, name="NREP", &
111 : description="Number of repeated calculation of each integral. ", &
112 1386 : default_i_val=1)
113 1386 : CALL section_add_keyword(section, keyword)
114 1386 : CALL keyword_release(keyword)
115 :
116 : CALL keyword_create(keyword, __LOCATION__, name="CHECK_ACCURACY", &
117 : description="Causes abortion when SHG and OS integrals differ "// &
118 : "more what's given by ACCURACY_LEVEL.", &
119 1386 : default_l_val=.TRUE., lone_keyword_l_val=.TRUE.)
120 1386 : CALL section_add_keyword(section, keyword)
121 1386 : CALL keyword_release(keyword)
122 :
123 : CALL keyword_create(keyword, __LOCATION__, name="ACCURACY_LEVEL", &
124 : description="Level of accuracy for comparison of SHG and OS "// &
125 : "integrals.", &
126 1386 : default_r_val=1.0E-8_dp)
127 1386 : CALL section_add_keyword(section, keyword)
128 1386 : CALL keyword_release(keyword)
129 :
130 : CALL keyword_create(keyword, __LOCATION__, name="CALCULATE_DERIVATIVES", &
131 : description="Calculates also the derivatives of the integrals.", &
132 1386 : default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
133 1386 : CALL section_add_keyword(section, keyword)
134 1386 : CALL keyword_release(keyword)
135 :
136 : CALL keyword_create(keyword, __LOCATION__, name="TEST_OVERLAP", &
137 : description="Calculates the integrals (a|b).", &
138 1386 : default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
139 1386 : CALL section_add_keyword(section, keyword)
140 1386 : CALL keyword_release(keyword)
141 :
142 1386 : CALL keyword_release(keyword)
143 : CALL keyword_create(keyword, __LOCATION__, name="TEST_COULOMB", &
144 : description="Calculates the integrals (a|1/r12|b).", &
145 1386 : default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
146 1386 : CALL section_add_keyword(section, keyword)
147 1386 : CALL keyword_release(keyword)
148 :
149 : CALL keyword_create(keyword, __LOCATION__, name="TEST_VERF", &
150 : description="Calculates the integrals (a|erf(omega*r12)/r12|b).", &
151 1386 : default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
152 1386 : CALL section_add_keyword(section, keyword)
153 1386 : CALL keyword_release(keyword)
154 :
155 : CALL keyword_create(keyword, __LOCATION__, name="TEST_VERFC", &
156 : description="Calculates the integrals (a|erfc(omega*r12)/r12|b).", &
157 1386 : default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
158 1386 : CALL section_add_keyword(section, keyword)
159 1386 : CALL keyword_release(keyword)
160 :
161 : CALL keyword_create(keyword, __LOCATION__, name="TEST_VGAUSS", &
162 : description="Calculates the integrals (a|exp(omega*r12^2)/r12|b).", &
163 1386 : default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
164 1386 : CALL section_add_keyword(section, keyword)
165 1386 : CALL keyword_release(keyword)
166 :
167 : CALL keyword_create(keyword, __LOCATION__, name="TEST_GAUSS", &
168 : description="Calculates the integrals (a|exp(omega*r12^2)|b).", &
169 1386 : default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
170 1386 : CALL section_add_keyword(section, keyword)
171 1386 : CALL keyword_release(keyword)
172 :
173 : CALL keyword_create(keyword, __LOCATION__, name="TEST_RA2M", &
174 : description="Calculates the integrals (a|(r-Ra)^(2m)|b).", &
175 1386 : default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
176 1386 : CALL section_add_keyword(section, keyword)
177 1386 : CALL keyword_release(keyword)
178 :
179 : CALL keyword_create(keyword, __LOCATION__, name="M", &
180 : description="Exponent in integral (a|(r-Ra)^(2m)|b).", &
181 1386 : default_i_val=1)
182 1386 : CALL section_add_keyword(section, keyword)
183 1386 : CALL keyword_release(keyword)
184 :
185 : CALL keyword_create(keyword, __LOCATION__, name="TEST_OVERLAP_ABA", &
186 : description="Calculates the integrals (a|b|b).", &
187 1386 : default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
188 1386 : CALL section_add_keyword(section, keyword)
189 1386 : CALL keyword_release(keyword)
190 :
191 : CALL keyword_create(keyword, __LOCATION__, name="TEST_OVERLAP_ABB", &
192 : description="Calculates the integrals (a|b|b).", &
193 1386 : default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
194 1386 : CALL section_add_keyword(section, keyword)
195 1386 : CALL keyword_release(keyword)
196 :
197 1386 : END SUBROUTINE create_shg_integrals_test_section
198 :
199 : ! **************************************************************************************************
200 : !> \brief Unit test for performance and accuracy of the SHG integrals
201 : !> \param iw output unit
202 : !> \param shg_integrals_test_section ...
203 : ! **************************************************************************************************
204 4 : SUBROUTINE shg_integrals_perf_acc_test(iw, shg_integrals_test_section)
205 : INTEGER, INTENT(IN) :: iw
206 : TYPE(section_vals_type), INTENT(INOUT), POINTER :: shg_integrals_test_section
207 :
208 : CHARACTER(len=*), PARAMETER :: routineN = 'shg_integrals_perf_acc_test'
209 : REAL(KIND=dp), PARAMETER :: d_nab = 1.0E-06_dp
210 :
211 : CHARACTER(LEN=default_string_length) :: basis_type
212 : INTEGER :: count_ab, handle, iab, jab, kab, lamax, &
213 : lbmax, lcamax, lcbmax, lmax, nab, &
214 : nab_min, nab_xyz, nrep, nrep_bas
215 : LOGICAL :: acc_check, calc_derivatives, &
216 : test_overlap_aba, test_overlap_abb
217 : REAL(KIND=dp) :: acc_param
218 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rab
219 4 : REAL(KIND=dp), DIMENSION(:), POINTER :: cell_par
220 4 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scona_shg, sconb_shg
221 : TYPE(gto_basis_set_type), POINTER :: fba, fbb, oba, obb
222 : TYPE(section_vals_type), POINTER :: basis_section
223 :
224 4 : CALL timeset(routineN, handle)
225 4 : NULLIFY (oba, obb, fba, fbb, basis_section, cell_par)
226 4 : CALL section_vals_val_get(shg_integrals_test_section, "ABC", r_vals=cell_par)
227 4 : CALL section_vals_val_get(shg_integrals_test_section, "NAB_MIN", i_val=nab_min)
228 4 : CALL section_vals_val_get(shg_integrals_test_section, "NREP", i_val=nrep)
229 4 : CALL section_vals_val_get(shg_integrals_test_section, "CHECK_ACCURACY", l_val=acc_check)
230 4 : CALL section_vals_val_get(shg_integrals_test_section, "ACCURACY_LEVEL", r_val=acc_param)
231 4 : CALL section_vals_val_get(shg_integrals_test_section, "CALCULATE_DERIVATIVES", l_val=calc_derivatives)
232 4 : CALL section_vals_val_get(shg_integrals_test_section, "TEST_OVERLAP_ABA", l_val=test_overlap_aba)
233 4 : CALL section_vals_val_get(shg_integrals_test_section, "TEST_OVERLAP_ABB", l_val=test_overlap_abb)
234 :
235 : !*** Read the basis set information
236 4 : basis_section => section_vals_get_subs_vals(shg_integrals_test_section, "BASIS")
237 4 : CALL section_vals_get(basis_section, n_repetition=nrep_bas)
238 4 : IF (.NOT. (nrep_bas == 2 .OR. nrep_bas == 3)) THEN
239 : CALL cp_abort(__LOCATION__, &
240 0 : "Provide basis sets")
241 : END IF
242 4 : CALL allocate_gto_basis_set(oba)
243 4 : CALL read_gto_basis_set(TRIM("A"), basis_type, oba, basis_section, irep=1)
244 26 : lamax = MAXVAL(oba%lmax)
245 4 : CALL allocate_gto_basis_set(obb)
246 4 : CALL read_gto_basis_set(TRIM("B"), basis_type, obb, basis_section, irep=2)
247 36 : lbmax = MAXVAL(obb%lmax)
248 4 : lmax = MAX(lamax, lbmax)
249 4 : IF (test_overlap_aba) THEN
250 2 : CALL allocate_gto_basis_set(fba)
251 2 : CALL read_gto_basis_set(TRIM("CA"), basis_type, fba, basis_section, irep=3)
252 32 : lcamax = MAXVAL(fba%lmax)
253 2 : lmax = MAX(lamax + lcamax, lbmax)
254 : END IF
255 4 : IF (test_overlap_abb) THEN
256 2 : CALL allocate_gto_basis_set(fbb)
257 2 : CALL read_gto_basis_set(TRIM("CB"), basis_type, fbb, basis_section, irep=3)
258 32 : lcbmax = MAXVAL(fbb%lmax)
259 2 : lmax = MAX(lamax, lbmax + lcbmax)
260 : END IF
261 4 : IF (test_overlap_aba .AND. test_overlap_abb) THEN
262 2 : lmax = MAX(MAX(lamax + lcamax, lbmax), MAX(lamax, lbmax + lcbmax))
263 : END IF
264 : !*** Initialize basis set information
265 4 : CALL init_orbital_pointers(lmax + 1)
266 4 : CALL init_spherical_harmonics(lmax, output_unit=-100)
267 4 : oba%norm_type = 2
268 4 : CALL init_orb_basis_set(oba)
269 4 : obb%norm_type = 2
270 4 : CALL init_orb_basis_set(obb)
271 4 : IF (test_overlap_aba) THEN
272 2 : fba%norm_type = 2
273 2 : CALL init_orb_basis_set(fba)
274 : END IF
275 4 : IF (test_overlap_abb) THEN
276 2 : fbb%norm_type = 2
277 2 : CALL init_orb_basis_set(fbb)
278 : END IF
279 : ! if shg integrals are later actually used in the code, contraction_matrix_shg should be
280 : ! moved to init_orb_basis_set and scon_shg should become an element of gto_basis_set_type
281 4 : CALL contraction_matrix_shg(oba, scona_shg)
282 4 : CALL contraction_matrix_shg(obb, sconb_shg)
283 :
284 : !*** Create range of rab (atomic distances) to be tested
285 4 : nab_xyz = CEILING(REAL(nab_min, KIND=dp)**(1.0_dp/3.0_dp) - d_nab)
286 4 : nab = nab_xyz**3
287 :
288 12 : ALLOCATE (rab(3, nab))
289 4 : count_ab = 0
290 12 : DO iab = 1, nab_xyz
291 28 : DO jab = 1, nab_xyz
292 56 : DO kab = 1, nab_xyz
293 32 : count_ab = count_ab + 1
294 144 : rab(:, count_ab) = [iab*ABS(cell_par(1)), jab*ABS(cell_par(2)), kab*ABS(cell_par(3))]/nab_xyz
295 : END DO
296 : END DO
297 : END DO
298 :
299 : !*** Calculate the SHG integrals
300 :
301 : CALL test_shg_operator12_integrals(oba, obb, rab, nrep, scona_shg, sconb_shg, &
302 : shg_integrals_test_section, acc_check, &
303 4 : acc_param, calc_derivatives, iw)
304 :
305 : CALL test_shg_overlap_integrals(oba, obb, rab, nrep, scona_shg, sconb_shg, &
306 : shg_integrals_test_section, acc_check, &
307 4 : acc_param, calc_derivatives, iw)
308 : CALL test_shg_ra2m_integrals(oba, obb, rab, nrep, scona_shg, sconb_shg, &
309 : shg_integrals_test_section, acc_check, &
310 4 : acc_param, calc_derivatives, iw)
311 :
312 : CALL test_shg_overlap_aba_integrals(oba, obb, fba, fbb, rab, nrep, scona_shg, sconb_shg, &
313 : shg_integrals_test_section, acc_check, &
314 4 : acc_param, calc_derivatives, iw)
315 :
316 4 : DEALLOCATE (scona_shg, sconb_shg, rab)
317 4 : CALL deallocate_gto_basis_set(oba)
318 4 : CALL deallocate_gto_basis_set(obb)
319 4 : IF (test_overlap_aba) CALL deallocate_gto_basis_set(fba)
320 4 : IF (test_overlap_abb) CALL deallocate_gto_basis_set(fbb)
321 :
322 4 : CALL timestop(handle)
323 :
324 12 : END SUBROUTINE shg_integrals_perf_acc_test
325 :
326 : ! **************************************************************************************************
327 : !> \brief tests two-center integrals of the type [a|O(r12)|b]
328 : !> \param oba basis set on a
329 : !> \param obb basis set on b
330 : !> \param rab distance between a and b
331 : !> \param nrep ...
332 : !> \param scona_shg SHG contraction matrix for a
333 : !> \param sconb_shg SHG contraction matrix for b
334 : !> \param shg_integrals_test_section ...
335 : !> \param acc_check if accuracy is checked
336 : !> \param acc_param accuracy level, if deviation larger abort
337 : !> \param calc_derivatives ...
338 : !> \param iw ...
339 : ! **************************************************************************************************
340 4 : SUBROUTINE test_shg_operator12_integrals(oba, obb, rab, nrep, scona_shg, sconb_shg, &
341 : shg_integrals_test_section, acc_check, &
342 : acc_param, calc_derivatives, iw)
343 : TYPE(gto_basis_set_type), POINTER :: oba, obb
344 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: rab
345 : INTEGER, INTENT(IN) :: nrep
346 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scona_shg, sconb_shg
347 : TYPE(section_vals_type), INTENT(IN), POINTER :: shg_integrals_test_section
348 : LOGICAL, INTENT(IN) :: acc_check
349 : REAL(KIND=dp), INTENT(IN) :: acc_param
350 : LOGICAL, INTENT(IN) :: calc_derivatives
351 : INTEGER, INTENT(IN) :: iw
352 :
353 : INTEGER :: iab, irep, nab, nfa, nfb
354 : LOGICAL :: test_any, test_coulomb, test_gauss, &
355 : test_verf, test_verfc, test_vgauss
356 : REAL(KIND=dp) :: ddmax_coulomb, ddmax_gauss, ddmax_verf, ddmax_verfc, ddmax_vgauss, ddtemp, &
357 : dmax_coulomb, dmax_gauss, dmax_verf, dmax_verfc, dmax_vgauss, dtemp, omega
358 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: vab_os, vab_shg
359 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dvab_os, dvab_shg
360 :
361 4 : CALL section_vals_val_get(shg_integrals_test_section, "TEST_COULOMB", l_val=test_coulomb)
362 4 : CALL section_vals_val_get(shg_integrals_test_section, "TEST_VERF", l_val=test_verf)
363 4 : CALL section_vals_val_get(shg_integrals_test_section, "TEST_VERFC", l_val=test_verfc)
364 4 : CALL section_vals_val_get(shg_integrals_test_section, "TEST_VGAUSS", l_val=test_vgauss)
365 4 : CALL section_vals_val_get(shg_integrals_test_section, "TEST_GAUSS", l_val=test_gauss)
366 :
367 4 : test_any = (test_coulomb .OR. test_verf .OR. test_verfc .OR. test_vgauss .OR. test_gauss)
368 :
369 : IF (test_any) THEN
370 2 : nfa = oba%nsgf
371 2 : nfb = obb%nsgf
372 16 : ALLOCATE (vab_shg(nfa, nfb), dvab_shg(nfa, nfb, 3))
373 10 : ALLOCATE (vab_os(nfa, nfb), dvab_os(nfa, nfb, 3))
374 2 : omega = 2.3_dp
375 2 : dmax_coulomb = 0.0_dp
376 2 : ddmax_coulomb = 0.0_dp
377 2 : dmax_verf = 0.0_dp
378 2 : ddmax_verf = 0.0_dp
379 2 : dmax_verfc = 0.0_dp
380 2 : ddmax_verfc = 0.0_dp
381 2 : dmax_vgauss = 0.0_dp
382 2 : ddmax_vgauss = 0.0_dp
383 2 : dmax_gauss = 0.0_dp
384 2 : ddmax_gauss = 0.0_dp
385 :
386 2 : nab = SIZE(rab, 2)
387 6 : DO irep = 1, nrep
388 38 : DO iab = 1, nab
389 : !*** Coulomb: (a|1/r12|b)
390 32 : IF (test_coulomb) THEN
391 : CALL int_operators_r12_ab_shg(operator_coulomb, vab_shg, dvab_shg, rab(:, iab), &
392 : oba, obb, scona_shg, sconb_shg, &
393 32 : calculate_forces=calc_derivatives)
394 : CALL int_operators_r12_ab_os(operator_coulomb, vab_os, dvab_os, rab(:, iab), &
395 32 : oba, obb, calculate_forces=calc_derivatives)
396 32 : CALL calculate_deviation_ab(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp)
397 32 : dmax_coulomb = MAX(dmax_coulomb, dtemp)
398 32 : ddmax_coulomb = MAX(ddmax_coulomb, ddtemp)
399 : END IF
400 : !*** verf: (a|erf(omega*r12)/r12|b)
401 32 : IF (test_verf) THEN
402 : CALL int_operators_r12_ab_shg(operator_verf, vab_shg, dvab_shg, rab(:, iab), &
403 : oba, obb, scona_shg, sconb_shg, omega, &
404 32 : calc_derivatives)
405 : CALL int_operators_r12_ab_os(operator_verf, vab_os, dvab_os, rab(:, iab), &
406 32 : oba, obb, omega=omega, calculate_forces=calc_derivatives)
407 32 : CALL calculate_deviation_ab(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp)
408 32 : dmax_verf = MAX(dmax_verf, dtemp)
409 32 : ddmax_verf = MAX(ddmax_verf, ddtemp)
410 : END IF
411 : !*** verfc: (a|erfc(omega*r12)/r12|b)
412 32 : IF (test_verfc) THEN
413 : CALL int_operators_r12_ab_shg(operator_verfc, vab_shg, dvab_shg, rab(:, iab), &
414 : oba, obb, scona_shg, sconb_shg, omega, &
415 32 : calc_derivatives)
416 : CALL int_operators_r12_ab_os(operator_verfc, vab_os, dvab_os, rab(:, iab), &
417 32 : oba, obb, omega=omega, calculate_forces=calc_derivatives)
418 32 : CALL calculate_deviation_ab(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp)
419 32 : dmax_verfc = MAX(dmax_verfc, dtemp)
420 32 : ddmax_verfc = MAX(ddmax_verfc, ddtemp)
421 : END IF
422 : !*** vgauss: (a|exp(omega*r12^2)/r12|b)
423 32 : IF (test_vgauss) THEN
424 : CALL int_operators_r12_ab_shg(operator_vgauss, vab_shg, dvab_shg, rab(:, iab), &
425 : oba, obb, scona_shg, sconb_shg, omega, &
426 32 : calc_derivatives)
427 : CALL int_operators_r12_ab_os(operator_vgauss, vab_os, dvab_os, rab(:, iab), &
428 32 : oba, obb, omega=omega, calculate_forces=calc_derivatives)
429 32 : CALL calculate_deviation_ab(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp)
430 32 : dmax_vgauss = MAX(dmax_vgauss, dtemp)
431 32 : ddmax_vgauss = MAX(ddmax_vgauss, ddtemp)
432 : END IF
433 : !*** gauss: (a|exp(omega*r12^2)|b)
434 36 : IF (test_gauss) THEN
435 : CALL int_operators_r12_ab_shg(operator_gauss, vab_shg, dvab_shg, rab(:, iab), &
436 : oba, obb, scona_shg, sconb_shg, omega, &
437 32 : calc_derivatives)
438 : CALL int_operators_r12_ab_os(operator_gauss, vab_os, dvab_os, rab(:, iab), &
439 32 : oba, obb, omega=omega, calculate_forces=calc_derivatives)
440 32 : CALL calculate_deviation_ab(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp)
441 32 : dmax_gauss = MAX(dmax_gauss, dtemp)
442 32 : ddmax_gauss = MAX(ddmax_gauss, ddtemp)
443 : END IF
444 : END DO
445 : END DO
446 :
447 2 : IF (iw > 0) THEN
448 1 : WRITE (iw, FMT="(/,T2,A)") "TEST INFO FOR 2-CENTER SHG and OS INTEGRALS:"
449 1 : WRITE (iw, FMT="(T2,A)") "Maximal deviation between SHG and OS integrals and their derivatives"
450 1 : IF (test_coulomb) THEN
451 1 : WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|1/r12|b]", &
452 2 : dmax_coulomb, ddmax_coulomb
453 : END IF
454 1 : IF (test_verf) THEN
455 1 : WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|erf(omega*r12)/r12|b]", &
456 2 : dmax_verf, ddmax_verf
457 : END IF
458 1 : IF (test_verfc) THEN
459 1 : WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|erfc(omega*r12)/r12|b]", &
460 2 : dmax_verfc, ddmax_verfc
461 : END IF
462 1 : IF (test_vgauss) THEN
463 1 : WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|exp(-omega*r12^2)/r12|b]", &
464 2 : dmax_vgauss, ddmax_vgauss
465 : END IF
466 1 : IF (test_gauss) THEN
467 1 : WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|exp(-omega*r12^2)|b]", &
468 2 : dmax_gauss, ddmax_gauss
469 : END IF
470 :
471 1 : IF (acc_check) THEN
472 1 : IF ((dmax_coulomb >= acc_param) .OR. (ddmax_coulomb >= acc_param)) THEN
473 0 : CPABORT("[a|1/r12|b]: Dev. larger than"//cp_to_string(acc_param))
474 : END IF
475 1 : IF ((dmax_verf >= acc_param) .OR. (ddmax_verf >= acc_param)) THEN
476 0 : CPABORT("[a|erf(omega*r12)/r12|b]: Dev. larger than"//cp_to_string(acc_param))
477 : END IF
478 1 : IF ((dmax_verfc >= acc_param) .OR. (ddmax_verfc >= acc_param)) THEN
479 0 : CPABORT("[a|erfc(omega*r12)/r12|b]: Dev. larger than"//cp_to_string(acc_param))
480 : END IF
481 1 : IF ((dmax_vgauss >= acc_param) .OR. (ddmax_vgauss >= acc_param)) THEN
482 0 : CPABORT("[a|exp(-omega*r12^2)/r12|b]: Dev. larger than"//cp_to_string(acc_param))
483 : END IF
484 1 : IF ((dmax_gauss >= acc_param) .OR. (ddmax_gauss >= acc_param)) THEN
485 0 : CPABORT("[a|exp(-omega*r12^2)|b]: Dev. larger than"//cp_to_string(acc_param))
486 : END IF
487 : END IF
488 : END IF
489 2 : DEALLOCATE (vab_shg, vab_os, dvab_shg, dvab_os)
490 : END IF
491 :
492 4 : END SUBROUTINE test_shg_operator12_integrals
493 :
494 : ! **************************************************************************************************
495 : !> \brief tests two center overlap integrals [a|b]
496 : !> \param oba ...
497 : !> \param obb ...
498 : !> \param rab ...
499 : !> \param nrep ...
500 : !> \param scona_shg ...
501 : !> \param sconb_shg ...
502 : !> \param shg_integrals_test_section ...
503 : !> \param acc_check ...
504 : !> \param acc_param ...
505 : !> \param calc_derivatives ...
506 : !> \param iw ...
507 : ! **************************************************************************************************
508 4 : SUBROUTINE test_shg_overlap_integrals(oba, obb, rab, nrep, scona_shg, sconb_shg, &
509 : shg_integrals_test_section, acc_check, &
510 : acc_param, calc_derivatives, iw)
511 : TYPE(gto_basis_set_type), POINTER :: oba, obb
512 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: rab
513 : INTEGER, INTENT(IN) :: nrep
514 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: scona_shg, sconb_shg
515 : TYPE(section_vals_type), INTENT(IN), POINTER :: shg_integrals_test_section
516 : LOGICAL, INTENT(IN) :: acc_check
517 : REAL(KIND=dp), INTENT(IN) :: acc_param
518 : LOGICAL, INTENT(IN) :: calc_derivatives
519 : INTEGER, INTENT(IN) :: iw
520 :
521 : INTEGER :: iab, irep, nab, nfa, nfb
522 : LOGICAL :: test_overlap
523 : REAL(KIND=dp) :: ddmax_overlap, ddtemp, dmax_overlap, &
524 : dtemp, dummy
525 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: sab_os, sab_shg
526 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dsab_os, dsab_shg
527 :
528 : CALL section_vals_val_get(shg_integrals_test_section, "TEST_OVERLAP", &
529 4 : l_val=test_overlap)
530 4 : IF (test_overlap) THEN
531 : !effectively switch off screening; makes no sense for the tests
532 22 : oba%set_radius(:) = 1.0E+09_dp
533 32 : obb%set_radius(:) = 1.0E+09_dp
534 42 : oba%pgf_radius(:, :) = 1.0E+09_dp
535 62 : obb%pgf_radius(:, :) = 1.0E+09_dp
536 2 : nfa = oba%nsgf
537 2 : nfb = obb%nsgf
538 2 : dummy = 0.0_dp
539 2 : dmax_overlap = 0.0_dp
540 2 : ddmax_overlap = 0.0_dp
541 16 : ALLOCATE (sab_shg(nfa, nfb), dsab_shg(nfa, nfb, 3))
542 10 : ALLOCATE (sab_os(nfa, nfb), dsab_os(nfa, nfb, 3))
543 2 : nab = SIZE(rab, 2)
544 6 : DO irep = 1, nrep
545 38 : DO iab = 1, nab
546 : CALL int_overlap_ab_shg(sab_shg, dsab_shg, rab(:, iab), oba, obb, &
547 32 : scona_shg, sconb_shg, calc_derivatives)
548 : CALL int_overlap_ab_os(sab_os, dsab_os, rab(:, iab), oba, obb, &
549 32 : calc_derivatives, debug=.FALSE., dmax=dummy)
550 32 : CALL calculate_deviation_ab(sab_shg, sab_os, dsab_shg, dsab_os, dtemp, ddtemp)
551 32 : dmax_overlap = MAX(dmax_overlap, dtemp)
552 36 : ddmax_overlap = MAX(ddmax_overlap, ddtemp)
553 : END DO
554 : END DO
555 :
556 2 : IF (iw > 0) THEN
557 1 : WRITE (iw, FMT="(/,T2,A)") "TEST INFO FOR 2-CENTER OVERLAP SHG and OS INTEGRALS:"
558 1 : WRITE (iw, FMT="(T2,A)") "Maximal deviation between SHG and OS integrals and their derivatives"
559 1 : WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|b]", &
560 2 : dmax_overlap, ddmax_overlap
561 : END IF
562 2 : IF (acc_check) THEN
563 2 : IF ((dmax_overlap >= acc_param) .OR. (ddmax_overlap >= acc_param)) THEN
564 0 : CPABORT("[a|b]: Deviation larger than"//cp_to_string(acc_param))
565 : END IF
566 : END IF
567 2 : DEALLOCATE (sab_shg, sab_os, dsab_shg, dsab_os)
568 : END IF
569 :
570 4 : END SUBROUTINE test_shg_overlap_integrals
571 :
572 : ! **************************************************************************************************
573 : !> \brief tests two-center integrals of the type [a|(r-Ra)^(2m)|b]
574 : !> \param oba ...
575 : !> \param obb ...
576 : !> \param rab ...
577 : !> \param nrep ...
578 : !> \param scona_shg ...
579 : !> \param sconb_shg ...
580 : !> \param shg_integrals_test_section ...
581 : !> \param acc_check ...
582 : !> \param acc_param ...
583 : !> \param calc_derivatives ...
584 : !> \param iw ...
585 : ! **************************************************************************************************
586 4 : SUBROUTINE test_shg_ra2m_integrals(oba, obb, rab, nrep, scona_shg, sconb_shg, &
587 : shg_integrals_test_section, acc_check, &
588 : acc_param, calc_derivatives, iw)
589 : TYPE(gto_basis_set_type), POINTER :: oba, obb
590 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: rab
591 : INTEGER, INTENT(IN) :: nrep
592 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: scona_shg, sconb_shg
593 : TYPE(section_vals_type), INTENT(IN), POINTER :: shg_integrals_test_section
594 : LOGICAL, INTENT(IN) :: acc_check
595 : REAL(KIND=dp), INTENT(IN) :: acc_param
596 : LOGICAL, INTENT(IN) :: calc_derivatives
597 : INTEGER, INTENT(IN) :: iw
598 :
599 : INTEGER :: iab, irep, m, nab, nfa, nfb
600 : LOGICAL :: test_ra2m
601 : REAL(KIND=dp) :: ddmax, ddtemp, dmax, dtemp
602 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: vab_os, vab_shg
603 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dvab_os, dvab_shg
604 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: scon_ra2m
605 :
606 : CALL section_vals_val_get(shg_integrals_test_section, "TEST_RA2M", &
607 4 : l_val=test_ra2m)
608 4 : IF (test_ra2m) THEN
609 : CALL section_vals_val_get(shg_integrals_test_section, "M", &
610 2 : i_val=m)
611 2 : nfa = oba%nsgf
612 2 : nfb = obb%nsgf
613 2 : dmax = 0.0_dp
614 2 : ddmax = 0.0_dp
615 2 : CALL contraction_matrix_shg_rx2m(oba, m, scona_shg, scon_ra2m)
616 16 : ALLOCATE (vab_shg(nfa, nfb), dvab_shg(nfa, nfb, 3))
617 10 : ALLOCATE (vab_os(nfa, nfb), dvab_os(nfa, nfb, 3))
618 2 : nab = SIZE(rab, 2)
619 6 : DO irep = 1, nrep
620 38 : DO iab = 1, nab
621 : CALL int_ra2m_ab_shg(vab_shg, dvab_shg, rab(:, iab), oba, obb, &
622 32 : scon_ra2m, sconb_shg, m, calc_derivatives)
623 32 : CALL int_ra2m_ab_os(vab_os, dvab_os, rab(:, iab), oba, obb, m, calc_derivatives)
624 32 : CALL calculate_deviation_ab(vab_shg, vab_os, dvab_shg, dvab_os, dtemp, ddtemp)
625 32 : dmax = MAX(dmax, dtemp)
626 36 : ddmax = MAX(ddmax, ddtemp)
627 : END DO
628 : END DO
629 2 : IF (iw > 0) THEN
630 1 : WRITE (iw, FMT="(/,T2,A)") "TEST INFO FOR 2-CENTER RA2m SHG and OS INTEGRALS:"
631 1 : WRITE (iw, FMT="(T2,A)") "Maximal deviation between SHG and OS integrals and their derivatives"
632 1 : WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|(r-Ra)^(2m)|b]", &
633 2 : dmax, ddmax
634 : END IF
635 2 : IF (acc_check) THEN
636 2 : IF ((dmax >= acc_param) .OR. (ddmax >= acc_param)) THEN
637 0 : CPABORT("[a|ra^(2m)|b]: Deviation larger than"//cp_to_string(acc_param))
638 : END IF
639 : END IF
640 2 : DEALLOCATE (scon_ra2m)
641 4 : DEALLOCATE (vab_shg, vab_os, dvab_shg, dvab_os)
642 : END IF
643 4 : END SUBROUTINE test_shg_ra2m_integrals
644 :
645 : ! **************************************************************************************************
646 : !> \brief test overlap integrals [a|b|a] and [a|b|b]
647 : !> \param oba ...
648 : !> \param obb ...
649 : !> \param fba ...
650 : !> \param fbb ...
651 : !> \param rab ...
652 : !> \param nrep ...
653 : !> \param scon_oba ...
654 : !> \param scon_obb ...
655 : !> \param shg_integrals_test_section ...
656 : !> \param acc_check ...
657 : !> \param acc_param ...
658 : !> \param calc_derivatives ...
659 : !> \param iw ...
660 : ! **************************************************************************************************
661 4 : SUBROUTINE test_shg_overlap_aba_integrals(oba, obb, fba, fbb, rab, nrep, scon_oba, scon_obb, &
662 : shg_integrals_test_section, acc_check, &
663 : acc_param, calc_derivatives, iw)
664 : TYPE(gto_basis_set_type), POINTER :: oba, obb, fba, fbb
665 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: rab
666 : INTEGER, INTENT(IN) :: nrep
667 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scon_oba, scon_obb
668 : TYPE(section_vals_type), INTENT(IN), POINTER :: shg_integrals_test_section
669 : LOGICAL, INTENT(IN) :: acc_check
670 : REAL(KIND=dp), INTENT(IN) :: acc_param
671 : LOGICAL, INTENT(IN) :: calc_derivatives
672 : INTEGER, INTENT(IN) :: iw
673 :
674 : INTEGER :: iab, irep, la_max, lb_max, lbb_max, &
675 : maxl_orb, maxl_ri, nab, nba, nbb, nfa, &
676 : nfb
677 4 : INTEGER, DIMENSION(:, :), POINTER :: ncg_none0
678 4 : INTEGER, DIMENSION(:, :, :), POINTER :: cg_none0_list, fba_index, fbb_index, &
679 4 : oba_index, obb_index
680 : LOGICAL :: test_overlap_aba, test_overlap_abb
681 : REAL(KIND=dp) :: ddmax_overlap_aba, ddmax_overlap_abb, &
682 : ddtemp, dmax_overlap_aba, &
683 : dmax_overlap_abb, dtemp, dummy
684 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: saba_os, saba_shg, sabb_os, sabb_shg
685 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dsaba_os, dsaba_shg, dsabb_os, dsabb_shg
686 4 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: cg_coeff
687 4 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: scona_mix, sconb_mix
688 :
689 : CALL section_vals_val_get(shg_integrals_test_section, "TEST_OVERLAP_ABA", &
690 4 : l_val=test_overlap_aba)
691 : CALL section_vals_val_get(shg_integrals_test_section, "TEST_OVERLAP_ABB", &
692 4 : l_val=test_overlap_abb)
693 4 : IF (test_overlap_aba .OR. test_overlap_abb) THEN
694 : !effectively switch off screening; makes no sense for the tests
695 4 : oba%set_radius(:) = 1.0E+09_dp
696 4 : obb%set_radius(:) = 1.0E+09_dp
697 18 : oba%pgf_radius(:, :) = 1.0E+09_dp
698 18 : obb%pgf_radius(:, :) = 1.0E+09_dp
699 2 : nba = oba%nsgf
700 2 : nbb = obb%nsgf
701 6 : maxl_orb = MAX(MAXVAL(oba%lmax), MAXVAL(obb%lmax))
702 4 : la_max = MAXVAL(oba%lmax)
703 4 : lb_max = MAXVAL(obb%lmax)
704 2 : IF (test_overlap_aba) THEN
705 32 : fba%set_radius(:) = 1.0E+09_dp
706 62 : fba%pgf_radius(:, :) = 1.0E+09_dp
707 2 : nfa = fba%nsgf
708 32 : maxl_ri = MAXVAL(fba%lmax) + 1 ! + 1 to avoid fail for l=0
709 20 : ALLOCATE (saba_shg(nba, nbb, nfa), dsaba_shg(nba, nbb, nfa, 3))
710 14 : ALLOCATE (saba_os(nba, nbb, nfa), dsaba_os(nba, nbb, nfa, 3))
711 2 : CALL contraction_matrix_shg_mix(oba, fba, oba_index, fba_index, scona_mix)
712 : END IF
713 2 : IF (test_overlap_abb) THEN
714 32 : fbb%set_radius(:) = 1.0E+09_dp
715 62 : fbb%pgf_radius(:, :) = 1.0E+09_dp
716 2 : nfb = fbb%nsgf
717 32 : maxl_ri = MAXVAL(fbb%lmax) + 1
718 36 : lbb_max = MAXVAL(obb%lmax) + MAXVAL(fbb%lmax)
719 20 : ALLOCATE (sabb_shg(nba, nbb, nfb), dsabb_shg(nba, nbb, nfb, 3))
720 14 : ALLOCATE (sabb_os(nba, nbb, nfb), dsabb_os(nba, nbb, nfb, 3))
721 2 : CALL contraction_matrix_shg_mix(obb, fbb, obb_index, fbb_index, sconb_mix)
722 : END IF
723 2 : dummy = 0.0_dp
724 2 : dmax_overlap_aba = 0.0_dp
725 2 : ddmax_overlap_aba = 0.0_dp
726 2 : dmax_overlap_abb = 0.0_dp
727 2 : ddmax_overlap_abb = 0.0_dp
728 2 : CALL get_clebsch_gordon_coefficients(cg_coeff, cg_none0_list, ncg_none0, maxl_orb, maxl_ri)
729 2 : nab = SIZE(rab, 2)
730 2 : IF (test_overlap_aba) THEN
731 4 : DO irep = 1, nrep
732 20 : DO iab = 1, nab
733 : CALL int_overlap_aba_shg(saba_shg, dsaba_shg, rab(:, iab), oba, obb, fba, &
734 : scon_obb, scona_mix, oba_index, fba_index, &
735 : cg_coeff, cg_none0_list, ncg_none0, &
736 16 : calc_derivatives)
737 : CALL int_overlap_aba_os(saba_os, dsaba_os, rab(:, iab), oba, obb, fba, &
738 16 : calc_derivatives, debug=.FALSE., dmax=dummy)
739 16 : CALL calculate_deviation_abx(saba_shg, saba_os, dsaba_shg, dsaba_os, dtemp, ddtemp)
740 16 : dmax_overlap_aba = MAX(dmax_overlap_aba, dtemp)
741 18 : ddmax_overlap_aba = MAX(ddmax_overlap_aba, ddtemp)
742 : END DO
743 : END DO
744 2 : DEALLOCATE (oba_index, fba_index, scona_mix)
745 2 : DEALLOCATE (saba_shg, saba_os, dsaba_shg, dsaba_os)
746 : END IF
747 2 : IF (test_overlap_abb) THEN
748 4 : DO irep = 1, nrep
749 20 : DO iab = 1, nab
750 : CALL int_overlap_abb_shg(sabb_shg, dsabb_shg, rab(:, iab), oba, obb, fbb, &
751 : scon_oba, sconb_mix, obb_index, fbb_index, &
752 : cg_coeff, cg_none0_list, ncg_none0, &
753 16 : calc_derivatives)
754 : CALL int_overlap_abb_os(sabb_os, dsabb_os, rab(:, iab), oba, obb, fbb, &
755 16 : calc_derivatives, debug=.FALSE., dmax=dummy)
756 16 : CALL calculate_deviation_abx(sabb_shg, sabb_os, dsabb_shg, dsabb_os, dtemp, ddtemp)
757 16 : dmax_overlap_abb = MAX(dmax_overlap_abb, dtemp)
758 18 : ddmax_overlap_abb = MAX(ddmax_overlap_abb, ddtemp)
759 : END DO
760 : END DO
761 2 : DEALLOCATE (obb_index, fbb_index, sconb_mix)
762 2 : DEALLOCATE (sabb_shg, sabb_os, dsabb_shg, dsabb_os)
763 : END IF
764 2 : IF (iw > 0) THEN
765 1 : WRITE (iw, FMT="(/,T2,A)") "TEST INFO [a|b|x] OVERLAP SHG and OS INTEGRALS:"
766 1 : WRITE (iw, FMT="(T2,A)") "Maximal deviation between SHG and OS integrals and their derivatives"
767 1 : IF (test_overlap_aba) THEN
768 1 : WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|b|a]", &
769 2 : dmax_overlap_aba, ddmax_overlap_aba
770 : END IF
771 1 : IF (test_overlap_abb) THEN
772 1 : WRITE (iw, FMT="(T2,A,T53,ES12.5,4X,ES12.5)") "SHG_INTEGRALS | [a|b|b]", &
773 2 : dmax_overlap_abb, ddmax_overlap_abb
774 : END IF
775 : END IF
776 2 : IF (acc_check) THEN
777 2 : IF ((dmax_overlap_aba >= acc_param) .OR. (ddmax_overlap_aba >= acc_param)) THEN
778 0 : CPABORT("[a|b|a]: Dev. larger than"//cp_to_string(acc_param))
779 : END IF
780 2 : IF ((dmax_overlap_abb >= acc_param) .OR. (ddmax_overlap_abb >= acc_param)) THEN
781 0 : CPABORT("[a|b|b]: Dev. larger than"//cp_to_string(acc_param))
782 : END IF
783 : END IF
784 2 : DEALLOCATE (cg_coeff, cg_none0_list, ncg_none0)
785 : END IF
786 :
787 8 : END SUBROUTINE test_shg_overlap_aba_integrals
788 :
789 : ! **************************************************************************************************
790 : !> \brief Calculation of the deviation between SHG and OS integrals
791 : !> \param vab_shg integral matrix obtained from the SHG scheme
792 : !> \param vab_os integral matrix obtained from the OS scheme
793 : !> \param dvab_shg derivative of the integrals, SHG
794 : !> \param dvab_os derivative of the integrals, OS
795 : !> \param dmax maximal deviation of vab matrices
796 : !> \param ddmax maximal deviation of dvab matrices
797 : ! **************************************************************************************************
798 224 : SUBROUTINE calculate_deviation_ab(vab_shg, vab_os, dvab_shg, dvab_os, dmax, ddmax)
799 :
800 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: vab_shg, vab_os
801 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: dvab_shg, dvab_os
802 : REAL(KIND=dp), INTENT(OUT) :: dmax, ddmax
803 :
804 : INTEGER :: i, j, k
805 : REAL(KIND=dp) :: diff
806 :
807 224 : dmax = 0.0_dp
808 224 : ddmax = 0.0_dp
809 :
810 : ! integrals vab
811 56448 : DO j = 1, SIZE(vab_shg, 2)
812 6747104 : DO i = 1, SIZE(vab_shg, 1)
813 6690656 : diff = ABS(vab_shg(i, j) - vab_os(i, j))
814 6746880 : dmax = MAX(dmax, diff)
815 : END DO
816 : END DO
817 :
818 : ! derivatives dvab
819 896 : DO k = 1, 3
820 169568 : DO j = 1, SIZE(dvab_shg, 2)
821 20241312 : DO i = 1, SIZE(dvab_shg, 1)
822 20071968 : diff = ABS(dvab_shg(i, j, k) - dvab_os(i, j, k))
823 20240640 : ddmax = MAX(ddmax, diff)
824 : END DO
825 : END DO
826 : END DO
827 :
828 224 : END SUBROUTINE calculate_deviation_ab
829 :
830 : ! **************************************************************************************************
831 : !> \brief Calculation of the deviation between SHG and OS integrals
832 : !> \param vab_shg integral matrix obtained from the SHG scheme
833 : !> \param vab_os integral matrix obtained from the OS scheme
834 : !> \param dvab_shg derivative of the integrals, SHG
835 : !> \param dvab_os derivative of the integrals, OS
836 : !> \param dmax maximal deviation of vab matrices
837 : !> \param ddmax maximal deviation of dvab matrices
838 : ! **************************************************************************************************
839 32 : SUBROUTINE calculate_deviation_abx(vab_shg, vab_os, dvab_shg, dvab_os, dmax, ddmax)
840 :
841 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: vab_shg, vab_os
842 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: dvab_shg, dvab_os
843 : REAL(KIND=dp), INTENT(OUT) :: dmax, ddmax
844 :
845 : INTEGER :: i, j, k, l
846 : REAL(KIND=dp) :: diff
847 :
848 32 : dmax = 0.0_dp
849 32 : ddmax = 0.0_dp
850 :
851 : ! integrals vab
852 8064 : DO k = 1, SIZE(vab_shg, 3)
853 112480 : DO j = 1, SIZE(vab_shg, 2)
854 634528 : DO i = 1, SIZE(vab_shg, 1)
855 522080 : diff = ABS(vab_shg(i, j, k) - vab_os(i, j, k))
856 626496 : dmax = MAX(dmax, diff)
857 : END DO
858 : END DO
859 : END DO
860 :
861 : ! derivatives dvab
862 128 : DO l = 1, 3
863 24224 : DO k = 1, SIZE(dvab_shg, 3)
864 337440 : DO j = 1, SIZE(dvab_shg, 2)
865 1903584 : DO i = 1, SIZE(dvab_shg, 1)
866 1566240 : diff = ABS(dvab_shg(i, j, k, l) - dvab_os(i, j, k, l))
867 1879488 : ddmax = MAX(ddmax, diff)
868 : END DO
869 : END DO
870 : END DO
871 : END DO
872 :
873 32 : END SUBROUTINE calculate_deviation_abx
874 :
875 : END MODULE shg_integrals_test
|
