27 use,
intrinsic :: iso_fortran_env
81 real(real64),
allocatable :: ddf(:, :, :)
82 real(real64),
allocatable :: ddv(:, :, :)
83 real(real64),
allocatable :: df_old(:, :)
84 real(real64),
allocatable :: dvin_old(:, :)
85 real(real64),
allocatable :: dvin(:, :)
86 real(real64),
allocatable :: dvout(:, :)
87 real(real64),
allocatable :: dvnew(:, :)
88 real(real64),
allocatable :: dresidual(:, :)
91 complex(real64),
allocatable :: zdf(:, :, :)
92 complex(real64),
allocatable :: zdv(:, :, :)
93 complex(real64),
allocatable :: zf_old(:, :)
94 complex(real64),
allocatable :: zvin_old(:, :)
95 complex(real64),
allocatable :: zvin(:, :)
96 complex(real64),
allocatable :: zvout(:, :)
97 complex(real64),
allocatable :: zvnew(:, :)
98 complex(real64),
allocatable :: zresidual(:, :)
100 type(type_t) :: func_type
101 integer :: d1, d2, d3
108 type(mixfield_t),
pointer :: p
111 integer,
parameter :: MAX_AUXMIXFIELD = 5
122 type(derivatives_t),
pointer :: der
125 real(real64) :: coeff
126 real(real64) :: residual_coeff
129 integer,
public :: ns
130 integer,
public :: ns_restart
137 type(mixfield_t) :: mixfield
138 integer :: nauxmixfield
139 type(mixfield_ptr_t) :: auxmixfield(MAX_AUXMIXFIELD)
145 real(real64) :: kerker_factor
146 logical :: precondition
147 type(nl_operator_t) :: preconditioner
151 procedure :: compute_residuals_aux_field
171 subroutine mix_init(smix, namespace, space, der, d1, d2, def_, func_type_, prefix_)
172 type(mix_t),
intent(out) :: smix
173 type(namespace_t),
intent(in) :: namespace
176 integer,
intent(in) :: d1
177 integer,
intent(in) :: d2
178 integer,
optional,
intent(in) :: def_
179 type(
type_t),
optional,
intent(in) :: func_type_
180 character(len=*),
optional,
intent(in) :: prefix_
183 character(len=32) :: prefix
190 def = option__mixingscheme__broyden
191 if (
present(def_)) def = def_
192 if (
present(func_type_))
then
193 func_type = func_type_
198 if (
present(prefix_)) prefix = prefix_
239 call parse_variable(namespace, trim(prefix)+
'MixingPreconditioner', .false., smix%precondition)
240 if (der%dim /= 3) smix%precondition = .false.
251 call parse_variable(namespace, trim(prefix)+
'Mixing', 0.3_real64, smix%coeff)
264 call parse_variable(namespace, trim(prefix)+
'MixingResidual', 0.05_real64, smix%residual_coeff)
265 if (smix%residual_coeff <=
m_zero .or. smix%residual_coeff >
m_one)
then
266 call messages_input_error(namespace,
'MixingResidual',
'Value should be positive and smaller than one.')
278 if (smix%scheme /= option__mixingscheme__linear)
then
279 call parse_variable(namespace, trim(prefix)//
'MixNumberSteps', 4, smix%ns)
294 if (smix%scheme /= option__mixingscheme__linear)
then
295 call parse_variable(namespace, trim(prefix)//
'MixingRestart', 20, smix%ns_restart)
301 write(
message(1),
'(A,I4,A,I4,A)')
"Info: Mixing uses ", smix%ns,
" steps and restarts after ", &
302 smix%ns_restart,
" steps."
315 call parse_variable(namespace, trim(prefix)//
'MixInterval', 1, smix%interval)
316 if (smix%interval < 1)
call messages_input_error(namespace,
'MixInterval',
'MixInterval must be larger or equal than 1')
331 call parse_variable(namespace, trim(prefix)+
'MixingKerker', .false., smix%kerker)
344 call parse_variable(namespace, trim(prefix)+
'MixingKerkerFactor', 1._real64, smix%kerker_factor)
346 smix%nauxmixfield = 0
347 do ii = 1,max_auxmixfield
348 nullify(smix%auxmixfield(ii)%p)
351 call mixfield_init(smix, smix%mixfield, d1, d2, smix%ns, func_type)
361 integer :: ns, maxp, ip, is
362 real(real64),
parameter :: weight = 50.0_real64
369 assert(.not. der%mesh%use_curvilinear)
374 call nl_operator_build(space, der%mesh, smix%preconditioner, der%mesh%np, const_w = .not. der%mesh%use_curvilinear)
376 ns = smix%preconditioner%stencil%size
378 if (smix%preconditioner%const_w)
then
387 select case (sum(abs(smix%preconditioner%stencil%points(1:der%dim, is))))
389 smix%preconditioner%w(is, ip) =
m_one + weight/8.0_real64
391 smix%preconditioner%w(is, ip) = weight/16.0_real64
393 smix%preconditioner%w(is, ip) = weight/32.0_real64
395 smix%preconditioner%w(is, ip) = weight/64.0_real64
414 type(
mix_t),
intent(inout) :: smix
429 type(
mix_t),
intent(inout) :: smix
439 smix%nauxmixfield = 0
440 do ii = 1,max_auxmixfield
441 nullify(smix%auxmixfield(ii)%p)
449 subroutine mix_dump(namespace, restart, smix, space, mesh, ierr)
453 class(
space_t),
intent(in) :: space
454 class(
mesh_t),
intent(in) :: mesh
455 integer,
intent(out) :: ierr
457 integer :: iunit, id2, id3, err, err2(4)
458 character(len=40) :: lines(7)
459 character(len=80) :: filename
470 message(1) =
"Debug: Writing mixing restart."
474 assert(mesh%np == smix%mixfield%d1)
478 write(lines(1),
'(a11,i1)')
'scheme= ', smix%scheme
480 write(lines(2),
'(a11,i10)')
'd1= ', mesh%np_global
481 write(lines(3),
'(a11,i10)')
'd2= ', smix%mixfield%d2
482 write(lines(4),
'(a11,i10)')
'd3= ', smix%mixfield%d3
483 write(lines(5),
'(a11,i10)')
'iter= ', smix%iter
484 write(lines(6),
'(a11,i10)')
'ns= ', smix%ns
485 write(lines(7),
'(a11,i10)')
'last_ipos= ', smix%last_ipos
493 if (smix%scheme /= option__mixingscheme__linear)
then
494 do id2 = 1, smix%mixfield%d2
495 do id3 = 1, smix%mixfield%d3
497 write(filename,
'(a3,i2.2,i2.2)')
'df_', id2, id3
498 if (smix%mixfield%func_type ==
type_float)
then
503 err2(1) = err2(1) + err
505 write(filename,
'(a3,i2.2,i2.2)')
'dv_', id2, id3
511 err2(2) = err2(2) + err
515 write(filename,
'(a6,i2.2)')
'f_old_', id2
516 if (smix%mixfield%func_type ==
type_float)
then
521 err2(3) = err2(3) + err
523 write(filename,
'(a8,i2.2)')
'vin_old_', id2
524 if (smix%mixfield%func_type ==
type_float)
then
529 err2(4) = err2(4) + err
533 if (err2(1) /= 0) ierr = ierr + 2
534 if (err2(2) /= 0) ierr = ierr + 4
535 if (err2(3) /= 0) ierr = ierr + 8
536 if (err2(4) /= 0) ierr = ierr + 16
539 message(1) =
"Debug: Writing mixing restart done."
547 subroutine mix_load(namespace, restart, smix, space, mesh, ierr)
550 type(
mix_t),
intent(inout) :: smix
551 class(
space_t),
intent(in) :: space
552 class(
mesh_t),
intent(in) :: mesh
553 integer,
intent(out) :: ierr
555 integer :: iunit, err, err2(4)
556 integer :: scheme, d1, d2, d3, ns
558 character(len=11) :: str
559 character(len=80) :: filename
560 character(len=256) :: lines(7)
572 message(1) =
"Debug: Reading mixing restart."
581 read(lines(1), *) str, scheme
582 read(lines(2), *) str, d1
583 read(lines(3), *) str, d2
584 read(lines(4), *) str, d3
585 read(lines(5), *) str, smix%iter
586 read(lines(6), *) str, ns
587 read(lines(7), *) str, smix%last_ipos
594 if (scheme /= smix%scheme .or. ns /= smix%ns)
then
595 message(1) =
"The mixing scheme from the restart data is not the same as the one used in the current calculation."
601 if (mesh%np_global /= d1 .or. mesh%np /= smix%mixfield%d1 .or. d2 /= smix%mixfield%d2)
then
602 message(1) =
"The dimensions of the arrays from the mixing restart data"
603 message(2) =
"are not the same as the ones used in this calculation."
613 if (smix%scheme /= option__mixingscheme__linear)
then
615 do id2 = 1, smix%mixfield%d2
616 do id3 = 1, smix%mixfield%d3
618 write(filename,
'(a3,i2.2,i2.2)')
'df_', id2, id3
619 if (smix%mixfield%func_type ==
type_float)
then
624 if (err /= 0) err2(1) = err2(1) + 1
626 write(filename,
'(a3,i2.2,i2.2)')
'dv_', id2, id3
627 if (smix%mixfield%func_type ==
type_float)
then
632 if (err /= 0) err2(2) = err2(2) + 1
636 write(filename,
'(a6,i2.2)')
'f_old_', id2
637 if (smix%mixfield%func_type ==
type_float)
then
642 if (err /= 0) err2(3) = err2(3) + 1
644 write(filename,
'(a8,i2.2)')
'vin_old_', id2
645 if (smix%mixfield%func_type ==
type_float)
then
650 if (err /= 0) err2(4) = err2(4) + 1
654 if (err2(1) /= 0) ierr = ierr + 8
655 if (err2(2) /= 0) ierr = ierr + 16
656 if (err2(3) /= 0) ierr = ierr + 32
657 if (err2(4) /= 0) ierr = ierr + 64
666 message(1) =
"Debug: Reading mixing restart done."
673 type(
mix_t),
intent(in) :: this
675 coefficient = this%coeff
678 integer pure function
mix_scheme(this) result(scheme)
679 type(
mix_t),
intent(in) :: this
684 integer pure function
mix_d3(this)
685 type(
mix_t),
intent(in) :: this
691 type(
mix_t),
target,
intent(in) :: this
692 type(
mixfield_t),
pointer,
intent(out) :: mixfield
694 mixfield => this%mixfield
698 subroutine mixing(namespace, smix)
699 type(namespace_t),
intent(in) :: namespace
700 type(
mix_t),
intent(inout) :: smix
704 if (smix%mixfield%func_type == type_float)
then
705 call dmixing(namespace, smix, smix%mixfield%dvin, smix%mixfield%dvout, smix%mixfield%dvnew)
707 call zmixing(namespace, smix, smix%mixfield%zvin, smix%mixfield%zvout, smix%mixfield%zvnew)
714 type(namespace_t),
intent(in) :: namespace
715 type(
mix_t),
intent(inout) :: smix
716 type(
mixfield_t),
target,
intent(in) :: mixfield
720 smix%nauxmixfield = smix%nauxmixfield + 1
721 smix%auxmixfield(smix%nauxmixfield)%p => mixfield
723 if (smix%scheme == option__mixingscheme__diis)
then
724 message(1) =
'Mixing scheme DIIS is not implemented for auxiliary mixing fields'
725 call messages_fatal(1, namespace=namespace)
732 subroutine mixfield_init(smix, mixfield, d1, d2, d3, func_type)
733 type(
mix_t),
intent(in) :: smix
735 integer,
intent(in) :: d1, d2, d3
736 type(type_t),
intent(in) :: func_type
744 mixfield%func_type = func_type
746 if (smix%scheme /= option__mixingscheme__linear)
then
747 if (mixfield%func_type == type_float)
then
748 safe_allocate( mixfield%ddf(1:d1, 1:d2, 1:d3))
749 safe_allocate( mixfield%ddv(1:d1, 1:d2, 1:d3))
750 safe_allocate(mixfield%dvin_old(1:d1, 1:d2))
751 safe_allocate( mixfield%df_old(1:d1, 1:d2))
753 safe_allocate( mixfield%zdf(1:d1, 1:d2, 1:d3))
754 safe_allocate( mixfield%zdv(1:d1, 1:d2, 1:d3))
755 safe_allocate(mixfield%zvin_old(1:d1, 1:d2))
756 safe_allocate( mixfield%zf_old(1:d1, 1:d2))
760 if (mixfield%func_type == type_float)
then
761 safe_allocate(mixfield%dvin(1:d1, 1:d2))
762 safe_allocate(mixfield%dvout(1:d1, 1:d2))
763 safe_allocate(mixfield%dvnew(1:d1, 1:d2))
764 safe_allocate(mixfield%dresidual(1:d1, 1:d2))
766 safe_allocate(mixfield%zvin(1:d1, 1:d2))
767 safe_allocate(mixfield%zvout(1:d1, 1:d2))
768 safe_allocate(mixfield%zvnew(1:d1, 1:d2))
769 safe_allocate(mixfield%zresidual(1:d1, 1:d2))
783 if (smix%scheme /= option__mixingscheme__linear)
then
784 if (mixfield%func_type == type_float)
then
785 safe_deallocate_a(mixfield%ddf)
786 safe_deallocate_a(mixfield%ddv)
787 safe_deallocate_a(mixfield%dvin_old)
788 safe_deallocate_a(mixfield%df_old)
790 safe_deallocate_a(mixfield%zdf)
791 safe_deallocate_a(mixfield%zdv)
792 safe_deallocate_a(mixfield%zvin_old)
793 safe_deallocate_a(mixfield%zf_old)
797 if (mixfield%func_type == type_float)
then
798 safe_deallocate_a(mixfield%dvin)
799 safe_deallocate_a(mixfield%dvout)
800 safe_deallocate_a(mixfield%dvnew)
801 safe_deallocate_a(mixfield%dresidual)
803 safe_deallocate_a(mixfield%zvin)
804 safe_deallocate_a(mixfield%zvout)
805 safe_deallocate_a(mixfield%zvnew)
806 safe_deallocate_a(mixfield%zresidual)
814 integer,
intent(in) :: scheme
816 integer :: d1, d2, d3
824 if (scheme /= option__mixingscheme__linear)
then
825 if (mixfield%func_type == type_float)
then
826 assert(
allocated(mixfield%ddf))
827 mixfield%ddf(1:d1, 1:d2, 1:d3) = m_zero
828 mixfield%ddv(1:d1, 1:d2, 1:d3) = m_zero
829 mixfield%dvin_old(1:d1, 1:d2) = m_zero
830 mixfield%df_old(1:d1, 1:d2) = m_zero
832 assert(
allocated(mixfield%zdf))
833 mixfield%zdf(1:d1, 1:d2, 1:d3) = m_z0
834 mixfield%zdv(1:d1, 1:d2, 1:d3) = m_z0
835 mixfield%zvin_old(1:d1, 1:d2) = m_z0
836 mixfield%zf_old(1:d1, 1:d2) = m_z0
840 if (mixfield%func_type == type_float)
then
841 mixfield%dvin(1:d1, 1:d2) = m_zero
842 mixfield%dvout(1:d1, 1:d2) = m_zero
843 mixfield%dvnew(1:d1, 1:d2) = m_zero
844 mixfield%dresidual(1:d1, 1:d2) = m_zero
846 mixfield%zvin(1:d1, 1:d2) = m_z0
847 mixfield%zvout(1:d1, 1:d2) = m_z0
848 mixfield%zvnew(1:d1, 1:d2) = m_z0
849 mixfield%zresidual(1:d1, 1:d2) = m_z0
861 real(real64),
intent(in) :: vin_re(:,:), vin_im(:,:)
865 mixfield%zvin(1:mixfield%d1, 1:mixfield%d2) = vin_re(1:mixfield%d1, 1:mixfield%d2) &
866 + m_zi * vin_im(1:mixfield%d1, 1:mixfield%d2)
873 real(real64),
intent(in) :: vout_re(:,:), vout_im(:,:)
877 mixfield%zvout(1:mixfield%d1, 1:mixfield%d2) = vout_re(1:mixfield%d1, 1:mixfield%d2) &
878 + m_zi * vout_im(1:mixfield%d1, 1:mixfield%d2)
885 real(real64),
intent(inout) :: re(:,:), im(:,:)
887 push_sub(mixfield_get_ddvnew)
889 re(1:mixfield%d1, 1:mixfield%d2) = real(mixfield%zvnew(1:mixfield%d1, 1:mixfield%d2), real64)
890 im(1:mixfield%d1, 1:mixfield%d2) = aimag(mixfield%zvnew(1:mixfield%d1, 1:mixfield%d2))
892 pop_sub(mixfield_get_ddvnew)
904 class(
mix_t),
intent(inout) :: this
912 do i = 1, this%nauxmixfield
913 aux_field => this%auxmixfield(i)%p
914 dims = [aux_field%d1, aux_field%d2]
916 if (aux_field%func_type == type_float)
then
917 call dcompute_residual(this, dims, aux_field%dvin, aux_field%dvout, aux_field%dresidual)
919 call zcompute_residual(this, dims, aux_field%zvin, aux_field%zvout, aux_field%zresidual)
941 type(
mix_t),
intent(inout) :: smix
948 call smix%compute_residuals_aux_field()
950 do i = 1, smix%nauxmixfield
951 field => smix%auxmixfield(i)%p
953 if (field%func_type == type_float)
then
954 call dmixing_linear(field%d1, field%d2, smix%coeff, field%dvin, field%dresidual, field%dvnew)
956 call zmixing_linear(field%d1, field%d2, smix%coeff, field%zvin, field%zresidual, field%zvnew)
969#include "mix_inc.F90"
972#include "complex.F90"
974#include "mix_inc.F90"
subroutine init_preconditioner()
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
real(real64), parameter, public m_zero
real(real64), parameter, public m_one
This module is intended to contain "only mathematical" functions and procedures.
This module defines various routines, operating on mesh functions.
This module defines the meshes, which are used in Octopus.
subroutine, public messages_warning(no_lines, all_nodes, namespace)
subroutine, public messages_obsolete_variable(namespace, name, rep)
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
subroutine, public messages_input_error(namespace, var, details, row, column)
subroutine, public messages_experimental(name, namespace)
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
subroutine, public mixfield_end(smix, mixfield)
Deallocate all arrays of a mixfield instance.
integer pure function, public mix_scheme(this)
subroutine dmixfield_set_vin(mixfield, vin)
real(real64) pure function, public mix_coefficient(this)
subroutine, public mixing(namespace, smix)
Main entry-point to SCF mixer.
subroutine, public mixfield_init(smix, mixfield, d1, d2, d3, func_type)
Initialise all attributes of a mixfield instance.
subroutine ddmixfield_get_vnew(mixfield, re, im)
subroutine zmixfield_get_vnew(mixfield, vnew)
subroutine dcompute_residual(smix, dims, vin, vout, residual)
Compute the residual of the potential (or density) for SCF mixing.
subroutine, public mix_get_field(this, mixfield)
subroutine compute_residuals_aux_field(this)
Compute the residuals of the auxilliary fields.
subroutine zmixfield_set_vout(mixfield, vout)
subroutine, public mix_add_auxmixfield(namespace, smix, mixfield)
subroutine dmixfield_get_vnew(mixfield, vnew)
subroutine dmixing_linear(d1, d2, coeff, vin, residual, vnew)
Linear mixing of the input and output potentials.
subroutine, public mix_dump(namespace, restart, smix, space, mesh, ierr)
subroutine, public mixfield_clear(scheme, mixfield)
Zero all potential and field attributes of a mixfield instance.
subroutine, public mix_load(namespace, restart, smix, space, mesh, ierr)
subroutine, public mix_init(smix, namespace, space, der, d1, d2, def_, func_type_, prefix_)
Initialise mix_t instance.
subroutine ddmixfield_set_vin(mixfield, vin_re, vin_im)
subroutine zmixfield_set_vin(mixfield, vin)
subroutine, public zmixing(namespace, smix, vin, vout, vnew)
Mix the input and output potentials (or densities) according to some scheme.
integer pure function, public mix_d3(this)
subroutine dmixfield_set_vout(mixfield, vout)
subroutine zcompute_residual(smix, dims, vin, vout, residual)
Compute the residual of the potential (or density) for SCF mixing.
subroutine ddmixfield_set_vout(mixfield, vout_re, vout_im)
subroutine, public dmixing(namespace, smix, vin, vout, vnew)
Mix the input and output potentials (or densities) according to some scheme.
subroutine linear_mixing_aux_field(smix)
Linear mixing of the auxilliary fields.
subroutine zmixing_linear(d1, d2, coeff, vin, residual, vnew)
Linear mixing of the input and output potentials.
subroutine, public mix_end(smix)
subroutine, public mix_clear(smix)
This module defines non-local operators.
subroutine, public nl_operator_init(op, label)
initialize an instance of a non-local operator by setting the label
subroutine, public nl_operator_update_gpu_buffers(op)
subroutine, public nl_operator_output_weights(this)
subroutine, public nl_operator_end(op)
subroutine, public nl_operator_build(space, mesh, op, np, const_w, regenerate)
subroutine, public nl_operator_allocate_gpu_buffers(op)
subroutine, public restart_read(restart, iunit, lines, nlines, ierr)
subroutine, public restart_close(restart, iunit)
Close a file previously opened with restart_open.
subroutine, public restart_write(restart, iunit, lines, nlines, ierr)
logical pure function, public restart_skip(restart)
Returns true if the restart information should neither be read nor written. This might happen because...
integer function, public restart_open(restart, filename, status, position, silent)
Open file "filename" found inside the current restart directory. Depending on the type of restart,...
This module is intended to contain "only mathematical" functions and procedures.
This module defines routines, generating operators for a cubic stencil.
subroutine, public stencil_cube_get_lapl(this, dim, order)
type(type_t), public type_float
class representing derivatives
Describes mesh distribution to nodes.
Quantities used in mixing: Input, output and new potentials, and the residuals.