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