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