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