Octopus
test.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
20#include "global.h"
21
22module test_oct_m
23 use accel_oct_m
24 use batch_oct_m
27 use box_oct_m
32 use clock_oct_m
34 use debug_oct_m
41 use global_oct_m
42 use grid_oct_m
48 use iihash_oct_m
50 use iso_c_binding
51 use io_csv_oct_m
52 use io_oct_m
54 use, intrinsic :: iso_fortran_env
59 use mesh_oct_m
67 use mpi_oct_m
74 use parser_oct_m
80 use sihash_oct_m
82 use space_oct_m
83 use sphash_oct_m
92 use types_oct_m
93 use v_ks_oct_m
97
98 ! TODO(Alex) Issue 951. Migrate the test implementations to separate modules in $OCTOPUS_ROOT/test_integrations/
99 ! Test implementation modules
102
103 implicit none
104
105 !Auxiliary quantities needed by the linear solver
106 real(real64) :: shift_aux
107 type(derivatives_t), pointer :: der_aux => null()
108 type(preconditioner_t) :: prec_aux
109
110 public :: test_run
111
112contains
113
117 subroutine test_run(namespace)
118 type(namespace_t), intent(in) :: namespace
119
120 type(test_parameters_t) :: param
121 integer :: test_mode
122
123 push_sub(test_run)
124
125 call messages_obsolete_variable(namespace, 'WhichTest', 'TestMode')
126
127 call param%init_from_file(namespace)
128
129 !%Variable TestMode
130 !%Type integer
131 !%Default hartree
132 !%Section Calculation Modes::Test
133 !%Description
134 !% Decides what kind of test should be performed.
135 !%Option hartree 1
136 !% Tests the Poisson solvers used to calculate the Hartree potential.
137 !%Option derivatives 2
138 !% Tests and benchmarks the implementation of the finite-difference operators, used to calculate derivatives.
139 !%Option orthogonalization 3
140 !% Tests the implementation of the orthogonalization routines.
141 !%Option interpolation 4
142 !% Test the interpolation routines.
143 !%Option ion_interaction 5
144 !% Tests the ion-ion interaction routines.
145 !%Option projector 6
146 !% Tests the code that applies the nonlocal part of the pseudopotentials
147 !% in case of spin-orbit coupling
148 !%Option dft_u 7
149 !% Tests the DFT+U part of the code for projections on the basis.
150 !%Option hamiltonian_apply 8
151 !% Tests the application of the Hamiltonian, or a part of it
152 !%Option density_calc 9
153 !% Calculation of the density.
154 !%Option exp_apply 10
155 !% Tests the exponential of the Hamiltonian
156 !%Option boundaries 11
157 !% Tests the boundaries conditions
158 !%Option subspace_diag 12
159 !% Tests the subspace diagonalization
160 !%Option batch_ops 13
161 !% Tests the batch operations
162 !%Option clock 18
163 !% Tests for clock
164 !%Option linear_solver 19
165 !% Tests the linear solvers
166 !%Option cgal 20
167 !% Tests for cgal interface
168 !%Option dense_eigensolver 21
169 !% Tests for dense eigensolvers (especially parallel ones)
170 !%Option grid_interpolation 22
171 !% Tests for grid interpolation and multigrid methods.
172 !%Option iihash 23
173 !% Tests for the integer-integer hash table.
174 !%Option sihash 24
175 !% Tests for the string-integer hash table.
176 !%Option sphash 25
177 !% Tests for the string-polymorphic hash table.
178 !%Option mpiwrappers 26
179 !% Tests for the MPI wrappers with large integer displacements.
180 !%Option regridding 27
181 !% Tests the regridding between two different grids.
182 !%Option helmholtz_decomposition 28
183 !% Test for the Helmholtz decomposition subroutines
184 !%Option vecpot_analytical 29
185 !% Tests analytically the vector potential from B field.
186 !%Option current_density 30
187 !% Tests the different contributions to the total electronic current density
188 !%Option mixing_tests 31
189 !% Unit test for the mixing
190 !%Option optimizers 32
191 !% Tests for the optimizers, using standard test functions
192 !%Option weighted_kmeans 33
193 !% Tests weighted kmeans algorithm, for assignment of centroids at clustered
194 !% around regions of high weight.
195 !%Option csv_input 34
196 !% Tests the reading of csv files in various dimensions
197 !%Option composition_chebyshev 35
198 !% Tests the chebyshev filtering and the related composition rule
199 !%Option lalg_adv 36
200 !% Tests for the advanced linear algebra routines
201 !%Option isdf_serial 37
202 !% Tests construction of interoperable separable density fit vectors using the
203 !% serial reference implementation.
204 !%Option isdf 38
205 !% Tests construction of interoperable separable density fit vectors.
206 !%Option centroids_t 39
207 !% Test centroids type methods
208 !%End
209
210 call parse_variable(namespace, 'TestMode', option__testmode__hartree, test_mode)
211
212 call messages_print_with_emphasis(msg="Test mode", namespace=namespace)
213 call messages_print_var_option("TestMode", test_mode, namespace=namespace)
214 call messages_print_var_option("TestType", param%type, namespace=namespace)
215 call messages_print_var_value("TestRepetitions", param%repetitions, namespace=namespace)
216 call messages_print_var_value("TestMinBlockSize", param%min_blocksize, namespace=namespace)
217 call messages_print_var_value("TestMaxBlockSize", param%max_blocksize, namespace=namespace)
218 call messages_print_with_emphasis(namespace=namespace)
219
220 select case (test_mode)
221 case (option__testmode__hartree)
222 call test_hartree(param, namespace)
223 case (option__testmode__derivatives)
224 call test_derivatives(param, namespace)
225 case (option__testmode__orthogonalization)
226 call test_orthogonalization(param, namespace)
227 case (option__testmode__interpolation)
228 call test_interpolation(param, namespace)
229 case (option__testmode__ion_interaction)
230 call test_ion_interaction(namespace)
231 case (option__testmode__projector)
232 call test_projector(param, namespace)
233 case (option__testmode__dft_u)
234 call test_dft_u(param, namespace)
235 case (option__testmode__hamiltonian_apply)
236 call test_hamiltonian(param, namespace)
237 case (option__testmode__density_calc)
238 call test_density_calc(param, namespace)
239 case (option__testmode__exp_apply)
240 call test_exponential(param, namespace)
241 case (option__testmode__boundaries)
242 call test_boundaries(param, namespace)
243 case (option__testmode__subspace_diag)
244 call test_subspace_diagonalization(param, namespace)
245 case (option__testmode__batch_ops)
246 call test_batch_ops(param, namespace)
247 case (option__testmode__clock)
248 call test_clock()
249 case (option__testmode__linear_solver)
250 call test_linear_solver(namespace)
251 case (option__testmode__cgal)
252 call test_cgal()
253 case (option__testmode__dense_eigensolver)
255 case (option__testmode__grid_interpolation)
257 case (option__testmode__iihash)
258 call test_iihash()
259 case (option__testmode__sihash)
260 call test_sihash()
261 case (option__testmode__sphash)
262 call test_sphash(namespace)
263 case (option__testmode__mpiwrappers)
264 call test_mpiwrappers()
265 case (option__testmode__regridding)
266 call test_regridding(namespace)
267 case (option__testmode__helmholtz_decomposition)
268 call test_helmholtz_decomposition(namespace)
269 case (option__testmode__vecpot_analytical)
270 call test_vecpot_analytical(namespace)
271 case (option__testmode__current_density)
272 call test_current_density(namespace)
273 case (option__testmode__mixing_tests)
274 call mix_tests_run()
275 case (option__testmode__optimizers)
276 call test_optimizers(namespace)
277 case (option__testmode__weighted_kmeans)
278 call test_weighted_kmeans(namespace)
279 case (option__testmode__csv_input)
280 call test_csv_input(namespace)
281 case (option__testmode__composition_chebyshev)
282 call test_composition_chebyshev(namespace)
283 case (option__testmode__lalg_adv)
284 call test_exponential_matrix(namespace)
285 case (option__testmode__isdf_serial)
286 call test_isdf(namespace, serial=.true.)
287 case (option__testmode__isdf)
288 call test_isdf(namespace)
289 end select
290
291 pop_sub(test_run)
292 end subroutine test_run
293
294 ! ---------------------------------------------------------
295 subroutine test_hartree(param, namespace)
296 type(test_parameters_t), intent(in) :: param
297 type(namespace_t), intent(in) :: namespace
298
299 type(electrons_t), pointer :: sys
300
301 push_sub(test_hartree)
302
303 call calc_mode_par%set_parallelization(p_strategy_states, default = .false.)
304
305 sys => electrons_t(namespace, generate_epot=.false.)
306 call sys%init_parallelization(mpi_world)
307 call poisson_test(sys%hm%psolver, sys%space, sys%gr, sys%ions%latt, namespace, param%repetitions)
308 safe_deallocate_p(sys)
309
310 pop_sub(test_hartree)
311 end subroutine test_hartree
312
313 ! ---------------------------------------------------------
314 subroutine test_helmholtz_decomposition(namespace)
315 type(namespace_t), intent(in) :: namespace
316
317 type(electrons_t), pointer :: sys
318 type(helmholtz_decomposition_t) :: helmholtz
319
321
322 ! First of all we have to initialize the grid and the poisson solver
323 call calc_mode_par%set_parallelization(p_strategy_states, default = .false.)
324
325 sys => electrons_t(namespace, generate_epot=.false.)
326 call sys%init_parallelization(mpi_world)
327
328 call helmholtz%init(namespace, sys%gr, sys%mc, sys%space)
329
330 ! Then we have to initialize the exact fields
331 write(message(1),'(a)') "Helmholtz decomposition: beginning Hertzian dipole test"
332 call messages_info(1, namespace=namespace)
333 call hertzian_dipole_test(helmholtz, sys%gr, sys%namespace, sys%space)
334
335 write(message(1),'(a)') "Helmholtz decomposition: beginning Gaussian test"
336 call messages_info(1, namespace=namespace)
337 call gaussian_test(helmholtz, sys%gr, sys%namespace, sys%space)
338
339 safe_deallocate_p(sys)
340
342 end subroutine test_helmholtz_decomposition
343
344 ! ---------------------------------------------------------
345 subroutine test_linear_solver(namespace)
346 use, intrinsic :: iso_c_binding, only: c_ptr
347 type(namespace_t), intent(in) :: namespace
348
349 type(electrons_t), pointer :: sys
350 real(real64), allocatable :: rho(:), x(:), center(:)
351 real(real64) :: rr, alpha, beta, res
352 integer :: ip, iter
353
354 real(real64), parameter :: threshold = 1.e-7_real64
355 type(c_ptr) :: userdata(0) ! workaround, literal [c_ptr::] crashes ifort
356
357 push_sub(test_linear_solver)
358
359 call calc_mode_par%set_parallelization(p_strategy_states, default = .false.)
360
361 sys => electrons_t(namespace, generate_epot=.false.)
362 call sys%init_parallelization(mpi_world)
363
364 ! We need to set up some auxiliary quantities called by the linear solver
365 call mesh_init_mesh_aux(sys%gr)
366 ! Shift of the linear equation
367 shift_aux = 0.25_real64
368 ! Preconditioner used for the QMR algorithm
369 call preconditioner_init(prec_aux, namespace, sys%gr, sys%mc, sys%space)
370 ! Derivative object needed
371 call set_der_aux(sys%gr%der)
372
373 ! Here we put a Gaussian as the right-hand side of the linear solver
374 ! Values are taken from the poisson_test routine
375 alpha = m_four * sys%gr%spacing(1)
376 beta = m_one / (alpha**sys%space%dim * sqrt(m_pi)**sys%space%dim)
377 ! The Gaussian is centered around the origin
378 safe_allocate(center(1:sys%space%dim))
379 center = m_zero
380
381 safe_allocate(rho(1:sys%gr%np))
382 rho = m_zero
383 do ip = 1, sys%gr%np
384 call mesh_r(sys%gr, ip, rr, origin = center(:))
385 rho(ip) = beta*exp(-(rr/alpha)**2)
386 end do
387
388 safe_allocate(x(1:sys%gr%np))
389
390 !Test the CG linear solver
391 x = m_zero
392 iter = 1000
393 call dconjugate_gradients(sys%gr%np, x, rho, laplacian_op, dmf_dotp_aux, iter, res, threshold, userdata)
394 write(message(1),'(a,i6,a)') "Info: CG converged with ", iter, " iterations."
395 write(message(2),'(a,e14.6)') "Info: The residue is ", res
396 write(message(3),'(a,e14.6)') "Info: Norm solution CG ", dmf_nrm2(sys%gr, x)
397 call messages_info(3, namespace=namespace)
398
399 call preconditioner_end(prec_aux)
400 safe_deallocate_a(x)
401 safe_deallocate_a(rho)
402 safe_deallocate_p(sys)
403
404 pop_sub(test_linear_solver)
405 end subroutine test_linear_solver
406
407 subroutine set_der_aux(der)
408 type(derivatives_t), target, intent(in) :: der
410 der_aux => der
412 end subroutine set_der_aux
413
414 ! ---------------------------------------------------------
416 subroutine laplacian_op(x, hx, userdata)
417 use, intrinsic :: iso_c_binding, only: c_ptr
418 real(real64), contiguous, intent(in) :: x(:)
419 real(real64), contiguous, intent(out) :: Hx(:)
420 type(c_ptr), intent(in) :: userdata(:)
421
422 real(real64), allocatable :: tmpx(:)
423
424 assert(associated(mesh_aux))
425
426 safe_allocate(tmpx(1:mesh_aux%np_part))
427 call lalg_copy(mesh_aux%np, x, tmpx)
428 call dderivatives_lapl(der_aux, tmpx, hx)
429 call lalg_scal(mesh_aux%np, -m_one, hx)
430 call lalg_axpy(mesh_aux%np, shift_aux, x, hx)
431 safe_deallocate_a(tmpx)
432
433 end subroutine laplacian_op
434
435 ! ---------------------------------------------------------
436 subroutine test_projector(param, namespace)
437 type(test_parameters_t), intent(in) :: param
438 type(namespace_t), intent(in) :: namespace
439
440 type(electrons_t), pointer :: sys
441 type(wfs_elec_t) :: epsib
442 integer :: itime
443 complex(real64), allocatable :: psi(:, :)
444
445 push_sub(test_projector)
446
447 call calc_mode_par%set_parallelization(p_strategy_states, default = .false.)
448
449 call messages_write('Info: Testing the nonlocal part of the pseudopotential with SOC')
450 call messages_new_line()
451 call messages_new_line()
452 call messages_info(namespace=namespace)
453
454 sys => electrons_t(namespace, generate_epot=.false.)
455 call sys%init_parallelization(mpi_world)
456
457 call states_elec_allocate_wfns(sys%st, sys%gr, wfs_type = type_cmplx)
458 call test_batch_set_gaussian(sys%st%group%psib(1, 1), sys%gr)
459
460 ! Initialize external potential
461 call hamiltonian_elec_epot_generate(sys%hm, sys%namespace, sys%space, sys%gr, sys%ions, sys%ext_partners, sys%st)
462
463 call sys%st%group%psib(1, 1)%copy_to(epsib)
464
465 call batch_set_zero(epsib)
466
467 do itime = 1, param%repetitions
468 call zproject_psi_batch(sys%gr, sys%gr%der%boundaries, sys%hm%ep%proj, &
469 sys%hm%ep%natoms, 2, sys%st%group%psib(1, 1), epsib)
470 end do
471
472 safe_allocate(psi(1:sys%gr%np, 1:sys%st%d%dim))
473 do itime = 1, epsib%nst
474 call batch_get_state(epsib, itime, sys%gr%np, psi)
475 write(message(1),'(a,i1,3x, f12.6)') "Norm state ", itime, zmf_nrm2(sys%gr, 2, psi)
476 call messages_info(1, namespace=sys%namespace)
477 end do
478 safe_deallocate_a(psi)
479
480 call epsib%end()
481 call states_elec_deallocate_wfns(sys%st)
482 safe_deallocate_p(sys)
483
484 pop_sub(test_projector)
485 end subroutine test_projector
486
487 ! ---------------------------------------------------------
488 subroutine test_dft_u(param, namespace)
489 type(test_parameters_t), intent(in) :: param
490 type(namespace_t), intent(in) :: namespace
491
492 type(electrons_t), pointer :: sys
493 type(wfs_elec_t) :: epsib, epsib2
494 integer :: itime, ist
495 type(orbitalbasis_t) :: basis
496 real(real64), allocatable :: ddot(:,:,:), dweight(:,:)
497 complex(real64), allocatable :: zdot(:,:,:), zweight(:,:)
498
499 push_sub(test_dft_u)
500
501 call calc_mode_par%unset_parallelization(p_strategy_states)
502 call calc_mode_par%unset_parallelization(p_strategy_kpoints)
503 call calc_mode_par%set_parallelization(p_strategy_domains, default = .true.)
504
505 call messages_write('Info: Testing some DFT+U routines')
506 call messages_new_line()
507 call messages_new_line()
508 call messages_info(namespace=namespace)
509
510 sys => electrons_t(namespace, generate_epot=.false.)
511 call sys%init_parallelization(mpi_world)
512
513 call states_elec_allocate_wfns(sys%st, sys%gr)
514 call test_batch_set_gaussian(sys%st%group%psib(1, 1), sys%gr)
515 if (sys%st%pack_states) call sys%st%pack()
516
517 call sys%st%group%psib(1, 1)%copy_to(epsib2, copy_data = .true.)
518
519 ! We set the phase of the batch if needed
520 if (.not. sys%hm%phase%is_allocated()) then
521 call sys%st%group%psib(1, 1)%copy_to(epsib, copy_data = .true.)
522 else
523 call sys%st%group%psib(1, 1)%copy_to(epsib)
524 call sys%hm%phase%apply_to(sys%gr, sys%gr%np, &
525 .false., epsib, src=sys%st%group%psib(1, 1))
526 epsib2%has_phase = .true.
527 end if
528
529 ! Initialize the orbital basis
530 call orbitalbasis_init(basis, sys%namespace, sys%space%periodic_dim)
531 if (states_are_real(sys%st)) then
532 call dorbitalbasis_build(basis, sys%namespace, sys%ions, sys%gr, sys%st%d%kpt, sys%st%d%dim, &
533 sys%hm%phase%is_allocated(), .false., .false.)
534 safe_allocate(dweight(1:basis%orbsets(1)%norbs, 1:epsib%nst_linear))
535 safe_allocate(ddot(1:sys%st%d%dim, 1:basis%orbsets(1)%norbs, 1:epsib%nst))
536 else
537 call zorbitalbasis_build(basis, sys%namespace, sys%ions, sys%gr, sys%st%d%kpt, sys%st%d%dim, &
538 sys%hm%phase%is_allocated(), .false., .false.)
539 call orbitalset_update_phase(basis%orbsets(1), sys%space%dim, sys%st%d%kpt, sys%kpoints, (sys%st%d%ispin==spin_polarized))
540 safe_allocate(zweight(1:basis%orbsets(1)%norbs, 1:epsib%nst_linear))
541 safe_allocate(zdot(1:sys%st%d%dim, 1:basis%orbsets(1)%norbs, 1:epsib%nst))
542
543 !We set the phase of the orbitals if needed
544 if (sys%hm%phase%is_allocated()) then
545 call orbitalset_update_phase(basis%orbsets(1), sys%space%dim, sys%st%d%kpt, sys%kpoints, &
546 (sys%st%d%ispin==spin_polarized))
547 end if
548 end if
549
550 do itime = 1, param%repetitions
551 call batch_set_zero(epsib2)
552 if (states_are_real(sys%st)) then
553 dweight = m_one
554 ddot = m_zero
555 call dorbitalset_get_coeff_batch(basis%orbsets(1), sys%st%d%dim, sys%st%group%psib(1, 1), ddot)
556 call dorbitalset_add_to_batch(basis%orbsets(1), sys%st%d%dim, epsib2, dweight)
557 else
558 zweight = m_one
559 zdot = m_zero
560 call zorbitalset_get_coeff_batch(basis%orbsets(1), sys%st%d%dim, epsib, zdot)
561 call zorbitalset_add_to_batch(basis%orbsets(1), sys%st%d%dim, epsib2, zweight)
562 end if
563 end do
564
565 if (epsib%is_packed()) then
566 call epsib%do_unpack(force = .true.)
567 end if
568
569 do ist = 1, epsib%nst
570 if (states_are_real(sys%st)) then
571 write(message(1),'(a,i2,3x,e13.6)') "Dotp state ", ist, ddot(1,1,ist)
572 else
573 write(message(1),'(a,i2,2(3x,e13.6))') "Dotp state ", ist, zdot(1,1,ist)
574 end if
575 call messages_info(1)
576 end do
577
578
579 call test_prints_info_batch(sys%st, sys%gr, epsib2)
580
581 safe_deallocate_a(dweight)
582 safe_deallocate_a(zweight)
583 safe_deallocate_a(ddot)
584 safe_deallocate_a(zdot)
585
586 call epsib%end()
587 call epsib2%end()
588 call orbitalbasis_end(basis)
589 call states_elec_deallocate_wfns(sys%st)
590 safe_deallocate_p(sys)
591
592 pop_sub(test_dft_u)
593 end subroutine test_dft_u
594
595 ! ---------------------------------------------------------
596 subroutine test_hamiltonian(param, namespace)
597 type(test_parameters_t), intent(in) :: param
598 type(namespace_t), intent(in) :: namespace
599
600 type(electrons_t), pointer :: sys
601 type(wfs_elec_t) :: hpsib
602 integer :: itime, terms
603
604 push_sub(test_hamiltonian)
605
606 !%Variable TestHamiltonianApply
607 !%Type integer
608 !%Default term_all
609 !%Section Calculation Modes::Test
610 !%Description
611 !% Decides which part of the Hamiltonian is applied.
612 !%Option term_all 0
613 !% Apply the full Hamiltonian.
614 !%Option term_kinetic 1
615 !% Apply only the kinetic operator
616 !%Option term_local_potential 2
617 !% Apply only the local potential.
618 !%Option term_non_local_potential 4
619 !% Apply only the non_local potential.
620 !%End
621 call parse_variable(namespace, 'TestHamiltonianApply', option__testhamiltonianapply__term_all, terms)
622 if (terms == 0) terms = huge(1)
623
624
625 call calc_mode_par%set_parallelization(p_strategy_states, default = .false.)
626
627 call messages_write('Info: Testing the application of the Hamiltonian')
628 call messages_new_line()
629 call messages_new_line()
630 call messages_info(namespace=namespace)
631
632 sys => electrons_t(namespace, generate_epot=.false.)
633 call sys%init_parallelization(mpi_world)
634
635 call states_elec_allocate_wfns(sys%st, sys%gr)
636 call test_batch_set_gaussian(sys%st%group%psib(1, 1), sys%gr)
637
638 ! Initialize external potential
639 if (sys%st%pack_states .and. sys%hm%apply_packed()) call sys%st%pack()
640 call hamiltonian_elec_epot_generate(sys%hm, sys%namespace, sys%space, sys%gr, sys%ions, sys%ext_partners, sys%st)
641 call density_calc(sys%st, sys%gr, sys%st%rho)
642 call v_ks_calc(sys%ks, sys%namespace, sys%space, sys%hm, sys%st, sys%ions, sys%ext_partners)
643
644 call boundaries_set(sys%gr%der%boundaries, sys%gr, sys%st%group%psib(1, 1))
645
646 call sys%st%group%psib(1, 1)%copy_to(hpsib)
647
648 if (sys%hm%apply_packed()) then
649 call sys%st%group%psib(1, 1)%do_pack()
650 call hpsib%do_pack(copy = .false.)
651 end if
652
653 do itime = 1, param%repetitions
654 if (states_are_real(sys%st)) then
655 call dhamiltonian_elec_apply_batch(sys%hm, sys%namespace, sys%gr, sys%st%group%psib(1, 1), hpsib, terms = terms, &
656 set_bc = .false.)
657 else
658 call zhamiltonian_elec_apply_batch(sys%hm, sys%namespace, sys%gr, sys%st%group%psib(1, 1), hpsib, terms = terms, &
659 set_bc = .false.)
660 end if
661 end do
662
663 if (hpsib%is_packed()) then
664 call hpsib%do_unpack(force = .true.)
665 end if
666
667 call test_prints_info_batch(sys%st, sys%gr, hpsib)
668
669 call hpsib%end(copy = .false.)
670 call states_elec_deallocate_wfns(sys%st)
671 safe_deallocate_p(sys)
672
673 pop_sub(test_hamiltonian)
674 end subroutine test_hamiltonian
675
676
677 ! ---------------------------------------------------------
678 subroutine test_density_calc(param, namespace)
679 type(test_parameters_t), intent(in) :: param
680 type(namespace_t), intent(in) :: namespace
681
682 type(electrons_t), pointer :: sys
683 integer :: itime
684
685 push_sub(test_density_calc)
686
687 call calc_mode_par%set_parallelization(p_strategy_states, default = .false.)
688
689 call messages_write('Info: Testing density calculation')
690 call messages_new_line()
692 call messages_info(namespace=namespace)
693
694 sys => electrons_t(namespace, generate_epot=.false.)
695 call sys%init_parallelization(mpi_world)
696
697 call states_elec_allocate_wfns(sys%st, sys%gr)
698 call test_batch_set_gaussian(sys%st%group%psib(1, 1), sys%gr)
699 if (sys%st%pack_states) call sys%st%pack()
700
701 do itime = 1, param%repetitions
702 call density_calc(sys%st, sys%gr, sys%st%rho)
703 end do
704
705 write(message(1),'(a,3x, f12.6)') "Norm density ", dmf_nrm2(sys%gr, sys%st%rho(:,1))
706 call messages_info(1, namespace=namespace)
707
708 call states_elec_deallocate_wfns(sys%st)
709 safe_deallocate_p(sys)
710
711 pop_sub(test_density_calc)
712 end subroutine test_density_calc
713
714
715 ! ---------------------------------------------------------
716 subroutine test_boundaries(param, namespace)
717 type(test_parameters_t), intent(in) :: param
718 type(namespace_t), intent(in) :: namespace
719
720 type(electrons_t), pointer :: sys
721 integer :: itime
722
723 push_sub(test_density_calc)
724
725 call calc_mode_par%set_parallelization(p_strategy_states, default = .false.)
726
727 call messages_write('Info: Testing boundary conditions')
728 call messages_new_line()
729 call messages_new_line()
730 call messages_info(namespace=namespace)
731
732 sys => electrons_t(namespace, generate_epot=.false.)
733 call sys%init_parallelization(mpi_world)
734
735 call states_elec_allocate_wfns(sys%st, sys%gr)
736 call test_batch_set_gaussian(sys%st%group%psib(1, 1), sys%gr)
737 if (sys%st%pack_states) call sys%st%pack()
738
739 do itime = 1, param%repetitions
740 call boundaries_set(sys%gr%der%boundaries, sys%gr, sys%st%group%psib(1, 1))
741 end do
742
743 call test_prints_info_batch(sys%st, sys%gr, sys%st%group%psib(1, 1))
744
745 call states_elec_deallocate_wfns(sys%st)
746 safe_deallocate_p(sys)
747
748 pop_sub(test_density_calc)
749 end subroutine test_boundaries
750
751 ! ---------------------------------------------------------
755 subroutine test_composition_chebyshev(namespace)
756 type(namespace_t), intent(in) :: namespace
757
758 type(electrons_t), pointer :: sys
759 type(chebyshev_filter_bounds_t), pointer :: bounds
760 integer, allocatable :: degree(:)
761 integer :: degree_n
762 type(wfs_elec_t) :: psib
763
765
766 call calc_mode_par%set_parallelization(p_strategy_states, default = .false.)
767
768 call messages_write('Info: Testing Chebyshev filtering and its composition rule')
769 call messages_new_line()
770 call messages_new_line()
771 call messages_info(namespace=namespace)
772
773 sys => electrons_t(namespace, generate_epot=.false.)
774 call sys%init_parallelization(mpi_world)
775
776 call states_elec_allocate_wfns(sys%st, sys%gr)
777
778 if (sys%st%pack_states .and. sys%hm%apply_packed()) call sys%st%pack()
779 call hamiltonian_elec_epot_generate(sys%hm, sys%namespace, sys%space, sys%gr, sys%ions, sys%ext_partners, sys%st)
780 call density_calc(sys%st, sys%gr, sys%st%rho)
781 call v_ks_calc(sys%ks, sys%namespace, sys%space, sys%hm, sys%st, sys%ions, sys%ext_partners)
782
783 bounds => chebyshev_filter_bounds_t(-1.0_real64, 1.0_real64)
784 safe_allocate(degree(1:sys%st%group%nblocks))
785
786 !%Variable TestCompositionOrder
787 !%Type integer
788 !%Section Calculation Modes::Test
789 !%Description
790 !% Specifies the value of m for the test of the composition rule
791 !% $T_{2n}(x) = 2*T_n(x)^2-1$.
792 !%
793 !%End
794 call parse_variable(namespace, 'TestCompositionOrder', 2, degree_n)
795
796 !Initialization
797 call test_batch_set_gaussian(sys%st%group%psib(1, 1), sys%gr)
798 call messages_new_line()
799 call messages_write('Info: Result for the original vector')
800 call messages_info(namespace=namespace)
801 call test_prints_info_batch(sys%st, sys%gr, sys%st%group%psib(1, 1))
802
803 ! Direct evaluation of mxn-th order Chebyshev filtering
804 degree = 2*degree_n
805 if (states_are_real(sys%st)) then
806 call dchebyshev_filter(namespace, sys%gr, sys%st, sys%hm, degree, bounds, 1, normalize=.false.)
807 else
808 call zchebyshev_filter(namespace, sys%gr, sys%st, sys%hm, degree, bounds, 1, normalize=.false.)
809 end if
810
812 call messages_write('Info: Result after calling 2n-th order filtering')
813 call messages_info(namespace=namespace)
814 call test_prints_info_batch(sys%st, sys%gr, sys%st%group%psib(1, 1))
815
816 ! Composition rule: m-th order Chebyshev filtering followed by n-th order filtering
817 call test_batch_set_gaussian(sys%st%group%psib(1, 1), sys%gr)
818 call sys%st%group%psib(1, 1)%copy_to(psib, copy_data=.true.)
819 degree = degree_n
820 if (states_are_real(sys%st)) then
821 call dchebyshev_filter(namespace, sys%gr, sys%st, sys%hm, degree, bounds, 1, normalize=.false.)
822 call dchebyshev_filter(namespace, sys%gr, sys%st, sys%hm, degree, bounds, 1, normalize=.false.)
823 else
824 call zchebyshev_filter(namespace, sys%gr, sys%st, sys%hm, degree, bounds, 1, normalize=.false.)
825 call zchebyshev_filter(namespace, sys%gr, sys%st, sys%hm, degree, bounds, 1, normalize=.false.)
826 end if
827 call batch_scal(sys%gr%np, m_two, sys%st%group%psib(1, 1))
828 call batch_axpy(sys%gr%np, -m_one, psib, sys%st%group%psib(1, 1))
829
830 call messages_new_line()
831 call messages_write('Info: Result using the composition rule')
832 call messages_info(namespace=namespace)
833 call test_prints_info_batch(sys%st, sys%gr, sys%st%group%psib(1, 1))
834
835
836 safe_deallocate_a(degree)
837 safe_deallocate_p(bounds)
838
839 call psib%end()
840 call states_elec_deallocate_wfns(sys%st)
841 safe_deallocate_p(sys)
842
844 end subroutine test_composition_chebyshev
845
846
847
848 ! ---------------------------------------------------------
849 subroutine test_exponential(param, namespace)
850 type(test_parameters_t), intent(in) :: param
851 type(namespace_t), intent(in) :: namespace
852
853 type(electrons_t), pointer :: sys
854 type(exponential_t) :: te
855 integer :: itime
856
857 push_sub(test_exponential)
858
859 call calc_mode_par%set_parallelization(p_strategy_states, default = .false.)
860
861 call messages_write('Info: Testing exponential')
862 call messages_new_line()
863 call messages_new_line()
864 call messages_info(namespace=namespace)
865
866 sys => electrons_t(namespace, generate_epot=.false.)
867 call sys%init_parallelization(mpi_world)
868
869 call states_elec_allocate_wfns(sys%st, sys%gr, wfs_type=type_cmplx)
870 call test_batch_set_gaussian(sys%st%group%psib(1, 1), sys%gr)
871
872 ! Initialize external potential
873 if (sys%st%pack_states .and. sys%hm%apply_packed()) call sys%st%pack()
874 call hamiltonian_elec_epot_generate(sys%hm, sys%namespace, sys%space, sys%gr, sys%ions, sys%ext_partners, sys%st)
875 call density_calc(sys%st, sys%gr, sys%st%rho)
876 call v_ks_calc(sys%ks, sys%namespace, sys%space, sys%hm, sys%st, sys%ions, sys%ext_partners)
877
878 call exponential_init(te, namespace)
879
880 if (sys%hm%apply_packed()) then
881 call sys%st%group%psib(1, 1)%do_pack()
882 end if
883
884 do itime = 1, param%repetitions
885 call te%apply_batch(sys%namespace, sys%gr, sys%hm, sys%st%group%psib(1, 1), 1e-3_real64)
886 end do
887
888 call test_prints_info_batch(sys%st, sys%gr, sys%st%group%psib(1, 1))
889
890 call states_elec_deallocate_wfns(sys%st)
891 safe_deallocate_p(sys)
892
893 pop_sub(test_exponential)
894 end subroutine test_exponential
895
896
897 ! ---------------------------------------------------------
898 subroutine test_subspace_diagonalization(param, namespace)
899 type(test_parameters_t), intent(in) :: param
900 type(namespace_t), intent(in) :: namespace
901
902 type(electrons_t), pointer :: sys
903 integer :: itime
904 type(subspace_t) :: sdiag
905 real(real64), allocatable :: diff(:)
906
908
909 call calc_mode_par%set_parallelization(p_strategy_states, default = .false.)
910
911 call messages_write('Info: Testing subspace diagonalization')
912 call messages_new_line()
913 call messages_new_line()
914 call messages_info(namespace=namespace)
915
916 sys => electrons_t(namespace, generate_epot=.false.)
917 call sys%init_parallelization(mpi_world)
918
919 call states_elec_allocate_wfns(sys%st, sys%gr)
920 call test_batch_set_gaussian(sys%st%group%psib(1, 1), sys%gr)
921
922 if (sys%st%pack_states .and. sys%hm%apply_packed()) call sys%st%pack()
923 call hamiltonian_elec_epot_generate(sys%hm, sys%namespace, sys%space, sys%gr, sys%ions, sys%ext_partners, sys%st)
924 call density_calc(sys%st, sys%gr, sys%st%rho)
925 call v_ks_calc(sys%ks, sys%namespace, sys%space, sys%hm, sys%st, sys%ions, sys%ext_partners)
926
927 call sdiag%init(sys%namespace, sys%st)
928
929 safe_allocate(diff(1:sys%st%nst))
930
931 do itime = 1, param%repetitions
932 if (states_are_real(sys%st)) then
933 call dsubspace_diag(sdiag, sys%namespace, sys%gr, sys%st, sys%hm, 1, sys%st%eigenval(:, 1), diff)
934 else
935 call zsubspace_diag(sdiag, sys%namespace, sys%gr, sys%st, sys%hm, 1, sys%st%eigenval(:, 1), diff)
936 end if
937 end do
938
939 safe_deallocate_a(diff)
940
941 call test_prints_info_batch(sys%st, sys%gr, sys%st%group%psib(1, 1))
942
943 call states_elec_deallocate_wfns(sys%st)
944 safe_deallocate_p(sys)
945
947 end subroutine test_subspace_diagonalization
948
949
950 ! ---------------------------------------------------------
951 subroutine test_batch_ops(param, namespace)
952 type(test_parameters_t), intent(in) :: param
953 type(namespace_t), intent(in) :: namespace
954
955 type(electrons_t), pointer :: sys
956 integer :: itime, ops, ops_default, ist, jst, nst
957 type(wfs_elec_t) :: xx, yy
958 real(real64), allocatable :: tmp(:)
959 real(real64), allocatable :: ddotv(:)
960 complex(real64), allocatable :: zdotv(:)
961 real(real64), allocatable :: ddot(:,:), df(:,:), dweight(:), dpoints(:,:,:)
962 complex(real64), allocatable :: zdot(:,:), zf(:,:), zweight(:), zpoints(:,:,:)
963 integer :: sp, block_size, size
964
965
966 push_sub(test_batch_ops)
967
968 !%Variable TestBatchOps
969 !%Type flag
970 !%Default all
971 !%Section Calculation Modes::Test
972 !%Description
973 !% Decides which part of the Hamiltonian is applied.
974 !%Option ops_axpy bit(1)
975 !% Tests batch_axpy operation
976 !%Option ops_scal bit(2)
977 !% Tests batch_scal operation
978 !%Option ops_nrm2 bit(3)
979 !% Tests batch_nrm2 operation
980 !%Option ops_dotp_matrix bit(4)
981 !% Tests X(mesh_batch_dotp_matrix)
982 !%Option ops_dotp_self bit(5)
983 !% Tests X(mesh_batch_dotp_self)
984 !%Option ops_dotp_vector bit(6)
985 !% Tests X(mesh_batch_dotp_vector)
986 !%Option ops_ax_function_py bit(7)
987 !% Tests X(batch_ax_function_py)
988 !%Option ops_get_points bit(8)
989 !% Tests batch_get_points and batch_set_points
990 !%End
991 ops_default = &
992 option__testbatchops__ops_axpy + &
993 option__testbatchops__ops_scal + &
994 option__testbatchops__ops_nrm2 + &
995 option__testbatchops__ops_dotp_matrix + &
996 option__testbatchops__ops_dotp_self + &
997 option__testbatchops__ops_dotp_vector + &
998 option__testbatchops__ops_ax_function_py + &
999 option__testbatchops__ops_get_points
1000
1001 call parse_variable(namespace, 'TestBatchOps', ops_default, ops)
1002
1003 call calc_mode_par%set_parallelization(p_strategy_states, default = .false.)
1004
1005 call messages_write('Info: Testing batch operations')
1006 call messages_new_line()
1007 call messages_new_line()
1008 call messages_info(namespace=namespace)
1009
1010 sys => electrons_t(namespace, generate_epot=.false.)
1011 call sys%init_parallelization(mpi_world)
1012
1013 call states_elec_allocate_wfns(sys%st, sys%gr)
1014 call test_batch_set_gaussian(sys%st%group%psib(1, 1), sys%gr)
1015 if (sys%st%pack_states) call sys%st%pack()
1016
1017 if (bitand(ops, option__testbatchops__ops_axpy) /= 0) then
1018 message(1) = 'Info: Testing axpy'
1019 call messages_info(1, namespace=sys%namespace)
1020
1021 call sys%st%group%psib(1, 1)%copy_to(xx, copy_data = .true.)
1022 call sys%st%group%psib(1, 1)%copy_to(yy, copy_data = .true.)
1023
1024 do itime = 1, param%repetitions
1025 call batch_axpy(sys%gr%np, 0.1_real64, xx, yy)
1026 end do
1027 call test_prints_info_batch(sys%st, sys%gr, yy, string = "axpy")
1028
1029 call xx%end()
1030 call yy%end()
1031 end if
1032
1033 if (bitand(ops, option__testbatchops__ops_scal) /= 0) then
1034 message(1) = 'Info: Testing scal'
1035 call messages_info(1, namespace=sys%namespace)
1036
1037 call sys%st%group%psib(1, 1)%copy_to(xx, copy_data = .true.)
1038 call sys%st%group%psib(1, 1)%copy_to(yy, copy_data = .true.)
1039
1040 do itime = 1, param%repetitions
1041 call batch_scal(sys%gr%np, 0.1_real64, yy)
1042 end do
1043 call test_prints_info_batch(sys%st, sys%gr, yy, string="scal")
1044
1045 call xx%end()
1046 call yy%end()
1047 end if
1048
1049 if (bitand(ops, option__testbatchops__ops_nrm2) /= 0) then
1050 message(1) = 'Info: Testing nrm2'
1051 call messages_info(1, namespace=sys%namespace)
1052
1053 call sys%st%group%psib(1, 1)%copy_to(xx, copy_data = .true.)
1054 call sys%st%group%psib(1, 1)%copy_to(yy, copy_data = .true.)
1055
1056 safe_allocate(tmp(1:xx%nst))
1057
1058 do itime = 1, param%repetitions
1059 call mesh_batch_nrm2(sys%gr, yy, tmp)
1060 end do
1061 do itime = 1, xx%nst
1062 write(message(1),'(a,i1,3x,e23.16)') "Nrm2 norm state ", itime, tmp(itime)
1063 call messages_info(1, namespace=sys%namespace)
1064 end do
1065
1066 safe_deallocate_a(tmp)
1067
1068 call xx%end()
1069 call yy%end()
1070 end if
1071
1072 if (bitand(ops, option__testbatchops__ops_dotp_matrix) /= 0) then
1073
1074 message(1) = 'Info: Testing dotp_matrix'
1075 call messages_info(1, namespace=sys%namespace)
1076
1077 call sys%st%group%psib(1, 1)%copy_to(xx, copy_data = .true.)
1078 call sys%st%group%psib(1, 1)%copy_to(yy, copy_data = .true.)
1079
1080 nst = sys%st%group%psib(1, 1)%nst
1081
1082 if (states_are_real(sys%st)) then
1083 safe_allocate(ddot(1:nst, 1:nst))
1084 do itime = 1, param%repetitions
1085 call dmesh_batch_dotp_matrix(sys%gr, xx, yy, ddot)
1086 end do
1087
1088 do ist = 1, nst
1089 do jst = 1, nst
1090 write(message(jst), '(a,2i3,3x,e23.16)') 'Dotp_matrix states', ist, jst, ddot(ist,jst)
1091 end do
1092 call messages_info(nst, namespace=sys%namespace)
1093 end do
1094 safe_deallocate_a(ddot)
1095 else
1096 safe_allocate(zdot(1:nst, 1:nst))
1097 do itime = 1, param%repetitions
1098 call zmesh_batch_dotp_matrix(sys%gr, xx, yy, zdot)
1099 end do
1100
1101 do ist = 1, nst
1102 do jst = 1, nst
1103 write(message(jst), '(a,2i3,3x,2(e23.16,1x))') 'Dotp_matrix states', ist, jst, zdot(ist,jst)
1104 end do
1105 call messages_info(nst, namespace=sys%namespace)
1106 end do
1107 safe_deallocate_a(zdot)
1108 end if
1109
1110 call xx%end()
1111 call yy%end()
1112 end if
1113
1114 if (bitand(ops, option__testbatchops__ops_dotp_vector) /= 0) then
1115
1116 message(1) = 'Info: Testing dotp_vector'
1117 call messages_info(1, namespace=sys%namespace)
1118
1119 call sys%st%group%psib(1, 1)%copy_to(xx, copy_data = .true.)
1120 call sys%st%group%psib(1, 1)%copy_to(yy, copy_data = .true.)
1121
1122 nst = sys%st%group%psib(1, 1)%nst
1123
1124 if (states_are_real(sys%st)) then
1125 safe_allocate(ddotv(1:nst))
1126 do itime = 1, param%repetitions
1127 call dmesh_batch_dotp_vector(sys%gr, xx, yy, ddotv)
1128 end do
1129
1130 do ist = 1, nst
1131 write(message(ist), '(a,i3,3x,e23.16)') 'Dotp_vector state', ist, ddotv(ist)
1132 end do
1133 call messages_info(nst, namespace=sys%namespace)
1134 safe_deallocate_a(ddotv)
1135 else
1136 safe_allocate(zdotv(1:nst))
1137 do itime = 1, param%repetitions
1138 call zmesh_batch_dotp_vector(sys%gr, xx, yy, zdotv)
1139 end do
1140
1141 do ist = 1, nst
1142 write(message(ist), '(a,i3,3x,2(e23.16,1x))') 'Dotp_vector state', ist, zdotv(ist)
1143 end do
1144 call messages_info(nst, namespace=sys%namespace)
1145 safe_deallocate_a(zdotv)
1146 end if
1147
1148 call xx%end()
1149 call yy%end()
1150 end if
1151
1152 if (bitand(ops, option__testbatchops__ops_dotp_self) /= 0) then
1153
1154 message(1) = 'Info: Testing dotp_self'
1155 call messages_info(1, namespace=sys%namespace)
1156
1157 call sys%st%group%psib(1, 1)%copy_to(xx, copy_data = .true.)
1158
1159 nst = sys%st%group%psib(1, 1)%nst
1160
1161 if (states_are_real(sys%st)) then
1162 safe_allocate(ddot(1:nst, 1:nst))
1163 do itime = 1, param%repetitions
1164 call dmesh_batch_dotp_self(sys%gr, xx, ddot)
1165 end do
1166
1167 do ist = 1, nst
1168 do jst = 1, nst
1169 write(message(jst), '(a,2i3,3x,e23.16)') 'Dotp_self states', ist, jst, ddot(ist,jst)
1170 end do
1171 call messages_info(nst, namespace=sys%namespace)
1172 end do
1173 safe_deallocate_a(ddot)
1174 else
1175 safe_allocate(zdot(1:nst, 1:nst))
1176 do itime = 1, param%repetitions
1177 call zmesh_batch_dotp_self(sys%gr, xx, zdot)
1178 end do
1179
1180 do ist = 1, nst
1181 do jst = 1, nst
1182 write(message(jst), '(a,2i3,3x,2(e23.16,1x))') 'Dotp_self states', ist, jst, zdot(ist,jst)
1183 end do
1184 call messages_info(nst, namespace=sys%namespace)
1185 end do
1186 safe_deallocate_a(zdot)
1187 end if
1188
1189 call xx%end()
1190 end if
1191
1192 if (bitand(ops, option__testbatchops__ops_ax_function_py) /= 0) then
1193 message(1) = 'Info: Testing ax_function_py'
1194 call messages_info(1, namespace=sys%namespace)
1195
1196 call sys%st%group%psib(1, 1)%copy_to(yy, copy_data = .true.)
1197
1198 if (states_are_real(sys%st)) then
1199 safe_allocate(df(sys%gr%np, sys%st%d%dim))
1200 safe_allocate(dweight(1:sys%st%group%psib(1, 1)%nst_linear))
1201 dweight = 0.1_real64
1202 call batch_get_state(sys%st%group%psib(1, 1), 1, sys%gr%np, df)
1203
1204 do itime = 1, param%repetitions
1205 call dbatch_ax_function_py(sys%gr%np, dweight, df, yy)
1206 end do
1207 safe_deallocate_a(df)
1208 safe_deallocate_a(dweight)
1209 else
1210 safe_allocate(zf(sys%gr%np, sys%st%d%dim))
1211 safe_allocate(zweight(1:sys%st%group%psib(1, 1)%nst_linear))
1212 zweight = cmplx(0.1_real64, m_zero, real64)
1213 call batch_get_state(sys%st%group%psib(1, 1), 1, sys%gr%np, zf)
1214
1215 do itime = 1, param%repetitions
1216 call zbatch_ax_function_py(sys%gr%np, zweight, zf, yy)
1217 end do
1218 safe_deallocate_a(zf)
1219 end if
1220 call test_prints_info_batch(sys%st, sys%gr, yy, string = "ax_function_py")
1221
1222 call yy%end()
1223 end if
1224
1225 if (bitand(ops, option__testbatchops__ops_get_points) /= 0) then
1226 ! The blocking over grid points mimics the actual use in the code
1227
1228 if (.not. accel_is_enabled()) then
1229
1230 call sys%st%group%psib(1, 1)%copy_to(yy)
1231
1232 if (states_are_real(sys%st)) then
1233
1234 block_size = max(40, cpu_hardware%l2%size/(2 * sizeof_real64 * sys%st%nst))
1235
1236 safe_allocate(dpoints(1:sys%st%nst, 1:sys%st%d%dim, 1:block_size))
1237
1238 do itime = 1, param%repetitions
1239 do sp = 1, sys%gr%np, block_size
1240 size = min(block_size, sys%gr%np - sp + 1)
1241 call batch_get_points(sys%st%group%psib(1, 1), sp, sp + size - 1, dpoints)
1242 call batch_set_points(yy, sp, sp + size - 1, dpoints)
1243 end do
1244 end do
1245 safe_deallocate_a(dpoints)
1246
1247 else
1248
1249 block_size = max(20, cpu_hardware%l2%size/(2 * sizeof_complex64 * sys%st%nst))
1250
1251 safe_allocate(zpoints(1:sys%st%nst, 1:sys%st%d%dim, 1:block_size))
1252
1253 do itime = 1, param%repetitions
1254 do sp = 1, sys%gr%np, block_size
1255 size = min(block_size, sys%gr%np - sp + 1)
1256 call batch_get_points(sys%st%group%psib(1, 1), sp, sp + size - 1, zpoints)
1257 call batch_set_points(yy, sp, sp + size - 1, zpoints)
1258 end do
1259 end do
1260 safe_deallocate_a(dpoints)
1261
1262 end if
1263
1264 call test_prints_info_batch(sys%st, sys%gr, yy, string = "get_points")
1265
1266 call yy%end()
1267 end if
1268 end if
1269
1270 call states_elec_deallocate_wfns(sys%st)
1271 safe_deallocate_p(sys)
1272
1273 pop_sub(test_batch_ops)
1274 end subroutine test_batch_ops
1275
1276! ---------------------------------------------------------
1277 subroutine test_derivatives(param, namespace)
1278 type(test_parameters_t), intent(in) :: param
1279 type(namespace_t), intent(in) :: namespace
1280
1281 type(electrons_t), pointer :: sys
1282
1283 push_sub(test_derivatives)
1284
1285 sys => electrons_t(namespace, generate_epot=.false.)
1286 call sys%init_parallelization(mpi_world)
1287
1288 message(1) = 'Info: Testing the finite-differences derivatives.'
1289 message(2) = ''
1290 call messages_info(2, namespace=namespace)
1291
1292 if (param%type == option__testtype__all .or. param%type == option__testtype__real) then
1293 call dderivatives_test(sys%gr%der, sys%namespace, param%repetitions, param%min_blocksize, param%max_blocksize)
1294 end if
1295
1296 if (param%type == option__testtype__all .or. param%type == option__testtype__complex) then
1297 call zderivatives_test(sys%gr%der, sys%namespace, param%repetitions, param%min_blocksize, param%max_blocksize)
1298 end if
1299
1300 safe_deallocate_p(sys)
1301
1302 pop_sub(test_derivatives)
1303 end subroutine test_derivatives
1304
1305 ! ---------------------------------------------------------
1306
1307 subroutine test_orthogonalization(param, namespace)
1308 type(test_parameters_t), intent(in) :: param
1309 type(namespace_t), intent(in) :: namespace
1310
1311 type(electrons_t), pointer :: sys
1312 integer :: itime
1313
1314 push_sub(test_orthogonalization)
1315
1316 call calc_mode_par%set_parallelization(p_strategy_states, default = .false.)
1317 call calc_mode_par%set_scalapack_compat()
1318
1319 sys => electrons_t(namespace, generate_epot=.false.)
1320 call sys%init_parallelization(mpi_world)
1321
1322 message(1) = 'Info: Testing orthogonalization.'
1323 message(2) = ''
1324 call messages_info(2, namespace=namespace)
1325
1326 if (param%type == option__testtype__all .or. param%type == option__testtype__real) then
1327 message(1) = 'Info: Real wave-functions.'
1328 call messages_info(1, namespace=namespace)
1329 do itime = 1, param%repetitions
1330 call dstates_elec_calc_orth_test(sys%st, sys%namespace, sys%gr, sys%kpoints)
1331 end do
1332 end if
1333
1334 if (param%type == option__testtype__all .or. param%type == option__testtype__complex) then
1335 message(1) = 'Info: Complex wave-functions.'
1336 call messages_info(1, namespace=namespace)
1337 do itime = 1, param%repetitions
1338 call zstates_elec_calc_orth_test(sys%st, sys%namespace, sys%gr, sys%kpoints)
1339 end do
1340 end if
1341
1342 safe_deallocate_p(sys)
1343
1344 pop_sub(test_orthogonalization)
1345 end subroutine test_orthogonalization
1346
1347 ! ---------------------------------------------------------
1348
1349 subroutine test_interpolation(param, namespace)
1350 type(test_parameters_t), intent(in) :: param
1351 type(namespace_t), intent(in) :: namespace
1352
1353 type(electrons_t), pointer :: sys
1354
1355 push_sub(test_interpolation)
1356
1357 sys => electrons_t(namespace, generate_epot=.false.)
1358 call sys%init_parallelization(mpi_world)
1359
1360 if (param%type == option__testtype__all .or. param%type == option__testtype__real) then
1361 call messages_write('Info: Testing real interpolation routines')
1362 call messages_new_line()
1363 call messages_new_line()
1364 call messages_info(namespace=namespace)
1365
1366 call dmesh_interpolation_test(sys%gr)
1367 end if
1368
1369 if (param%type == option__testtype__all .or. param%type == option__testtype__complex) then
1370 call messages_new_line()
1371 call messages_write('Info: Testing complex interpolation routines')
1373 call messages_new_line()
1374 call messages_info(namespace=namespace)
1375
1376 call zmesh_interpolation_test(sys%gr)
1377 end if
1378
1379 safe_deallocate_p(sys)
1380
1381 pop_sub(test_interpolation)
1382 end subroutine test_interpolation
1383
1384
1385 ! ---------------------------------------------------------
1386
1387 subroutine test_ion_interaction(namespace)
1388 type(namespace_t), intent(in) :: namespace
1389
1390 type(electrons_t), pointer :: sys
1391
1392 push_sub(test_ion_interaction)
1393
1394 sys => electrons_t(namespace, generate_epot=.false.)
1395 call sys%init_parallelization(mpi_world)
1396
1397 call ion_interaction_test(sys%space, sys%ions%latt, sys%ions%atom, sys%ions%natoms, sys%ions%pos, &
1398 sys%gr%box%bounding_box_l, namespace, sys%mc)
1399
1400 safe_deallocate_p(sys)
1401
1403 end subroutine test_ion_interaction
1404
1405 ! ---------------------------------------------------------
1406
1407 subroutine test_prints_info_batch(st, gr, psib, string)
1408 type(states_elec_t), intent(in) :: st
1409 type(grid_t), intent(in) :: gr
1410 class(batch_t), intent(inout) :: psib
1411 character(*), optional, intent(in) :: string
1412
1413 integer :: itime
1414 complex(real64), allocatable :: zpsi(:, :)
1415 real(real64), allocatable :: dpsi(:, :)
1416
1417 character(80) :: string_
1418
1419 string_ = optional_default(string, "")
1420
1421 push_sub(test_prints_info_batch)
1422
1423 if (states_are_real(st)) then
1424 safe_allocate(dpsi(1:gr%np, 1:st%d%dim))
1425 else
1426 safe_allocate(zpsi(1:gr%np, 1:st%d%dim))
1427 end if
1428
1429 do itime = 1, psib%nst
1430 if (states_are_real(st)) then
1431 call batch_get_state(psib, itime, gr%np, dpsi)
1432 write(message(1),'(a,i2,3x,e23.16)') "Norm state "//trim(string_)//" ", itime, dmf_nrm2(gr, st%d%dim, dpsi)
1433 else
1434 call batch_get_state(psib, itime, gr%np, zpsi)
1435 write(message(1),'(a,i2,3x,e23.16)') "Norm state "//trim(string_)//" ", itime, zmf_nrm2(gr, st%d%dim, zpsi)
1436 end if
1437 call messages_info(1)
1438 end do
1439
1440 if (states_are_real(st)) then
1441 safe_deallocate_a(dpsi)
1442 else
1443 safe_deallocate_a(zpsi)
1444 end if
1445
1446 pop_sub(test_prints_info_batch)
1447
1448 end subroutine test_prints_info_batch
1449
1450
1451 ! ---------------------------------------------------------
1452 subroutine test_clock()
1453
1454 type(iteration_counter_t) :: clock, other_clock
1455
1456 push_sub(test_clock)
1457
1458 ! Clock operations
1459 write(message(1),'(a)') " Operation Counter Time Global step"
1460 call messages_info(1)
1461
1462 clock = clock_t(time_step=m_two, initial_iteration=2)
1463 call write_clock("initialization")
1464
1465 clock = clock + 5
1466 call write_clock("addition (+5)")
1467
1468 clock = clock - 3
1469 call write_clock("subtraction (-3)")
1470
1471 call clock%reset()
1472 call write_clock("reset (=0)")
1473
1474 call clock%set(clock_t(time_step=m_four, initial_iteration=1))
1475 call write_clock("set (=2)")
1476
1477 ! Clock comparisons:
1478 write(message(1),'(a)') " Clock comparisons:"
1479 call messages_info(1)
1480
1481 other_clock = clock_t(time_step=m_one, initial_iteration=5)
1483 call write_condition_result("4 < 5", clock < other_clock)
1484 call write_condition_result("4 <= 5", clock <= other_clock)
1485 call write_condition_result("4 > 5", clock > other_clock)
1486 call write_condition_result("4 >= 5", clock >= other_clock)
1487 call write_condition_result("4 == 5", clock == other_clock)
1488 call write_condition_result("4 /= 5", clock /= other_clock)
1489
1490 pop_sub(test_clock)
1491 contains
1492 subroutine write_clock(operation)
1493 character(len=*), intent(in) :: operation
1494
1495 write(message(1),'(a17,1x,i6,1x,f10.1,1x,i16)') operation, clock%counter(), clock%value(), clock%global_step()
1496 call messages_info(1)
1497 end subroutine write_clock
1498
1499 subroutine write_condition_result(condition, result)
1500 character(len=*), intent(in) :: condition
1501 logical, intent(in) :: result
1503 write(message(1),'(a10," = ",i1," (",l1,")")') condition, abs(transfer(result, 0)), result
1504 call messages_info(1)
1505 end subroutine write_condition_result
1506
1507 end subroutine test_clock
1508
1509
1510 ! ---------------------------------------------------------
1511 subroutine test_cgal()
1512
1513 type(cgal_polyhedra_t) :: cgal_poly
1514
1515 push_sub(test_cgal)
1516
1517 call cgal_polyhedron_init(cgal_poly, "28-cgal.02-X.off", verbose = .true.)
1518
1519 if (cgal_polyhedron_point_inside(cgal_poly, 30._real64, 10._real64, 30._real64)) then
1520 message(1) = "cgal_polyhedron_point_inside"
1521 call messages_info(1)
1522 end if
1523
1524 call cgal_polyhedron_end(cgal_poly)
1525
1526 pop_sub(test_cgal)
1527 end subroutine test_cgal
1528
1529
1530 ! ---------------------------------------------------------
1531 subroutine test_dense_eigensolver()
1532 integer :: N, ii, jj, N_list(4), i_N
1533 real(real64), allocatable :: matrix(:, :), eigenvectors(:, :), eigenvalues(:), test(:)
1534 real(real64), allocatable :: differences(:)
1535
1536 push_sub(test_dense_eigensolver)
1537
1538 n_list = [15, 32, 100, 500]
1539
1540 do i_n = 1, 4
1541 n = n_list(i_n)
1542 safe_allocate(matrix(1:n, 1:n))
1543 safe_allocate(eigenvectors(1:n, 1:n))
1544 safe_allocate(eigenvalues(1:n))
1545 safe_allocate(test(1:n))
1546 safe_allocate(differences(1:n))
1548
1549 ! set up symmetrix matrix
1550 do jj = 1, n
1551 do ii = 1, n
1552 matrix(ii, jj) = ii * jj
1553 end do
1554 end do
1555
1556 ! parallel solver
1557 eigenvectors(1:n, 1:n) = matrix(1:n, 1:n)
1558 call lalg_eigensolve_parallel(n, eigenvectors, eigenvalues)
1559
1560 do ii = 1, n
1561 test(:) = matmul(matrix, eigenvectors(:, ii)) - eigenvalues(ii) * eigenvectors(:, ii)
1562 differences(ii) = sum(abs(test)) / sum(abs(eigenvectors(:, ii)))
1563 end do
1564 write(message(1), "(A, I3, A, E13.6)") "Parallel solver - N: ", n, &
1565 ", average difference: ", sum(differences)/n
1566 call messages_info(1)
1567
1568 ! serial solver
1569 eigenvectors(1:n, 1:n) = matrix(1:n, 1:n)
1570 call lalg_eigensolve(n, eigenvectors, eigenvalues)
1571
1572 do ii = 1, n
1573 test(:) = matmul(matrix, eigenvectors(:, ii)) - eigenvalues(ii) * eigenvectors(:, ii)
1574 differences(ii) = sum(abs(test)) / sum(abs(eigenvectors(:, ii)))
1575 end do
1576 write(message(1), "(A, I3, A, E13.6)") "Serial solver - N: ", n, &
1577 ", average difference: ", sum(differences)/n
1578 call messages_info(1)
1579
1580 safe_deallocate_a(matrix)
1581 safe_deallocate_a(eigenvectors)
1582 safe_deallocate_a(eigenvalues)
1583 safe_deallocate_a(test)
1584 safe_deallocate_a(differences)
1585 end do
1586
1588 end subroutine test_dense_eigensolver
1589
1590 subroutine test_batch_set_gaussian(psib, mesh)
1591 class(batch_t), intent(inout) :: psib
1592 class(mesh_t), intent(in) :: mesh
1593
1594 real(real64), allocatable :: dff(:)
1595 complex(real64), allocatable :: zff(:)
1596 integer :: ist, ip
1597 real(real64) :: da, db, dc
1598 complex(real64) :: za, zb, zc
1599
1600 push_sub(test_batch_set_gaussian)
1601
1602 ! use a similar function as in the derivatives test
1603 da = m_one/mesh%box%bounding_box_l(1)
1604 db = 10.0_real64
1605 dc = 100.0_real64
1607 if (type_is_complex(psib%type())) then
1608 ! we make things more "complex"
1609 za = da + m_zi*0.01_real64
1610 zb = db*exp(m_zi*0.345_real64)
1611 zc = dc - m_zi*50.0_real64
1612
1613 safe_allocate(zff(1:mesh%np))
1614 do ist = 1, psib%nst_linear
1615 za = za * ist
1616 zb = zb / ist
1617 do ip = 1, mesh%np
1618 zff(ip) = zb*exp(-za*sum(mesh%x(ip, :)**2)) + zc
1619 end do
1620 call batch_set_state(psib, ist, mesh%np, zff)
1621 end do
1622 call zmesh_batch_normalize(mesh, psib)
1623 safe_deallocate_a(zff)
1624 else
1625 safe_allocate(dff(1:mesh%np))
1626 do ist = 1, psib%nst_linear
1627 da = da * ist
1628 db = db / ist
1629 do ip = 1, mesh%np
1630 dff(ip) = db*exp(-da*sum(mesh%x(ip, :)**2)) + dc
1631 end do
1632 call batch_set_state(psib, ist, mesh%np, dff)
1633 end do
1634 call dmesh_batch_normalize(mesh, psib)
1635 safe_deallocate_a(dff)
1636 end if
1637
1639 end subroutine test_batch_set_gaussian
1640
1641 ! ---------------------------------------------------------
1642 subroutine test_grid_interpolation()
1643 type(electrons_t), pointer :: sys
1644 type(multigrid_t) :: mgrid
1645
1646 push_sub(test_grid_interpolation)
1647
1648 sys => electrons_t(global_namespace, generate_epot=.false.)
1649 call sys%init_parallelization(mpi_world)
1650
1651 call multigrid_init(mgrid, global_namespace, sys%space, sys%gr, sys%gr%der, &
1652 sys%gr%stencil, sys%mc, nlevels=3)
1653
1654 call multigrid_test_interpolation(mgrid, sys%space)
1655
1656 call multigrid_end(mgrid)
1657
1658 safe_deallocate_p(sys)
1659
1661 end subroutine test_grid_interpolation
1662
1663 subroutine test_iihash()
1664
1665 type(iihash_t) :: h
1666 integer :: value
1667 logical :: found
1668
1669 call iihash_init(h)
1670 call iihash_insert(h, 1,5)
1671 call iihash_insert(h, 3,7)
1672
1673 value = iihash_lookup(h, 1, found)
1674 write(message(1),*) 'hash[1] :', found, value
1675 call messages_info()
1676
1677 value = iihash_lookup(h, 2, found)
1678 write(message(1),*) 'hash[2] :', found, value
1679 call messages_info()
1680
1681 value = iihash_lookup(h, 3, found)
1682 write(message(1),*) 'hash[3] :', found, value
1683 call messages_info()
1684
1686
1687 end subroutine test_iihash
1688
1689 subroutine test_sihash()
1690 type(sihash_t) :: h
1691 type(sihash_iterator_t) :: it
1692
1693 integer :: counter
1694 integer :: value, sum
1695 logical :: found
1696
1697 call sihash_init(h)
1698 call sihash_insert(h, "one",5)
1699 call sihash_insert(h, "three",7)
1700 call sihash_insert(h, "the answer", 42)
1701
1702 value = sihash_lookup(h, "one", found)
1703 write(message(1),*) 'hash["one"]: ', found, value
1704 call messages_info()
1705
1706 value = sihash_lookup(h, "two", found)
1707 write(message(1),*) 'hash["two"]: ', found, value
1708 call messages_info()
1709
1710 value = sihash_lookup(h, "three", found)
1711 write(message(1),*) 'hash["three"]: ', found, value
1712 call messages_info()
1713
1714 sum = 0
1715 counter = 1
1716 call it%start(h)
1717
1718 do while (it%has_next())
1719 value = it%get_next()
1720 sum = sum + value
1721 write(message(1),'(I3,A,I5)') counter,': hash[...] = ',value
1722 call messages_info()
1723 counter = counter + 1
1724 end do
1725 write(message(1),*) 'counter = ', counter
1726 write(message(2),*) 'sum = ', sum
1727 call messages_info(2)
1728
1729
1730 call sihash_end(h)
1731
1732 end subroutine test_sihash
1733
1734
1735 subroutine test_sphash(namespace)
1736 type(namespace_t), intent(in) :: namespace
1738 type(sphash_t) :: h
1739 type(sphash_iterator_t) :: it
1740
1741 logical :: found
1742
1743
1744 type(iteration_counter_t) :: clock_1
1745 type(iteration_counter_t), allocatable :: clock_2
1746 type(space_t) :: space_1
1747
1748 class(*), pointer :: value
1749
1750 integer :: count_clock, count_space
1751
1752 safe_allocate(clock_2)
1753
1754 clock_1 = clock_t(1.d-5, 1)
1755 clock_2 = clock_t(2.d-5, 2)
1756 space_1 = space_t(namespace)
1757
1759 call sphash_insert(h, "one", clock_1)
1760 call sphash_insert(h, "two", space_1)
1761 call sphash_insert(h, "three", clock_2, clone=.true.)
1762
1763 value => sphash_lookup(h, "one", found)
1764 select type(value)
1765 type is (iteration_counter_t)
1766 write(message(1),*) 'hash["one"]: ', found, value%counter()
1767 call messages_info()
1768 type is (space_t)
1769 write(message(1),*) 'hash["one"]: ', found, value%short_info()
1770 call messages_info()
1771 class default
1772 write(message(1),*) 'wrong type. found = ', found
1773 call messages_info()
1774 end select
1775
1776 value => sphash_lookup(h, "two", found)
1777 select type(value)
1778 type is (iteration_counter_t)
1779 write(message(1),*) 'hash["two"]: ', found, value%counter()
1780 call messages_info()
1781 type is (space_t)
1782 write(message(1),*) 'hash["two"]: ', found, value%short_info()
1783 call messages_info()
1784 class default
1785 write(message(1),*) 'wrong type. found = ',found
1786 call messages_info()
1787 end select
1788
1789 safe_deallocate_a(clock_2)
1790
1791 value => sphash_lookup(h, "three", found)
1792 select type(value)
1793 type is (iteration_counter_t)
1794 write(message(1),*) 'hash["three"]: ', found, value%counter()
1795 call messages_info()
1796 type is (space_t)
1797 write(message(1),*) 'hash["three"]: ', found, value%short_info()
1798 call messages_info()
1799 class default
1800 write(message(1),*) 'wrong type. found = ',found
1801 call messages_info()
1802 end select
1803
1804 count_clock = 0
1805 count_space = 0
1806
1807 call it%start(h)
1808
1809 do while (it%has_next())
1810 value => it%get_next()
1811 select type(value)
1812 type is (iteration_counter_t)
1813 count_clock = count_clock + 1
1814 type is (space_t)
1815 count_space = count_space + 1
1816 end select
1817 end do
1818
1819 write(message(1), *) 'Count_clock = ', count_clock
1820 write(message(2), *) 'Count_space = ', count_space
1821 call messages_info(2)
1822
1823 call sphash_end(h)
1824
1825 end subroutine test_sphash
1826
1827 ! ---------------------------------------------------------
1828 subroutine test_regridding(namespace)
1829 type(namespace_t), intent(in) :: namespace
1831 type(electrons_t), pointer :: sysA, sysB
1832 type(regridding_t), pointer :: regridding
1833 real(real64), allocatable :: ff_A(:), ff_A_reference(:), ff_B(:), ff_B_reference(:), diff_A(:), diff_B(:)
1834 real(real64) :: norm_ff, norm_diff
1835 integer :: ip, ierr
1836
1837 push_sub(test_regridding)
1838
1839 sysa => electrons_t(namespace_t("A", namespace), generate_epot=.false.)
1840 sysb => electrons_t(namespace_t("B", namespace), generate_epot=.false.)
1841
1842 call calc_mode_par%set_parallelization(p_strategy_states, default=.true.)
1843 call sysa%init_parallelization(mpi_world)
1844 call sysb%init_parallelization(mpi_world)
1845
1846
1847 safe_allocate(ff_a(1:sysa%gr%np))
1848 safe_allocate(ff_a_reference(1:sysa%gr%np))
1849 safe_allocate(diff_a(1:sysa%gr%np))
1850 safe_allocate(ff_b(1:sysb%gr%np))
1851 safe_allocate(ff_b_reference(1:sysb%gr%np))
1852 safe_allocate(diff_b(1:sysb%gr%np))
1853
1854 do ip = 1, sysa%gr%np
1855 ff_a_reference(ip) = values(sysa%gr%x(ip, :))
1856 end do
1857 do ip = 1, sysb%gr%np
1858 ff_b_reference(ip) = values(sysb%gr%x(ip, :))
1859 end do
1860
1861 ! forward mapping A->B
1862 regridding => regridding_t(sysb%gr, sysa%gr, sysa%space, sysa%namespace)
1863 call regridding%do_transfer(ff_b, ff_a_reference)
1864 safe_deallocate_p(regridding)
1865
1866 ! check that mapped function is exactly the same for the points that are available on both grids
1867 do ip = 1, sysb%gr%np
1868 if (abs(ff_b(ip)) <= m_epsilon) then
1869 ff_b_reference(ip) = m_zero
1870 diff_b(ip) = m_zero
1871 else
1872 diff_b(ip) = abs(ff_b_reference(ip) - ff_b(ip))
1873 end if
1874 end do
1875 norm_ff = dmf_nrm2(sysb%gr, ff_b_reference)
1876 norm_diff = dmf_nrm2(sysb%gr, diff_b)
1877
1878 write(message(1),'(a, E14.6)') "Forward: difference of reference to mapped function (rel.): ", &
1879 norm_diff/norm_ff
1880 call messages_info(1, namespace=namespace)
1881
1882 call dio_function_output(io_function_fill_how('AxisX'), ".", "forward_reference", namespace, sysb%space, &
1883 sysb%gr, ff_b_reference, unit_one, ierr)
1884 call dio_function_output(io_function_fill_how('AxisX'), ".", "forward_mapped", namespace, sysb%space, &
1885 sysb%gr, ff_b, unit_one, ierr)
1886 call dio_function_output(io_function_fill_how('AxisX'), ".", "forward_original", namespace, sysa%space, &
1887 sysa%gr, ff_a_reference, unit_one, ierr)
1888
1889 ! backward mapping B->A
1890 regridding => regridding_t(sysa%gr, sysb%gr, sysb%space, sysb%namespace)
1891 call regridding%do_transfer(ff_a, ff_b_reference)
1892 safe_deallocate_p(regridding)
1893
1894 do ip = 1, sysa%gr%np
1895 if (abs(ff_a(ip)) <= m_epsilon) then
1896 ff_a_reference(ip) = m_zero
1897 diff_a(ip) = m_zero
1898 else
1899 diff_a(ip) = abs(ff_a_reference(ip) - ff_a(ip))
1900 end if
1901 end do
1902 norm_ff = dmf_nrm2(sysa%gr, ff_a_reference)
1903 norm_diff = dmf_nrm2(sysa%gr, diff_a)
1904
1905 write(message(1),'(a, E14.6)') "Backward: difference of reference to mapped function (rel.): ", &
1906 norm_diff/norm_ff
1907 call messages_info(1, namespace=namespace)
1908
1909 call dio_function_output(io_function_fill_how('AxisX'), ".", "backward_reference", namespace, sysa%space, &
1910 sysa%gr, ff_a_reference, unit_one, ierr)
1911 call dio_function_output(io_function_fill_how('AxisX'), ".", "backward_mapped", namespace, sysa%space, &
1912 sysa%gr, ff_a, unit_one, ierr)
1913 call dio_function_output(io_function_fill_how('AxisX'), ".", "backward_original", namespace, sysb%space, &
1914 sysb%gr, ff_b_reference, unit_one, ierr)
1915
1916 safe_deallocate_a(ff_a)
1917 safe_deallocate_a(ff_a_reference)
1918 safe_deallocate_a(ff_b)
1919 safe_deallocate_a(ff_b_reference)
1920 safe_deallocate_a(diff_a)
1921 safe_deallocate_a(diff_b)
1922 safe_deallocate_p(sysa)
1923 safe_deallocate_p(sysb)
1924
1925 call messages_info(1, namespace=namespace)
1926
1927 pop_sub(test_regridding)
1928 contains
1929 real(real64) function values(xx)
1930 real(real64), intent(in) :: xx(:)
1931 real(real64) :: xx0(1:size(xx, dim=1))
1932 real(real64), parameter :: aa = m_half
1933 real(real64), parameter :: bb = m_four
1934
1935 ! no push_sub/pop_sub because this function is called often
1936 xx0(:) = m_one
1937 values = bb * exp(-aa*sum((xx-xx0)**2))
1938 end function values
1939 end subroutine test_regridding
1940
1946 subroutine test_vecpot_analytical(namespace)
1947 type(namespace_t), intent(in) :: namespace
1948
1949 class(maxwell_t), pointer :: maxwell_system
1950
1951 real(real64), allocatable :: magnetic_field(:,:)
1952 real(real64), allocatable :: vector_potential_mag(:,:)
1953 real(real64), allocatable :: vector_potential_analytical(:,:)
1954 real(real64), allocatable :: delta(:,:)
1955 real(real64) :: exp_factor
1956 real(real64) :: xx
1957 real(real64) :: yy
1958 real(real64) :: zz
1959 real(real64) :: sigma
1960 integer :: ip, j, ierr, nn
1961 integer(int64) :: out_how
1962 character(len=MAX_PATH_LEN) :: fname, fname2, fname3
1963
1964 out_how = 32
1965 maxwell_system => maxwell_t(namespace)
1966 sigma = maxwell_system%gr%box%bounding_box_l(1)/10_real64 ! this is exponential width
1967 call maxwell_system%init_parallelization(mpi_world)
1968
1969 safe_allocate(magnetic_field(1:maxwell_system%gr%np_part, 1:3))
1970 safe_allocate(vector_potential_mag(1:maxwell_system%gr%np_part, 1:3))
1971 safe_allocate(vector_potential_analytical(1:maxwell_system%gr%np_part, 1:3))
1972 safe_allocate(delta(1:maxwell_system%gr%np, 1:3))
1973
1974 !%Variable TestVectorPotentialType
1975 !%Type integer
1976 !%Default bounded
1977 !%Section Calculation Modes::Test
1978 !%Description
1979 !% Select whether bounded or unbounded type will be used for vector potential tests
1980 !%Option bounded 1
1981 !% Analytical Vector Potential formulation is bounded by spatial gaussian
1982 !%Option unbounded 2
1983 !% Analytical Vector Potential is not bounded
1984 !%End
1985 call parse_variable(namespace, 'TestVectorPotentialType', option__testvectorpotentialtype__bounded, nn)
1986
1987 select case (nn)
1988 case (option__testvectorpotentialtype__bounded) ! bounded input
1989 do ip = 1, maxwell_system%gr%np_part
1990 xx = maxwell_system%gr%x(ip, 1)
1991 yy = maxwell_system%gr%x(ip, 2)
1992 zz = maxwell_system%gr%x(ip, 3)
1993 exp_factor = exp((-xx**2 - yy**2 - zz**2)*1/(2*sigma**2))
1994 magnetic_field(ip, 1) = exp_factor*yy*(1 - (-xx**2 + yy**2)/(3*sigma**2) - zz**2/(3*sigma**2))
1995 magnetic_field(ip, 2) = exp_factor * xx * (1 + (-xx**2 + yy**2)/(3*sigma**2) - zz**2/(3*sigma**2))
1996 magnetic_field(ip, 3) = exp_factor * 2 * xx * yy * zz * 1/(3*sigma**2)
1997
1998 vector_potential_analytical(ip, 1) = m_third * xx * zz * exp_factor
1999 vector_potential_analytical(ip, 2) = - m_third * yy * zz * exp_factor
2000 vector_potential_analytical(ip, 3) = m_third * (-xx**2 + yy**2) * exp_factor
2001 end do
2002 case (option__testvectorpotentialtype__unbounded) ! unbounded input,
2003 ! TODO this unit test requires implementation of BCs for Helmholtz decomposition
2004 do ip = 1, maxwell_system%gr%np_part
2005 magnetic_field(ip, 1) = maxwell_system%gr%x(ip, 2)
2006 magnetic_field(ip, 2) = maxwell_system%gr%x(ip, 1)
2007 magnetic_field(ip, 3) = m_zero
2008
2009 vector_potential_analytical(ip, 1) = m_third * maxwell_system%gr%x(ip, 1) * maxwell_system%gr%x(ip, 3)
2010 vector_potential_analytical(ip, 2) = - m_third * maxwell_system%gr%x(ip, 2) * maxwell_system%gr%x(ip, 3)
2011 vector_potential_analytical(ip, 3) = - m_third * (maxwell_system%gr%x(ip, 1)**2 - maxwell_system%gr%x(ip, 2)**2)
2012 end do
2013 end select
2014 call maxwell_system%helmholtz%get_vector_potential(namespace, vector_potential_mag, magnetic_field)
2015
2016 do j = 1, 3
2017 delta(:,:) = m_zero
2018 do ip = 1, maxwell_system%gr%np
2019 delta(ip,j) = vector_potential_analytical(ip, j) - vector_potential_mag(ip, j)
2020 end do
2021 end do
2022
2023 do j = 1, 3
2024 write(message(j),*) 'j, norm2(delta)', j, norm2(delta(:,j))
2025 end do
2026 call messages_info(3)
2027
2028 write(fname, '(a)') 'deviation_from_analytical_formulation' ! use messages later
2029 call io_function_output_vector(out_how , './' , trim(fname), namespace, maxwell_system%space, maxwell_system%gr, &
2030 delta, unit_one, ierr)
2031 write(fname2, '(a)') 'vector_potential_analytical'
2032 call io_function_output_vector(out_how , './' , trim(fname2), namespace, maxwell_system%space, maxwell_system%gr, &
2033 vector_potential_analytical, unit_one, ierr)
2034 write(fname3, '(a)') 'vector_potential_mag'
2035 call io_function_output_vector(out_how , './' , trim(fname3), namespace, maxwell_system%space, maxwell_system%gr, &
2036 vector_potential_mag, unit_one, ierr)
2037
2038 safe_deallocate_a(magnetic_field)
2039 safe_deallocate_a(vector_potential_mag)
2040 safe_deallocate_a(vector_potential_analytical)
2041 safe_deallocate_a(delta)
2042 safe_deallocate_p(maxwell_system)
2043
2044 end subroutine test_vecpot_analytical
2045
2046 ! ---------------------------------------------------------
2047 subroutine multigrid_test_interpolation(mgrid, space)
2048 type(multigrid_t), intent(in) :: mgrid
2049 class(space_t), intent(in) :: space
2050
2051 real(real64), allocatable :: guess0(:), res0(:), guess1(:)
2052 type(mesh_t), pointer :: mesh0, mesh1
2053 real(real64) :: delta, xx(3,2), alpha, beta, rr
2054 integer :: nn, ip, ierr
2055
2057
2058 message(1) = 'Info: Testing the grid interpolation.'
2059 message(2) = ''
2060 call messages_info(2)
2061
2062 mesh0 => mgrid%level(0)%mesh
2063 mesh1 => mgrid%level(1)%mesh
2064
2065 safe_allocate(guess0(1:mesh0%np_part))
2066 safe_allocate(res0(1:mesh0%np_part))
2067 safe_allocate(guess1(1:mesh1%np_part))
2068
2069 alpha = m_four*mesh0%spacing(1)
2070 beta = m_one / (alpha**space%dim * sqrt(m_pi)**space%dim)
2071
2072 ! Set the centers of the Gaussians by hand
2073 xx(1, 1) = m_one
2074 xx(2, 1) = -m_half
2075 xx(3, 1) = m_two
2076 xx(1, 2) = -m_two
2077 xx(2, 2) = m_zero
2078 xx(3, 2) = -m_one
2079 xx = xx * alpha
2080
2081 ! Density as sum of Gaussians
2082 guess0 = m_zero
2083 do nn = 1, 2
2084 do ip = 1, mesh0%np
2085 call mesh_r(mesh0, ip, rr, origin = xx(:, nn))
2086 guess0(ip) = guess0(ip) + (-1)**nn * beta*exp(-(rr/alpha)**2)
2087 end do
2088 end do
2089
2090 call dio_function_output (io_function_fill_how('AxisX'), ".", "interpolation_target", global_namespace, &
2091 space, mesh0, guess0, unit_one, ierr)
2092 call dio_function_output (io_function_fill_how('AxisZ'), ".", "interpolation_target", global_namespace, &
2093 space, mesh0, guess0, unit_one, ierr)
2094 call dio_function_output (io_function_fill_how('PlaneZ'), ".", "interpolation_target", global_namespace, &
2095 space, mesh0, guess0, unit_one, ierr)
2096
2097 ! We start by testing the interpolation scheme. For this, we generate a function on the fine grid
2098 ! and we inject it on the coarse grid. Then we interpolate it back to the fine grid and we compare
2099 ! This allows for testing the quality of the interpolation scheme
2100
2101 ! move to level 1
2102 call dmultigrid_fine2coarse(mgrid%level(1)%tt, mgrid%level(0)%der, mesh1, guess0, guess1, injection)
2103 ! back to level 0
2104 call dmultigrid_coarse2fine(mgrid%level(1)%tt, mgrid%level(1)%der, mesh0, guess1, res0)
2105
2106 call dio_function_output (io_function_fill_how('AxisX'), ".", "interpolation_result", global_namespace, &
2107 space, mesh0, res0, unit_one, ierr)
2108 call dio_function_output (io_function_fill_how('AxisZ'), ".", "interpolation_result", global_namespace, &
2109 space, mesh0, res0, unit_one, ierr)
2110 call dio_function_output (io_function_fill_how('PlaneZ'), ".", "interpolation_result", global_namespace, &
2111 space, mesh0, res0, unit_one, ierr)
2112
2113 delta = dmf_nrm2(mesh0, guess0(1:mesh0%np)-res0(1:mesh0%np))
2114 write(message(1),'(a,e13.6)') 'Interpolation test (abs.) = ', delta
2115
2116 ! Now we test if the restriction+interpolation combination returns the original result or not
2117 ! This allows to test if restriction and interpolation are adjoint operators or not.
2118
2119 ! move to level 1
2120 call dmultigrid_fine2coarse(mgrid%level(1)%tt, mgrid%level(0)%der, mesh1, guess0, guess1, fullweight)
2121 ! back to level 0
2122 call dmultigrid_coarse2fine(mgrid%level(1)%tt, mgrid%level(1)%der, mesh0, guess1, res0)
2123
2124 call dio_function_output (io_function_fill_how('AxisX'), ".", "restriction_result", global_namespace, &
2125 space, mesh0, res0, unit_one, ierr)
2126 call dio_function_output (io_function_fill_how('AxisZ'), ".", "restriction_result", global_namespace, &
2127 space, mesh0, res0, unit_one, ierr)
2128 call dio_function_output (io_function_fill_how('PlaneZ'), ".", "restriction_result", global_namespace, &
2129 space, mesh0, res0, unit_one, ierr)
2130
2131 delta = dmf_nrm2(mesh0, guess0(1:mesh0%np)-res0(1:mesh0%np))
2132 write(message(2),'(a,e13.6)') 'Restriction test (abs.) = ', delta
2133 call messages_info(2)
2134
2135 safe_deallocate_a(guess0)
2136 safe_deallocate_a(res0)
2137 safe_deallocate_a(guess1)
2138
2140 end subroutine multigrid_test_interpolation
2141
2144 subroutine test_current_density(namespace)
2145 type(namespace_t), intent(in) :: namespace
2146 type(electrons_t), pointer :: sys
2147
2148 type(current_t) :: current
2149 character(len=MAX_PATH_LEN) :: fname
2150 integer :: ierr, ip, idir
2151 integer(int64) :: out_how
2152 real(real64), allocatable :: current_para_ref(:,:), current_dia_ref(:,:), current_mag_ref(:,:), delta(:)
2153 real(real64) :: xx(3), rr, a0, mag_curr, sin_thet, sin_phi, cos_phi, vec_pot_slope
2154 complex(real64) :: alpha
2155
2156 sys => electrons_t(namespace, generate_epot=.false.)
2157 call sys%init_parallelization(mpi_world)
2158
2159 alpha = (0.0_real64, 0.5_real64)
2160 a0 = m_one
2161 vec_pot_slope = 0.4_real64 ! units of B field
2162
2163 call states_elec_allocate_wfns(sys%st, sys%gr, wfs_type = type_cmplx)
2164 call set_hydrogen_states(sys%st%group%psib(1, 1), sys%gr, namespace, alpha, a0)
2165
2166 call current_init(current, namespace)
2167
2168 call hamiltonian_elec_epot_generate(sys%hm, sys%namespace, sys%space, sys%gr, sys%ions, sys%ext_partners, sys%st)
2169
2170 safe_allocate(sys%hm%hm_base%vector_potential(1:3, 1:sys%gr%np))
2171
2172 sys%hm%hm_base%vector_potential = m_zero
2173 do ip = 1, sys%gr%np
2174 xx = sys%gr%x(ip, 1:3)
2175 sys%hm%hm_base%vector_potential(2, ip) = vec_pot_slope * xx(1) / p_c ! vector potential is here devided by c_0
2176 end do
2177
2178 call states_elec_allocate_current(sys%st, sys%space, sys%gr)
2179 call density_calc(sys%st, sys%gr, sys%st%rho)
2180 ! paramagnetic + diamagnetic current densities
2181 call current_calculate(current, namespace, sys%gr, sys%hm, sys%space, sys%st)
2182
2183 safe_allocate(current_para_ref(1:sys%gr%np,1:3))
2184 safe_allocate(current_dia_ref(1:sys%gr%np,1:3))
2185 safe_allocate(current_mag_ref(1:sys%gr%np,1:3))
2186 safe_allocate(delta(1:sys%gr%np))
2187
2188 ! analytic paramagnetic
2189 current_para_ref(:,:) = m_zero
2190 do ip = 1, sys%gr%np
2191 call mesh_r(sys%gr, ip, rr)
2192 xx = sys%gr%x(ip, 1:3)
2193 if (rr > r_small) then
2194 current_para_ref(ip,1:3) = - (psi_2s(rr, a0) * dr_psi_1s(rr, a0) - &
2195 psi_1s(rr, a0) * dr_psi_2s(rr, a0) ) * aimag(alpha) / (1 + abs(alpha)**2) * xx(1:3) / rr
2196 end if
2197 end do
2198
2199 write(fname, '(a)') 'current_para'
2200 out_how = io_function_fill_how("PlaneZ")
2201 call io_function_output_vector(out_how , './' , trim(fname), namespace, sys%space, sys%gr, &
2202 sys%st%current_para(:,:,1), unit_one, ierr)
2203
2204 write(fname, '(a)') 'current_para-ref'
2205 out_how = io_function_fill_how("PlaneZ")
2206 call io_function_output_vector(out_how , './' , trim(fname), namespace, sys%space, sys%gr, &
2207 current_para_ref(:,:), unit_one, ierr)
2208
2209 do idir = 1, 3
2210 delta = m_zero
2211 delta(:) = current_para_ref(1:sys%gr%np, idir) - sys%st%current_para(1:sys%gr%np, idir, 1)
2212 write(message(idir),*) 'idir =',idir,', norm2(delta paramagnetic)',norm2(delta)
2213 end do
2214 call messages_info(3)
2215
2216 ! analytic diamagnetic
2217 current_dia_ref(:,:) = m_zero
2218 do ip = 1, sys%gr%np
2219 call mesh_r(sys%gr, ip, rr)
2220 current_dia_ref(ip,1:3) = - sys%hm%hm_base%vector_potential(1:3,ip) *&
2221 m_one / (1 + abs(alpha)**2) * abs(lc_hydrogen_state(rr, alpha, a0))**2
2222 end do
2223
2224 write(fname, '(a)') 'current_dia'
2225 out_how = io_function_fill_how("PlaneZ")
2226 call io_function_output_vector(out_how , './' , trim(fname), namespace, sys%space, sys%gr, &
2227 sys%st%current_dia(:,:,1), unit_one, ierr)
2228
2229 write(fname, '(a)') 'current_dia-ref'
2230 out_how = io_function_fill_how("PlaneZ")
2231 call io_function_output_vector(out_how , './' , trim(fname), namespace, sys%space, sys%gr, &
2232 current_dia_ref(:,:), unit_one, ierr)
2233
2234 do idir = 1, 3
2235 delta = m_zero
2236 delta(:) = current_dia_ref(1:sys%gr%np, idir) - sys%st%current_dia(1:sys%gr%np, idir, 1)
2237 write(message(idir),*) 'idir =',idir,', norm2(delta diamagnetic)',norm2(delta)
2238 end do
2239 call messages_info(3)
2240
2241 ! magnetization current
2242 call current_calculate_mag(sys%gr%der, sys%st)
2243
2244 ! analytic magnetization
2245 current_mag_ref(:,:) = m_zero
2246 do ip = 1, sys%gr%np
2247 call mesh_r(sys%gr, ip, rr)
2248 xx = sys%gr%x(ip, 1:3)
2249 if (norm2(xx(1:2)) > r_small .and. rr > r_small) then
2250 sin_thet = norm2(xx(1:2)) / rr
2251 sin_phi = xx(2) / norm2(xx(1:2))
2252 cos_phi = xx(1) / norm2(xx(1:2))
2253 mag_curr = m_two * (psi_1s(rr,a0)*dr_psi_1s(rr,a0) + real(alpha)*(psi_1s(rr,a0)*dr_psi_2s(rr,a0) &
2254 + dr_psi_1s(rr,a0)*psi_2s(rr,a0)) + abs(alpha)**2 * psi_2s(rr,a0)*dr_psi_2s(rr,a0))
2255 ! minus signs are reversed because of the charge of the electron
2256 current_mag_ref(ip,1) = m_half * mag_curr * sin_thet * sin_phi / (1+abs(alpha)**2)
2257 current_mag_ref(ip,2) = -m_half * mag_curr * sin_thet * cos_phi / (1+abs(alpha)**2)
2258 endif
2259 end do
2260
2261 write(fname, '(a)') 'current_mag'
2262 out_how = io_function_fill_how("PlaneZ")
2263 call io_function_output_vector(out_how , './' , trim(fname), namespace, sys%space, sys%gr, &
2264 sys%st%current_mag(:,:,1), unit_one, ierr)
2265
2266 write(fname, '(a)') 'current_mag-ref'
2267 out_how = io_function_fill_how("PlaneZ")
2268 call io_function_output_vector(out_how , './' , trim(fname), namespace, sys%space, sys%gr, &
2269 current_mag_ref(:,:), unit_one, ierr)
2270
2271 do idir = 1, 3
2272 delta = m_zero
2273 delta(:) = current_mag_ref(1:sys%gr%np, idir) - sys%st%current_mag(1:sys%gr%np, idir, 1)
2274 write(message(idir),*) 'idir =',idir,', norm2(delta magnetization)',norm2(delta)
2275 end do
2276 call messages_info(3)
2277
2278 safe_deallocate_a(current_para_ref)
2279 safe_deallocate_a(current_dia_ref)
2280 safe_deallocate_a(current_mag_ref)
2281 safe_deallocate_a(delta)
2282 safe_deallocate_p(sys)
2283
2284 end subroutine test_current_density
2285
2286
2287 subroutine set_hydrogen_states(psib, mesh, namespace, alpha, a0)
2288 class(batch_t), intent(inout) :: psib
2289 class(mesh_t), intent(in) :: mesh
2290 type(namespace_t), intent(in) :: namespace
2291 complex(real64), intent(in) :: alpha
2292 real(real64), intent(in) :: a0
2293
2294 complex(real64), allocatable :: zff(:)
2295 integer :: ip
2296 real(real64) :: rr
2297
2298 push_sub(set_hydrogen_states)
2299
2300 ! psi = phi_{1s} + \alpha \phi_{2s}
2301 safe_allocate(zff(1:mesh%np))
2302 if (type_is_complex(psib%type())) then
2303 do ip = 1, mesh%np
2304 call mesh_r(mesh, ip, rr)
2305 zff(ip) = m_one / sqrt(1 + abs(alpha)**2) * lc_hydrogen_state(rr, alpha, a0)
2306 call batch_set_state(psib, 1, mesh%np, zff)
2307 end do
2308 safe_deallocate_a(zff)
2309 else
2310 write(message(1),*) "States should be complex for the linear combination of hydrogenic states to work"
2311 call messages_info(1, namespace=namespace)
2312 end if
2313
2314 pop_sub(set_hydrogen_states)
2315 end subroutine set_hydrogen_states
2316
2317 complex(real64) function lc_hydrogen_state(rr, alpha, a0)
2318 real(real64), intent(in) :: a0, rr
2319 complex(real64), intent(in) :: alpha
2320
2321 lc_hydrogen_state = psi_1s(rr, a0) + alpha * psi_2s(rr, a0)
2322 end function lc_hydrogen_state
2323
2324 ! phi_{1s} = R_{10}(r)Y_{00}(\theta, \phi) = 1/(2\sqrt(\pi)) * 2 a0^{-3/2} * \exp(-r/a0)
2325 real(real64) function psi_1s(rr, a0)
2326 real(real64), intent(in) :: a0, rr
2327
2328 psi_1s = a0**(-m_three/m_two) * exp(-rr / a0) / sqrt(m_pi)
2329 end function psi_1s
2330
2331 ! phi_{2s} = R_{20}(r)Y_{00}(\theta, \phi) = 1/(2\sqrt(\pi)) * \sqrt(2) a0^{-3/2} (1 - r/(2a0)) \exp(-r/(2a0))
2332 real(real64) function psi_2s(rr, a0)
2333 real(real64), intent(in) :: a0, rr
2334
2335 psi_2s = sqrt(m_two) * a0**(-m_three/m_two) * (m_one - rr/(m_two * a0)) * exp(-rr/(m_two * a0)) &
2336 / (m_two * sqrt(m_pi))
2337 end function psi_2s
2338
2339 real(real64) function dr_psi_1s(rr, a0)
2340 real(real64), intent(in) :: a0, rr
2341
2342 dr_psi_1s = -(m_one / a0) * psi_1s(rr, a0)
2343 end function dr_psi_1s
2344
2345 real(real64) function dr_psi_2s(rr, a0)
2346 real(real64), intent(in) :: a0, rr
2347
2348 dr_psi_2s = -(m_half / a0) * psi_2s(rr, a0) - &
2349 a0**(-m_five/m_two) * exp(-rr/(m_two * a0)) / (sqrt(m_two) * m_two * sqrt(m_pi))
2350 end function dr_psi_2s
2351
2352
2353 subroutine test_csv_input(namespace)
2354 type(namespace_t), intent(in) :: namespace
2355
2356 integer :: ierr
2357 integer(int64) :: i, j, k
2358 integer(int64) :: dims(3)
2359 character(len=MAX_PATH_LEN) :: fname
2360
2361 real(real64), allocatable :: ff(:)
2362
2363 push_sub(test_csv_input)
2364
2365 !%Variable TestCSVFileName
2366 !%Type string
2367 !%Section Calculation Modes::Test
2368 !%Description
2369 !% The file name of the csv file to be read
2370 !%End
2371 call parse_variable(namespace, 'TestCSVFileName', "", fname)
2372
2373 message(1) = "Attempting to probe "//trim(fname)
2374 call messages_info(1, namespace=namespace)
2375
2376 call io_csv_get_info(fname, dims, ierr)
2377 if (ierr == 0) then
2378 message(1) = "Probing successful."
2379 write(message(2), '("found dimensions: ",3I20)') dims
2380 call messages_info(2, namespace=namespace)
2381 else
2382 write(message(1),'("Probing failed. ierr = ",I5)') ierr
2383 call messages_fatal(1, namespace=namespace)
2384 end if
2385
2386 safe_allocate(ff(1:dims(1)*dims(2)*dims(3)))
2387
2388 message(1) = "Attempting to read "//trim(fname)
2389 call messages_info(1, namespace=namespace)
2390
2391 call dread_csv(fname, dims(1)*dims(2)*dims(3), ff, ierr)
2392 if (ierr == 0) then
2393 message(1) = "Reading successful."
2394 call messages_info(1, namespace=namespace)
2395
2396 do k=1, min(4_int64, dims(3))
2397 do j=1, min(4_int64, dims(2))
2398 write(message(j), '("data ",2I5, 1X, 4F8.2)') k, j, &
2399 (ff(i + dims(1)*(j-1) + dims(1)* dims(2)*(k-1)), i=1, min(4_int64, dims(1)))
2400 end do
2401 write(message(int(j, int32)), '("")')
2402 call messages_info(int(j, int32), namespace=namespace)
2403 end do
2404
2405 else
2406 message(1) = "Reading failed."
2407 call messages_fatal(1, namespace=namespace)
2408 end if
2409
2410 safe_deallocate_a(ff)
2411
2413 end subroutine test_csv_input
2414
2415end module test_oct_m
2416
2417!! Local Variables:
2418!! mode: f90
2419!! coding: utf-8
2420!! End:
batchified version of the BLAS axpy routine:
Definition: batch_ops.F90:156
scale a batch by a constant or vector
Definition: batch_ops.F90:164
There are several ways how to call batch_set_state and batch_get_state:
Definition: batch_ops.F90:203
constant times a vector plus a vector
Definition: lalg_basic.F90:173
Copies a vector x, to a vector y.
Definition: lalg_basic.F90:188
scales a vector by a constant
Definition: lalg_basic.F90:159
Prints out to iunit a message in the form: ["InputVariable" = value] where "InputVariable" is given b...
Definition: messages.F90:182
double exp(double __x) __attribute__((__nothrow__
double sqrt(double __x) __attribute__((__nothrow__
pure logical function, public accel_is_enabled()
Definition: accel.F90:418
This module implements batches of mesh functions.
Definition: batch.F90:135
This module implements common operations on batches of mesh functions.
Definition: batch_ops.F90:118
subroutine, public dbatch_ax_function_py(np, aa, psi, yy)
This routine performs a set of axpy operations adding the same function psi to all functions of a bat...
Definition: batch_ops.F90:1040
subroutine, public batch_set_zero(this, np, async)
fill all mesh functions of the batch with zero
Definition: batch_ops.F90:244
subroutine, public zbatch_ax_function_py(np, aa, psi, yy)
This routine performs a set of axpy operations adding the same function psi to all functions of a bat...
Definition: batch_ops.F90:2569
Module implementing boundary conditions in Octopus.
Definition: boundaries.F90:124
Module, implementing a factory for boxes.
This module handles the calculation mode.
integer, parameter, public p_strategy_kpoints
parallelization in k-points
integer, parameter, public p_strategy_domains
parallelization in domains
type(calc_mode_par_t), public calc_mode_par
Singleton instance of parallel calculation mode.
integer, parameter, public p_strategy_states
parallelization in states
subroutine, public cgal_polyhedron_init(cgal_poly, fname, verbose)
logical function, public cgal_polyhedron_point_inside(cgal_poly, xq, yq, zq)
subroutine, public cgal_polyhedron_end(cgal_poly)
pure real(real64) function center(this)
Center of the filter interval.
This module implements a calculator for the density and defines related functions.
Definition: density.F90:122
subroutine, public density_calc(st, gr, density, istin)
Computes the density from the orbitals in st.
Definition: density.F90:612
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
subroutine, public dderivatives_test(this, namespace, repetitions, min_blocksize, max_blocksize)
unit test for derivatives
subroutine, public dderivatives_lapl(der, ff, op_ff, ghost_update, set_bc, factor)
apply the Laplacian to a mesh function
subroutine, public zderivatives_test(this, namespace, repetitions, min_blocksize, max_blocksize)
unit test for derivatives
subroutine, public zchebyshev_filter(namespace, mesh, st, hm, degree, bounds, ik, normalize)
Chebyshev Filter.
subroutine, public dchebyshev_filter(namespace, mesh, st, hm, degree, bounds, ik, normalize)
Chebyshev Filter.
integer, parameter, public spin_polarized
subroutine, public exponential_init(te, namespace, full_batch)
real(real64), parameter, public m_two
Definition: global.F90:192
real(real64), parameter, public m_zero
Definition: global.F90:190
real(real64), parameter, public m_four
Definition: global.F90:194
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:188
complex(real64), parameter, public m_zi
Definition: global.F90:204
real(real64), parameter, public m_epsilon
Definition: global.F90:206
real(real64), parameter, public m_half
Definition: global.F90:196
real(real64), parameter, public m_one
Definition: global.F90:191
This module implements the underlying real-space grid.
Definition: grid.F90:119
subroutine, public zhamiltonian_elec_apply_batch(hm, namespace, mesh, psib, hpsib, terms, set_bc)
subroutine, public hamiltonian_elec_epot_generate(this, namespace, space, gr, ions, ext_partners, st, time)
subroutine, public dhamiltonian_elec_apply_batch(hm, namespace, mesh, psib, hpsib, terms, set_bc)
type(hardware_t), public cpu_hardware
Global instance of CPU hardware specification.
Definition: hardware.F90:61
integer, parameter, public sizeof_real64
Number of bytes to store a variable of type real(real64)
Definition: hardware.F90:55
integer, parameter, public sizeof_complex64
Number of bytes to store a variable of type complex(real64)
Definition: hardware.F90:58
The Helmholtz decomposition is intended to contain "only mathematical" functions and procedures to co...
Test suit for the Helmholtz decomposition module.
subroutine, public gaussian_test(helmholtz, sys_grid, namespace, space)
subroutine, public hertzian_dipole_test(helmholtz, sys_grid, namespace, space)
This module implements a simple hash table for non-negative integer keys and integer values.
Definition: iihash.F90:127
subroutine, public iihash_end(h)
Free a hash table.
Definition: iihash.F90:186
subroutine, public iihash_insert(h, key, val)
Insert a (key, val) pair into the hash table h.
Definition: iihash.F90:208
integer function, public iihash_lookup(h, key, found)
Look up a value in the hash table h. If found is present, it indicates if key could be found in the t...
Definition: iihash.F90:233
subroutine, public iihash_init(h)
Initialize a hash table h.
Definition: iihash.F90:163
subroutine, public dio_function_output(how, dir, fname, namespace, space, mesh, ff, unit, ierr, pos, atoms, grp, root)
Top-level IO routine for functions defined on the mesh.
integer(int64) function, public io_function_fill_how(where)
Use this function to quickly plot functions for debugging purposes: call dio_function_output(io_funct...
Definition: io.F90:116
subroutine, public ion_interaction_test(space, latt, atom, natoms, pos, lsize, namespace, mc)
This modules takes care of testing some linear algebra routines.
subroutine, public test_exponential_matrix(namespace)
Unit tests for the exponential of a matrix.
This module defines functions over batches of mesh functions.
Definition: mesh_batch.F90:118
subroutine, public dmesh_batch_dotp_matrix(mesh, aa, bb, dot, reduce)
Calculate the overlap matrix of two batches.
Definition: mesh_batch.F90:273
subroutine, public dmesh_batch_dotp_self(mesh, aa, dot, reduce)
calculate the overlap matrix of a batch with itself
Definition: mesh_batch.F90:535
subroutine, public mesh_batch_nrm2(mesh, aa, nrm2, reduce)
Calculate the norms (norm2, not the square!) of a batch of mesh functions.
Definition: mesh_batch.F90:179
subroutine, public zmesh_batch_dotp_vector(mesh, aa, bb, dot, reduce, cproduct)
calculate the vector of dot-products of mesh functions between two batches
subroutine, public zmesh_batch_dotp_matrix(mesh, aa, bb, dot, reduce)
Calculate the overlap matrix of two batches.
subroutine, public zmesh_batch_normalize(mesh, psib, norm)
Normalize a batch.
subroutine, public zmesh_batch_dotp_self(mesh, aa, dot, reduce)
calculate the overlap matrix of a batch with itself
subroutine, public dmesh_batch_normalize(mesh, psib, norm)
Normalize a batch.
subroutine, public dmesh_batch_dotp_vector(mesh, aa, bb, dot, reduce, cproduct)
calculate the vector of dot-products of mesh functions between two batches
Definition: mesh_batch.F90:649
This module defines various routines, operating on mesh functions.
class(mesh_t), pointer, public mesh_aux
Globally-scoped pointer to the mesh instance.
real(real64) function, public dmf_dotp_aux(f1, f2)
dot product between two vectors (mesh functions)
subroutine, public mesh_init_mesh_aux(mesh)
Initialise a pointer to the grid/mesh, that is globally exposed, such that low level mesh operations ...
subroutine, public zmesh_interpolation_test(mesh)
subroutine, public dmesh_interpolation_test(mesh)
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:120
pure subroutine, public mesh_r(mesh, ip, rr, origin, coords)
return the distance to the origin for a given grid point
Definition: mesh.F90:339
subroutine, public messages_print_with_emphasis(msg, iunit, namespace)
Definition: messages.F90:904
character(len=512), private msg
Definition: messages.F90:167
subroutine, public messages_obsolete_variable(namespace, name, rep)
Definition: messages.F90:1029
subroutine, public messages_new_line()
Definition: messages.F90:1118
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:162
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:600
This modules takes care of testing optimizers using standard test functions.
subroutine, public test_optimizers(namespace)
Unit tests for different optimizers.
This module implements unit tests for the mixing methods.
Definition: mix_tests.F90:117
subroutine, public mix_tests_run()
Definition: mix_tests.F90:139
type(mpi_grp_t), public mpi_world
Definition: mpi.F90:272
subroutine, public test_mpiwrappers
Definition: mpi_test.F90:130
This module handles the communicators for the various parallelization strategies.
Definition: multicomm.F90:147
subroutine, public multigrid_end(mgrid)
Definition: multigrid.F90:516
subroutine, public multigrid_init(mgrid, namespace, space, mesh, der, stencil, mc, nlevels)
Definition: multigrid.F90:187
type(namespace_t), public global_namespace
Definition: namespace.F90:134
subroutine, public orbitalbasis_end(this)
subroutine, public zorbitalbasis_build(this, namespace, ions, mesh, kpt, ndim, skip_s_orb, use_all_orb, verbose)
This routine is an interface for constructing the orbital basis.
subroutine, public dorbitalbasis_build(this, namespace, ions, mesh, kpt, ndim, skip_s_orb, use_all_orb, verbose)
This routine is an interface for constructing the orbital basis.
subroutine, public orbitalbasis_init(this, namespace, periodic_dim)
subroutine, public dorbitalset_add_to_batch(os, ndim, psib, weight)
Definition: orbitalset.F90:843
subroutine, public zorbitalset_add_to_batch(os, ndim, psib, weight)
subroutine, public orbitalset_update_phase(os, dim, kpt, kpoints, spin_polarized, vec_pot, vec_pot_var, kpt_max)
Build the phase correction to the global phase in case the orbital crosses the border of the simulato...
Definition: orbitalset.F90:281
subroutine, public dorbitalset_get_coeff_batch(os, ndim, psib, dot, reduce)
Definition: orbitalset.F90:597
subroutine, public zorbitalset_get_coeff_batch(os, ndim, psib, dot, reduce)
subroutine, public poisson_test(this, space, mesh, latt, namespace, repetitions)
This routine checks the Hartree solver selected in the input file by calculating numerically and anal...
subroutine, public preconditioner_end(this)
subroutine, public preconditioner_init(this, namespace, gr, mc, space)
subroutine, public zproject_psi_batch(mesh, bnd, pj, npj, dim, psib, ppsib)
To optimize the application of the non-local operator in parallel, the projectors are applied in step...
Definition: projector.F90:1256
Implementation details for regridding.
Definition: regridding.F90:172
This module implements a simple hash table for string valued keys and integer values using the C++ ST...
Definition: sihash.F90:120
subroutine, public sihash_insert(h, key, val)
Insert a (key, val) pair into the hash table h.
Definition: sihash.F90:205
subroutine, public sihash_init(h)
Initialize a hash table h with size entries. Since we use separate chaining, the number of entries in...
Definition: sihash.F90:164
integer function, public sihash_lookup(h, key, found)
Look up a value in the hash table h. If found is present, it indicates if key could be found in the t...
Definition: sihash.F90:231
subroutine, public sihash_end(h)
Free a hash table.
Definition: sihash.F90:185
This module is intended to contain "only mathematical" functions and procedures.
Definition: solvers.F90:117
This module implements a simple hash table for string valued keys and integer values using the C++ ST...
Definition: sphash.F90:120
subroutine, public sphash_init(h)
Initialize a hash table h with size entries. Since we use separate chaining, the number of entries in...
Definition: sphash.F90:224
subroutine, public sphash_insert(h, key, val, clone)
Insert a (key, val) pair into the hash table h. If clone=.true., the object will be copied.
Definition: sphash.F90:293
subroutine, public sphash_end(h)
Free a hash table.
Definition: sphash.F90:250
class(*) function, pointer, public sphash_lookup(h, key, found)
Look up a value in the hash table h. If found is present, it indicates if key could be found in the t...
Definition: sphash.F90:325
pure logical function, public states_are_real(st)
subroutine, public zstates_elec_calc_orth_test(st, namespace, mesh, kpoints)
subroutine, public dstates_elec_calc_orth_test(st, namespace, mesh, kpoints)
This module handles spin dimensions of the states and the k-point distribution.
subroutine, public states_elec_deallocate_wfns(st)
Deallocates the KS wavefunctions defined within a states_elec_t structure.
subroutine, public states_elec_allocate_wfns(st, mesh, wfs_type, skip, packed)
Allocates the KS wavefunctions defined within a states_elec_t structure.
subroutine, public dsubspace_diag(this, namespace, mesh, st, hm, ik, eigenval, diff, nonortho)
Diagonalises the Hamiltonian in the subspace defined by the states.
Definition: subspace.F90:299
subroutine, public zsubspace_diag(this, namespace, mesh, st, hm, ik, eigenval, diff, nonortho)
Diagonalises the Hamiltonian in the subspace defined by the states.
Definition: subspace.F90:767
Integration tests for ISDF.
Definition: test_isdf.F90:116
subroutine, public test_isdf(namespace, serial)
Set up an electron system, compute some optimal centroid positions, and use these to build a set of I...
Definition: test_isdf.F90:158
subroutine, public test_weighted_kmeans(namespace)
Test weighted kmeans algorithm for a finite system.
This module implements a unit-test like runmode for Octopus.
Definition: test.F90:117
subroutine set_der_aux(der)
Definition: test.F90:503
real(real64) function psi_2s(rr, a0)
Definition: test.F90:2428
subroutine test_sihash()
Definition: test.F90:1785
subroutine test_cgal()
Definition: test.F90:1607
subroutine test_helmholtz_decomposition(namespace)
Definition: test.F90:410
subroutine test_hartree(param, namespace)
Definition: test.F90:391
subroutine test_derivatives(param, namespace)
Definition: test.F90:1373
subroutine laplacian_op(x, hx, userdata)
Computes .
Definition: test.F90:512
subroutine test_subspace_diagonalization(param, namespace)
Definition: test.F90:994
subroutine, public test_run(namespace)
Components and integration test runner.
Definition: test.F90:213
subroutine test_density_calc(param, namespace)
Definition: test.F90:774
subroutine test_iihash()
Definition: test.F90:1759
real(real64) function psi_1s(rr, a0)
Definition: test.F90:2421
subroutine test_current_density(namespace)
Here we test the different contributions to the total electronic current density.
Definition: test.F90:2240
subroutine test_batch_set_gaussian(psib, mesh)
Definition: test.F90:1686
subroutine test_regridding(namespace)
Definition: test.F90:1924
subroutine test_sphash(namespace)
Definition: test.F90:1831
subroutine test_hamiltonian(param, namespace)
Definition: test.F90:692
subroutine test_dense_eigensolver()
Definition: test.F90:1627
subroutine test_interpolation(param, namespace)
Definition: test.F90:1445
subroutine multigrid_test_interpolation(mgrid, space)
Definition: test.F90:2143
subroutine test_ion_interaction(namespace)
Definition: test.F90:1483
subroutine test_boundaries(param, namespace)
Definition: test.F90:812
subroutine test_orthogonalization(param, namespace)
Definition: test.F90:1403
real(real64) function dr_psi_2s(rr, a0)
Definition: test.F90:2441
subroutine test_batch_ops(param, namespace)
Definition: test.F90:1047
complex(real64) function lc_hydrogen_state(rr, alpha, a0)
Definition: test.F90:2413
subroutine test_projector(param, namespace)
Definition: test.F90:532
subroutine test_dft_u(param, namespace)
Definition: test.F90:584
subroutine test_prints_info_batch(st, gr, psib, string)
Definition: test.F90:1503
subroutine test_composition_chebyshev(namespace)
Test the composition rule for Chebyshev polynomials.
Definition: test.F90:851
real(real64) function dr_psi_1s(rr, a0)
Definition: test.F90:2435
subroutine test_linear_solver(namespace)
Definition: test.F90:441
subroutine test_grid_interpolation()
Definition: test.F90:1738
subroutine set_hydrogen_states(psib, mesh, namespace, alpha, a0)
Definition: test.F90:2383
subroutine test_exponential(param, namespace)
Definition: test.F90:945
subroutine test_csv_input(namespace)
Definition: test.F90:2449
subroutine test_vecpot_analytical(namespace)
Here, analytical formulation for vector potential and B field are used. Ref: Sangita Sen and Erik I....
Definition: test.F90:2042
subroutine test_clock()
Definition: test.F90:1548
type(type_t), public type_cmplx
Definition: types.F90:136
logical pure function, public type_is_complex(this)
Definition: types.F90:182
This module defines the unit system, used for input and output.
type(unit_t), public unit_one
some special units required for particular quantities
subroutine, public v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval, time, calc_energy, calc_current, force_semilocal)
Definition: v_ks.F90:749
Class defining batches of mesh functions.
Definition: batch.F90:161
class representing derivatives
Class describing the electron system.
Definition: electrons.F90:220
Description of the grid, containing information on derivatives, stencil, and symmetries.
Definition: grid.F90:171
This class implements the iteration counter used by the multisystem algorithms. As any iteration coun...
Describes mesh distribution to nodes.
Definition: mesh.F90:187
contains the information of the meshes and provides the transfer functions
Definition: regridding.F90:198
The states_elec_t class contains all electronic wave functions.
batches of electronic states
Definition: wfs_elec.F90:141
int true(void)
subroutine write_clock(operation)
Definition: test.F90:1588
subroutine write_condition_result(condition, result)
Definition: test.F90:1595
real(real64) function values(xx)
Definition: test.F90:2025