(git:e7e05ae)
optimize_embedding_potential.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 MODULE optimize_embedding_potential
9 
10  USE atomic_kind_types, ONLY: atomic_kind_type,&
11  get_atomic_kind,&
12  get_atomic_kind_set
13  USE cell_types, ONLY: cell_type
14  USE cp_blacs_env, ONLY: cp_blacs_env_create,&
15  cp_blacs_env_release,&
16  cp_blacs_env_type
17  USE cp_control_types, ONLY: dft_control_type
18  USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
19  dbcsr_deallocate_matrix_set
20  USE cp_files, ONLY: close_file,&
21  open_file
22  USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
23  cp_fm_scale,&
24  cp_fm_scale_and_add,&
25  cp_fm_trace
26  USE cp_fm_diag, ONLY: choose_eigv_solver
27  USE cp_fm_struct, ONLY: cp_fm_struct_create,&
28  cp_fm_struct_release,&
29  cp_fm_struct_type
30  USE cp_fm_types, ONLY: &
31  cp_fm_copy_general, cp_fm_create, cp_fm_get_element, cp_fm_get_info, cp_fm_release, &
32  cp_fm_set_all, cp_fm_to_fm, cp_fm_to_fm_submat, cp_fm_type, cp_fm_write_unformatted
33  USE cp_log_handling, ONLY: cp_get_default_logger,&
34  cp_logger_type
35  USE cp_output_handling, ONLY: cp_p_file,&
36  cp_print_key_finished_output,&
37  cp_print_key_should_output,&
38  cp_print_key_unit_nr
39  USE cp_realspace_grid_cube, ONLY: cp_cube_to_pw,&
40  cp_pw_to_cube,&
41  cp_pw_to_simple_volumetric
42  USE dbcsr_api, ONLY: dbcsr_p_type
43  USE embed_types, ONLY: opt_embed_pot_type
44  USE force_env_types, ONLY: force_env_type
45  USE input_constants, ONLY: &
46  embed_diff, embed_fa, embed_grid_angstrom, embed_grid_bohr, embed_level_shift, embed_none, &
47  embed_quasi_newton, embed_resp, embed_steep_desc
48  USE input_section_types, ONLY: section_get_ival,&
49  section_get_ivals,&
50  section_get_rval,&
51  section_vals_get_subs_vals,&
52  section_vals_type,&
53  section_vals_val_get
54  USE kinds, ONLY: default_path_length,&
55  dp
56  USE lri_environment_types, ONLY: lri_kind_type
57  USE mathconstants, ONLY: pi
58  USE message_passing, ONLY: mp_para_env_type
59  USE mixed_environment_utils, ONLY: get_subsys_map_index
60  USE parallel_gemm_api, ONLY: parallel_gemm
61  USE particle_list_types, ONLY: particle_list_type
62  USE particle_types, ONLY: particle_type
63  USE pw_env_types, ONLY: pw_env_get,&
64  pw_env_type
65  USE pw_methods, ONLY: &
66  pw_axpy, pw_copy, pw_derive, pw_dr2, pw_integral_ab, pw_integrate_function, pw_scale, &
67  pw_transfer, pw_zero
68  USE pw_poisson_methods, ONLY: pw_poisson_solve
69  USE pw_poisson_types, ONLY: pw_poisson_type
70  USE pw_pool_types, ONLY: pw_pool_type
71  USE pw_types, ONLY: pw_c1d_gs_type,&
72  pw_r3d_rs_type
73  USE qs_collocate_density, ONLY: calculate_rho_resp_all,&
74  calculate_wavefunction
75  USE qs_environment_types, ONLY: get_qs_env,&
76  qs_environment_type,&
77  set_qs_env
78  USE qs_integrate_potential_single, ONLY: integrate_v_rspace_one_center
79  USE qs_kind_types, ONLY: get_qs_kind,&
80  qs_kind_type
81  USE qs_kinetic, ONLY: build_kinetic_matrix
82  USE qs_ks_types, ONLY: qs_ks_env_type
83  USE qs_mo_types, ONLY: get_mo_set,&
84  mo_set_type
85  USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
86  USE qs_rho_types, ONLY: qs_rho_get,&
87  qs_rho_type
88  USE qs_subsys_types, ONLY: qs_subsys_get,&
89  qs_subsys_type
90  USE xc, ONLY: smooth_cutoff
91  USE xc_rho_cflags_types, ONLY: xc_rho_cflags_setall,&
92  xc_rho_cflags_type
93  USE xc_rho_set_types, ONLY: xc_rho_set_create,&
94  xc_rho_set_release,&
95  xc_rho_set_type,&
96  xc_rho_set_update
97 #include "./base/base_uses.f90"
98 
99  IMPLICIT NONE
100 
101  PRIVATE
102 
103  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'optimize_embedding_potential'
104 
105  PUBLIC :: prepare_embed_opt, init_embed_pot, release_opt_embed, calculate_embed_pot_grad, &
106  opt_embed_step, print_rho_diff, step_control, max_dens_diff, print_emb_opt_info, &
107  conv_check_embed, make_subsys_embed_pot, print_embed_restart, find_aux_dimen, &
108  read_embed_pot, understand_spin_states, given_embed_pot, print_rho_spin_diff, &
109  print_pot_simple_grid, get_prev_density, get_max_subsys_diff, coulomb_guess
110 
111 CONTAINS
112 
113 ! **************************************************************************************************
114 !> \brief Find out whether we need to swap alpha- and beta- spind densities in the second subsystem
115 !> \brief It's only needed because by default alpha-spins go first in a subsystem.
116 !> \brief By swapping we impose the constraint:
117 !> \brief rho_1(alpha) + rho_2(alpha) = rho_total(alpha)
118 !> \brief rho_1(beta) + rho_2(beta) = rho_total(beta)
119 !> \param force_env ...
120 !> \param ref_subsys_number ...
121 !> \param change_spin ...
122 !> \param open_shell_embed ...
123 !> \param all_nspins ...
124 !> \return ...
125 !> \author Vladimir Rybkin
126 ! **************************************************************************************************
127  SUBROUTINE understand_spin_states(force_env, ref_subsys_number, change_spin, open_shell_embed, all_nspins)
128  TYPE(force_env_type), POINTER :: force_env
129  INTEGER :: ref_subsys_number
130  LOGICAL :: change_spin, open_shell_embed
131  INTEGER, ALLOCATABLE, DIMENSION(:) :: all_nspins
132 
133  INTEGER :: i_force_eval, nspins, sub_spin_1, &
134  sub_spin_2, total_spin
135  INTEGER, DIMENSION(2) :: nelectron_spin
136  INTEGER, DIMENSION(2, 3) :: all_spins
137  TYPE(dft_control_type), POINTER :: dft_control
138 
139  change_spin = .false.
140  open_shell_embed = .false.
141  ALLOCATE (all_nspins(ref_subsys_number))
142  IF (ref_subsys_number .EQ. 3) THEN
143  all_spins = 0
144  DO i_force_eval = 1, ref_subsys_number
145  CALL get_qs_env(qs_env=force_env%sub_force_env(i_force_eval)%force_env%qs_env, &
146  nelectron_spin=nelectron_spin, dft_control=dft_control)
147  all_spins(:, i_force_eval) = nelectron_spin
148  nspins = dft_control%nspins
149  all_nspins(i_force_eval) = nspins
150  END DO
151 
152  ! Find out whether we need a spin-dependend embedding potential
153  IF (.NOT. ((all_nspins(1) .EQ. 1) .AND. (all_nspins(2) .EQ. 1) .AND. (all_nspins(3) .EQ. 1))) &
154  open_shell_embed = .true.
155 
156  ! If it's open shell, we need to check spin states
157  IF (open_shell_embed) THEN
158 
159  IF (all_nspins(3) .EQ. 1) THEN
160  total_spin = 0
161  ELSE
162  total_spin = all_spins(1, 3) - all_spins(2, 3)
163  END IF
164  IF (all_nspins(1) .EQ. 1) THEN
165  sub_spin_1 = 0
166  ELSE
167  sub_spin_1 = all_spins(1, 1) - all_spins(2, 1)
168  END IF
169  IF (all_nspins(2) .EQ. 1) THEN
170  sub_spin_2 = 0
171  ELSE
172  sub_spin_2 = all_spins(1, 2) - all_spins(2, 2)
173  END IF
174  IF ((sub_spin_1 + sub_spin_2) .EQ. total_spin) THEN
175  change_spin = .false.
176  ELSE
177  IF (abs(sub_spin_1 - sub_spin_2) .EQ. total_spin) THEN
178  change_spin = .true.
179  ELSE
180  cpabort("Spin states of subsystems are not compatible.")
181  END IF
182  END IF
183 
184  END IF ! not open_shell
185  ELSE
186  cpabort("Reference subsystem must be the third FORCE_EVAL.")
187  END IF
188 
189  END SUBROUTINE understand_spin_states
190 
191 ! **************************************************************************************************
192 !> \brief ...
193 !> \param qs_env ...
194 !> \param embed_pot ...
195 !> \param add_const_pot ...
196 !> \param Fermi_Amaldi ...
197 !> \param const_pot ...
198 !> \param open_shell_embed ...
199 !> \param spin_embed_pot ...
200 !> \param pot_diff ...
201 !> \param Coulomb_guess ...
202 !> \param grid_opt ...
203 ! **************************************************************************************************
204  SUBROUTINE init_embed_pot(qs_env, embed_pot, add_const_pot, Fermi_Amaldi, const_pot, open_shell_embed, &
205  spin_embed_pot, pot_diff, Coulomb_guess, grid_opt)
206  TYPE(qs_environment_type), POINTER :: qs_env
207  TYPE(pw_r3d_rs_type), POINTER :: embed_pot
208  LOGICAL :: add_const_pot, fermi_amaldi
209  TYPE(pw_r3d_rs_type), POINTER :: const_pot
210  LOGICAL :: open_shell_embed
211  TYPE(pw_r3d_rs_type), POINTER :: spin_embed_pot, pot_diff
212  LOGICAL :: coulomb_guess, grid_opt
213 
214  INTEGER :: nelectrons
215  INTEGER, DIMENSION(2) :: nelectron_spin
216  REAL(kind=dp) :: factor
217  TYPE(pw_env_type), POINTER :: pw_env
218  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
219  TYPE(pw_r3d_rs_type), POINTER :: v_hartree_r_space
220 
221  ! Extract plane waves environment
222  CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, &
223  nelectron_spin=nelectron_spin, &
224  v_hartree_rspace=v_hartree_r_space)
225 
226  ! Prepare plane-waves pool
227  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
228 
229  ! Create embedding potential and set to zero
230  NULLIFY (embed_pot)
231  ALLOCATE (embed_pot)
232  CALL auxbas_pw_pool%create_pw(embed_pot)
233  CALL pw_zero(embed_pot)
234 
235  ! Spin embedding potential if asked
236  IF (open_shell_embed) THEN
237  NULLIFY (spin_embed_pot)
238  ALLOCATE (spin_embed_pot)
239  CALL auxbas_pw_pool%create_pw(spin_embed_pot)
240  CALL pw_zero(spin_embed_pot)
241  END IF
242 
243  ! Coulomb guess/constant potential
244  IF (coulomb_guess) THEN
245  NULLIFY (pot_diff)
246  ALLOCATE (pot_diff)
247  CALL auxbas_pw_pool%create_pw(pot_diff)
248  CALL pw_zero(pot_diff)
249  END IF
250 
251  ! Initialize constant part of the embedding potential
252  IF (add_const_pot .AND. (.NOT. grid_opt)) THEN
253  ! Now the constant potential is the Coulomb one
254  NULLIFY (const_pot)
255  ALLOCATE (const_pot)
256  CALL auxbas_pw_pool%create_pw(const_pot)
257  CALL pw_zero(const_pot)
258  END IF
259 
260  ! Add Fermi-Amaldi potential if requested
261  IF (fermi_amaldi) THEN
262 
263  ! Extract Hartree potential
264  NULLIFY (v_hartree_r_space)
265  CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, &
266  v_hartree_rspace=v_hartree_r_space)
267  CALL pw_copy(v_hartree_r_space, embed_pot)
268 
269  ! Calculate the number of electrons
270  nelectrons = nelectron_spin(1) + nelectron_spin(2)
271  factor = (real(nelectrons, dp) - 1.0_dp)/(real(nelectrons, dp))
272 
273  ! Scale the Hartree potential to get Fermi-Amaldi
274  CALL pw_scale(embed_pot, a=factor)
275 
276  ! Copy Fermi-Amaldi to embedding potential for basis-based optimization
277  IF (.NOT. grid_opt) CALL pw_copy(embed_pot, embed_pot)
278 
279  END IF
280 
281  END SUBROUTINE init_embed_pot
282 
283 ! **************************************************************************************************
284 !> \brief Creates and allocates objects for optimization of embedding potential
285 !> \param qs_env ...
286 !> \param opt_embed ...
287 !> \param opt_embed_section ...
288 !> \author Vladimir Rybkin
289 ! **************************************************************************************************
290  SUBROUTINE prepare_embed_opt(qs_env, opt_embed, opt_embed_section)
291  TYPE(qs_environment_type), POINTER :: qs_env
292  TYPE(opt_embed_pot_type) :: opt_embed
293  TYPE(section_vals_type), POINTER :: opt_embed_section
294 
295  INTEGER :: diff_size, i_dens, size_prev_dens
296  TYPE(cp_blacs_env_type), POINTER :: blacs_env
297  TYPE(cp_fm_struct_type), POINTER :: fm_struct
298  TYPE(mp_para_env_type), POINTER :: para_env
299  TYPE(pw_env_type), POINTER :: pw_env
300  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
301 
302  !TYPE(pw_env_type), POINTER :: pw_env
303  !TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
304 
305  ! First, read the input
306 
307  CALL read_opt_embed_section(opt_embed, opt_embed_section)
308 
309  ! All these are needed for optimization in a finite Gaussian basis
310  IF (.NOT. opt_embed%grid_opt) THEN
311  ! Create blacs environment
312  CALL get_qs_env(qs_env=qs_env, &
313  para_env=para_env)
314  NULLIFY (blacs_env)
315  CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env)
316 
317  ! Reveal the dimension of the RI basis
318  CALL find_aux_dimen(qs_env, opt_embed%dimen_aux)
319 
320  ! Prepare the object for integrals
321  CALL make_lri_object(qs_env, opt_embed%lri)
322 
323  ! In case if spin embedding potential has to be optimized,
324  ! the dimension of variational space is two times larger
325  IF (opt_embed%open_shell_embed) THEN
326  opt_embed%dimen_var_aux = 2*opt_embed%dimen_aux
327  ELSE
328  opt_embed%dimen_var_aux = opt_embed%dimen_aux
329  END IF
330 
331  ! Allocate expansion coefficients and gradient
332  NULLIFY (opt_embed%embed_pot_grad, opt_embed%embed_pot_coef, opt_embed%step, fm_struct)
333 
334  NULLIFY (opt_embed%prev_embed_pot_grad, opt_embed%prev_embed_pot_coef, opt_embed%prev_step)
335  CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, &
336  nrow_global=opt_embed%dimen_var_aux, ncol_global=1)
337  ALLOCATE (opt_embed%embed_pot_grad, opt_embed%embed_pot_coef, &
338  opt_embed%prev_embed_pot_grad, opt_embed%prev_embed_pot_coef, &
339  opt_embed%step, opt_embed%prev_step)
340  CALL cp_fm_create(opt_embed%embed_pot_grad, fm_struct, name="pot_grad")
341  CALL cp_fm_create(opt_embed%embed_pot_coef, fm_struct, name="pot_coef")
342  CALL cp_fm_create(opt_embed%prev_embed_pot_grad, fm_struct, name="prev_pot_grad")
343  CALL cp_fm_create(opt_embed%prev_embed_pot_coef, fm_struct, name="prev_pot_coef")
344  CALL cp_fm_create(opt_embed%step, fm_struct, name="step")
345  CALL cp_fm_create(opt_embed%prev_step, fm_struct, name="prev_step")
346 
347  CALL cp_fm_struct_release(fm_struct)
348  CALL cp_fm_set_all(opt_embed%embed_pot_grad, 0.0_dp)
349  CALL cp_fm_set_all(opt_embed%prev_embed_pot_grad, 0.0_dp)
350  CALL cp_fm_set_all(opt_embed%embed_pot_coef, 0.0_dp)
351  CALL cp_fm_set_all(opt_embed%prev_embed_pot_coef, 0.0_dp)
352  CALL cp_fm_set_all(opt_embed%step, 0.0_dp)
353 
354  CALL cp_fm_set_all(opt_embed%prev_step, 0.0_dp)
355 
356  ! Allocate Hessian
357  NULLIFY (opt_embed%embed_pot_hess, opt_embed%prev_embed_pot_hess, fm_struct)
358  CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, &
359  nrow_global=opt_embed%dimen_var_aux, ncol_global=opt_embed%dimen_var_aux)
360  ALLOCATE (opt_embed%embed_pot_hess, opt_embed%prev_embed_pot_hess)
361  CALL cp_fm_create(opt_embed%embed_pot_hess, fm_struct, name="pot_Hess")
362  CALL cp_fm_create(opt_embed%prev_embed_pot_hess, fm_struct, name="prev_pot_Hess")
363  CALL cp_fm_struct_release(fm_struct)
364 
365  ! Special structure for the kinetic energy matrix
366  NULLIFY (fm_struct, opt_embed%kinetic_mat)
367  CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, &
368  nrow_global=opt_embed%dimen_aux, ncol_global=opt_embed%dimen_aux)
369  ALLOCATE (opt_embed%kinetic_mat)
370  CALL cp_fm_create(opt_embed%kinetic_mat, fm_struct, name="kinetic_mat")
371  CALL cp_fm_struct_release(fm_struct)
372  CALL cp_fm_set_all(opt_embed%kinetic_mat, 0.0_dp)
373 
374  ! Hessian is set as a unit matrix
375  CALL cp_fm_set_all(opt_embed%embed_pot_hess, 0.0_dp, -1.0_dp)
376  CALL cp_fm_set_all(opt_embed%prev_embed_pot_hess, 0.0_dp, -1.0_dp)
377 
378  ! Release blacs environment
379  CALL cp_blacs_env_release(blacs_env)
380 
381  END IF
382 
383  CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
384  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
385  NULLIFY (opt_embed%prev_subsys_dens)
386  size_prev_dens = sum(opt_embed%all_nspins(1:(SIZE(opt_embed%all_nspins) - 1)))
387  ALLOCATE (opt_embed%prev_subsys_dens(size_prev_dens))
388  DO i_dens = 1, size_prev_dens
389  CALL auxbas_pw_pool%create_pw(opt_embed%prev_subsys_dens(i_dens))
390  CALL pw_zero(opt_embed%prev_subsys_dens(i_dens))
391  END DO
392  ALLOCATE (opt_embed%max_subsys_dens_diff(size_prev_dens))
393 
394  ! Array to store functional values
395  ALLOCATE (opt_embed%w_func(opt_embed%n_iter))
396  opt_embed%w_func = 0.0_dp
397 
398  ! Allocate max_diff and int_diff
399  diff_size = 1
400  IF (opt_embed%open_shell_embed) diff_size = 2
401  ALLOCATE (opt_embed%max_diff(diff_size))
402  ALLOCATE (opt_embed%int_diff(diff_size))
403  ALLOCATE (opt_embed%int_diff_square(diff_size))
404 
405  ! FAB update
406  IF (opt_embed%fab) THEN
407  NULLIFY (opt_embed%prev_embed_pot)
408  ALLOCATE (opt_embed%prev_embed_pot)
409  CALL auxbas_pw_pool%create_pw(opt_embed%prev_embed_pot)
410  CALL pw_zero(opt_embed%prev_embed_pot)
411  IF (opt_embed%open_shell_embed) THEN
412  NULLIFY (opt_embed%prev_spin_embed_pot)
413  ALLOCATE (opt_embed%prev_spin_embed_pot)
414  CALL auxbas_pw_pool%create_pw(opt_embed%prev_spin_embed_pot)
415  CALL pw_zero(opt_embed%prev_spin_embed_pot)
416  END IF
417  END IF
418 
419  ! Set allowed energy decrease parameter
420  opt_embed%allowed_decrease = 0.0001_dp
421 
422  ! Regularization contribution is set to zero
423  opt_embed%reg_term = 0.0_dp
424 
425  ! Step is accepted in the beginning
426  opt_embed%accept_step = .true.
427  opt_embed%newton_step = .false.
428  opt_embed%last_accepted = 1
429 
430  ! Set maximum and minimum trust radii
431  opt_embed%max_trad = opt_embed%trust_rad*7.900_dp
432  opt_embed%min_trad = opt_embed%trust_rad*0.125*0.065_dp
433 
434  END SUBROUTINE prepare_embed_opt
435 
436 ! **************************************************************************************************
437 !> \brief ...
438 !> \param opt_embed ...
439 !> \param opt_embed_section ...
440 ! **************************************************************************************************
441  SUBROUTINE read_opt_embed_section(opt_embed, opt_embed_section)
442  TYPE(opt_embed_pot_type) :: opt_embed
443  TYPE(section_vals_type), POINTER :: opt_embed_section
444 
445  INTEGER :: embed_guess, embed_optimizer
446 
447  ! Read keywords
448  CALL section_vals_val_get(opt_embed_section, "REG_LAMBDA", &
449  r_val=opt_embed%lambda)
450 
451  CALL section_vals_val_get(opt_embed_section, "N_ITER", &
452  i_val=opt_embed%n_iter)
453 
454  CALL section_vals_val_get(opt_embed_section, "TRUST_RAD", &
455  r_val=opt_embed%trust_rad)
456 
457  CALL section_vals_val_get(opt_embed_section, "DENS_CONV_MAX", &
458  r_val=opt_embed%conv_max)
459 
460  CALL section_vals_val_get(opt_embed_section, "DENS_CONV_INT", &
461  r_val=opt_embed%conv_int)
462 
463  CALL section_vals_val_get(opt_embed_section, "SPIN_DENS_CONV_MAX", &
464  r_val=opt_embed%conv_max_spin)
465 
466  CALL section_vals_val_get(opt_embed_section, "SPIN_DENS_CONV_INT", &
467  r_val=opt_embed%conv_int_spin)
468 
469  CALL section_vals_val_get(opt_embed_section, "CHARGE_DISTR_WIDTH", &
470  r_val=opt_embed%eta)
471 
472  CALL section_vals_val_get(opt_embed_section, "READ_EMBED_POT", &
473  l_val=opt_embed%read_embed_pot)
474 
475  CALL section_vals_val_get(opt_embed_section, "READ_EMBED_POT_CUBE", &
476  l_val=opt_embed%read_embed_pot_cube)
477 
478  CALL section_vals_val_get(opt_embed_section, "GRID_OPT", &
479  l_val=opt_embed%grid_opt)
480 
481  CALL section_vals_val_get(opt_embed_section, "LEEUWEN-BAERENDS", &
482  l_val=opt_embed%leeuwen)
483 
484  CALL section_vals_val_get(opt_embed_section, "FAB", &
485  l_val=opt_embed%fab)
486 
487  CALL section_vals_val_get(opt_embed_section, "VW_CUTOFF", &
488  r_val=opt_embed%vw_cutoff)
489 
490  CALL section_vals_val_get(opt_embed_section, "VW_SMOOTH_CUT_RANGE", &
491  r_val=opt_embed%vw_smooth_cutoff_range)
492 
493  CALL section_vals_val_get(opt_embed_section, "OPTIMIZER", i_val=embed_optimizer)
494  SELECT CASE (embed_optimizer)
495  CASE (embed_steep_desc)
496  opt_embed%steep_desc = .true.
497  CASE (embed_quasi_newton)
498  opt_embed%steep_desc = .false.
499  opt_embed%level_shift = .false.
500  CASE (embed_level_shift)
501  opt_embed%steep_desc = .false.
502  opt_embed%level_shift = .true.
503  CASE DEFAULT
504  opt_embed%steep_desc = .true.
505  END SELECT
506 
507  CALL section_vals_val_get(opt_embed_section, "POT_GUESS", i_val=embed_guess)
508  SELECT CASE (embed_guess)
509  CASE (embed_none)
510  opt_embed%add_const_pot = .false.
511  opt_embed%Fermi_Amaldi = .false.
512  opt_embed%Coulomb_guess = .false.
513  opt_embed%diff_guess = .false.
514  CASE (embed_diff)
515  opt_embed%add_const_pot = .true.
516  opt_embed%Fermi_Amaldi = .false.
517  opt_embed%Coulomb_guess = .false.
518  opt_embed%diff_guess = .true.
519  CASE (embed_fa)
520  opt_embed%add_const_pot = .true.
521  opt_embed%Fermi_Amaldi = .true.
522  opt_embed%Coulomb_guess = .false.
523  opt_embed%diff_guess = .false.
524  CASE (embed_resp)
525  opt_embed%add_const_pot = .true.
526  opt_embed%Fermi_Amaldi = .true.
527  opt_embed%Coulomb_guess = .true.
528  opt_embed%diff_guess = .false.
529  CASE DEFAULT
530  opt_embed%add_const_pot = .false.
531  opt_embed%Fermi_Amaldi = .false.
532  opt_embed%Coulomb_guess = .false.
533  opt_embed%diff_guess = .false.
534  END SELECT
535 
536  END SUBROUTINE read_opt_embed_section
537 
538 ! **************************************************************************************************
539 !> \brief Find the dimension of the auxiliary basis for the expansion of the embedding potential
540 !> \param qs_env ...
541 !> \param dimen_aux ...
542 ! **************************************************************************************************
543  SUBROUTINE find_aux_dimen(qs_env, dimen_aux)
544  TYPE(qs_environment_type), POINTER :: qs_env
545  INTEGER :: dimen_aux
546 
547  INTEGER :: iatom, ikind, nsgf
548  INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
549  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
550  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
551  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
552 
553  ! First, reveal the dimension of the RI basis
554  CALL get_qs_env(qs_env=qs_env, &
555  particle_set=particle_set, &
556  qs_kind_set=qs_kind_set, &
557  atomic_kind_set=atomic_kind_set)
558 
559  CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
560 
561  dimen_aux = 0
562  DO iatom = 1, SIZE(particle_set)
563  ikind = kind_of(iatom)
564  CALL get_qs_kind(qs_kind=qs_kind_set(ikind), nsgf=nsgf, basis_type="RI_AUX")
565  dimen_aux = dimen_aux + nsgf
566  END DO
567 
568  END SUBROUTINE find_aux_dimen
569 
570 ! **************************************************************************************************
571 !> \brief Prepare the lri_kind_type object for integrals between density and aux. basis functions
572 !> \param qs_env ...
573 !> \param lri ...
574 ! **************************************************************************************************
575  SUBROUTINE make_lri_object(qs_env, lri)
576  TYPE(qs_environment_type), POINTER :: qs_env
577  TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri
578 
579  INTEGER :: ikind, natom, nkind, nsgf
580  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
581  TYPE(atomic_kind_type), POINTER :: atomic_kind
582  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
583 
584  NULLIFY (atomic_kind, lri)
585  CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, &
586  qs_kind_set=qs_kind_set)
587  nkind = SIZE(atomic_kind_set)
588 
589  ALLOCATE (lri(nkind))
590  ! Here we need only v_int and acoef (the latter as dummies)
591  DO ikind = 1, nkind
592  NULLIFY (lri(ikind)%acoef)
593  NULLIFY (lri(ikind)%v_int)
594  atomic_kind => atomic_kind_set(ikind)
595  CALL get_atomic_kind(atomic_kind=atomic_kind, natom=natom)
596  CALL get_qs_kind(qs_kind=qs_kind_set(ikind), nsgf=nsgf, basis_type="RI_AUX")
597  ALLOCATE (lri(ikind)%acoef(natom, nsgf))
598  lri(ikind)%acoef = 0._dp
599  ALLOCATE (lri(ikind)%v_int(natom, nsgf))
600  lri(ikind)%v_int = 0._dp
601  END DO
602 
603  END SUBROUTINE make_lri_object
604 
605 ! **************************************************************************************************
606 !> \brief Read the external embedding potential, not to be optimized
607 !> \param qs_env ...
608 ! **************************************************************************************************
609  SUBROUTINE given_embed_pot(qs_env)
610  TYPE(qs_environment_type), POINTER :: qs_env
611 
612  LOGICAL :: open_shell_embed
613  TYPE(dft_control_type), POINTER :: dft_control
614  TYPE(pw_env_type), POINTER :: pw_env
615  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool_subsys
616  TYPE(pw_r3d_rs_type), POINTER :: embed_pot, spin_embed_pot
617  TYPE(section_vals_type), POINTER :: input, qs_section
618 
619  qs_env%given_embed_pot = .true.
620  NULLIFY (input, dft_control, embed_pot, spin_embed_pot, embed_pot, spin_embed_pot, &
621  qs_section)
622  CALL get_qs_env(qs_env=qs_env, &
623  input=input, &
624  dft_control=dft_control, &
625  pw_env=pw_env)
626  qs_section => section_vals_get_subs_vals(input, "DFT%QS")
627  open_shell_embed = .false.
628  IF (dft_control%nspins .EQ. 2) open_shell_embed = .true.
629 
630  ! Prepare plane-waves pool
631  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool_subsys)
632 
633  ! Create embedding potential
634  !CALL get_qs_env(qs_env=qs_env, &
635  ! embed_pot=embed_pot)
636  ALLOCATE (embed_pot)
637  CALL auxbas_pw_pool_subsys%create_pw(embed_pot)
638  IF (open_shell_embed) THEN
639  ! Create spin embedding potential
640  ALLOCATE (spin_embed_pot)
641  CALL auxbas_pw_pool_subsys%create_pw(spin_embed_pot)
642  END IF
643  ! Read the cubes
644  CALL read_embed_pot_cube(embed_pot, spin_embed_pot, qs_section, open_shell_embed)
645 
646  IF (.NOT. open_shell_embed) THEN
647  CALL set_qs_env(qs_env=qs_env, embed_pot=embed_pot)
648  ELSE
649  CALL set_qs_env(qs_env=qs_env, embed_pot=embed_pot, spin_embed_pot=spin_embed_pot)
650  END IF
651 
652  END SUBROUTINE given_embed_pot
653 
654 ! **************************************************************************************************
655 !> \brief ...
656 !> \param qs_env ...
657 !> \param embed_pot ...
658 !> \param spin_embed_pot ...
659 !> \param section ...
660 !> \param opt_embed ...
661 ! **************************************************************************************************
662  SUBROUTINE read_embed_pot(qs_env, embed_pot, spin_embed_pot, section, opt_embed)
663  TYPE(qs_environment_type), POINTER :: qs_env
664  TYPE(pw_r3d_rs_type), POINTER :: embed_pot, spin_embed_pot
665  TYPE(section_vals_type), POINTER :: section
666  TYPE(opt_embed_pot_type) :: opt_embed
667 
668  ! Read the potential as a vector in the auxiliary basis
669  IF (opt_embed%read_embed_pot) &
670  CALL read_embed_pot_vector(qs_env, embed_pot, spin_embed_pot, section, &
671  opt_embed%embed_pot_coef, opt_embed%open_shell_embed)
672  ! Read the potential as a cube (two cubes for open shell)
673  IF (opt_embed%read_embed_pot_cube) &
674  CALL read_embed_pot_cube(embed_pot, spin_embed_pot, section, opt_embed%open_shell_embed)
675 
676  END SUBROUTINE read_embed_pot
677 
678 ! **************************************************************************************************
679 !> \brief ...
680 !> \param embed_pot ...
681 !> \param spin_embed_pot ...
682 !> \param section ...
683 !> \param open_shell_embed ...
684 ! **************************************************************************************************
685  SUBROUTINE read_embed_pot_cube(embed_pot, spin_embed_pot, section, open_shell_embed)
686  TYPE(pw_r3d_rs_type), INTENT(IN) :: embed_pot, spin_embed_pot
687  TYPE(section_vals_type), POINTER :: section
688  LOGICAL :: open_shell_embed
689 
690  CHARACTER(LEN=default_path_length) :: filename
691  LOGICAL :: exist
692  REAL(kind=dp) :: scaling_factor
693 
694  exist = .false.
695  CALL section_vals_val_get(section, "EMBED_CUBE_FILE_NAME", c_val=filename)
696  INQUIRE (file=filename, exist=exist)
697  IF (.NOT. exist) &
698  cpabort("Embedding cube file not found. ")
699 
700  scaling_factor = 1.0_dp
701  CALL cp_cube_to_pw(embed_pot, filename, scaling_factor)
702 
703  ! Spin-dependent part of the potential
704  IF (open_shell_embed) THEN
705  exist = .false.
706  CALL section_vals_val_get(section, "EMBED_SPIN_CUBE_FILE_NAME", c_val=filename)
707  INQUIRE (file=filename, exist=exist)
708  IF (.NOT. exist) &
709  cpabort("Embedding spin cube file not found. ")
710 
711  scaling_factor = 1.0_dp
712  CALL cp_cube_to_pw(spin_embed_pot, filename, scaling_factor)
713  END IF
714 
715  END SUBROUTINE read_embed_pot_cube
716 
717 ! **************************************************************************************************
718 !> \brief Read the embedding potential from the binary file
719 !> \param qs_env ...
720 !> \param embed_pot ...
721 !> \param spin_embed_pot ...
722 !> \param section ...
723 !> \param embed_pot_coef ...
724 !> \param open_shell_embed ...
725 ! **************************************************************************************************
726  SUBROUTINE read_embed_pot_vector(qs_env, embed_pot, spin_embed_pot, section, embed_pot_coef, open_shell_embed)
727  TYPE(qs_environment_type), POINTER :: qs_env
728  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: embed_pot
729  TYPE(pw_r3d_rs_type), INTENT(IN), POINTER :: spin_embed_pot
730  TYPE(section_vals_type), POINTER :: section
731  TYPE(cp_fm_type), INTENT(IN) :: embed_pot_coef
732  LOGICAL, INTENT(IN) :: open_shell_embed
733 
734  CHARACTER(LEN=default_path_length) :: filename
735  INTEGER :: dimen_aux, dimen_restart_basis, &
736  dimen_var_aux, l_global, lll, &
737  nrow_local, restart_unit
738  INTEGER, DIMENSION(:), POINTER :: row_indices
739  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coef, coef_read
740  TYPE(cp_blacs_env_type), POINTER :: blacs_env
741  TYPE(cp_fm_struct_type), POINTER :: fm_struct
742  TYPE(cp_fm_type) :: my_embed_pot_coef
743  TYPE(mp_para_env_type), POINTER :: para_env
744 
745  ! Get the vector dimension
746  CALL find_aux_dimen(qs_env, dimen_aux)
747  IF (open_shell_embed) THEN
748  dimen_var_aux = dimen_aux*2
749  ELSE
750  dimen_var_aux = dimen_aux
751  END IF
752 
753  ! We need a temporary vector of coefficients
754  CALL get_qs_env(qs_env=qs_env, &
755  para_env=para_env)
756  NULLIFY (blacs_env)
757  NULLIFY (fm_struct)
758  CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env)
759  CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, &
760  nrow_global=dimen_var_aux, ncol_global=1)
761  CALL cp_fm_create(my_embed_pot_coef, fm_struct, name="my_pot_coef")
762 
763  CALL cp_fm_struct_release(fm_struct)
764  CALL cp_fm_set_all(my_embed_pot_coef, 0.0_dp)
765 
766  ! Read the coefficients vector
767  restart_unit = -1
768 
769  ! Allocate the attay to read the coefficients
770  ALLOCATE (coef(dimen_var_aux))
771  coef = 0.0_dp
772 
773  IF (para_env%is_source()) THEN
774 
775  ! Get the restart file name
776  CALL embed_restart_file_name(filename, section)
777 
778  CALL open_file(file_name=filename, &
779  file_action="READ", &
780  file_form="UNFORMATTED", &
781  file_status="OLD", &
782  unit_number=restart_unit)
783 
784  READ (restart_unit) dimen_restart_basis
785  ! Check the dimensions of the bases: the actual and the restart one
786  IF (.NOT. (dimen_restart_basis == dimen_aux)) &
787  cpabort("Wrong dimension of the embedding basis in the restart file.")
788 
789  ALLOCATE (coef_read(dimen_var_aux))
790  coef_read = 0.0_dp
791 
792  READ (restart_unit) coef_read
793  coef(:) = coef_read(:)
794  DEALLOCATE (coef_read)
795 
796  ! Close restart file
797  CALL close_file(unit_number=restart_unit)
798 
799  END IF
800 
801  ! Broadcast the coefficients on all processes
802  CALL para_env%bcast(coef)
803 
804  ! Copy to fm_type structure
805  ! Information about full matrix gradient
806  CALL cp_fm_get_info(matrix=my_embed_pot_coef, &
807  nrow_local=nrow_local, &
808  row_indices=row_indices)
809 
810  DO lll = 1, nrow_local
811  l_global = row_indices(lll)
812  my_embed_pot_coef%local_data(lll, 1) = coef(l_global)
813  END DO
814 
815  DEALLOCATE (coef)
816 
817  ! Copy to the my_embed_pot_coef to embed_pot_coef
818  CALL cp_fm_copy_general(my_embed_pot_coef, embed_pot_coef, para_env)
819 
820  ! Build the embedding potential on the grid
821  CALL update_embed_pot(embed_pot_coef, dimen_aux, embed_pot, spin_embed_pot, &
822  qs_env, .false., open_shell_embed)
823 
824  ! Release my_embed_pot_coef
825  CALL cp_fm_release(my_embed_pot_coef)
826 
827  ! Release blacs environment
828  CALL cp_blacs_env_release(blacs_env)
829 
830  END SUBROUTINE read_embed_pot_vector
831 
832 ! **************************************************************************************************
833 !> \brief Find the embedding restart file name
834 !> \param filename ...
835 !> \param section ...
836 ! **************************************************************************************************
837  SUBROUTINE embed_restart_file_name(filename, section)
838  CHARACTER(LEN=default_path_length), INTENT(OUT) :: filename
839  TYPE(section_vals_type), POINTER :: section
840 
841  LOGICAL :: exist
842 
843  exist = .false.
844  CALL section_vals_val_get(section, "EMBED_RESTART_FILE_NAME", c_val=filename)
845  INQUIRE (file=filename, exist=exist)
846  IF (.NOT. exist) &
847  cpabort("Embedding restart file not found. ")
848 
849  END SUBROUTINE embed_restart_file_name
850 
851 ! **************************************************************************************************
852 !> \brief Deallocate stuff for optimizing embedding potential
853 !> \param opt_embed ...
854 ! **************************************************************************************************
855  SUBROUTINE release_opt_embed(opt_embed)
856  TYPE(opt_embed_pot_type) :: opt_embed
857 
858  INTEGER :: i_dens, i_spin, ikind
859 
860  IF (.NOT. opt_embed%grid_opt) THEN
861  CALL cp_fm_release(opt_embed%embed_pot_grad)
862  CALL cp_fm_release(opt_embed%embed_pot_coef)
863  CALL cp_fm_release(opt_embed%step)
864  CALL cp_fm_release(opt_embed%prev_step)
865  CALL cp_fm_release(opt_embed%embed_pot_hess)
866  CALL cp_fm_release(opt_embed%prev_embed_pot_grad)
867  CALL cp_fm_release(opt_embed%prev_embed_pot_coef)
868  CALL cp_fm_release(opt_embed%prev_embed_pot_hess)
869  CALL cp_fm_release(opt_embed%kinetic_mat)
870  DEALLOCATE (opt_embed%embed_pot_grad, opt_embed%embed_pot_coef, &
871  opt_embed%step, opt_embed%prev_step, opt_embed%embed_pot_hess, &
872  opt_embed%prev_embed_pot_grad, opt_embed%prev_embed_pot_coef, &
873  opt_embed%prev_embed_pot_hess, opt_embed%kinetic_mat)
874  DEALLOCATE (opt_embed%w_func)
875  DEALLOCATE (opt_embed%max_diff)
876  DEALLOCATE (opt_embed%int_diff)
877 
878  DO ikind = 1, SIZE(opt_embed%lri)
879  DEALLOCATE (opt_embed%lri(ikind)%v_int)
880  DEALLOCATE (opt_embed%lri(ikind)%acoef)
881  END DO
882  DEALLOCATE (opt_embed%lri)
883  END IF
884 
885  IF (ASSOCIATED(opt_embed%prev_subsys_dens)) THEN
886  DO i_dens = 1, SIZE(opt_embed%prev_subsys_dens)
887  CALL opt_embed%prev_subsys_dens(i_dens)%release()
888  END DO
889  DEALLOCATE (opt_embed%prev_subsys_dens)
890  END IF
891  DEALLOCATE (opt_embed%max_subsys_dens_diff)
892 
893  DEALLOCATE (opt_embed%all_nspins)
894 
895  IF (ASSOCIATED(opt_embed%const_pot)) THEN
896  CALL opt_embed%const_pot%release()
897  DEALLOCATE (opt_embed%const_pot)
898  END IF
899 
900  IF (ASSOCIATED(opt_embed%pot_diff)) THEN
901  CALL opt_embed%pot_diff%release()
902  DEALLOCATE (opt_embed%pot_diff)
903  END IF
904 
905  IF (ASSOCIATED(opt_embed%prev_embed_pot)) THEN
906  CALL opt_embed%prev_embed_pot%release()
907  DEALLOCATE (opt_embed%prev_embed_pot)
908  END IF
909  IF (ASSOCIATED(opt_embed%prev_spin_embed_pot)) THEN
910  CALL opt_embed%prev_spin_embed_pot%release()
911  DEALLOCATE (opt_embed%prev_spin_embed_pot)
912  END IF
913  IF (ASSOCIATED(opt_embed%v_w)) THEN
914  DO i_spin = 1, SIZE(opt_embed%v_w)
915  CALL opt_embed%v_w(i_spin)%release()
916  END DO
917  DEALLOCATE (opt_embed%v_w)
918  END IF
919 
920  END SUBROUTINE release_opt_embed
921 
922 ! **************************************************************************************************
923 !> \brief Calculates subsystem Coulomb potential from the RESP charges of the total system
924 !> \param v_rspace ...
925 !> \param rhs ...
926 !> \param mapping_section ...
927 !> \param qs_env ...
928 !> \param nforce_eval ...
929 !> \param iforce_eval ...
930 !> \param eta ...
931 ! **************************************************************************************************
932  SUBROUTINE coulomb_guess(v_rspace, rhs, mapping_section, qs_env, nforce_eval, iforce_eval, eta)
933  TYPE(pw_r3d_rs_type) :: v_rspace
934  REAL(kind=dp), DIMENSION(:), POINTER :: rhs
935  TYPE(section_vals_type), POINTER :: mapping_section
936  TYPE(qs_environment_type), POINTER :: qs_env
937  INTEGER :: nforce_eval, iforce_eval
938  REAL(kind=dp) :: eta
939 
940  INTEGER :: iparticle, jparticle, natom
941  INTEGER, DIMENSION(:), POINTER :: map_index
942  REAL(kind=dp) :: dvol, normalize_factor
943  REAL(kind=dp), DIMENSION(:), POINTER :: rhs_subsys
944  TYPE(particle_list_type), POINTER :: particles
945  TYPE(pw_c1d_gs_type) :: v_resp_gspace
946  TYPE(pw_env_type), POINTER :: pw_env
947  TYPE(pw_poisson_type), POINTER :: poisson_env
948  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
949  TYPE(pw_r3d_rs_type) :: rho_resp, v_resp_rspace
950  TYPE(qs_subsys_type), POINTER :: subsys
951 
952  ! Get available particles
953  NULLIFY (subsys)
954  CALL get_qs_env(qs_env=qs_env, subsys=subsys, pw_env=pw_env)
955  CALL qs_subsys_get(subsys, particles=particles)
956  natom = particles%n_els
957 
958  ALLOCATE (rhs_subsys(natom))
959 
960  NULLIFY (map_index)
961  CALL get_subsys_map_index(mapping_section, natom, iforce_eval, nforce_eval, &
962  map_index, .true.)
963 
964  ! Mapping particles from iforce_eval environment to the embed env
965  DO iparticle = 1, natom
966  jparticle = map_index(iparticle)
967  rhs_subsys(iparticle) = rhs(jparticle)
968  END DO
969 
970  ! Prepare plane waves
971  NULLIFY (auxbas_pw_pool)
972 
973  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
974  poisson_env=poisson_env)
975 
976  CALL auxbas_pw_pool%create_pw(v_resp_gspace)
977 
978  CALL auxbas_pw_pool%create_pw(v_resp_rspace)
979 
980  CALL auxbas_pw_pool%create_pw(rho_resp)
981 
982  ! Calculate charge density
983  CALL pw_zero(rho_resp)
984  CALL calculate_rho_resp_all(rho_resp, rhs_subsys, natom, eta, qs_env)
985 
986  ! Calculate potential
987  CALL pw_poisson_solve(poisson_env, rho_resp, &
988  vhartree=v_resp_rspace)
989  dvol = v_resp_rspace%pw_grid%dvol
990  CALL pw_scale(v_resp_rspace, dvol)
991  normalize_factor = sqrt((eta/pi)**3)
992  !normalize_factor = -2.0_dp
993  CALL pw_scale(v_resp_rspace, normalize_factor)
994 
995  ! Hard copy potential
996  CALL pw_copy(v_resp_rspace, v_rspace)
997 
998  ! Release plane waves
999  CALL v_resp_gspace%release()
1000  CALL v_resp_rspace%release()
1001  CALL rho_resp%release()
1002 
1003  ! Deallocate map_index array
1004  DEALLOCATE (map_index)
1005  ! Deallocate charges
1006  DEALLOCATE (rhs_subsys)
1007 
1008  END SUBROUTINE coulomb_guess
1009 
1010 ! **************************************************************************************************
1011 !> \brief Creates a subsystem embedding potential
1012 !> \param qs_env ...
1013 !> \param embed_pot ...
1014 !> \param embed_pot_subsys ...
1015 !> \param spin_embed_pot ...
1016 !> \param spin_embed_pot_subsys ...
1017 !> \param open_shell_embed ...
1018 !> \param change_spin_sign ...
1019 !> \author Vladimir Rybkin
1020 ! **************************************************************************************************
1021  SUBROUTINE make_subsys_embed_pot(qs_env, embed_pot, embed_pot_subsys, &
1022  spin_embed_pot, spin_embed_pot_subsys, open_shell_embed, &
1023  change_spin_sign)
1024  TYPE(qs_environment_type), POINTER :: qs_env
1025  TYPE(pw_r3d_rs_type), INTENT(IN) :: embed_pot
1026  TYPE(pw_r3d_rs_type), POINTER :: embed_pot_subsys
1027  TYPE(pw_r3d_rs_type), INTENT(IN), POINTER :: spin_embed_pot
1028  TYPE(pw_r3d_rs_type), POINTER :: spin_embed_pot_subsys
1029  LOGICAL :: open_shell_embed, change_spin_sign
1030 
1031  TYPE(pw_env_type), POINTER :: pw_env
1032  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool_subsys
1033 
1034  ! Extract plane waves environment
1035  CALL get_qs_env(qs_env, pw_env=pw_env)
1036 
1037  ! Prepare plane-waves pool
1038  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool_subsys)
1039 
1040  ! Create embedding potential and set to zero
1041  NULLIFY (embed_pot_subsys)
1042  ALLOCATE (embed_pot_subsys)
1043  CALL auxbas_pw_pool_subsys%create_pw(embed_pot_subsys)
1044 
1045  ! Hard copy the grid
1046  CALL pw_copy(embed_pot, embed_pot_subsys)
1047 
1048  IF (open_shell_embed) THEN
1049  NULLIFY (spin_embed_pot_subsys)
1050  ALLOCATE (spin_embed_pot_subsys)
1051  CALL auxbas_pw_pool_subsys%create_pw(spin_embed_pot_subsys)
1052  ! Hard copy the grid
1053  IF (change_spin_sign) THEN
1054  CALL pw_axpy(spin_embed_pot, spin_embed_pot_subsys, -1.0_dp, 0.0_dp, allow_noncompatible_grids=.true.)
1055  ELSE
1056  CALL pw_copy(spin_embed_pot, spin_embed_pot_subsys)
1057  END IF
1058  END IF
1059 
1060  END SUBROUTINE make_subsys_embed_pot
1061 
1062 ! **************************************************************************************************
1063 !> \brief Calculates the derivative of the embedding potential wrt to the expansion coefficients
1064 !> \param qs_env ...
1065 !> \param diff_rho_r ...
1066 !> \param diff_rho_spin ...
1067 !> \param opt_embed ...
1068 !> \author Vladimir Rybkin
1069 ! **************************************************************************************************
1070 
1071  SUBROUTINE calculate_embed_pot_grad(qs_env, diff_rho_r, diff_rho_spin, opt_embed)
1072  TYPE(qs_environment_type), POINTER :: qs_env
1073  TYPE(pw_r3d_rs_type), INTENT(IN) :: diff_rho_r, diff_rho_spin
1074  TYPE(opt_embed_pot_type) :: opt_embed
1075 
1076  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_embed_pot_grad'
1077 
1078  INTEGER :: handle
1079  TYPE(cp_blacs_env_type), POINTER :: blacs_env
1080  TYPE(cp_fm_struct_type), POINTER :: fm_struct
1081  TYPE(cp_fm_type) :: embed_pot_coeff_spin, &
1082  embed_pot_coeff_spinless, &
1083  regular_term, spin_reg, spinless_reg
1084  TYPE(mp_para_env_type), POINTER :: para_env
1085  TYPE(pw_env_type), POINTER :: pw_env
1086  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1087 
1088  CALL timeset(routinen, handle)
1089 
1090  ! We destroy the previous gradient and Hessian:
1091  ! current data are now previous data
1092  CALL cp_fm_to_fm(opt_embed%embed_pot_grad, opt_embed%prev_embed_pot_grad)
1093  CALL cp_fm_to_fm(opt_embed%embed_pot_Hess, opt_embed%prev_embed_pot_Hess)
1094 
1095  NULLIFY (pw_env)
1096 
1097  CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, para_env=para_env)
1098 
1099  ! Get plane waves pool
1100  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1101 
1102  ! Calculate potential gradient coefficients
1103  CALL calculate_embed_pot_grad_inner(qs_env, opt_embed%dimen_aux, diff_rho_r, diff_rho_spin, &
1104  opt_embed%embed_pot_grad, &
1105  opt_embed%open_shell_embed, opt_embed%lri)
1106 
1107  ! Add regularization with kinetic matrix
1108  IF (opt_embed%i_iter .EQ. 1) THEN ! Else it is kept in memory
1109  CALL compute_kinetic_mat(qs_env, opt_embed%kinetic_mat)
1110  END IF
1111 
1112  CALL cp_fm_get_info(matrix=opt_embed%embed_pot_grad, &
1113  matrix_struct=fm_struct)
1114  CALL cp_fm_create(regular_term, fm_struct, name="regular_term")
1115  CALL cp_fm_set_all(regular_term, 0.0_dp)
1116 
1117  ! In case of open shell embedding we need two terms of dimen_aux=dimen_var_aux/2 for
1118  ! the spinless and the spin parts
1119  IF (opt_embed%open_shell_embed) THEN
1120  ! Prepare auxiliary full matrices
1121  NULLIFY (fm_struct, blacs_env)
1122 
1123  !CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env)
1124 
1125  CALL cp_fm_get_info(matrix=opt_embed%embed_pot_coef, context=blacs_env)
1126  CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, &
1127  nrow_global=opt_embed%dimen_aux, ncol_global=1)
1128  CALL cp_fm_create(embed_pot_coeff_spinless, fm_struct, name="pot_coeff_spinless")
1129  CALL cp_fm_create(embed_pot_coeff_spin, fm_struct, name="pot_coeff_spin")
1130  CALL cp_fm_create(spinless_reg, fm_struct, name="spinless_reg")
1131  CALL cp_fm_create(spin_reg, fm_struct, name="spin_reg")
1132  CALL cp_fm_set_all(embed_pot_coeff_spinless, 0.0_dp)
1133  CALL cp_fm_set_all(embed_pot_coeff_spin, 0.0_dp)
1134  CALL cp_fm_set_all(spinless_reg, 0.0_dp)
1135  CALL cp_fm_set_all(spin_reg, 0.0_dp)
1136  CALL cp_fm_struct_release(fm_struct)
1137 
1138  ! Copy coefficients to the auxiliary structures
1139  CALL cp_fm_to_fm_submat(msource=opt_embed%embed_pot_coef, &
1140  mtarget=embed_pot_coeff_spinless, &
1141  nrow=opt_embed%dimen_aux, ncol=1, &
1142  s_firstrow=1, s_firstcol=1, &
1143  t_firstrow=1, t_firstcol=1)
1144  CALL cp_fm_to_fm_submat(msource=opt_embed%embed_pot_coef, &
1145  mtarget=embed_pot_coeff_spin, &
1146  nrow=opt_embed%dimen_aux, ncol=1, &
1147  s_firstrow=opt_embed%dimen_aux + 1, s_firstcol=1, &
1148  t_firstrow=1, t_firstcol=1)
1149  ! Multiply
1150  CALL parallel_gemm(transa="N", transb="N", m=opt_embed%dimen_aux, n=1, &
1151  k=opt_embed%dimen_aux, alpha=1.0_dp, &
1152  matrix_a=opt_embed%kinetic_mat, matrix_b=embed_pot_coeff_spinless, &
1153  beta=0.0_dp, matrix_c=spinless_reg)
1154  CALL parallel_gemm(transa="N", transb="N", m=opt_embed%dimen_aux, n=1, &
1155  k=opt_embed%dimen_aux, alpha=1.0_dp, &
1156  matrix_a=opt_embed%kinetic_mat, matrix_b=embed_pot_coeff_spin, &
1157  beta=0.0_dp, matrix_c=spin_reg)
1158  ! Copy from the auxiliary structures to the full regularization term
1159  CALL cp_fm_to_fm_submat(msource=spinless_reg, &
1160  mtarget=regular_term, &
1161  nrow=opt_embed%dimen_aux, ncol=1, &
1162  s_firstrow=1, s_firstcol=1, &
1163  t_firstrow=1, t_firstcol=1)
1164  CALL cp_fm_to_fm_submat(msource=spin_reg, &
1165  mtarget=regular_term, &
1166  nrow=opt_embed%dimen_aux, ncol=1, &
1167  s_firstrow=1, s_firstcol=1, &
1168  t_firstrow=opt_embed%dimen_aux + 1, t_firstcol=1)
1169  ! Release internally used auxiliary structures
1170  CALL cp_fm_release(embed_pot_coeff_spinless)
1171  CALL cp_fm_release(embed_pot_coeff_spin)
1172  CALL cp_fm_release(spin_reg)
1173  CALL cp_fm_release(spinless_reg)
1174 
1175  ELSE ! Simply multiply
1176  CALL parallel_gemm(transa="N", transb="N", m=opt_embed%dimen_var_aux, n=1, &
1177  k=opt_embed%dimen_var_aux, alpha=1.0_dp, &
1178  matrix_a=opt_embed%kinetic_mat, matrix_b=opt_embed%embed_pot_coef, &
1179  beta=0.0_dp, matrix_c=regular_term)
1180  END IF
1181 
1182  ! Scale by the regularization parameter and add to the gradient
1183  CALL cp_fm_scale_and_add(1.0_dp, opt_embed%embed_pot_grad, 4.0_dp*opt_embed%lambda, regular_term)
1184 
1185  ! Calculate the regularization contribution to the energy functional
1186  CALL cp_fm_trace(opt_embed%embed_pot_coef, regular_term, opt_embed%reg_term)
1187  opt_embed%reg_term = 2.0_dp*opt_embed%lambda*opt_embed%reg_term
1188 
1189  ! Deallocate regular term
1190  CALL cp_fm_release(regular_term)
1191 
1192  CALL timestop(handle)
1193 
1194  END SUBROUTINE calculate_embed_pot_grad
1195 
1196 ! **************************************************************************************************
1197 !> \brief Performs integration for the embedding potential gradient
1198 !> \param qs_env ...
1199 !> \param dimen_aux ...
1200 !> \param rho_r ...
1201 !> \param rho_spin ...
1202 !> \param embed_pot_grad ...
1203 !> \param open_shell_embed ...
1204 !> \param lri ...
1205 !> \author Vladimir Rybkin
1206 ! **************************************************************************************************
1207  SUBROUTINE calculate_embed_pot_grad_inner(qs_env, dimen_aux, rho_r, rho_spin, embed_pot_grad, &
1208  open_shell_embed, lri)
1209  TYPE(qs_environment_type), POINTER :: qs_env
1210  INTEGER :: dimen_aux
1211  TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_r, rho_spin
1212  TYPE(cp_fm_type), INTENT(IN) :: embed_pot_grad
1213  LOGICAL :: open_shell_embed
1214  TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri
1215 
1216  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_embed_pot_grad_inner'
1217 
1218  INTEGER :: handle, iatom, ikind, l_global, lll, &
1219  nrow_local, nsgf, start_pos
1220  INTEGER, DIMENSION(:), POINTER :: row_indices
1221  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: pot_grad
1222  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1223  TYPE(cell_type), POINTER :: cell
1224  TYPE(dft_control_type), POINTER :: dft_control
1225  TYPE(mp_para_env_type), POINTER :: para_env
1226  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1227  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1228 
1229 ! Needed to store integrals
1230 
1231  CALL timeset(routinen, handle)
1232 
1233  CALL get_qs_env(qs_env=qs_env, &
1234  particle_set=particle_set, &
1235  qs_kind_set=qs_kind_set, &
1236  dft_control=dft_control, &
1237  cell=cell, &
1238  atomic_kind_set=atomic_kind_set, &
1239  para_env=para_env)
1240 
1241  ! Create wf_vector and gradient
1242  IF (open_shell_embed) THEN
1243  ALLOCATE (pot_grad(dimen_aux*2))
1244  ELSE
1245  ALLOCATE (pot_grad(dimen_aux))
1246  END IF
1247 
1248  ! Use lri subroutine
1249  DO ikind = 1, SIZE(lri)
1250  lri(ikind)%v_int = 0.0_dp
1251  END DO
1252 
1253  CALL integrate_v_rspace_one_center(rho_r, qs_env, lri, &
1254  .false., "RI_AUX")
1255  DO ikind = 1, SIZE(lri)
1256  CALL para_env%sum(lri(ikind)%v_int)
1257  END DO
1258 
1259  pot_grad = 0.0_dp
1260  start_pos = 1
1261  DO ikind = 1, SIZE(lri)
1262  DO iatom = 1, SIZE(lri(ikind)%v_int, dim=1)
1263  nsgf = SIZE(lri(ikind)%v_int(iatom, :))
1264  pot_grad(start_pos:start_pos + nsgf - 1) = lri(ikind)%v_int(iatom, :)
1265  start_pos = start_pos + nsgf
1266  END DO
1267  END DO
1268 
1269  ! Open-shell embedding
1270  IF (open_shell_embed) THEN
1271  DO ikind = 1, SIZE(lri)
1272  lri(ikind)%v_int = 0.0_dp
1273  END DO
1274 
1275  CALL integrate_v_rspace_one_center(rho_spin, qs_env, lri, &
1276  .false., "RI_AUX")
1277  DO ikind = 1, SIZE(lri)
1278  CALL para_env%sum(lri(ikind)%v_int)
1279  END DO
1280 
1281  start_pos = dimen_aux + 1
1282  DO ikind = 1, SIZE(lri)
1283  DO iatom = 1, SIZE(lri(ikind)%v_int, dim=1)
1284  nsgf = SIZE(lri(ikind)%v_int(iatom, :))
1285  pot_grad(start_pos:start_pos + nsgf - 1) = lri(ikind)%v_int(iatom, :)
1286  start_pos = start_pos + nsgf
1287  END DO
1288  END DO
1289  END IF
1290 
1291  ! Scale by the cell volume
1292  pot_grad = pot_grad*rho_r%pw_grid%dvol
1293 
1294  ! Information about full matrix gradient
1295  CALL cp_fm_get_info(matrix=embed_pot_grad, &
1296  nrow_local=nrow_local, &
1297  row_indices=row_indices)
1298 
1299  ! Copy the gradient into the full matrix
1300  DO lll = 1, nrow_local
1301  l_global = row_indices(lll)
1302  embed_pot_grad%local_data(lll, 1) = pot_grad(l_global)
1303  END DO
1304 
1305  DEALLOCATE (pot_grad)
1306 
1307  CALL timestop(handle)
1308 
1309  END SUBROUTINE calculate_embed_pot_grad_inner
1310 
1311 ! **************************************************************************************************
1312 !> \brief Calculates kinetic energy matrix in auxiliary basis in the fm format
1313 !> \param qs_env ...
1314 !> \param kinetic_mat ...
1315 !> \author Vladimir Rybkin
1316 ! **************************************************************************************************
1317  SUBROUTINE compute_kinetic_mat(qs_env, kinetic_mat)
1318  TYPE(qs_environment_type), POINTER :: qs_env
1319  TYPE(cp_fm_type), INTENT(IN) :: kinetic_mat
1320 
1321  CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_kinetic_mat'
1322 
1323  INTEGER :: handle
1324  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_t
1325  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1326  POINTER :: sab_orb
1327  TYPE(qs_ks_env_type), POINTER :: ks_env
1328 
1329  CALL timeset(routinen, handle)
1330 
1331  NULLIFY (ks_env, sab_orb, matrix_t)
1332 
1333  ! First, get the dbcsr structure from the overlap matrix
1334  CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, sab_orb=sab_orb)
1335 
1336  ! Calculate kinetic matrix
1337  CALL build_kinetic_matrix(ks_env, matrix_t=matrix_t, &
1338  matrix_name="KINETIC ENERGY MATRIX", &
1339  basis_type="RI_AUX", &
1340  sab_nl=sab_orb, calculate_forces=.false.)
1341 
1342  ! Change to the fm format
1343  CALL copy_dbcsr_to_fm(matrix_t(1)%matrix, kinetic_mat)
1344 
1345  ! Release memory
1346  CALL dbcsr_deallocate_matrix_set(matrix_t)
1347 
1348  CALL timestop(handle)
1349 
1350  END SUBROUTINE compute_kinetic_mat
1351 
1352 ! **************************************************************************************************
1353 !> \brief Regularizes the Wu-Yang potential on the grid
1354 !> \param potential ...
1355 !> \param pw_env ...
1356 !> \param lambda ...
1357 !> \param reg_term ...
1358 ! **************************************************************************************************
1359  SUBROUTINE grid_regularize(potential, pw_env, lambda, reg_term)
1360 
1361  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: potential
1362  TYPE(pw_env_type), POINTER :: pw_env
1363  REAL(kind=dp) :: lambda, reg_term
1364 
1365  INTEGER :: i, j, k
1366  INTEGER, DIMENSION(3) :: lb, n, ub
1367  TYPE(pw_c1d_gs_type) :: dr2_pot, grid_reg_g, potential_g
1368  TYPE(pw_c1d_gs_type), DIMENSION(3) :: dpot_g
1369  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1370  TYPE(pw_r3d_rs_type) :: grid_reg, square_norm_dpot
1371  TYPE(pw_r3d_rs_type), DIMENSION(3) :: dpot
1372 
1373  !
1374  ! First, the contribution to the gradient
1375  !
1376 
1377  ! Get some of the grids ready
1378  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1379 
1380  CALL auxbas_pw_pool%create_pw(potential_g)
1381 
1382  CALL auxbas_pw_pool%create_pw(dr2_pot)
1383 
1384  CALL auxbas_pw_pool%create_pw(grid_reg)
1385 
1386  CALL auxbas_pw_pool%create_pw(grid_reg_g)
1387  CALL pw_zero(grid_reg_g)
1388 
1389  ! Transfer potential to the reciprocal space
1390  CALL pw_transfer(potential, potential_g)
1391 
1392  ! Calculate second derivatives: dx^2, dy^2, dz^2
1393  DO i = 1, 3
1394  CALL pw_dr2(potential_g, dr2_pot, i, i)
1395  CALL pw_axpy(dr2_pot, grid_reg_g, 1.0_dp)
1396  END DO
1397  ! Transfer potential to the real space
1398  CALL pw_transfer(grid_reg_g, grid_reg)
1399 
1400  ! Update the potential with a regularization term
1401  CALL pw_axpy(grid_reg, potential, -4.0_dp*lambda)
1402 
1403  !
1404  ! Second, the contribution to the functional
1405  !
1406  DO i = 1, 3
1407  CALL auxbas_pw_pool%create_pw(dpot(i))
1408  CALL auxbas_pw_pool%create_pw(dpot_g(i))
1409  END DO
1410 
1411  CALL auxbas_pw_pool%create_pw(square_norm_dpot)
1412 
1413  DO i = 1, 3
1414  n(:) = 0
1415  n(i) = 1
1416  CALL pw_copy(potential_g, dpot_g(i))
1417  CALL pw_derive(dpot_g(i), n(:))
1418  CALL pw_transfer(dpot_g(i), dpot(i))
1419  END DO
1420 
1421  lb(1:3) = square_norm_dpot%pw_grid%bounds_local(1, 1:3)
1422  ub(1:3) = square_norm_dpot%pw_grid%bounds_local(2, 1:3)
1423 !$OMP PARALLEL DO DEFAULT(NONE) &
1424 !$OMP PRIVATE(i,j,k) &
1425 !$OMP SHARED(dpot, lb, square_norm_dpot, ub)
1426  DO k = lb(3), ub(3)
1427  DO j = lb(2), ub(2)
1428  DO i = lb(1), ub(1)
1429  square_norm_dpot%array(i, j, k) = (dpot(1)%array(i, j, k)* &
1430  dpot(1)%array(i, j, k) + &
1431  dpot(2)%array(i, j, k)* &
1432  dpot(2)%array(i, j, k) + &
1433  dpot(3)%array(i, j, k)* &
1434  dpot(3)%array(i, j, k))
1435  END DO
1436  END DO
1437  END DO
1438 !$OMP END PARALLEL DO
1439 
1440  reg_term = 2*lambda*pw_integrate_function(fun=square_norm_dpot)
1441 
1442  ! Release
1443  CALL auxbas_pw_pool%give_back_pw(potential_g)
1444  CALL auxbas_pw_pool%give_back_pw(dr2_pot)
1445  CALL auxbas_pw_pool%give_back_pw(grid_reg)
1446  CALL auxbas_pw_pool%give_back_pw(grid_reg_g)
1447  CALL auxbas_pw_pool%give_back_pw(square_norm_dpot)
1448  DO i = 1, 3
1449  CALL auxbas_pw_pool%give_back_pw(dpot(i))
1450  CALL auxbas_pw_pool%give_back_pw(dpot_g(i))
1451  END DO
1452 
1453  END SUBROUTINE grid_regularize
1454 
1455 ! **************************************************************************************************
1456 !> \brief Takes maximization step in embedding potential optimization
1457 !> \param diff_rho_r ...
1458 !> \param diff_rho_spin ...
1459 !> \param opt_embed ...
1460 !> \param embed_pot ...
1461 !> \param spin_embed_pot ...
1462 !> \param rho_r_ref ...
1463 !> \param qs_env ...
1464 !> \author Vladimir Rybkin
1465 ! **************************************************************************************************
1466  SUBROUTINE opt_embed_step(diff_rho_r, diff_rho_spin, opt_embed, embed_pot, spin_embed_pot, rho_r_ref, qs_env)
1467 
1468  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: diff_rho_r, diff_rho_spin
1469  TYPE(opt_embed_pot_type) :: opt_embed
1470  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: embed_pot
1471  TYPE(pw_r3d_rs_type), INTENT(IN), POINTER :: spin_embed_pot
1472  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_ref
1473  TYPE(qs_environment_type), POINTER :: qs_env
1474 
1475  CHARACTER(LEN=*), PARAMETER :: routinen = 'opt_embed_step'
1476  REAL(kind=dp), PARAMETER :: thresh = 0.000001_dp
1477 
1478  INTEGER :: handle, l_global, lll, nrow_local
1479  INTEGER, DIMENSION(:), POINTER :: row_indices
1480  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenval
1481  TYPE(cp_fm_struct_type), POINTER :: fm_struct
1482  TYPE(cp_fm_type) :: diag_grad, diag_step, fm_u, fm_u_scale
1483  TYPE(pw_env_type), POINTER :: pw_env
1484 
1485  CALL timeset(routinen, handle)
1486 
1487  IF (opt_embed%grid_opt) THEN ! Grid based optimization
1488 
1489  opt_embed%step_len = opt_embed%trust_rad
1490  CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1491  IF (opt_embed%leeuwen) THEN
1492  CALL leeuwen_baerends_potential_update(pw_env, embed_pot, spin_embed_pot, diff_rho_r, diff_rho_spin, &
1493  rho_r_ref, opt_embed%open_shell_embed, opt_embed%trust_rad)
1494  ELSE
1495  IF (opt_embed%fab) THEN
1496  CALL fab_update(qs_env, rho_r_ref, opt_embed%prev_embed_pot, opt_embed%prev_spin_embed_pot, &
1497  embed_pot, spin_embed_pot, &
1498  diff_rho_r, diff_rho_spin, opt_embed%v_w, opt_embed%i_iter, opt_embed%trust_rad, &
1499  opt_embed%open_shell_embed, opt_embed%vw_cutoff, opt_embed%vw_smooth_cutoff_range)
1500  ELSE
1501  CALL grid_based_step(diff_rho_r, diff_rho_spin, pw_env, opt_embed, embed_pot, spin_embed_pot)
1502  END IF
1503  END IF
1504 
1505  ELSE ! Finite basis optimization
1506  ! If the previous step has been rejected, we go back to the previous expansion coefficients
1507  IF (.NOT. opt_embed%accept_step) &
1508  CALL cp_fm_scale_and_add(1.0_dp, opt_embed%embed_pot_coef, -1.0_dp, opt_embed%step)
1509 
1510  ! Do a simple steepest descent
1511  IF (opt_embed%steep_desc) THEN
1512  IF (opt_embed%i_iter .GT. 2) &
1513  opt_embed%trust_rad = barzilai_borwein(opt_embed%step, opt_embed%prev_step, &
1514  opt_embed%embed_pot_grad, opt_embed%prev_embed_pot_grad)
1515  IF (abs(opt_embed%trust_rad) .GT. opt_embed%max_trad) THEN
1516  IF (opt_embed%trust_rad .GT. 0.0_dp) THEN
1517  opt_embed%trust_rad = opt_embed%max_trad
1518  ELSE
1519  opt_embed%trust_rad = -opt_embed%max_trad
1520  END IF
1521  END IF
1522 
1523  CALL cp_fm_to_fm(opt_embed%step, opt_embed%prev_step)
1524  CALL cp_fm_scale_and_add(0.0_dp, opt_embed%prev_step, 1.0_dp, opt_embed%step)
1525  CALL cp_fm_set_all(opt_embed%step, 0.0_dp)
1526  CALL cp_fm_scale_and_add(1.0_dp, opt_embed%step, opt_embed%trust_rad, opt_embed%embed_pot_grad)
1527  opt_embed%step_len = opt_embed%trust_rad
1528  ELSE
1529 
1530  ! First, update the Hessian inverse if needed
1531  IF (opt_embed%i_iter > 1) THEN
1532  IF (opt_embed%accept_step) & ! We don't update Hessian if the step has been rejected
1533  CALL symm_rank_one_update(opt_embed%embed_pot_grad, opt_embed%prev_embed_pot_grad, &
1534  opt_embed%step, opt_embed%prev_embed_pot_Hess, opt_embed%embed_pot_Hess)
1535  END IF
1536 
1537  ! Add regularization term to the Hessian
1538  !CALL cp_fm_scale_and_add(1.0_dp, opt_embed%embed_pot_Hess, 4.0_dp*opt_embed%lambda, &
1539  ! opt_embed%kinetic_mat)
1540 
1541  ! Else use the first initial Hessian. Now it's just the unit matrix: embed_pot_hess
1542  ! Second, invert the Hessian
1543  ALLOCATE (eigenval(opt_embed%dimen_var_aux))
1544  eigenval = 0.0_dp
1545  CALL cp_fm_get_info(matrix=opt_embed%embed_pot_hess, &
1546  matrix_struct=fm_struct)
1547  CALL cp_fm_create(fm_u, fm_struct, name="fm_U")
1548  CALL cp_fm_create(fm_u_scale, fm_struct, name="fm_U")
1549  CALL cp_fm_set_all(fm_u, 0.0_dp)
1550  CALL cp_fm_set_all(fm_u_scale, 0.0_dp)
1551  CALL cp_fm_get_info(matrix=opt_embed%embed_pot_grad, &
1552  matrix_struct=fm_struct)
1553  CALL cp_fm_create(diag_grad, fm_struct, name="diag_grad")
1554  CALL cp_fm_set_all(diag_grad, 0.0_dp)
1555  CALL cp_fm_create(diag_step, fm_struct, name="diag_step")
1556  CALL cp_fm_set_all(diag_step, 0.0_dp)
1557 
1558  ! Store the Hessian as it will be destroyed in diagonalization: use fm_U_scal for it
1559  CALL cp_fm_to_fm(opt_embed%embed_pot_hess, fm_u_scale)
1560 
1561  ! Diagonalize Hessian
1562  CALL choose_eigv_solver(opt_embed%embed_pot_hess, fm_u, eigenval)
1563 
1564  ! Copy the Hessian back
1565  CALL cp_fm_to_fm(fm_u_scale, opt_embed%embed_pot_hess)
1566 
1567  ! Find the step in diagonal representation, begin with gradient
1568  CALL parallel_gemm(transa="T", transb="N", m=opt_embed%dimen_var_aux, n=1, &
1569  k=opt_embed%dimen_var_aux, alpha=1.0_dp, &
1570  matrix_a=fm_u, matrix_b=opt_embed%embed_pot_grad, beta=0.0_dp, &
1571  matrix_c=diag_grad)
1572 
1573  CALL cp_fm_get_info(matrix=opt_embed%embed_pot_coef, &
1574  nrow_local=nrow_local, &
1575  row_indices=row_indices)
1576 
1577  DO lll = 1, nrow_local
1578  l_global = row_indices(lll)
1579  IF (abs(eigenval(l_global)) .GE. thresh) THEN
1580  diag_step%local_data(lll, 1) = &
1581  -diag_grad%local_data(lll, 1)/(eigenval(l_global))
1582  ELSE
1583  diag_step%local_data(lll, 1) = 0.0_dp
1584  END IF
1585  END DO
1586  CALL cp_fm_trace(diag_step, diag_step, opt_embed%step_len)
1587 
1588  ! Transform step to a non-diagonal representation
1589  CALL parallel_gemm(transa="N", transb="N", m=opt_embed%dimen_var_aux, n=1, &
1590  k=opt_embed%dimen_var_aux, alpha=1.0_dp, &
1591  matrix_a=fm_u, matrix_b=diag_step, beta=0.0_dp, &
1592  matrix_c=opt_embed%step)
1593 
1594  ! Now use fm_U_scale for scaled eigenvectors
1595  CALL cp_fm_to_fm(fm_u, fm_u_scale)
1596  CALL cp_fm_column_scale(fm_u_scale, eigenval)
1597 
1598  CALL cp_fm_release(fm_u_scale)
1599 
1600  ! Scale the step to fit within the trust radius: it it's less already,
1601  ! then take the Newton step
1602  CALL cp_fm_trace(opt_embed%step, opt_embed%step, opt_embed%step_len)
1603  IF (opt_embed%step_len .GT. opt_embed%trust_rad) THEN
1604 
1605  IF (opt_embed%level_shift) THEN
1606  ! Find a level shift parameter and apply it
1607  CALL level_shift(opt_embed, diag_grad, eigenval, diag_step)
1608  ELSE ! Just scale
1609  CALL cp_fm_trace(diag_step, diag_step, opt_embed%step_len)
1610  CALL cp_fm_scale(opt_embed%trust_rad/opt_embed%step_len, diag_step)
1611  END IF
1612  CALL cp_fm_trace(diag_step, diag_step, opt_embed%step_len)
1613  ! Transform step to a non-diagonal representation
1614  CALL parallel_gemm(transa="N", transb="N", m=opt_embed%dimen_var_aux, n=1, &
1615  k=opt_embed%dimen_var_aux, alpha=1.0_dp, &
1616  matrix_a=fm_u, matrix_b=diag_step, beta=0.0_dp, &
1617  matrix_c=opt_embed%step)
1618  CALL cp_fm_trace(opt_embed%step, opt_embed%step, opt_embed%step_len)
1619 
1620  ! Recalculate step in diagonal representation
1621  opt_embed%newton_step = .false.
1622  ELSE
1623  opt_embed%newton_step = .true.
1624  END IF
1625 
1626  ! Release some memory
1627  DEALLOCATE (eigenval)
1628  ! Release more memory
1629  CALL cp_fm_release(diag_grad)
1630  CALL cp_fm_release(diag_step)
1631  CALL cp_fm_release(fm_u)
1632 
1633  END IF ! grad_descent
1634 
1635  ! Update the coefficients
1636  CALL cp_fm_scale_and_add(1.0_dp, opt_embed%embed_pot_coef, 1.0_dp, opt_embed%step)
1637 
1638  ! Update the embedding potential
1639  CALL update_embed_pot(opt_embed%embed_pot_coef, opt_embed%dimen_aux, embed_pot, &
1640  spin_embed_pot, qs_env, opt_embed%add_const_pot, &
1641  opt_embed%open_shell_embed, opt_embed%const_pot)
1642  END IF ! Grid-based optimization
1643 
1644  CALL timestop(handle)
1645 
1646  END SUBROUTINE opt_embed_step
1647 
1648 !
1649 ! **************************************************************************************************
1650 !> \brief ...
1651 !> \param diff_rho_r ...
1652 !> \param diff_rho_spin ...
1653 !> \param pw_env ...
1654 !> \param opt_embed ...
1655 !> \param embed_pot ...
1656 !> \param spin_embed_pot ...
1657 ! **************************************************************************************************
1658  SUBROUTINE grid_based_step(diff_rho_r, diff_rho_spin, pw_env, opt_embed, embed_pot, spin_embed_pot)
1659 
1660  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: diff_rho_r, diff_rho_spin
1661  TYPE(pw_env_type), POINTER :: pw_env
1662  TYPE(opt_embed_pot_type) :: opt_embed
1663  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: embed_pot
1664  TYPE(pw_r3d_rs_type), POINTER :: spin_embed_pot
1665 
1666  CHARACTER(LEN=*), PARAMETER :: routinen = 'grid_based_step'
1667 
1668  INTEGER :: handle
1669  REAL(kind=dp) :: my_reg_term
1670 
1671  CALL timeset(routinen, handle)
1672 
1673  ! Take the step for spin-free part
1674  CALL pw_axpy(diff_rho_r, embed_pot, opt_embed%step_len)
1675  ! Regularize
1676  CALL grid_regularize(embed_pot, pw_env, opt_embed%lambda, my_reg_term)
1677  opt_embed%reg_term = opt_embed%reg_term + my_reg_term
1678 
1679  IF (opt_embed%open_shell_embed) THEN
1680  CALL pw_axpy(diff_rho_spin, spin_embed_pot, opt_embed%step_len)
1681  CALL grid_regularize(spin_embed_pot, pw_env, opt_embed%lambda, my_reg_term)
1682  opt_embed%reg_term = opt_embed%reg_term + my_reg_term
1683  END IF
1684 
1685  CALL timestop(handle)
1686 
1687  END SUBROUTINE grid_based_step
1688 
1689 ! **************************************************************************************************
1690 !> \brief ... Adds variable part of to the embedding potential
1691 !> \param embed_pot_coef ...
1692 !> \param dimen_aux ...
1693 !> \param embed_pot ...
1694 !> \param spin_embed_pot ...
1695 !> \param qs_env ...
1696 !> \param add_const_pot ...
1697 !> \param open_shell_embed ...
1698 !> \param const_pot ...
1699 !> \author Vladimir Rybkin
1700 ! **************************************************************************************************
1701 
1702  SUBROUTINE update_embed_pot(embed_pot_coef, dimen_aux, embed_pot, spin_embed_pot, &
1703  qs_env, add_const_pot, open_shell_embed, const_pot)
1704  TYPE(cp_fm_type), INTENT(IN) :: embed_pot_coef
1705  INTEGER :: dimen_aux
1706  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: embed_pot
1707  TYPE(pw_r3d_rs_type), INTENT(IN), POINTER :: spin_embed_pot
1708  TYPE(qs_environment_type), POINTER :: qs_env
1709  LOGICAL :: add_const_pot, open_shell_embed
1710  TYPE(pw_r3d_rs_type), INTENT(IN), OPTIONAL :: const_pot
1711 
1712  CHARACTER(LEN=*), PARAMETER :: routinen = 'update_embed_pot'
1713 
1714  INTEGER :: handle, l_global, lll, nrow_local
1715  INTEGER, DIMENSION(:), POINTER :: row_indices
1716  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: wf_vector
1717  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1718  TYPE(cell_type), POINTER :: cell
1719  TYPE(cp_blacs_env_type), POINTER :: blacs_env
1720  TYPE(cp_fm_struct_type), POINTER :: fm_struct
1721  TYPE(cp_fm_type) :: embed_pot_coef_spin, &
1722  embed_pot_coef_spinless
1723  TYPE(cp_fm_type), POINTER :: mo_coeff
1724  TYPE(dft_control_type), POINTER :: dft_control
1725  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1726  TYPE(mp_para_env_type), POINTER :: para_env
1727  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1728  TYPE(pw_c1d_gs_type) :: rho_g
1729  TYPE(pw_env_type), POINTER :: pw_env
1730  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1731  TYPE(pw_r3d_rs_type) :: psi_l
1732  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1733 
1734  CALL timeset(routinen, handle)
1735  ! Get MO coefficients: we need only the structure, therefore don't care about the spin
1736  CALL get_qs_env(qs_env=qs_env, &
1737  particle_set=particle_set, &
1738  qs_kind_set=qs_kind_set, &
1739  dft_control=dft_control, &
1740  cell=cell, &
1741  atomic_kind_set=atomic_kind_set, &
1742  pw_env=pw_env, mos=mos, para_env=para_env)
1743  CALL get_mo_set(mo_set=mos(1), mo_coeff=mo_coeff)
1744 
1745  ! Get plane waves pool
1746  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1747 
1748  ! get some of the grids ready
1749  CALL auxbas_pw_pool%create_pw(rho_g)
1750 
1751  CALL auxbas_pw_pool%create_pw(psi_l)
1752 
1753  ! Create wf_vector and auxiliary wave functions
1754  ALLOCATE (wf_vector(dimen_aux))
1755  wf_vector = 0.0_dp
1756 
1757  ! Create auxiliary full matrices for open-shell case
1758  IF (open_shell_embed) THEN
1759  NULLIFY (blacs_env)
1760  CALL cp_fm_get_info(matrix=embed_pot_coef, context=blacs_env)
1761  CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, &
1762  nrow_global=dimen_aux, ncol_global=1)
1763  CALL cp_fm_create(embed_pot_coef_spinless, fm_struct, name="pot_coeff_spinless")
1764  CALL cp_fm_create(embed_pot_coef_spin, fm_struct, name="pot_coeff_spin")
1765  CALL cp_fm_set_all(embed_pot_coef_spinless, 0.0_dp)
1766  CALL cp_fm_set_all(embed_pot_coef_spin, 0.0_dp)
1767  CALL cp_fm_struct_release(fm_struct)
1768 
1769  ! Copy coefficients to the auxiliary structures
1770  CALL cp_fm_to_fm_submat(embed_pot_coef, &
1771  mtarget=embed_pot_coef_spinless, &
1772  nrow=dimen_aux, ncol=1, &
1773  s_firstrow=1, s_firstcol=1, &
1774  t_firstrow=1, t_firstcol=1)
1775  CALL cp_fm_to_fm_submat(embed_pot_coef, &
1776  mtarget=embed_pot_coef_spin, &
1777  nrow=dimen_aux, ncol=1, &
1778  s_firstrow=dimen_aux + 1, s_firstcol=1, &
1779  t_firstrow=1, t_firstcol=1)
1780 
1781  ! Spinless potential
1782  CALL cp_fm_get_info(matrix=embed_pot_coef_spinless, &
1783  nrow_local=nrow_local, &
1784  row_indices=row_indices)
1785 
1786  ! Copy fm_coeff to an array
1787  DO lll = 1, nrow_local
1788  l_global = row_indices(lll)
1789  wf_vector(l_global) = embed_pot_coef_spinless%local_data(lll, 1)
1790  END DO
1791  CALL para_env%sum(wf_vector)
1792 
1793  ! Calculate the variable part of the embedding potential
1794  CALL calculate_wavefunction(mo_coeff, 1, psi_l, rho_g, atomic_kind_set, &
1795  qs_kind_set, cell, dft_control, particle_set, pw_env, &
1796  basis_type="RI_AUX", &
1797  external_vector=wf_vector)
1798  ! Update the full embedding potential
1799  IF (add_const_pot) THEN
1800  CALL pw_copy(const_pot, embed_pot)
1801  ELSE
1802  CALL pw_zero(embed_pot)
1803  END IF
1804 
1805  CALL pw_axpy(psi_l, embed_pot)
1806 
1807  ! Spin-dependent potential
1808  wf_vector = 0.0_dp
1809  CALL cp_fm_get_info(matrix=embed_pot_coef_spin, &
1810  nrow_local=nrow_local, &
1811  row_indices=row_indices)
1812 
1813  ! Copy fm_coeff to an array
1814  DO lll = 1, nrow_local
1815  l_global = row_indices(lll)
1816  wf_vector(l_global) = embed_pot_coef_spin%local_data(lll, 1)
1817  END DO
1818  CALL para_env%sum(wf_vector)
1819 
1820  ! Calculate the variable part of the embedding potential
1821  CALL calculate_wavefunction(mo_coeff, 1, psi_l, rho_g, atomic_kind_set, &
1822  qs_kind_set, cell, dft_control, particle_set, pw_env, &
1823  basis_type="RI_AUX", &
1824  external_vector=wf_vector)
1825  ! No constant potential for spin-dependent potential
1826  CALL pw_zero(spin_embed_pot)
1827  CALL pw_axpy(psi_l, spin_embed_pot)
1828 
1829  ELSE ! Closed shell
1830 
1831  CALL cp_fm_get_info(matrix=embed_pot_coef, &
1832  nrow_local=nrow_local, &
1833  row_indices=row_indices)
1834 
1835  ! Copy fm_coeff to an array
1836  DO lll = 1, nrow_local
1837  l_global = row_indices(lll)
1838  wf_vector(l_global) = embed_pot_coef%local_data(lll, 1)
1839  END DO
1840  CALL para_env%sum(wf_vector)
1841 
1842  ! Calculate the variable part of the embedding potential
1843  CALL calculate_wavefunction(mo_coeff, 1, psi_l, rho_g, atomic_kind_set, &
1844  qs_kind_set, cell, dft_control, particle_set, pw_env)
1845 
1846  CALL calculate_wavefunction(mo_coeff, 1, psi_l, rho_g, atomic_kind_set, &
1847  qs_kind_set, cell, dft_control, particle_set, pw_env, &
1848  basis_type="RI_AUX", &
1849  external_vector=wf_vector)
1850 
1851  ! Update the full embedding potential
1852  IF (add_const_pot) THEN
1853  CALL pw_copy(const_pot, embed_pot)
1854  ELSE
1855  CALL pw_zero(embed_pot)
1856  END IF
1857 
1858  CALL pw_axpy(psi_l, embed_pot)
1859  END IF ! Open/closed shell
1860 
1861  ! Deallocate memory and release objects
1862  DEALLOCATE (wf_vector)
1863  CALL auxbas_pw_pool%give_back_pw(psi_l)
1864  CALL auxbas_pw_pool%give_back_pw(rho_g)
1865 
1866  IF (open_shell_embed) THEN
1867  CALL cp_fm_release(embed_pot_coef_spin)
1868  CALL cp_fm_release(embed_pot_coef_spinless)
1869  END IF
1870 
1871  CALL timestop(handle)
1872 
1873  END SUBROUTINE update_embed_pot
1874 
1875 ! **************************************************************************************************
1876 !> \brief BFGS update of the inverse Hessian in the full matrix format
1877 !> \param grad ...
1878 !> \param prev_grad ...
1879 !> \param step ...
1880 !> \param prev_inv_Hess ...
1881 !> \param inv_Hess ...
1882 !> \author Vladimir Rybkin
1883 ! **************************************************************************************************
1884  SUBROUTINE inv_hessian_update(grad, prev_grad, step, prev_inv_Hess, inv_Hess)
1885  TYPE(cp_fm_type), INTENT(IN) :: grad, prev_grad, step, prev_inv_hess, &
1886  inv_hess
1887 
1888  INTEGER :: mat_size
1889  REAL(kind=dp) :: factor1, s_dot_y, y_dot_b_inv_y
1890  TYPE(cp_fm_struct_type), POINTER :: fm_struct_mat, fm_struct_vec
1891  TYPE(cp_fm_type) :: b_inv_y, b_inv_y_s, s_s, s_y, s_y_b_inv, &
1892  y
1893 
1894  ! Recover the dimension
1895  CALL cp_fm_get_info(matrix=inv_hess, &
1896  nrow_global=mat_size)
1897 
1898  CALL cp_fm_set_all(inv_hess, 0.0_dp)
1899  CALL cp_fm_to_fm(prev_inv_hess, inv_hess)
1900 
1901  ! Get full matrix structures
1902  NULLIFY (fm_struct_mat, fm_struct_vec)
1903 
1904  CALL cp_fm_get_info(matrix=prev_inv_hess, &
1905  matrix_struct=fm_struct_mat)
1906  CALL cp_fm_get_info(matrix=grad, &
1907  matrix_struct=fm_struct_vec)
1908 
1909  ! Allocate intermediates
1910  CALL cp_fm_create(b_inv_y, fm_struct_vec, name="B_inv_y")
1911  CALL cp_fm_create(y, fm_struct_vec, name="y")
1912 
1913  CALL cp_fm_create(s_s, fm_struct_mat, name="s_s")
1914  CALL cp_fm_create(s_y, fm_struct_mat, name="s_y")
1915  CALL cp_fm_create(b_inv_y_s, fm_struct_mat, name="B_inv_y_s")
1916  CALL cp_fm_create(s_y_b_inv, fm_struct_mat, name="s_y_B_inv")
1917 
1918  CALL cp_fm_set_all(b_inv_y, 0.0_dp)
1919  CALL cp_fm_set_all(s_s, 0.0_dp)
1920  CALL cp_fm_set_all(s_y, 0.0_dp)
1921  CALL cp_fm_set_all(b_inv_y_s, 0.0_dp)
1922  CALL cp_fm_set_all(s_y_b_inv, 0.0_dp)
1923 
1924  ! Calculate intermediates
1925  ! y the is gradient difference
1926  CALL cp_fm_get_info(matrix=grad)
1927  CALL cp_fm_to_fm(grad, y)
1928  CALL cp_fm_scale_and_add(1.0_dp, y, -1.0_dp, prev_grad)
1929 
1930  ! First term
1931  CALL parallel_gemm(transa="N", transb="N", m=mat_size, n=1, &
1932  k=mat_size, alpha=1.0_dp, &
1933  matrix_a=prev_inv_hess, matrix_b=y, beta=0.0_dp, &
1934  matrix_c=b_inv_y)
1935 
1936  CALL parallel_gemm(transa="N", transb="T", m=mat_size, n=mat_size, &
1937  k=1, alpha=1.0_dp, &
1938  matrix_a=step, matrix_b=step, beta=0.0_dp, &
1939  matrix_c=s_s)
1940 
1941  CALL parallel_gemm(transa="N", transb="T", m=mat_size, n=mat_size, &
1942  k=1, alpha=1.0_dp, &
1943  matrix_a=step, matrix_b=y, beta=0.0_dp, &
1944  matrix_c=s_y)
1945 
1946  CALL cp_fm_trace(step, y, s_dot_y)
1947 
1948  CALL cp_fm_trace(y, y, s_dot_y)
1949  CALL cp_fm_trace(step, step, s_dot_y)
1950 
1951  CALL cp_fm_trace(y, b_inv_y, y_dot_b_inv_y)
1952 
1953  factor1 = (s_dot_y + y_dot_b_inv_y)/(s_dot_y)**2
1954 
1955  CALL cp_fm_scale_and_add(1.0_dp, inv_hess, factor1, s_s)
1956 
1957  ! Second term
1958  CALL parallel_gemm(transa="N", transb="T", m=mat_size, n=mat_size, &
1959  k=1, alpha=1.0_dp, &
1960  matrix_a=b_inv_y, matrix_b=step, beta=0.0_dp, &
1961  matrix_c=b_inv_y_s)
1962 
1963  CALL parallel_gemm(transa="N", transb="N", m=mat_size, n=mat_size, &
1964  k=mat_size, alpha=1.0_dp, &
1965  matrix_a=s_y, matrix_b=prev_inv_hess, beta=0.0_dp, &
1966  matrix_c=s_y_b_inv)
1967 
1968  CALL cp_fm_scale_and_add(1.0_dp, b_inv_y_s, 1.0_dp, s_y_b_inv)
1969 
1970  ! Assemble the new inverse Hessian
1971  CALL cp_fm_scale_and_add(1.0_dp, inv_hess, -s_dot_y, b_inv_y_s)
1972 
1973  ! Deallocate intermediates
1974  CALL cp_fm_release(y)
1975  CALL cp_fm_release(b_inv_y)
1976  CALL cp_fm_release(s_s)
1977  CALL cp_fm_release(s_y)
1978  CALL cp_fm_release(b_inv_y_s)
1979  CALL cp_fm_release(s_y_b_inv)
1980 
1981  END SUBROUTINE inv_hessian_update
1982 
1983 ! **************************************************************************************************
1984 !> \brief ...
1985 !> \param grad ...
1986 !> \param prev_grad ...
1987 !> \param step ...
1988 !> \param prev_Hess ...
1989 !> \param Hess ...
1990 ! **************************************************************************************************
1991  SUBROUTINE hessian_update(grad, prev_grad, step, prev_Hess, Hess)
1992  TYPE(cp_fm_type), INTENT(IN) :: grad, prev_grad, step, prev_hess, hess
1993 
1994  INTEGER :: mat_size
1995  REAL(kind=dp) :: s_b_s, y_t_s
1996  TYPE(cp_blacs_env_type), POINTER :: blacs_env
1997  TYPE(cp_fm_struct_type), POINTER :: fm_struct_mat, fm_struct_vec, &
1998  fm_struct_vec_t
1999  TYPE(cp_fm_type) :: b_s, b_s_s_b, s_t_b, y, y_y_t
2000  TYPE(mp_para_env_type), POINTER :: para_env
2001 
2002  ! Recover the dimension
2003  CALL cp_fm_get_info(matrix=hess, &
2004  nrow_global=mat_size, para_env=para_env)
2005 
2006  CALL cp_fm_set_all(hess, 0.0_dp)
2007  CALL cp_fm_to_fm(prev_hess, hess)
2008 
2009  ! WARNING: our Hessian must be negative-definite, whereas BFGS makes it positive-definite!
2010  ! Therefore, we change sign in the beginning and in the end.
2011  CALL cp_fm_scale(-1.0_dp, hess)
2012 
2013  ! Create blacs environment
2014  NULLIFY (blacs_env)
2015  CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env)
2016 
2017  ! Get full matrix structures
2018  NULLIFY (fm_struct_mat, fm_struct_vec, fm_struct_vec_t)
2019 
2020  CALL cp_fm_get_info(matrix=prev_hess, &
2021  matrix_struct=fm_struct_mat)
2022  CALL cp_fm_get_info(matrix=grad, &
2023  matrix_struct=fm_struct_vec)
2024  CALL cp_fm_struct_create(fm_struct_vec_t, para_env=para_env, context=blacs_env, &
2025  nrow_global=1, ncol_global=mat_size)
2026 
2027  ! Allocate intermediates
2028  CALL cp_fm_create(b_s, fm_struct_vec, name="B_s")
2029  CALL cp_fm_create(s_t_b, fm_struct_vec_t, name="s_t_B")
2030  CALL cp_fm_create(y, fm_struct_vec, name="y")
2031 
2032  CALL cp_fm_create(y_y_t, fm_struct_mat, name="y_y_t")
2033  CALL cp_fm_create(b_s_s_b, fm_struct_mat, name="B_s_s_B")
2034 
2035  CALL cp_fm_set_all(y_y_t, 0.0_dp)
2036  CALL cp_fm_set_all(y, 0.0_dp)
2037  CALL cp_fm_set_all(b_s_s_b, 0.0_dp)
2038  CALL cp_fm_set_all(b_s, 0.0_dp)
2039  CALL cp_fm_set_all(s_t_b, 0.0_dp)
2040 
2041  ! Release the structure created only here
2042  CALL cp_fm_struct_release(fm_struct_vec_t)
2043 
2044  ! Calculate intermediates
2045  ! y the is gradient difference
2046  CALL cp_fm_to_fm(grad, y)
2047  CALL cp_fm_scale_and_add(1.0_dp, y, -1.0_dp, prev_grad)
2048 
2049  ! First term
2050  CALL parallel_gemm(transa="N", transb="T", m=mat_size, n=mat_size, &
2051  k=1, alpha=1.0_dp, &
2052  matrix_a=y, matrix_b=y, beta=0.0_dp, &
2053  matrix_c=y_y_t)
2054 
2055  CALL cp_fm_trace(y, step, y_t_s)
2056 
2057  CALL cp_fm_scale_and_add(1.0_dp, hess, (1.0_dp/y_t_s), y_y_t)
2058 
2059  ! Second term
2060  CALL parallel_gemm(transa="N", transb="N", m=mat_size, n=1, &
2061  k=mat_size, alpha=1.0_dp, &
2062  matrix_a=hess, matrix_b=step, beta=0.0_dp, &
2063  matrix_c=b_s)
2064 
2065  CALL cp_fm_trace(b_s, step, s_b_s)
2066 
2067  CALL parallel_gemm(transa="T", transb="N", m=1, n=mat_size, &
2068  k=mat_size, alpha=1.0_dp, &
2069  matrix_a=step, matrix_b=hess, beta=0.0_dp, &
2070  matrix_c=s_t_b)
2071 
2072  CALL parallel_gemm(transa="N", transb="N", m=mat_size, n=mat_size, &
2073  k=1, alpha=1.0_dp, &
2074  matrix_a=b_s, matrix_b=s_t_b, beta=0.0_dp, &
2075  matrix_c=b_s_s_b)
2076 
2077  CALL cp_fm_scale_and_add(1.0_dp, hess, -(1.0_dp/s_b_s), b_s_s_b)
2078 
2079  ! WARNING: our Hessian must be negative-definite, whereas BFGS makes it positive-definite!
2080  ! Therefore, we change sign in the beginning and in the end.
2081  CALL cp_fm_scale(-1.0_dp, hess)
2082 
2083  ! Release blacs environment
2084  CALL cp_blacs_env_release(blacs_env)
2085 
2086  ! Deallocate intermediates
2087  CALL cp_fm_release(y_y_t)
2088  CALL cp_fm_release(b_s_s_b)
2089  CALL cp_fm_release(b_s)
2090  CALL cp_fm_release(s_t_b)
2091  CALL cp_fm_release(y)
2092 
2093  END SUBROUTINE hessian_update
2094 
2095 ! **************************************************************************************************
2096 !> \brief ...
2097 !> \param grad ...
2098 !> \param prev_grad ...
2099 !> \param step ...
2100 !> \param prev_Hess ...
2101 !> \param Hess ...
2102 ! **************************************************************************************************
2103  SUBROUTINE symm_rank_one_update(grad, prev_grad, step, prev_Hess, Hess)
2104  TYPE(cp_fm_type), INTENT(IN) :: grad, prev_grad, step, prev_hess, hess
2105 
2106  INTEGER :: mat_size
2107  REAL(kind=dp) :: factor
2108  TYPE(cp_fm_struct_type), POINTER :: fm_struct_mat, fm_struct_vec
2109  TYPE(cp_fm_type) :: b_x, y, y_b_x_y_b_x
2110 
2111  ! Recover the dimension
2112  CALL cp_fm_get_info(matrix=hess, nrow_global=mat_size)
2113 
2114  CALL cp_fm_set_all(hess, 0.0_dp)
2115  CALL cp_fm_to_fm(prev_hess, hess)
2116 
2117  ! Get full matrix structures
2118  NULLIFY (fm_struct_mat, fm_struct_vec)
2119 
2120  CALL cp_fm_get_info(matrix=prev_hess, &
2121  matrix_struct=fm_struct_mat)
2122  CALL cp_fm_get_info(matrix=grad, &
2123  matrix_struct=fm_struct_vec)
2124 
2125  ! Allocate intermediates
2126  CALL cp_fm_create(y, fm_struct_vec, name="y")
2127  CALL cp_fm_create(b_x, fm_struct_vec, name="B_x")
2128  CALL cp_fm_create(y_b_x_y_b_x, fm_struct_mat, name="y_B_x_y_B_x")
2129 
2130  CALL cp_fm_set_all(y, 0.0_dp)
2131  CALL cp_fm_set_all(b_x, 0.0_dp)
2132  CALL cp_fm_set_all(y_b_x_y_b_x, 0.0_dp)
2133 
2134  ! Calculate intermediates
2135  ! y the is gradient difference
2136  CALL cp_fm_to_fm(grad, y)
2137  CALL cp_fm_scale_and_add(1.0_dp, y, -1.0_dp, prev_grad)
2138 
2139  CALL parallel_gemm(transa="N", transb="N", m=mat_size, n=1, &
2140  k=mat_size, alpha=1.0_dp, &
2141  matrix_a=hess, matrix_b=step, beta=0.0_dp, &
2142  matrix_c=b_x)
2143 
2144  CALL cp_fm_scale_and_add(1.0_dp, y, -1.0_dp, b_x)
2145 
2146  CALL parallel_gemm(transa="N", transb="T", m=mat_size, n=mat_size, &
2147  k=1, alpha=1.0_dp, &
2148  matrix_a=y, matrix_b=y, beta=0.0_dp, &
2149  matrix_c=y_b_x_y_b_x)
2150 
2151  ! Scaling factor
2152  CALL cp_fm_trace(y, step, factor)
2153 
2154  ! Assemble the Hessian
2155  CALL cp_fm_scale_and_add(1.0_dp, hess, (1.0_dp/factor), y_b_x_y_b_x)
2156 
2157  ! Deallocate intermediates
2158  CALL cp_fm_release(y)
2159  CALL cp_fm_release(b_x)
2160  CALL cp_fm_release(y_b_x_y_b_x)
2161 
2162  END SUBROUTINE symm_rank_one_update
2163 
2164 ! **************************************************************************************************
2165 !> \brief Controls the step, changes the trust radius if needed in maximization of the V_emb
2166 !> \param opt_embed ...
2167 !> \author Vladimir Rybkin
2168 ! **************************************************************************************************
2169  SUBROUTINE step_control(opt_embed)
2170  TYPE(opt_embed_pot_type) :: opt_embed
2171 
2172  CHARACTER(LEN=*), PARAMETER :: routinen = 'step_control'
2173 
2174  INTEGER :: handle
2175  REAL(kind=dp) :: actual_ener_change, ener_ratio, &
2176  lin_term, pred_ener_change, quad_term
2177  TYPE(cp_fm_struct_type), POINTER :: fm_struct
2178  TYPE(cp_fm_type) :: h_b
2179 
2180  CALL timeset(routinen, handle)
2181 
2182  NULLIFY (fm_struct)
2183  CALL cp_fm_get_info(matrix=opt_embed%embed_pot_grad, &
2184  matrix_struct=fm_struct)
2185  CALL cp_fm_create(h_b, fm_struct, name="H_b")
2186  CALL cp_fm_set_all(h_b, 0.0_dp)
2187 
2188  ! Calculate the quadratic estimate for the energy
2189  ! Linear term
2190  CALL cp_fm_trace(opt_embed%step, opt_embed%embed_pot_grad, lin_term)
2191 
2192  ! Quadratic term
2193  CALL parallel_gemm(transa="N", transb="N", m=opt_embed%dimen_aux, n=1, &
2194  k=opt_embed%dimen_aux, alpha=1.0_dp, &
2195  matrix_a=opt_embed%embed_pot_Hess, matrix_b=opt_embed%step, &
2196  beta=0.0_dp, matrix_c=h_b)
2197  CALL cp_fm_trace(opt_embed%step, h_b, quad_term)
2198 
2199  pred_ener_change = lin_term + 0.5_dp*quad_term
2200 
2201  ! Reveal actual energy change
2202  actual_ener_change = opt_embed%w_func(opt_embed%i_iter) - &
2203  opt_embed%w_func(opt_embed%last_accepted)
2204 
2205  ener_ratio = actual_ener_change/pred_ener_change
2206 
2207  CALL cp_fm_release(h_b)
2208 
2209  IF (actual_ener_change .GT. 0.0_dp) THEN ! If energy increases
2210  ! We accept step
2211  opt_embed%accept_step = .true.
2212  ! If energy change is larger than the predicted one, increase trust radius twice
2213  ! Else (between 0 and 1) leave as it is, unless Newton step has been taken and if the step is less than max
2214  IF ((ener_ratio .GT. 1.0_dp) .AND. (.NOT. opt_embed%newton_step) .AND. &
2215  (opt_embed%trust_rad .LT. opt_embed%max_trad)) &
2216  opt_embed%trust_rad = 2.0_dp*opt_embed%trust_rad
2217  ELSE ! Energy decreases
2218  ! If the decrease is not too large we allow this step to be taken
2219  ! Otherwise, the step is rejected
2220  IF (abs(actual_ener_change) .GE. opt_embed%allowed_decrease) THEN
2221  opt_embed%accept_step = .false.
2222  END IF
2223  ! Trust radius is decreased 4 times unless it's smaller than the minimal allowed value
2224  IF (opt_embed%trust_rad .GE. opt_embed%min_trad) &
2225  opt_embed%trust_rad = 0.25_dp*opt_embed%trust_rad
2226  END IF
2227 
2228  IF (opt_embed%accept_step) opt_embed%last_accepted = opt_embed%i_iter
2229 
2230  CALL timestop(handle)
2231 
2232  END SUBROUTINE step_control
2233 
2234 ! **************************************************************************************************
2235 !> \brief ...
2236 !> \param opt_embed ...
2237 !> \param diag_grad ...
2238 !> \param eigenval ...
2239 !> \param diag_step ...
2240 ! **************************************************************************************************
2241  SUBROUTINE level_shift(opt_embed, diag_grad, eigenval, diag_step)
2242  TYPE(opt_embed_pot_type) :: opt_embed
2243  TYPE(cp_fm_type), INTENT(IN) :: diag_grad
2244  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenval
2245  TYPE(cp_fm_type), INTENT(IN) :: diag_step
2246 
2247  CHARACTER(LEN=*), PARAMETER :: routinen = 'level_shift'
2248  INTEGER, PARAMETER :: max_iter = 25
2249  REAL(kind=dp), PARAMETER :: thresh = 0.00001_dp
2250 
2251  INTEGER :: handle, i_iter, l_global, lll, &
2252  min_index, nrow_local
2253  INTEGER, ALLOCATABLE, DIMENSION(:) :: red_eigenval_map
2254  INTEGER, DIMENSION(:), POINTER :: row_indices
2255  LOGICAL :: converged, do_shift
2256  REAL(kind=dp) :: diag_grad_norm, grad_min, hess_min, shift, shift_max, shift_min, step_len, &
2257  step_minus_trad, step_minus_trad_first, step_minus_trad_max, step_minus_trad_min
2258  TYPE(mp_para_env_type), POINTER :: para_env
2259 
2260  CALL timeset(routinen, handle)
2261 
2262  ! Array properties
2263  CALL cp_fm_get_info(matrix=opt_embed%embed_pot_coef, &
2264  nrow_local=nrow_local, &
2265  row_indices=row_indices, &
2266  para_env=para_env)
2267 
2268  min_index = minloc(abs(eigenval), dim=1)
2269  hess_min = eigenval(min_index)
2270  CALL cp_fm_get_element(diag_grad, min_index, 1, grad_min)
2271 
2272  CALL cp_fm_trace(diag_grad, diag_grad, diag_grad_norm)
2273 
2274  IF (hess_min .LT. 0.0_dp) THEN
2275  !shift_min = -2.0_dp*(diag_grad_norm/opt_embed%trust_rad - min(hess_min, 0.0_dp))
2276  !shift_max = max(0.0_dp, -hess_min + 0.5_dp*grad_min/opt_embed%trust_rad)
2277  !shift_max = MIN(-hess_min+0.5_dp*grad_min/opt_embed%trust_rad, 0.0_dp)
2278  shift_max = hess_min + 0.1
2279  shift_min = diag_grad_norm/opt_embed%trust_rad
2280  shift_min = 10.0_dp
2281  !If (abs(shift_max) .LE. thresh) then
2282  ! shift_min = -20.0_dp*(diag_grad_norm/opt_embed%trust_rad - min(hess_min, 0.0_dp))
2283  !Else
2284  ! shift_min = 20.0_dp*shift_max
2285  !Endif
2286 
2287  ! The boundary values
2288  step_minus_trad_max = shifted_step(diag_grad, eigenval, shift_max, opt_embed%trust_rad)
2289  step_minus_trad_min = shifted_step(diag_grad, eigenval, shift_min, opt_embed%trust_rad)
2290 
2291  ! Find zero by bisection
2292  converged = .false.
2293  do_shift = .false.
2294  IF (abs(step_minus_trad_max) .LE. thresh) THEN
2295  shift = shift_max
2296  ELSE
2297  IF (abs(step_minus_trad_min) .LE. thresh) THEN
2298  shift = shift_min
2299  ELSE
2300  DO i_iter = 1, max_iter
2301  shift = 0.5_dp*(shift_max + shift_min)
2302  step_minus_trad = shifted_step(diag_grad, eigenval, shift, opt_embed%trust_rad)
2303  IF (i_iter .EQ. 1) step_minus_trad_first = step_minus_trad
2304  IF (step_minus_trad .GT. 0.0_dp) shift_max = shift
2305  IF (step_minus_trad .LT. 0.0_dp) shift_min = shift
2306  !IF (ABS(shift_max-shift_min) .LT. thresh) converged = .TRUE.
2307  IF (abs(step_minus_trad) .LT. thresh) converged = .true.
2308  IF (converged) EXIT
2309  END DO
2310  IF (abs(step_minus_trad) .LT. abs(step_minus_trad_first)) do_shift = .true.
2311  END IF
2312  END IF
2313  ! Apply level-shifting
2314  IF (converged .OR. do_shift) THEN
2315  DO lll = 1, nrow_local
2316  l_global = row_indices(lll)
2317  IF (abs(eigenval(l_global)) .GE. thresh) THEN
2318  diag_step%local_data(lll, 1) = &
2319  -diag_grad%local_data(lll, 1)/(eigenval(l_global) - shift)
2320  ELSE
2321  diag_step%local_data(lll, 1) = 0.0_dp
2322  END IF
2323  END DO
2324  END IF
2325  IF (.NOT. converged) THEN ! Scale if shift has not been found
2326  CALL cp_fm_trace(diag_step, diag_step, step_len)
2327  CALL cp_fm_scale(opt_embed%trust_rad/step_len, diag_step)
2328  END IF
2329 
2330  ! Special case
2331  ELSE ! Hess min .LT. 0.0_dp
2332  ! First, find all positive eigenvalues
2333  ALLOCATE (red_eigenval_map(opt_embed%dimen_var_aux))
2334  red_eigenval_map = 0
2335  DO lll = 1, nrow_local
2336  l_global = row_indices(lll)
2337  IF (eigenval(l_global) .GE. 0.0_dp) THEN
2338  red_eigenval_map(l_global) = 1
2339  END IF
2340  END DO
2341  CALL para_env%sum(red_eigenval_map)
2342 
2343  ! Set shift as -hess_min and find step on the reduced space of negative-value
2344  ! eigenvectors
2345  shift = -hess_min
2346  DO lll = 1, nrow_local
2347  l_global = row_indices(lll)
2348  IF (red_eigenval_map(l_global) .EQ. 0) THEN
2349  IF (abs(eigenval(l_global)) .GE. thresh) THEN
2350  diag_step%local_data(lll, 1) = &
2351  -diag_grad%local_data(lll, 1)/(eigenval(l_global) - shift)
2352  ELSE
2353  diag_step%local_data(lll, 1) = 0.0_dp
2354  END IF
2355  ELSE
2356  diag_step%local_data(lll, 1) = 0.0_dp
2357  END IF
2358  END DO
2359 
2360  ! Find the step length of such a step
2361  CALL cp_fm_trace(diag_step, diag_step, step_len)
2362 
2363  END IF
2364 
2365  CALL timestop(handle)
2366 
2367  END SUBROUTINE level_shift
2368 
2369 ! **************************************************************************************************
2370 !> \brief ...
2371 !> \param diag_grad ...
2372 !> \param eigenval ...
2373 !> \param shift ...
2374 !> \param trust_rad ...
2375 !> \return ...
2376 ! **************************************************************************************************
2377  FUNCTION shifted_step(diag_grad, eigenval, shift, trust_rad) RESULT(step_minus_trad)
2378  TYPE(cp_fm_type), INTENT(IN) :: diag_grad
2379  REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
2380  INTENT(IN) :: eigenval
2381  REAL(kind=dp), INTENT(IN) :: shift, trust_rad
2382  REAL(kind=dp) :: step_minus_trad
2383 
2384  REAL(kind=dp), PARAMETER :: thresh = 0.000001_dp
2385 
2386  INTEGER :: l_global, lll, nrow_local
2387  INTEGER, DIMENSION(:), POINTER :: row_indices
2388  REAL(kind=dp) :: step, step_1d
2389  TYPE(mp_para_env_type), POINTER :: para_env
2390 
2391  CALL cp_fm_get_info(matrix=diag_grad, &
2392  nrow_local=nrow_local, &
2393  row_indices=row_indices, &
2394  para_env=para_env)
2395 
2396  step = 0.0_dp
2397  DO lll = 1, nrow_local
2398  l_global = row_indices(lll)
2399  IF ((abs(eigenval(l_global)) .GE. thresh) .AND. (abs(diag_grad%local_data(lll, 1)) .GE. thresh)) THEN
2400  step_1d = -diag_grad%local_data(lll, 1)/(eigenval(l_global) + shift)
2401  step = step + step_1d**2
2402  END IF
2403  END DO
2404 
2405  CALL para_env%sum(step)
2406 
2407  step_minus_trad = sqrt(step) - trust_rad
2408 
2409  END FUNCTION shifted_step
2410 
2411 ! **************************************************************************************************
2412 !> \brief ...
2413 !> \param step ...
2414 !> \param prev_step ...
2415 !> \param grad ...
2416 !> \param prev_grad ...
2417 !> \return ...
2418 !> \retval length ...
2419 ! **************************************************************************************************
2420  FUNCTION barzilai_borwein(step, prev_step, grad, prev_grad) RESULT(length)
2421  TYPE(cp_fm_type), INTENT(IN) :: step, prev_step, grad, prev_grad
2422  REAL(kind=dp) :: length
2423 
2424  REAL(kind=dp) :: denominator, numerator
2425  TYPE(cp_fm_struct_type), POINTER :: fm_struct
2426  TYPE(cp_fm_type) :: grad_diff, step_diff
2427 
2428  ! Get full matrix structures
2429  NULLIFY (fm_struct)
2430 
2431  CALL cp_fm_get_info(matrix=grad, &
2432  matrix_struct=fm_struct)
2433 
2434  ! Allocate intermediates
2435  CALL cp_fm_create(grad_diff, fm_struct, name="grad_diff")
2436  CALL cp_fm_create(step_diff, fm_struct, name="step_diff")
2437 
2438  ! Calculate intermediates
2439  CALL cp_fm_to_fm(grad, grad_diff)
2440  CALL cp_fm_to_fm(step, step_diff)
2441 
2442  CALL cp_fm_scale_and_add(1.0_dp, grad_diff, -1.0_dp, prev_grad)
2443  CALL cp_fm_scale_and_add(1.0_dp, step_diff, -1.0_dp, prev_step)
2444 
2445  CALL cp_fm_trace(step_diff, grad_diff, numerator)
2446  CALL cp_fm_trace(grad_diff, grad_diff, denominator)
2447 
2448  ! Release intermediates
2449  CALL cp_fm_release(grad_diff)
2450  CALL cp_fm_release(step_diff)
2451 
2452  length = numerator/denominator
2453 
2454  END FUNCTION barzilai_borwein
2455 
2456 ! **************************************************************************************************
2457 !> \brief ...
2458 !> \param pw_env ...
2459 !> \param embed_pot ...
2460 !> \param spin_embed_pot ...
2461 !> \param diff_rho_r ...
2462 !> \param diff_rho_spin ...
2463 !> \param rho_r_ref ...
2464 !> \param open_shell_embed ...
2465 !> \param step_len ...
2466 ! **************************************************************************************************
2467  SUBROUTINE leeuwen_baerends_potential_update(pw_env, embed_pot, spin_embed_pot, diff_rho_r, diff_rho_spin, &
2468  rho_r_ref, open_shell_embed, step_len)
2469  TYPE(pw_env_type), POINTER :: pw_env
2470  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: embed_pot
2471  TYPE(pw_r3d_rs_type), INTENT(IN), POINTER :: spin_embed_pot
2472  TYPE(pw_r3d_rs_type), INTENT(IN) :: diff_rho_r, diff_rho_spin
2473  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_ref
2474  LOGICAL, INTENT(IN) :: open_shell_embed
2475  REAL(kind=dp), INTENT(IN) :: step_len
2476 
2477  CHARACTER(LEN=*), PARAMETER :: routinen = 'Leeuwen_Baerends_potential_update'
2478 
2479  INTEGER :: handle, i, i_spin, j, k, nspins
2480  INTEGER, DIMENSION(3) :: lb, ub
2481  REAL(kind=dp) :: my_rho, rho_cutoff
2482  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
2483  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: new_embed_pot, rho_n_1, temp_embed_pot
2484 
2485  CALL timeset(routinen, handle)
2486 
2487  rho_cutoff = epsilon(0.0_dp)
2488 
2489  ! Prepare plane-waves pool
2490  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
2491  NULLIFY (new_embed_pot)
2492 
2493  nspins = 1
2494  IF (open_shell_embed) nspins = 2
2495  NULLIFY (new_embed_pot)
2496  ALLOCATE (new_embed_pot(nspins))
2497  DO i_spin = 1, nspins
2498  CALL auxbas_pw_pool%create_pw(new_embed_pot(i_spin))
2499  CALL pw_zero(new_embed_pot(i_spin))
2500  END DO
2501 
2502  lb(1:3) = embed_pot%pw_grid%bounds_local(1, 1:3)
2503  ub(1:3) = embed_pot%pw_grid%bounds_local(2, 1:3)
2504 
2505  IF (.NOT. open_shell_embed) THEN
2506 !$OMP PARALLEL DO DEFAULT(NONE) &
2507 !$OMP PRIVATE(i,j,k, my_rho) &
2508 !$OMP SHARED(new_embed_pot, embed_pot, diff_rho_r, rho_r_ref, lb, ub, rho_cutoff, step_len)
2509  DO k = lb(3), ub(3)
2510  DO j = lb(2), ub(2)
2511  DO i = lb(1), ub(1)
2512  IF (rho_r_ref(1)%array(i, j, k) .GT. rho_cutoff) THEN
2513  my_rho = rho_r_ref(1)%array(i, j, k)
2514  ELSE
2515  my_rho = rho_cutoff
2516  END IF
2517  new_embed_pot(1)%array(i, j, k) = step_len*embed_pot%array(i, j, k)* &
2518  (diff_rho_r%array(i, j, k) + rho_r_ref(1)%array(i, j, k))/my_rho
2519  END DO
2520  END DO
2521  END DO
2522 !$OMP END PARALLEL DO
2523  CALL pw_copy(new_embed_pot(1), embed_pot)
2524 
2525  ELSE
2526  ! One has to work with spin components rather than with total and spin density
2527  NULLIFY (rho_n_1)
2528  ALLOCATE (rho_n_1(nspins))
2529  NULLIFY (temp_embed_pot)
2530  ALLOCATE (temp_embed_pot(nspins))
2531  DO i_spin = 1, nspins
2532  CALL auxbas_pw_pool%create_pw(rho_n_1(i_spin))
2533  CALL pw_zero(rho_n_1(i_spin))
2534  CALL auxbas_pw_pool%create_pw(temp_embed_pot(i_spin))
2535  CALL pw_zero(temp_embed_pot(i_spin))
2536  END DO
2537  CALL pw_copy(diff_rho_r, rho_n_1(1))
2538  CALL pw_copy(diff_rho_r, rho_n_1(2))
2539  CALL pw_axpy(diff_rho_spin, rho_n_1(1), 1.0_dp)
2540  CALL pw_axpy(diff_rho_spin, rho_n_1(2), -1.0_dp)
2541  CALL pw_scale(rho_n_1(1), a=0.5_dp)
2542  CALL pw_scale(rho_n_1(2), a=0.5_dp)
2543 
2544  CALL pw_copy(embed_pot, temp_embed_pot(1))
2545  CALL pw_copy(embed_pot, temp_embed_pot(2))
2546  CALL pw_axpy(spin_embed_pot, temp_embed_pot(1), 1.0_dp)
2547  CALL pw_axpy(spin_embed_pot, temp_embed_pot(2), -1.0_dp)
2548 
2549  IF (SIZE(rho_r_ref) .EQ. 2) THEN
2550  CALL pw_axpy(rho_r_ref(1), rho_n_1(1), 1.0_dp)
2551  CALL pw_axpy(rho_r_ref(2), rho_n_1(2), 1.0_dp)
2552 
2553 !$OMP PARALLEL DO DEFAULT(NONE) &
2554 !$OMP PRIVATE(i,j,k, my_rho) &
2555 !$OMP SHARED(new_embed_pot, temp_embed_pot, rho_r_ref, rho_n_1, lb, ub, rho_cutoff, step_len)
2556  DO k = lb(3), ub(3)
2557  DO j = lb(2), ub(2)
2558  DO i = lb(1), ub(1)
2559  IF (rho_r_ref(1)%array(i, j, k) .GT. rho_cutoff) THEN
2560  my_rho = rho_r_ref(1)%array(i, j, k)
2561  ELSE
2562  my_rho = rho_cutoff
2563  END IF
2564  new_embed_pot(1)%array(i, j, k) = step_len*temp_embed_pot(1)%array(i, j, k)* &
2565  (rho_n_1(1)%array(i, j, k))/my_rho
2566  IF (rho_r_ref(2)%array(i, j, k) .GT. rho_cutoff) THEN
2567  my_rho = rho_r_ref(2)%array(i, j, k)
2568  ELSE
2569  my_rho = rho_cutoff
2570  END IF
2571  new_embed_pot(2)%array(i, j, k) = step_len*temp_embed_pot(2)%array(i, j, k)* &
2572  (rho_n_1(2)%array(i, j, k))/my_rho
2573  END DO
2574  END DO
2575  END DO
2576 !$OMP END PARALLEL DO
2577 
2578  ELSE ! Reference system is closed-shell
2579  CALL pw_axpy(rho_r_ref(1), rho_n_1(1), 1.0_dp)
2580  ! The beta spin component is here equal to the difference: nothing to do
2581 
2582 !$OMP PARALLEL DO DEFAULT(NONE) &
2583 !$OMP PRIVATE(i,j,k, my_rho) &
2584 !$OMP SHARED(new_embed_pot, rho_r_ref, temp_embed_pot, rho_n_1, lb, ub, rho_cutoff, step_len)
2585  DO k = lb(3), ub(3)
2586  DO j = lb(2), ub(2)
2587  DO i = lb(1), ub(1)
2588  IF (rho_r_ref(1)%array(i, j, k) .GT. rho_cutoff) THEN
2589  my_rho = 0.5_dp*rho_r_ref(1)%array(i, j, k)
2590  ELSE
2591  my_rho = rho_cutoff
2592  END IF
2593  new_embed_pot(1)%array(i, j, k) = step_len*temp_embed_pot(1)%array(i, j, k)* &
2594  (rho_n_1(1)%array(i, j, k))/my_rho
2595  new_embed_pot(2)%array(i, j, k) = step_len*temp_embed_pot(2)%array(i, j, k)* &
2596  (rho_n_1(2)%array(i, j, k))/my_rho
2597  END DO
2598  END DO
2599  END DO
2600 !$OMP END PARALLEL DO
2601  END IF
2602 
2603  CALL pw_copy(new_embed_pot(1), embed_pot)
2604  CALL pw_axpy(new_embed_pot(2), embed_pot, 1.0_dp)
2605  CALL pw_scale(embed_pot, a=0.5_dp)
2606  CALL pw_copy(new_embed_pot(1), spin_embed_pot)
2607  CALL pw_axpy(new_embed_pot(2), spin_embed_pot, -1.0_dp)
2608  CALL pw_scale(spin_embed_pot, a=0.5_dp)
2609 
2610  DO i_spin = 1, nspins
2611  CALL rho_n_1(i_spin)%release()
2612  CALL temp_embed_pot(i_spin)%release()
2613  END DO
2614  DEALLOCATE (rho_n_1)
2615  DEALLOCATE (temp_embed_pot)
2616  END IF
2617 
2618  DO i_spin = 1, nspins
2619  CALL new_embed_pot(i_spin)%release()
2620  END DO
2621 
2622  DEALLOCATE (new_embed_pot)
2623 
2624  CALL timestop(handle)
2625 
2626  END SUBROUTINE leeuwen_baerends_potential_update
2627 
2628 ! **************************************************************************************************
2629 !> \brief ...
2630 !> \param qs_env ...
2631 !> \param rho_r_ref ...
2632 !> \param prev_embed_pot ...
2633 !> \param prev_spin_embed_pot ...
2634 !> \param embed_pot ...
2635 !> \param spin_embed_pot ...
2636 !> \param diff_rho_r ...
2637 !> \param diff_rho_spin ...
2638 !> \param v_w_ref ...
2639 !> \param i_iter ...
2640 !> \param step_len ...
2641 !> \param open_shell_embed ...
2642 !> \param vw_cutoff ...
2643 !> \param vw_smooth_cutoff_range ...
2644 ! **************************************************************************************************
2645  SUBROUTINE fab_update(qs_env, rho_r_ref, prev_embed_pot, prev_spin_embed_pot, embed_pot, spin_embed_pot, &
2646  diff_rho_r, diff_rho_spin, v_w_ref, i_iter, step_len, open_shell_embed, &
2647  vw_cutoff, vw_smooth_cutoff_range)
2648  TYPE(qs_environment_type), POINTER :: qs_env
2649  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_ref
2650  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: prev_embed_pot
2651  TYPE(pw_r3d_rs_type), INTENT(IN), POINTER :: prev_spin_embed_pot
2652  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: embed_pot
2653  TYPE(pw_r3d_rs_type), INTENT(IN), POINTER :: spin_embed_pot
2654  TYPE(pw_r3d_rs_type), INTENT(IN) :: diff_rho_r, diff_rho_spin
2655  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_w_ref
2656  INTEGER, INTENT(IN) :: i_iter
2657  REAL(kind=dp) :: step_len
2658  LOGICAL :: open_shell_embed
2659  REAL(kind=dp) :: vw_cutoff, vw_smooth_cutoff_range
2660 
2661  CHARACTER(LEN=*), PARAMETER :: routinen = 'FAB_update'
2662 
2663  INTEGER :: handle, i_spin, nspins
2664  TYPE(pw_env_type), POINTER :: pw_env
2665  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
2666  TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: new_embed_pot, temp_embed_pot, v_w
2667  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: curr_rho
2668 
2669  CALL timeset(routinen, handle)
2670 
2671  ! Update formula: v(n+1) = v(n-1) - v_w(ref) + v_w(n)
2672 
2673  CALL get_qs_env(qs_env=qs_env, &
2674  pw_env=pw_env)
2675  ! Get plane waves pool
2676  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
2677 
2678  ! We calculate von Weizsaecker potential for the reference density
2679  ! only at the first iteration
2680  IF (i_iter .LE. 1) THEN
2681  nspins = SIZE(rho_r_ref)
2682  NULLIFY (v_w_ref)
2683  ALLOCATE (v_w_ref(nspins))
2684  DO i_spin = 1, nspins
2685  CALL auxbas_pw_pool%create_pw(v_w_ref(i_spin))
2686  END DO
2687  CALL von_weizsacker(rho_r_ref, v_w_ref, qs_env, vw_cutoff, vw_smooth_cutoff_range)
2688  ! For the first step previous are set to current
2689  CALL pw_copy(embed_pot, prev_embed_pot)
2690  CALL pw_axpy(diff_rho_r, embed_pot, 0.5_dp)
2691  IF (open_shell_embed) THEN
2692  CALL pw_copy(spin_embed_pot, prev_spin_embed_pot)
2693  CALL pw_axpy(diff_rho_r, embed_pot, 0.5_dp)
2694  END IF
2695 
2696  ELSE
2697 
2698  ! Reference can be closed shell, but total embedding - open shell:
2699  ! redefine nspins
2700  nspins = 1
2701  IF (open_shell_embed) nspins = 2
2702  ALLOCATE (new_embed_pot(nspins))
2703  ALLOCATE (v_w(nspins))
2704  NULLIFY (curr_rho)
2705  ALLOCATE (curr_rho(nspins))
2706  DO i_spin = 1, nspins
2707  CALL auxbas_pw_pool%create_pw(new_embed_pot(i_spin))
2708  CALL pw_zero(new_embed_pot(i_spin))
2709 
2710  CALL auxbas_pw_pool%create_pw(v_w(i_spin))
2711  CALL pw_zero(v_w(i_spin))
2712 
2713  CALL auxbas_pw_pool%create_pw(curr_rho(i_spin))
2714  CALL pw_zero(curr_rho(i_spin))
2715  END DO
2716 
2717  ! Now, deal with the current density
2718 
2719  IF (.NOT. open_shell_embed) THEN
2720  ! Reconstruct current density
2721  CALL pw_copy(diff_rho_r, curr_rho(1))
2722  CALL pw_axpy(rho_r_ref(1), curr_rho(1), 1.0_dp)
2723  ! Compute von Weizsaecker potential
2724  CALL von_weizsacker(curr_rho, v_w, qs_env, vw_cutoff, vw_smooth_cutoff_range)
2725  ! Compute new embedding potential
2726  CALL pw_copy(prev_embed_pot, new_embed_pot(1))
2727  CALL pw_axpy(v_w(1), new_embed_pot(1), step_len)
2728  CALL pw_axpy(v_w_ref(1), new_embed_pot(1), -step_len)
2729  ! Copy the potentials
2730 
2731  CALL pw_copy(embed_pot, prev_embed_pot)
2732  CALL pw_copy(new_embed_pot(1), embed_pot)
2733 
2734  ELSE
2735  ! Reconstruct current density
2736  CALL pw_copy(diff_rho_r, curr_rho(1))
2737  CALL pw_copy(diff_rho_r, curr_rho(2))
2738  CALL pw_axpy(diff_rho_spin, curr_rho(1), 1.0_dp)
2739  CALL pw_axpy(diff_rho_spin, curr_rho(2), -1.0_dp)
2740  CALL pw_scale(curr_rho(1), a=0.5_dp)
2741  CALL pw_scale(curr_rho(2), a=0.5_dp)
2742 
2743  IF (SIZE(rho_r_ref) .EQ. 1) THEN ! If reference system is closed-shell
2744  CALL pw_axpy(rho_r_ref(1), curr_rho(1), 0.5_dp)
2745  CALL pw_axpy(rho_r_ref(1), curr_rho(2), 0.5_dp)
2746  ELSE ! If reference system is open-shell
2747  CALL pw_axpy(rho_r_ref(1), curr_rho(1), 1.0_dp)
2748  CALL pw_axpy(rho_r_ref(2), curr_rho(2), 1.0_dp)
2749  END IF
2750 
2751  ! Compute von Weizsaecker potential
2752  CALL von_weizsacker(curr_rho, v_w, qs_env, vw_cutoff, vw_smooth_cutoff_range)
2753 
2754  ! Reconstruct corrent spin components of the potential
2755  ALLOCATE (temp_embed_pot(nspins))
2756  DO i_spin = 1, nspins
2757  CALL auxbas_pw_pool%create_pw(temp_embed_pot(i_spin))
2758  CALL pw_zero(temp_embed_pot(i_spin))
2759  END DO
2760  CALL pw_copy(embed_pot, temp_embed_pot(1))
2761  CALL pw_copy(embed_pot, temp_embed_pot(2))
2762  CALL pw_axpy(spin_embed_pot, temp_embed_pot(1), 1.0_dp)
2763  CALL pw_axpy(spin_embed_pot, temp_embed_pot(2), -1.0_dp)
2764 
2765  ! Compute new embedding potential
2766  IF (SIZE(v_w_ref) .EQ. 1) THEN ! Reference system is closed-shell
2767  CALL pw_copy(temp_embed_pot(1), new_embed_pot(1))
2768  CALL pw_axpy(v_w(1), new_embed_pot(1), 0.5_dp*step_len)
2769  CALL pw_axpy(v_w_ref(1), new_embed_pot(1), -0.5_dp*step_len)
2770 
2771  CALL pw_copy(temp_embed_pot(2), new_embed_pot(2))
2772  CALL pw_axpy(v_w(2), new_embed_pot(2), 0.5_dp)
2773  CALL pw_axpy(v_w_ref(1), new_embed_pot(2), -0.5_dp)
2774 
2775  ELSE ! Reference system is open-shell
2776 
2777  DO i_spin = 1, nspins
2778  CALL pw_copy(temp_embed_pot(i_spin), new_embed_pot(i_spin))
2779  CALL pw_axpy(v_w(1), new_embed_pot(i_spin), step_len)
2780  CALL pw_axpy(v_w_ref(i_spin), new_embed_pot(i_spin), -step_len)
2781  END DO
2782  END IF
2783 
2784  ! Update embedding potentials
2785  CALL pw_copy(embed_pot, prev_embed_pot)
2786  CALL pw_copy(spin_embed_pot, prev_spin_embed_pot)
2787 
2788  CALL pw_copy(new_embed_pot(1), embed_pot)
2789  CALL pw_axpy(new_embed_pot(2), embed_pot, 1.0_dp)
2790  CALL pw_scale(embed_pot, a=0.5_dp)
2791  CALL pw_copy(new_embed_pot(1), spin_embed_pot)
2792  CALL pw_axpy(new_embed_pot(2), spin_embed_pot, -1.0_dp)
2793  CALL pw_scale(spin_embed_pot, a=0.5_dp)
2794 
2795  DO i_spin = 1, nspins
2796  CALL temp_embed_pot(i_spin)%release()
2797  END DO
2798  DEALLOCATE (temp_embed_pot)
2799 
2800  END IF
2801 
2802  DO i_spin = 1, nspins
2803  CALL curr_rho(i_spin)%release()
2804  CALL new_embed_pot(i_spin)%release()
2805  CALL v_w(i_spin)%release()
2806  END DO
2807 
2808  DEALLOCATE (new_embed_pot)
2809  DEALLOCATE (v_w)
2810  DEALLOCATE (curr_rho)
2811 
2812  END IF
2813 
2814  CALL timestop(handle)
2815 
2816  END SUBROUTINE fab_update
2817 
2818 ! **************************************************************************************************
2819 !> \brief ...
2820 !> \param rho_r ...
2821 !> \param v_w ...
2822 !> \param qs_env ...
2823 !> \param vw_cutoff ...
2824 !> \param vw_smooth_cutoff_range ...
2825 ! **************************************************************************************************
2826  SUBROUTINE von_weizsacker(rho_r, v_w, qs_env, vw_cutoff, vw_smooth_cutoff_range)
2827  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
2828  TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN) :: v_w
2829  TYPE(qs_environment_type), POINTER :: qs_env
2830  REAL(kind=dp), INTENT(IN) :: vw_cutoff, vw_smooth_cutoff_range
2831 
2832  REAL(kind=dp), PARAMETER :: one_4 = 0.25_dp, one_8 = 0.125_dp
2833 
2834  INTEGER :: i, i_spin, j, k, nspins
2835  INTEGER, DIMENSION(3) :: lb, ub
2836  REAL(kind=dp) :: density_smooth_cut_range, my_rho, &
2837  rho_cutoff
2838  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: rhoa, rhob
2839  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
2840  TYPE(pw_env_type), POINTER :: pw_env
2841  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
2842  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: tau
2843  TYPE(section_vals_type), POINTER :: input, xc_section
2844  TYPE(xc_rho_cflags_type) :: needs
2845  TYPE(xc_rho_set_type) :: rho_set
2846 
2847  rho_cutoff = epsilon(0.0_dp)
2848 
2849  nspins = SIZE(rho_r)
2850 
2851  NULLIFY (xc_section)
2852 
2853  CALL get_qs_env(qs_env=qs_env, &
2854  pw_env=pw_env, &
2855  input=input)
2856 
2857  ! Get plane waves pool
2858  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
2859 
2860  ! get some of the grids ready
2861  NULLIFY (rho_g)
2862  ALLOCATE (rho_g(nspins))
2863  DO i_spin = 1, nspins
2864  CALL auxbas_pw_pool%create_pw(rho_g(i_spin))
2865  CALL pw_transfer(rho_r(i_spin), rho_g(i_spin))
2866  END DO
2867 
2868  xc_section => section_vals_get_subs_vals(input, "DFT%XC")
2869 
2870  CALL xc_rho_set_create(rho_set, &
2871  rho_r(1)%pw_grid%bounds_local, &
2872  rho_cutoff=section_get_rval(xc_section, "density_cutoff"), &
2873  drho_cutoff=section_get_rval(xc_section, "gradient_cutoff"), &
2874  tau_cutoff=section_get_rval(xc_section, "tau_cutoff"))
2875 
2876  CALL xc_rho_cflags_setall(needs, .false.)
2877 
2878  IF (nspins .EQ. 2) THEN
2879  needs%rho_spin = .true.
2880  needs%norm_drho_spin = .true.
2881  needs%laplace_rho_spin = .true.
2882  ELSE
2883  needs%rho = .true.
2884  needs%norm_drho = .true.
2885  needs%laplace_rho = .true.
2886  END IF
2887 
2888  CALL xc_rho_set_update(rho_set, rho_r, rho_g, tau, needs, &
2889  section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
2890  section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
2891  auxbas_pw_pool)
2892 
2893  CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", &
2894  r_val=rho_cutoff)
2895  CALL section_vals_val_get(xc_section, "DENSITY_SMOOTH_CUTOFF_RANGE", &
2896  r_val=density_smooth_cut_range)
2897 
2898  lb(1:3) = rho_r(1)%pw_grid%bounds_local(1, 1:3)
2899  ub(1:3) = rho_r(1)%pw_grid%bounds_local(2, 1:3)
2900 
2901  IF (nspins .EQ. 2) THEN
2902 !$OMP PARALLEL DO DEFAULT(NONE) &
2903 !$OMP PRIVATE(i,j,k, my_rho) &
2904 !$OMP SHARED(v_w, rho_r, rho_set, lb, ub, rho_cutoff)
2905  DO k = lb(3), ub(3)
2906  DO j = lb(2), ub(2)
2907  DO i = lb(1), ub(1)
2908  IF (rho_r(1)%array(i, j, k) .GT. rho_cutoff) THEN
2909  my_rho = rho_r(1)%array(i, j, k)
2910  ELSE
2911  my_rho = rho_cutoff
2912  END IF
2913  v_w(1)%array(i, j, k) = one_8*rho_set%norm_drhoa(i, j, k)**2/my_rho**2 - &
2914  one_4*rho_set%laplace_rhoa(i, j, k)/my_rho
2915 
2916  IF (rho_r(2)%array(i, j, k) .GT. rho_cutoff) THEN
2917  my_rho = rho_r(2)%array(i, j, k)
2918  ELSE
2919  my_rho = rho_cutoff
2920  END IF
2921  v_w(2)%array(i, j, k) = one_8*rho_set%norm_drhob(i, j, k)**2/my_rho**2 - &
2922  one_4*rho_set%laplace_rhob(i, j, k)/my_rho
2923  END DO
2924  END DO
2925  END DO
2926 !$OMP END PARALLEL DO
2927  ELSE
2928 !$OMP PARALLEL DO DEFAULT(NONE) &
2929 !$OMP PRIVATE(i,j,k, my_rho) &
2930 !$OMP SHARED(v_w, rho_r, rho_set, lb, ub, rho_cutoff)
2931  DO k = lb(3), ub(3)
2932  DO j = lb(2), ub(2)
2933  DO i = lb(1), ub(1)
2934  IF (rho_r(1)%array(i, j, k) .GT. rho_cutoff) THEN
2935  my_rho = rho_r(1)%array(i, j, k)
2936  v_w(1)%array(i, j, k) = one_8*rho_set%norm_drho(i, j, k)**2/my_rho**2 - &
2937  one_4*rho_set%laplace_rho(i, j, k)/my_rho
2938  ELSE
2939  v_w(1)%array(i, j, k) = 0.0_dp
2940  END IF
2941  END DO
2942  END DO
2943  END DO
2944 !$OMP END PARALLEL DO
2945 
2946  END IF
2947 
2948  ! Smoothen the von Weizsaecker potential
2949  IF (nspins == 2) THEN
2950  density_smooth_cut_range = 0.5_dp*density_smooth_cut_range
2951  rho_cutoff = 0.5_dp*rho_cutoff
2952  END IF
2953  DO i_spin = 1, nspins
2954  CALL smooth_cutoff(pot=v_w(i_spin)%array, rho=rho_r(i_spin)%array, rhoa=rhoa, rhob=rhob, &
2955  rho_cutoff=vw_cutoff, &
2956  rho_smooth_cutoff_range=vw_smooth_cutoff_range)
2957  END DO
2958 
2959  CALL xc_rho_set_release(rho_set, pw_pool=auxbas_pw_pool)
2960 
2961  DO i_spin = 1, nspins
2962  CALL rho_g(i_spin)%release()
2963  END DO
2964  DEALLOCATE (rho_g)
2965 
2966  END SUBROUTINE von_weizsacker
2967 
2968 ! **************************************************************************************************
2969 !> \brief ...
2970 !> \param diff_rho_r ...
2971 !> \return ...
2972 ! **************************************************************************************************
2973  FUNCTION max_dens_diff(diff_rho_r) RESULT(total_max_diff)
2974  TYPE(pw_r3d_rs_type), INTENT(IN) :: diff_rho_r
2975  REAL(kind=dp) :: total_max_diff
2976 
2977  INTEGER :: size_x, size_y, size_z
2978  REAL(kind=dp) :: max_diff
2979  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: grid_3d
2980 
2981  !, i_x, i_y, i_z
2982 
2983  ! Get the sizes
2984  size_x = SIZE(diff_rho_r%array, 1)
2985  size_y = SIZE(diff_rho_r%array, 2)
2986  size_z = SIZE(diff_rho_r%array, 3)
2987 
2988  ! Allocate the density
2989  ALLOCATE (grid_3d(size_x, size_y, size_z))
2990 
2991  ! Copy density
2992  grid_3d(:, :, :) = diff_rho_r%array(:, :, :)
2993 
2994  ! Find the maximum absolute value
2995  max_diff = maxval(abs(grid_3d))
2996  total_max_diff = max_diff
2997  CALL diff_rho_r%pw_grid%para%group%max(total_max_diff)
2998 
2999  ! Deallocate the density
3000  DEALLOCATE (grid_3d)
3001 
3002  END FUNCTION max_dens_diff
3003 
3004 ! **************************************************************************************************
3005 !> \brief Prints a cube for the (rho_A + rho_B - rho_ref) to be minimized in embedding
3006 !> \param diff_rho_r ...
3007 !> \param i_iter ...
3008 !> \param qs_env ...
3009 !> \param final_one ...
3010 !> \author Vladimir Rybkin
3011 ! **************************************************************************************************
3012  SUBROUTINE print_rho_diff(diff_rho_r, i_iter, qs_env, final_one)
3013  TYPE(pw_r3d_rs_type), INTENT(IN) :: diff_rho_r
3014  INTEGER, INTENT(IN) :: i_iter
3015  TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
3016  LOGICAL, INTENT(IN) :: final_one
3017 
3018  CHARACTER(LEN=default_path_length) :: filename, my_pos_cube, title
3019  INTEGER :: unit_nr
3020  TYPE(cp_logger_type), POINTER :: logger
3021  TYPE(particle_list_type), POINTER :: particles
3022  TYPE(qs_subsys_type), POINTER :: subsys
3023  TYPE(section_vals_type), POINTER :: dft_section, input
3024 
3025  NULLIFY (subsys, input)
3026 
3027  CALL get_qs_env(qs_env=qs_env, &
3028  subsys=subsys, &
3029  input=input)
3030  dft_section => section_vals_get_subs_vals(input, "DFT")
3031  CALL qs_subsys_get(subsys, particles=particles)
3032 
3033  logger => cp_get_default_logger()
3034  IF (btest(cp_print_key_should_output(logger%iter_info, input, &
3035  "DFT%QS%OPT_EMBED%EMBED_DENS_DIFF"), cp_p_file)) THEN
3036  my_pos_cube = "REWIND"
3037  IF (.NOT. final_one) THEN
3038  WRITE (filename, '(a5,I3.3,a1,I1.1)') "DIFF_", i_iter
3039  ELSE
3040  WRITE (filename, '(a5,I3.3,a1,I1.1)') "DIFF"
3041  END IF
3042  unit_nr = cp_print_key_unit_nr(logger, input, "DFT%QS%OPT_EMBED%EMBED_DENS_DIFF", &
3043  extension=".cube", middle_name=trim(filename), file_position=my_pos_cube, &
3044  log_filename=.false.)
3045 
3046  WRITE (title, *) "EMBEDDING DENSITY DIFFERENCE ", " optimization step ", i_iter
3047  CALL cp_pw_to_cube(diff_rho_r, unit_nr, title, particles=particles, &
3048  stride=section_get_ivals(dft_section, "QS%OPT_EMBED%EMBED_DENS_DIFF%STRIDE"))
3049  CALL cp_print_key_finished_output(unit_nr, logger, input, &
3050  "DFT%QS%OPT_EMBED%EMBED_DENS_DIFF")
3051  END IF
3052 
3053  END SUBROUTINE print_rho_diff
3054 
3055 ! **************************************************************************************************
3056 !> \brief Prints a cube for the (spin_rho_A + spin_rho_B - spin_rho_ref) to be minimized in embedding
3057 !> \param spin_diff_rho_r ...
3058 !> \param i_iter ...
3059 !> \param qs_env ...
3060 !> \param final_one ...
3061 !> \author Vladimir Rybkin
3062 ! **************************************************************************************************
3063  SUBROUTINE print_rho_spin_diff(spin_diff_rho_r, i_iter, qs_env, final_one)
3064  TYPE(pw_r3d_rs_type), INTENT(IN) :: spin_diff_rho_r
3065  INTEGER, INTENT(IN) :: i_iter
3066  TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
3067  LOGICAL, INTENT(IN) :: final_one
3068 
3069  CHARACTER(LEN=default_path_length) :: filename, my_pos_cube, title
3070  INTEGER :: unit_nr
3071  TYPE(cp_logger_type), POINTER :: logger
3072  TYPE(particle_list_type), POINTER :: particles
3073  TYPE(qs_subsys_type), POINTER :: subsys
3074  TYPE(section_vals_type), POINTER :: dft_section, input
3075 
3076  NULLIFY (subsys, input)
3077 
3078  CALL get_qs_env(qs_env=qs_env, &
3079  subsys=subsys, &
3080  input=input)
3081  dft_section => section_vals_get_subs_vals(input, "DFT")
3082  CALL qs_subsys_get(subsys, particles=particles)
3083 
3084  logger => cp_get_default_logger()
3085  IF (btest(cp_print_key_should_output(logger%iter_info, input, &
3086  "DFT%QS%OPT_EMBED%EMBED_DENS_DIFF"), cp_p_file)) THEN
3087  my_pos_cube = "REWIND"
3088  IF (.NOT. final_one) THEN
3089  WRITE (filename, '(a5,I3.3,a1,I1.1)') "SPIN_DIFF_", i_iter
3090  ELSE
3091  WRITE (filename, '(a9,I3.3,a1,I1.1)') "SPIN_DIFF"
3092  END IF
3093  unit_nr = cp_print_key_unit_nr(logger, input, "DFT%QS%OPT_EMBED%EMBED_DENS_DIFF", &
3094  extension=".cube", middle_name=trim(filename), file_position=my_pos_cube, &
3095  log_filename=.false.)
3096 
3097  WRITE (title, *) "EMBEDDING SPIN DENSITY DIFFERENCE ", " optimization step ", i_iter
3098  CALL cp_pw_to_cube(spin_diff_rho_r, unit_nr, title, particles=particles, &
3099  stride=section_get_ivals(dft_section, "QS%OPT_EMBED%EMBED_DENS_DIFF%STRIDE"))
3100  CALL cp_print_key_finished_output(unit_nr, logger, input, &
3101  "DFT%QS%OPT_EMBED%EMBED_DENS_DIFF")
3102  END IF
3103 
3104  END SUBROUTINE print_rho_spin_diff
3105 ! **************************************************************************************************
3106 !> \brief Print embedding potential as a cube and as a binary (for restarting)
3107 !> \param qs_env ...
3108 !> \param dimen_aux ...
3109 !> \param embed_pot_coef ...
3110 !> \param embed_pot ...
3111 !> \param i_iter ...
3112 !> \param embed_pot_spin ...
3113 !> \param open_shell_embed ...
3114 !> \param grid_opt ...
3115 !> \param final_one ...
3116 ! **************************************************************************************************
3117  SUBROUTINE print_embed_restart(qs_env, dimen_aux, embed_pot_coef, embed_pot, i_iter, &
3118  embed_pot_spin, open_shell_embed, grid_opt, final_one)
3119  TYPE(qs_environment_type), POINTER :: qs_env
3120  INTEGER :: dimen_aux
3121  TYPE(cp_fm_type), INTENT(IN), POINTER :: embed_pot_coef
3122  TYPE(pw_r3d_rs_type), INTENT(IN) :: embed_pot
3123  INTEGER :: i_iter
3124  TYPE(pw_r3d_rs_type), INTENT(IN), POINTER :: embed_pot_spin
3125  LOGICAL :: open_shell_embed, grid_opt, final_one
3126 
3127  CHARACTER(LEN=default_path_length) :: filename, my_pos_cube, title
3128  INTEGER :: unit_nr
3129  TYPE(cp_logger_type), POINTER :: logger
3130  TYPE(particle_list_type), POINTER :: particles
3131  TYPE(qs_subsys_type), POINTER :: subsys
3132  TYPE(section_vals_type), POINTER :: dft_section, input
3133 
3134  NULLIFY (input)
3135  CALL get_qs_env(qs_env=qs_env, subsys=subsys, &
3136  input=input)
3137 
3138  ! First we print an unformatted file
3139  IF (.NOT. grid_opt) THEN ! Only for finite basis optimization
3140  logger => cp_get_default_logger()
3141  IF (btest(cp_print_key_should_output(logger%iter_info, input, &
3142  "DFT%QS%OPT_EMBED%EMBED_POT_VECTOR"), cp_p_file)) THEN
3143  IF (.NOT. final_one) THEN
3144  WRITE (filename, '(a10,I3.3)') "embed_pot_", i_iter
3145  ELSE
3146  WRITE (filename, '(a10,I3.3)') "embed_pot"
3147  END IF
3148  unit_nr = cp_print_key_unit_nr(logger, input, "DFT%QS%OPT_EMBED%EMBED_POT_VECTOR", extension=".wfn", &
3149  file_form="UNFORMATTED", middle_name=trim(filename), file_position="REWIND")
3150  IF (unit_nr > 0) THEN
3151  WRITE (unit_nr) dimen_aux
3152  END IF
3153  CALL cp_fm_write_unformatted(embed_pot_coef, unit_nr)
3154  IF (unit_nr > 0) THEN
3155  CALL close_file(unit_nr)
3156  END IF
3157  END IF
3158  END IF
3159 
3160  ! Second, cube files
3161  dft_section => section_vals_get_subs_vals(input, "DFT")
3162  CALL qs_subsys_get(subsys, particles=particles)
3163 
3164  logger => cp_get_default_logger()
3165  IF (btest(cp_print_key_should_output(logger%iter_info, input, &
3166  "DFT%QS%OPT_EMBED%EMBED_POT_CUBE"), cp_p_file)) THEN
3167  my_pos_cube = "REWIND"
3168  IF (.NOT. final_one) THEN
3169  WRITE (filename, '(a10,I3.3)') "embed_pot_", i_iter
3170  ELSE
3171  WRITE (filename, '(a10,I3.3)') "embed_pot"
3172  END IF
3173  unit_nr = cp_print_key_unit_nr(logger, input, "DFT%QS%OPT_EMBED%EMBED_POT_CUBE", &
3174  extension=".cube", middle_name=trim(filename), file_position=my_pos_cube, &
3175  log_filename=.false.)
3176 
3177  WRITE (title, *) "EMBEDDING POTENTIAL at optimization step ", i_iter
3178  CALL cp_pw_to_cube(embed_pot, unit_nr, title, particles=particles)
3179 !, &
3180 ! stride=section_get_ivals(dft_section, "QS%OPT_EMBED%EMBED_POT_CUBE%STRIDE"))
3181  CALL cp_print_key_finished_output(unit_nr, logger, input, &
3182  "DFT%QS%OPT_EMBED%EMBED_POT_CUBE")
3183  IF (open_shell_embed) THEN ! Print spin part of the embedding potential
3184  my_pos_cube = "REWIND"
3185  IF (.NOT. final_one) THEN
3186  WRITE (filename, '(a15,I3.3)') "spin_embed_pot_", i_iter
3187  ELSE
3188  WRITE (filename, '(a15,I3.3)') "spin_embed_pot"
3189  END IF
3190  unit_nr = cp_print_key_unit_nr(logger, input, "DFT%QS%OPT_EMBED%EMBED_POT_CUBE", &
3191  extension=".cube", middle_name=trim(filename), file_position=my_pos_cube, &
3192  log_filename=.false.)
3193 
3194  WRITE (title, *) "SPIN EMBEDDING POTENTIAL at optimization step ", i_iter
3195  CALL cp_pw_to_cube(embed_pot_spin, unit_nr, title, particles=particles)
3196 !, &
3197 ! stride=section_get_ivals(dft_section, "QS%OPT_EMBED%EMBED_POT_CUBE%STRIDE"))
3198  CALL cp_print_key_finished_output(unit_nr, logger, input, &
3199  "DFT%QS%OPT_EMBED%EMBED_POT_CUBE")
3200  END IF
3201  END IF
3202 
3203  END SUBROUTINE print_embed_restart
3204 
3205 ! **************************************************************************************************
3206 !> \brief Prints a volumetric file: X Y Z value for interfacing with external programs.
3207 !> \param qs_env ...
3208 !> \param embed_pot ...
3209 !> \param embed_pot_spin ...
3210 !> \param i_iter ...
3211 !> \param open_shell_embed ...
3212 !> \param final_one ...
3213 !> \param qs_env_cluster ...
3214 ! **************************************************************************************************
3215  SUBROUTINE print_pot_simple_grid(qs_env, embed_pot, embed_pot_spin, i_iter, open_shell_embed, &
3216  final_one, qs_env_cluster)
3217  TYPE(qs_environment_type), POINTER :: qs_env
3218  TYPE(pw_r3d_rs_type), INTENT(IN) :: embed_pot
3219  TYPE(pw_r3d_rs_type), INTENT(IN), POINTER :: embed_pot_spin
3220  INTEGER :: i_iter
3221  LOGICAL :: open_shell_embed, final_one
3222  TYPE(qs_environment_type), POINTER :: qs_env_cluster
3223 
3224  CHARACTER(LEN=default_path_length) :: filename
3225  INTEGER :: my_units, unit_nr
3226  LOGICAL :: angstrom, bohr
3227  TYPE(cp_logger_type), POINTER :: logger
3228  TYPE(pw_env_type), POINTER :: pw_env
3229  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
3230  TYPE(pw_r3d_rs_type) :: pot_alpha, pot_beta
3231  TYPE(section_vals_type), POINTER :: dft_section, input
3232 
3233  NULLIFY (input)
3234  CALL get_qs_env(qs_env=qs_env, input=input, pw_env=pw_env)
3235 
3236  ! Second, cube files
3237  dft_section => section_vals_get_subs_vals(input, "DFT")
3238 
3239  NULLIFY (logger)
3240  logger => cp_get_default_logger()
3241  IF (btest(cp_print_key_should_output(logger%iter_info, input, &
3242  "DFT%QS%OPT_EMBED%WRITE_SIMPLE_GRID"), cp_p_file)) THEN
3243 
3244  ! Figure out the units
3245  angstrom = .false.
3246  bohr = .true.
3247  CALL section_vals_val_get(dft_section, "QS%OPT_EMBED%WRITE_SIMPLE_GRID%UNITS", i_val=my_units)
3248  SELECT CASE (my_units)
3249  CASE (embed_grid_bohr)
3250  bohr = .true.
3251  angstrom = .false.
3252  CASE (embed_grid_angstrom)
3253  bohr = .false.
3254  angstrom = .true.
3255  CASE DEFAULT
3256  bohr = .true.
3257  angstrom = .false.
3258  END SELECT
3259 
3260  ! Get alpha and beta potentials
3261  ! Prepare plane-waves pool
3262  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
3263 
3264  ! Create embedding potential and set to zero
3265  CALL auxbas_pw_pool%create_pw(pot_alpha)
3266  CALL pw_zero(pot_alpha)
3267 
3268  CALL pw_copy(embed_pot, pot_alpha)
3269 
3270  IF (open_shell_embed) THEN
3271  CALL auxbas_pw_pool%create_pw(pot_beta)
3272  CALL pw_copy(embed_pot, pot_beta)
3273  ! Add spin potential to the alpha, and subtract from the beta part
3274  CALL pw_axpy(embed_pot_spin, pot_alpha, 1.0_dp)
3275  CALL pw_axpy(embed_pot_spin, pot_beta, -1.0_dp)
3276  END IF
3277 
3278  IF (.NOT. final_one) THEN
3279  WRITE (filename, '(a10,I3.3)') "embed_pot_", i_iter
3280  ELSE
3281  WRITE (filename, '(a10,I3.3)') "embed_pot"
3282  END IF
3283  unit_nr = cp_print_key_unit_nr(logger, input, "DFT%QS%OPT_EMBED%WRITE_SIMPLE_GRID", extension=".dat", &
3284  middle_name=trim(filename), file_form="FORMATTED", file_position="REWIND")
3285 
3286  IF (open_shell_embed) THEN ! Print spin part of the embedding potential
3287  CALL cp_pw_to_simple_volumetric(pw=pot_alpha, unit_nr=unit_nr, &
3288  stride=section_get_ivals(dft_section, "QS%OPT_EMBED%WRITE_SIMPLE_GRID%STRIDE"), &
3289  pw2=pot_beta)
3290  ELSE
3291  CALL cp_pw_to_simple_volumetric(pot_alpha, unit_nr, &
3292  stride=section_get_ivals(dft_section, "QS%OPT_EMBED%WRITE_SIMPLE_GRID%STRIDE"))
3293  END IF
3294 
3295  CALL cp_print_key_finished_output(unit_nr, logger, input, &
3296  "DFT%QS%OPT_EMBED%WRITE_SIMPLE_GRID")
3297  ! Release structures
3298  CALL pot_alpha%release()
3299  IF (open_shell_embed) THEN
3300  CALL pot_beta%release()
3301  END IF
3302 
3303  END IF
3304 
3305  ! Fold the coordinates and write into separate file: needed to have the grid correspond to coordinates
3306  ! Needed for external software.
3307  CALL print_folded_coordinates(qs_env_cluster, input)
3308 
3309  END SUBROUTINE print_pot_simple_grid
3310 
3311 ! **************************************************************************************************
3312 !> \brief ...
3313 !> \param qs_env ...
3314 !> \param input ...
3315 ! **************************************************************************************************
3316  SUBROUTINE print_folded_coordinates(qs_env, input)
3317  TYPE(qs_environment_type), POINTER :: qs_env
3318  TYPE(section_vals_type), POINTER :: input
3319 
3320  CHARACTER(LEN=2), ALLOCATABLE, DIMENSION(:) :: particles_el
3321  CHARACTER(LEN=default_path_length) :: filename
3322  INTEGER :: iat, n, unit_nr
3323  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: particles_r
3324  REAL(kind=dp), DIMENSION(3) :: center, r_pbc, s
3325  TYPE(cell_type), POINTER :: cell
3326  TYPE(cp_logger_type), POINTER :: logger
3327  TYPE(particle_list_type), POINTER :: particles
3328  TYPE(qs_subsys_type), POINTER :: subsys
3329 
3330  NULLIFY (logger)
3331  logger => cp_get_default_logger()
3332  IF (btest(cp_print_key_should_output(logger%iter_info, input, &
3333  "DFT%QS%OPT_EMBED%WRITE_SIMPLE_GRID/FOLD_COORD"), cp_p_file)) THEN
3334  CALL get_qs_env(qs_env=qs_env, cell=cell, subsys=subsys)
3335  CALL qs_subsys_get(subsys=subsys, particles=particles)
3336 
3337  ! Prepare the file
3338  WRITE (filename, '(a14)') "folded_cluster"
3339  unit_nr = cp_print_key_unit_nr(logger, input, &
3340  "DFT%QS%OPT_EMBED%WRITE_SIMPLE_GRID/FOLD_COORD", extension=".dat", &
3341  middle_name=trim(filename), file_form="FORMATTED", file_position="REWIND")
3342  IF (unit_nr > 0) THEN
3343 
3344  n = particles%n_els
3345  ALLOCATE (particles_el(n))
3346  ALLOCATE (particles_r(3, n))
3347  DO iat = 1, n
3348  CALL get_atomic_kind(particles%els(iat)%atomic_kind, element_symbol=particles_el(iat))
3349  particles_r(:, iat) = particles%els(iat)%r(:)
3350  END DO
3351 
3352  ! Fold the coordinates
3353  center(:) = cell%hmat(:, 1)/2.0_dp + cell%hmat(:, 2)/2.0_dp + cell%hmat(:, 3)/2.0_dp
3354 
3355  ! Print folded coordinates to file
3356  DO iat = 1, SIZE(particles_el)
3357  r_pbc(:) = particles_r(:, iat) - center
3358  s = matmul(cell%h_inv, r_pbc)
3359  s = s - anint(s)
3360  r_pbc = matmul(cell%hmat, s)
3361  r_pbc = r_pbc + center
3362  WRITE (unit_nr, '(a4,4f12.6)') particles_el(iat), r_pbc(:)
3363  END DO
3364 
3365  CALL cp_print_key_finished_output(unit_nr, logger, input, &
3366  "DFT%QS%OPT_EMBED%WRITE_SIMPLE_GRID/FOLD_COORD")
3367 
3368  DEALLOCATE (particles_el)
3369  DEALLOCATE (particles_r)
3370  END IF
3371 
3372  END IF ! Should output
3373 
3374  END SUBROUTINE print_folded_coordinates
3375 
3376 ! **************************************************************************************************
3377 !> \brief ...
3378 !> \param output_unit ...
3379 !> \param step_num ...
3380 !> \param opt_embed ...
3381 ! **************************************************************************************************
3382  SUBROUTINE print_emb_opt_info(output_unit, step_num, opt_embed)
3383  INTEGER :: output_unit, step_num
3384  TYPE(opt_embed_pot_type) :: opt_embed
3385 
3386  IF (output_unit > 0) THEN
3387  WRITE (unit=output_unit, fmt="(/,T2,8('-'),A,I5,1X,12('-'))") &
3388  " Optimize embedding potential info at step = ", step_num
3389  WRITE (unit=output_unit, fmt="(T2,A,F20.10)") &
3390  " Functional value = ", opt_embed%w_func(step_num)
3391  IF (step_num .GT. 1) THEN
3392  WRITE (unit=output_unit, fmt="(T2,A,F20.10)") &
3393  " Real energy change = ", opt_embed%w_func(step_num) - &
3394  opt_embed%w_func(step_num - 1)
3395 
3396  WRITE (unit=output_unit, fmt="(T2,A,F20.10)") &
3397  " Step size = ", opt_embed%step_len
3398 
3399  END IF
3400 
3401  WRITE (unit=output_unit, fmt="(T2,A,F20.10)") &
3402  " Trust radius = ", opt_embed%trust_rad
3403 
3404  WRITE (unit=output_unit, fmt="(T2,51('-'))")
3405  END IF
3406 
3407  END SUBROUTINE print_emb_opt_info
3408 
3409 ! **************************************************************************************************
3410 !> \brief ...
3411 !> \param opt_embed ...
3412 !> \param force_env ...
3413 !> \param subsys_num ...
3414 ! **************************************************************************************************
3415  SUBROUTINE get_prev_density(opt_embed, force_env, subsys_num)
3416  TYPE(opt_embed_pot_type) :: opt_embed
3417  TYPE(force_env_type), POINTER :: force_env
3418  INTEGER :: subsys_num
3419 
3420  INTEGER :: i_dens_start, i_spin, nspins
3421  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
3422  TYPE(qs_rho_type), POINTER :: rho
3423 
3424  NULLIFY (rho_r, rho)
3425  CALL get_qs_env(force_env%qs_env, rho=rho)
3426  CALL qs_rho_get(rho_struct=rho, rho_r=rho_r)
3427 
3428  nspins = opt_embed%all_nspins(subsys_num)
3429 
3430  i_dens_start = sum(opt_embed%all_nspins(1:subsys_num)) - nspins + 1
3431 
3432  DO i_spin = 1, nspins
3433  opt_embed%prev_subsys_dens(i_dens_start + i_spin - 1)%array(:, :, :) = &
3434  rho_r(i_spin)%array(:, :, :)
3435  END DO
3436 
3437  END SUBROUTINE get_prev_density
3438 
3439 ! **************************************************************************************************
3440 !> \brief ...
3441 !> \param opt_embed ...
3442 !> \param force_env ...
3443 !> \param subsys_num ...
3444 ! **************************************************************************************************
3445  SUBROUTINE get_max_subsys_diff(opt_embed, force_env, subsys_num)
3446  TYPE(opt_embed_pot_type) :: opt_embed
3447  TYPE(force_env_type), POINTER :: force_env
3448  INTEGER :: subsys_num
3449 
3450  INTEGER :: i_dens_start, i_spin, nspins
3451  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
3452  TYPE(qs_rho_type), POINTER :: rho
3453 
3454  NULLIFY (rho_r, rho)
3455  CALL get_qs_env(force_env%qs_env, rho=rho)
3456  CALL qs_rho_get(rho_struct=rho, rho_r=rho_r)
3457 
3458  nspins = opt_embed%all_nspins(subsys_num)
3459 
3460  i_dens_start = sum(opt_embed%all_nspins(1:subsys_num)) - nspins + 1
3461 
3462  DO i_spin = 1, nspins
3463  CALL pw_axpy(rho_r(i_spin), opt_embed%prev_subsys_dens(i_dens_start + i_spin - 1), 1.0_dp, -1.0_dp, &
3464  allow_noncompatible_grids=.true.)
3465  opt_embed%max_subsys_dens_diff(i_dens_start + i_spin - 1) = &
3466  max_dens_diff(opt_embed%prev_subsys_dens(i_dens_start + i_spin - 1))
3467  END DO
3468 
3469  END SUBROUTINE get_max_subsys_diff
3470 
3471 ! **************************************************************************************************
3472 !> \brief ...
3473 !> \param opt_embed ...
3474 !> \param diff_rho_r ...
3475 !> \param diff_rho_spin ...
3476 !> \param output_unit ...
3477 ! **************************************************************************************************
3478  SUBROUTINE conv_check_embed(opt_embed, diff_rho_r, diff_rho_spin, output_unit)
3479  TYPE(opt_embed_pot_type) :: opt_embed
3480  TYPE(pw_r3d_rs_type), INTENT(IN) :: diff_rho_r, diff_rho_spin
3481  INTEGER :: output_unit
3482 
3483  INTEGER :: i_dens, i_dens_start, i_spin
3484  LOGICAL :: conv_int_diff, conv_max_diff
3485  REAL(kind=dp) :: int_diff, int_diff_spin, &
3486  int_diff_square, int_diff_square_spin, &
3487  max_diff, max_diff_spin
3488 
3489  ! Calculate the convergence target values
3490  opt_embed%max_diff(1) = max_dens_diff(diff_rho_r)
3491  opt_embed%int_diff(1) = pw_integrate_function(fun=diff_rho_r, oprt='ABS')
3492  opt_embed%int_diff_square(1) = pw_integral_ab(diff_rho_r, diff_rho_r)
3493  IF (opt_embed%open_shell_embed) THEN
3494  opt_embed%max_diff(2) = max_dens_diff(diff_rho_spin)
3495  opt_embed%int_diff(2) = pw_integrate_function(fun=diff_rho_spin, oprt='ABS')
3496  opt_embed%int_diff_square(2) = pw_integral_ab(diff_rho_spin, diff_rho_spin)
3497  END IF
3498 
3499  ! Find out the convergence
3500  max_diff = opt_embed%max_diff(1)
3501 
3502  ! Maximum value criterium
3503  ! Open shell
3504  IF (opt_embed%open_shell_embed) THEN
3505  max_diff_spin = opt_embed%max_diff(2)
3506  IF ((max_diff .LE. opt_embed%conv_max) .AND. (max_diff_spin .LE. opt_embed%conv_max_spin)) THEN
3507  conv_max_diff = .true.
3508  ELSE
3509  conv_max_diff = .false.
3510  END IF
3511  ELSE
3512  ! Closed shell
3513  IF (max_diff .LE. opt_embed%conv_max) THEN
3514  conv_max_diff = .true.
3515  ELSE
3516  conv_max_diff = .false.
3517  END IF
3518  END IF
3519 
3520  ! Integrated value criterium
3521  int_diff = opt_embed%int_diff(1)
3522  ! Open shell
3523  IF (opt_embed%open_shell_embed) THEN
3524  int_diff_spin = opt_embed%int_diff(2)
3525  IF ((int_diff .LE. opt_embed%conv_int) .AND. (int_diff_spin .LE. opt_embed%conv_int_spin)) THEN
3526  conv_int_diff = .true.
3527  ELSE
3528  conv_int_diff = .false.
3529  END IF
3530  ELSE
3531  ! Closed shell
3532  IF (int_diff .LE. opt_embed%conv_int) THEN
3533  conv_int_diff = .true.
3534  ELSE
3535  conv_int_diff = .false.
3536  END IF
3537  END IF
3538 
3539  ! Integrated squared value criterium
3540  int_diff_square = opt_embed%int_diff_square(1)
3541  ! Open shell
3542  IF (opt_embed%open_shell_embed) int_diff_square_spin = opt_embed%int_diff_square(2)
3543 
3544  IF ((conv_max_diff) .AND. (conv_int_diff)) THEN
3545  opt_embed%converged = .true.
3546  ELSE
3547  opt_embed%converged = .false.
3548  END IF
3549 
3550  ! Print the information
3551  IF (output_unit > 0) THEN
3552  WRITE (unit=output_unit, fmt="(/,T2,A)") &
3553  " Convergence check :"
3554 
3555  ! Maximum value of density
3556  WRITE (unit=output_unit, fmt="(T2,A,F20.10)") &
3557  " Maximum density difference = ", max_diff
3558  WRITE (unit=output_unit, fmt="(T2,A,F20.10)") &
3559  " Convergence limit for max. density diff. = ", opt_embed%conv_max
3560 
3561  IF (opt_embed%open_shell_embed) THEN
3562 
3563  WRITE (unit=output_unit, fmt="(T2,A,F20.10)") &
3564  " Maximum spin density difference = ", max_diff_spin
3565  WRITE (unit=output_unit, fmt="(T2,A,F20.10)") &
3566  " Convergence limit for max. spin dens.diff.= ", opt_embed%conv_max_spin
3567 
3568  END IF
3569 
3570  IF (conv_max_diff) THEN
3571  WRITE (unit=output_unit, fmt="(T2,2A)") &
3572  " Convergence in max. density diff. = ", &
3573  " YES"
3574  ELSE
3575  WRITE (unit=output_unit, fmt="(T2,2A)") &
3576  " Convergence in max. density diff. = ", &
3577  " NO"
3578  END IF
3579 
3580  ! Integrated abs. value of density
3581  WRITE (unit=output_unit, fmt="(T2,A,F20.10)") &
3582  " Integrated density difference = ", int_diff
3583  WRITE (unit=output_unit, fmt="(T2,A,F20.10)") &
3584  " Conv. limit for integrated density diff. = ", opt_embed%conv_int
3585  IF (opt_embed%open_shell_embed) THEN
3586  WRITE (unit=output_unit, fmt="(T2,A,F20.10)") &
3587  " Integrated spin density difference = ", int_diff_spin
3588  WRITE (unit=output_unit, fmt="(T2,A,F20.10)") &
3589  " Conv. limit for integrated spin dens.diff.= ", opt_embed%conv_int_spin
3590  END IF
3591 
3592  IF (conv_int_diff) THEN
3593  WRITE (unit=output_unit, fmt="(T2,2A)") &
3594  " Convergence in integrated density diff. = ", &
3595  " YES"
3596  ELSE
3597  WRITE (unit=output_unit, fmt="(T2,2A)") &
3598  " Convergence in integrated density diff. = ", &
3599  " NO"
3600  END IF
3601 
3602  ! Integrated squared value of density
3603  WRITE (unit=output_unit, fmt="(T2,A,F20.10)") &
3604  " Integrated squared density difference = ", int_diff_square
3605  IF (opt_embed%open_shell_embed) THEN
3606  WRITE (unit=output_unit, fmt="(T2,A,F20.10)") &
3607  " Integrated squared spin density difference= ", int_diff_square_spin
3608  END IF
3609 
3610  ! Maximum subsystem density change
3611  WRITE (unit=output_unit, fmt="(/,T2,A)") &
3612  " Maximum density change in:"
3613  DO i_dens = 1, (SIZE(opt_embed%all_nspins) - 1)
3614  i_dens_start = sum(opt_embed%all_nspins(1:i_dens)) - opt_embed%all_nspins(i_dens) + 1
3615  DO i_spin = 1, opt_embed%all_nspins(i_dens)
3616  WRITE (unit=output_unit, fmt="(T4,A10,I3,A6,I3,A1,F20.10)") &
3617  " subsystem ", i_dens, ', spin', i_spin, ":", &
3618  opt_embed%max_subsys_dens_diff(i_dens_start + i_spin - 1)
3619  END DO
3620  END DO
3621 
3622  END IF
3623 
3624  IF ((opt_embed%converged) .AND. (output_unit > 0)) THEN
3625  WRITE (unit=output_unit, fmt="(/,T2,A)") repeat("*", 79)
3626  WRITE (unit=output_unit, fmt="(T2,A,T25,A,T78,A)") &
3627  "***", "EMBEDDING POTENTIAL OPTIMIZATION COMPLETED", "***"
3628  WRITE (unit=output_unit, fmt="(T2,A)") repeat("*", 79)
3629  END IF
3630 
3631  END SUBROUTINE conv_check_embed
3632 
3633 END MODULE optimize_embedding_potential
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition: cp_dbcsr_operations.F:89
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition: atomic_kind_types.F:223
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition: atomic_kind_types.F:138
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cp_blacs_env
methods related to the blacs parallel environment
Definition: cp_blacs_env.F:15
cp_blacs_env::cp_blacs_env_release
subroutine, public cp_blacs_env_release(blacs_env)
releases the given blacs_env
Definition: cp_blacs_env.F:282
cp_blacs_env::cp_blacs_env_create
subroutine, public cp_blacs_env_create(blacs_env, para_env, blacs_grid_layout, blacs_repeatable, row_major, grid_2d)
allocates and initializes a type that represent a blacs context
Definition: cp_blacs_env.F:123
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition: cp_dbcsr_operations.F:208
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition: cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition: cp_files.F:308
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition: cp_files.F:119
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition: cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_column_scale
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
Definition: cp_fm_basic_linalg.F:1816
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition: cp_fm_basic_linalg.F:101
cp_fm_basic_linalg::cp_fm_scale
subroutine, public cp_fm_scale(alpha, matrix_a)
scales a matrix matrix_a = alpha * matrix_b
Definition: cp_fm_basic_linalg.F:1614
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition: cp_fm_diag.F:17
cp_fm_diag::choose_eigv_solver
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition: cp_fm_diag.F:208
cp_fm_struct
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition: cp_fm_struct.F:125
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
cp_fm_types
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
cp_fm_types::cp_fm_copy_general
subroutine, public cp_fm_copy_general(source, destination, para_env)
General copy of a fm matrix to another fm matrix. Uses non-blocking MPI rather than ScaLAPACK.
Definition: cp_fm_types.F:1538
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition: cp_fm_types.F:1016
cp_fm_types::cp_fm_write_unformatted
subroutine, public cp_fm_write_unformatted(fm, unit)
...
Definition: cp_fm_types.F:2147
cp_fm_types::cp_fm_get_element
subroutine, public cp_fm_get_element(matrix, irow_global, icol_global, alpha, local)
returns an element of a fm this value is valid on every cpu using this call is expensive
Definition: cp_fm_types.F:643
cp_fm_types::cp_fm_to_fm_submat
subroutine, public cp_fm_to_fm_submat(msource, mtarget, nrow, ncol, s_firstrow, s_firstcol, t_firstrow, t_firstcol)
copy just a part ot the matrix
Definition: cp_fm_types.F:1473
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition: cp_fm_types.F:535
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition: cp_output_handling.F:1013
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition: cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition: cp_output_handling.F:345
cp_realspace_grid_cube
A wrapper around pw_to_cube() which accepts particle_list_type.
Definition: cp_realspace_grid_cube.F:12
cp_realspace_grid_cube::cp_pw_to_cube
subroutine, public cp_pw_to_cube(pw, unit_nr, title, particles, stride, zero_tails, silent, mpi_io)
...
Definition: cp_realspace_grid_cube.F:45
cp_realspace_grid_cube::cp_cube_to_pw
subroutine, public cp_cube_to_pw(grid, filename, scaling, silent)
Thin wrapper around routine cube_to_pw.
Definition: cp_realspace_grid_cube.F:113
cp_realspace_grid_cube::cp_pw_to_simple_volumetric
subroutine, public cp_pw_to_simple_volumetric(pw, unit_nr, stride, pw2)
Prints grid in a simple format: X Y Z value.
Definition: cp_realspace_grid_cube.F:90
embed_types
Definition: embed_types.F:8
force_env_types
Interface for the force calculations.
Definition: force_env_types.F:18
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::embed_grid_angstrom
integer, parameter, public embed_grid_angstrom
Definition: input_constants.F:1253
input_constants::embed_steep_desc
integer, parameter, public embed_steep_desc
Definition: input_constants.F:1250
input_constants::embed_level_shift
integer, parameter, public embed_level_shift
Definition: input_constants.F:1250
input_constants::embed_resp
integer, parameter, public embed_resp
Definition: input_constants.F:1255
input_constants::embed_quasi_newton
integer, parameter, public embed_quasi_newton
Definition: input_constants.F:1250
input_constants::embed_none
integer, parameter, public embed_none
Definition: input_constants.F:1255
input_constants::embed_fa
integer, parameter, public embed_fa
Definition: input_constants.F:1255
input_constants::embed_diff
integer, parameter, public embed_diff
Definition: input_constants.F:1255
input_constants::embed_grid_bohr
integer, parameter, public embed_grid_bohr
Definition: input_constants.F:1253
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_get_rval
real(kind=dp) function, public section_get_rval(section_vals, keyword_name)
...
Definition: input_section_types.F:941
input_section_types::section_get_ivals
integer function, dimension(:), pointer, public section_get_ivals(section_vals, keyword_name)
...
Definition: input_section_types.F:989
input_section_types::section_get_ival
integer function, public section_get_ival(section_vals, keyword_name)
...
Definition: input_section_types.F:973
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition: input_section_types.F:724
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition: input_section_types.F:1040
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_path_length
integer, parameter, public default_path_length
Definition: kinds.F:58
lri_environment_types
contains the types and subroutines for dealing with the lri_env lri : local resolution of the identit...
Definition: lri_environment_types.F:18
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition: mathconstants.F:110
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
mixed_environment_utils
Util mixed_environment.
Definition: mixed_environment_utils.F:12
mixed_environment_utils::get_subsys_map_index
subroutine, public get_subsys_map_index(mapping_section, natom, iforce_eval, nforce_eval, map_index, force_eval_embed)
performs mapping of the subsystems of different force_eval
Definition: mixed_environment_utils.F:141
optimize_embedding_potential
Definition: optimize_embedding_potential.F:8
optimize_embedding_potential::print_embed_restart
subroutine, public print_embed_restart(qs_env, dimen_aux, embed_pot_coef, embed_pot, i_iter, embed_pot_spin, open_shell_embed, grid_opt, final_one)
Print embedding potential as a cube and as a binary (for restarting)
Definition: optimize_embedding_potential.F:3119
optimize_embedding_potential::make_subsys_embed_pot
subroutine, public make_subsys_embed_pot(qs_env, embed_pot, embed_pot_subsys, spin_embed_pot, spin_embed_pot_subsys, open_shell_embed, change_spin_sign)
Creates a subsystem embedding potential.
Definition: optimize_embedding_potential.F:1024
optimize_embedding_potential::print_rho_spin_diff
subroutine, public print_rho_spin_diff(spin_diff_rho_r, i_iter, qs_env, final_one)
Prints a cube for the (spin_rho_A + spin_rho_B - spin_rho_ref) to be minimized in embedding.
Definition: optimize_embedding_potential.F:3064
optimize_embedding_potential::get_max_subsys_diff
subroutine, public get_max_subsys_diff(opt_embed, force_env, subsys_num)
...
Definition: optimize_embedding_potential.F:3446
optimize_embedding_potential::print_emb_opt_info
subroutine, public print_emb_opt_info(output_unit, step_num, opt_embed)
...
Definition: optimize_embedding_potential.F:3383
optimize_embedding_potential::understand_spin_states
subroutine, public understand_spin_states(force_env, ref_subsys_number, change_spin, open_shell_embed, all_nspins)
Find out whether we need to swap alpha- and beta- spind densities in the second subsystem.
Definition: optimize_embedding_potential.F:128
optimize_embedding_potential::read_embed_pot
subroutine, public read_embed_pot(qs_env, embed_pot, spin_embed_pot, section, opt_embed)
...
Definition: optimize_embedding_potential.F:663
optimize_embedding_potential::given_embed_pot
subroutine, public given_embed_pot(qs_env)
Read the external embedding potential, not to be optimized.
Definition: optimize_embedding_potential.F:610
optimize_embedding_potential::find_aux_dimen
subroutine, public find_aux_dimen(qs_env, dimen_aux)
Find the dimension of the auxiliary basis for the expansion of the embedding potential.
Definition: optimize_embedding_potential.F:544
optimize_embedding_potential::get_prev_density
subroutine, public get_prev_density(opt_embed, force_env, subsys_num)
...
Definition: optimize_embedding_potential.F:3416
optimize_embedding_potential::coulomb_guess
subroutine, public coulomb_guess(v_rspace, rhs, mapping_section, qs_env, nforce_eval, iforce_eval, eta)
Calculates subsystem Coulomb potential from the RESP charges of the total system.
Definition: optimize_embedding_potential.F:933
optimize_embedding_potential::print_rho_diff
subroutine, public print_rho_diff(diff_rho_r, i_iter, qs_env, final_one)
Prints a cube for the (rho_A + rho_B - rho_ref) to be minimized in embedding.
Definition: optimize_embedding_potential.F:3013
optimize_embedding_potential::conv_check_embed
subroutine, public conv_check_embed(opt_embed, diff_rho_r, diff_rho_spin, output_unit)
...
Definition: optimize_embedding_potential.F:3479
optimize_embedding_potential::max_dens_diff
real(kind=dp) function, public max_dens_diff(diff_rho_r)
...
Definition: optimize_embedding_potential.F:2974
optimize_embedding_potential::step_control
subroutine, public step_control(opt_embed)
Controls the step, changes the trust radius if needed in maximization of the V_emb.
Definition: optimize_embedding_potential.F:2170
optimize_embedding_potential::opt_embed_step
subroutine, public opt_embed_step(diff_rho_r, diff_rho_spin, opt_embed, embed_pot, spin_embed_pot, rho_r_ref, qs_env)
Takes maximization step in embedding potential optimization.
Definition: optimize_embedding_potential.F:1467
optimize_embedding_potential::calculate_embed_pot_grad
subroutine, public calculate_embed_pot_grad(qs_env, diff_rho_r, diff_rho_spin, opt_embed)
Calculates the derivative of the embedding potential wrt to the expansion coefficients.
Definition: optimize_embedding_potential.F:1072
optimize_embedding_potential::print_pot_simple_grid
subroutine, public print_pot_simple_grid(qs_env, embed_pot, embed_pot_spin, i_iter, open_shell_embed, final_one, qs_env_cluster)
Prints a volumetric file: X Y Z value for interfacing with external programs.
Definition: optimize_embedding_potential.F:3217
optimize_embedding_potential::prepare_embed_opt
subroutine, public prepare_embed_opt(qs_env, opt_embed, opt_embed_section)
Creates and allocates objects for optimization of embedding potential.
Definition: optimize_embedding_potential.F:291
optimize_embedding_potential::release_opt_embed
subroutine, public release_opt_embed(opt_embed)
Deallocate stuff for optimizing embedding potential.
Definition: optimize_embedding_potential.F:856
optimize_embedding_potential::init_embed_pot
subroutine, public init_embed_pot(qs_env, embed_pot, add_const_pot, Fermi_Amaldi, const_pot, open_shell_embed, spin_embed_pot, pot_diff, Coulomb_guess, grid_opt)
...
Definition: optimize_embedding_potential.F:206
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition: parallel_gemm_api.F:14
particle_list_types
represent a simple array based list of the given type
Definition: particle_list_types.F:15
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
pw_env_types
container for various plainwaves related things
Definition: pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition: pw_env_types.F:113
pw_methods
Definition: pw_methods.F:25
pw_methods::pw_derive
subroutine, public pw_derive(pw, n)
Calculate the derivative of a plane wave vector.
Definition: pw_methods.F:10106
pw_methods::pw_dr2
subroutine, public pw_dr2(pw, pwdr2, i, j)
Calculate the tensorial 2nd derivative of a plane wave vector.
Definition: pw_methods.F:10204
pw_poisson_methods
Definition: pw_poisson_methods.F:13
pw_poisson_types
functions related to the poisson solver on regular grids
Definition: pw_poisson_types.F:22
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
pw_types
Definition: pw_types.F:24
qs_collocate_density
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
Definition: qs_collocate_density.F:38
qs_collocate_density::calculate_wavefunction
subroutine, public calculate_wavefunction(mo_vectors, ivector, rho, rho_gspace, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, pw_env, basis_type, external_vector)
maps a given wavefunction on the grid
Definition: qs_collocate_density.F:2860
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_environment_types::set_qs_env
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, WannierCentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Set the QUICKSTEP environment.
Definition: qs_environment_types.F:1035
qs_integrate_potential_single
Build up the plane wave density by collocating the primitive Gaussian functions (pgf).
Definition: qs_integrate_potential_single.F:30
qs_integrate_potential_single::integrate_v_rspace_one_center
subroutine, public integrate_v_rspace_one_center(v_rspace, qs_env, int_res, calculate_forces, basis_type, atomlist)
computes integrals of product of v_rspace times a one-center function required for LRIGPW
Definition: qs_integrate_potential_single.F:763
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451
qs_kinetic
Calculation of kinetic energy matrix and forces.
Definition: qs_kinetic.F:15
qs_kinetic::build_kinetic_matrix
subroutine, public build_kinetic_matrix(ks_env, matrix_t, matrixkp_t, matrix_name, basis_type, sab_nl, calculate_forces, matrix_p, matrixkp_p, eps_filter, nderivative)
Calculation of the kinetic energy matrix over Cartesian Gaussian functions.
Definition: qs_kinetic.F:101
qs_ks_types
Definition: qs_ks_types.F:15
qs_mo_types
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kTS, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition: qs_mo_types.F:397
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition: qs_neighbor_list_types.F:17
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229
qs_subsys_types
types that represent a quickstep subsys
Definition: qs_subsys_types.F:12
qs_subsys_types::qs_subsys_get
subroutine, public qs_subsys_get(subsys, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, cp_subsys, nelectron_total, nelectron_spin)
...
Definition: qs_subsys_types.F:134
xc_rho_cflags_types
contains the structure
Definition: xc_rho_cflags_types.F:14
xc_rho_cflags_types::xc_rho_cflags_setall
elemental subroutine, public xc_rho_cflags_setall(cflags, value)
sets all the flags to the given value
Definition: xc_rho_cflags_types.F:63
xc_rho_set_types
contains the structure
Definition: xc_rho_set_types.F:14
xc_rho_set_types::xc_rho_set_update
subroutine, public xc_rho_set_update(rho_set, rho_r, rho_g, tau, needs, xc_deriv_method_id, xc_rho_smooth_id, pw_pool)
updates the given rho set with the density given by rho_r (and rho_g). The rho set will contain the c...
Definition: xc_rho_set_types.F:691
xc_rho_set_types::xc_rho_set_create
subroutine, public xc_rho_set_create(rho_set, local_bounds, rho_cutoff, drho_cutoff, tau_cutoff)
allocates and does (minimal) initialization of a rho_set
Definition: xc_rho_set_types.F:111
xc_rho_set_types::xc_rho_set_release
subroutine, public xc_rho_set_release(rho_set, pw_pool)
releases the given rho_set
Definition: xc_rho_set_types.F:129
xc
Exchange and Correlation functional calculations.
Definition: xc.F:17
xc::smooth_cutoff
subroutine, public smooth_cutoff(pot, rho, rhoa, rhob, rho_cutoff, rho_smooth_cutoff_range, e_0, e_0_scale_factor)
smooths the cutoff on rho with a function smoothderiv_rho that is 0 for rho<rho_cutoff and 1 for rho>...
Definition: xc.F:910