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