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 Calculation of the Hamiltonian integral matrix <a|H|b> for
10 : !> semi-empirical methods
11 : !> \author JGH
12 : ! **************************************************************************************************
13 : MODULE se_core_matrix
14 : USE atomic_kind_types, ONLY: atomic_kind_type,&
15 : get_atomic_kind,&
16 : get_atomic_kind_set
17 : USE cp_control_types, ONLY: dft_control_type
18 : USE cp_dbcsr_api, ONLY: &
19 : dbcsr_add, dbcsr_copy, dbcsr_deallocate_matrix, dbcsr_distribute, dbcsr_get_block_p, &
20 : dbcsr_p_type, dbcsr_replicate_all, dbcsr_set, dbcsr_sum_replicated, dbcsr_type
21 : USE cp_dbcsr_contrib, ONLY: dbcsr_get_block_diag
22 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
23 : USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
24 : USE cp_log_handling, ONLY: cp_get_default_logger,&
25 : cp_logger_type
26 : USE cp_output_handling, ONLY: cp_p_file,&
27 : cp_print_key_finished_output,&
28 : cp_print_key_should_output,&
29 : cp_print_key_unit_nr
30 : USE input_constants, ONLY: &
31 : do_method_am1, do_method_mndo, do_method_mndod, do_method_pdg, do_method_pm3, &
32 : do_method_pm6, do_method_pm6fm, do_method_pnnl, do_method_rm1
33 : USE input_section_types, ONLY: section_vals_val_get
34 : USE kinds, ONLY: dp
35 : USE message_passing, ONLY: mp_para_env_type
36 : USE particle_types, ONLY: particle_type
37 : USE physcon, ONLY: evolt
38 : USE qs_energy_types, ONLY: qs_energy_type
39 : USE qs_environment_types, ONLY: get_qs_env,&
40 : qs_environment_type
41 : USE qs_force_types, ONLY: qs_force_type
42 : USE qs_kind_types, ONLY: get_qs_kind,&
43 : qs_kind_type
44 : USE qs_ks_types, ONLY: qs_ks_env_type,&
45 : set_ks_env
46 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
47 : neighbor_list_iterate,&
48 : neighbor_list_iterator_create,&
49 : neighbor_list_iterator_p_type,&
50 : neighbor_list_iterator_release,&
51 : neighbor_list_set_p_type
52 : USE qs_overlap, ONLY: build_overlap_matrix
53 : USE qs_rho_types, ONLY: qs_rho_get,&
54 : qs_rho_type
55 : USE semi_empirical_int_arrays, ONLY: rij_threshold
56 : USE semi_empirical_types, ONLY: get_se_param,&
57 : semi_empirical_type
58 : USE semi_empirical_utils, ONLY: get_se_type
59 : USE virial_methods, ONLY: virial_pair_force
60 : USE virial_types, ONLY: virial_type
61 : #include "./base/base_uses.f90"
62 :
63 : IMPLICIT NONE
64 :
65 : PRIVATE
66 :
67 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'se_core_matrix'
68 :
69 : PUBLIC :: build_se_core_matrix
70 :
71 : CONTAINS
72 :
73 : ! **************************************************************************************************
74 : !> \brief ...
75 : !> \param qs_env ...
76 : !> \param para_env ...
77 : !> \param calculate_forces ...
78 : ! **************************************************************************************************
79 7382 : SUBROUTINE build_se_core_matrix(qs_env, para_env, calculate_forces)
80 :
81 : TYPE(qs_environment_type), POINTER :: qs_env
82 : TYPE(mp_para_env_type), POINTER :: para_env
83 : LOGICAL, INTENT(IN) :: calculate_forces
84 :
85 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_se_core_matrix'
86 :
87 : INTEGER :: after, atom_a, atom_b, handle, i, iatom, icol, icor, ikind, inode, irow, itype, &
88 : iw, j, jatom, jkind, natom, natorb_a, nkind, nr_a, nra, nrb
89 7382 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, nrt
90 : LOGICAL :: defined, found, omit_headers, use_virial
91 7382 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: se_defined
92 : REAL(KIND=dp) :: delta, dr, econst, eheat, eisol, kh, &
93 : udd, uff, upp, uss, ZPa, ZPb, ZSa, ZSb
94 7382 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ZPt, ZSt
95 7382 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: hmt, umt
96 : REAL(KIND=dp), DIMENSION(16) :: ha, hb, ua
97 : REAL(KIND=dp), DIMENSION(3) :: force_ab, rij
98 7382 : REAL(KIND=dp), DIMENSION(:), POINTER :: beta_a, sto_exponents_a
99 7382 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dsmat, h_block, h_blocka, pabmat, pamat, &
100 7382 : s_block
101 7382 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
102 : TYPE(cp_logger_type), POINTER :: logger
103 7382 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_h, matrix_p, matrix_s
104 : TYPE(dbcsr_type), POINTER :: diagmat_h, diagmat_p
105 : TYPE(dft_control_type), POINTER :: dft_control
106 : TYPE(neighbor_list_iterator_p_type), &
107 7382 : DIMENSION(:), POINTER :: nl_iterator
108 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
109 7382 : POINTER :: sab_orb
110 7382 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
111 : TYPE(qs_energy_type), POINTER :: energy
112 7382 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
113 7382 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
114 : TYPE(qs_ks_env_type), POINTER :: ks_env
115 : TYPE(qs_rho_type), POINTER :: rho
116 : TYPE(semi_empirical_type), POINTER :: se_kind_a
117 : TYPE(virial_type), POINTER :: virial
118 :
119 : ! REAL(KIND=dp), DIMENSION(3) :: R
120 :
121 7382 : CALL timeset(routineN, handle)
122 :
123 7382 : NULLIFY (logger, energy)
124 7382 : logger => cp_get_default_logger()
125 :
126 7382 : NULLIFY (rho, force, atomic_kind_set, qs_kind_set, sab_orb, &
127 7382 : diagmat_h, diagmat_p, particle_set, matrix_p, ks_env)
128 :
129 : CALL get_qs_env(qs_env, &
130 : matrix_s=matrix_s, &
131 : matrix_h=matrix_h, &
132 : ks_env=ks_env, &
133 : particle_set=particle_set, &
134 : atomic_kind_set=atomic_kind_set, &
135 : qs_kind_set=qs_kind_set, &
136 : dft_control=dft_control, &
137 : energy=energy, &
138 : force=force, &
139 : virial=virial, &
140 : rho=rho, &
141 7382 : sab_orb=sab_orb)
142 :
143 : ! calculate overlap matrix
144 7382 : IF (calculate_forces) THEN
145 : CALL build_overlap_matrix(ks_env, nderivative=1, matrix_s=matrix_s, &
146 : matrix_name="OVERLAP", &
147 : basis_type_a="ORB", &
148 : basis_type_b="ORB", &
149 3024 : sab_nl=sab_orb)
150 3024 : CALL set_ks_env(ks_env, matrix_s=matrix_s)
151 3024 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
152 : ELSE
153 : CALL build_overlap_matrix(ks_env, matrix_s=matrix_s, &
154 : matrix_name="OVERLAP", &
155 : basis_type_a="ORB", &
156 : basis_type_b="ORB", &
157 4358 : sab_nl=sab_orb)
158 4358 : CALL set_ks_env(ks_env, matrix_s=matrix_s)
159 4358 : use_virial = .FALSE.
160 : END IF
161 :
162 7382 : IF (calculate_forces) THEN
163 3024 : CALL qs_rho_get(rho, rho_ao=matrix_p)
164 :
165 3024 : IF (SIZE(matrix_p) == 2) THEN
166 140 : CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
167 : END IF
168 3024 : delta = dft_control%qs_control%se_control%delta
169 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
170 3024 : atom_of_kind=atom_of_kind)
171 3024 : ALLOCATE (diagmat_p)
172 3024 : CALL dbcsr_get_block_diag(matrix_p(1)%matrix, diagmat_p)
173 3024 : CALL dbcsr_replicate_all(diagmat_p)
174 : END IF
175 :
176 : ! Allocate the core Hamiltonian matrix
177 7382 : CALL dbcsr_allocate_matrix_set(matrix_h, 1)
178 7382 : ALLOCATE (matrix_h(1)%matrix)
179 7382 : CALL dbcsr_copy(matrix_h(1)%matrix, matrix_s(1)%matrix, "CORE HAMILTONIAN MATRIX")
180 7382 : CALL dbcsr_set(matrix_h(1)%matrix, 0.0_dp)
181 :
182 : ! Allocate a diagonal block matrix
183 7382 : ALLOCATE (diagmat_h)
184 7382 : CALL dbcsr_get_block_diag(matrix_s(1)%matrix, diagmat_h)
185 7382 : CALL dbcsr_set(diagmat_h, 0.0_dp)
186 7382 : CALL dbcsr_replicate_all(diagmat_h)
187 :
188 : ! kh might be set in qs_control
189 : itype = get_se_type(dft_control%qs_control%method_id)
190 7382 : kh = 0.5_dp
191 :
192 7382 : nkind = SIZE(atomic_kind_set)
193 :
194 22146 : ALLOCATE (se_defined(nkind))
195 22146 : ALLOCATE (hmt(16, nkind))
196 14764 : ALLOCATE (umt(16, nkind))
197 :
198 22146 : ALLOCATE (ZSt(nkind))
199 14764 : ALLOCATE (ZPt(nkind))
200 14764 : ALLOCATE (nrt(nkind))
201 :
202 7382 : econst = 0.0_dp
203 23774 : DO ikind = 1, nkind
204 16392 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom)
205 16392 : CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a)
206 : CALL get_se_param(se_kind_a, defined=defined, natorb=natorb_a, &
207 : beta=beta_a, uss=uss, upp=upp, udd=udd, uff=uff, eisol=eisol, eheat=eheat, &
208 16392 : nr=nr_a, sto_exponents=sto_exponents_a)
209 16392 : econst = econst - (eisol - eheat)*REAL(natom, dp)
210 16392 : se_defined(ikind) = (defined .AND. natorb_a >= 1)
211 16392 : hmt(1, ikind) = beta_a(0)
212 65568 : hmt(2:4, ikind) = beta_a(1)
213 98352 : hmt(5:9, ikind) = beta_a(2)
214 131136 : hmt(10:16, ikind) = beta_a(3)
215 16392 : umt(1, ikind) = uss
216 65568 : umt(2:4, ikind) = upp
217 98352 : umt(5:9, ikind) = udd
218 131136 : umt(10:16, ikind) = uff
219 :
220 16392 : ZSt(ikind) = sto_exponents_a(0)
221 16392 : ZPt(ikind) = sto_exponents_a(1)
222 56558 : nrt(ikind) = nr_a
223 :
224 : END DO
225 7382 : energy%core_self = econst
226 :
227 7382 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
228 432736 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
229 425354 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, inode=inode, r=rij)
230 425354 : IF (.NOT. se_defined(ikind)) CYCLE
231 425354 : IF (.NOT. se_defined(jkind)) CYCLE
232 7231018 : ha(1:16) = hmt(1:16, ikind)
233 7231018 : ua(1:16) = umt(1:16, ikind)
234 7231018 : hb(1:16) = hmt(1:16, jkind)
235 :
236 425354 : nra = nrt(ikind)
237 425354 : nrb = nrt(jkind)
238 425354 : ZSa = ZSt(ikind)
239 425354 : ZSb = ZSt(jkind)
240 425354 : ZPa = ZPt(ikind)
241 425354 : ZPb = ZPt(jkind)
242 :
243 425354 : IF (inode == 1) THEN
244 161534 : SELECT CASE (dft_control%qs_control%method_id)
245 : CASE (do_method_am1, do_method_rm1, do_method_mndo, do_method_pdg, &
246 : do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_mndod, do_method_pnnl)
247 80767 : NULLIFY (h_blocka)
248 80767 : CALL dbcsr_get_block_p(diagmat_h, iatom, iatom, h_blocka, found)
249 80767 : CPASSERT(ASSOCIATED(h_blocka))
250 161534 : IF (calculate_forces) THEN
251 34572 : CALL dbcsr_get_block_p(diagmat_p, iatom, iatom, pamat, found)
252 34572 : CPASSERT(ASSOCIATED(pamat))
253 : END IF
254 : END SELECT
255 : END IF
256 1701416 : dr = SUM(rij(:)**2)
257 425354 : IF (iatom == jatom .AND. dr < rij_threshold) THEN
258 :
259 28364 : SELECT CASE (dft_control%qs_control%method_id)
260 : CASE DEFAULT
261 0 : CPABORT("Unknown method for semi-empirical calculation")
262 : CASE (do_method_am1, do_method_rm1, do_method_mndo, do_method_pdg, &
263 : do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_mndod, do_method_pnnl)
264 115992 : DO i = 1, SIZE(h_blocka, 1)
265 115992 : h_blocka(i, i) = h_blocka(i, i) + ua(i)
266 : END DO
267 : END SELECT
268 :
269 : ELSE
270 396990 : IF (iatom <= jatom) THEN
271 191958 : irow = iatom
272 191958 : icol = jatom
273 : ELSE
274 205032 : irow = jatom
275 205032 : icol = iatom
276 : END IF
277 396990 : NULLIFY (h_block)
278 : CALL dbcsr_get_block_p(matrix_h(1)%matrix, &
279 396990 : irow, icol, h_block, found)
280 396990 : CPASSERT(ASSOCIATED(h_block))
281 : ! two-centre one-electron term
282 396990 : NULLIFY (s_block)
283 :
284 : CALL dbcsr_get_block_p(matrix_s(1)%matrix, &
285 396990 : irow, icol, s_block, found)
286 396990 : CPASSERT(ASSOCIATED(s_block))
287 396990 : IF (irow == iatom) THEN
288 808558 : DO i = 1, SIZE(h_block, 1)
289 2988620 : DO j = 1, SIZE(h_block, 2)
290 2796662 : h_block(i, j) = h_block(i, j) + kh*(ha(i) + hb(j))*s_block(i, j)
291 : END DO
292 : END DO
293 : ELSE
294 859632 : DO i = 1, SIZE(h_block, 1)
295 3201963 : DO j = 1, SIZE(h_block, 2)
296 2996931 : h_block(i, j) = h_block(i, j) + kh*(ha(j) + hb(i))*s_block(i, j)
297 : END DO
298 : END DO
299 : END IF
300 396990 : IF (calculate_forces) THEN
301 173955 : atom_a = atom_of_kind(iatom)
302 173955 : atom_b = atom_of_kind(jatom)
303 :
304 173955 : CALL dbcsr_get_block_p(matrix_p(1)%matrix, irow, icol, pabmat, found)
305 173955 : CPASSERT(ASSOCIATED(pabmat))
306 695820 : DO icor = 1, 3
307 521865 : force_ab(icor) = 0._dp
308 :
309 521865 : CALL dbcsr_get_block_p(matrix_s(icor + 1)%matrix, irow, icol, dsmat, found)
310 521865 : CPASSERT(ASSOCIATED(dsmat))
311 12004470 : dsmat = 2._dp*kh*dsmat*pabmat
312 1217685 : IF (irow == iatom) THEN
313 983133 : DO i = 1, SIZE(h_block, 1)
314 2966799 : DO j = 1, SIZE(h_block, 2)
315 2713941 : force_ab(icor) = force_ab(icor) + (ha(i) + hb(j))*dsmat(i, j)
316 : END DO
317 : END DO
318 : ELSE
319 1039803 : DO i = 1, SIZE(h_block, 1)
320 3157656 : DO j = 1, SIZE(h_block, 2)
321 2888649 : force_ab(icor) = force_ab(icor) + (ha(j) + hb(i))*dsmat(i, j)
322 : END DO
323 : END DO
324 : END IF
325 : END DO
326 : END IF
327 :
328 : END IF
329 :
330 432736 : IF (calculate_forces .AND. (iatom /= jatom .OR. dr > rij_threshold)) THEN
331 511099 : IF (irow == iatom) force_ab = -force_ab
332 : force(ikind)%all_potential(:, atom_a) = &
333 695820 : force(ikind)%all_potential(:, atom_a) - force_ab(:)
334 : force(jkind)%all_potential(:, atom_b) = &
335 695820 : force(jkind)%all_potential(:, atom_b) + force_ab(:)
336 173955 : IF (use_virial) THEN
337 0 : CALL virial_pair_force(virial%pv_virial, -1.0_dp, force_ab, rij)
338 : END IF
339 : END IF
340 :
341 : END DO
342 7382 : CALL neighbor_list_iterator_release(nl_iterator)
343 :
344 7382 : DEALLOCATE (se_defined, hmt, umt, ZSt, ZPt, nrt)
345 :
346 7382 : CALL dbcsr_sum_replicated(diagmat_h)
347 7382 : CALL dbcsr_distribute(diagmat_h)
348 7382 : CALL dbcsr_add(matrix_h(1)%matrix, diagmat_h, 1.0_dp, 1.0_dp)
349 7382 : CALL set_ks_env(ks_env, matrix_h=matrix_h)
350 :
351 7382 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
352 : qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN"), cp_p_file)) THEN
353 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN", &
354 1264 : extension=".Log")
355 1264 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
356 1264 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
357 1264 : after = MIN(MAX(after, 1), 16)
358 : CALL cp_dbcsr_write_sparse_matrix(matrix_h(1)%matrix, 4, after, qs_env, para_env, &
359 1264 : scale=evolt, output_unit=iw, omit_headers=omit_headers)
360 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
361 1264 : "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN")
362 : END IF
363 :
364 7382 : IF (calculate_forces) THEN
365 3024 : IF (SIZE(matrix_p) == 2) THEN
366 140 : CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
367 : END IF
368 3024 : DEALLOCATE (atom_of_kind)
369 3024 : CALL dbcsr_deallocate_matrix(diagmat_p)
370 : END IF
371 :
372 7382 : CALL dbcsr_deallocate_matrix(diagmat_h)
373 :
374 7382 : CALL timestop(handle)
375 :
376 14764 : END SUBROUTINE build_se_core_matrix
377 :
378 : ! **************************************************************************************************
379 : !> \brief ...
380 : !> \param R ...
381 : !> \param nra ...
382 : !> \param nrb ...
383 : !> \param ZSA ...
384 : !> \param ZSB ...
385 : !> \param ZPA ...
386 : !> \param ZPB ...
387 : !> \param S ...
388 : ! **************************************************************************************************
389 0 : SUBROUTINE makeS(R, nra, nrb, ZSA, ZSB, ZPA, ZPB, S)
390 :
391 : REAL(kind=dp), DIMENSION(3) :: R
392 : INTEGER :: nra, nrb
393 : REAL(kind=dp) :: ZSA, ZSB, ZPA, ZPB
394 : REAL(kind=dp), DIMENSION(4, 4) :: S
395 :
396 : INTEGER, DIMENSION(4, 4), PARAMETER :: &
397 : nc1 = RESHAPE([2, 4, 4, 6, 4, 3, 6, 7, 4, 6, 4, 8, 6, 7, 8, 5], [4, 4]), &
398 : nc2 = RESHAPE([4, 4, 8, 8, 6, 8, 6, 12, 8, 8, 12, 8, 10, 12, 14, 16], [4, 4]), &
399 : nc3 = RESHAPE([4, 6, 8, 10, 4, 8, 8, 12, 8, 6, 12, 14, 8, 12, 8, 16], [4, 4]), &
400 : nc4 = RESHAPE([4, 8, 11, 14, 8, 6, 12, 14, 11, 12, 10, 20, 14, 14, 20, 12], [4, 4]), &
401 : nc5 = RESHAPE([2, 4, 6, 8, 4, 4, 8, 8, 6, 8, 6, 12, 8, 8, 12, 8], [4, 4])
402 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: c1 = RESHAPE([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, &
403 : 1, 1, 1, 1, -1, 1, 2, 3, -1, -2, 1, 2, -2, -1, -3, 1, -3, -2, -1, -4, 0, -1, -2, 2, -1, 1,&
404 : -2, -1, 2, -2, 3, -3, 2, -1, -3, 6, 0, -1, -1, -2, 1, 0, -2, -4, -1, 2, -1, -3, 2, 4, 3, &
405 : -4, 0, 0, 0, -3, 0, 0, 1, -1, 0, 1, 0, 3, -3, -1, 3, 1, 0, 0, 0, -1, 0, 0, 1, 2, 0, -1, 0,&
406 : 3, 1, -2, -3, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, -1, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0 &
407 : , 0, 0, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
408 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
409 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
410 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
411 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
412 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
413 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
414 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: c2 = RESHAPE([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, &
415 : 1, 1, 1, 1, -1, 1, 1, 2, -2, -1, 1, 1, -3, -2, -1, 1, -4, -3, -2, -1, 1, -1, 1, 1, 1, 1, &
416 : -2, 1, 1, 1, 1, -3, 1, 1, 1, 1, -1, -1, -1, 2, 1, -1, -2, -2, 3, -2, -2, -3, 6, 2, -1, -3,&
417 : 0, 0, 1, -2, -2, -1, 1, 1, -3, 2, -1, 3, -4, -3, -2, -1, 0, 0, -1, -1, 1, 1, 1, -2, -1, -1&
418 : , 2, 3, -4, 2, 4, 3, 0, 0, -1, -2, 0, -1, 0, -2, 3, 2, -2, -1, 6, 2, -1, -3, 0, 0, -1, -1,&
419 : 0, 1, 0, 1, -1, -1, 1, -1, 1, -3, -1, 3, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 2, 0, -4, 2, 4, 3 &
420 : , 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 1, 1, -2, -3, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0&
421 : , -3, -1, 3, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 0, 1, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
422 : 0, 0, 0, 0, 0, -2, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0&
423 : , 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, &
424 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
425 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],&
426 : [4, 4, 20])
427 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: c3 = RESHAPE([-1, -1, -1, -1, -1, -1, -1, -1, -1, &
428 : -1, -1, -1, -1, -1, -1, -1, -1, -2, -3, -4, 1, -1, -2, -3, 1, 1, -1, -2, 2, 1, 1, -1, 1, 1&
429 : , 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 3, 1, 1, -1, -3, -6, -1, 1, 2, -2, 1, -2, 2, 1, -2, &
430 : 2, -3, 3, 0, 2, 3, 4, 0, 1, 2, 3, -1, -1, 1, 2, -2, -1, -3, 1, 0, 1, -1, -4, 0, 1, 1, 2, &
431 : -1, 1, 2, 4, 1, -2, 3, 3, 0, 0, 3, 6, 0, -1, -2, 2, -1, 0, -2, -1, 2, -2, 1, -3, 0, 0, 1, &
432 : -1, 0, -1, -1, 3, 1, 0, -1, 1, -1, -1, -1, -3, 0, 0, 0, 4, 0, 0, 0, -2, 0, 0, -2, -4, 0, 2&
433 : , 0, -3, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, -1, -2, 0, 1, 0, -3, 0, 0, 0, 0, 0, 0, 0, -3, 0, 0,&
434 : 1, -1, 0, 1, 0, 3, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, -1, 0, -1, 0, 1, 0, 0, 0, 0, 0, 0, 0,&
435 : 0, 0, 0, 0, 2, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, &
436 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0&
437 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
438 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
439 : , 0], [4, 4, 20])
440 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: c4 = RESHAPE([-1, -1, -1, -1, -1, -1, -1, -1, -1, &
441 : -1, -1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -2, -3,&
442 : 1, 1, -1, -2, 2, 1, 2, -1, 3, 2, 1, 3, -1, 1, 2, 3, -1, -1, 1, 2, -2, -1, -1, 1, -3, -2, &
443 : -1, -2, 0, 1, -1, -3, 1, -1, 1, 1, -1, 1, -1, 2, -3, 1, 2, -1, 0, -1, 2, 4, -1, 1, -1, -1 &
444 : , 2, -1, -1, -1, 4, -1, -1, -3, 0, 1, -1, -1, -1, 0, 1, 2, -1, -1, -1, -1, -1, -2, -1, 3, &
445 : 0, -1, 2, -1, 1, 0, -1, -2, -2, 1, 2, 2, 1, 2, -2, 1, 0, 0, -2, 4, 0, 0, -1, 1, 2, -1, 1, &
446 : -1, -4, 1, 1, 2, 0, 0, 1, -3, 0, 0, 1, -1, 1, 1, -1, -1, 3, -1, 1, -3, 0, 0, -1, 3, 0, 0, &
447 : -1, -2, -1, 1, 0, -1, 3, 2, -1, -1, 0, 0, 0, -3, 0, 0, 1, 2, 0, -1, 0, -1, -3, -2, -1, 1, &
448 : 0, 0, 0, 1, 0, 0, 0, -1, 0, 0, 0, 2, -1, -1, 2, 0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, 0, -1, 1&
449 : , 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
450 : 0, 2, 0, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,&
451 : 0, 0, 0, -1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 0, &
452 : 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0], [4, 4, 20])
453 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: c5 = RESHAPE([-1, -1, -1, -1, -1, -1, -1, -1, -1, &
454 : -1, -1, -1, -1, -1, -1, -1, 1, -1, -2, -3, 1, 1, -1, -2, 2, 1, 2, -1, 3, 2, 1, 3, 0, 1, -1&
455 : , -3, 1, 1, 1, 1, -1, 1, 1, 2, -3, 1, 2, 1, 0, 1, 1, 1, -1, -1, 1, 2, 1, 1, -1, 1, 1, -2, &
456 : 1, -3, 0, 0, 2, -1, 0, 0, 1, 2, -2, -1, -2, 2, 1, -2, -2, -3, 0, 0, 1, 3, 0, 0, 1, 1, 1, &
457 : -1, 1, 1, -3, 1, -1, 1, 0, 0, 0, 3, 0, 0, -1, -2, 0, -1, 0, -1, 3, 2, -1, 3, 0, 0, 0, 1, 0&
458 : , 0, -1, -1, 0, 1, 0, -2, -1, -1, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0,&
459 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,&
460 : 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 &
461 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
462 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
463 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
464 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20&
465 : ])
466 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: ma1 = RESHAPE([2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7 &
467 : , 5, 6, 7, 8, 0, 2, 3, 4, 2, 2, 4, 5, 3, 4, 4, 6, 4, 5, 6, 6, 0, 1, 1, 3, 1, 0, 3, 4, 1, 3&
468 : , 2, 5, 3, 4, 5, 4, 0, 0, 0, 2, 0, 0, 2, 3, 0, 2, 0, 4, 2, 3, 4, 2, 0, 0, 0, 1, 0, 0, 1, 2&
469 : , 0, 1, 0, 3, 1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 2, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0&
470 : , 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
471 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
472 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
473 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
474 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
475 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
476 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
477 : , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
478 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: ma2 = RESHAPE([2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7 &
479 : , 5, 6, 7, 8, 1, 2, 3, 4, 2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7, 1, 1, 3, 4, 2, 3, 3, 5, 3, 4&
480 : , 5, 5, 4, 5, 6, 7, 0, 0, 2, 3, 1, 2, 2, 4, 2, 3, 4, 4, 3, 4, 5, 6, 0, 0, 2, 2, 1, 2, 1, 4&
481 : , 2, 2, 4, 3, 3, 4, 5, 6, 0, 0, 1, 1, 0, 1, 0, 3, 1, 1, 3, 2, 2, 3, 4, 5, 0, 0, 1, 1, 0, 1&
482 : , 0, 3, 1, 1, 3, 1, 2, 3, 4, 5, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 2, 0, 1, 2, 3, 4, 0, 0, 0, 0&
483 : , 0, 0, 0, 2, 0, 0, 2, 0, 1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 2, 3, 0, 0&
484 : , 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2&
485 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
486 : , 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
487 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
488 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
489 : , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
490 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: ma3 = RESHAPE([2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7 &
491 : , 5, 6, 7, 8, 1, 2, 3, 4, 2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7, 1, 2, 3, 4, 1, 3, 4, 5, 3, 3&
492 : , 5, 6, 4, 5, 5, 7, 0, 1, 2, 3, 0, 2, 3, 4, 2, 2, 4, 5, 3, 4, 4, 6, 0, 1, 2, 3, 0, 2, 2, 4&
493 : , 2, 1, 4, 5, 2, 4, 3, 6, 0, 0, 1, 2, 0, 1, 1, 3, 1, 0, 3, 4, 1, 3, 2, 5, 0, 0, 1, 2, 0, 1&
494 : , 1, 3, 1, 0, 3, 4, 1, 3, 1, 5, 0, 0, 0, 1, 0, 0, 0, 2, 0, 0, 2, 3, 0, 2, 0, 4, 0, 0, 0, 1&
495 : , 0, 0, 0, 2, 0, 0, 2, 3, 0, 2, 0, 4, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 2, 0, 1, 0, 3, 0, 0&
496 : , 0, 0, 0, 0, 0, 1, 0, 0, 1, 2, 0, 1, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 2&
497 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
498 : , 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
499 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
500 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
501 : , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
502 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: ma4 = RESHAPE([2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7 &
503 : , 5, 6, 7, 8, 2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7, 5, 6, 7, 8, 0, 2, 3, 4, 2, 2, 4, 5, 3, 4&
504 : , 4, 6, 4, 5, 6, 6, 0, 2, 3, 4, 2, 2, 4, 5, 3, 4, 4, 6, 4, 5, 6, 6, 0, 1, 2, 3, 1, 0, 3, 4&
505 : , 2, 3, 4, 5, 3, 4, 5, 6, 0, 1, 2, 3, 1, 0, 3, 4, 2, 3, 2, 5, 3, 4, 5, 4, 0, 0, 2, 3, 0, 0&
506 : , 2, 3, 2, 2, 2, 5, 3, 3, 5, 4, 0, 0, 1, 2, 0, 0, 2, 3, 1, 2, 2, 4, 2, 3, 4, 2, 0, 0, 1, 2&
507 : , 0, 0, 1, 2, 1, 1, 0, 4, 2, 2, 4, 2, 0, 0, 0, 2, 0, 0, 1, 2, 0, 1, 0, 4, 2, 2, 4, 2, 0, 0&
508 : , 0, 1, 0, 0, 0, 1, 0, 0, 0, 3, 1, 1, 3, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 3, 1, 1, 3, 0&
509 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0&
510 : , 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2&
511 : , 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
512 : , 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
513 : , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
514 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: ma5 = RESHAPE([2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7 &
515 : , 5, 6, 7, 8, 0, 2, 3, 4, 2, 2, 4, 5, 3, 4, 4, 6, 4, 5, 6, 6, 0, 1, 2, 3, 1, 2, 3, 4, 2, 3&
516 : , 4, 5, 3, 4, 5, 6, 0, 0, 2, 3, 0, 0, 3, 3, 2, 3, 2, 5, 3, 3, 5, 4, 0, 0, 1, 2, 0, 0, 2, 3&
517 : , 1, 2, 2, 4, 2, 3, 4, 4, 0, 0, 0, 2, 0, 0, 2, 2, 0, 2, 0, 4, 2, 2, 4, 2, 0, 0, 0, 1, 0, 0&
518 : , 1, 1, 0, 1, 0, 3, 1, 1, 3, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 3, 0, 0, 0, 0, 0&
519 : , 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 2, 0, 0, 0&
520 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
521 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
522 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
523 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
524 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
525 : , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
526 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: mb1 = RESHAPE([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 &
527 : , 0, 0, 0, 0, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 0, 2, 3, 2, 2, 4, 2, 2, 3, 2&
528 : , 4, 2, 2, 2, 2, 4, 0, 3, 4, 3, 3, 0, 3, 3, 4, 3, 6, 3, 3, 3, 3, 6, 0, 0, 0, 4, 0, 0, 4, 4&
529 : , 0, 4, 0, 4, 4, 4, 4, 8, 0, 0, 0, 5, 0, 0, 5, 5, 0, 5, 0, 5, 5, 5, 5, 0, 0, 0, 0, 0, 0, 0&
530 : , 0, 6, 0, 0, 0, 6, 0, 6, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 0, 0, 7, 0, 0, 0, 0, 0&
531 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
532 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
533 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
534 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
535 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
536 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
537 : , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
538 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: mb2 = RESHAPE([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 &
539 : , 1, 1, 1, 1, 2, 0, 2, 2, 2, 2, 0, 2, 2, 2, 2, 0, 2, 2, 2, 2, 0, 3, 0, 0, 0, 0, 3, 0, 0, 0&
540 : , 0, 3, 0, 0, 0, 0, 1, 2, 3, 1, 3, 3, 2, 3, 3, 1, 3, 2, 3, 3, 3, 3, 0, 0, 1, 4, 1, 1, 5, 1&
541 : , 1, 4, 1, 5, 1, 1, 1, 1, 0, 0, 4, 5, 2, 4, 4, 4, 4, 5, 4, 4, 4, 4, 4, 4, 0, 0, 2, 3, 0, 2&
542 : , 0, 2, 2, 3, 2, 7, 2, 2, 2, 2, 0, 0, 3, 4, 0, 3, 0, 5, 3, 4, 5, 6, 5, 5, 5, 5, 0, 0, 0, 0&
543 : , 0, 0, 0, 3, 0, 0, 3, 0, 3, 3, 3, 3, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 6, 0, 4, 6, 6, 6, 0, 0&
544 : , 0, 0, 0, 0, 0, 4, 0, 0, 4, 0, 0, 4, 4, 4, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 5, 0, 0, 5, 7, 7&
545 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
546 : , 6, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
547 : , 0, 0, 0, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
548 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
549 : , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
550 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: mb3 = RESHAPE([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 &
551 : , 1, 1, 1, 1, 2, 2, 2, 2, 0, 2, 2, 2, 2, 0, 2, 2, 2, 2, 0, 2, 0, 0, 0, 0, 3, 0, 0, 0, 0, 3&
552 : , 0, 0, 0, 0, 3, 0, 1, 3, 3, 3, 2, 3, 1, 3, 3, 2, 3, 3, 1, 3, 2, 3, 0, 1, 1, 1, 0, 1, 4, 1&
553 : , 1, 5, 1, 1, 4, 1, 5, 1, 0, 2, 4, 4, 0, 4, 5, 4, 4, 4, 4, 4, 5, 4, 4, 4, 0, 0, 2, 2, 0, 2&
554 : , 3, 2, 2, 0, 2, 2, 3, 2, 7, 2, 0, 0, 3, 5, 0, 3, 4, 5, 3, 0, 5, 5, 4, 5, 6, 5, 0, 0, 0, 3&
555 : , 0, 0, 0, 3, 0, 0, 3, 3, 0, 3, 0, 3, 0, 0, 0, 4, 0, 0, 0, 6, 0, 0, 6, 6, 0, 6, 0, 6, 0, 0&
556 : , 0, 0, 0, 0, 0, 4, 0, 0, 4, 4, 0, 4, 0, 4, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 5, 7, 0, 5, 0, 7&
557 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 0, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0&
558 : , 0, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
559 : , 0, 0, 0, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
560 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
561 : , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
562 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: mb4 = RESHAPE([2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 &
563 : , 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 3, 3, 3, 4, 3, 3, 3, 3&
564 : , 4, 3, 3, 3, 3, 4, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 0, 2, 4, 4, 2, 4, 4, 2&
565 : , 4, 4, 0, 4, 4, 2, 4, 0, 0, 0, 2, 2, 0, 2, 0, 0, 2, 0, 6, 2, 2, 0, 2, 6, 0, 3, 0, 0, 3, 0&
566 : , 5, 5, 0, 5, 4, 0, 0, 5, 0, 2, 0, 1, 3, 5, 1, 0, 1, 1, 3, 1, 2, 5, 5, 1, 5, 8, 0, 0, 1, 3&
567 : , 0, 0, 4, 6, 1, 4, 6, 3, 3, 6, 3, 6, 0, 0, 4, 1, 0, 0, 2, 4, 4, 2, 4, 1, 1, 4, 1, 4, 0, 0&
568 : , 2, 4, 0, 0, 5, 5, 2, 5, 0, 6, 4, 5, 6, 8, 0, 0, 0, 2, 0, 0, 3, 3, 0, 3, 0, 4, 2, 3, 4, 6&
569 : , 0, 0, 0, 5, 0, 0, 0, 6, 0, 0, 0, 2, 5, 6, 2, 0, 0, 0, 0, 3, 0, 0, 0, 4, 0, 0, 0, 7, 3, 4&
570 : , 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3&
571 : , 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
572 : , 0, 4, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 0, 0, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0&
573 : , 0, 0, 0, 5, 0, 0, 5, 0], [4, 4, 20])
574 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: mb5 = RESHAPE([2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 &
575 : , 2, 2, 2, 2, 0, 3, 3, 3, 3, 4, 3, 3, 3, 3, 4, 3, 3, 3, 3, 4, 0, 0, 4, 4, 0, 0, 4, 0, 4, 4&
576 : , 0, 4, 4, 0, 4, 0, 0, 1, 0, 0, 1, 2, 0, 5, 0, 0, 6, 0, 0, 5, 0, 6, 0, 0, 1, 5, 0, 0, 5, 1&
577 : , 1, 5, 2, 5, 5, 1, 5, 2, 0, 0, 2, 1, 0, 0, 1, 6, 2, 1, 4, 1, 1, 6, 1, 8, 0, 0, 0, 2, 0, 0&
578 : , 2, 3, 0, 2, 0, 6, 2, 3, 6, 4, 0, 0, 0, 3, 0, 0, 3, 4, 0, 3, 0, 2, 3, 4, 2, 6, 0, 0, 0, 0&
579 : , 0, 0, 0, 0, 0, 0, 0, 7, 0, 0, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 3, 0, 0, 0&
580 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 5, 0&
581 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
582 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
583 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
584 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
585 : , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
586 :
587 : INTEGER :: k, k1, k2, mu
588 : REAL(kind=dp) :: cp, ct, fac1, fac2, J, Jc, Jcc, Jss, rr, &
589 : sp, st, xx, yy, za, zb
590 : REAL(kind=dp), DIMENSION(3) :: v
591 : REAL(kind=dp), DIMENSION(3, 3) :: Arot
592 :
593 0 : S(:, :) = 0.0_dp
594 :
595 0 : v(:) = R(:)
596 0 : rr = NORM2(v)
597 :
598 0 : IF (rr < 1.0e-20_dp) THEN
599 :
600 0 : DO mu = 1, 4
601 0 : S(mu, mu) = 1.0_dp
602 : END DO
603 :
604 : ELSE
605 :
606 0 : fac1 = 1.0_dp
607 0 : IF (nra == 1) THEN
608 : fac1 = fac1*2.0_dp
609 : ELSE
610 : IF (nra == 2) THEN
611 : fac1 = fac1*SQRT(4.0_dp/3.0_dp)
612 : ELSE
613 : IF (nra == 3) THEN
614 : fac1 = fac1*SQRT(8.0_dp/45.0_dp)
615 : ELSE
616 : IF (nra == 4) THEN
617 : fac1 = fac1*SQRT(4.0_dp/315.0_dp)
618 : ELSE
619 0 : WRITE (*, *) 'nra= ', nra
620 0 : RETURN
621 : END IF
622 : END IF
623 : END IF
624 : END IF
625 0 : IF (nrb == 1) THEN
626 0 : fac1 = fac1*2.0_dp
627 : ELSE
628 0 : IF (nrb == 2) THEN
629 0 : fac1 = fac1*SQRT(4.0_dp/3.0_dp)
630 : ELSE
631 0 : IF (nrb == 3) THEN
632 0 : fac1 = fac1*SQRT(8.0_dp/45.0_dp)
633 : ELSE
634 0 : IF (nrb == 4) THEN
635 0 : fac1 = fac1*SQRT(4.0_dp/315.0_dp)
636 : ELSE
637 0 : WRITE (*, *) 'nrb= ', nrb
638 0 : RETURN
639 : END IF
640 : END IF
641 : END IF
642 : END IF
643 :
644 0 : ct = -v(3)/rr
645 0 : IF (ABS(ct) < 1.0_dp) THEN
646 0 : st = SQRT(1.0_dp - ct**2)
647 0 : cp = -v(1)/(rr*st)
648 0 : sp = -v(2)/(rr*st)
649 0 : Arot(1, 1) = ct*cp
650 0 : Arot(1, 2) = -sp
651 0 : Arot(1, 3) = st*cp
652 0 : Arot(2, 1) = ct*sp
653 0 : Arot(2, 2) = cp
654 0 : Arot(2, 3) = st*sp
655 0 : Arot(3, 1) = -st
656 0 : Arot(3, 2) = 0.0_dp
657 0 : Arot(3, 3) = ct
658 : ELSE
659 0 : Arot(1, 1) = ct
660 0 : Arot(1, 2) = 0.0_dp
661 0 : Arot(1, 3) = 0.0_dp
662 0 : Arot(2, 1) = 0.0_dp
663 0 : Arot(2, 2) = 1.0_dp
664 0 : Arot(2, 3) = 0.0_dp
665 0 : Arot(3, 1) = 0.0_dp
666 0 : Arot(3, 2) = 0.0_dp
667 0 : Arot(3, 3) = ct
668 : END IF
669 :
670 0 : za = ZSA
671 0 : zb = ZSB
672 0 : fac2 = SQRT(za**(2*nra + 1)*zb**(2*nrb + 1))
673 0 : xx = 0.5_dp*rr*(za + zb)
674 0 : yy = 0.5_dp*rr*(za - zb)
675 :
676 0 : J = 0.0_dp
677 0 : DO k = 1, nc1(nra, nrb)
678 0 : J = J + REAL(c1(nra, nrb, k), dp)*AA(ma1(nra, nrb, k), xx)*BB(mb1(nra, nrb, k), yy)
679 : END DO
680 0 : J = J*rr**(nra + nrb + 1)
681 0 : J = J/2.0_dp**(nra + nrb + 2)
682 :
683 0 : S(1, 1) = S(1, 1) + fac1*fac2*J
684 :
685 0 : za = ZPA
686 0 : zb = ZSB
687 0 : fac2 = SQRT(za**(2*nra + 1)*zb**(2*nrb + 1))
688 0 : xx = 0.5_dp*rr*(za + zb)
689 0 : yy = 0.5_dp*rr*(za - zb)
690 :
691 0 : Jc = 0.0_dp
692 0 : DO k = 1, nc2(nra, nrb)
693 0 : Jc = Jc + REAL(c2(nra, nrb, k), dp)*AA(ma2(nra, nrb, k), xx)*BB(mb2(nra, nrb, k), yy)
694 : END DO
695 0 : Jc = Jc*rr**(nra + nrb + 1)
696 0 : Jc = Jc/2.0_dp**(nra + nrb + 2)
697 :
698 0 : DO k1 = 1, 3
699 : S(k1 + 1, 1) = S(k1 + 1, 1) &
700 0 : & + SQRT(3.0_dp)*Arot(k1, 3)*fac1*fac2*Jc
701 : END DO
702 :
703 0 : za = ZSA
704 0 : zb = ZPB
705 0 : fac2 = SQRT(za**(2*nra + 1)*zb**(2*nrb + 1))
706 0 : xx = 0.5_dp*rr*(za + zb)
707 0 : yy = 0.5_dp*rr*(za - zb)
708 :
709 0 : Jc = 0.0_dp
710 0 : DO k = 1, nc3(nra, nrb)
711 0 : Jc = Jc + REAL(c3(nra, nrb, k), dp)*AA(ma3(nra, nrb, k), xx)*BB(mb3(nra, nrb, k), yy)
712 : END DO
713 0 : Jc = Jc*rr**(nra + nrb + 1)
714 0 : Jc = Jc/2.0_dp**(nra + nrb + 2)
715 :
716 0 : DO k1 = 1, 3
717 : S(1, k1 + 1) = S(1, k1 + 1) &
718 0 : & - SQRT(3.0_dp)*Arot(k1, 3)*fac1*fac2*Jc
719 : END DO
720 :
721 0 : za = ZPA
722 0 : zb = ZPB
723 0 : fac2 = SQRT(za**(2*nra + 1)*zb**(2*nrb + 1))
724 0 : xx = 0.5_dp*rr*(za + zb)
725 0 : yy = 0.5_dp*rr*(za - zb)
726 :
727 0 : Jss = 0.0_dp
728 0 : DO k = 1, nc4(nra, nrb)
729 0 : Jss = Jss + REAL(c4(nra, nrb, k), dp)*AA(ma4(nra, nrb, k), xx)*BB(mb4(nra, nrb, k), yy)
730 : END DO
731 0 : Jss = Jss*rr**(nra + nrb + 1)
732 0 : Jss = Jss/2.0_dp**(nra + nrb + 2)
733 :
734 0 : Jcc = 0.0_dp
735 0 : DO k = 1, nc5(nra, nrb)
736 0 : Jcc = Jcc + REAL(c5(nra, nrb, k), dp)*AA(ma5(nra, nrb, k), xx)*BB(mb5(nra, nrb, k), yy)
737 : END DO
738 0 : Jcc = Jcc*rr**(nra + nrb + 1)
739 0 : Jcc = Jcc/2.0_dp**(nra + nrb + 2)
740 :
741 0 : DO k1 = 1, 3
742 0 : DO k2 = 1, 3
743 : S(k1 + 1, k2 + 1) = S(k1 + 1, k2 + 1) &
744 : & + 1.5_dp*Arot(k1, 1)*Arot(k2, 1)*fac1*fac2*Jss &
745 : & + 1.5_dp*Arot(k1, 2)*Arot(k2, 2)*fac1*fac2*Jss &
746 0 : & - 3.0_dp*Arot(k1, 3)*Arot(k2, 3)*fac1*fac2*Jcc
747 : END DO
748 : END DO
749 :
750 : END IF
751 :
752 : END SUBROUTINE makeS
753 :
754 : ! **************************************************************************************************
755 : !> \brief ...
756 : !> \param R ...
757 : !> \param nra ...
758 : !> \param nrb ...
759 : !> \param ZSA ...
760 : !> \param ZSB ...
761 : !> \param ZPA ...
762 : !> \param ZPB ...
763 : !> \param dS ...
764 : ! **************************************************************************************************
765 0 : SUBROUTINE makedS(R, nra, nrb, ZSA, ZSB, ZPA, ZPB, dS)
766 :
767 : REAL(kind=dp), DIMENSION(3) :: R
768 : INTEGER :: nra, nrb
769 : REAL(kind=dp) :: ZSA, ZSB, ZPA, ZPB
770 : REAL(kind=dp), DIMENSION(4, 4, 3) :: dS
771 :
772 : INTEGER, DIMENSION(4, 4), PARAMETER :: &
773 : nc1 = RESHAPE([2, 4, 4, 6, 4, 3, 6, 7, 4, 6, 4, 8, 6, 7, 8, 5], [4, 4]), &
774 : nc2 = RESHAPE([4, 4, 8, 8, 6, 8, 6, 12, 8, 8, 12, 8, 10, 12, 14, 16], [4, 4]), &
775 : nc3 = RESHAPE([4, 6, 8, 10, 4, 8, 8, 12, 8, 6, 12, 14, 8, 12, 8, 16], [4, 4]), &
776 : nc4 = RESHAPE([4, 8, 11, 14, 8, 6, 12, 14, 11, 12, 10, 20, 14, 14, 20, 12], [4, 4]), &
777 : nc5 = RESHAPE([2, 4, 6, 8, 4, 4, 8, 8, 6, 8, 6, 12, 8, 8, 12, 8], [4, 4])
778 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: c1 = RESHAPE([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, &
779 : 1, 1, 1, 1, -1, 1, 2, 3, -1, -2, 1, 2, -2, -1, -3, 1, -3, -2, -1, -4, 0, -1, -2, 2, -1, 1,&
780 : -2, -1, 2, -2, 3, -3, 2, -1, -3, 6, 0, -1, -1, -2, 1, 0, -2, -4, -1, 2, -1, -3, 2, 4, 3, &
781 : -4, 0, 0, 0, -3, 0, 0, 1, -1, 0, 1, 0, 3, -3, -1, 3, 1, 0, 0, 0, -1, 0, 0, 1, 2, 0, -1, 0,&
782 : 3, 1, -2, -3, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, -1, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0 &
783 : , 0, 0, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
784 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
785 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
786 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
787 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
788 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
789 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
790 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: c2 = RESHAPE([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, &
791 : 1, 1, 1, 1, -1, 1, 1, 2, -2, -1, 1, 1, -3, -2, -1, 1, -4, -3, -2, -1, 1, -1, 1, 1, 1, 1, &
792 : -2, 1, 1, 1, 1, -3, 1, 1, 1, 1, -1, -1, -1, 2, 1, -1, -2, -2, 3, -2, -2, -3, 6, 2, -1, -3,&
793 : 0, 0, 1, -2, -2, -1, 1, 1, -3, 2, -1, 3, -4, -3, -2, -1, 0, 0, -1, -1, 1, 1, 1, -2, -1, -1&
794 : , 2, 3, -4, 2, 4, 3, 0, 0, -1, -2, 0, -1, 0, -2, 3, 2, -2, -1, 6, 2, -1, -3, 0, 0, -1, -1,&
795 : 0, 1, 0, 1, -1, -1, 1, -1, 1, -3, -1, 3, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 2, 0, -4, 2, 4, 3 &
796 : , 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 1, 1, -2, -3, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0&
797 : , -3, -1, 3, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 0, 1, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
798 : 0, 0, 0, 0, 0, -2, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0&
799 : , 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, &
800 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
801 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],&
802 : [4, 4, 20])
803 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: c3 = RESHAPE([-1, -1, -1, -1, -1, -1, -1, -1, -1, &
804 : -1, -1, -1, -1, -1, -1, -1, -1, -2, -3, -4, 1, -1, -2, -3, 1, 1, -1, -2, 2, 1, 1, -1, 1, 1&
805 : , 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 3, 1, 1, -1, -3, -6, -1, 1, 2, -2, 1, -2, 2, 1, -2, &
806 : 2, -3, 3, 0, 2, 3, 4, 0, 1, 2, 3, -1, -1, 1, 2, -2, -1, -3, 1, 0, 1, -1, -4, 0, 1, 1, 2, &
807 : -1, 1, 2, 4, 1, -2, 3, 3, 0, 0, 3, 6, 0, -1, -2, 2, -1, 0, -2, -1, 2, -2, 1, -3, 0, 0, 1, &
808 : -1, 0, -1, -1, 3, 1, 0, -1, 1, -1, -1, -1, -3, 0, 0, 0, 4, 0, 0, 0, -2, 0, 0, -2, -4, 0, 2&
809 : , 0, -3, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, -1, -2, 0, 1, 0, -3, 0, 0, 0, 0, 0, 0, 0, -3, 0, 0,&
810 : 1, -1, 0, 1, 0, 3, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, -1, 0, -1, 0, 1, 0, 0, 0, 0, 0, 0, 0,&
811 : 0, 0, 0, 0, 2, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, &
812 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0&
813 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
814 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
815 : , 0], [4, 4, 20])
816 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: c4 = RESHAPE([-1, -1, -1, -1, -1, -1, -1, -1, -1, &
817 : -1, -1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -2, -3,&
818 : 1, 1, -1, -2, 2, 1, 2, -1, 3, 2, 1, 3, -1, 1, 2, 3, -1, -1, 1, 2, -2, -1, -1, 1, -3, -2, &
819 : -1, -2, 0, 1, -1, -3, 1, -1, 1, 1, -1, 1, -1, 2, -3, 1, 2, -1, 0, -1, 2, 4, -1, 1, -1, -1 &
820 : , 2, -1, -1, -1, 4, -1, -1, -3, 0, 1, -1, -1, -1, 0, 1, 2, -1, -1, -1, -1, -1, -2, -1, 3, &
821 : 0, -1, 2, -1, 1, 0, -1, -2, -2, 1, 2, 2, 1, 2, -2, 1, 0, 0, -2, 4, 0, 0, -1, 1, 2, -1, 1, &
822 : -1, -4, 1, 1, 2, 0, 0, 1, -3, 0, 0, 1, -1, 1, 1, -1, -1, 3, -1, 1, -3, 0, 0, -1, 3, 0, 0, &
823 : -1, -2, -1, 1, 0, -1, 3, 2, -1, -1, 0, 0, 0, -3, 0, 0, 1, 2, 0, -1, 0, -1, -3, -2, -1, 1, &
824 : 0, 0, 0, 1, 0, 0, 0, -1, 0, 0, 0, 2, -1, -1, 2, 0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, 0, -1, 1&
825 : , 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
826 : 0, 2, 0, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,&
827 : 0, 0, 0, -1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 0, &
828 : 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0], [4, 4, 20])
829 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: c5 = RESHAPE([-1, -1, -1, -1, -1, -1, -1, -1, -1, &
830 : -1, -1, -1, -1, -1, -1, -1, 1, -1, -2, -3, 1, 1, -1, -2, 2, 1, 2, -1, 3, 2, 1, 3, 0, 1, -1&
831 : , -3, 1, 1, 1, 1, -1, 1, 1, 2, -3, 1, 2, 1, 0, 1, 1, 1, -1, -1, 1, 2, 1, 1, -1, 1, 1, -2, &
832 : 1, -3, 0, 0, 2, -1, 0, 0, 1, 2, -2, -1, -2, 2, 1, -2, -2, -3, 0, 0, 1, 3, 0, 0, 1, 1, 1, &
833 : -1, 1, 1, -3, 1, -1, 1, 0, 0, 0, 3, 0, 0, -1, -2, 0, -1, 0, -1, 3, 2, -1, 3, 0, 0, 0, 1, 0&
834 : , 0, -1, -1, 0, 1, 0, -2, -1, -1, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0,&
835 : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,&
836 : 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 &
837 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
838 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
839 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
840 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20&
841 : ])
842 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: ma1 = RESHAPE([2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7 &
843 : , 5, 6, 7, 8, 0, 2, 3, 4, 2, 2, 4, 5, 3, 4, 4, 6, 4, 5, 6, 6, 0, 1, 1, 3, 1, 0, 3, 4, 1, 3&
844 : , 2, 5, 3, 4, 5, 4, 0, 0, 0, 2, 0, 0, 2, 3, 0, 2, 0, 4, 2, 3, 4, 2, 0, 0, 0, 1, 0, 0, 1, 2&
845 : , 0, 1, 0, 3, 1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 2, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0&
846 : , 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
847 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
848 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
849 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
850 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
851 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
852 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
853 : , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
854 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: ma2 = RESHAPE([2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7 &
855 : , 5, 6, 7, 8, 1, 2, 3, 4, 2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7, 1, 1, 3, 4, 2, 3, 3, 5, 3, 4&
856 : , 5, 5, 4, 5, 6, 7, 0, 0, 2, 3, 1, 2, 2, 4, 2, 3, 4, 4, 3, 4, 5, 6, 0, 0, 2, 2, 1, 2, 1, 4&
857 : , 2, 2, 4, 3, 3, 4, 5, 6, 0, 0, 1, 1, 0, 1, 0, 3, 1, 1, 3, 2, 2, 3, 4, 5, 0, 0, 1, 1, 0, 1&
858 : , 0, 3, 1, 1, 3, 1, 2, 3, 4, 5, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 2, 0, 1, 2, 3, 4, 0, 0, 0, 0&
859 : , 0, 0, 0, 2, 0, 0, 2, 0, 1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 2, 3, 0, 0&
860 : , 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2&
861 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
862 : , 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
863 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
864 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
865 : , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
866 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: ma3 = RESHAPE([2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7 &
867 : , 5, 6, 7, 8, 1, 2, 3, 4, 2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7, 1, 2, 3, 4, 1, 3, 4, 5, 3, 3&
868 : , 5, 6, 4, 5, 5, 7, 0, 1, 2, 3, 0, 2, 3, 4, 2, 2, 4, 5, 3, 4, 4, 6, 0, 1, 2, 3, 0, 2, 2, 4&
869 : , 2, 1, 4, 5, 2, 4, 3, 6, 0, 0, 1, 2, 0, 1, 1, 3, 1, 0, 3, 4, 1, 3, 2, 5, 0, 0, 1, 2, 0, 1&
870 : , 1, 3, 1, 0, 3, 4, 1, 3, 1, 5, 0, 0, 0, 1, 0, 0, 0, 2, 0, 0, 2, 3, 0, 2, 0, 4, 0, 0, 0, 1&
871 : , 0, 0, 0, 2, 0, 0, 2, 3, 0, 2, 0, 4, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 2, 0, 1, 0, 3, 0, 0&
872 : , 0, 0, 0, 0, 0, 1, 0, 0, 1, 2, 0, 1, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 2&
873 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
874 : , 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
875 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
876 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
877 : , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
878 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: ma4 = RESHAPE([2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7 &
879 : , 5, 6, 7, 8, 2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7, 5, 6, 7, 8, 0, 2, 3, 4, 2, 2, 4, 5, 3, 4&
880 : , 4, 6, 4, 5, 6, 6, 0, 2, 3, 4, 2, 2, 4, 5, 3, 4, 4, 6, 4, 5, 6, 6, 0, 1, 2, 3, 1, 0, 3, 4&
881 : , 2, 3, 4, 5, 3, 4, 5, 6, 0, 1, 2, 3, 1, 0, 3, 4, 2, 3, 2, 5, 3, 4, 5, 4, 0, 0, 2, 3, 0, 0&
882 : , 2, 3, 2, 2, 2, 5, 3, 3, 5, 4, 0, 0, 1, 2, 0, 0, 2, 3, 1, 2, 2, 4, 2, 3, 4, 2, 0, 0, 1, 2&
883 : , 0, 0, 1, 2, 1, 1, 0, 4, 2, 2, 4, 2, 0, 0, 0, 2, 0, 0, 1, 2, 0, 1, 0, 4, 2, 2, 4, 2, 0, 0&
884 : , 0, 1, 0, 0, 0, 1, 0, 0, 0, 3, 1, 1, 3, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 3, 1, 1, 3, 0&
885 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0&
886 : , 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2&
887 : , 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
888 : , 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
889 : , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
890 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: ma5 = RESHAPE([2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7 &
891 : , 5, 6, 7, 8, 0, 2, 3, 4, 2, 2, 4, 5, 3, 4, 4, 6, 4, 5, 6, 6, 0, 1, 2, 3, 1, 2, 3, 4, 2, 3&
892 : , 4, 5, 3, 4, 5, 6, 0, 0, 2, 3, 0, 0, 3, 3, 2, 3, 2, 5, 3, 3, 5, 4, 0, 0, 1, 2, 0, 0, 2, 3&
893 : , 1, 2, 2, 4, 2, 3, 4, 4, 0, 0, 0, 2, 0, 0, 2, 2, 0, 2, 0, 4, 2, 2, 4, 2, 0, 0, 0, 1, 0, 0&
894 : , 1, 1, 0, 1, 0, 3, 1, 1, 3, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 3, 0, 0, 0, 0, 0&
895 : , 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 2, 0, 0, 0&
896 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
897 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
898 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
899 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
900 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
901 : , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
902 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: mb1 = RESHAPE([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 &
903 : , 0, 0, 0, 0, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 0, 2, 3, 2, 2, 4, 2, 2, 3, 2&
904 : , 4, 2, 2, 2, 2, 4, 0, 3, 4, 3, 3, 0, 3, 3, 4, 3, 6, 3, 3, 3, 3, 6, 0, 0, 0, 4, 0, 0, 4, 4&
905 : , 0, 4, 0, 4, 4, 4, 4, 8, 0, 0, 0, 5, 0, 0, 5, 5, 0, 5, 0, 5, 5, 5, 5, 0, 0, 0, 0, 0, 0, 0&
906 : , 0, 6, 0, 0, 0, 6, 0, 6, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 0, 0, 7, 0, 0, 0, 0, 0&
907 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
908 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
909 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
910 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
911 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
912 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
913 : , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
914 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: mb2 = RESHAPE([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 &
915 : , 1, 1, 1, 1, 2, 0, 2, 2, 2, 2, 0, 2, 2, 2, 2, 0, 2, 2, 2, 2, 0, 3, 0, 0, 0, 0, 3, 0, 0, 0&
916 : , 0, 3, 0, 0, 0, 0, 1, 2, 3, 1, 3, 3, 2, 3, 3, 1, 3, 2, 3, 3, 3, 3, 0, 0, 1, 4, 1, 1, 5, 1&
917 : , 1, 4, 1, 5, 1, 1, 1, 1, 0, 0, 4, 5, 2, 4, 4, 4, 4, 5, 4, 4, 4, 4, 4, 4, 0, 0, 2, 3, 0, 2&
918 : , 0, 2, 2, 3, 2, 7, 2, 2, 2, 2, 0, 0, 3, 4, 0, 3, 0, 5, 3, 4, 5, 6, 5, 5, 5, 5, 0, 0, 0, 0&
919 : , 0, 0, 0, 3, 0, 0, 3, 0, 3, 3, 3, 3, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 6, 0, 4, 6, 6, 6, 0, 0&
920 : , 0, 0, 0, 0, 0, 4, 0, 0, 4, 0, 0, 4, 4, 4, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 5, 0, 0, 5, 7, 7&
921 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
922 : , 6, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
923 : , 0, 0, 0, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
924 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
925 : , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
926 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: mb3 = RESHAPE([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 &
927 : , 1, 1, 1, 1, 2, 2, 2, 2, 0, 2, 2, 2, 2, 0, 2, 2, 2, 2, 0, 2, 0, 0, 0, 0, 3, 0, 0, 0, 0, 3&
928 : , 0, 0, 0, 0, 3, 0, 1, 3, 3, 3, 2, 3, 1, 3, 3, 2, 3, 3, 1, 3, 2, 3, 0, 1, 1, 1, 0, 1, 4, 1&
929 : , 1, 5, 1, 1, 4, 1, 5, 1, 0, 2, 4, 4, 0, 4, 5, 4, 4, 4, 4, 4, 5, 4, 4, 4, 0, 0, 2, 2, 0, 2&
930 : , 3, 2, 2, 0, 2, 2, 3, 2, 7, 2, 0, 0, 3, 5, 0, 3, 4, 5, 3, 0, 5, 5, 4, 5, 6, 5, 0, 0, 0, 3&
931 : , 0, 0, 0, 3, 0, 0, 3, 3, 0, 3, 0, 3, 0, 0, 0, 4, 0, 0, 0, 6, 0, 0, 6, 6, 0, 6, 0, 6, 0, 0&
932 : , 0, 0, 0, 0, 0, 4, 0, 0, 4, 4, 0, 4, 0, 4, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 5, 7, 0, 5, 0, 7&
933 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 0, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0&
934 : , 0, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
935 : , 0, 0, 0, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
936 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
937 : , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
938 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: mb4 = RESHAPE([2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 &
939 : , 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 3, 3, 3, 4, 3, 3, 3, 3&
940 : , 4, 3, 3, 3, 3, 4, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 0, 2, 4, 4, 2, 4, 4, 2&
941 : , 4, 4, 0, 4, 4, 2, 4, 0, 0, 0, 2, 2, 0, 2, 0, 0, 2, 0, 6, 2, 2, 0, 2, 6, 0, 3, 0, 0, 3, 0&
942 : , 5, 5, 0, 5, 4, 0, 0, 5, 0, 2, 0, 1, 3, 5, 1, 0, 1, 1, 3, 1, 2, 5, 5, 1, 5, 8, 0, 0, 1, 3&
943 : , 0, 0, 4, 6, 1, 4, 6, 3, 3, 6, 3, 6, 0, 0, 4, 1, 0, 0, 2, 4, 4, 2, 4, 1, 1, 4, 1, 4, 0, 0&
944 : , 2, 4, 0, 0, 5, 5, 2, 5, 0, 6, 4, 5, 6, 8, 0, 0, 0, 2, 0, 0, 3, 3, 0, 3, 0, 4, 2, 3, 4, 6&
945 : , 0, 0, 0, 5, 0, 0, 0, 6, 0, 0, 0, 2, 5, 6, 2, 0, 0, 0, 0, 3, 0, 0, 0, 4, 0, 0, 0, 7, 3, 4&
946 : , 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3&
947 : , 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
948 : , 0, 4, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 0, 0, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0&
949 : , 0, 0, 0, 5, 0, 0, 5, 0], [4, 4, 20])
950 : INTEGER, DIMENSION(4, 4, 20), PARAMETER :: mb5 = RESHAPE([2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 &
951 : , 2, 2, 2, 2, 0, 3, 3, 3, 3, 4, 3, 3, 3, 3, 4, 3, 3, 3, 3, 4, 0, 0, 4, 4, 0, 0, 4, 0, 4, 4&
952 : , 0, 4, 4, 0, 4, 0, 0, 1, 0, 0, 1, 2, 0, 5, 0, 0, 6, 0, 0, 5, 0, 6, 0, 0, 1, 5, 0, 0, 5, 1&
953 : , 1, 5, 2, 5, 5, 1, 5, 2, 0, 0, 2, 1, 0, 0, 1, 6, 2, 1, 4, 1, 1, 6, 1, 8, 0, 0, 0, 2, 0, 0&
954 : , 2, 3, 0, 2, 0, 6, 2, 3, 6, 4, 0, 0, 0, 3, 0, 0, 3, 4, 0, 3, 0, 2, 3, 4, 2, 6, 0, 0, 0, 0&
955 : , 0, 0, 0, 0, 0, 0, 0, 7, 0, 0, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 3, 0, 0, 0&
956 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 5, 0&
957 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
958 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
959 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
960 : , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
961 : , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
962 :
963 : INTEGER :: k, k1, k2, mu
964 : REAL(kind=dp) :: cp, ct, dJ, dJc, dJcc, dJss, dxx, dyy, &
965 : f, fac1, fac2, J, Jc, Jcc, Jss, rr, &
966 : sp, st, w, w1, w2, xx, yy, za, zb
967 : REAL(kind=dp), DIMENSION(3) :: dcp, dct, dsp, dst, v
968 : REAL(kind=dp), DIMENSION(3, 3) :: Arot
969 : REAL(kind=dp), DIMENSION(3, 3, 3) :: dArot
970 :
971 0 : dS(:, :, :) = 0.0_dp
972 :
973 0 : v(:) = R(:)
974 0 : rr = NORM2(v)
975 :
976 0 : IF (rr < 1.0e-20_dp) THEN
977 :
978 0 : DO mu = 1, 4
979 0 : dS(mu, mu, :) = 0.0_dp
980 : END DO
981 :
982 : ELSE
983 :
984 0 : fac1 = 1.0_dp
985 0 : IF (nra == 1) THEN
986 : fac1 = fac1*2.0_dp
987 : ELSE
988 : IF (nra == 2) THEN
989 : fac1 = fac1*SQRT(4.0_dp/3.0_dp)
990 : ELSE
991 : IF (nra == 3) THEN
992 : fac1 = fac1*SQRT(8.0_dp/45.0_dp)
993 : ELSE
994 : IF (nra == 4) THEN
995 : fac1 = fac1*SQRT(4.0_dp/315.0_dp)
996 : ELSE
997 0 : WRITE (*, *) 'nra= ', nra
998 0 : RETURN
999 : END IF
1000 : END IF
1001 : END IF
1002 : END IF
1003 0 : IF (nrb == 1) THEN
1004 0 : fac1 = fac1*2.0_dp
1005 : ELSE
1006 0 : IF (nrb == 2) THEN
1007 0 : fac1 = fac1*SQRT(4.0_dp/3.0_dp)
1008 : ELSE
1009 0 : IF (nrb == 3) THEN
1010 0 : fac1 = fac1*SQRT(8.0_dp/45.0_dp)
1011 : ELSE
1012 0 : IF (nrb == 4) THEN
1013 0 : fac1 = fac1*SQRT(4.0_dp/315.0_dp)
1014 : ELSE
1015 0 : WRITE (*, *) 'nrb= ', nrb
1016 0 : RETURN
1017 : END IF
1018 : END IF
1019 : END IF
1020 : END IF
1021 :
1022 0 : ct = -v(3)/rr
1023 0 : IF (ABS(ct) >= 1.0_dp) THEN
1024 :
1025 0 : dct(:) = v(:)*v(3)/rr**3
1026 0 : dct(3) = dct(3) - 1.0_dp/rr
1027 :
1028 0 : Arot(1, 1) = ct
1029 0 : Arot(1, 2) = 0.0_dp
1030 0 : Arot(1, 3) = 0.0_dp
1031 0 : Arot(2, 1) = 0.0_dp
1032 0 : Arot(2, 2) = 1.0_dp
1033 0 : Arot(2, 3) = 0.0_dp
1034 0 : Arot(3, 1) = 0.0_dp
1035 0 : Arot(3, 2) = 0.0_dp
1036 0 : Arot(3, 3) = ct
1037 :
1038 0 : dArot(1, 1, :) = dct(:)
1039 0 : dArot(1, 2, :) = 0.0_dp
1040 0 : dArot(1, 3, :) = 0.0_dp
1041 0 : dArot(2, 1, :) = 0.0_dp
1042 0 : dArot(2, 2, :) = 0.0_dp
1043 0 : dArot(2, 3, :) = 0.0_dp
1044 0 : dArot(3, 1, :) = 0.0_dp
1045 0 : dArot(3, 2, :) = 0.0_dp
1046 0 : dArot(3, 3, :) = dct(:)
1047 :
1048 : ELSE
1049 :
1050 0 : xx = SQRT(v(1)**2 + v(2)**2)
1051 0 : st = xx/rr
1052 0 : cp = -v(1)/xx
1053 0 : sp = -v(2)/xx
1054 :
1055 0 : dct(:) = v(:)*v(3)/rr**3
1056 0 : dct(3) = dct(3) - 1.0_dp/rr
1057 0 : dst(:) = -ct*dct(:)/st
1058 0 : dcp(:) = v(:)*v(1)/(rr**3*st)
1059 0 : dcp(:) = dcp(:) + v(1)*dst(:)/(rr*st**2)
1060 0 : dcp(1) = dcp(1) - 1.0_dp/(rr*st)
1061 0 : dsp(:) = v(:)*v(2)/(rr**3*st)
1062 0 : dsp(:) = dsp(:) + v(2)*dst(:)/(rr*st**2)
1063 0 : dsp(2) = dsp(2) - 1.0_dp/(rr*st)
1064 :
1065 0 : Arot(1, 1) = ct*cp
1066 0 : Arot(1, 2) = -sp
1067 0 : Arot(1, 3) = st*cp
1068 0 : Arot(2, 1) = ct*sp
1069 0 : Arot(2, 2) = cp
1070 0 : Arot(2, 3) = st*sp
1071 0 : Arot(3, 1) = -st
1072 0 : Arot(3, 2) = 0.0_dp
1073 0 : Arot(3, 3) = ct
1074 :
1075 0 : dArot(1, 1, :) = dct(:)*cp + ct*dcp(:)
1076 0 : dArot(1, 2, :) = -dsp(:)
1077 0 : dArot(1, 3, :) = dst(:)*cp + st*dcp(:)
1078 0 : dArot(2, 1, :) = dct(:)*sp + ct*dsp(:)
1079 0 : dArot(2, 2, :) = dcp(:)
1080 0 : dArot(2, 3, :) = dst(:)*sp + st*dsp(:)
1081 0 : dArot(3, 1, :) = -dst(:)
1082 0 : dArot(3, 2, :) = 0.0_dp
1083 0 : dArot(3, 3, :) = dct(:)
1084 :
1085 : END IF
1086 :
1087 0 : za = ZSA
1088 0 : zb = ZSB
1089 0 : fac2 = SQRT(za**(2*nra + 1)*zb**(2*nrb + 1))
1090 0 : xx = 0.5_dp*rr*(za + zb)
1091 0 : yy = 0.5_dp*rr*(za - zb)
1092 0 : dxx = 0.5_dp*(za + zb)
1093 0 : dyy = 0.5_dp*(za - zb)
1094 :
1095 0 : w = 0.0_dp
1096 0 : w1 = 0.0_dp
1097 0 : w2 = 0.0_dp
1098 0 : f = rr**(nra + nrb + 1)/2.0_dp**(nra + nrb + 2)
1099 0 : DO k = 1, nc1(nra, nrb)
1100 0 : w = w + REAL(c1(nra, nrb, k), dp)*AA(ma1(nra, nrb, k), xx)*BB(mb1(nra, nrb, k), yy)
1101 0 : w1 = w1 + REAL(c1(nra, nrb, k), dp)*AA(ma1(nra, nrb, k) + 1, xx)*BB(mb1(nra, nrb, k), yy)
1102 0 : w2 = w2 + REAL(c1(nra, nrb, k), dp)*AA(ma1(nra, nrb, k), xx)*BB(mb1(nra, nrb, k) + 1, yy)
1103 : END DO
1104 0 : J = f*w
1105 0 : dJ = f*REAL(nra + nrb + 1, dp)*w/rr
1106 0 : dJ = dJ - dxx*f*w1
1107 0 : dJ = dJ - dyy*f*w2
1108 :
1109 0 : dS(1, 1, :) = dS(1, 1, :) + fac1*fac2*dJ*v(:)/rr
1110 :
1111 0 : za = ZPA
1112 0 : zb = ZSB
1113 0 : fac2 = SQRT(za**(2*nra + 1)*zb**(2*nrb + 1))
1114 0 : xx = 0.5_dp*rr*(za + zb)
1115 0 : yy = 0.5_dp*rr*(za - zb)
1116 0 : dxx = 0.5_dp*(za + zb)
1117 0 : dyy = 0.5_dp*(za - zb)
1118 :
1119 0 : w = 0.0_dp
1120 0 : w1 = 0.0_dp
1121 0 : w2 = 0.0_dp
1122 0 : f = rr**(nra + nrb + 1)/2.0_dp**(nra + nrb + 2)
1123 0 : DO k = 1, nc2(nra, nrb)
1124 0 : w = w + REAL(c2(nra, nrb, k), dp)*AA(ma2(nra, nrb, k), xx)*BB(mb2(nra, nrb, k), yy)
1125 0 : w1 = w1 + REAL(c2(nra, nrb, k), dp)*AA(ma2(nra, nrb, k) + 1, xx)*BB(mb2(nra, nrb, k), yy)
1126 0 : w2 = w2 + REAL(c2(nra, nrb, k), dp)*AA(ma2(nra, nrb, k), xx)*BB(mb2(nra, nrb, k) + 1, yy)
1127 : END DO
1128 0 : Jc = f*w
1129 0 : dJc = f*REAL(nra + nrb + 1, dp)*w/rr
1130 0 : dJc = dJc - dxx*f*w1
1131 0 : dJc = dJc - dyy*f*w2
1132 :
1133 0 : DO k1 = 1, 3
1134 : dS(k1 + 1, 1, :) = dS(k1 + 1, 1, :) &
1135 : & + SQRT(3.0_dp)*Arot(k1, 3)*fac1*fac2*dJc*v(:)/rr &
1136 0 : & + SQRT(3.0_dp)*dArot(k1, 3, :)*fac1*fac2*Jc
1137 : END DO
1138 :
1139 0 : za = ZSA
1140 0 : zb = ZPB
1141 0 : fac2 = SQRT(za**(2*nra + 1)*zb**(2*nrb + 1))
1142 0 : xx = 0.5_dp*rr*(za + zb)
1143 0 : yy = 0.5_dp*rr*(za - zb)
1144 0 : dxx = 0.5_dp*(za + zb)
1145 0 : dyy = 0.5_dp*(za - zb)
1146 :
1147 0 : w = 0.0_dp
1148 0 : w1 = 0.0_dp
1149 0 : w2 = 0.0_dp
1150 0 : f = rr**(nra + nrb + 1)/2.0_dp**(nra + nrb + 2)
1151 0 : DO k = 1, nc3(nra, nrb)
1152 0 : w = w + REAL(c3(nra, nrb, k), dp)*AA(ma3(nra, nrb, k), xx)*BB(mb3(nra, nrb, k), yy)
1153 0 : w1 = w1 + REAL(c3(nra, nrb, k), dp)*AA(ma3(nra, nrb, k) + 1, xx)*BB(mb3(nra, nrb, k), yy)
1154 0 : w2 = w2 + REAL(c3(nra, nrb, k), dp)*AA(ma3(nra, nrb, k), xx)*BB(mb3(nra, nrb, k) + 1, yy)
1155 : END DO
1156 0 : Jc = f*w
1157 0 : dJc = f*REAL(nra + nrb + 1, dp)*w/rr
1158 0 : dJc = dJc - dxx*f*w1
1159 0 : dJc = dJc - dyy*f*w2
1160 :
1161 0 : DO k1 = 1, 3
1162 : dS(1, k1 + 1, :) = dS(1, k1 + 1, :) &
1163 : & - SQRT(3.0_dp)*Arot(k1, 3)*fac1*fac2*dJc*v(:)/rr &
1164 0 : & - SQRT(3.0_dp)*dArot(k1, 3, :)*fac1*fac2*Jc
1165 : END DO
1166 :
1167 0 : za = ZPA
1168 0 : zb = ZPB
1169 0 : fac2 = SQRT(za**(2*nra + 1)*zb**(2*nrb + 1))
1170 0 : xx = 0.5_dp*rr*(za + zb)
1171 0 : yy = 0.5_dp*rr*(za - zb)
1172 0 : dxx = 0.5_dp*(za + zb)
1173 0 : dyy = 0.5_dp*(za - zb)
1174 :
1175 0 : w = 0.0_dp
1176 0 : w1 = 0.0_dp
1177 0 : w2 = 0.0_dp
1178 0 : f = rr**(nra + nrb + 1)/2.0_dp**(nra + nrb + 2)
1179 0 : DO k = 1, nc4(nra, nrb)
1180 0 : w = w + REAL(c4(nra, nrb, k), dp)*AA(ma4(nra, nrb, k), xx)*BB(mb4(nra, nrb, k), yy)
1181 0 : w1 = w1 + REAL(c4(nra, nrb, k), dp)*AA(ma4(nra, nrb, k) + 1, xx)*BB(mb4(nra, nrb, k), yy)
1182 0 : w2 = w2 + REAL(c4(nra, nrb, k), dp)*AA(ma4(nra, nrb, k), xx)*BB(mb4(nra, nrb, k) + 1, yy)
1183 : END DO
1184 0 : Jss = f*w
1185 0 : dJss = f*REAL(nra + nrb + 1, dp)*w/rr
1186 0 : dJss = dJss - dxx*f*w1
1187 0 : dJss = dJss - dyy*f*w2
1188 :
1189 0 : w = 0.0_dp
1190 0 : w1 = 0.0_dp
1191 0 : w2 = 0.0_dp
1192 0 : f = rr**(nra + nrb + 1)/2.0_dp**(nra + nrb + 2)
1193 0 : DO k = 1, nc5(nra, nrb)
1194 0 : w = w + REAL(c5(nra, nrb, k), dp)*AA(ma5(nra, nrb, k), xx)*BB(mb5(nra, nrb, k), yy)
1195 0 : w1 = w1 + REAL(c5(nra, nrb, k), dp)*AA(ma5(nra, nrb, k) + 1, xx)*BB(mb5(nra, nrb, k), yy)
1196 0 : w2 = w2 + REAL(c5(nra, nrb, k), dp)*AA(ma5(nra, nrb, k), xx)*BB(mb5(nra, nrb, k) + 1, yy)
1197 : END DO
1198 0 : Jcc = f*w
1199 0 : dJcc = f*REAL(nra + nrb + 1, dp)*w/rr
1200 0 : dJcc = dJcc - dxx*f*w1
1201 0 : dJcc = dJcc - dyy*f*w2
1202 :
1203 0 : DO k1 = 1, 3
1204 0 : DO k2 = 1, 3
1205 : dS(k1 + 1, k2 + 1, :) = dS(k1 + 1, k2 + 1, :) &
1206 : & + 1.5_dp*Arot(k1, 1)*Arot(k2, 1)*fac1*fac2*dJss*v(:)/rr &
1207 : & + 1.5_dp*dArot(k1, 1, :)*Arot(k2, 1)*fac1*fac2*Jss &
1208 : & + 1.5_dp*Arot(k1, 1)*dArot(k2, 1, :)*fac1*fac2*Jss &
1209 : & + 1.5_dp*Arot(k1, 2)*Arot(k2, 2)*fac1*fac2*dJss*v(:)/rr &
1210 : & + 1.5_dp*dArot(k1, 2, :)*Arot(k2, 2)*fac1*fac2*Jss &
1211 : & + 1.5_dp*Arot(k1, 2)*dArot(k2, 2, :)*fac1*fac2*Jss &
1212 : & - 3.0_dp*Arot(k1, 3)*Arot(k2, 3)*fac1*fac2*dJcc*v(:)/rr &
1213 : & - 3.0_dp*dArot(k1, 3, :)*Arot(k2, 3)*fac1*fac2*Jcc &
1214 0 : & - 3.0_dp*Arot(k1, 3)*dArot(k2, 3, :)*fac1*fac2*Jcc
1215 : END DO
1216 : END DO
1217 :
1218 : END IF
1219 :
1220 : END SUBROUTINE makedS
1221 :
1222 : ! **************************************************************************************************
1223 : !> \brief ...
1224 : !> \param n ...
1225 : !> \param x ...
1226 : !> \return ...
1227 : ! **************************************************************************************************
1228 0 : FUNCTION AA(n, x)
1229 :
1230 : INTEGER :: n
1231 : REAL(kind=dp) :: x, AA
1232 :
1233 : REAL(kind=dp) :: p
1234 :
1235 0 : IF (n == 0) THEN
1236 : p = 1.0_dp
1237 : ELSE
1238 : IF (n == 1) THEN
1239 0 : p = 1.0_dp + x
1240 : ELSE
1241 : IF (n == 2) THEN
1242 : p = 2.0_dp + x*( &
1243 0 : 2.0_dp + x)
1244 : ELSE
1245 : IF (n == 3) THEN
1246 : p = 6.0_dp + x*( &
1247 : 6.0_dp + x*( &
1248 0 : 3.0_dp + x))
1249 : ELSE
1250 : IF (n == 4) THEN
1251 : p = 24.0_dp + x*( &
1252 : 24.0_dp + x*( &
1253 : 12.0_dp + x*( &
1254 0 : 4.0_dp + x)))
1255 : ELSE
1256 : IF (n == 5) THEN
1257 : p = 120.0_dp + x*( &
1258 : 120.0_dp + x*( &
1259 : 60.0_dp + x*( &
1260 : 20.0_dp + x*( &
1261 0 : 5.0_dp + x))))
1262 : ELSE
1263 : IF (n == 6) THEN
1264 : p = 720.0_dp + x*( &
1265 : 720.0_dp + x*( &
1266 : 360.0_dp + x*( &
1267 : 120.0_dp + x*( &
1268 : 30.0_dp + x*( &
1269 0 : 6.0_dp + x)))))
1270 : ELSE
1271 : IF (n == 7) THEN
1272 : p = 5040.0_dp + x*( &
1273 : 5040.0_dp + x*( &
1274 : 2520.0_dp + x*( &
1275 : 840.0_dp + x*( &
1276 : 210.0_dp + x*( &
1277 : 42.0_dp + x*( &
1278 0 : 7.0_dp + x))))))
1279 : ELSE
1280 : IF (n == 8) THEN
1281 : p = 40320.0_dp + x*( &
1282 : 40320.0_dp + x*( &
1283 : 20160.0_dp + x*( &
1284 : 6720.0_dp + x*( &
1285 : 1680.0_dp + x*( &
1286 : 336.0_dp + x*( &
1287 : 56.0_dp + x*( &
1288 0 : 8.0_dp + x)))))))
1289 : ELSE
1290 : IF (n == 9) THEN
1291 : p = 362880.0_dp + x*( &
1292 : 362880.0_dp + x*( &
1293 : 181440.0_dp + x*( &
1294 : 60480.0_dp + x*( &
1295 : 15120.0_dp + x*( &
1296 : 3024.0_dp + x*( &
1297 : 504.0_dp + x*( &
1298 : 72.0_dp + x*( &
1299 0 : 9.0_dp + x))))))))
1300 : ELSE
1301 : IF (n == 10) THEN
1302 : p = 3628800.0_dp + x*( &
1303 : 3628800.0_dp + x*( &
1304 : 1814400.0_dp + x*( &
1305 : 604800.0_dp + x*( &
1306 : 151200.0_dp + x*( &
1307 : 30240.0_dp + x*( &
1308 : 5040.0_dp + x*( &
1309 : 720.0_dp + x*( &
1310 : 90.0_dp + x*( &
1311 0 : 10.0_dp + x)))))))))
1312 : ELSE
1313 0 : p = 1.0_dp
1314 0 : WRITE (*, *) ' n= ', n, ' in AA(n,x) '
1315 : END IF
1316 : END IF
1317 : END IF
1318 : END IF
1319 : END IF
1320 : END IF
1321 : END IF
1322 : END IF
1323 : END IF
1324 : END IF
1325 : END IF
1326 :
1327 0 : AA = EXP(-x)*p/x**(n + 1)
1328 :
1329 0 : END FUNCTION AA
1330 :
1331 : ! **************************************************************************************************
1332 : !> \brief ...
1333 : !> \param n ...
1334 : !> \param y ...
1335 : !> \return ...
1336 : ! **************************************************************************************************
1337 0 : FUNCTION BB(n, y)
1338 :
1339 : INTEGER :: n
1340 : REAL(kind=dp) :: y, BB
1341 :
1342 0 : IF (ABS(y) > 1.0e-20_dp) THEN
1343 0 : BB = REAL((-1)**(n + 1), dp)*AA(n, -y) - AA(n, y)
1344 : ELSE
1345 0 : IF (MOD(n, 2) == 0) THEN
1346 0 : BB = 2.0_dp/REAL(n + 1, dp)
1347 : ELSE
1348 0 : BB = -y*2.0_dp/REAL(n + 2, dp)
1349 : END IF
1350 : END IF
1351 :
1352 0 : END FUNCTION BB
1353 :
1354 : END MODULE se_core_matrix
1355 :
|
