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