Octopus
opt_control.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
19#include "global.h"
20
26 use debug_oct_m
30 use filter_oct_m
31 use global_oct_m
32 use grid_oct_m
34 use initst_oct_m
35 use ions_oct_m
36 use iso_c_binding
37 use, intrinsic :: iso_fortran_env
38 use io_oct_m
40 use lasers_oct_m
41 use loct_oct_m
42 use math_oct_m
50 use output_oct_m
51 use parser_oct_m
60 use target_oct_m
62 use td_oct_m
63
64 implicit none
65
66 private
67 public :: &
74
75
77 type(filter_t), save :: filter
78 type(oct_t), save :: oct
79 type(oct_iterator_t), save :: iterator
80 type(target_t), save :: oct_target
81 type(opt_control_state_t), save :: initial_st
82
83
85 type(controlfunction_t), save :: par_
86 type(electrons_t), pointer :: sys_
87 type(hamiltonian_elec_t), pointer :: hm_
88 type(td_t), pointer :: td_
89 real(real64), allocatable :: x_(:)
90 integer :: index_
91
92contains
93
94 ! ---------------------------------------------------------
95 subroutine opt_control_run(system)
96 class(*), intent(inout) :: system
97
98 push_sub(opt_control_run)
99
100 select type (system)
101 class is (multisystem_basic_t)
102 message(1) = "CalculationMode = opt_control not implemented for multi-system calculations"
103 call messages_fatal(1, namespace=system%namespace)
104 type is (electrons_t)
105 call opt_control_run_legacy(system)
106 end select
107
108 pop_sub(opt_control_run)
109 end subroutine opt_control_run
110
113 subroutine opt_control_run_legacy(sys)
114 type(electrons_t), target, intent(inout) :: sys
115
116 type(td_t), target :: td
117 type(controlfunction_t) :: par, par_new, par_prev
118 logical :: stop_loop
119 real(real64) :: j1
120 type(oct_prop_t) :: prop_chi, prop_psi
121 type(states_elec_t) :: psi
122
123 type(lasers_t), pointer :: lasers
124
125 push_sub(opt_control_run_legacy)
126
127 if (sys%hm%pcm%run_pcm) then
128 call messages_not_implemented("PCM for CalculationMode /= gs or td", namespace=sys%namespace)
129 end if
130
131 if (sys%kpoints%use_symmetries) then
132 call messages_experimental("KPoints symmetries with CalculationMode = opt_control", namespace=sys%namespace)
133 end if
134
135 ! Creates a directory where the optimal control stuff will be written. The name of the directory
136 ! is stored in the preprocessor macro OCT_DIR, which should be defined in src/include/global.h
137 call io_mkdir(oct_dir, sys%namespace)
138
139 ! Initializes the time propagator. Then, it forces the propagation to be self consistent, in case
140 ! the theory level is not "independent particles".
141 call td_init(td, sys%namespace, sys%space, sys%gr, sys%ions, sys%st, sys%ks, sys%hm, sys%ext_partners, sys%outp)
142 if (sys%hm%theory_level /= independent_particles) call propagator_elec_set_scf_prop(td%tr, threshold = 1.0e-14_real64)
143
144 ! Read general information about how the OCT run will be made, from inp file. "oct_read_inp" is
145 ! in the opt_control_global_oct_m module (like the definition of the oct_t data type)
146 call oct_read_inp(oct, sys%namespace)
147
148 ! Read info about, and prepare, the control functions
149 call controlfunction_mod_init(sys%ext_partners, sys%namespace, td%dt, td%max_iter, oct%mode_fixed_fluence)
150 call controlfunction_init(par, td%dt, td%max_iter)
151 call controlfunction_set(par, sys%ext_partners)
152 ! This prints the initial control parameters, exactly as described in the inp file,
153 ! that is, without applying any envelope or filter.
154 call controlfunction_write(oct_dir//'initial_laser_inp', par, sys%namespace)
156 call controlfunction_to_h(par, sys%ext_partners)
157 call messages_print_with_emphasis(msg="TD ext. fields after processing", namespace=sys%namespace)
158 lasers => list_get_lasers(sys%ext_partners)
159 if(associated(lasers)) call laser_write_info(lasers%lasers, namespace=sys%namespace)
160 call messages_print_with_emphasis(namespace=sys%namespace)
161 call controlfunction_write(oct_dir//'initial_laser', par, sys%namespace)
162
163
164 ! Startup of the iterator data type (takes care of counting iterations, stopping, etc).
165 call oct_iterator_init(iterator, sys%namespace, par)
166
167
168 ! Initialization of the propagation_oct_m module.
169 call propagation_mod_init(td%max_iter, oct%eta, oct%delta, oct%number_checkpoints, &
170 (oct%algorithm == option__octscheme__oct_zbr98), &
171 (oct%algorithm == option__octscheme__oct_cg) .or. &
172 (oct%algorithm == option__octscheme__oct_bfgs) .or. &
173 (oct%algorithm == option__octscheme__oct_nlopt_lbfgs))
175
176 ! If filters are to be used, they also have to be initialized.
177 call filter_init(td%max_iter, sys%namespace, td%dt, filter)
178 call filter_write(filter, sys%namespace)
181 ! Figure out the starting wavefunction(s), and the target.
182 call initial_state_init(sys, initial_st)
183 call target_init(sys%gr, sys%kpoints, sys%namespace, sys%space, sys%ions, initial_st, td, &
184 controlfunction_w0(par), oct_target, oct, sys%hm%ep, sys%mc)
185
186 ! Sanity checks.
187 call check_faulty_runmodes(sys, td%tr)
189
190 ! Informative output.
191 call opt_control_get_qs(psi, initial_st)
192 call output_states(sys%outp, sys%namespace, sys%space, oct_dir//'initial', psi, sys%gr, sys%ions, sys%hm, -1)
193 call target_output(oct_target, sys%namespace, sys%space, sys%gr, oct_dir//'target', sys%ions, sys%hm, sys%outp)
194 call states_elec_end(psi)
195
196
197 ! mode switcher; here is where the real run is made.
198 select case (oct%algorithm)
199 case (option__octscheme__oct_zbr98)
200 message(1) = "Info: Starting OCT iteration using scheme: ZBR98"
201 call messages_info(1, namespace=sys%namespace)
202 call scheme_zbr98()
203 case (option__octscheme__oct_wg05)
204 message(1) = "Info: Starting OCT iteration using scheme: WG05"
205 call messages_info(1, namespace=sys%namespace)
207 case (option__octscheme__oct_zr98)
208 message(1) = "Info: Starting OCT iteration using scheme: ZR98"
209 call messages_info(1, namespace=sys%namespace)
210 call scheme_mt03()
211 case (option__octscheme__oct_mt03)
212 message(1) = "Info: Starting OCT iteration using scheme: MT03"
213 call messages_info(1, namespace=sys%namespace)
214 call scheme_mt03()
215 case (option__octscheme__oct_krotov)
216 message(1) = "Info: Starting OCT iteration using scheme: KROTOV"
217 call messages_info(1, namespace=sys%namespace)
218 call scheme_mt03()
219 case (option__octscheme__oct_straight_iteration)
220 message(1) = "Info: Starting OCT iterations using scheme: STRAIGHT ITERATION"
221 call messages_info(1, namespace=sys%namespace)
223 case (option__octscheme__oct_cg)
224 message(1) = "Info: Starting OCT iterations using scheme: CONJUGATE GRADIENTS"
225 call messages_info(1, namespace=sys%namespace)
226 call scheme_cg()
227 case (option__octscheme__oct_bfgs)
228 message(1) = "Info: Starting OCT iterations using scheme: BFGS"
229 call messages_info(1, namespace=sys%namespace)
230 call scheme_cg()
231 case (option__octscheme__oct_direct)
232 message(1) = "Info: Starting OCT iterations using scheme: DIRECT OPTIMIZATION (NELDER-MEAD)"
233 call messages_info(1, namespace=sys%namespace)
234 call scheme_direct()
235 case (option__octscheme__oct_nlopt_bobyqa)
236 message(1) = "Info: Starting OCT iterations using scheme: DIRECT OPTIMIZATION (NLOPT - BOBYQA)"
237 call messages_info(1, namespace=sys%namespace)
238 call scheme_nlopt()
239 case (option__octscheme__oct_nlopt_lbfgs)
240 message(1) = "Info: Starting OCT iterations using scheme: DIRECT OPTIMIZATION (NLOPT - LBFGS)"
241 call messages_info(1, namespace=sys%namespace)
242 call scheme_nlopt()
243 case default
244 call messages_input_error(sys%namespace, 'OCTScheme')
245 end select
246
247 ! do final test run: propagate initial state with optimal field
248 call oct_finalcheck(sys, td)
249
250 ! clean up
251 call controlfunction_end(par)
252 call oct_iterator_end(iterator, sys%namespace)
253 call filter_end(filter)
254 call td_end(td)
255 call opt_control_state_end(initial_st)
256 call target_end(oct_target, oct)
258
260
261 contains
262
263
264 ! ---------------------------------------------------------
265 subroutine scheme_straight_iteration()
267
269 call controlfunction_copy(par_new, par)
270 ctr_loop: do
271 call controlfunction_copy(par_prev, par)
272 call f_striter(sys, td, par, j1)
273 stop_loop = iteration_manager(sys%namespace, j1, par_prev, par, iterator)
274 if (clean_stop(sys%mc%master_comm) .or. stop_loop) exit ctr_loop
275 end do ctr_loop
276
277 call controlfunction_end(par_new)
278 call controlfunction_end(par_prev)
280 end subroutine scheme_straight_iteration
281 ! ---------------------------------------------------------
282
283
284 ! ---------------------------------------------------------
285 subroutine scheme_mt03()
286 type(opt_control_state_t) :: psi
288
289 call opt_control_state_null(psi)
290 call opt_control_state_copy(psi, initial_st)
291 call oct_prop_init(prop_chi, sys%namespace, "chi", sys%gr, sys%mc)
292 call oct_prop_init(prop_psi, sys%namespace, "psi", sys%gr, sys%mc)
293
294 call controlfunction_copy(par_new, par)
295 ctr_loop: do
296 call controlfunction_copy(par_prev, par)
297 call f_iter(sys, td, psi, par, prop_psi, prop_chi, j1)
298 stop_loop = iteration_manager(sys%namespace, j1, par, par_prev, iterator)
299 if (clean_stop(sys%mc%master_comm) .or. stop_loop) exit ctr_loop
300 end do ctr_loop
301
302 call opt_control_state_end(psi)
303 call oct_prop_end(prop_chi)
304 call oct_prop_end(prop_psi)
305 call controlfunction_end(par_new)
306 call controlfunction_end(par_prev)
308 end subroutine scheme_mt03
309 ! ---------------------------------------------------------
310
311
312 ! ---------------------------------------------------------
313 subroutine scheme_wg05()
314 type(opt_control_state_t) :: psi
316
317 call oct_prop_init(prop_chi, sys%namespace, "chi", sys%gr, sys%mc)
318 call oct_prop_init(prop_psi, sys%namespace, "psi", sys%gr, sys%mc)
319
320 if (oct%mode_fixed_fluence) then
322 end if
323
324 call opt_control_state_null(psi)
325
326 call controlfunction_copy(par_new, par)
327 ctr_loop: do
328 call controlfunction_copy(par_prev, par)
329 call f_wg05(sys, td, psi, par, prop_psi, prop_chi, j1)
330 stop_loop = iteration_manager(sys%namespace, j1, par, par_prev, iterator)
331 if (clean_stop(sys%mc%master_comm) .or. stop_loop) exit ctr_loop
332 end do ctr_loop
333
334 call opt_control_state_end(psi)
335 call oct_prop_end(prop_chi)
336 call oct_prop_end(prop_psi)
337 call controlfunction_end(par_new)
338 call controlfunction_end(par_prev)
340 end subroutine scheme_wg05
341 ! ---------------------------------------------------------
342
343
344 ! ---------------------------------------------------------
345 subroutine scheme_zbr98()
346 type(opt_control_state_t) :: qcpsi
348
349 call opt_control_state_null(qcpsi)
350 call opt_control_state_copy(qcpsi, initial_st)
351 call oct_prop_init(prop_chi, sys%namespace, "chi", sys%gr, sys%mc)
352 call oct_prop_init(prop_psi, sys%namespace, "psi", sys%gr, sys%mc)
353
354 call controlfunction_copy(par_prev, par)
355 call propagate_forward(sys, td, par, oct_target, qcpsi, prop_psi)
356 j1 = target_j1(oct_target, sys%namespace, sys%gr, sys%kpoints, qcpsi)
357 stop_loop = iteration_manager(sys%namespace, j1, par, par_prev, iterator)
358 if (clean_stop(sys%mc%master_comm) .or. stop_loop) then
359 call opt_control_state_end(qcpsi)
360 call oct_prop_end(prop_chi)
361 call oct_prop_end(prop_psi)
363 return
364 end if
365
366 call controlfunction_copy(par_new, par)
367 ctr_loop: do
368 call controlfunction_copy(par_prev, par)
369 call f_zbr98(sys, td, qcpsi, prop_psi, prop_chi, par)
370 j1 = target_j1(oct_target, sys%namespace, sys%gr, sys%kpoints, qcpsi)
371 stop_loop = iteration_manager(sys%namespace, j1, par, par_prev, iterator)
372 if (clean_stop(sys%mc%master_comm) .or. stop_loop) exit ctr_loop
373 end do ctr_loop
374
375 call opt_control_state_end(qcpsi)
376 call oct_prop_end(prop_chi)
377 call oct_prop_end(prop_psi)
379 call controlfunction_end(par_prev)
381 end subroutine scheme_zbr98
382 ! ---------------------------------------------------------
383
384
385 ! ---------------------------------------------------------
386 subroutine scheme_cg()
387 integer :: dof, ierr, maxiter
388 real(real64):: step, minvalue
389 real(real64), allocatable :: theta(:)
390 real(real64), allocatable :: x(:)
391 real(real64) :: f
392 type(opt_control_state_t) :: qcpsi
394
396
397 call opt_control_state_null(qcpsi)
398 call opt_control_state_copy(qcpsi, initial_st)
399 call propagate_forward(sys, td, par, oct_target, qcpsi)
400 f = - target_j1(oct_target, sys%namespace, sys%gr, sys%kpoints, qcpsi, sys%ions) - controlfunction_j2(par)
401 call opt_control_state_end(qcpsi)
402 call iteration_manager_direct(-f, par, iterator, sys)
403 if (oct_iterator_maxiter(iterator) == 0) then
404 ! Nothing to do.
406 return
407 end if
408
409 if (oct%random_initial_guess) call controlfunction_randomize(par)
410
411 ! Set the module pointers, so that the opt_control_cg_calc and opt_control_cg_write_info routines
412 ! can use them.
413 call controlfunction_copy(par_, par)
414 sys_ => sys
415 hm_ => sys%hm
416 td_ => td
417
418 dof = controlfunction_dof(par)
419 safe_allocate(x(1:dof))
420 safe_allocate(theta(1:dof))
421 call controlfunction_get_theta(par, theta)
422 x = theta
423
424 step = oct%direct_step * m_pi
425 maxiter = oct_iterator_maxiter(iterator) - 1
426
427 select case (oct%algorithm)
428 case (option__octscheme__oct_bfgs)
429 call minimize_multidim(minmethod_bfgs2, dof, x, step, 0.1_real64, &
430 real(oct_iterator_tolerance(iterator), real64), real(oct_iterator_tolerance(iterator), real64), &
431 maxiter, opt_control_cg_calc, opt_control_cg_write_info, minvalue, ierr)
432 case (option__octscheme__oct_cg)
433 call minimize_multidim(minmethod_fr_cg, dof, x, step, 0.1_real64, &
434 real(oct_iterator_tolerance(iterator), real64), real(oct_iterator_tolerance(iterator), real64), &
435 maxiter, opt_control_cg_calc, opt_control_cg_write_info, minvalue, ierr)
436 end select
437
438 if (ierr /= 0) then
439 if (ierr <= 1024) then
440 message(1) = "Error occurred during the GSL minimization procedure:"
441 call loct_strerror(ierr, message(2))
442 call messages_fatal(2, namespace=sys%namespace)
443 else
444 message(1) = "The optimization did not meet the convergence criterion."
445 call messages_info(1, namespace=sys%namespace)
446 end if
447 end if
448
449 call controlfunction_end(par_)
450 safe_deallocate_a(x)
451 safe_deallocate_a(theta)
453 end subroutine scheme_cg
454 ! ---------------------------------------------------------
455
456
457 ! ---------------------------------------------------------
458 subroutine scheme_direct()
459 integer :: ierr, maxiter
460 real(real64):: minvalue, step
461 real(real64), allocatable :: theta(:)
462 real(real64), allocatable :: x(:)
463 real(real64) :: f
464 integer :: dim
465 type(opt_control_state_t) :: qcpsi
466
468
470 dim = controlfunction_dof(par)
471
472 call opt_control_state_null(qcpsi)
473 call opt_control_state_copy(qcpsi, initial_st)
474 call propagate_forward(sys, td, par, oct_target, qcpsi)
475 f = - target_j1(oct_target, sys%namespace, sys%gr, sys%kpoints, qcpsi, sys%ions) - controlfunction_j2(par)
476 call opt_control_state_end(qcpsi)
477 call iteration_manager_direct(-f, par, iterator, sys)
478 if (oct_iterator_maxiter(iterator) == 0) then
479 ! Nothing to do.
481 return
482 end if
483
484 safe_allocate(x(1:dim))
485 safe_allocate(theta(1:dim))
486
487 if (oct%random_initial_guess) call controlfunction_randomize(par)
488
489 ! Set the module pointers, so that the opt_control_direct_calc and opt_control_direct_message_info routines
490 ! can use them.
491 call controlfunction_copy(par_, par)
492 sys_ => sys
493 hm_ => sys%hm
494 td_ => td
495
496 ! theta may be in single precision, whereas x is always double precision.
497 call controlfunction_get_theta(par, theta)
498 x = theta
499
500 step = oct%direct_step * m_pi
501 maxiter = oct_iterator_maxiter(iterator)
502
503 call minimize_multidim_nograd(minmethod_nmsimplex, dim, x, step, &
504 real(oct_iterator_tolerance(iterator), real64) , maxiter, &
505 opt_control_direct_calc, opt_control_direct_message_info, minvalue, ierr)
506
507 if (ierr /= 0) then
508 if (ierr <= 1024) then
509 message(1) = "Error occurred during the GSL minimization procedure:"
510 call loct_strerror(ierr, message(2))
511 call messages_fatal(2, namespace=sys%namespace)
512 else
513 message(1) = "The OCT direct optimization did not meet the convergence criterion."
514 call messages_info(1, namespace=sys%namespace)
515 end if
516 end if
517
518 call controlfunction_end(par_)
519 safe_deallocate_a(x)
520 safe_deallocate_a(theta)
522 end subroutine scheme_direct
523 ! ---------------------------------------------------------
524
525
526 ! ---------------------------------------------------------
527 subroutine scheme_nlopt()
528#if defined(HAVE_NLOPT)
529 integer :: method, dim, maxiter, ierr
530 real(real64), allocatable :: x(:), xl(:), xu(:)
531 real(real64) :: step, toldr, minimum, f
532 type(opt_control_state_t) :: qcpsi
534
536
537 call opt_control_state_null(qcpsi)
538 call opt_control_state_copy(qcpsi, initial_st)
539 call propagate_forward(sys, td, par, oct_target, qcpsi)
540 f = - target_j1(oct_target, sys%namespace, sys%gr, sys%kpoints, qcpsi, sys%ions) - controlfunction_j2(par)
541 call opt_control_state_end(qcpsi)
542 call iteration_manager_direct(-f, par, iterator, sys)
543 if (oct_iterator_maxiter(iterator) == 0) then
544 ! Nothing to do.
546 return
547 end if
548
549 dim = controlfunction_dof(par)
550 safe_allocate(x(dim))
551 safe_allocate(xl(1:dim))
552 safe_allocate(xu(1:dim))
553 call controlfunction_bounds(par, xl, xu)
554
555 if (oct%random_initial_guess) call controlfunction_randomize(par)
556
557 ! Set the module pointers, so that the calc_point and write_iter_info routines
558 ! can use them.
559 call controlfunction_copy(par_, par)
560 sys_ => sys
561 hm_ => sys%hm
562 td_ => td
563
564 call controlfunction_get_theta(par, x)
565
566 maxiter = oct_iterator_maxiter(iterator)
567 step = oct%direct_step
568 select case (oct%algorithm)
569 case (option__octscheme__oct_nlopt_bobyqa)
571 case (option__octscheme__oct_nlopt_lbfgs)
572 method = minmethod_nlopt_lbfgs
573 end select
574 toldr = oct_iterator_tolerance(iterator)
575
576 call minimize_multidim_nlopt(ierr, method, dim, x, step, toldr, maxiter, opt_control_nlopt_func, minimum, &
577 xl, xu)
578 if (ierr < 1 .or. ierr > 4) then
579 message(1) = "The nlopt minimization procedure did not find convergence, or found an error"
580 write(message(2),'(a,i5)') "Error code =", ierr
581 call messages_info(2, namespace=sys%namespace)
582 end if
583
584 call controlfunction_end(par_)
585 safe_deallocate_a(xl)
586 safe_deallocate_a(xu)
587 safe_deallocate_a(x)
589#endif
590 end subroutine scheme_nlopt
591 ! ---------------------------------------------------------
592
593 end subroutine opt_control_run_legacy
594 ! ---------------------------------------------------------
595
596
597 ! ---------------------------------------------------------
598 subroutine f_zbr98(sys, td, qcpsi, prop_psi, prop_chi, par)
599 type(electrons_t), intent(inout) :: sys
600 type(td_t), intent(inout) :: td
601 type(opt_control_state_t), intent(inout) :: qcpsi
602 type(oct_prop_t), intent(inout) :: prop_psi
603 type(oct_prop_t), intent(inout) :: prop_chi
604 type(controlfunction_t), intent(inout) :: par
605
606 type(states_elec_t) :: chi
607 type(opt_control_state_t) :: qcchi
608 type(controlfunction_t) :: par_chi
609
610 push_sub(f_zbr98)
611
612 call controlfunction_copy(par_chi, par)
613
614 call target_get_state(oct_target, chi)
615 call opt_control_state_init(qcchi, chi, sys%ions)
616 call bwd_step(sys, td, oct_target, par, par_chi, qcchi, prop_chi, prop_psi)
617 call opt_control_state_copy(qcpsi, initial_st)
618 call fwd_step(sys, td, oct_target, par, par_chi, qcpsi, prop_chi, prop_psi)
619
621 call opt_control_state_end(qcchi)
622 call controlfunction_end(par_chi)
623 pop_sub(f_zbr98)
624 end subroutine f_zbr98
625
626
627 ! ---------------------------------------------------------
628 subroutine f_wg05(sys, td, qcpsi, par, prop_psi, prop_chi, j1)
629 type(electrons_t), intent(inout) :: sys
630 type(td_t), intent(inout) :: td
631 type(opt_control_state_t), intent(inout) :: qcpsi
632 type(controlfunction_t), intent(inout) :: par
633 type(oct_prop_t), intent(inout) :: prop_psi
634 type(oct_prop_t), intent(inout) :: prop_chi
635 real(real64), intent(out) :: j1
636
637 real(real64) :: new_penalty
638 type(opt_control_state_t) :: qcchi
639 type(controlfunction_t) :: parp
640
641 push_sub(f_wg05)
642
643 if (oct_iterator_current(iterator) == 0) then
644 call opt_control_state_copy(qcpsi, initial_st)
645 call propagate_forward(sys, td, par, oct_target, qcpsi, prop_psi)
646 j1 = target_j1(oct_target, sys%namespace, sys%gr, sys%kpoints, qcpsi)
647 pop_sub(f_wg05)
648 return
649 end if
650
651 call controlfunction_copy(parp, par)
652
653 call opt_control_state_null(qcchi)
654 call opt_control_state_copy(qcchi, qcpsi)
655 call target_chi(oct_target, sys%namespace, sys%gr, sys%kpoints, qcpsi, qcchi, sys%ions)
656 call bwd_step(sys, td, oct_target, par, parp, qcchi, prop_chi, prop_psi)
657
658 call controlfunction_filter(parp, filter)
659
660 ! recalc field
661 if (oct%mode_fixed_fluence) then
663 call controlfunction_set_alpha(parp, new_penalty)
665 end if
666
667 call controlfunction_copy(par, parp)
669
670 call opt_control_state_copy(qcpsi, initial_st)
671 call propagate_forward(sys, td, par, oct_target, qcpsi, prop_psi)
672
673 j1 = target_j1(oct_target, sys%namespace, sys%gr, sys%kpoints, qcpsi)
674
675 call opt_control_state_end(qcchi)
676 call controlfunction_end(parp)
677 pop_sub(f_wg05)
678 end subroutine f_wg05
679 ! ---------------------------------------------------------
680
681
682 ! ---------------------------------------------------------
683 subroutine f_striter(sys, td, par, j1)
684 type(electrons_t), intent(inout) :: sys
685 type(td_t), intent(inout) :: td
686 type(controlfunction_t), intent(inout) :: par
687 real(real64), intent(out) :: j1
688
689 type(opt_control_state_t) :: qcpsi
690 type(opt_control_state_t) :: qcchi
691 type(controlfunction_t) :: par_chi
692 type(oct_prop_t) :: prop_chi, prop_psi
693
694 push_sub(f_striter)
695
696 call oct_prop_init(prop_chi, sys%namespace, "chi", sys%gr, sys%mc)
697 call oct_prop_init(prop_psi, sys%namespace, "psi", sys%gr, sys%mc)
698
700
701 call controlfunction_copy(par_chi, par)
702
703 ! First, a forward propagation with the input field.
704 call opt_control_state_null(qcpsi)
705 call opt_control_state_copy(qcpsi, initial_st)
706 call propagate_forward(sys, td, par, oct_target, qcpsi, prop_psi)
707
708 ! Check the performance.
709 j1 = target_j1(oct_target, sys%namespace, sys%gr, sys%kpoints, qcpsi, sys%ions)
710
711 ! Set the boundary condition for the backward propagation.
712 call opt_control_state_null(qcchi)
713 call opt_control_state_copy(qcchi, qcpsi)
714 call target_chi(oct_target, sys%namespace, sys%gr, sys%kpoints, qcpsi, qcchi, sys%ions)
715
716 ! Backward propagation, while at the same time finding the output field,
717 ! which is placed at par_chi
718 call bwd_step_2(sys, td, oct_target, par, par_chi, qcchi, prop_chi, prop_psi)
719 !if (oct%mode_fixed_fluence) call controlfunction_set_fluence(par_chi)
720
721 ! Copy par_chi to par
722 call controlfunction_copy(par, par_chi)
723
724 call opt_control_state_end(qcpsi)
725 call opt_control_state_end(qcchi)
726 call controlfunction_end(par_chi)
727 call oct_prop_end(prop_chi)
728 call oct_prop_end(prop_psi)
729
730 pop_sub(f_striter)
731 end subroutine f_striter
732 ! ---------------------------------------------------------
733
734
735 ! ---------------------------------------------------------
736 subroutine f_iter(sys, td, qcpsi, par, prop_psi, prop_chi, j1)
737 type(electrons_t), intent(inout) :: sys
738 type(td_t), intent(inout) :: td
739 type(opt_control_state_t), intent(inout) :: qcpsi
740 type(controlfunction_t), intent(inout) :: par
741 type(oct_prop_t), intent(inout) :: prop_psi
742 type(oct_prop_t), intent(inout) :: prop_chi
743 real(real64), intent(out) :: j1
744
745 type(opt_control_state_t) :: qcchi
746 type(controlfunction_t) :: par_chi
747
748 push_sub(f_iter)
749
750 if (oct_iterator_current(iterator) == 0) then
751 call opt_control_state_copy(qcpsi, initial_st)
752 call propagate_forward(sys, td, par, oct_target, qcpsi, prop_psi)
753 j1 = target_j1(oct_target, sys%namespace, sys%gr, sys%kpoints, qcpsi)
754 pop_sub(f_iter)
755 return
756 end if
757
758 call controlfunction_copy(par_chi, par)
759
760 call opt_control_state_null(qcchi)
761 call opt_control_state_copy(qcchi, qcpsi)
762 call target_chi(oct_target, sys%namespace, sys%gr, sys%kpoints, qcpsi, qcchi, sys%ions)
763 call bwd_step(sys, td, oct_target, par, par_chi, qcchi, prop_chi, prop_psi)
764
765 call opt_control_state_copy(qcpsi, initial_st)
766 call fwd_step(sys, td, oct_target, par, par_chi, qcpsi, prop_chi, prop_psi)
767
768 j1 = target_j1(oct_target, sys%namespace, sys%gr, sys%kpoints, qcpsi)
769
770 call opt_control_state_end(qcchi)
771 call controlfunction_end(par_chi)
772 pop_sub(f_iter)
773 end subroutine f_iter
774 ! ---------------------------------------------------------
775
776#include "opt_control_c_inc.F90"
777#include "check_input_inc.F90"
778#include "finalcheck_inc.F90"
779
780
781end module opt_control_oct_m
782
783!! Local Variables:
784!! mode: f90
785!! coding: utf-8
786!! End:
Define which routines can be seen from the outside.
Definition: loct.F90:147
double sqrt(double __x) __attribute__((__nothrow__
subroutine minimize_multidim(method, dim, x, step, line_tol, tolgrad, toldr, maxiter, f, write_iter_info, minimum, ierr)
Definition: minimizer.F90:403
This module contains the definition of the data type that holds a "control function" used for OCT run...
subroutine, public controlfunction_filter(par, filter)
subroutine, public controlfunction_apply_envelope(cp)
subroutine, public controlfunction_set_fluence(par)
subroutine, public controlfunction_bounds(par, lower_bounds, upper_bounds)
subroutine, public controlfunction_to_realtime(par)
integer pure function, public controlfunction_dof(par)
real(real64) function, public controlfunction_fluence(par)
subroutine, public controlfunction_set_alpha(par, alpha)
subroutine, public controlfunction_write(filename, cp, namespace)
subroutine, public controlfunction_prepare_initial(par)
"Prepares" the initial guess control field: maybe it has to be normalized to a certain fluence,...
real(real64) pure function, public controlfunction_w0(par)
subroutine, public controlfunction_end(cp)
subroutine, public controlfunction_get_theta(par, theta)
real(real64) function, public controlfunction_j2(par)
subroutine, public controlfunction_randomize(par)
real(real64) pure function, public controlfunction_alpha(par, ipar)
real(real64) pure function, public controlfunction_targetfluence()
subroutine, public controlfunction_copy(cp_out, cp_in)
subroutine, public controlfunction_set_rep(par)
Transforms the control function to frequency space, if this is the space in which the functions are d...
subroutine, public controlfunction_mod_close()
subroutine, public controlfunction_to_h(cp, ext_partners)
subroutine, public controlfunction_init(cp, dt, ntiter)
Before using an controlfunction_t variable, it needs to be initialized, either by calling controlfunc...
subroutine, public controlfunction_mod_init(ext_partners, namespace, dt, max_iter, mode_fixed_fluence)
Initializes the module, should be the first subroutine to be called (the last one should be controlfu...
subroutine, public controlfunction_set(cp, ext_partners)
The external fields defined in epot_t "ep" are transferred to the control functions described in "cp"...
type(lasers_t) function, pointer, public list_get_lasers(partners)
subroutine, public filter_init(steps, namespace, dt, filter)
Definition: filter.F90:154
subroutine, public filter_write(filter, namespace)
Definition: filter.F90:333
subroutine, public filter_end(filter)
Definition: filter.F90:367
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:186
character(len= *), parameter, public oct_dir
Definition: global.F90:259
This module implements the underlying real-space grid.
Definition: grid.F90:117
subroutine, public initial_state_init(sys, qcstate)
Definition: initst.F90:154
Definition: io.F90:114
subroutine, public io_mkdir(fname, namespace, parents)
Definition: io.F90:311
A module to handle KS potential, without the external potential.
integer, parameter, public independent_particles
subroutine, public laser_write_info(lasers, namespace, dt, max_iter, iunit)
Definition: lasers.F90:930
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:115
subroutine, public messages_print_with_emphasis(msg, iunit, namespace)
Definition: messages.F90:920
subroutine, public messages_not_implemented(feature, namespace)
Definition: messages.F90:1113
character(len=512), private msg
Definition: messages.F90:165
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:160
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:414
subroutine, public messages_input_error(namespace, var, details, row, column)
Definition: messages.F90:713
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1085
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:616
integer minmethod_fr_cg
Definition: minimizer.F90:134
integer minmethod_nmsimplex
Definition: minimizer.F90:134
integer minmethod_bfgs2
Definition: minimizer.F90:134
integer minmethod_nlopt_lbfgs
Definition: minimizer.F90:134
integer minmethod_nlopt_bobyqa
Definition: minimizer.F90:134
This module implements the basic mulsisystem class, a container system for other systems.
This module contains the definition of the oct_t data type, which contains some of the basic informat...
subroutine, public oct_read_inp(oct, namespace)
Reads, from the inp file, some global information about how the QOCT run should be.
logical function, public iteration_manager(namespace, j1, par, par_prev, iterator)
subroutine, public oct_iterator_init(iterator, namespace, par)
integer pure function, public oct_iterator_maxiter(iterator)
subroutine, public oct_iterator_end(iterator, namespace)
integer pure function, public oct_iterator_current(iterator)
subroutine, public iteration_manager_direct(j, par, iterator, sys, dx)
real(real64) pure function, public oct_iterator_tolerance(iterator)
This module contains the main procedure ("opt_control_run") that is used when optimal control runs ar...
subroutine, public opt_control_cg_calc(n, x, f, getgrad, df)
subroutine check_faulty_runmodes(sys, tr)
subroutine, public opt_control_run(system)
subroutine oct_finalcheck(sys, td)
subroutine, public opt_control_cg_write_info(iter, n, val, maxdx, maxdf, x)
interface is required by its being passed as dummy routine to minimize_multidim
subroutine f_striter(sys, td, par, j1)
subroutine, public opt_control_direct_message_info(iter, n, val, maxdx, x)
subroutine f_zbr98(sys, td, qcpsi, prop_psi, prop_chi, par)
subroutine f_iter(sys, td, qcpsi, par, prop_psi, prop_chi, j1)
subroutine opt_control_run_legacy(sys)
This is the main procedure for all types of optimal control runs. It is called from the "run" procedu...
subroutine f_wg05(sys, td, qcpsi, par, prop_psi, prop_chi, j1)
subroutine, public opt_control_function_forward(x, f)
subroutine opt_control_nlopt_func(val, n, x, grad, need_gradient, f_data)
subroutine, public opt_control_direct_calc(n, x, f)
No intents here is unfortunately required because this will be passed to newuoa routines as a dummy f...
This module holds the "opt_control_state_t" datatype, which contains a quantum-classical state.
subroutine, public opt_control_state_end(ocs)
subroutine, public opt_control_state_null(ocs)
subroutine, public opt_control_state_copy(ocsout, ocsin)
subroutine, public opt_control_state_init(ocs, qstate, ions)
subroutine, public opt_control_get_qs(qstate, ocs)
this module contains the output system
Definition: output.F90:115
subroutine, public output_states(outp, namespace, space, dir, st, gr, ions, hm, iter)
Definition: output.F90:1888
subroutine, public propagation_mod_init(niter, eta, delta, number_checkpoints, zbr98, gradients)
This subroutine must be called before any QOCT propagations are done. It simply stores in the module ...
subroutine, public oct_prop_end(prop)
subroutine, public bwd_step(sys, td, tg, par, par_chi, qcchi, prop_chi, prop_psi)
subroutine, public propagate_forward(sys, td, par, tg, qcpsi, prop, write_iter)
subroutine, public fwd_step(sys, td, tg, par, par_chi, qcpsi, prop_chi, prop_psi)
subroutine, public oct_prop_init(prop, namespace, dirname, mesh, mc)
subroutine, public bwd_step_2(sys, td, tg, par, par_chi, qcchi, prop_chi, prop_psi)
subroutine, public propagator_elec_set_scf_prop(tr, threshold)
logical function, public clean_stop(comm)
returns true if a file named stop exists
Definition: restart.F90:276
This module handles spin dimensions of the states and the k-point distribution.
subroutine, public states_elec_end(st)
finalize the states_elec_t object
real(real64) function, public target_j1(tg, namespace, gr, kpoints, qcpsi, ions)
Calculates the J1 functional, i.e.: in the time-independent case, or else in the time-dependent cas...
Definition: target.F90:578
subroutine, public target_init(gr, kpoints, namespace, space, ions, qcs, td, w0, tg, oct, ep, mc)
The target is initialized, mainly by reading from the inp file.
Definition: target.F90:217
subroutine, public target_get_state(tg, st)
This just copies the states_elec_t variable present in target, into st.
Definition: target.F90:203
subroutine, public target_output(tg, namespace, space, gr, dir, ions, hm, outp)
Definition: target.F90:396
subroutine, public target_chi(tg, namespace, gr, kpoints, qcpsi_in, qcchi_out, ions)
Calculate .
Definition: target.F90:630
subroutine, public target_end(tg, oct)
Definition: target.F90:363
Definition: td.F90:114
subroutine, public td_end(td)
Definition: td.F90:603
subroutine, public td_init(td, namespace, space, gr, ions, st, ks, hm, ext_partners, outp)
Definition: td.F90:247
subroutine scheme_zbr98()
subroutine scheme_cg()
subroutine scheme_direct()
subroutine scheme_straight_iteration()
subroutine scheme_mt03()
subroutine scheme_wg05()
subroutine scheme_nlopt()
This is the data type used to hold a control function.
Class describing the electron system.
Definition: electrons.F90:218
Container class for lists of system_oct_m::system_t.
This is the datatype that contains the objects that are propagated: in principle this could be both t...
The states_elec_t class contains all electronic wave functions.