1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Density Derived atomic point charges from a QM calculation 8!> (see Bloechl, J. Chem. Phys. Vol. 103 pp. 7422-7428) 9!> \par History 10!> 08.2005 created [tlaino] 11!> \author Teodoro Laino 12! ************************************************************************************************** 13MODULE cp_ddapc_util 14 15 USE atomic_charges, ONLY: print_atomic_charges 16 USE cell_types, ONLY: cell_type 17 USE cp_control_types, ONLY: ddapc_restraint_type,& 18 dft_control_type 19 USE cp_ddapc_forces, ONLY: evaluate_restraint_functional 20 USE cp_ddapc_methods, ONLY: build_A_matrix,& 21 build_b_vector,& 22 build_der_A_matrix_rows,& 23 build_der_b_vector,& 24 cleanup_g_dot_rvec_sin_cos,& 25 prep_g_dot_rvec_sin_cos 26 USE cp_ddapc_types, ONLY: cp_ddapc_create,& 27 cp_ddapc_type 28 USE cp_log_handling, ONLY: cp_get_default_logger,& 29 cp_logger_type 30 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 31 cp_print_key_unit_nr 32 USE cp_para_types, ONLY: cp_para_env_type 33 USE input_constants, ONLY: do_full_density,& 34 do_spin_density 35 USE input_section_types, ONLY: section_vals_get_subs_vals,& 36 section_vals_type,& 37 section_vals_val_get 38 USE kinds, ONLY: default_string_length,& 39 dp 40 USE mathconstants, ONLY: pi 41 USE message_passing, ONLY: mp_sum 42 USE particle_types, ONLY: particle_type 43 USE pw_env_types, ONLY: pw_env_get,& 44 pw_env_type 45 USE pw_methods, ONLY: pw_axpy,& 46 pw_copy,& 47 pw_transfer 48 USE pw_pool_types, ONLY: pw_pool_create_pw,& 49 pw_pool_give_back_pw,& 50 pw_pool_type 51 USE pw_types, ONLY: COMPLEXDATA1D,& 52 RECIPROCALSPACE,& 53 pw_p_type,& 54 pw_type 55 USE qs_charges_types, ONLY: qs_charges_type 56 USE qs_environment_types, ONLY: get_qs_env,& 57 qs_environment_type 58 USE qs_kind_types, ONLY: qs_kind_type 59 USE qs_rho_types, ONLY: qs_rho_get,& 60 qs_rho_type 61#include "./base/base_uses.f90" 62 63 IMPLICIT NONE 64 PRIVATE 65 66 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE. 67 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_ddapc_util' 68 PUBLIC :: get_ddapc, & 69 restraint_functional_potential, & 70 modify_hartree_pot, & 71 cp_ddapc_init 72 73CONTAINS 74 75! ************************************************************************************************** 76!> \brief Initialize the cp_ddapc_environment 77!> \param qs_env ... 78!> \par History 79!> 08.2005 created [tlaino] 80!> \author Teodoro Laino 81! ************************************************************************************************** 82 SUBROUTINE cp_ddapc_init(qs_env) 83 TYPE(qs_environment_type), POINTER :: qs_env 84 85 CHARACTER(len=*), PARAMETER :: routineN = 'cp_ddapc_init', routineP = moduleN//':'//routineN 86 87 INTEGER :: handle, i, iw, iw2, n_rep_val, num_gauss 88 LOGICAL :: allocate_ddapc_env, unimplemented 89 REAL(KIND=dp) :: gcut, pfact, rcmin, Vol 90 REAL(KIND=dp), DIMENSION(:), POINTER :: inp_radii, radii 91 TYPE(cell_type), POINTER :: cell, super_cell 92 TYPE(cp_logger_type), POINTER :: logger 93 TYPE(cp_para_env_type), POINTER :: para_env 94 TYPE(dft_control_type), POINTER :: dft_control 95 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 96 TYPE(pw_env_type), POINTER :: pw_env 97 TYPE(pw_p_type), POINTER :: rho0_s_gs, rho_core 98 TYPE(pw_pool_type), POINTER :: auxbas_pool 99 TYPE(pw_type), POINTER :: rho_tot_g 100 TYPE(qs_charges_type), POINTER :: qs_charges 101 TYPE(qs_rho_type), POINTER :: rho 102 TYPE(section_vals_type), POINTER :: density_fit_section 103 104 CALL timeset(routineN, handle) 105 logger => cp_get_default_logger() 106 NULLIFY (dft_control, rho, rho_tot_g, rho_core, rho0_s_gs, pw_env, & 107 radii, inp_radii, particle_set, qs_charges, para_env) 108 109 CALL get_qs_env(qs_env, dft_control=dft_control) 110 allocate_ddapc_env = qs_env%cp_ddapc_ewald%do_solvation .OR. & 111 qs_env%cp_ddapc_ewald%do_qmmm_periodic_decpl .OR. & 112 qs_env%cp_ddapc_ewald%do_decoupling .OR. & 113 qs_env%cp_ddapc_ewald%do_restraint 114 unimplemented = dft_control%qs_control%semi_empirical .OR. & 115 dft_control%qs_control%dftb .OR. & 116 dft_control%qs_control%xtb 117 IF (allocate_ddapc_env .AND. unimplemented) THEN 118 CPABORT("DDAP charges work only with GPW/GAPW code.") 119 END IF 120 allocate_ddapc_env = allocate_ddapc_env .OR. & 121 qs_env%cp_ddapc_ewald%do_property 122 allocate_ddapc_env = allocate_ddapc_env .AND. (.NOT. unimplemented) 123 IF (allocate_ddapc_env) THEN 124 CALL get_qs_env(qs_env=qs_env, & 125 dft_control=dft_control, & 126 rho=rho, & 127 rho_core=rho_core, & 128 rho0_s_gs=rho0_s_gs, & 129 pw_env=pw_env, & 130 qs_charges=qs_charges, & 131 particle_set=particle_set, & 132 cell=cell, & 133 super_cell=super_cell, & 134 para_env=para_env) 135 density_fit_section => section_vals_get_subs_vals(qs_env%input, "DFT%DENSITY_FITTING") 136 iw = cp_print_key_unit_nr(logger, density_fit_section, & 137 "PROGRAM_RUN_INFO", ".FitCharge") 138 IF (iw > 0) THEN 139 WRITE (iw, '(/,A)') " Initializing the DDAPC Environment" 140 END IF 141 CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pool) 142 CALL pw_pool_create_pw(auxbas_pool, rho_tot_g, in_space=RECIPROCALSPACE, & 143 use_data=COMPLEXDATA1D) 144 Vol = rho_tot_g%pw_grid%vol 145 ! 146 ! Get Input Parameters 147 ! 148 CALL section_vals_val_get(density_fit_section, "RADII", n_rep_val=n_rep_val) 149 IF (n_rep_val /= 0) THEN 150 CALL section_vals_val_get(density_fit_section, "RADII", r_vals=inp_radii) 151 num_gauss = SIZE(inp_radii) 152 ALLOCATE (radii(num_gauss)) 153 radii = inp_radii 154 ELSE 155 CALL section_vals_val_get(density_fit_section, "NUM_GAUSS", i_val=num_gauss) 156 CALL section_vals_val_get(density_fit_section, "MIN_RADIUS", r_val=rcmin) 157 CALL section_vals_val_get(density_fit_section, "PFACTOR", r_val=pfact) 158 ALLOCATE (radii(num_gauss)) 159 DO i = 1, num_gauss 160 radii(i) = rcmin*pfact**(i - 1) 161 END DO 162 END IF 163 CALL section_vals_val_get(density_fit_section, "GCUT", r_val=gcut) 164 ! Create DDAPC environment 165 iw2 = cp_print_key_unit_nr(logger, density_fit_section, & 166 "PROGRAM_RUN_INFO/CONDITION_NUMBER", ".FitCharge") 167 ! Initialization of the cp_ddapc_env and of the cp_ddapc_ewald environment 168 !NB pass qs_env%para_env for parallelization of ewald_ddapc_pot() 169 CALL cp_ddapc_create(para_env, & 170 qs_env%cp_ddapc_env, & 171 qs_env%cp_ddapc_ewald, & 172 particle_set, & 173 radii, & 174 cell, & 175 super_cell, & 176 rho_tot_g, & 177 gcut, & 178 iw2, & 179 Vol, & 180 qs_env%input) 181 CALL cp_print_key_finished_output(iw2, logger, density_fit_section, & 182 "PROGRAM_RUN_INFO/CONDITION_NUMBER") 183 DEALLOCATE (radii) 184 CALL pw_pool_give_back_pw(auxbas_pool, rho_tot_g, & 185 accept_non_compatible=.TRUE.) 186 END IF 187 CALL timestop(handle) 188 END SUBROUTINE cp_ddapc_init 189 190! ************************************************************************************************** 191!> \brief Computes the Density Derived Atomic Point Charges 192!> \param qs_env ... 193!> \param calc_force ... 194!> \param density_fit_section ... 195!> \param density_type ... 196!> \param qout1 ... 197!> \param qout2 ... 198!> \param out_radii ... 199!> \param dq_out ... 200!> \param ext_rho_tot_g ... 201!> \param Itype_of_density ... 202!> \param iwc ... 203!> \par History 204!> 08.2005 created [tlaino] 205!> \author Teodoro Laino 206! ************************************************************************************************** 207 RECURSIVE SUBROUTINE get_ddapc(qs_env, calc_force, density_fit_section, & 208 density_type, qout1, qout2, out_radii, dq_out, ext_rho_tot_g, & 209 Itype_of_density, iwc) 210 TYPE(qs_environment_type), POINTER :: qs_env 211 LOGICAL, INTENT(in), OPTIONAL :: calc_force 212 TYPE(section_vals_type), POINTER :: density_fit_section 213 INTEGER, OPTIONAL :: density_type 214 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: qout1, qout2, out_radii 215 REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, & 216 POINTER :: dq_out 217 TYPE(pw_type), OPTIONAL, POINTER :: ext_rho_tot_g 218 CHARACTER(LEN=*), OPTIONAL :: Itype_of_density 219 INTEGER, INTENT(IN), OPTIONAL :: iwc 220 221 CHARACTER(len=*), PARAMETER :: routineN = 'get_ddapc', routineP = moduleN//':'//routineN 222 223 CHARACTER(LEN=default_string_length) :: type_of_density 224 INTEGER :: handle, handle2, handle3, i, ii, & 225 iparticle, iparticle0, ispin, iw, j, & 226 myid, n_rep_val, ndim, nparticles, & 227 num_gauss, pmax, pmin 228 LOGICAL :: need_f 229 REAL(KIND=dp) :: c1, c3, ch_dens, gcut, pfact, rcmin, Vol 230 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: AmI_bv, AmI_cv, bv, cv, cvT_AmI, & 231 cvT_AmI_dAmj, dAmj_qv, qtot, qv 232 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dbv, g_dot_rvec_cos, g_dot_rvec_sin 233 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dAm, dqv, tv 234 REAL(KIND=dp), DIMENSION(:), POINTER :: inp_radii, radii 235 TYPE(cell_type), POINTER :: cell, super_cell 236 TYPE(cp_ddapc_type), POINTER :: cp_ddapc_env 237 TYPE(cp_logger_type), POINTER :: logger 238 TYPE(dft_control_type), POINTER :: dft_control 239 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 240 TYPE(pw_env_type), POINTER :: pw_env 241 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_g, rho_r 242 TYPE(pw_p_type), POINTER :: rho0_s_gs, rho_core 243 TYPE(pw_pool_type), POINTER :: auxbas_pool 244 TYPE(pw_type), POINTER :: rho_tot_g 245 TYPE(qs_charges_type), POINTER :: qs_charges 246 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 247 TYPE(qs_rho_type), POINTER :: rho 248 249!NB variables for doing build_der_A_matrix_rows in blocks 250!NB refactor math in inner loop - no need for dqv0 251!!NB refactor math in inner loop - new temporaries 252 253 EXTERNAL dgemv, dgemm 254 255 CALL timeset(routineN, handle) 256 need_f = .FALSE. 257 IF (PRESENT(calc_force)) need_f = calc_force 258 logger => cp_get_default_logger() 259 NULLIFY (dft_control, rho, rho_tot_g, rho_core, rho0_s_gs, pw_env, rho_g, rho_r, & 260 radii, inp_radii, particle_set, qs_kind_set, qs_charges, cp_ddapc_env) 261 CALL get_qs_env(qs_env=qs_env, & 262 dft_control=dft_control, & 263 rho=rho, & 264 rho_core=rho_core, & 265 rho0_s_gs=rho0_s_gs, & 266 pw_env=pw_env, & 267 qs_charges=qs_charges, & 268 particle_set=particle_set, & 269 qs_kind_set=qs_kind_set, & 270 cell=cell, & 271 super_cell=super_cell) 272 273 CALL qs_rho_get(rho, rho_r=rho_r, rho_g=rho_g) 274 275 IF (PRESENT(iwc)) THEN 276 iw = iwc 277 ELSE 278 iw = cp_print_key_unit_nr(logger, density_fit_section, & 279 "PROGRAM_RUN_INFO", ".FitCharge") 280 END IF 281 CALL pw_env_get(pw_env=pw_env, & 282 auxbas_pw_pool=auxbas_pool) 283 CALL pw_pool_create_pw(auxbas_pool, rho_tot_g, in_space=RECIPROCALSPACE, & 284 use_data=COMPLEXDATA1D) 285 IF (PRESENT(ext_rho_tot_g)) THEN 286 ! If provided use the input density in g-space 287 CALL pw_transfer(ext_rho_tot_g, rho_tot_g) 288 type_of_density = Itype_of_density 289 ELSE 290 IF (PRESENT(density_type)) THEN 291 myid = density_type 292 ELSE 293 CALL section_vals_val_get(qs_env%input, & 294 "PROPERTIES%FIT_CHARGE%TYPE_OF_DENSITY", i_val=myid) 295 END IF 296 SELECT CASE (myid) 297 CASE (do_full_density) 298 ! Otherwise build the total QS density (electron+nuclei) in G-space 299 IF (dft_control%qs_control%gapw) THEN 300 CALL pw_transfer(rho0_s_gs%pw, rho_tot_g) 301 ELSE 302 CALL pw_transfer(rho_core%pw, rho_tot_g) 303 END IF 304 DO ispin = 1, SIZE(rho_g) 305 CALL pw_axpy(rho_g(ispin)%pw, rho_tot_g) 306 END DO 307 type_of_density = "FULL DENSITY" 308 CASE (do_spin_density) 309 CALL pw_copy(rho_g(1)%pw, rho_tot_g) 310 CALL pw_axpy(rho_g(2)%pw, rho_tot_g, alpha=-1._dp) 311 type_of_density = "SPIN DENSITY" 312 END SELECT 313 END IF 314 Vol = rho_r(1)%pw%pw_grid%vol 315 ch_dens = 0.0_dp 316 ! should use pw_integrate 317 IF (rho_tot_g%pw_grid%have_g0) ch_dens = REAL(rho_tot_g%cc(1), KIND=dp) 318 CALL mp_sum(ch_dens, logger%para_env%group) 319 ! 320 ! Get Input Parameters 321 ! 322 CALL section_vals_val_get(density_fit_section, "RADII", n_rep_val=n_rep_val) 323 IF (n_rep_val /= 0) THEN 324 CALL section_vals_val_get(density_fit_section, "RADII", r_vals=inp_radii) 325 num_gauss = SIZE(inp_radii) 326 ALLOCATE (radii(num_gauss)) 327 radii = inp_radii 328 ELSE 329 CALL section_vals_val_get(density_fit_section, "NUM_GAUSS", i_val=num_gauss) 330 CALL section_vals_val_get(density_fit_section, "MIN_RADIUS", r_val=rcmin) 331 CALL section_vals_val_get(density_fit_section, "PFACTOR", r_val=pfact) 332 ALLOCATE (radii(num_gauss)) 333 DO i = 1, num_gauss 334 radii(i) = rcmin*pfact**(i - 1) 335 END DO 336 END IF 337 IF (PRESENT(out_radii)) THEN 338 IF (ASSOCIATED(out_radii)) THEN 339 DEALLOCATE (out_radii) 340 END IF 341 ALLOCATE (out_radii(SIZE(radii))) 342 out_radii = radii 343 END IF 344 CALL section_vals_val_get(density_fit_section, "GCUT", r_val=gcut) 345 cp_ddapc_env => qs_env%cp_ddapc_env 346 ! 347 ! Start with the linear system 348 ! 349 ndim = SIZE(particle_set)*SIZE(radii) 350 ALLOCATE (bv(ndim)) 351 ALLOCATE (qv(ndim)) 352 ALLOCATE (qtot(SIZE(particle_set))) 353 ALLOCATE (cv(ndim)) 354 CALL timeset(routineN//"-charges", handle2) 355 bv(:) = 0.0_dp 356 cv(:) = 1.0_dp/Vol 357 CALL build_b_vector(bv, cp_ddapc_env%gfunc, cp_ddapc_env%w, & 358 particle_set, radii, rho_tot_g, gcut); bv(:) = bv(:)/Vol 359 CALL mp_sum(bv, rho_tot_g%pw_grid%para%group) 360 c1 = DOT_PRODUCT(cv, MATMUL(cp_ddapc_env%AmI, bv)) - ch_dens 361 c1 = c1/cp_ddapc_env%c0 362 qv(:) = -MATMUL(cp_ddapc_env%AmI, (bv - c1*cv)) 363 j = 0 364 qtot = 0.0_dp 365 DO i = 1, ndim, num_gauss 366 j = j + 1 367 DO ii = 1, num_gauss 368 qtot(j) = qtot(j) + qv((i - 1) + ii) 369 END DO 370 END DO 371 IF (PRESENT(qout1)) THEN 372 IF (ASSOCIATED(qout1)) THEN 373 CPASSERT(SIZE(qout1) == SIZE(qv)) 374 ELSE 375 ALLOCATE (qout1(SIZE(qv))) 376 END IF 377 qout1 = qv 378 END IF 379 IF (PRESENT(qout2)) THEN 380 IF (ASSOCIATED(qout2)) THEN 381 CPASSERT(SIZE(qout2) == SIZE(qtot)) 382 ELSE 383 ALLOCATE (qout2(SIZE(qtot))) 384 END IF 385 qout2 = qtot 386 END IF 387 CALL print_atomic_charges(particle_set, qs_kind_set, iw, title=" DDAP "// & 388 TRIM(type_of_density)//" charges:", atomic_charges=qtot) 389 CALL timestop(handle2) 390 ! 391 ! If requested evaluate also the correction to derivatives due to Pulay Forces 392 ! 393 IF (need_f) THEN 394 CALL timeset(routineN//"-forces", handle3) 395 IF (iw > 0) THEN 396 WRITE (iw, '(T3,A)') " Evaluating DDAPC atomic derivatives .." 397 END IF 398 ALLOCATE (dAm(ndim, ndim, 3)) 399 ALLOCATE (dbv(ndim, 3)) 400 ALLOCATE (dqv(ndim, SIZE(particle_set), 3)) 401 !NB refactor math in inner loop - no more dqv0, but new temporaries instead 402 ALLOCATE (cvT_AmI(ndim)) 403 ALLOCATE (cvT_AmI_dAmj(ndim)) 404 ALLOCATE (tv(ndim, SIZE(particle_set), 3)) 405 ALLOCATE (AmI_cv(ndim)) 406 cvT_AmI(:) = MATMUL(cv, cp_ddapc_env%AmI) 407 AmI_cv(:) = MATMUL(cp_ddapc_env%AmI, cv) 408 ALLOCATE (dAmj_qv(ndim)) 409 ALLOCATE (AmI_bv(ndim)) 410 AmI_bv(:) = MATMUL(cp_ddapc_env%AmI, bv) 411 412 !NB call routine to precompute sin(g.r) and cos(g.r), 413 ! so it doesn't have to be done for each r_i-r_j pair in build_der_A_matrix_rows() 414 CALL prep_g_dot_rvec_sin_cos(rho_tot_g, particle_set, gcut, g_dot_rvec_sin, g_dot_rvec_cos) 415 !NB do build_der_A_matrix_rows in blocks, for more efficient use of DGEMM 416#define NPSET 100 417 DO iparticle0 = 1, SIZE(particle_set), NPSET 418 nparticles = MIN(NPSET, SIZE(particle_set) - iparticle0 + 1) 419 !NB each dAm is supposed to have one block of rows and one block of columns 420 !NB for derivatives with respect to each atom. build_der_A_matrix_rows() 421 !NB just returns rows, since dAm is symmetric, and missing columns can be 422 !NB reconstructed with a simple transpose, as below 423 CALL build_der_A_matrix_rows(dAm, cp_ddapc_env%gfunc, cp_ddapc_env%w, & 424 particle_set, radii, rho_tot_g, gcut, iparticle0, & 425 nparticles, g_dot_rvec_sin, g_dot_rvec_cos) 426 !NB no more reduction of dbv and dAm - instead we go through with each node's contribution 427 !NB and reduce resulting charges/forces once, at the end. Intermediate speedup can be 428 !NB had by reducing dqv after the inner loop, and then other routines don't need to know 429 !NB that contributions to dqv are distributed over the nodes. 430 !NB also get rid of zeroing of dAm and division by Vol**2 - it's slow, and can be done 431 !NB more quickly later, to a scalar or vector rather than a matrix 432 DO iparticle = iparticle0, iparticle0 + nparticles - 1 433 IF (debug_this_module) THEN 434 CALL debug_der_A_matrix(dAm, particle_set, radii, rho_tot_g, & 435 gcut, iparticle, Vol, qs_env) 436 cp_ddapc_env => qs_env%cp_ddapc_env 437 END IF 438 dbv(:, :) = 0.0_dp 439 CALL build_der_b_vector(dbv, cp_ddapc_env%gfunc, cp_ddapc_env%w, & 440 particle_set, radii, rho_tot_g, gcut, iparticle) 441 dbv(:, :) = dbv(:, :)/Vol 442 IF (debug_this_module) THEN 443 CALL debug_der_b_vector(dbv, particle_set, radii, rho_tot_g, & 444 gcut, iparticle, Vol, qs_env) 445 cp_ddapc_env => qs_env%cp_ddapc_env 446 END IF 447 DO j = 1, 3 448 !NB dAmj is actually pretty sparse - one block of cols + one block of rows - use this here: 449 pmin = (iparticle - 1)*SIZE(radii) + 1 450 pmax = iparticle*SIZE(radii) 451 !NB multiply by block of columns that aren't explicitly in dAm, but can be reconstructured 452 !NB as transpose of relevant block of rows 453 IF (pmin > 1) THEN 454 dAmj_qv(:pmin - 1) = MATMUL(TRANSPOSE(dAm(pmin:pmax, :pmin - 1, j)), qv(pmin:pmax)) 455 cvT_AmI_dAmj(:pmin - 1) = MATMUL(TRANSPOSE(dAm(pmin:pmax, :pmin - 1, j)), cvT_AmI(pmin:pmax)) 456 ENDIF 457 !NB multiply by block of rows that are explicitly in dAm 458 dAmj_qv(pmin:pmax) = MATMUL(dAm(pmin:pmax, :, j), qv(:)) 459 cvT_AmI_dAmj(pmin:pmax) = MATMUL(dAm(pmin:pmax, :, j), cvT_AmI(:)) 460 !NB multiply by block of columns that aren't explicitly in dAm, but can be reconstructured 461 !NB as transpose of relevant block of rows 462 IF (pmax < SIZE(particle_set)*SIZE(radii)) THEN 463 dAmj_qv(pmax + 1:) = MATMUL(TRANSPOSE(dAm(pmin:pmax, pmax + 1:, j)), qv(pmin:pmax)) 464 cvT_AmI_dAmj(pmax + 1:) = MATMUL(TRANSPOSE(dAm(pmin:pmax, pmax + 1:, j)), cvT_AmI(pmin:pmax)) 465 ENDIF 466 dAmj_qv(:) = dAmj_qv(:)/(Vol*Vol) 467 cvT_AmI_dAmj(:) = cvT_AmI_dAmj(:)/(Vol*Vol) 468 c3 = DOT_PRODUCT(cvT_AmI_dAmj, AmI_bv) - DOT_PRODUCT(cvT_AmI, dbv(:, j)) - c1*DOT_PRODUCT(cvT_AmI_dAmj, AmI_cv) 469 tv(:, iparticle, j) = -(dAmj_qv(:) + dbv(:, j) + c3/cp_ddapc_env%c0*cv) 470 END DO ! j 471 !NB zero relevant parts of dAm here 472 dAm((iparticle - 1)*SIZE(radii) + 1:iparticle*SIZE(radii), :, :) = 0.0_dp 473 !! dAm(:,(iparticle-1)*SIZE(radii)+1:iparticle*SIZE(radii),:) = 0.0_dp 474 END DO ! iparticle 475 END DO ! iparticle0 476 !NB final part of refactoring of math - one dgemm is faster than many dgemv 477 CALL dgemm('N', 'N', SIZE(dqv, 1), SIZE(dqv, 2)*SIZE(dqv, 3), SIZE(cp_ddapc_env%AmI, 2), 1.0_dp, & 478 cp_ddapc_env%AmI, SIZE(cp_ddapc_env%AmI, 1), tv, SIZE(tv, 1), 0.0_dp, dqv, SIZE(dqv, 1)) 479 !NB deallocate g_dot_rvec_sin and g_dot_rvec_cos 480 CALL cleanup_g_dot_rvec_sin_cos(g_dot_rvec_sin, g_dot_rvec_cos) 481 !NB moved reduction out to where dqv is used to compute 482 !NB a force contribution (smaller array to reduce, just size(particle_set) x 3) 483 !NB namely ewald_ddapc_force(), cp_decl_ddapc_forces(), restraint_functional_force() 484 CPASSERT(PRESENT(dq_out)) 485 IF (.NOT. ASSOCIATED(dq_out)) THEN 486 ALLOCATE (dq_out(SIZE(dqv, 1), SIZE(dqv, 2), SIZE(dqv, 3))) 487 ELSE 488 CPASSERT(SIZE(dqv, 1) == SIZE(dq_out, 1)) 489 CPASSERT(SIZE(dqv, 2) == SIZE(dq_out, 2)) 490 CPASSERT(SIZE(dqv, 3) == SIZE(dq_out, 3)) 491 END IF 492 dq_out = dqv 493 IF (debug_this_module) THEN 494 CALL debug_charge(dqv, qs_env, density_fit_section, & 495 particle_set, radii, rho_tot_g, type_of_density) 496 cp_ddapc_env => qs_env%cp_ddapc_env 497 END IF 498 DEALLOCATE (dqv) 499 DEALLOCATE (dAm) 500 DEALLOCATE (dbv) 501 !NB deallocate new temporaries 502 DEALLOCATE (cvT_AmI) 503 DEALLOCATE (cvT_AmI_dAmj) 504 DEALLOCATE (AmI_cv) 505 DEALLOCATE (tv) 506 DEALLOCATE (dAmj_qv) 507 DEALLOCATE (AmI_bv) 508 CALL timestop(handle3) 509 END IF 510 ! 511 ! End of charge fit 512 ! 513 DEALLOCATE (radii) 514 DEALLOCATE (bv) 515 DEALLOCATE (cv) 516 DEALLOCATE (qv) 517 DEALLOCATE (qtot) 518 IF (.NOT. PRESENT(iwc)) THEN 519 CALL cp_print_key_finished_output(iw, logger, density_fit_section, & 520 "PROGRAM_RUN_INFO") 521 END IF 522 CALL pw_pool_give_back_pw(auxbas_pool, rho_tot_g, & 523 accept_non_compatible=.TRUE.) 524 CALL timestop(handle) 525 END SUBROUTINE get_ddapc 526 527! ************************************************************************************************** 528!> \brief modify hartree potential to handle restraints in DDAPC scheme 529!> \param v_hartree_gspace ... 530!> \param density_fit_section ... 531!> \param particle_set ... 532!> \param AmI ... 533!> \param radii ... 534!> \param charges ... 535!> \param ddapc_restraint_control ... 536!> \param energy_res ... 537!> \par History 538!> 02.2006 modified [Teo] 539! ************************************************************************************************** 540 SUBROUTINE restraint_functional_potential(v_hartree_gspace, & 541 density_fit_section, particle_set, AmI, radii, charges, & 542 ddapc_restraint_control, energy_res) 543 TYPE(pw_p_type) :: v_hartree_gspace 544 TYPE(section_vals_type), POINTER :: density_fit_section 545 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 546 REAL(KIND=dp), DIMENSION(:, :), POINTER :: AmI 547 REAL(KIND=dp), DIMENSION(:), POINTER :: radii, charges 548 TYPE(ddapc_restraint_type), INTENT(INOUT) :: ddapc_restraint_control 549 REAL(KIND=dp), INTENT(INOUT) :: energy_res 550 551 CHARACTER(len=*), PARAMETER :: routineN = 'restraint_functional_potential', & 552 routineP = moduleN//':'//routineN 553 554 COMPLEX(KIND=dp) :: g_corr, phase 555 INTEGER :: handle, idim, ig, igauss, iparticle, & 556 n_gauss 557 REAL(KIND=dp) :: arg, fac, fac2, g2, gcut, gcut2, gfunc, & 558 gvec(3), rc, rc2, rvec(3), sfac, Vol, w 559 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cv, uv 560 TYPE(pw_type), POINTER :: g_hartree 561 562 CALL timeset(routineN, handle) 563 NULLIFY (g_hartree) 564 n_gauss = SIZE(radii) 565 ALLOCATE (cv(n_gauss*SIZE(particle_set))) 566 ALLOCATE (uv(n_gauss*SIZE(particle_set))) 567 uv = 0.0_dp 568 CALL evaluate_restraint_functional(ddapc_restraint_control, n_gauss, uv, & 569 charges, energy_res) 570 ! 571 CALL section_vals_val_get(density_fit_section, "GCUT", r_val=gcut) 572 gcut2 = gcut*gcut 573 g_hartree => v_hartree_gspace%pw 574 Vol = g_hartree%pw_grid%vol 575 cv = 1.0_dp/Vol 576 sfac = -1.0_dp/Vol 577 fac = DOT_PRODUCT(cv, MATMUL(AmI, cv)) 578 fac2 = DOT_PRODUCT(cv, MATMUL(AmI, uv)) 579 cv(:) = uv - cv*fac2/fac 580 cv(:) = MATMUL(AmI, cv) 581 IF (g_hartree%pw_grid%have_g0) g_hartree%cc(1) = g_hartree%cc(1) + sfac*fac2/fac 582 DO ig = g_hartree%pw_grid%first_gne0, g_hartree%pw_grid%ngpts_cut_local 583 g2 = g_hartree%pw_grid%gsq(ig) 584 w = 4.0_dp*pi*(g2 - gcut2)**2.0_dp/(g2*gcut2) 585 IF (g2 > gcut2) EXIT 586 gvec = g_hartree%pw_grid%g(:, ig) 587 g_corr = 0.0_dp 588 idim = 0 589 DO iparticle = 1, SIZE(particle_set) 590 DO igauss = 1, SIZE(radii) 591 idim = idim + 1 592 rc = radii(igauss) 593 rc2 = rc*rc 594 rvec = particle_set(iparticle)%r 595 arg = DOT_PRODUCT(gvec, rvec) 596 phase = CMPLX(COS(arg), -SIN(arg), KIND=dp) 597 gfunc = EXP(-g2*rc2/4.0_dp) 598 g_corr = g_corr + gfunc*cv(idim)*phase 599 END DO 600 END DO 601 g_corr = g_corr*w 602 g_hartree%cc(ig) = g_hartree%cc(ig) + sfac*g_corr/Vol 603 END DO 604 CALL timestop(handle) 605 END SUBROUTINE restraint_functional_potential 606 607! ************************************************************************************************** 608!> \brief Modify the Hartree potential 609!> \param v_hartree_gspace ... 610!> \param density_fit_section ... 611!> \param particle_set ... 612!> \param M ... 613!> \param AmI ... 614!> \param radii ... 615!> \param charges ... 616!> \par History 617!> 08.2005 created [tlaino] 618!> \author Teodoro Laino 619! ************************************************************************************************** 620 SUBROUTINE modify_hartree_pot(v_hartree_gspace, density_fit_section, & 621 particle_set, M, AmI, radii, charges) 622 TYPE(pw_p_type) :: v_hartree_gspace 623 TYPE(section_vals_type), POINTER :: density_fit_section 624 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 625 REAL(KIND=dp), DIMENSION(:, :), POINTER :: M, AmI 626 REAL(KIND=dp), DIMENSION(:), POINTER :: radii, charges 627 628 CHARACTER(len=*), PARAMETER :: routineN = 'modify_hartree_pot', & 629 routineP = moduleN//':'//routineN 630 631 COMPLEX(KIND=dp) :: g_corr, phase 632 INTEGER :: handle, idim, ig, igauss, iparticle 633 REAL(kind=dp) :: arg, fac, fac2, g2, gcut, gcut2, gfunc, & 634 gvec(3), rc, rc2, rvec(3), sfac, Vol, w 635 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cv, uv 636 TYPE(pw_type), POINTER :: g_hartree 637 638 CALL timeset(routineN, handle) 639 NULLIFY (g_hartree) 640 CALL section_vals_val_get(density_fit_section, "GCUT", r_val=gcut) 641 gcut2 = gcut*gcut 642 g_hartree => v_hartree_gspace%pw 643 Vol = g_hartree%pw_grid%vol 644 ALLOCATE (cv(SIZE(M, 1))) 645 ALLOCATE (uv(SIZE(M, 1))) 646 cv = 1.0_dp/Vol 647 uv(:) = MATMUL(M, charges) 648 sfac = -1.0_dp/Vol 649 fac = DOT_PRODUCT(cv, MATMUL(AmI, cv)) 650 fac2 = DOT_PRODUCT(cv, MATMUL(AmI, uv)) 651 cv(:) = uv - cv*fac2/fac 652 cv(:) = MATMUL(AmI, cv) 653 IF (g_hartree%pw_grid%have_g0) g_hartree%cc(1) = g_hartree%cc(1) + sfac*fac2/fac 654 DO ig = g_hartree%pw_grid%first_gne0, g_hartree%pw_grid%ngpts_cut_local 655 g2 = g_hartree%pw_grid%gsq(ig) 656 w = 4.0_dp*pi*(g2 - gcut2)**2.0_dp/(g2*gcut2) 657 IF (g2 > gcut2) EXIT 658 gvec = g_hartree%pw_grid%g(:, ig) 659 g_corr = 0.0_dp 660 idim = 0 661 DO iparticle = 1, SIZE(particle_set) 662 DO igauss = 1, SIZE(radii) 663 idim = idim + 1 664 rc = radii(igauss) 665 rc2 = rc*rc 666 rvec = particle_set(iparticle)%r 667 arg = DOT_PRODUCT(gvec, rvec) 668 phase = CMPLX(COS(arg), -SIN(arg), KIND=dp) 669 gfunc = EXP(-g2*rc2/4.0_dp) 670 g_corr = g_corr + gfunc*cv(idim)*phase 671 END DO 672 END DO 673 g_corr = g_corr*w 674 g_hartree%cc(ig) = g_hartree%cc(ig) + sfac*g_corr/Vol 675 END DO 676 CALL timestop(handle) 677 END SUBROUTINE modify_hartree_pot 678 679! ************************************************************************************************** 680!> \brief To Debug the derivative of the B vector for the solution of the 681!> linear system 682!> \param dbv ... 683!> \param particle_set ... 684!> \param radii ... 685!> \param rho_tot_g ... 686!> \param gcut ... 687!> \param iparticle ... 688!> \param Vol ... 689!> \param qs_env ... 690!> \par History 691!> 08.2005 created [tlaino] 692!> \author Teodoro Laino 693! ************************************************************************************************** 694 SUBROUTINE debug_der_b_vector(dbv, particle_set, radii, & 695 rho_tot_g, gcut, iparticle, Vol, qs_env) 696 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: dbv 697 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 698 REAL(KIND=dp), DIMENSION(:), POINTER :: radii 699 TYPE(pw_type), POINTER :: rho_tot_g 700 REAL(KIND=dp), INTENT(IN) :: gcut 701 INTEGER, INTENT(in) :: iparticle 702 REAL(KIND=dp), INTENT(IN) :: Vol 703 TYPE(qs_environment_type), POINTER :: qs_env 704 705 CHARACTER(len=*), PARAMETER :: routineN = 'debug_der_b_vector', & 706 routineP = moduleN//':'//routineN 707 708 INTEGER :: handle, i, kk, ndim 709 REAL(KIND=dp) :: dx, rvec(3), v0 710 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: bv1, bv2, ddbv 711 TYPE(cp_ddapc_type), POINTER :: cp_ddapc_env 712 713 NULLIFY (cp_ddapc_env) 714 CALL timeset(routineN, handle) 715 dx = 0.01_dp 716 ndim = SIZE(particle_set)*SIZE(radii) 717 ALLOCATE (bv1(ndim)) 718 ALLOCATE (bv2(ndim)) 719 ALLOCATE (ddbv(ndim)) 720 rvec = particle_set(iparticle)%r 721 cp_ddapc_env => qs_env%cp_ddapc_env 722 DO i = 1, 3 723 bv1(:) = 0.0_dp 724 bv2(:) = 0.0_dp 725 particle_set(iparticle)%r(i) = rvec(i) + dx 726 CALL build_b_vector(bv1, cp_ddapc_env%gfunc, cp_ddapc_env%w, & 727 particle_set, radii, rho_tot_g, gcut); bv1(:) = bv1(:)/Vol 728 CALL mp_sum(bv1, rho_tot_g%pw_grid%para%group) 729 particle_set(iparticle)%r(i) = rvec(i) - dx 730 CALL build_b_vector(bv2, cp_ddapc_env%gfunc, cp_ddapc_env%w, & 731 particle_set, radii, rho_tot_g, gcut); bv2(:) = bv2(:)/Vol 732 CALL mp_sum(bv2, rho_tot_g%pw_grid%para%group) 733 ddbv(:) = (bv1(:) - bv2(:))/(2.0_dp*dx) 734 DO kk = 1, SIZE(ddbv) 735 IF (ddbv(kk) .GT. 1.0E-8_dp) THEN 736 v0 = ABS(dbv(kk, i) - ddbv(kk))/ddbv(kk)*100.0_dp 737 WRITE (*, *) "Error % on B ::", v0 738 IF (v0 .GT. 0.1_dp) THEN 739 WRITE (*, '(A,2I5,2F15.9)') "ERROR IN DERIVATIVE OF B VECTOR, IPARTICLE, ICOORD:", iparticle, i, & 740 dbv(kk, i), ddbv(kk) 741 CPABORT("") 742 END IF 743 END IF 744 END DO 745 particle_set(iparticle)%r = rvec 746 END DO 747 DEALLOCATE (bv1) 748 DEALLOCATE (bv2) 749 DEALLOCATE (ddbv) 750 CALL timestop(handle) 751 END SUBROUTINE debug_der_b_vector 752 753! ************************************************************************************************** 754!> \brief To Debug the derivative of the A matrix for the solution of the 755!> linear system 756!> \param dAm ... 757!> \param particle_set ... 758!> \param radii ... 759!> \param rho_tot_g ... 760!> \param gcut ... 761!> \param iparticle ... 762!> \param Vol ... 763!> \param qs_env ... 764!> \par History 765!> 08.2005 created [tlaino] 766!> \author Teodoro Laino 767! ************************************************************************************************** 768 SUBROUTINE debug_der_A_matrix(dAm, particle_set, radii, & 769 rho_tot_g, gcut, iparticle, Vol, qs_env) 770 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: dAm 771 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 772 REAL(KIND=dp), DIMENSION(:), POINTER :: radii 773 TYPE(pw_type), POINTER :: rho_tot_g 774 REAL(KIND=dp), INTENT(IN) :: gcut 775 INTEGER, INTENT(in) :: iparticle 776 REAL(KIND=dp), INTENT(IN) :: Vol 777 TYPE(qs_environment_type), POINTER :: qs_env 778 779 CHARACTER(len=*), PARAMETER :: routineN = 'debug_der_A_matrix', & 780 routineP = moduleN//':'//routineN 781 782 INTEGER :: handle, i, kk, ll, ndim 783 REAL(KIND=dp) :: dx, rvec(3), v0 784 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Am1, Am2, ddAm, g_dot_rvec_cos, & 785 g_dot_rvec_sin 786 TYPE(cp_ddapc_type), POINTER :: cp_ddapc_env 787 788!NB new temporaries sin(g.r) and cos(g.r), as used in get_ddapc, to speed up build_der_A_matrix() 789 790 NULLIFY (cp_ddapc_env) 791 CALL timeset(routineN, handle) 792 dx = 0.01_dp 793 ndim = SIZE(particle_set)*SIZE(radii) 794 ALLOCATE (Am1(ndim, ndim)) 795 ALLOCATE (Am2(ndim, ndim)) 796 ALLOCATE (ddAm(ndim, ndim)) 797 rvec = particle_set(iparticle)%r 798 cp_ddapc_env => qs_env%cp_ddapc_env 799 CALL prep_g_dot_rvec_sin_cos(rho_tot_g, particle_set, gcut, g_dot_rvec_sin, g_dot_rvec_cos) 800 DO i = 1, 3 801 Am1 = 0.0_dp 802 Am2 = 0.0_dp 803 particle_set(iparticle)%r(i) = rvec(i) + dx 804 CALL build_A_matrix(Am1, cp_ddapc_env%gfunc, cp_ddapc_env%w, & 805 particle_set, radii, rho_tot_g, gcut, g_dot_rvec_sin, g_dot_rvec_cos) 806 Am1(:, :) = Am1(:, :)/(Vol*Vol) 807 CALL mp_sum(Am1, rho_tot_g%pw_grid%para%group) 808 particle_set(iparticle)%r(i) = rvec(i) - dx 809 CALL build_A_matrix(Am2, cp_ddapc_env%gfunc, cp_ddapc_env%w, & 810 particle_set, radii, rho_tot_g, gcut, g_dot_rvec_sin, g_dot_rvec_cos) 811 Am2(:, :) = Am2(:, :)/(Vol*Vol) 812 CALL mp_sum(Am2, rho_tot_g%pw_grid%para%group) 813 ddAm(:, :) = (Am1 - Am2)/(2.0_dp*dx) 814 DO kk = 1, SIZE(ddAm, 1) 815 DO ll = 1, SIZE(ddAm, 2) 816 IF (ddAm(kk, ll) .GT. 1.0E-8_dp) THEN 817 v0 = ABS(dAm(kk, ll, i) - ddAm(kk, ll))/ddAm(kk, ll)*100.0_dp 818 WRITE (*, *) "Error % on A ::", v0, Am1(kk, ll), Am2(kk, ll), iparticle, i, kk, ll 819 IF (v0 .GT. 0.1_dp) THEN 820 WRITE (*, '(A,4I5,2F15.9)') "ERROR IN DERIVATIVE OF A MATRIX, IPARTICLE, ICOORD:", iparticle, i, kk, ll, & 821 dAm(kk, ll, i), ddAm(kk, ll) 822 CPABORT("") 823 END IF 824 END IF 825 END DO 826 END DO 827 particle_set(iparticle)%r = rvec 828 END DO 829 CALL cleanup_g_dot_rvec_sin_cos(g_dot_rvec_sin, g_dot_rvec_cos) 830 DEALLOCATE (Am1) 831 DEALLOCATE (Am2) 832 DEALLOCATE (ddAm) 833 CALL timestop(handle) 834 END SUBROUTINE debug_der_A_matrix 835 836! ************************************************************************************************** 837!> \brief To Debug the fitted charges 838!> \param dqv ... 839!> \param qs_env ... 840!> \param density_fit_section ... 841!> \param particle_set ... 842!> \param radii ... 843!> \param rho_tot_g ... 844!> \param type_of_density ... 845!> \par History 846!> 08.2005 created [tlaino] 847!> \author Teodoro Laino 848! ************************************************************************************************** 849 SUBROUTINE debug_charge(dqv, qs_env, density_fit_section, & 850 particle_set, radii, rho_tot_g, type_of_density) 851 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: dqv 852 TYPE(qs_environment_type), POINTER :: qs_env 853 TYPE(section_vals_type), POINTER :: density_fit_section 854 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 855 REAL(KIND=dp), DIMENSION(:), POINTER :: radii 856 TYPE(pw_type), POINTER :: rho_tot_g 857 CHARACTER(LEN=*) :: type_of_density 858 859 CHARACTER(len=*), PARAMETER :: routineN = 'debug_charge', routineP = moduleN//':'//routineN 860 861 INTEGER :: handle, i, iparticle, kk, ndim 862 REAL(KIND=dp) :: dx, rvec(3) 863 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ddqv 864 REAL(KIND=dp), DIMENSION(:), POINTER :: qtot1, qtot2 865 866 CALL timeset(routineN, handle) 867 WRITE (*, *) "DEBUG_CHARGE_ROUTINE" 868 ndim = SIZE(particle_set)*SIZE(radii) 869 NULLIFY (qtot1, qtot2) 870 ALLOCATE (qtot1(ndim)) 871 ALLOCATE (qtot2(ndim)) 872 ALLOCATE (ddqv(ndim)) 873 ! 874 dx = 0.001_dp 875 DO iparticle = 1, SIZE(particle_set) 876 rvec = particle_set(iparticle)%r 877 DO i = 1, 3 878 particle_set(iparticle)%r(i) = rvec(i) + dx 879 CALL get_ddapc(qs_env, .FALSE., density_fit_section, qout1=qtot1, & 880 ext_rho_tot_g=rho_tot_g, Itype_of_density=type_of_density) 881 particle_set(iparticle)%r(i) = rvec(i) - dx 882 CALL get_ddapc(qs_env, .FALSE., density_fit_section, qout1=qtot2, & 883 ext_rho_tot_g=rho_tot_g, Itype_of_density=type_of_density) 884 ddqv(:) = (qtot1 - qtot2)/(2.0_dp*dx) 885 DO kk = 1, SIZE(qtot1) - 1, SIZE(radii) 886 IF (ANY(ddqv(kk:kk + 2) .GT. 1.0E-8_dp)) THEN 887 WRITE (*, '(A,2F12.6,F12.2)') "Error :", SUM(dqv(kk:kk + 2, iparticle, i)), SUM(ddqv(kk:kk + 2)), & 888 ABS((SUM(ddqv(kk:kk + 2)) - SUM(dqv(kk:kk + 2, iparticle, i)))/SUM(ddqv(kk:kk + 2))*100.0_dp) 889 END IF 890 END DO 891 particle_set(iparticle)%r = rvec 892 END DO 893 END DO 894 ! 895 DEALLOCATE (qtot1) 896 DEALLOCATE (qtot2) 897 DEALLOCATE (ddqv) 898 CALL timestop(handle) 899 END SUBROUTINE debug_charge 900 901END MODULE cp_ddapc_util 902