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