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