(git:374b731)
Loading...
Searching...
No Matches
se_core_matrix.F
Go to the documentation of this file.
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! **************************************************************************************************
13MODULE 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_operations, ONLY: dbcsr_allocate_matrix_set
19 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
20 USE cp_log_handling, ONLY: cp_get_default_logger,&
21 cp_logger_type
22 USE cp_output_handling, ONLY: cp_p_file,&
23 cp_print_key_finished_output,&
24 cp_print_key_should_output,&
25 cp_print_key_unit_nr
26 USE dbcsr_api, ONLY: &
27 dbcsr_add, dbcsr_copy, dbcsr_deallocate_matrix, dbcsr_distribute, dbcsr_get_block_diag, &
28 dbcsr_get_block_p, dbcsr_p_type, dbcsr_replicate_all, dbcsr_set, dbcsr_sum_replicated, &
29 dbcsr_type
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
71CONTAINS
72
73! **************************************************************************************************
74!> \brief ...
75!> \param qs_env ...
76!> \param para_env ...
77!> \param calculate_forces ...
78! **************************************************************************************************
79 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, nrt
90 LOGICAL :: defined, found, omit_headers, use_virial
91 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: zpt, zst
95 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 REAL(kind=dp), DIMENSION(:), POINTER :: beta_a, sto_exponents_a
99 REAL(kind=dp), DIMENSION(:, :), POINTER :: dsmat, h_block, h_blocka, pabmat, pamat, &
100 s_block
101 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
102 TYPE(cp_logger_type), POINTER :: logger
103 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 DIMENSION(:), POINTER :: nl_iterator
108 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
109 POINTER :: sab_orb
110 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
111 TYPE(qs_energy_type), POINTER :: energy
112 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
113 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 CALL timeset(routinen, handle)
122
123 NULLIFY (logger, energy)
124 logger => cp_get_default_logger()
125
126 NULLIFY (rho, force, atomic_kind_set, qs_kind_set, sab_orb, &
127 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 sab_orb=sab_orb)
142
143 ! calculate overlap matrix
144 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 sab_nl=sab_orb)
150 CALL set_ks_env(ks_env, matrix_s=matrix_s)
151 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 sab_nl=sab_orb)
158 CALL set_ks_env(ks_env, matrix_s=matrix_s)
159 use_virial = .false.
160 END IF
161
162 IF (calculate_forces) THEN
163 CALL qs_rho_get(rho, rho_ao=matrix_p)
164
165 IF (SIZE(matrix_p) == 2) THEN
166 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 delta = dft_control%qs_control%se_control%delta
169 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
170 atom_of_kind=atom_of_kind)
171 ALLOCATE (diagmat_p)
172 CALL dbcsr_get_block_diag(matrix_p(1)%matrix, diagmat_p)
173 CALL dbcsr_replicate_all(diagmat_p)
174 END IF
175
176 ! Allocate the core Hamiltonian matrix
177 CALL dbcsr_allocate_matrix_set(matrix_h, 1)
178 ALLOCATE (matrix_h(1)%matrix)
179 CALL dbcsr_copy(matrix_h(1)%matrix, matrix_s(1)%matrix, "CORE HAMILTONIAN MATRIX")
180 CALL dbcsr_set(matrix_h(1)%matrix, 0.0_dp)
181
182 ! Allocate a diagonal block matrix
183 ALLOCATE (diagmat_h)
184 CALL dbcsr_get_block_diag(matrix_s(1)%matrix, diagmat_h)
185 CALL dbcsr_set(diagmat_h, 0.0_dp)
186 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 kh = 0.5_dp
191
192 nkind = SIZE(atomic_kind_set)
193
194 ALLOCATE (se_defined(nkind))
195 ALLOCATE (hmt(16, nkind))
196 ALLOCATE (umt(16, nkind))
197
198 ALLOCATE (zst(nkind))
199 ALLOCATE (zpt(nkind))
200 ALLOCATE (nrt(nkind))
201
202 econst = 0.0_dp
203 DO ikind = 1, nkind
204 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom)
205 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 nr=nr_a, sto_exponents=sto_exponents_a)
209 econst = econst - (eisol - eheat)*real(natom, dp)
210 se_defined(ikind) = (defined .AND. natorb_a >= 1)
211 hmt(1, ikind) = beta_a(0)
212 hmt(2:4, ikind) = beta_a(1)
213 hmt(5:9, ikind) = beta_a(2)
214 hmt(10:16, ikind) = beta_a(3)
215 umt(1, ikind) = uss
216 umt(2:4, ikind) = upp
217 umt(5:9, ikind) = udd
218 umt(10:16, ikind) = uff
219
220 zst(ikind) = sto_exponents_a(0)
221 zpt(ikind) = sto_exponents_a(1)
222 nrt(ikind) = nr_a
223
224 END DO
225 energy%core_self = econst
226
227 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
228 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
229 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, inode=inode, r=rij)
230 IF (.NOT. se_defined(ikind)) cycle
231 IF (.NOT. se_defined(jkind)) cycle
232 ha(1:16) = hmt(1:16, ikind)
233 ua(1:16) = umt(1:16, ikind)
234 hb(1:16) = hmt(1:16, jkind)
235
236 nra = nrt(ikind)
237 nrb = nrt(jkind)
238 zsa = zst(ikind)
239 zsb = zst(jkind)
240 zpa = zpt(ikind)
241 zpb = zpt(jkind)
242
243 IF (inode == 1) THEN
244 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 NULLIFY (h_blocka)
248 CALL dbcsr_get_block_p(diagmat_h, iatom, iatom, h_blocka, found)
249 cpassert(ASSOCIATED(h_blocka))
250 IF (calculate_forces) THEN
251 CALL dbcsr_get_block_p(diagmat_p, iatom, iatom, pamat, found)
252 cpassert(ASSOCIATED(pamat))
253 END IF
254 END SELECT
255 END IF
256 dr = sum(rij(:)**2)
257 IF (iatom == jatom .AND. dr < rij_threshold) THEN
258
259 SELECT CASE (dft_control%qs_control%method_id)
260 CASE DEFAULT
261 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 DO i = 1, SIZE(h_blocka, 1)
265 h_blocka(i, i) = h_blocka(i, i) + ua(i)
266 END DO
267 END SELECT
268
269 ELSE
270 IF (iatom <= jatom) THEN
271 irow = iatom
272 icol = jatom
273 ELSE
274 irow = jatom
275 icol = iatom
276 END IF
277 NULLIFY (h_block)
278 CALL dbcsr_get_block_p(matrix_h(1)%matrix, &
279 irow, icol, h_block, found)
280 cpassert(ASSOCIATED(h_block))
281 ! two-centre one-electron term
282 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 irow, icol, s_block, found)
304 cpassert(ASSOCIATED(s_block))
305 IF (irow == iatom) THEN
306 DO i = 1, SIZE(h_block, 1)
307 DO j = 1, SIZE(h_block, 2)
308 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 DO i = 1, SIZE(h_block, 1)
313 DO j = 1, SIZE(h_block, 2)
314 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 IF (calculate_forces) THEN
319 atom_a = atom_of_kind(iatom)
320 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 CALL dbcsr_get_block_p(matrix_p(1)%matrix, irow, icol, pabmat, found)
331 cpassert(ASSOCIATED(pabmat))
332 DO icor = 1, 3
333 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 CALL dbcsr_get_block_p(matrix_s(icor + 1)%matrix, irow, icol, dsmat, found)
345 cpassert(ASSOCIATED(dsmat))
346 dsmat = 2._dp*kh*dsmat*pabmat
347 IF (irow == iatom) THEN
348 DO i = 1, SIZE(h_block, 1)
349 DO j = 1, SIZE(h_block, 2)
350 force_ab(icor) = force_ab(icor) + (ha(i) + hb(j))*dsmat(i, j)
351 END DO
352 END DO
353 ELSE
354 DO i = 1, SIZE(h_block, 1)
355 DO j = 1, SIZE(h_block, 2)
356 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 IF (calculate_forces .AND. (iatom /= jatom .OR. dr > rij_threshold)) THEN
366 IF (irow == iatom) force_ab = -force_ab
367 force(ikind)%all_potential(:, atom_a) = &
368 force(ikind)%all_potential(:, atom_a) - force_ab(:)
369 force(jkind)%all_potential(:, atom_b) = &
370 force(jkind)%all_potential(:, atom_b) + force_ab(:)
371 IF (use_virial) THEN
372 CALL virial_pair_force(virial%pv_virial, -1.0_dp, force_ab, rij)
373 END IF
374 END IF
375
376 END DO
377 CALL neighbor_list_iterator_release(nl_iterator)
378
379 DEALLOCATE (se_defined, hmt, umt, zst, zpt, nrt)
380
381 CALL dbcsr_sum_replicated(diagmat_h)
382 CALL dbcsr_distribute(diagmat_h)
383 CALL dbcsr_add(matrix_h(1)%matrix, diagmat_h, 1.0_dp, 1.0_dp)
384 CALL set_ks_env(ks_env, matrix_h=matrix_h)
385
386 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 extension=".Log")
390 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
391 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
392 after = min(max(after, 1), 16)
393 CALL cp_dbcsr_write_sparse_matrix(matrix_h(1)%matrix, 4, after, qs_env, para_env, &
394 scale=evolt, output_unit=iw, omit_headers=omit_headers)
395 CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
396 "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN")
397 END IF
398
399 IF (calculate_forces) THEN
400 IF (SIZE(matrix_p) == 2) THEN
401 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 DEALLOCATE (atom_of_kind)
404 CALL dbcsr_deallocate_matrix(diagmat_p)
405 END IF
406
407 CALL dbcsr_deallocate_matrix(diagmat_h)
408
409 CALL timestop(handle)
410
411 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 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 s(:, :) = 0.0_dp
629
630 v(:) = r(:)
631 rr = sqrt(dot_product(v, v))
632
633 IF (rr < 1.0e-20_dp) THEN
634
635 DO mu = 1, 4
636 s(mu, mu) = 1.0_dp
637 END DO
638
639 ELSE
640
641 fac1 = 1.0_dp
642 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 WRITE (*, *) 'nra= ', nra
655 RETURN
656 END IF
657 END IF
658 END IF
659 END IF
660 IF (nrb == 1) THEN
661 fac1 = fac1*2.0_dp
662 ELSE
663 IF (nrb == 2) THEN
664 fac1 = fac1*sqrt(4.0_dp/3.0_dp)
665 ELSE
666 IF (nrb == 3) THEN
667 fac1 = fac1*sqrt(8.0_dp/45.0_dp)
668 ELSE
669 IF (nrb == 4) THEN
670 fac1 = fac1*sqrt(4.0_dp/315.0_dp)
671 ELSE
672 WRITE (*, *) 'nrb= ', nrb
673 RETURN
674 END IF
675 END IF
676 END IF
677 END IF
678
679 ct = -v(3)/rr
680 IF (abs(ct) < 1.0_dp) THEN
681 st = sqrt(1.0_dp - ct**2)
682 cp = -v(1)/(rr*st)
683 sp = -v(2)/(rr*st)
684 arot(1, 1) = ct*cp
685 arot(1, 2) = -sp
686 arot(1, 3) = st*cp
687 arot(2, 1) = ct*sp
688 arot(2, 2) = cp
689 arot(2, 3) = st*sp
690 arot(3, 1) = -st
691 arot(3, 2) = 0.0_dp
692 arot(3, 3) = ct
693 ELSE
694 arot(1, 1) = ct
695 arot(1, 2) = 0.0_dp
696 arot(1, 3) = 0.0_dp
697 arot(2, 1) = 0.0_dp
698 arot(2, 2) = 1.0_dp
699 arot(2, 3) = 0.0_dp
700 arot(3, 1) = 0.0_dp
701 arot(3, 2) = 0.0_dp
702 arot(3, 3) = ct
703 END IF
704
705 za = zsa
706 zb = zsb
707 fac2 = sqrt(za**(2*nra + 1)*zb**(2*nrb + 1))
708 xx = 0.5_dp*rr*(za + zb)
709 yy = 0.5_dp*rr*(za - zb)
710
711 j = 0.0_dp
712 DO k = 1, nc1(nra, nrb)
713 j = j + real(c1(nra, nrb, k), dp)*aa(ma1(nra, nrb, k), xx)*bb(mb1(nra, nrb, k), yy)
714 END DO
715 j = j*rr**(nra + nrb + 1)
716 j = j/2.0_dp**(nra + nrb + 2)
717
718 s(1, 1) = s(1, 1) + fac1*fac2*j
719
720 za = zpa
721 zb = zsb
722 fac2 = sqrt(za**(2*nra + 1)*zb**(2*nrb + 1))
723 xx = 0.5_dp*rr*(za + zb)
724 yy = 0.5_dp*rr*(za - zb)
725
726 jc = 0.0_dp
727 DO k = 1, nc2(nra, nrb)
728 jc = jc + real(c2(nra, nrb, k), dp)*aa(ma2(nra, nrb, k), xx)*bb(mb2(nra, nrb, k), yy)
729 END DO
730 jc = jc*rr**(nra + nrb + 1)
731 jc = jc/2.0_dp**(nra + nrb + 2)
732
733 DO k1 = 1, 3
734 s(k1 + 1, 1) = s(k1 + 1, 1) &
735 & + sqrt(3.0_dp)*arot(k1, 3)*fac1*fac2*jc
736 END DO
737
738 za = zsa
739 zb = zpb
740 fac2 = sqrt(za**(2*nra + 1)*zb**(2*nrb + 1))
741 xx = 0.5_dp*rr*(za + zb)
742 yy = 0.5_dp*rr*(za - zb)
743
744 jc = 0.0_dp
745 DO k = 1, nc3(nra, nrb)
746 jc = jc + real(c3(nra, nrb, k), dp)*aa(ma3(nra, nrb, k), xx)*bb(mb3(nra, nrb, k), yy)
747 END DO
748 jc = jc*rr**(nra + nrb + 1)
749 jc = jc/2.0_dp**(nra + nrb + 2)
750
751 DO k1 = 1, 3
752 s(1, k1 + 1) = s(1, k1 + 1) &
753 & - sqrt(3.0_dp)*arot(k1, 3)*fac1*fac2*jc
754 END DO
755
756 za = zpa
757 zb = zpb
758 fac2 = sqrt(za**(2*nra + 1)*zb**(2*nrb + 1))
759 xx = 0.5_dp*rr*(za + zb)
760 yy = 0.5_dp*rr*(za - zb)
761
762 jss = 0.0_dp
763 DO k = 1, nc4(nra, nrb)
764 jss = jss + real(c4(nra, nrb, k), dp)*aa(ma4(nra, nrb, k), xx)*bb(mb4(nra, nrb, k), yy)
765 END DO
766 jss = jss*rr**(nra + nrb + 1)
767 jss = jss/2.0_dp**(nra + nrb + 2)
768
769 jcc = 0.0_dp
770 DO k = 1, nc5(nra, nrb)
771 jcc = jcc + real(c5(nra, nrb, k), dp)*aa(ma5(nra, nrb, k), xx)*bb(mb5(nra, nrb, k), yy)
772 END DO
773 jcc = jcc*rr**(nra + nrb + 1)
774 jcc = jcc/2.0_dp**(nra + nrb + 2)
775
776 DO k1 = 1, 3
777 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 & - 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 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 ds(:, :, :) = 0.0_dp
1007
1008 v(:) = r(:)
1009 rr = sqrt(dot_product(v, v))
1010
1011 IF (rr < 1.0e-20_dp) THEN
1012
1013 DO mu = 1, 4
1014 ds(mu, mu, :) = 0.0_dp
1015 END DO
1016
1017 ELSE
1018
1019 fac1 = 1.0_dp
1020 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 WRITE (*, *) 'nra= ', nra
1033 RETURN
1034 END IF
1035 END IF
1036 END IF
1037 END IF
1038 IF (nrb == 1) THEN
1039 fac1 = fac1*2.0_dp
1040 ELSE
1041 IF (nrb == 2) THEN
1042 fac1 = fac1*sqrt(4.0_dp/3.0_dp)
1043 ELSE
1044 IF (nrb == 3) THEN
1045 fac1 = fac1*sqrt(8.0_dp/45.0_dp)
1046 ELSE
1047 IF (nrb == 4) THEN
1048 fac1 = fac1*sqrt(4.0_dp/315.0_dp)
1049 ELSE
1050 WRITE (*, *) 'nrb= ', nrb
1051 RETURN
1052 END IF
1053 END IF
1054 END IF
1055 END IF
1056
1057 ct = -v(3)/rr
1058 IF (abs(ct) >= 1.0_dp) THEN
1059
1060 dct(:) = v(:)*v(3)/rr**3
1061 dct(3) = dct(3) - 1.0_dp/rr
1062
1063 arot(1, 1) = ct
1064 arot(1, 2) = 0.0_dp
1065 arot(1, 3) = 0.0_dp
1066 arot(2, 1) = 0.0_dp
1067 arot(2, 2) = 1.0_dp
1068 arot(2, 3) = 0.0_dp
1069 arot(3, 1) = 0.0_dp
1070 arot(3, 2) = 0.0_dp
1071 arot(3, 3) = ct
1072
1073 darot(1, 1, :) = dct(:)
1074 darot(1, 2, :) = 0.0_dp
1075 darot(1, 3, :) = 0.0_dp
1076 darot(2, 1, :) = 0.0_dp
1077 darot(2, 2, :) = 0.0_dp
1078 darot(2, 3, :) = 0.0_dp
1079 darot(3, 1, :) = 0.0_dp
1080 darot(3, 2, :) = 0.0_dp
1081 darot(3, 3, :) = dct(:)
1082
1083 ELSE
1084
1085 xx = sqrt(v(1)**2 + v(2)**2)
1086 st = xx/rr
1087 cp = -v(1)/xx
1088 sp = -v(2)/xx
1089
1090 dct(:) = v(:)*v(3)/rr**3
1091 dct(3) = dct(3) - 1.0_dp/rr
1092 dst(:) = -ct*dct(:)/st
1093 dcp(:) = v(:)*v(1)/(rr**3*st)
1094 dcp(:) = dcp(:) + v(1)*dst(:)/(rr*st**2)
1095 dcp(1) = dcp(1) - 1.0_dp/(rr*st)
1096 dsp(:) = v(:)*v(2)/(rr**3*st)
1097 dsp(:) = dsp(:) + v(2)*dst(:)/(rr*st**2)
1098 dsp(2) = dsp(2) - 1.0_dp/(rr*st)
1099
1100 arot(1, 1) = ct*cp
1101 arot(1, 2) = -sp
1102 arot(1, 3) = st*cp
1103 arot(2, 1) = ct*sp
1104 arot(2, 2) = cp
1105 arot(2, 3) = st*sp
1106 arot(3, 1) = -st
1107 arot(3, 2) = 0.0_dp
1108 arot(3, 3) = ct
1109
1110 darot(1, 1, :) = dct(:)*cp + ct*dcp(:)
1111 darot(1, 2, :) = -dsp(:)
1112 darot(1, 3, :) = dst(:)*cp + st*dcp(:)
1113 darot(2, 1, :) = dct(:)*sp + ct*dsp(:)
1114 darot(2, 2, :) = dcp(:)
1115 darot(2, 3, :) = dst(:)*sp + st*dsp(:)
1116 darot(3, 1, :) = -dst(:)
1117 darot(3, 2, :) = 0.0_dp
1118 darot(3, 3, :) = dct(:)
1119
1120 END IF
1121
1122 za = zsa
1123 zb = zsb
1124 fac2 = sqrt(za**(2*nra + 1)*zb**(2*nrb + 1))
1125 xx = 0.5_dp*rr*(za + zb)
1126 yy = 0.5_dp*rr*(za - zb)
1127 dxx = 0.5_dp*(za + zb)
1128 dyy = 0.5_dp*(za - zb)
1129
1130 w = 0.0_dp
1131 w1 = 0.0_dp
1132 w2 = 0.0_dp
1133 f = rr**(nra + nrb + 1)/2.0_dp**(nra + nrb + 2)
1134 DO k = 1, nc1(nra, nrb)
1135 w = w + real(c1(nra, nrb, k), dp)*aa(ma1(nra, nrb, k), xx)*bb(mb1(nra, nrb, k), yy)
1136 w1 = w1 + real(c1(nra, nrb, k), dp)*aa(ma1(nra, nrb, k) + 1, xx)*bb(mb1(nra, nrb, k), yy)
1137 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 j = f*w
1140 dj = f*real(nra + nrb + 1, dp)*w/rr
1141 dj = dj - dxx*f*w1
1142 dj = dj - dyy*f*w2
1143
1144 ds(1, 1, :) = ds(1, 1, :) + fac1*fac2*dj*v(:)/rr
1145
1146 za = zpa
1147 zb = zsb
1148 fac2 = sqrt(za**(2*nra + 1)*zb**(2*nrb + 1))
1149 xx = 0.5_dp*rr*(za + zb)
1150 yy = 0.5_dp*rr*(za - zb)
1151 dxx = 0.5_dp*(za + zb)
1152 dyy = 0.5_dp*(za - zb)
1153
1154 w = 0.0_dp
1155 w1 = 0.0_dp
1156 w2 = 0.0_dp
1157 f = rr**(nra + nrb + 1)/2.0_dp**(nra + nrb + 2)
1158 DO k = 1, nc2(nra, nrb)
1159 w = w + real(c2(nra, nrb, k), dp)*aa(ma2(nra, nrb, k), xx)*bb(mb2(nra, nrb, k), yy)
1160 w1 = w1 + real(c2(nra, nrb, k), dp)*aa(ma2(nra, nrb, k) + 1, xx)*bb(mb2(nra, nrb, k), yy)
1161 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 jc = f*w
1164 djc = f*real(nra + nrb + 1, dp)*w/rr
1165 djc = djc - dxx*f*w1
1166 djc = djc - dyy*f*w2
1167
1168 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 & + sqrt(3.0_dp)*darot(k1, 3, :)*fac1*fac2*jc
1172 END DO
1173
1174 za = zsa
1175 zb = zpb
1176 fac2 = sqrt(za**(2*nra + 1)*zb**(2*nrb + 1))
1177 xx = 0.5_dp*rr*(za + zb)
1178 yy = 0.5_dp*rr*(za - zb)
1179 dxx = 0.5_dp*(za + zb)
1180 dyy = 0.5_dp*(za - zb)
1181
1182 w = 0.0_dp
1183 w1 = 0.0_dp
1184 w2 = 0.0_dp
1185 f = rr**(nra + nrb + 1)/2.0_dp**(nra + nrb + 2)
1186 DO k = 1, nc3(nra, nrb)
1187 w = w + real(c3(nra, nrb, k), dp)*aa(ma3(nra, nrb, k), xx)*bb(mb3(nra, nrb, k), yy)
1188 w1 = w1 + real(c3(nra, nrb, k), dp)*aa(ma3(nra, nrb, k) + 1, xx)*bb(mb3(nra, nrb, k), yy)
1189 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 jc = f*w
1192 djc = f*real(nra + nrb + 1, dp)*w/rr
1193 djc = djc - dxx*f*w1
1194 djc = djc - dyy*f*w2
1195
1196 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 & - sqrt(3.0_dp)*darot(k1, 3, :)*fac1*fac2*jc
1200 END DO
1201
1202 za = zpa
1203 zb = zpb
1204 fac2 = sqrt(za**(2*nra + 1)*zb**(2*nrb + 1))
1205 xx = 0.5_dp*rr*(za + zb)
1206 yy = 0.5_dp*rr*(za - zb)
1207 dxx = 0.5_dp*(za + zb)
1208 dyy = 0.5_dp*(za - zb)
1209
1210 w = 0.0_dp
1211 w1 = 0.0_dp
1212 w2 = 0.0_dp
1213 f = rr**(nra + nrb + 1)/2.0_dp**(nra + nrb + 2)
1214 DO k = 1, nc4(nra, nrb)
1215 w = w + real(c4(nra, nrb, k), dp)*aa(ma4(nra, nrb, k), xx)*bb(mb4(nra, nrb, k), yy)
1216 w1 = w1 + real(c4(nra, nrb, k), dp)*aa(ma4(nra, nrb, k) + 1, xx)*bb(mb4(nra, nrb, k), yy)
1217 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 jss = f*w
1220 djss = f*real(nra + nrb + 1, dp)*w/rr
1221 djss = djss - dxx*f*w1
1222 djss = djss - dyy*f*w2
1223
1224 w = 0.0_dp
1225 w1 = 0.0_dp
1226 w2 = 0.0_dp
1227 f = rr**(nra + nrb + 1)/2.0_dp**(nra + nrb + 2)
1228 DO k = 1, nc5(nra, nrb)
1229 w = w + real(c5(nra, nrb, k), dp)*aa(ma5(nra, nrb, k), xx)*bb(mb5(nra, nrb, k), yy)
1230 w1 = w1 + real(c5(nra, nrb, k), dp)*aa(ma5(nra, nrb, k) + 1, xx)*bb(mb5(nra, nrb, k), yy)
1231 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 jcc = f*w
1234 djcc = f*real(nra + nrb + 1, dp)*w/rr
1235 djcc = djcc - dxx*f*w1
1236 djcc = djcc - dyy*f*w2
1237
1238 DO k1 = 1, 3
1239 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 & - 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 FUNCTION aa(n, x)
1264
1265 INTEGER :: n
1266 REAL(kind=dp) :: x, aa
1267
1268 REAL(kind=dp) :: p
1269
1270 IF (n == 0) THEN
1271 p = 1.0_dp
1272 ELSE
1273 IF (n == 1) THEN
1274 p = 1.0_dp + x
1275 ELSE
1276 IF (n == 2) THEN
1277 p = 2.0_dp + x*( &
1278 2.0_dp + x)
1279 ELSE
1280 IF (n == 3) THEN
1281 p = 6.0_dp + x*( &
1282 6.0_dp + x*( &
1283 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 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 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 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 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 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 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 10.0_dp + x)))))))))
1347 ELSE
1348 p = 1.0_dp
1349 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 aa = exp(-x)*p/x**(n + 1)
1363
1364 END FUNCTION aa
1365
1366! **************************************************************************************************
1367!> \brief ...
1368!> \param n ...
1369!> \param y ...
1370!> \return ...
1371! **************************************************************************************************
1372 FUNCTION bb(n, y)
1373
1374 INTEGER :: n
1375 REAL(kind=dp) :: y, bb
1376
1377 IF (abs(y) > 1.0e-20_dp) THEN
1378 bb = real((-1)**(n + 1), dp)*aa(n, -y) - aa(n, y)
1379 ELSE
1380 IF (mod(n, 2) == 0) THEN
1381 bb = 2.0_dp/real(n + 1, dp)
1382 ELSE
1383 bb = -y*2.0_dp/real(n + 2, dp)
1384 END IF
1385 END IF
1386
1387 END FUNCTION bb
1388
1389END MODULE se_core_matrix
1390
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:81
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_output
DBCSR output in CP2K.
Definition cp_dbcsr_output.F:17
cp_dbcsr_output::cp_dbcsr_write_sparse_matrix
subroutine, public cp_dbcsr_write_sparse_matrix(sparse_matrix, before, after, qs_env, para_env, first_row, last_row, first_col, last_col, scale, output_unit, omit_headers)
...
Definition cp_dbcsr_output.F:163
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1013
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
dct
the type I Discrete Cosine Transform (DCT-I)
Definition dct.F:16
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::do_method_pdg
integer, parameter, public do_method_pdg
Definition input_constants.F:250
input_constants::do_method_pnnl
integer, parameter, public do_method_pnnl
Definition input_constants.F:250
input_constants::do_method_rm1
integer, parameter, public do_method_rm1
Definition input_constants.F:250
input_constants::do_method_pm3
integer, parameter, public do_method_pm3
Definition input_constants.F:250
input_constants::do_method_mndo
integer, parameter, public do_method_mndo
Definition input_constants.F:250
input_constants::do_method_mndod
integer, parameter, public do_method_mndod
Definition input_constants.F:250
input_constants::do_method_am1
integer, parameter, public do_method_am1
Definition input_constants.F:250
input_constants::do_method_pm6fm
integer, parameter, public do_method_pm6fm
Definition input_constants.F:250
input_constants::do_method_pm6
integer, parameter, public do_method_pm6
Definition input_constants.F:250
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1040
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::sp
integer, parameter, public sp
Definition kinds.F:33
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::evolt
real(kind=dp), parameter, public evolt
Definition physcon.F:183
qs_energy_types
Definition qs_energy_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:502
qs_force_types
Definition qs_force_types.F:13
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:451
qs_ks_types
Definition qs_ks_types.F:15
qs_ks_types::set_ks_env
subroutine, public set_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, kpoints, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
Definition qs_ks_types.F:567
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition qs_neighbor_list_types.F:165
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition qs_neighbor_list_types.F:251
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition qs_neighbor_list_types.F:549
qs_overlap
Calculation of overlap matrix, its derivatives and forces.
Definition qs_overlap.F:19
qs_overlap::build_overlap_matrix
subroutine, public build_overlap_matrix(ks_env, matrix_s, matrixkp_s, matrix_name, nderivative, basis_type_a, basis_type_b, sab_nl, calculate_forces, matrix_p, matrixkp_p)
Calculation of the overlap matrix over Cartesian Gaussian functions.
Definition qs_overlap.F:120
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
se_core_matrix
Calculation of the Hamiltonian integral matrix <a|H|b> for semi-empirical methods.
Definition se_core_matrix.F:13
se_core_matrix::build_se_core_matrix
subroutine, public build_se_core_matrix(qs_env, para_env, calculate_forces)
...
Definition se_core_matrix.F:80
semi_empirical_int_arrays
Arrays of parameters used in the semi-empirical calculations \References Everywhere in this module TC...
Definition semi_empirical_int_arrays.F:17
semi_empirical_int_arrays::rij_threshold
real(kind=dp), parameter, public rij_threshold
Definition semi_empirical_int_arrays.F:27
semi_empirical_types
Definition of the semi empirical parameter types.
Definition semi_empirical_types.F:12
semi_empirical_types::get_se_param
subroutine, public get_se_param(sep, name, typ, defined, z, zeff, natorb, eheat, beta, sto_exponents, uss, upp, udd, uff, alp, eisol, gss, gsp, gpp, gp2, acoul, nr, de, ass, asp, app, hsp, gsd, gpd, gdd, ppddg, dpddg, ngauss)
Get info from the semi-empirical type.
Definition semi_empirical_types.F:363
semi_empirical_utils
Working with the semi empirical parameter types.
Definition semi_empirical_utils.F:12
semi_empirical_utils::get_se_type
integer function, public get_se_type(se_method)
Gives back the unique semi_empirical METHOD type.
Definition semi_empirical_utils.F:336
virial_methods
Definition virial_methods.F:12
virial_methods::virial_pair_force
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
Definition virial_methods.F:139
virial_types
Definition virial_types.F:13
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cp_control_types::dft_control_type
Definition cp_control_types.F:559
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:215
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:135
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
semi_empirical_types::semi_empirical_type
Semi-empirical type.
Definition semi_empirical_types.F:57
virial_types::virial_type
Definition virial_types.F:26