Octopus
mix.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2014 M. Marques, A. Castro, A. Rubio, G. Bertsch, M. Oliveira
2!! Copyright (C) 2024 A. Buccheri
3!!
4!! This program is free software; you can redistribute it and/or modify
5!! it under the terms of the GNU General Public License as published by
6!! the Free Software Foundation; either version 2, or (at your option)
7!! any later version.
8!!
9!! This program is distributed in the hope that it will be useful,
10!! but WITHOUT ANY WARRANTY; without even the implied warranty of
11!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12!! GNU General Public License for more details.
13!!
14!! You should have received a copy of the GNU General Public License
15!! along with this program; if not, write to the Free Software
16!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
17!! 02110-1301, USA.
18!!
19
20#include "global.h"
21
22module mix_oct_m
23 use comm_oct_m
24 use debug_oct_m
26 use global_oct_m
27 use, intrinsic :: iso_fortran_env
30 use math_oct_m
31 use mesh_oct_m
36 use parser_oct_m
41 use space_oct_m
43 use types_oct_m
45
46 implicit none
47
48 private
49 public :: &
50 mix_t, &
51 mix_init, &
52 mix_clear, &
53 mix_end, &
54 mix_dump, &
55 mix_load, &
56 mixing, &
57 dmixing, &
58 zmixing, &
60 mix_scheme, &
61 mix_d3, &
63 mixfield_t, &
71
72 ! A better design would be to have mixfield_t and mixfield_ptr_t separate from mix_t:
73 ! mix_t would contain constants, and should be immutable with only getters.
74 ! mixfield_type encapsulates the data used in mixing, and has bound procedures for
75 ! operating on this data.
76
78 type mixfield_t
79 private
80 ! Real versions of the data
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(:, :)
89
90 ! Complex versions of the data
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(:, :)
99
100 type(type_t) :: func_type
101 integer :: d1, d2, d3
102 ! d1 is mesh%np, d2 vary depending on where this called
103 ! d3 is for SCF history of a quantity.
104 end type mixfield_t
105
106 type mixfield_ptr_t
107 private
108 type(mixfield_t), pointer :: p
109 end type mixfield_ptr_t
110
111 integer, parameter :: MAX_AUXMIXFIELD = 5
112
120 type mix_t
121 private
122 type(derivatives_t), pointer :: der
123
124 integer :: scheme
125 real(real64) :: coeff
126 real(real64) :: residual_coeff
127 integer :: iter
128 ! !< include the iterations done in previous calculations.
129 integer, public :: ns
130 integer, public :: ns_restart
131 integer :: interval
132
133 integer :: last_ipos
134 ! This should have been an attribute of type(mixfield_t)
135
136 ! Fields
137 type(mixfield_t) :: mixfield
138 integer :: nauxmixfield
139 type(mixfield_ptr_t) :: auxmixfield(MAX_AUXMIXFIELD)
140 integer :: ipos
141 ! This should have been an attribute of type(mixfield_ptr_t)
142
143 ! Preconditioning
144 logical :: kerker
145 real(real64) :: kerker_factor
146 logical :: precondition
147 type(nl_operator_t) :: preconditioner
148 ! [GPAW](https://wiki.fysik.dtu.dk/gpaw/documentation/densitymix/densitymix.html)
149
150 contains
151 procedure :: compute_residuals_aux_field
152
153 end type mix_t
154
155 interface mixfield_set_vin
157 end interface mixfield_set_vin
158
159 interface mixfield_set_vout
161 end interface mixfield_set_vout
162
163 interface mixfield_get_vnew
165 end interface mixfield_get_vnew
166
167
168contains
169
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
174 class(space_t), intent(in) :: space
175 type(derivatives_t), target, intent(in) :: der
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_
182 integer :: def, ii
183 character(len=32) :: prefix
184 type(type_t) :: func_type
186 push_sub(mix_init)
188 smix%der => der
190 def = option__mixingscheme__broyden
191 if (present(def_)) def = def_
192 if (present(func_type_)) then
193 func_type = func_type_
194 else
195 func_type = type_float
196 end if
197 prefix = ''
198 if (present(prefix_)) prefix = prefix_
200 call messages_obsolete_variable(namespace, 'TypeOfMixing', 'MixingScheme')
202 !%Variable MixingScheme
203 !%Type integer
204 !%Default broyden
205 !%Section SCF::Mixing
206 !%Description
207 !% The scheme used to produce, at each iteration in the self-consistent cycle
208 !% that attempts to solve the Kohn-Sham equations, the input density from the value
209 !% of the input and output densities of previous iterations.
210 !%Option linear 0
211 !% Simple linear mixing.
212 !%Option broyden 2
213 !% Broyden scheme [C. G Broyden, <i>Math. Comp.</i> <b>19</b>, 577 (1965);
214 !% D. D. Johnson, <i>Phys. Rev. B</i> <b>38</b>, 12807 (1988)].
215 !% The scheme is slightly adapted, see the comments in the code.
216 !% For complex functions (e.g. Sternheimer with <tt>EMEta</tt> > 0), we use the generalization
217 !% with a complex dot product.
218 !%Option diis 9
219 !% Direct inversion in the iterative subspace (diis)
220 !% scheme [P. Pulay, <i>Chem. Phys. Lett.</i>, <b>73</b>, 393
221 !% (1980)] as described in [G. Kresse, and J. Furthmueller,
222 !% <i>Phys. Rev. B</i> <b>54</b>, 11169 (1996)].
223 !%End
224 call parse_variable(namespace, trim(prefix)//'MixingScheme', def, smix%scheme)
225 if (.not. varinfo_valid_option('MixingScheme', smix%scheme)) then
226 call messages_input_error(namespace, 'MixingScheme', 'invalid option')
227 end if
228 call messages_print_var_option("MixingScheme", smix%scheme, namespace=namespace)
229
230 !%Variable MixingPreconditioner
231 !%Type logical
232 !%Default false
233 !%Section SCF::Mixing
234 !%Description
235 !% (Experimental) If set to yes, Octopus will use a preconditioner
236 !% for the mixing operator.
237 !% This preconditioner is disabled for systems with dimension other than 3.
238 !%End
239 call parse_variable(namespace, trim(prefix)+'MixingPreconditioner', .false., smix%precondition)
240 if (der%dim /= 3) smix%precondition = .false.
241 if (smix%precondition) call messages_experimental('MixingPreconditioner', namespace=namespace)
242
243 !%Variable Mixing
244 !%Type float
245 !%Default 0.3
246 !%Section SCF::Mixing
247 !%Description
248 !% The linear, Broyden and DIIS scheme depend on a "mixing parameter", set by this variable.
249 !% Must be 0 < <tt>Mixing</tt> <= 1.
250 !%End
251 call parse_variable(namespace, trim(prefix)+'Mixing', 0.3_real64, smix%coeff)
252 if (smix%coeff <= m_zero .or. smix%coeff > m_one) then
253 call messages_input_error(namespace, 'Mixing', 'Value should be positive and smaller than one')
254 end if
255
256 !%Variable MixingResidual
257 !%Type float
258 !%Default 0.05
259 !%Section SCF::Mixing
260 !%Description
261 !% In the DIIS mixing it is benefitial to include a bit of
262 !% residual into the mixing. This parameter controls this amount.
263 !%End
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.')
267 end if
268
269 !%Variable MixNumberSteps
270 !%Type integer
271 !%Default 4
272 !%Section SCF::Mixing
273 !%Description
274 !% In the Broyden scheme, the new input density or potential is constructed
275 !% from the values of the densities/potentials of a given number of previous iterations.
276 !% This number is set by this variable. Must be greater than 1.
277 !%End
278 if (smix%scheme /= option__mixingscheme__linear) then
279 call parse_variable(namespace, trim(prefix)//'MixNumberSteps', 4, smix%ns)
280 if (smix%ns <= 1) call messages_input_error(namespace, 'MixNumberSteps')
281 else
282 smix%ns = 0
283 end if
284
285 !%Variable MixingRestart
286 !%Type integer
287 !%Default 20
288 !%Section SCF::Mixing
289 !%Description
290 !% In the Broyden scheme, the mixing is restarted after
291 !% the number of iterations given by this variable.
292 !% Set this to zero to disable restarting the mixing.
293 !%End
294 if (smix%scheme /= option__mixingscheme__linear) then
295 call parse_variable(namespace, trim(prefix)//'MixingRestart', 20, smix%ns_restart)
296 if (smix%ns_restart < 0) call messages_input_error(namespace, 'MixingRestart')
297 else
298 smix%ns_restart = 0
299 end if
300
301 write(message(1), '(A,I4,A,I4,A)') "Info: Mixing uses ", smix%ns, " steps and restarts after ", &
302 smix%ns_restart, " steps."
303 call messages_info(1, namespace=namespace)
304
305 !%Variable MixInterval
306 !%Type integer
307 !%Default 1
308 !%Section SCF::Mixing
309 !%Description
310 !% When this variable is set to a value different than 1 (the
311 !% default) a combined mixing scheme will be used, with MixInterval
312 !% - 1 steps of linear mixing followed by 1 step of the selected
313 !% mixing. For the moment this variable only works with DIIS mixing.
314 !%End
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')
317
318 smix%iter = 0
319
320 !%Variable MixingKerker
321 !%Type logical
322 !%Default false
323 !%Section SCF::Mixing
324 !%Description
325 !% (Experimental) If set to yes, Octopus will use the Kerker preconditioner
326 !% for the mixing operator applied in the linear and Broyden mixing schemes.
327 !% When using preconditioning with the Broyden scheme, the user should increase the <tt>Mixing</tt> from its default
328 !% to between 0.5 and 0.8, for optimal results.
329 !% The implementation follows Yoshinori Shiihara et al., Modelling Simul. Mater. Sci. Eng. 16 (2008), 035004
330 !%End
331 call parse_variable(namespace, trim(prefix)+'MixingKerker', .false., smix%kerker)
332 if(smix%kerker) call messages_experimental('MixingKerker')
333
334 !%Variable MixingKerkerFactor
335 !%Type float
336 !%Default 1.0
337 !%Section SCF::Mixing
338 !%Description
339 !% The screening factor, $q_0$ in the Kerker preconditioner in units of inverse length.
340 !% For small wave vectors, q, screening behaves like $\frac{Aq}{q_0}$, where A is the mixing coefficient.
341 !% An optimal value is determined empirically for a given system and mixing coeffient, however the default
342 !% is a reasonable choice in many cases.
343 !%End
344 call parse_variable(namespace, trim(prefix)+'MixingKerkerFactor', 1._real64, smix%kerker_factor)
345
346 smix%nauxmixfield = 0
347 do ii = 1,max_auxmixfield
348 nullify(smix%auxmixfield(ii)%p)
349 end do
350
351 call mixfield_init(smix, smix%mixfield, d1, d2, smix%ns, func_type)
352 call mix_clear(smix)
353 if (smix%precondition) call init_preconditioner()
354
355 pop_sub(mix_init)
356
357 contains
358
359 subroutine init_preconditioner()
360
361 integer :: ns, maxp, ip, is
362 real(real64), parameter :: weight = 50.0_real64
363
364 ! This the mixing preconditioner from GPAW:
365 !
366 ! https://wiki.fysik.dtu.dk/gpaw/documentation/densitymix/densitymix.html
367 !
368
369 assert(.not. der%mesh%use_curvilinear)
370 assert(der%dim == 3)
371
372 call nl_operator_init(smix%preconditioner, "Mixing preconditioner")
373 call stencil_cube_get_lapl(smix%preconditioner%stencil, der%dim, 1)
374 call nl_operator_build(space, der%mesh, smix%preconditioner, der%mesh%np, const_w = .not. der%mesh%use_curvilinear)
375
376 ns = smix%preconditioner%stencil%size
377
378 if (smix%preconditioner%const_w) then
379 maxp = 1
380 else
381 maxp = der%mesh%np
382 end if
383
384 do ip = 1, maxp
385
386 do is = 1, ns
387 select case (sum(abs(smix%preconditioner%stencil%points(1:der%dim, is))))
388 case (0)
389 smix%preconditioner%w(is, ip) = m_one + weight/8.0_real64
390 case (1)
391 smix%preconditioner%w(is, ip) = weight/16.0_real64
392 case (2)
393 smix%preconditioner%w(is, ip) = weight/32.0_real64
394 case (3)
395 smix%preconditioner%w(is, ip) = weight/64.0_real64
396 case default
397 assert(.false.)
398 end select
399
400 end do
401 end do
402
403 call nl_operator_allocate_gpu_buffers(smix%preconditioner)
404 call nl_operator_update_gpu_buffers(smix%preconditioner)
405 call nl_operator_output_weights(smix%preconditioner)
406
407 end subroutine init_preconditioner
408
409 end subroutine mix_init
410
411
412 ! ---------------------------------------------------------
413 subroutine mix_clear(smix)
414 type(mix_t), intent(inout) :: smix
415
416 push_sub(mix_clear)
417
418 call mixfield_clear(smix%scheme, smix%mixfield)
419
420 smix%iter = 0
421 smix%last_ipos = 0
422
423 pop_sub(mix_clear)
424 end subroutine mix_clear
425
426
427 ! ---------------------------------------------------------
428 subroutine mix_end(smix)
429 type(mix_t), intent(inout) :: smix
430
431 integer :: ii
432
433 push_sub(mix_end)
434
435 if (smix%precondition) call nl_operator_end(smix%preconditioner)
436
437 call mixfield_end(smix, smix%mixfield)
438
439 smix%nauxmixfield = 0
440 do ii = 1,max_auxmixfield
441 nullify(smix%auxmixfield(ii)%p)
442 end do
443
444 pop_sub(mix_end)
445 end subroutine mix_end
446
447
448 ! ---------------------------------------------------------
449 subroutine mix_dump(namespace, restart, smix, space, mesh, ierr)
450 type(namespace_t), intent(in) :: namespace
451 type(restart_t), intent(in) :: restart
452 type(mix_t), intent(in) :: smix
453 class(space_t), intent(in) :: space
454 class(mesh_t), intent(in) :: mesh
455 integer, intent(out) :: ierr
456
457 integer :: iunit, id2, id3, err, err2(4)
458 character(len=40) :: lines(7)
459 character(len=80) :: filename
460
461 push_sub(mix_dump)
462
463 ierr = 0
464
465 if (restart_skip(restart)) then
466 pop_sub(mix_dump)
467 return
468 end if
469
470 message(1) = "Debug: Writing mixing restart."
471 call messages_info(1, namespace=namespace, debug_only=.true.)
472
473 ! functions to be written need to be compatible with the mesh
474 assert(mesh%np == smix%mixfield%d1)
475
476 ! First we write some information about the mixing
477 iunit = restart_open(restart, 'mixing')
478 write(lines(1), '(a11,i1)') 'scheme= ', smix%scheme
479 ! Number of global mesh points have to be written, not only smix%d1
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
486 call restart_write(restart, iunit, lines, 7, err)
487 ierr = ierr + err
488 call restart_close(restart, iunit)
489
490 ! Now we write the different functions.
491 ! These are not needed when using linear mixing, so we will make sure we skip this step in that case.
492 err2 = 0
493 if (smix%scheme /= option__mixingscheme__linear) then
494 do id2 = 1, smix%mixfield%d2
495 do id3 = 1, smix%mixfield%d3
496
497 write(filename,'(a3,i2.2,i2.2)') 'df_', id2, id3
498 if (smix%mixfield%func_type == type_float) then
499 call restart_write_mesh_function(restart, space, filename, mesh, smix%mixfield%ddf(1:mesh%np, id2, id3), err)
500 else
501 call restart_write_mesh_function(restart, space, filename, mesh, smix%mixfield%zdf(1:mesh%np, id2, id3), err)
502 end if
503 err2(1) = err2(1) + err
504
505 write(filename,'(a3,i2.2,i2.2)') 'dv_', id2, id3
506 if (smix%mixfield%func_type == type_float) then
507 call restart_write_mesh_function(restart, space, filename, mesh, smix%mixfield%ddv(1:mesh%np, id2, id3), err)
508 else
509 call restart_write_mesh_function(restart, space, filename, mesh, smix%mixfield%zdv(1:mesh%np, id2, id3), err)
510 end if
511 err2(2) = err2(2) + err
512
513 end do
514
515 write(filename,'(a6,i2.2)') 'f_old_', id2
516 if (smix%mixfield%func_type == type_float) then
517 call restart_write_mesh_function(restart, space, filename, mesh, smix%mixfield%df_old(1:mesh%np, id2), err)
518 else
519 call restart_write_mesh_function(restart, space, filename, mesh, smix%mixfield%zf_old(1:mesh%np, id2), err)
520 end if
521 err2(3) = err2(3) + err
522
523 write(filename,'(a8,i2.2)') 'vin_old_', id2
524 if (smix%mixfield%func_type == type_float) then
525 call restart_write_mesh_function(restart, space, filename, mesh, smix%mixfield%dvin_old(1:mesh%np, id2), err)
526 else
527 call restart_write_mesh_function(restart, space, filename, mesh, smix%mixfield%zvin_old(1:mesh%np, id2), err)
528 end if
529 err2(4) = err2(4) + err
530
531 end do
532
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
537 end if
538
539 message(1) = "Debug: Writing mixing restart done."
540 call messages_info(1, namespace=namespace, debug_only=.true.)
541
542 pop_sub(mix_dump)
543 end subroutine mix_dump
544
545
546 !---------------------------------------------------------
547 subroutine mix_load(namespace, restart, smix, space, mesh, ierr)
548 type(namespace_t), intent(in) :: namespace
549 type(restart_t), intent(in) :: restart
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
554
555 integer :: iunit, err, err2(4)
556 integer :: scheme, d1, d2, d3, ns
557 integer :: id2, id3
558 character(len=11) :: str
559 character(len=80) :: filename
560 character(len=256) :: lines(7)
561
562 push_sub(mix_load)
563
564 ierr = 0
565
566 if (restart_skip(restart)) then
567 ierr = -1
568 pop_sub(mix_load)
569 return
570 end if
571
572 message(1) = "Debug: Reading mixing restart."
573 call messages_info(1, namespace=namespace, debug_only=.true.)
574
575 ! First we read some information about the mixing
576 iunit = restart_open(restart, 'mixing')
577 call restart_read(restart, iunit, lines, 7, err)
578 if (err /= 0) then
579 ierr = ierr + 1
580 else
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
588 end if
589 call restart_close(restart, iunit)
590
591
592 if (ierr == 0) then
593 ! We can only use the restart information if the mixing scheme and the number of steps used remained the same
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."
596 call messages_warning(1, namespace=namespace)
597 ierr = ierr + 2
598 end if
599
600 ! Check the dimensions of the arrays to be read
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."
604 call messages_warning(2, namespace=namespace)
605 ierr = ierr + 4
606 end if
607 end if
608
609
610 ! Now we read the different functions.
611 ! Note that we may have more or less functions than the ones needed (d3 /= smix%d3)
612 if (ierr == 0) then
613 if (smix%scheme /= option__mixingscheme__linear) then
614 err2 = 0
615 do id2 = 1, smix%mixfield%d2
616 do id3 = 1, smix%mixfield%d3
617
618 write(filename,'(a3,i2.2,i2.2)') 'df_', id2, id3
619 if (smix%mixfield%func_type == type_float) then
620 call restart_read_mesh_function(restart, space, filename, mesh, smix%mixfield%ddf(1:mesh%np, id2, id3), err)
621 else
622 call restart_read_mesh_function(restart, space, filename, mesh, smix%mixfield%zdf(1:mesh%np, id2, id3), err)
623 end if
624 if (err /= 0) err2(1) = err2(1) + 1
625
626 write(filename,'(a3,i2.2,i2.2)') 'dv_', id2, id3
627 if (smix%mixfield%func_type == type_float) then
628 call restart_read_mesh_function(restart, space, filename, mesh, smix%mixfield%ddv(1:mesh%np, id2, id3), err)
629 else
630 call restart_read_mesh_function(restart, space, filename, mesh, smix%mixfield%zdv(1:mesh%np, id2, id3), err)
631 end if
632 if (err /= 0) err2(2) = err2(2) + 1
633
634 end do
635
636 write(filename,'(a6,i2.2)') 'f_old_', id2
637 if (smix%mixfield%func_type == type_float) then
638 call restart_read_mesh_function(restart, space, filename, mesh, smix%mixfield%df_old(1:mesh%np, id2), err)
639 else
640 call restart_read_mesh_function(restart, space, filename, mesh, smix%mixfield%zf_old(1:mesh%np, id2), err)
641 end if
642 if (err /= 0) err2(3) = err2(3) + 1
643
644 write(filename,'(a8,i2.2)') 'vin_old_', id2
645 if (smix%mixfield%func_type == type_float) then
646 call restart_read_mesh_function(restart, space, filename, mesh, smix%mixfield%dvin_old(1:mesh%np, id2), err)
647 else
648 call restart_read_mesh_function(restart, space, filename, mesh, smix%mixfield%zvin_old(1:mesh%np, id2), err)
649 end if
650 if (err /= 0) err2(4) = err2(4) + 1
651
652 end do
653
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
658 end if
659 end if
660
661 if (ierr /= 0) then
662 ! Something went wront, so make sure we start from scratch
663 call mix_clear(smix)
664 end if
665
666 message(1) = "Debug: Reading mixing restart done."
667 call messages_info(1, namespace=namespace, debug_only=.true.)
668
669 pop_sub(mix_load)
670 end subroutine mix_load
671
672 real(real64) pure function mix_coefficient(this) result(coefficient)
673 type(mix_t), intent(in) :: this
674
675 coefficient = this%coeff
676 end function mix_coefficient
677
678 integer pure function mix_scheme(this) result(scheme)
679 type(mix_t), intent(in) :: this
680
681 scheme = this%scheme
682 end function mix_scheme
683
684 integer pure function mix_d3(this)
685 type(mix_t), intent(in) :: this
686
687 mix_d3 = this%mixfield%d3
688 end function mix_d3
689
690 subroutine mix_get_field(this, mixfield)
691 type(mix_t), target, intent(in) :: this
692 type(mixfield_t), pointer, intent(out) :: mixfield
693
694 mixfield => this%mixfield
695 end subroutine mix_get_field
696
698 subroutine mixing(namespace, smix)
699 type(namespace_t), intent(in) :: namespace
700 type(mix_t), intent(inout) :: smix
701
702 push_sub(mixing)
703
704 if (smix%mixfield%func_type == type_float) then
705 call dmixing(namespace, smix, smix%mixfield%dvin, smix%mixfield%dvout, smix%mixfield%dvnew)
706 else
707 call zmixing(namespace, smix, smix%mixfield%zvin, smix%mixfield%zvout, smix%mixfield%zvnew)
708 end if
709
710 pop_sub(mixing)
711 end subroutine mixing
712
713 subroutine mix_add_auxmixfield(namespace, smix, mixfield)
714 type(namespace_t), intent(in) :: namespace
715 type(mix_t), intent(inout) :: smix
716 type(mixfield_t), target, intent(in) :: mixfield
717
718 push_sub(mix_add_auxmixfield)
719
720 smix%nauxmixfield = smix%nauxmixfield + 1
721 smix%auxmixfield(smix%nauxmixfield)%p => mixfield
722
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)
726 end if
727
728 pop_sub(mix_add_auxmixfield)
729 end subroutine mix_add_auxmixfield
730
732 subroutine mixfield_init(smix, mixfield, d1, d2, d3, func_type)
733 type(mix_t), intent(in) :: smix
734 type(mixfield_t), intent(inout) :: mixfield
735 integer, intent(in) :: d1, d2, d3
736 type(type_t), intent(in) :: func_type
737
738 push_sub(mixfield_init)
739
740 mixfield%d1 = d1
741 mixfield%d2 = d2
742 mixfield%d3 = d3
743
744 mixfield%func_type = func_type
745
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))
752 else
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))
757 end if
758 end if
759
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))
765 else
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))
770 end if
772 pop_sub(mixfield_init)
773 end subroutine mixfield_init
774
776 subroutine mixfield_end(smix, mixfield)
777 type(mix_t), intent(inout) :: smix
778 type(mixfield_t), intent(inout) :: mixfield
779
780 push_sub(mixfield_end)
781
782 ! Arrays got allocated for all mixing schemes, except linear mixing
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)
789 else
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)
794 end if
795 end if
796
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)
802 else
803 safe_deallocate_a(mixfield%zvin)
804 safe_deallocate_a(mixfield%zvout)
805 safe_deallocate_a(mixfield%zvnew)
806 safe_deallocate_a(mixfield%zresidual)
807 end if
808
809 pop_sub(mixfield_end)
810 end subroutine mixfield_end
811
813 subroutine mixfield_clear(scheme, mixfield)
814 integer, intent(in) :: scheme
815 type(mixfield_t), intent(inout) :: mixfield
816 integer :: d1, d2, d3
817
818 push_sub(mixfield_clear)
819
820 d1 = mixfield%d1
821 d2 = mixfield%d2
822 d3 = mixfield%d3
823
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
831 else
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
837 end if
838 end if
839
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
845 else
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
850 end if
851
852 pop_sub(mixfield_clear)
853 end subroutine mixfield_clear
854
855 ! While it would be desirable to calve off auxillary fields and associated routines to separate modules
856 ! the Broyden mixing routines make calls to auxilliary fields. One would need to separate those
857 ! dependencies first.
858
859 subroutine ddmixfield_set_vin(mixfield, vin_re, vin_im)
860 type(mixfield_t), intent(inout) :: mixfield
861 real(real64), intent(in) :: vin_re(:,:), vin_im(:,:)
862
863 push_sub(ddmixfield_set_vin)
864
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)
867
868 pop_sub(ddmixfield_set_vin)
869 end subroutine ddmixfield_set_vin
870
871 subroutine ddmixfield_set_vout(mixfield, vout_re, vout_im)
872 type(mixfield_t), intent(inout) :: mixfield
873 real(real64), intent(in) :: vout_re(:,:), vout_im(:,:)
874
875 push_sub(ddmixfield_set_vout)
876
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)
879
880 pop_sub(ddmixfield_set_vout)
881 end subroutine ddmixfield_set_vout
882
883 subroutine ddmixfield_get_vnew(mixfield, re, im)
884 type(mixfield_t), intent(in) :: mixfield
885 real(real64), intent(inout) :: re(:,:), im(:,:)
886
887 push_sub(mixfield_get_ddvnew)
888
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))
891
892 pop_sub(mixfield_get_ddvnew)
893 end subroutine ddmixfield_get_vnew
894
903 subroutine compute_residuals_aux_field(this)
904 class(mix_t), intent(inout) :: this
905
906 type(mixfield_t), pointer :: aux_field
907 integer :: dims(2)
908 integer :: i
909
911
912 do i = 1, this%nauxmixfield
913 aux_field => this%auxmixfield(i)%p
914 dims = [aux_field%d1, aux_field%d2]
915
916 if (aux_field%func_type == type_float) then
917 call dcompute_residual(this, dims, aux_field%dvin, aux_field%dvout, aux_field%dresidual)
918 else
919 call zcompute_residual(this, dims, aux_field%zvin, aux_field%zvout, aux_field%zresidual)
920 end if
921
922 end do
923 nullify(aux_field)
924
926
927 end subroutine compute_residuals_aux_field
928
940 subroutine linear_mixing_aux_field(smix)
941 type(mix_t), intent(inout) :: smix
942
943 type(mixfield_t), pointer :: field
944 integer :: i
945
947
948 call smix%compute_residuals_aux_field()
949
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)
955 else
956 call zmixing_linear(field%d1, field%d2, smix%coeff, field%zvin, field%zresidual, field%zvnew)
957 end if
958
959 end do
960
961 nullify(field)
962
965
966#include "undef.F90"
967#include "real.F90"
968
969#include "mix_inc.F90"
970
971#include "undef.F90"
972#include "complex.F90"
973
974#include "mix_inc.F90"
975
976end module mix_oct_m
977
978!! Local Variables:
979!! mode: f90
980!! coding: utf-8
981!! End:
subroutine init_preconditioner()
Definition: mix.F90:453
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
real(real64), parameter, public m_zero
Definition: global.F90:187
real(real64), parameter, public m_one
Definition: global.F90:188
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:115
This module defines various routines, operating on mesh functions.
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:118
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:543
subroutine, public messages_obsolete_variable(namespace, name, rep)
Definition: messages.F90:1057
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:160
subroutine, public messages_input_error(namespace, var, details, row, column)
Definition: messages.F90:723
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1097
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:624
subroutine, public mixfield_end(smix, mixfield)
Deallocate all arrays of a mixfield instance.
Definition: mix.F90:870
integer pure function, public mix_scheme(this)
Definition: mix.F90:772
subroutine dmixfield_set_vin(mixfield, vin)
Definition: mix.F90:1671
real(real64) pure function, public mix_coefficient(this)
Definition: mix.F90:766
subroutine, public mixing(namespace, smix)
Main entry-point to SCF mixer.
Definition: mix.F90:792
subroutine, public mixfield_init(smix, mixfield, d1, d2, d3, func_type)
Initialise all attributes of a mixfield instance.
Definition: mix.F90:826
subroutine ddmixfield_get_vnew(mixfield, re, im)
Definition: mix.F90:977
subroutine zmixfield_get_vnew(mixfield, vnew)
Definition: mix.F90:2348
subroutine dcompute_residual(smix, dims, vin, vout, residual)
Compute the residual of the potential (or density) for SCF mixing.
Definition: mix.F90:1219
subroutine, public mix_get_field(this, mixfield)
Definition: mix.F90:784
subroutine compute_residuals_aux_field(this)
Compute the residuals of the auxilliary fields.
Definition: mix.F90:997
subroutine zmixfield_set_vout(mixfield, vout)
Definition: mix.F90:2336
subroutine, public mix_add_auxmixfield(namespace, smix, mixfield)
Definition: mix.F90:807
subroutine dmixfield_get_vnew(mixfield, vnew)
Definition: mix.F90:1695
subroutine dmixing_linear(d1, d2, coeff, vin, residual, vnew)
Linear mixing of the input and output potentials.
Definition: mix.F90:1187
subroutine, public mix_dump(namespace, restart, smix, space, mesh, ierr)
Definition: mix.F90:543
subroutine, public mixfield_clear(scheme, mixfield)
Zero all potential and field attributes of a mixfield instance.
Definition: mix.F90:907
subroutine, public mix_load(namespace, restart, smix, space, mesh, ierr)
Definition: mix.F90:641
subroutine, public mix_init(smix, namespace, space, der, d1, d2, def_, func_type_, prefix_)
Initialise mix_t instance.
Definition: mix.F90:265
subroutine ddmixfield_set_vin(mixfield, vin_re, vin_im)
Definition: mix.F90:953
subroutine zmixfield_set_vin(mixfield, vin)
Definition: mix.F90:2324
subroutine, public zmixing(namespace, smix, vin, vout, vnew)
Mix the input and output potentials (or densities) according to some scheme.
Definition: mix.F90:1786
integer pure function, public mix_d3(this)
Definition: mix.F90:778
subroutine dmixfield_set_vout(mixfield, vout)
Definition: mix.F90:1683
subroutine zcompute_residual(smix, dims, vin, vout, residual)
Compute the residual of the potential (or density) for SCF mixing.
Definition: mix.F90:1872
subroutine ddmixfield_set_vout(mixfield, vout_re, vout_im)
Definition: mix.F90:965
subroutine, public dmixing(namespace, smix, vin, vout, vnew)
Mix the input and output potentials (or densities) according to some scheme.
Definition: mix.F90:1133
subroutine linear_mixing_aux_field(smix)
Linear mixing of the auxilliary fields.
Definition: mix.F90:1034
subroutine zmixing_linear(d1, d2, coeff, vin, residual, vnew)
Linear mixing of the input and output potentials.
Definition: mix.F90:1840
subroutine, public mix_end(smix)
Definition: mix.F90:522
subroutine, public mix_clear(smix)
Definition: mix.F90:507
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)
Definition: restart.F90:935
subroutine, public restart_close(restart, iunit)
Close a file previously opened with restart_open.
Definition: restart.F90:952
subroutine, public restart_write(restart, iunit, lines, nlines, ierr)
Definition: restart.F90:908
logical pure function, public restart_skip(restart)
Returns true if the restart information should neither be read nor written. This might happen because...
Definition: restart.F90:969
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,...
Definition: restart.F90:861
This module is intended to contain "only mathematical" functions and procedures.
Definition: solvers.F90:115
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
Definition: types.F90:133
class representing derivatives
Describes mesh distribution to nodes.
Definition: mesh.F90:186
God class for mixing.
Definition: mix.F90:213
Quantities used in mixing: Input, output and new potentials, and the residuals.
Definition: mix.F90:171
int true(void)