(git:ccc2433)
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 ! **************************************************************************************************
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_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 
71 CONTAINS
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 
1389 END 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
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