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