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