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