(git:d3d49ac)
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-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Calculation of the Hamiltonian integral matrix <a|H|b> for
10!> semi-empirical methods
11!> \author JGH
12! **************************************************************************************************
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_api, ONLY: &
19 dbcsr_add, dbcsr_copy, dbcsr_deallocate_matrix, dbcsr_distribute, dbcsr_get_block_p, &
20 dbcsr_p_type, dbcsr_replicate_all, dbcsr_set, dbcsr_sum_replicated, dbcsr_type
21 USE cp_dbcsr_contrib, ONLY: dbcsr_get_block_diag
22 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
23 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
24 USE cp_log_handling, ONLY: cp_get_default_logger,&
25 cp_logger_type
26 USE cp_output_handling, ONLY: cp_p_file,&
27 cp_print_key_finished_output,&
28 cp_print_key_should_output,&
29 cp_print_key_unit_nr
30 USE input_constants, ONLY: &
31 do_method_am1, do_method_mndo, do_method_mndod, do_method_pdg, do_method_pm3, &
32 do_method_pm6, do_method_pm6fm, do_method_pnnl, do_method_rm1
33 USE input_section_types, ONLY: section_vals_val_get
34 USE kinds, ONLY: dp
35 USE message_passing, ONLY: mp_para_env_type
36 USE particle_types, ONLY: particle_type
37 USE physcon, ONLY: evolt
38 USE qs_energy_types, ONLY: qs_energy_type
39 USE qs_environment_types, ONLY: get_qs_env,&
40 qs_environment_type
41 USE qs_force_types, ONLY: qs_force_type
42 USE qs_kind_types, ONLY: get_qs_kind,&
43 qs_kind_type
44 USE qs_ks_types, ONLY: qs_ks_env_type,&
45 set_ks_env
46 USE qs_neighbor_list_types, ONLY: get_iterator_info,&
47 neighbor_list_iterate,&
48 neighbor_list_iterator_create,&
49 neighbor_list_iterator_p_type,&
50 neighbor_list_iterator_release,&
51 neighbor_list_set_p_type
52 USE qs_overlap, ONLY: build_overlap_matrix
53 USE qs_rho_types, ONLY: qs_rho_get,&
54 qs_rho_type
55 USE semi_empirical_int_arrays, ONLY: rij_threshold
56 USE semi_empirical_types, ONLY: get_se_param,&
57 semi_empirical_type
58 USE semi_empirical_utils, ONLY: get_se_type
59 USE virial_methods, ONLY: virial_pair_force
60 USE virial_types, ONLY: virial_type
61#include "./base/base_uses.f90"
62
63 IMPLICIT NONE
64
65 PRIVATE
66
67 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'se_core_matrix'
68
69 PUBLIC :: build_se_core_matrix
70
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("Unknown method for semi-empirical calculation")
262 CASE (do_method_am1, do_method_rm1, do_method_mndo, do_method_pdg, &
263 do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_mndod, do_method_pnnl)
264 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 cpassert(ASSOCIATED(s_block))
287 IF (irow == iatom) THEN
288 DO i = 1, SIZE(h_block, 1)
289 DO j = 1, SIZE(h_block, 2)
290 h_block(i, j) = h_block(i, j) + kh*(ha(i) + hb(j))*s_block(i, j)
291 END DO
292 END DO
293 ELSE
294 DO i = 1, SIZE(h_block, 1)
295 DO j = 1, SIZE(h_block, 2)
296 h_block(i, j) = h_block(i, j) + kh*(ha(j) + hb(i))*s_block(i, j)
297 END DO
298 END DO
299 END IF
300 IF (calculate_forces) THEN
301 atom_a = atom_of_kind(iatom)
302 atom_b = atom_of_kind(jatom)
303
304 CALL dbcsr_get_block_p(matrix_p(1)%matrix, irow, icol, pabmat, found)
305 cpassert(ASSOCIATED(pabmat))
306 DO icor = 1, 3
307 force_ab(icor) = 0._dp
308
309 CALL dbcsr_get_block_p(matrix_s(icor + 1)%matrix, irow, icol, dsmat, found)
310 cpassert(ASSOCIATED(dsmat))
311 dsmat = 2._dp*kh*dsmat*pabmat
312 IF (irow == iatom) THEN
313 DO i = 1, SIZE(h_block, 1)
314 DO j = 1, SIZE(h_block, 2)
315 force_ab(icor) = force_ab(icor) + (ha(i) + hb(j))*dsmat(i, j)
316 END DO
317 END DO
318 ELSE
319 DO i = 1, SIZE(h_block, 1)
320 DO j = 1, SIZE(h_block, 2)
321 force_ab(icor) = force_ab(icor) + (ha(j) + hb(i))*dsmat(i, j)
322 END DO
323 END DO
324 END IF
325 END DO
326 END IF
327
328 END IF
329
330 IF (calculate_forces .AND. (iatom /= jatom .OR. dr > rij_threshold)) THEN
331 IF (irow == iatom) force_ab = -force_ab
332 force(ikind)%all_potential(:, atom_a) = &
333 force(ikind)%all_potential(:, atom_a) - force_ab(:)
334 force(jkind)%all_potential(:, atom_b) = &
335 force(jkind)%all_potential(:, atom_b) + force_ab(:)
336 IF (use_virial) THEN
337 CALL virial_pair_force(virial%pv_virial, -1.0_dp, force_ab, rij)
338 END IF
339 END IF
340
341 END DO
342 CALL neighbor_list_iterator_release(nl_iterator)
343
344 DEALLOCATE (se_defined, hmt, umt, zst, zpt, nrt)
345
346 CALL dbcsr_sum_replicated(diagmat_h)
347 CALL dbcsr_distribute(diagmat_h)
348 CALL dbcsr_add(matrix_h(1)%matrix, diagmat_h, 1.0_dp, 1.0_dp)
349 CALL set_ks_env(ks_env, matrix_h=matrix_h)
350
351 IF (btest(cp_print_key_should_output(logger%iter_info, &
352 qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN"), cp_p_file)) THEN
353 iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN", &
354 extension=".Log")
355 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
356 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
357 after = min(max(after, 1), 16)
358 CALL cp_dbcsr_write_sparse_matrix(matrix_h(1)%matrix, 4, after, qs_env, para_env, &
359 scale=evolt, output_unit=iw, omit_headers=omit_headers)
360 CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
361 "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN")
362 END IF
363
364 IF (calculate_forces) THEN
365 IF (SIZE(matrix_p) == 2) THEN
366 CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
367 END IF
368 DEALLOCATE (atom_of_kind)
369 CALL dbcsr_deallocate_matrix(diagmat_p)
370 END IF
371
372 CALL dbcsr_deallocate_matrix(diagmat_h)
373
374 CALL timestop(handle)
375
376 END SUBROUTINE build_se_core_matrix
377
378! **************************************************************************************************
379!> \brief ...
380!> \param R ...
381!> \param nra ...
382!> \param nrb ...
383!> \param ZSA ...
384!> \param ZSB ...
385!> \param ZPA ...
386!> \param ZPB ...
387!> \param S ...
388! **************************************************************************************************
389 SUBROUTINE makes(R, nra, nrb, ZSA, ZSB, ZPA, ZPB, S)
390
391 REAL(kind=dp), DIMENSION(3) :: r
392 INTEGER :: nra, nrb
393 REAL(kind=dp) :: zsa, zsb, zpa, zpb
394 REAL(kind=dp), DIMENSION(4, 4) :: s
395
396 INTEGER, DIMENSION(4, 4), PARAMETER :: &
397 nc1 = reshape([2, 4, 4, 6, 4, 3, 6, 7, 4, 6, 4, 8, 6, 7, 8, 5], [4, 4]), &
398 nc2 = reshape([4, 4, 8, 8, 6, 8, 6, 12, 8, 8, 12, 8, 10, 12, 14, 16], [4, 4]), &
399 nc3 = reshape([4, 6, 8, 10, 4, 8, 8, 12, 8, 6, 12, 14, 8, 12, 8, 16], [4, 4]), &
400 nc4 = reshape([4, 8, 11, 14, 8, 6, 12, 14, 11, 12, 10, 20, 14, 14, 20, 12], [4, 4]), &
401 nc5 = reshape([2, 4, 6, 8, 4, 4, 8, 8, 6, 8, 6, 12, 8, 8, 12, 8], [4, 4])
402 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: c1 = reshape([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, &
403 1, 1, 1, 1, -1, 1, 2, 3, -1, -2, 1, 2, -2, -1, -3, 1, -3, -2, -1, -4, 0, -1, -2, 2, -1, 1,&
404 -2, -1, 2, -2, 3, -3, 2, -1, -3, 6, 0, -1, -1, -2, 1, 0, -2, -4, -1, 2, -1, -3, 2, 4, 3, &
405 -4, 0, 0, 0, -3, 0, 0, 1, -1, 0, 1, 0, 3, -3, -1, 3, 1, 0, 0, 0, -1, 0, 0, 1, 2, 0, -1, 0,&
406 3, 1, -2, -3, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, -1, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0 &
407 , 0, 0, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
408 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
409 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
410 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
411 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
412 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
413 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
414 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: c2 = reshape([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, &
415 1, 1, 1, 1, -1, 1, 1, 2, -2, -1, 1, 1, -3, -2, -1, 1, -4, -3, -2, -1, 1, -1, 1, 1, 1, 1, &
416 -2, 1, 1, 1, 1, -3, 1, 1, 1, 1, -1, -1, -1, 2, 1, -1, -2, -2, 3, -2, -2, -3, 6, 2, -1, -3,&
417 0, 0, 1, -2, -2, -1, 1, 1, -3, 2, -1, 3, -4, -3, -2, -1, 0, 0, -1, -1, 1, 1, 1, -2, -1, -1&
418 , 2, 3, -4, 2, 4, 3, 0, 0, -1, -2, 0, -1, 0, -2, 3, 2, -2, -1, 6, 2, -1, -3, 0, 0, -1, -1,&
419 0, 1, 0, 1, -1, -1, 1, -1, 1, -3, -1, 3, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 2, 0, -4, 2, 4, 3 &
420 , 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 1, 1, -2, -3, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0&
421 , -3, -1, 3, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 0, 1, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
422 0, 0, 0, 0, 0, -2, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0&
423 , 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, &
424 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
425 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],&
426 [4, 4, 20])
427 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: c3 = reshape([-1, -1, -1, -1, -1, -1, -1, -1, -1, &
428 -1, -1, -1, -1, -1, -1, -1, -1, -2, -3, -4, 1, -1, -2, -3, 1, 1, -1, -2, 2, 1, 1, -1, 1, 1&
429 , 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 3, 1, 1, -1, -3, -6, -1, 1, 2, -2, 1, -2, 2, 1, -2, &
430 2, -3, 3, 0, 2, 3, 4, 0, 1, 2, 3, -1, -1, 1, 2, -2, -1, -3, 1, 0, 1, -1, -4, 0, 1, 1, 2, &
431 -1, 1, 2, 4, 1, -2, 3, 3, 0, 0, 3, 6, 0, -1, -2, 2, -1, 0, -2, -1, 2, -2, 1, -3, 0, 0, 1, &
432 -1, 0, -1, -1, 3, 1, 0, -1, 1, -1, -1, -1, -3, 0, 0, 0, 4, 0, 0, 0, -2, 0, 0, -2, -4, 0, 2&
433 , 0, -3, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, -1, -2, 0, 1, 0, -3, 0, 0, 0, 0, 0, 0, 0, -3, 0, 0,&
434 1, -1, 0, 1, 0, 3, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, -1, 0, -1, 0, 1, 0, 0, 0, 0, 0, 0, 0,&
435 0, 0, 0, 0, 2, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, &
436 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0&
437 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
438 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
439 , 0], [4, 4, 20])
440 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: c4 = reshape([-1, -1, -1, -1, -1, -1, -1, -1, -1, &
441 -1, -1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -2, -3,&
442 1, 1, -1, -2, 2, 1, 2, -1, 3, 2, 1, 3, -1, 1, 2, 3, -1, -1, 1, 2, -2, -1, -1, 1, -3, -2, &
443 -1, -2, 0, 1, -1, -3, 1, -1, 1, 1, -1, 1, -1, 2, -3, 1, 2, -1, 0, -1, 2, 4, -1, 1, -1, -1 &
444 , 2, -1, -1, -1, 4, -1, -1, -3, 0, 1, -1, -1, -1, 0, 1, 2, -1, -1, -1, -1, -1, -2, -1, 3, &
445 0, -1, 2, -1, 1, 0, -1, -2, -2, 1, 2, 2, 1, 2, -2, 1, 0, 0, -2, 4, 0, 0, -1, 1, 2, -1, 1, &
446 -1, -4, 1, 1, 2, 0, 0, 1, -3, 0, 0, 1, -1, 1, 1, -1, -1, 3, -1, 1, -3, 0, 0, -1, 3, 0, 0, &
447 -1, -2, -1, 1, 0, -1, 3, 2, -1, -1, 0, 0, 0, -3, 0, 0, 1, 2, 0, -1, 0, -1, -3, -2, -1, 1, &
448 0, 0, 0, 1, 0, 0, 0, -1, 0, 0, 0, 2, -1, -1, 2, 0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, 0, -1, 1&
449 , 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
450 0, 2, 0, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,&
451 0, 0, 0, -1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 0, &
452 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0], [4, 4, 20])
453 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: c5 = reshape([-1, -1, -1, -1, -1, -1, -1, -1, -1, &
454 -1, -1, -1, -1, -1, -1, -1, 1, -1, -2, -3, 1, 1, -1, -2, 2, 1, 2, -1, 3, 2, 1, 3, 0, 1, -1&
455 , -3, 1, 1, 1, 1, -1, 1, 1, 2, -3, 1, 2, 1, 0, 1, 1, 1, -1, -1, 1, 2, 1, 1, -1, 1, 1, -2, &
456 1, -3, 0, 0, 2, -1, 0, 0, 1, 2, -2, -1, -2, 2, 1, -2, -2, -3, 0, 0, 1, 3, 0, 0, 1, 1, 1, &
457 -1, 1, 1, -3, 1, -1, 1, 0, 0, 0, 3, 0, 0, -1, -2, 0, -1, 0, -1, 3, 2, -1, 3, 0, 0, 0, 1, 0&
458 , 0, -1, -1, 0, 1, 0, -2, -1, -1, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0,&
459 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,&
460 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 &
461 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
462 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
463 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
464 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20&
465 ])
466 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: ma1 = reshape([2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7 &
467 , 5, 6, 7, 8, 0, 2, 3, 4, 2, 2, 4, 5, 3, 4, 4, 6, 4, 5, 6, 6, 0, 1, 1, 3, 1, 0, 3, 4, 1, 3&
468 , 2, 5, 3, 4, 5, 4, 0, 0, 0, 2, 0, 0, 2, 3, 0, 2, 0, 4, 2, 3, 4, 2, 0, 0, 0, 1, 0, 0, 1, 2&
469 , 0, 1, 0, 3, 1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 2, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0&
470 , 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
471 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
472 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
473 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
474 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
475 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
476 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
477 , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
478 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: ma2 = reshape([2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7 &
479 , 5, 6, 7, 8, 1, 2, 3, 4, 2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7, 1, 1, 3, 4, 2, 3, 3, 5, 3, 4&
480 , 5, 5, 4, 5, 6, 7, 0, 0, 2, 3, 1, 2, 2, 4, 2, 3, 4, 4, 3, 4, 5, 6, 0, 0, 2, 2, 1, 2, 1, 4&
481 , 2, 2, 4, 3, 3, 4, 5, 6, 0, 0, 1, 1, 0, 1, 0, 3, 1, 1, 3, 2, 2, 3, 4, 5, 0, 0, 1, 1, 0, 1&
482 , 0, 3, 1, 1, 3, 1, 2, 3, 4, 5, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 2, 0, 1, 2, 3, 4, 0, 0, 0, 0&
483 , 0, 0, 0, 2, 0, 0, 2, 0, 1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 2, 3, 0, 0&
484 , 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2&
485 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
486 , 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
487 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
488 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
489 , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
490 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: ma3 = reshape([2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7 &
491 , 5, 6, 7, 8, 1, 2, 3, 4, 2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7, 1, 2, 3, 4, 1, 3, 4, 5, 3, 3&
492 , 5, 6, 4, 5, 5, 7, 0, 1, 2, 3, 0, 2, 3, 4, 2, 2, 4, 5, 3, 4, 4, 6, 0, 1, 2, 3, 0, 2, 2, 4&
493 , 2, 1, 4, 5, 2, 4, 3, 6, 0, 0, 1, 2, 0, 1, 1, 3, 1, 0, 3, 4, 1, 3, 2, 5, 0, 0, 1, 2, 0, 1&
494 , 1, 3, 1, 0, 3, 4, 1, 3, 1, 5, 0, 0, 0, 1, 0, 0, 0, 2, 0, 0, 2, 3, 0, 2, 0, 4, 0, 0, 0, 1&
495 , 0, 0, 0, 2, 0, 0, 2, 3, 0, 2, 0, 4, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 2, 0, 1, 0, 3, 0, 0&
496 , 0, 0, 0, 0, 0, 1, 0, 0, 1, 2, 0, 1, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 2&
497 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
498 , 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
499 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
500 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
501 , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
502 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: ma4 = reshape([2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7 &
503 , 5, 6, 7, 8, 2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7, 5, 6, 7, 8, 0, 2, 3, 4, 2, 2, 4, 5, 3, 4&
504 , 4, 6, 4, 5, 6, 6, 0, 2, 3, 4, 2, 2, 4, 5, 3, 4, 4, 6, 4, 5, 6, 6, 0, 1, 2, 3, 1, 0, 3, 4&
505 , 2, 3, 4, 5, 3, 4, 5, 6, 0, 1, 2, 3, 1, 0, 3, 4, 2, 3, 2, 5, 3, 4, 5, 4, 0, 0, 2, 3, 0, 0&
506 , 2, 3, 2, 2, 2, 5, 3, 3, 5, 4, 0, 0, 1, 2, 0, 0, 2, 3, 1, 2, 2, 4, 2, 3, 4, 2, 0, 0, 1, 2&
507 , 0, 0, 1, 2, 1, 1, 0, 4, 2, 2, 4, 2, 0, 0, 0, 2, 0, 0, 1, 2, 0, 1, 0, 4, 2, 2, 4, 2, 0, 0&
508 , 0, 1, 0, 0, 0, 1, 0, 0, 0, 3, 1, 1, 3, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 3, 1, 1, 3, 0&
509 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0&
510 , 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2&
511 , 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
512 , 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
513 , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
514 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: ma5 = reshape([2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7 &
515 , 5, 6, 7, 8, 0, 2, 3, 4, 2, 2, 4, 5, 3, 4, 4, 6, 4, 5, 6, 6, 0, 1, 2, 3, 1, 2, 3, 4, 2, 3&
516 , 4, 5, 3, 4, 5, 6, 0, 0, 2, 3, 0, 0, 3, 3, 2, 3, 2, 5, 3, 3, 5, 4, 0, 0, 1, 2, 0, 0, 2, 3&
517 , 1, 2, 2, 4, 2, 3, 4, 4, 0, 0, 0, 2, 0, 0, 2, 2, 0, 2, 0, 4, 2, 2, 4, 2, 0, 0, 0, 1, 0, 0&
518 , 1, 1, 0, 1, 0, 3, 1, 1, 3, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 3, 0, 0, 0, 0, 0&
519 , 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 2, 0, 0, 0&
520 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
521 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
522 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
523 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
524 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
525 , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
526 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: mb1 = reshape([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 &
527 , 0, 0, 0, 0, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 0, 2, 3, 2, 2, 4, 2, 2, 3, 2&
528 , 4, 2, 2, 2, 2, 4, 0, 3, 4, 3, 3, 0, 3, 3, 4, 3, 6, 3, 3, 3, 3, 6, 0, 0, 0, 4, 0, 0, 4, 4&
529 , 0, 4, 0, 4, 4, 4, 4, 8, 0, 0, 0, 5, 0, 0, 5, 5, 0, 5, 0, 5, 5, 5, 5, 0, 0, 0, 0, 0, 0, 0&
530 , 0, 6, 0, 0, 0, 6, 0, 6, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 0, 0, 7, 0, 0, 0, 0, 0&
531 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
532 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
533 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
534 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
535 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
536 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
537 , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
538 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: mb2 = reshape([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 &
539 , 1, 1, 1, 1, 2, 0, 2, 2, 2, 2, 0, 2, 2, 2, 2, 0, 2, 2, 2, 2, 0, 3, 0, 0, 0, 0, 3, 0, 0, 0&
540 , 0, 3, 0, 0, 0, 0, 1, 2, 3, 1, 3, 3, 2, 3, 3, 1, 3, 2, 3, 3, 3, 3, 0, 0, 1, 4, 1, 1, 5, 1&
541 , 1, 4, 1, 5, 1, 1, 1, 1, 0, 0, 4, 5, 2, 4, 4, 4, 4, 5, 4, 4, 4, 4, 4, 4, 0, 0, 2, 3, 0, 2&
542 , 0, 2, 2, 3, 2, 7, 2, 2, 2, 2, 0, 0, 3, 4, 0, 3, 0, 5, 3, 4, 5, 6, 5, 5, 5, 5, 0, 0, 0, 0&
543 , 0, 0, 0, 3, 0, 0, 3, 0, 3, 3, 3, 3, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 6, 0, 4, 6, 6, 6, 0, 0&
544 , 0, 0, 0, 0, 0, 4, 0, 0, 4, 0, 0, 4, 4, 4, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 5, 0, 0, 5, 7, 7&
545 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
546 , 6, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
547 , 0, 0, 0, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
548 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
549 , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
550 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: mb3 = reshape([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 &
551 , 1, 1, 1, 1, 2, 2, 2, 2, 0, 2, 2, 2, 2, 0, 2, 2, 2, 2, 0, 2, 0, 0, 0, 0, 3, 0, 0, 0, 0, 3&
552 , 0, 0, 0, 0, 3, 0, 1, 3, 3, 3, 2, 3, 1, 3, 3, 2, 3, 3, 1, 3, 2, 3, 0, 1, 1, 1, 0, 1, 4, 1&
553 , 1, 5, 1, 1, 4, 1, 5, 1, 0, 2, 4, 4, 0, 4, 5, 4, 4, 4, 4, 4, 5, 4, 4, 4, 0, 0, 2, 2, 0, 2&
554 , 3, 2, 2, 0, 2, 2, 3, 2, 7, 2, 0, 0, 3, 5, 0, 3, 4, 5, 3, 0, 5, 5, 4, 5, 6, 5, 0, 0, 0, 3&
555 , 0, 0, 0, 3, 0, 0, 3, 3, 0, 3, 0, 3, 0, 0, 0, 4, 0, 0, 0, 6, 0, 0, 6, 6, 0, 6, 0, 6, 0, 0&
556 , 0, 0, 0, 0, 0, 4, 0, 0, 4, 4, 0, 4, 0, 4, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 5, 7, 0, 5, 0, 7&
557 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 0, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0&
558 , 0, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
559 , 0, 0, 0, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
560 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
561 , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
562 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: mb4 = reshape([2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 &
563 , 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 3, 3, 3, 4, 3, 3, 3, 3&
564 , 4, 3, 3, 3, 3, 4, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 0, 2, 4, 4, 2, 4, 4, 2&
565 , 4, 4, 0, 4, 4, 2, 4, 0, 0, 0, 2, 2, 0, 2, 0, 0, 2, 0, 6, 2, 2, 0, 2, 6, 0, 3, 0, 0, 3, 0&
566 , 5, 5, 0, 5, 4, 0, 0, 5, 0, 2, 0, 1, 3, 5, 1, 0, 1, 1, 3, 1, 2, 5, 5, 1, 5, 8, 0, 0, 1, 3&
567 , 0, 0, 4, 6, 1, 4, 6, 3, 3, 6, 3, 6, 0, 0, 4, 1, 0, 0, 2, 4, 4, 2, 4, 1, 1, 4, 1, 4, 0, 0&
568 , 2, 4, 0, 0, 5, 5, 2, 5, 0, 6, 4, 5, 6, 8, 0, 0, 0, 2, 0, 0, 3, 3, 0, 3, 0, 4, 2, 3, 4, 6&
569 , 0, 0, 0, 5, 0, 0, 0, 6, 0, 0, 0, 2, 5, 6, 2, 0, 0, 0, 0, 3, 0, 0, 0, 4, 0, 0, 0, 7, 3, 4&
570 , 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3&
571 , 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
572 , 0, 4, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 0, 0, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0&
573 , 0, 0, 0, 5, 0, 0, 5, 0], [4, 4, 20])
574 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: mb5 = reshape([2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 &
575 , 2, 2, 2, 2, 0, 3, 3, 3, 3, 4, 3, 3, 3, 3, 4, 3, 3, 3, 3, 4, 0, 0, 4, 4, 0, 0, 4, 0, 4, 4&
576 , 0, 4, 4, 0, 4, 0, 0, 1, 0, 0, 1, 2, 0, 5, 0, 0, 6, 0, 0, 5, 0, 6, 0, 0, 1, 5, 0, 0, 5, 1&
577 , 1, 5, 2, 5, 5, 1, 5, 2, 0, 0, 2, 1, 0, 0, 1, 6, 2, 1, 4, 1, 1, 6, 1, 8, 0, 0, 0, 2, 0, 0&
578 , 2, 3, 0, 2, 0, 6, 2, 3, 6, 4, 0, 0, 0, 3, 0, 0, 3, 4, 0, 3, 0, 2, 3, 4, 2, 6, 0, 0, 0, 0&
579 , 0, 0, 0, 0, 0, 0, 0, 7, 0, 0, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 3, 0, 0, 0&
580 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 5, 0&
581 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
582 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
583 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
584 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
585 , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
586
587 INTEGER :: k, k1, k2, mu
588 REAL(kind=dp) :: cp, ct, fac1, fac2, j, jc, jcc, jss, rr, &
589 sp, st, xx, yy, za, zb
590 REAL(kind=dp), DIMENSION(3) :: v
591 REAL(kind=dp), DIMENSION(3, 3) :: arot
592
593 s(:, :) = 0.0_dp
594
595 v(:) = r(:)
596 rr = norm2(v)
597
598 IF (rr < 1.0e-20_dp) THEN
599
600 DO mu = 1, 4
601 s(mu, mu) = 1.0_dp
602 END DO
603
604 ELSE
605
606 fac1 = 1.0_dp
607 IF (nra == 1) THEN
608 fac1 = fac1*2.0_dp
609 ELSE
610 IF (nra == 2) THEN
611 fac1 = fac1*sqrt(4.0_dp/3.0_dp)
612 ELSE
613 IF (nra == 3) THEN
614 fac1 = fac1*sqrt(8.0_dp/45.0_dp)
615 ELSE
616 IF (nra == 4) THEN
617 fac1 = fac1*sqrt(4.0_dp/315.0_dp)
618 ELSE
619 WRITE (*, *) 'nra= ', nra
620 RETURN
621 END IF
622 END IF
623 END IF
624 END IF
625 IF (nrb == 1) THEN
626 fac1 = fac1*2.0_dp
627 ELSE
628 IF (nrb == 2) THEN
629 fac1 = fac1*sqrt(4.0_dp/3.0_dp)
630 ELSE
631 IF (nrb == 3) THEN
632 fac1 = fac1*sqrt(8.0_dp/45.0_dp)
633 ELSE
634 IF (nrb == 4) THEN
635 fac1 = fac1*sqrt(4.0_dp/315.0_dp)
636 ELSE
637 WRITE (*, *) 'nrb= ', nrb
638 RETURN
639 END IF
640 END IF
641 END IF
642 END IF
643
644 ct = -v(3)/rr
645 IF (abs(ct) < 1.0_dp) THEN
646 st = sqrt(1.0_dp - ct**2)
647 cp = -v(1)/(rr*st)
648 sp = -v(2)/(rr*st)
649 arot(1, 1) = ct*cp
650 arot(1, 2) = -sp
651 arot(1, 3) = st*cp
652 arot(2, 1) = ct*sp
653 arot(2, 2) = cp
654 arot(2, 3) = st*sp
655 arot(3, 1) = -st
656 arot(3, 2) = 0.0_dp
657 arot(3, 3) = ct
658 ELSE
659 arot(1, 1) = ct
660 arot(1, 2) = 0.0_dp
661 arot(1, 3) = 0.0_dp
662 arot(2, 1) = 0.0_dp
663 arot(2, 2) = 1.0_dp
664 arot(2, 3) = 0.0_dp
665 arot(3, 1) = 0.0_dp
666 arot(3, 2) = 0.0_dp
667 arot(3, 3) = ct
668 END IF
669
670 za = zsa
671 zb = zsb
672 fac2 = sqrt(za**(2*nra + 1)*zb**(2*nrb + 1))
673 xx = 0.5_dp*rr*(za + zb)
674 yy = 0.5_dp*rr*(za - zb)
675
676 j = 0.0_dp
677 DO k = 1, nc1(nra, nrb)
678 j = j + real(c1(nra, nrb, k), dp)*aa(ma1(nra, nrb, k), xx)*bb(mb1(nra, nrb, k), yy)
679 END DO
680 j = j*rr**(nra + nrb + 1)
681 j = j/2.0_dp**(nra + nrb + 2)
682
683 s(1, 1) = s(1, 1) + fac1*fac2*j
684
685 za = zpa
686 zb = zsb
687 fac2 = sqrt(za**(2*nra + 1)*zb**(2*nrb + 1))
688 xx = 0.5_dp*rr*(za + zb)
689 yy = 0.5_dp*rr*(za - zb)
690
691 jc = 0.0_dp
692 DO k = 1, nc2(nra, nrb)
693 jc = jc + real(c2(nra, nrb, k), dp)*aa(ma2(nra, nrb, k), xx)*bb(mb2(nra, nrb, k), yy)
694 END DO
695 jc = jc*rr**(nra + nrb + 1)
696 jc = jc/2.0_dp**(nra + nrb + 2)
697
698 DO k1 = 1, 3
699 s(k1 + 1, 1) = s(k1 + 1, 1) &
700 & + sqrt(3.0_dp)*arot(k1, 3)*fac1*fac2*jc
701 END DO
702
703 za = zsa
704 zb = zpb
705 fac2 = sqrt(za**(2*nra + 1)*zb**(2*nrb + 1))
706 xx = 0.5_dp*rr*(za + zb)
707 yy = 0.5_dp*rr*(za - zb)
708
709 jc = 0.0_dp
710 DO k = 1, nc3(nra, nrb)
711 jc = jc + real(c3(nra, nrb, k), dp)*aa(ma3(nra, nrb, k), xx)*bb(mb3(nra, nrb, k), yy)
712 END DO
713 jc = jc*rr**(nra + nrb + 1)
714 jc = jc/2.0_dp**(nra + nrb + 2)
715
716 DO k1 = 1, 3
717 s(1, k1 + 1) = s(1, k1 + 1) &
718 & - sqrt(3.0_dp)*arot(k1, 3)*fac1*fac2*jc
719 END DO
720
721 za = zpa
722 zb = zpb
723 fac2 = sqrt(za**(2*nra + 1)*zb**(2*nrb + 1))
724 xx = 0.5_dp*rr*(za + zb)
725 yy = 0.5_dp*rr*(za - zb)
726
727 jss = 0.0_dp
728 DO k = 1, nc4(nra, nrb)
729 jss = jss + real(c4(nra, nrb, k), dp)*aa(ma4(nra, nrb, k), xx)*bb(mb4(nra, nrb, k), yy)
730 END DO
731 jss = jss*rr**(nra + nrb + 1)
732 jss = jss/2.0_dp**(nra + nrb + 2)
733
734 jcc = 0.0_dp
735 DO k = 1, nc5(nra, nrb)
736 jcc = jcc + real(c5(nra, nrb, k), dp)*aa(ma5(nra, nrb, k), xx)*bb(mb5(nra, nrb, k), yy)
737 END DO
738 jcc = jcc*rr**(nra + nrb + 1)
739 jcc = jcc/2.0_dp**(nra + nrb + 2)
740
741 DO k1 = 1, 3
742 DO k2 = 1, 3
743 s(k1 + 1, k2 + 1) = s(k1 + 1, k2 + 1) &
744 & + 1.5_dp*arot(k1, 1)*arot(k2, 1)*fac1*fac2*jss &
745 & + 1.5_dp*arot(k1, 2)*arot(k2, 2)*fac1*fac2*jss &
746 & - 3.0_dp*arot(k1, 3)*arot(k2, 3)*fac1*fac2*jcc
747 END DO
748 END DO
749
750 END IF
751
752 END SUBROUTINE makes
753
754! **************************************************************************************************
755!> \brief ...
756!> \param R ...
757!> \param nra ...
758!> \param nrb ...
759!> \param ZSA ...
760!> \param ZSB ...
761!> \param ZPA ...
762!> \param ZPB ...
763!> \param dS ...
764! **************************************************************************************************
765 SUBROUTINE makeds(R, nra, nrb, ZSA, ZSB, ZPA, ZPB, dS)
766
767 REAL(kind=dp), DIMENSION(3) :: r
768 INTEGER :: nra, nrb
769 REAL(kind=dp) :: zsa, zsb, zpa, zpb
770 REAL(kind=dp), DIMENSION(4, 4, 3) :: ds
771
772 INTEGER, DIMENSION(4, 4), PARAMETER :: &
773 nc1 = reshape([2, 4, 4, 6, 4, 3, 6, 7, 4, 6, 4, 8, 6, 7, 8, 5], [4, 4]), &
774 nc2 = reshape([4, 4, 8, 8, 6, 8, 6, 12, 8, 8, 12, 8, 10, 12, 14, 16], [4, 4]), &
775 nc3 = reshape([4, 6, 8, 10, 4, 8, 8, 12, 8, 6, 12, 14, 8, 12, 8, 16], [4, 4]), &
776 nc4 = reshape([4, 8, 11, 14, 8, 6, 12, 14, 11, 12, 10, 20, 14, 14, 20, 12], [4, 4]), &
777 nc5 = reshape([2, 4, 6, 8, 4, 4, 8, 8, 6, 8, 6, 12, 8, 8, 12, 8], [4, 4])
778 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: c1 = reshape([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, &
779 1, 1, 1, 1, -1, 1, 2, 3, -1, -2, 1, 2, -2, -1, -3, 1, -3, -2, -1, -4, 0, -1, -2, 2, -1, 1,&
780 -2, -1, 2, -2, 3, -3, 2, -1, -3, 6, 0, -1, -1, -2, 1, 0, -2, -4, -1, 2, -1, -3, 2, 4, 3, &
781 -4, 0, 0, 0, -3, 0, 0, 1, -1, 0, 1, 0, 3, -3, -1, 3, 1, 0, 0, 0, -1, 0, 0, 1, 2, 0, -1, 0,&
782 3, 1, -2, -3, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, -1, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0 &
783 , 0, 0, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
784 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
785 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
786 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
787 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
788 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
789 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
790 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: c2 = reshape([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, &
791 1, 1, 1, 1, -1, 1, 1, 2, -2, -1, 1, 1, -3, -2, -1, 1, -4, -3, -2, -1, 1, -1, 1, 1, 1, 1, &
792 -2, 1, 1, 1, 1, -3, 1, 1, 1, 1, -1, -1, -1, 2, 1, -1, -2, -2, 3, -2, -2, -3, 6, 2, -1, -3,&
793 0, 0, 1, -2, -2, -1, 1, 1, -3, 2, -1, 3, -4, -3, -2, -1, 0, 0, -1, -1, 1, 1, 1, -2, -1, -1&
794 , 2, 3, -4, 2, 4, 3, 0, 0, -1, -2, 0, -1, 0, -2, 3, 2, -2, -1, 6, 2, -1, -3, 0, 0, -1, -1,&
795 0, 1, 0, 1, -1, -1, 1, -1, 1, -3, -1, 3, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 2, 0, -4, 2, 4, 3 &
796 , 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 1, 1, -2, -3, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0&
797 , -3, -1, 3, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 0, 1, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
798 0, 0, 0, 0, 0, -2, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0&
799 , 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, &
800 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
801 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],&
802 [4, 4, 20])
803 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: c3 = reshape([-1, -1, -1, -1, -1, -1, -1, -1, -1, &
804 -1, -1, -1, -1, -1, -1, -1, -1, -2, -3, -4, 1, -1, -2, -3, 1, 1, -1, -2, 2, 1, 1, -1, 1, 1&
805 , 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 3, 1, 1, -1, -3, -6, -1, 1, 2, -2, 1, -2, 2, 1, -2, &
806 2, -3, 3, 0, 2, 3, 4, 0, 1, 2, 3, -1, -1, 1, 2, -2, -1, -3, 1, 0, 1, -1, -4, 0, 1, 1, 2, &
807 -1, 1, 2, 4, 1, -2, 3, 3, 0, 0, 3, 6, 0, -1, -2, 2, -1, 0, -2, -1, 2, -2, 1, -3, 0, 0, 1, &
808 -1, 0, -1, -1, 3, 1, 0, -1, 1, -1, -1, -1, -3, 0, 0, 0, 4, 0, 0, 0, -2, 0, 0, -2, -4, 0, 2&
809 , 0, -3, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, -1, -2, 0, 1, 0, -3, 0, 0, 0, 0, 0, 0, 0, -3, 0, 0,&
810 1, -1, 0, 1, 0, 3, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, -1, 0, -1, 0, 1, 0, 0, 0, 0, 0, 0, 0,&
811 0, 0, 0, 0, 2, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, &
812 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0&
813 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
814 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
815 , 0], [4, 4, 20])
816 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: c4 = reshape([-1, -1, -1, -1, -1, -1, -1, -1, -1, &
817 -1, -1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -2, -3,&
818 1, 1, -1, -2, 2, 1, 2, -1, 3, 2, 1, 3, -1, 1, 2, 3, -1, -1, 1, 2, -2, -1, -1, 1, -3, -2, &
819 -1, -2, 0, 1, -1, -3, 1, -1, 1, 1, -1, 1, -1, 2, -3, 1, 2, -1, 0, -1, 2, 4, -1, 1, -1, -1 &
820 , 2, -1, -1, -1, 4, -1, -1, -3, 0, 1, -1, -1, -1, 0, 1, 2, -1, -1, -1, -1, -1, -2, -1, 3, &
821 0, -1, 2, -1, 1, 0, -1, -2, -2, 1, 2, 2, 1, 2, -2, 1, 0, 0, -2, 4, 0, 0, -1, 1, 2, -1, 1, &
822 -1, -4, 1, 1, 2, 0, 0, 1, -3, 0, 0, 1, -1, 1, 1, -1, -1, 3, -1, 1, -3, 0, 0, -1, 3, 0, 0, &
823 -1, -2, -1, 1, 0, -1, 3, 2, -1, -1, 0, 0, 0, -3, 0, 0, 1, 2, 0, -1, 0, -1, -3, -2, -1, 1, &
824 0, 0, 0, 1, 0, 0, 0, -1, 0, 0, 0, 2, -1, -1, 2, 0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, 0, -1, 1&
825 , 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
826 0, 2, 0, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,&
827 0, 0, 0, -1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 0, &
828 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0], [4, 4, 20])
829 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: c5 = reshape([-1, -1, -1, -1, -1, -1, -1, -1, -1, &
830 -1, -1, -1, -1, -1, -1, -1, 1, -1, -2, -3, 1, 1, -1, -2, 2, 1, 2, -1, 3, 2, 1, 3, 0, 1, -1&
831 , -3, 1, 1, 1, 1, -1, 1, 1, 2, -3, 1, 2, 1, 0, 1, 1, 1, -1, -1, 1, 2, 1, 1, -1, 1, 1, -2, &
832 1, -3, 0, 0, 2, -1, 0, 0, 1, 2, -2, -1, -2, 2, 1, -2, -2, -3, 0, 0, 1, 3, 0, 0, 1, 1, 1, &
833 -1, 1, 1, -3, 1, -1, 1, 0, 0, 0, 3, 0, 0, -1, -2, 0, -1, 0, -1, 3, 2, -1, 3, 0, 0, 0, 1, 0&
834 , 0, -1, -1, 0, 1, 0, -2, -1, -1, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0,&
835 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,&
836 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 &
837 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
838 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
839 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
840 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20&
841 ])
842 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: ma1 = reshape([2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7 &
843 , 5, 6, 7, 8, 0, 2, 3, 4, 2, 2, 4, 5, 3, 4, 4, 6, 4, 5, 6, 6, 0, 1, 1, 3, 1, 0, 3, 4, 1, 3&
844 , 2, 5, 3, 4, 5, 4, 0, 0, 0, 2, 0, 0, 2, 3, 0, 2, 0, 4, 2, 3, 4, 2, 0, 0, 0, 1, 0, 0, 1, 2&
845 , 0, 1, 0, 3, 1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 2, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0&
846 , 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
847 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
848 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
849 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
850 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
851 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
852 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
853 , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
854 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: ma2 = reshape([2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7 &
855 , 5, 6, 7, 8, 1, 2, 3, 4, 2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7, 1, 1, 3, 4, 2, 3, 3, 5, 3, 4&
856 , 5, 5, 4, 5, 6, 7, 0, 0, 2, 3, 1, 2, 2, 4, 2, 3, 4, 4, 3, 4, 5, 6, 0, 0, 2, 2, 1, 2, 1, 4&
857 , 2, 2, 4, 3, 3, 4, 5, 6, 0, 0, 1, 1, 0, 1, 0, 3, 1, 1, 3, 2, 2, 3, 4, 5, 0, 0, 1, 1, 0, 1&
858 , 0, 3, 1, 1, 3, 1, 2, 3, 4, 5, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 2, 0, 1, 2, 3, 4, 0, 0, 0, 0&
859 , 0, 0, 0, 2, 0, 0, 2, 0, 1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 2, 3, 0, 0&
860 , 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2&
861 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
862 , 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
863 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
864 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
865 , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
866 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: ma3 = reshape([2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7 &
867 , 5, 6, 7, 8, 1, 2, 3, 4, 2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7, 1, 2, 3, 4, 1, 3, 4, 5, 3, 3&
868 , 5, 6, 4, 5, 5, 7, 0, 1, 2, 3, 0, 2, 3, 4, 2, 2, 4, 5, 3, 4, 4, 6, 0, 1, 2, 3, 0, 2, 2, 4&
869 , 2, 1, 4, 5, 2, 4, 3, 6, 0, 0, 1, 2, 0, 1, 1, 3, 1, 0, 3, 4, 1, 3, 2, 5, 0, 0, 1, 2, 0, 1&
870 , 1, 3, 1, 0, 3, 4, 1, 3, 1, 5, 0, 0, 0, 1, 0, 0, 0, 2, 0, 0, 2, 3, 0, 2, 0, 4, 0, 0, 0, 1&
871 , 0, 0, 0, 2, 0, 0, 2, 3, 0, 2, 0, 4, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 2, 0, 1, 0, 3, 0, 0&
872 , 0, 0, 0, 0, 0, 1, 0, 0, 1, 2, 0, 1, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 2&
873 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
874 , 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
875 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
876 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
877 , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
878 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: ma4 = reshape([2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7 &
879 , 5, 6, 7, 8, 2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7, 5, 6, 7, 8, 0, 2, 3, 4, 2, 2, 4, 5, 3, 4&
880 , 4, 6, 4, 5, 6, 6, 0, 2, 3, 4, 2, 2, 4, 5, 3, 4, 4, 6, 4, 5, 6, 6, 0, 1, 2, 3, 1, 0, 3, 4&
881 , 2, 3, 4, 5, 3, 4, 5, 6, 0, 1, 2, 3, 1, 0, 3, 4, 2, 3, 2, 5, 3, 4, 5, 4, 0, 0, 2, 3, 0, 0&
882 , 2, 3, 2, 2, 2, 5, 3, 3, 5, 4, 0, 0, 1, 2, 0, 0, 2, 3, 1, 2, 2, 4, 2, 3, 4, 2, 0, 0, 1, 2&
883 , 0, 0, 1, 2, 1, 1, 0, 4, 2, 2, 4, 2, 0, 0, 0, 2, 0, 0, 1, 2, 0, 1, 0, 4, 2, 2, 4, 2, 0, 0&
884 , 0, 1, 0, 0, 0, 1, 0, 0, 0, 3, 1, 1, 3, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 3, 1, 1, 3, 0&
885 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0&
886 , 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2&
887 , 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
888 , 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
889 , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
890 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: ma5 = reshape([2, 3, 4, 5, 3, 4, 5, 6, 4, 5, 6, 7 &
891 , 5, 6, 7, 8, 0, 2, 3, 4, 2, 2, 4, 5, 3, 4, 4, 6, 4, 5, 6, 6, 0, 1, 2, 3, 1, 2, 3, 4, 2, 3&
892 , 4, 5, 3, 4, 5, 6, 0, 0, 2, 3, 0, 0, 3, 3, 2, 3, 2, 5, 3, 3, 5, 4, 0, 0, 1, 2, 0, 0, 2, 3&
893 , 1, 2, 2, 4, 2, 3, 4, 4, 0, 0, 0, 2, 0, 0, 2, 2, 0, 2, 0, 4, 2, 2, 4, 2, 0, 0, 0, 1, 0, 0&
894 , 1, 1, 0, 1, 0, 3, 1, 1, 3, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 3, 0, 0, 0, 0, 0&
895 , 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 2, 0, 0, 0&
896 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
897 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
898 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
899 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
900 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
901 , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
902 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: mb1 = reshape([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 &
903 , 0, 0, 0, 0, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 0, 2, 3, 2, 2, 4, 2, 2, 3, 2&
904 , 4, 2, 2, 2, 2, 4, 0, 3, 4, 3, 3, 0, 3, 3, 4, 3, 6, 3, 3, 3, 3, 6, 0, 0, 0, 4, 0, 0, 4, 4&
905 , 0, 4, 0, 4, 4, 4, 4, 8, 0, 0, 0, 5, 0, 0, 5, 5, 0, 5, 0, 5, 5, 5, 5, 0, 0, 0, 0, 0, 0, 0&
906 , 0, 6, 0, 0, 0, 6, 0, 6, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 0, 0, 7, 0, 0, 0, 0, 0&
907 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
908 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
909 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
910 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
911 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
912 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
913 , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
914 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: mb2 = reshape([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 &
915 , 1, 1, 1, 1, 2, 0, 2, 2, 2, 2, 0, 2, 2, 2, 2, 0, 2, 2, 2, 2, 0, 3, 0, 0, 0, 0, 3, 0, 0, 0&
916 , 0, 3, 0, 0, 0, 0, 1, 2, 3, 1, 3, 3, 2, 3, 3, 1, 3, 2, 3, 3, 3, 3, 0, 0, 1, 4, 1, 1, 5, 1&
917 , 1, 4, 1, 5, 1, 1, 1, 1, 0, 0, 4, 5, 2, 4, 4, 4, 4, 5, 4, 4, 4, 4, 4, 4, 0, 0, 2, 3, 0, 2&
918 , 0, 2, 2, 3, 2, 7, 2, 2, 2, 2, 0, 0, 3, 4, 0, 3, 0, 5, 3, 4, 5, 6, 5, 5, 5, 5, 0, 0, 0, 0&
919 , 0, 0, 0, 3, 0, 0, 3, 0, 3, 3, 3, 3, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 6, 0, 4, 6, 6, 6, 0, 0&
920 , 0, 0, 0, 0, 0, 4, 0, 0, 4, 0, 0, 4, 4, 4, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 5, 0, 0, 5, 7, 7&
921 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
922 , 6, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
923 , 0, 0, 0, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
924 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
925 , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
926 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: mb3 = reshape([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 &
927 , 1, 1, 1, 1, 2, 2, 2, 2, 0, 2, 2, 2, 2, 0, 2, 2, 2, 2, 0, 2, 0, 0, 0, 0, 3, 0, 0, 0, 0, 3&
928 , 0, 0, 0, 0, 3, 0, 1, 3, 3, 3, 2, 3, 1, 3, 3, 2, 3, 3, 1, 3, 2, 3, 0, 1, 1, 1, 0, 1, 4, 1&
929 , 1, 5, 1, 1, 4, 1, 5, 1, 0, 2, 4, 4, 0, 4, 5, 4, 4, 4, 4, 4, 5, 4, 4, 4, 0, 0, 2, 2, 0, 2&
930 , 3, 2, 2, 0, 2, 2, 3, 2, 7, 2, 0, 0, 3, 5, 0, 3, 4, 5, 3, 0, 5, 5, 4, 5, 6, 5, 0, 0, 0, 3&
931 , 0, 0, 0, 3, 0, 0, 3, 3, 0, 3, 0, 3, 0, 0, 0, 4, 0, 0, 0, 6, 0, 0, 6, 6, 0, 6, 0, 6, 0, 0&
932 , 0, 0, 0, 0, 0, 4, 0, 0, 4, 4, 0, 4, 0, 4, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 5, 7, 0, 5, 0, 7&
933 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 0, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0&
934 , 0, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
935 , 0, 0, 0, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
936 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
937 , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
938 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: mb4 = reshape([2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 &
939 , 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 3, 3, 3, 4, 3, 3, 3, 3&
940 , 4, 3, 3, 3, 3, 4, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 0, 2, 4, 4, 2, 4, 4, 2&
941 , 4, 4, 0, 4, 4, 2, 4, 0, 0, 0, 2, 2, 0, 2, 0, 0, 2, 0, 6, 2, 2, 0, 2, 6, 0, 3, 0, 0, 3, 0&
942 , 5, 5, 0, 5, 4, 0, 0, 5, 0, 2, 0, 1, 3, 5, 1, 0, 1, 1, 3, 1, 2, 5, 5, 1, 5, 8, 0, 0, 1, 3&
943 , 0, 0, 4, 6, 1, 4, 6, 3, 3, 6, 3, 6, 0, 0, 4, 1, 0, 0, 2, 4, 4, 2, 4, 1, 1, 4, 1, 4, 0, 0&
944 , 2, 4, 0, 0, 5, 5, 2, 5, 0, 6, 4, 5, 6, 8, 0, 0, 0, 2, 0, 0, 3, 3, 0, 3, 0, 4, 2, 3, 4, 6&
945 , 0, 0, 0, 5, 0, 0, 0, 6, 0, 0, 0, 2, 5, 6, 2, 0, 0, 0, 0, 3, 0, 0, 0, 4, 0, 0, 0, 7, 3, 4&
946 , 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3&
947 , 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
948 , 0, 4, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 0, 0, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0&
949 , 0, 0, 0, 5, 0, 0, 5, 0], [4, 4, 20])
950 INTEGER, DIMENSION(4, 4, 20), PARAMETER :: mb5 = reshape([2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 &
951 , 2, 2, 2, 2, 0, 3, 3, 3, 3, 4, 3, 3, 3, 3, 4, 3, 3, 3, 3, 4, 0, 0, 4, 4, 0, 0, 4, 0, 4, 4&
952 , 0, 4, 4, 0, 4, 0, 0, 1, 0, 0, 1, 2, 0, 5, 0, 0, 6, 0, 0, 5, 0, 6, 0, 0, 1, 5, 0, 0, 5, 1&
953 , 1, 5, 2, 5, 5, 1, 5, 2, 0, 0, 2, 1, 0, 0, 1, 6, 2, 1, 4, 1, 1, 6, 1, 8, 0, 0, 0, 2, 0, 0&
954 , 2, 3, 0, 2, 0, 6, 2, 3, 6, 4, 0, 0, 0, 3, 0, 0, 3, 4, 0, 3, 0, 2, 3, 4, 2, 6, 0, 0, 0, 0&
955 , 0, 0, 0, 0, 0, 0, 0, 7, 0, 0, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 3, 0, 0, 0&
956 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 5, 0&
957 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
958 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
959 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
960 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0&
961 , 0, 0, 0, 0, 0, 0, 0, 0], [4, 4, 20])
962
963 INTEGER :: k, k1, k2, mu
964 REAL(kind=dp) :: cp, ct, dj, djc, djcc, djss, dxx, dyy, &
965 f, fac1, fac2, j, jc, jcc, jss, rr, &
966 sp, st, w, w1, w2, xx, yy, za, zb
967 REAL(kind=dp), DIMENSION(3) :: dcp, dct, dsp, dst, v
968 REAL(kind=dp), DIMENSION(3, 3) :: arot
969 REAL(kind=dp), DIMENSION(3, 3, 3) :: darot
970
971 ds(:, :, :) = 0.0_dp
972
973 v(:) = r(:)
974 rr = norm2(v)
975
976 IF (rr < 1.0e-20_dp) THEN
977
978 DO mu = 1, 4
979 ds(mu, mu, :) = 0.0_dp
980 END DO
981
982 ELSE
983
984 fac1 = 1.0_dp
985 IF (nra == 1) THEN
986 fac1 = fac1*2.0_dp
987 ELSE
988 IF (nra == 2) THEN
989 fac1 = fac1*sqrt(4.0_dp/3.0_dp)
990 ELSE
991 IF (nra == 3) THEN
992 fac1 = fac1*sqrt(8.0_dp/45.0_dp)
993 ELSE
994 IF (nra == 4) THEN
995 fac1 = fac1*sqrt(4.0_dp/315.0_dp)
996 ELSE
997 WRITE (*, *) 'nra= ', nra
998 RETURN
999 END IF
1000 END IF
1001 END IF
1002 END IF
1003 IF (nrb == 1) THEN
1004 fac1 = fac1*2.0_dp
1005 ELSE
1006 IF (nrb == 2) THEN
1007 fac1 = fac1*sqrt(4.0_dp/3.0_dp)
1008 ELSE
1009 IF (nrb == 3) THEN
1010 fac1 = fac1*sqrt(8.0_dp/45.0_dp)
1011 ELSE
1012 IF (nrb == 4) THEN
1013 fac1 = fac1*sqrt(4.0_dp/315.0_dp)
1014 ELSE
1015 WRITE (*, *) 'nrb= ', nrb
1016 RETURN
1017 END IF
1018 END IF
1019 END IF
1020 END IF
1021
1022 ct = -v(3)/rr
1023 IF (abs(ct) >= 1.0_dp) THEN
1024
1025 dct(:) = v(:)*v(3)/rr**3
1026 dct(3) = dct(3) - 1.0_dp/rr
1027
1028 arot(1, 1) = ct
1029 arot(1, 2) = 0.0_dp
1030 arot(1, 3) = 0.0_dp
1031 arot(2, 1) = 0.0_dp
1032 arot(2, 2) = 1.0_dp
1033 arot(2, 3) = 0.0_dp
1034 arot(3, 1) = 0.0_dp
1035 arot(3, 2) = 0.0_dp
1036 arot(3, 3) = ct
1037
1038 darot(1, 1, :) = dct(:)
1039 darot(1, 2, :) = 0.0_dp
1040 darot(1, 3, :) = 0.0_dp
1041 darot(2, 1, :) = 0.0_dp
1042 darot(2, 2, :) = 0.0_dp
1043 darot(2, 3, :) = 0.0_dp
1044 darot(3, 1, :) = 0.0_dp
1045 darot(3, 2, :) = 0.0_dp
1046 darot(3, 3, :) = dct(:)
1047
1048 ELSE
1049
1050 xx = sqrt(v(1)**2 + v(2)**2)
1051 st = xx/rr
1052 cp = -v(1)/xx
1053 sp = -v(2)/xx
1054
1055 dct(:) = v(:)*v(3)/rr**3
1056 dct(3) = dct(3) - 1.0_dp/rr
1057 dst(:) = -ct*dct(:)/st
1058 dcp(:) = v(:)*v(1)/(rr**3*st)
1059 dcp(:) = dcp(:) + v(1)*dst(:)/(rr*st**2)
1060 dcp(1) = dcp(1) - 1.0_dp/(rr*st)
1061 dsp(:) = v(:)*v(2)/(rr**3*st)
1062 dsp(:) = dsp(:) + v(2)*dst(:)/(rr*st**2)
1063 dsp(2) = dsp(2) - 1.0_dp/(rr*st)
1064
1065 arot(1, 1) = ct*cp
1066 arot(1, 2) = -sp
1067 arot(1, 3) = st*cp
1068 arot(2, 1) = ct*sp
1069 arot(2, 2) = cp
1070 arot(2, 3) = st*sp
1071 arot(3, 1) = -st
1072 arot(3, 2) = 0.0_dp
1073 arot(3, 3) = ct
1074
1075 darot(1, 1, :) = dct(:)*cp + ct*dcp(:)
1076 darot(1, 2, :) = -dsp(:)
1077 darot(1, 3, :) = dst(:)*cp + st*dcp(:)
1078 darot(2, 1, :) = dct(:)*sp + ct*dsp(:)
1079 darot(2, 2, :) = dcp(:)
1080 darot(2, 3, :) = dst(:)*sp + st*dsp(:)
1081 darot(3, 1, :) = -dst(:)
1082 darot(3, 2, :) = 0.0_dp
1083 darot(3, 3, :) = dct(:)
1084
1085 END IF
1086
1087 za = zsa
1088 zb = zsb
1089 fac2 = sqrt(za**(2*nra + 1)*zb**(2*nrb + 1))
1090 xx = 0.5_dp*rr*(za + zb)
1091 yy = 0.5_dp*rr*(za - zb)
1092 dxx = 0.5_dp*(za + zb)
1093 dyy = 0.5_dp*(za - zb)
1094
1095 w = 0.0_dp
1096 w1 = 0.0_dp
1097 w2 = 0.0_dp
1098 f = rr**(nra + nrb + 1)/2.0_dp**(nra + nrb + 2)
1099 DO k = 1, nc1(nra, nrb)
1100 w = w + real(c1(nra, nrb, k), dp)*aa(ma1(nra, nrb, k), xx)*bb(mb1(nra, nrb, k), yy)
1101 w1 = w1 + real(c1(nra, nrb, k), dp)*aa(ma1(nra, nrb, k) + 1, xx)*bb(mb1(nra, nrb, k), yy)
1102 w2 = w2 + real(c1(nra, nrb, k), dp)*aa(ma1(nra, nrb, k), xx)*bb(mb1(nra, nrb, k) + 1, yy)
1103 END DO
1104 j = f*w
1105 dj = f*real(nra + nrb + 1, dp)*w/rr
1106 dj = dj - dxx*f*w1
1107 dj = dj - dyy*f*w2
1108
1109 ds(1, 1, :) = ds(1, 1, :) + fac1*fac2*dj*v(:)/rr
1110
1111 za = zpa
1112 zb = zsb
1113 fac2 = sqrt(za**(2*nra + 1)*zb**(2*nrb + 1))
1114 xx = 0.5_dp*rr*(za + zb)
1115 yy = 0.5_dp*rr*(za - zb)
1116 dxx = 0.5_dp*(za + zb)
1117 dyy = 0.5_dp*(za - zb)
1118
1119 w = 0.0_dp
1120 w1 = 0.0_dp
1121 w2 = 0.0_dp
1122 f = rr**(nra + nrb + 1)/2.0_dp**(nra + nrb + 2)
1123 DO k = 1, nc2(nra, nrb)
1124 w = w + real(c2(nra, nrb, k), dp)*aa(ma2(nra, nrb, k), xx)*bb(mb2(nra, nrb, k), yy)
1125 w1 = w1 + real(c2(nra, nrb, k), dp)*aa(ma2(nra, nrb, k) + 1, xx)*bb(mb2(nra, nrb, k), yy)
1126 w2 = w2 + real(c2(nra, nrb, k), dp)*aa(ma2(nra, nrb, k), xx)*bb(mb2(nra, nrb, k) + 1, yy)
1127 END DO
1128 jc = f*w
1129 djc = f*real(nra + nrb + 1, dp)*w/rr
1130 djc = djc - dxx*f*w1
1131 djc = djc - dyy*f*w2
1132
1133 DO k1 = 1, 3
1134 ds(k1 + 1, 1, :) = ds(k1 + 1, 1, :) &
1135 & + sqrt(3.0_dp)*arot(k1, 3)*fac1*fac2*djc*v(:)/rr &
1136 & + sqrt(3.0_dp)*darot(k1, 3, :)*fac1*fac2*jc
1137 END DO
1138
1139 za = zsa
1140 zb = zpb
1141 fac2 = sqrt(za**(2*nra + 1)*zb**(2*nrb + 1))
1142 xx = 0.5_dp*rr*(za + zb)
1143 yy = 0.5_dp*rr*(za - zb)
1144 dxx = 0.5_dp*(za + zb)
1145 dyy = 0.5_dp*(za - zb)
1146
1147 w = 0.0_dp
1148 w1 = 0.0_dp
1149 w2 = 0.0_dp
1150 f = rr**(nra + nrb + 1)/2.0_dp**(nra + nrb + 2)
1151 DO k = 1, nc3(nra, nrb)
1152 w = w + real(c3(nra, nrb, k), dp)*aa(ma3(nra, nrb, k), xx)*bb(mb3(nra, nrb, k), yy)
1153 w1 = w1 + real(c3(nra, nrb, k), dp)*aa(ma3(nra, nrb, k) + 1, xx)*bb(mb3(nra, nrb, k), yy)
1154 w2 = w2 + real(c3(nra, nrb, k), dp)*aa(ma3(nra, nrb, k), xx)*bb(mb3(nra, nrb, k) + 1, yy)
1155 END DO
1156 jc = f*w
1157 djc = f*real(nra + nrb + 1, dp)*w/rr
1158 djc = djc - dxx*f*w1
1159 djc = djc - dyy*f*w2
1160
1161 DO k1 = 1, 3
1162 ds(1, k1 + 1, :) = ds(1, k1 + 1, :) &
1163 & - sqrt(3.0_dp)*arot(k1, 3)*fac1*fac2*djc*v(:)/rr &
1164 & - sqrt(3.0_dp)*darot(k1, 3, :)*fac1*fac2*jc
1165 END DO
1166
1167 za = zpa
1168 zb = zpb
1169 fac2 = sqrt(za**(2*nra + 1)*zb**(2*nrb + 1))
1170 xx = 0.5_dp*rr*(za + zb)
1171 yy = 0.5_dp*rr*(za - zb)
1172 dxx = 0.5_dp*(za + zb)
1173 dyy = 0.5_dp*(za - zb)
1174
1175 w = 0.0_dp
1176 w1 = 0.0_dp
1177 w2 = 0.0_dp
1178 f = rr**(nra + nrb + 1)/2.0_dp**(nra + nrb + 2)
1179 DO k = 1, nc4(nra, nrb)
1180 w = w + real(c4(nra, nrb, k), dp)*aa(ma4(nra, nrb, k), xx)*bb(mb4(nra, nrb, k), yy)
1181 w1 = w1 + real(c4(nra, nrb, k), dp)*aa(ma4(nra, nrb, k) + 1, xx)*bb(mb4(nra, nrb, k), yy)
1182 w2 = w2 + real(c4(nra, nrb, k), dp)*aa(ma4(nra, nrb, k), xx)*bb(mb4(nra, nrb, k) + 1, yy)
1183 END DO
1184 jss = f*w
1185 djss = f*real(nra + nrb + 1, dp)*w/rr
1186 djss = djss - dxx*f*w1
1187 djss = djss - dyy*f*w2
1188
1189 w = 0.0_dp
1190 w1 = 0.0_dp
1191 w2 = 0.0_dp
1192 f = rr**(nra + nrb + 1)/2.0_dp**(nra + nrb + 2)
1193 DO k = 1, nc5(nra, nrb)
1194 w = w + real(c5(nra, nrb, k), dp)*aa(ma5(nra, nrb, k), xx)*bb(mb5(nra, nrb, k), yy)
1195 w1 = w1 + real(c5(nra, nrb, k), dp)*aa(ma5(nra, nrb, k) + 1, xx)*bb(mb5(nra, nrb, k), yy)
1196 w2 = w2 + real(c5(nra, nrb, k), dp)*aa(ma5(nra, nrb, k), xx)*bb(mb5(nra, nrb, k) + 1, yy)
1197 END DO
1198 jcc = f*w
1199 djcc = f*real(nra + nrb + 1, dp)*w/rr
1200 djcc = djcc - dxx*f*w1
1201 djcc = djcc - dyy*f*w2
1202
1203 DO k1 = 1, 3
1204 DO k2 = 1, 3
1205 ds(k1 + 1, k2 + 1, :) = ds(k1 + 1, k2 + 1, :) &
1206 & + 1.5_dp*arot(k1, 1)*arot(k2, 1)*fac1*fac2*djss*v(:)/rr &
1207 & + 1.5_dp*darot(k1, 1, :)*arot(k2, 1)*fac1*fac2*jss &
1208 & + 1.5_dp*arot(k1, 1)*darot(k2, 1, :)*fac1*fac2*jss &
1209 & + 1.5_dp*arot(k1, 2)*arot(k2, 2)*fac1*fac2*djss*v(:)/rr &
1210 & + 1.5_dp*darot(k1, 2, :)*arot(k2, 2)*fac1*fac2*jss &
1211 & + 1.5_dp*arot(k1, 2)*darot(k2, 2, :)*fac1*fac2*jss &
1212 & - 3.0_dp*arot(k1, 3)*arot(k2, 3)*fac1*fac2*djcc*v(:)/rr &
1213 & - 3.0_dp*darot(k1, 3, :)*arot(k2, 3)*fac1*fac2*jcc &
1214 & - 3.0_dp*arot(k1, 3)*darot(k2, 3, :)*fac1*fac2*jcc
1215 END DO
1216 END DO
1217
1218 END IF
1219
1220 END SUBROUTINE makeds
1221
1222! **************************************************************************************************
1223!> \brief ...
1224!> \param n ...
1225!> \param x ...
1226!> \return ...
1227! **************************************************************************************************
1228 FUNCTION aa(n, x)
1229
1230 INTEGER :: n
1231 REAL(kind=dp) :: x, aa
1232
1233 REAL(kind=dp) :: p
1234
1235 IF (n == 0) THEN
1236 p = 1.0_dp
1237 ELSE
1238 IF (n == 1) THEN
1239 p = 1.0_dp + x
1240 ELSE
1241 IF (n == 2) THEN
1242 p = 2.0_dp + x*( &
1243 2.0_dp + x)
1244 ELSE
1245 IF (n == 3) THEN
1246 p = 6.0_dp + x*( &
1247 6.0_dp + x*( &
1248 3.0_dp + x))
1249 ELSE
1250 IF (n == 4) THEN
1251 p = 24.0_dp + x*( &
1252 24.0_dp + x*( &
1253 12.0_dp + x*( &
1254 4.0_dp + x)))
1255 ELSE
1256 IF (n == 5) THEN
1257 p = 120.0_dp + x*( &
1258 120.0_dp + x*( &
1259 60.0_dp + x*( &
1260 20.0_dp + x*( &
1261 5.0_dp + x))))
1262 ELSE
1263 IF (n == 6) THEN
1264 p = 720.0_dp + x*( &
1265 720.0_dp + x*( &
1266 360.0_dp + x*( &
1267 120.0_dp + x*( &
1268 30.0_dp + x*( &
1269 6.0_dp + x)))))
1270 ELSE
1271 IF (n == 7) THEN
1272 p = 5040.0_dp + x*( &
1273 5040.0_dp + x*( &
1274 2520.0_dp + x*( &
1275 840.0_dp + x*( &
1276 210.0_dp + x*( &
1277 42.0_dp + x*( &
1278 7.0_dp + x))))))
1279 ELSE
1280 IF (n == 8) THEN
1281 p = 40320.0_dp + x*( &
1282 40320.0_dp + x*( &
1283 20160.0_dp + x*( &
1284 6720.0_dp + x*( &
1285 1680.0_dp + x*( &
1286 336.0_dp + x*( &
1287 56.0_dp + x*( &
1288 8.0_dp + x)))))))
1289 ELSE
1290 IF (n == 9) THEN
1291 p = 362880.0_dp + x*( &
1292 362880.0_dp + x*( &
1293 181440.0_dp + x*( &
1294 60480.0_dp + x*( &
1295 15120.0_dp + x*( &
1296 3024.0_dp + x*( &
1297 504.0_dp + x*( &
1298 72.0_dp + x*( &
1299 9.0_dp + x))))))))
1300 ELSE
1301 IF (n == 10) THEN
1302 p = 3628800.0_dp + x*( &
1303 3628800.0_dp + x*( &
1304 1814400.0_dp + x*( &
1305 604800.0_dp + x*( &
1306 151200.0_dp + x*( &
1307 30240.0_dp + x*( &
1308 5040.0_dp + x*( &
1309 720.0_dp + x*( &
1310 90.0_dp + x*( &
1311 10.0_dp + x)))))))))
1312 ELSE
1313 p = 1.0_dp
1314 WRITE (*, *) ' n= ', n, ' in AA(n,x) '
1315 END IF
1316 END IF
1317 END IF
1318 END IF
1319 END IF
1320 END IF
1321 END IF
1322 END IF
1323 END IF
1324 END IF
1325 END IF
1326
1327 aa = exp(-x)*p/x**(n + 1)
1328
1329 END FUNCTION aa
1330
1331! **************************************************************************************************
1332!> \brief ...
1333!> \param n ...
1334!> \param y ...
1335!> \return ...
1336! **************************************************************************************************
1337 FUNCTION bb(n, y)
1338
1339 INTEGER :: n
1340 REAL(kind=dp) :: y, bb
1341
1342 IF (abs(y) > 1.0e-20_dp) THEN
1343 bb = real((-1)**(n + 1), dp)*aa(n, -y) - aa(n, y)
1344 ELSE
1345 IF (mod(n, 2) == 0) THEN
1346 bb = 2.0_dp/real(n + 1, dp)
1347 ELSE
1348 bb = -y*2.0_dp/real(n + 2, dp)
1349 END IF
1350 END IF
1351
1352 END FUNCTION bb
1353
1354END MODULE se_core_matrix
1355
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
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_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_deallocate_matrix
subroutine, public dbcsr_deallocate_matrix(matrix)
...
Definition cp_dbcsr_api.F:235
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:693
cp_dbcsr_api::dbcsr_replicate_all
subroutine, public dbcsr_replicate_all(matrix)
...
Definition cp_dbcsr_api.F:1156
cp_dbcsr_api::dbcsr_distribute
subroutine, public dbcsr_distribute(matrix)
...
Definition cp_dbcsr_api.F:540
cp_dbcsr_api::dbcsr_sum_replicated
subroutine, public dbcsr_sum_replicated(matrix)
...
Definition cp_dbcsr_api.F:1223
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:254
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_get_block_diag
subroutine, public dbcsr_get_block_diag(matrix, diag)
Copies the diagonal blocks of matrix into diag.
Definition cp_dbcsr_contrib.F:450
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, cartesian_basis)
...
Definition cp_dbcsr_output.F:165
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:853
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:1073
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:18
input_constants::do_method_pdg
integer, parameter, public do_method_pdg
Definition input_constants.F:263
input_constants::do_method_pnnl
integer, parameter, public do_method_pnnl
Definition input_constants.F:263
input_constants::do_method_rm1
integer, parameter, public do_method_rm1
Definition input_constants.F:263
input_constants::do_method_pm3
integer, parameter, public do_method_pm3
Definition input_constants.F:263
input_constants::do_method_mndo
integer, parameter, public do_method_mndo
Definition input_constants.F:263
input_constants::do_method_mndod
integer, parameter, public do_method_mndod
Definition input_constants.F:263
input_constants::do_method_am1
integer, parameter, public do_method_am1
Definition input_constants.F:263
input_constants::do_method_pm6fm
integer, parameter, public do_method_pm6fm
Definition input_constants.F:263
input_constants::do_method_pm6
integer, parameter, public do_method_pm6
Definition input_constants.F:263
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:1052
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, mimic, 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_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, 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, matrix_vhxc, 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, xcint_weights, 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, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_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, harris_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, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:534
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, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, 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_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, hund_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, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:491
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, exc_accint, 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_vhxc, 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, xcint_weights, 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_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
Definition qs_ks_types.F:591
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:550
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, ext_kpoints)
Calculation of the overlap matrix over Cartesian Gaussian functions.
Definition qs_overlap.F:121
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:364
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:337
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:142
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:757
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
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:750
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:222
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:190
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:137
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