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