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