1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \par History 8!> CJM APRIL-30-2009: only uses fist_env 9!> Teodoro Laino [tlaino] - 05.2009 : Generalization to different Ewald 10!> methods (initial framework) 11!> \author CJM 12! ************************************************************************************************** 13 14MODULE fist_pol_scf 15 USE atomic_kind_types, ONLY: atomic_kind_type,& 16 get_atomic_kind 17 USE cell_types, ONLY: cell_type 18 USE cp_log_handling, ONLY: cp_get_default_logger,& 19 cp_logger_type 20 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 21 cp_print_key_unit_nr 22 USE distribution_1d_types, ONLY: distribution_1d_type 23 USE ewald_environment_types, ONLY: ewald_env_get,& 24 ewald_environment_type 25 USE ewald_pw_types, ONLY: ewald_pw_type 26 USE ewalds_multipole, ONLY: ewald_multipole_evaluate 27 USE fist_energy_types, ONLY: fist_energy_type 28 USE fist_nonbond_env_types, ONLY: fist_nonbond_env_type 29 USE input_constants, ONLY: do_fist_pol_cg,& 30 do_fist_pol_sc 31 USE input_section_types, ONLY: section_vals_type 32 USE kinds, ONLY: dp 33 USE message_passing, ONLY: mp_sum 34 USE multipole_types, ONLY: multipole_type 35 USE particle_types, ONLY: particle_type 36 USE pw_poisson_types, ONLY: do_ewald_ewald,& 37 do_ewald_pme,& 38 do_ewald_spme 39#include "./base/base_uses.f90" 40 41 IMPLICIT NONE 42 PRIVATE 43 LOGICAL, PRIVATE :: debug_this_module = .FALSE. 44 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fist_pol_scf' 45 46 PUBLIC :: fist_pol_evaluate 47 48CONTAINS 49 50! ************************************************************************************************** 51!> \brief Main driver for evaluating energy/forces in a polarizable forcefield 52!> \param atomic_kind_set ... 53!> \param multipoles ... 54!> \param ewald_env ... 55!> \param ewald_pw ... 56!> \param fist_nonbond_env ... 57!> \param cell ... 58!> \param particle_set ... 59!> \param local_particles ... 60!> \param thermo ... 61!> \param vg_coulomb ... 62!> \param pot_nonbond ... 63!> \param f_nonbond ... 64!> \param fg_coulomb ... 65!> \param use_virial ... 66!> \param pv_g ... 67!> \param pv_nonbond ... 68!> \param mm_section ... 69!> \param do_ipol ... 70!> \author Toon.Verstraelen@gmail.com (2010-03-01) 71! ************************************************************************************************** 72 SUBROUTINE fist_pol_evaluate(atomic_kind_set, multipoles, ewald_env, & 73 ewald_pw, fist_nonbond_env, cell, particle_set, local_particles, & 74 thermo, vg_coulomb, pot_nonbond, f_nonbond, fg_coulomb, use_virial, & 75 pv_g, pv_nonbond, mm_section, do_ipol) 76 77 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 78 TYPE(multipole_type), POINTER :: multipoles 79 TYPE(ewald_environment_type), POINTER :: ewald_env 80 TYPE(ewald_pw_type), POINTER :: ewald_pw 81 TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env 82 TYPE(cell_type), POINTER :: cell 83 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 84 TYPE(distribution_1d_type), POINTER :: local_particles 85 TYPE(fist_energy_type), POINTER :: thermo 86 REAL(KIND=dp) :: vg_coulomb, pot_nonbond 87 REAL(KIND=dp), DIMENSION(:, :) :: f_nonbond, fg_coulomb 88 LOGICAL, INTENT(IN) :: use_virial 89 REAL(KIND=dp), DIMENSION(3, 3) :: pv_g, pv_nonbond 90 TYPE(section_vals_type), POINTER :: mm_section 91 INTEGER :: do_ipol 92 93 CHARACTER(len=*), PARAMETER :: routineN = 'fist_pol_evaluate', & 94 routineP = moduleN//':'//routineN 95 96 SELECT CASE (do_ipol) 97 CASE (do_fist_pol_sc) 98 CALL fist_pol_evaluate_sc(atomic_kind_set, multipoles, ewald_env, & 99 ewald_pw, fist_nonbond_env, cell, particle_set, local_particles, & 100 thermo, vg_coulomb, pot_nonbond, f_nonbond, fg_coulomb, use_virial, & 101 pv_g, pv_nonbond, mm_section) 102 CASE (do_fist_pol_cg) 103 CALL fist_pol_evaluate_cg(atomic_kind_set, multipoles, ewald_env, & 104 ewald_pw, fist_nonbond_env, cell, particle_set, local_particles, & 105 thermo, vg_coulomb, pot_nonbond, f_nonbond, fg_coulomb, use_virial, & 106 pv_g, pv_nonbond, mm_section) 107 END SELECT 108 109 END SUBROUTINE fist_pol_evaluate 110 111! ************************************************************************************************** 112!> \brief Self-consistent solver for a polarizable force-field 113!> \param atomic_kind_set ... 114!> \param multipoles ... 115!> \param ewald_env ... 116!> \param ewald_pw ... 117!> \param fist_nonbond_env ... 118!> \param cell ... 119!> \param particle_set ... 120!> \param local_particles ... 121!> \param thermo ... 122!> \param vg_coulomb ... 123!> \param pot_nonbond ... 124!> \param f_nonbond ... 125!> \param fg_coulomb ... 126!> \param use_virial ... 127!> \param pv_g ... 128!> \param pv_nonbond ... 129!> \param mm_section ... 130!> \author Toon.Verstraelen@gmail.com (2010-03-01) 131!> \note 132!> Method: Given an initial guess of the induced dipoles, the electrostatic 133!> field is computed at each dipole. Then new induced dipoles are computed 134!> following p = alpha x E. This is repeated until a convergence criterion is 135!> met. The convergence is measured as the RSMD of the derivatives of the 136!> electrostatic energy (including dipole self-energy) towards the components 137!> of the dipoles. 138! ************************************************************************************************** 139 SUBROUTINE fist_pol_evaluate_sc(atomic_kind_set, multipoles, ewald_env, ewald_pw, & 140 fist_nonbond_env, cell, particle_set, local_particles, thermo, vg_coulomb, & 141 pot_nonbond, f_nonbond, fg_coulomb, use_virial, pv_g, pv_nonbond, mm_section) 142 143 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 144 TYPE(multipole_type), POINTER :: multipoles 145 TYPE(ewald_environment_type), POINTER :: ewald_env 146 TYPE(ewald_pw_type), POINTER :: ewald_pw 147 TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env 148 TYPE(cell_type), POINTER :: cell 149 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 150 TYPE(distribution_1d_type), POINTER :: local_particles 151 TYPE(fist_energy_type), POINTER :: thermo 152 REAL(KIND=dp) :: vg_coulomb, pot_nonbond 153 REAL(KIND=dp), DIMENSION(:, :) :: f_nonbond, fg_coulomb 154 LOGICAL, INTENT(IN) :: use_virial 155 REAL(KIND=dp), DIMENSION(3, 3) :: pv_g, pv_nonbond 156 TYPE(section_vals_type), POINTER :: mm_section 157 158 CHARACTER(len=*), PARAMETER :: routineN = 'fist_pol_evaluate_sc', & 159 routineP = moduleN//':'//routineN 160 161 INTEGER :: ewald_type, handle, i, iatom, ii, ikind, & 162 iter, iw, iw2, j, max_ipol_iter, & 163 natom_of_kind, natoms, nkind, ntot 164 INTEGER, DIMENSION(:), POINTER :: atom_list 165 LOGICAL :: iwarn 166 REAL(KIND=dp) :: apol, cpol, eps_pol, pot_nonbond_local, & 167 rmsd, tmp_trace 168 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: efield1, efield2 169 TYPE(atomic_kind_type), POINTER :: atomic_kind 170 TYPE(cp_logger_type), POINTER :: logger 171 172 CALL timeset(routineN, handle) 173 NULLIFY (logger, atomic_kind) 174 logger => cp_get_default_logger() 175 176 iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%ITER_INFO", & 177 extension=".mmLog") 178 iw2 = cp_print_key_unit_nr(logger, mm_section, "PRINT%EWALD_INFO", & 179 extension=".mmLog") 180 181 CALL ewald_env_get(ewald_env, max_ipol_iter=max_ipol_iter, eps_pol=eps_pol, & 182 ewald_type=ewald_type) 183 184 natoms = SIZE(particle_set) 185 ALLOCATE (efield1(3, natoms)) 186 ALLOCATE (efield2(9, natoms)) 187 188 nkind = SIZE(atomic_kind_set) 189 IF (iw > 0) WRITE (iw, FMT='(/,T5,"POL_SCF|","Method: self-consistent")') 190 IF (iw > 0) WRITE (iw, FMT='(T5,"POL_SCF|"," Iteration",7X,"Conv.",T49,"Electrostatic & Induction Energy")') 191 pol_scf: DO iter = 1, max_ipol_iter 192 ! Evaluate the electrostatic with Ewald schemes 193 CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, & 194 particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, & 195 multipoles, do_correction_bonded=.TRUE., do_forces=.FALSE., & 196 do_stress=.FALSE., do_efield=.TRUE., iw2=iw2, do_debug=.FALSE., & 197 atomic_kind_set=atomic_kind_set, mm_section=mm_section, & 198 efield1=efield1, efield2=efield2) 199 CALL mp_sum(pot_nonbond_local, logger%para_env%group) 200 201 ! compute the new dipoles, qudrupoles, and check for convergence 202 ntot = 0 203 rmsd = 0.0_dp 204 thermo%e_induction = 0.0_dp 205 DO ikind = 1, nkind 206 atomic_kind => atomic_kind_set(ikind) 207 CALL get_atomic_kind(atomic_kind, apol=apol, cpol=cpol, natom=natom_of_kind, atom_list=atom_list) 208 ! ignore atoms with dipole and quadrupole polarizability zero 209 IF (apol == 0 .AND. cpol == 0) CYCLE 210 ! increment counter correctly 211 IF (apol /= 0) ntot = ntot + natom_of_kind 212 IF (cpol /= 0) ntot = ntot + natom_of_kind 213 214 DO iatom = 1, natom_of_kind 215 ii = atom_list(iatom) 216 IF (apol /= 0) THEN 217 DO i = 1, 3 218 ! the rmsd of the derivatives of the energy towards the 219 ! components of the atomic dipole moments 220 rmsd = rmsd + (multipoles%dipoles(i, ii)/apol - efield1(i, ii))**2 221 END DO 222 END IF 223 IF (cpol /= 0) THEN 224 rmsd = rmsd + (multipoles%quadrupoles(1, 1, ii)/cpol - efield2(1, ii))**2 225 rmsd = rmsd + (multipoles%quadrupoles(2, 1, ii)/cpol - efield2(2, ii))**2 226 rmsd = rmsd + (multipoles%quadrupoles(3, 1, ii)/cpol - efield2(3, ii))**2 227 rmsd = rmsd + (multipoles%quadrupoles(1, 2, ii)/cpol - efield2(4, ii))**2 228 rmsd = rmsd + (multipoles%quadrupoles(2, 2, ii)/cpol - efield2(5, ii))**2 229 rmsd = rmsd + (multipoles%quadrupoles(3, 2, ii)/cpol - efield2(6, ii))**2 230 rmsd = rmsd + (multipoles%quadrupoles(1, 3, ii)/cpol - efield2(7, ii))**2 231 rmsd = rmsd + (multipoles%quadrupoles(2, 3, ii)/cpol - efield2(8, ii))**2 232 rmsd = rmsd + (multipoles%quadrupoles(3, 3, ii)/cpol - efield2(9, ii))**2 233 END IF 234! compute dipole 235 multipoles%dipoles(:, ii) = apol*efield1(:, ii) 236! compute quadrupole 237 IF (cpol /= 0) THEN 238 multipoles%quadrupoles(1, 1, ii) = cpol*efield2(1, ii) 239 multipoles%quadrupoles(2, 1, ii) = cpol*efield2(2, ii) 240 multipoles%quadrupoles(3, 1, ii) = cpol*efield2(3, ii) 241 multipoles%quadrupoles(1, 2, ii) = cpol*efield2(4, ii) 242 multipoles%quadrupoles(2, 2, ii) = cpol*efield2(5, ii) 243 multipoles%quadrupoles(3, 2, ii) = cpol*efield2(6, ii) 244 multipoles%quadrupoles(1, 3, ii) = cpol*efield2(7, ii) 245 multipoles%quadrupoles(2, 3, ii) = cpol*efield2(8, ii) 246 multipoles%quadrupoles(3, 3, ii) = cpol*efield2(9, ii) 247 END IF 248 ! Compute the new induction term while we are here 249 IF (apol /= 0) THEN 250 thermo%e_induction = thermo%e_induction + & 251 DOT_PRODUCT(multipoles%dipoles(:, ii), & 252 multipoles%dipoles(:, ii))/apol/2.0_dp 253 END IF 254 IF (cpol /= 0) THEN 255 tmp_trace = 0._dp 256 DO i = 1, 3 257 DO j = 1, 3 258 tmp_trace = tmp_trace + & 259 multipoles%quadrupoles(i, j, ii)*multipoles%quadrupoles(i, j, ii) 260 END DO 261 END DO 262 thermo%e_induction = thermo%e_induction + tmp_trace/cpol/6.0_dp 263 END IF 264 END DO 265 END DO 266 rmsd = SQRT(rmsd/REAL(ntot, KIND=dp)) 267 IF (iw > 0) THEN 268 ! print the energy that is minimized (this is electrostatic + induction) 269 WRITE (iw, FMT='(T5,"POL_SCF|",5X,I5,5X,E12.6,T61,F20.10)') iter, & 270 rmsd, vg_coulomb + pot_nonbond_local + thermo%e_induction 271 END IF 272 IF (rmsd <= eps_pol) THEN 273 IF (iw > 0) WRITE (iw, FMT='(T5,"POL_SCF|",1X,"Self-consistent Polarization achieved.")') 274 EXIT pol_scf 275 END IF 276 277 iwarn = ((rmsd > eps_pol) .AND. (iter == max_ipol_iter)) 278 IF (iwarn .AND. iw > 0) WRITE (iw, FMT='(T5,"POL_SCF|",1X,"Self-consistent Polarization not converged!")') 279 IF (iwarn) & 280 CPWARN("Self-consistent Polarization not converged! ") 281 END DO pol_scf 282 283 ! Now evaluate after convergence to obtain forces and converged energies 284 CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, & 285 particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, & 286 multipoles, do_correction_bonded=.TRUE., do_forces=.TRUE., & 287 do_stress=use_virial, do_efield=.FALSE., iw2=iw2, do_debug=.FALSE., & 288 atomic_kind_set=atomic_kind_set, mm_section=mm_section, & 289 forces_local=fg_coulomb, forces_glob=f_nonbond, & 290 pv_local=pv_g, pv_glob=pv_nonbond) 291 pot_nonbond = pot_nonbond + pot_nonbond_local 292 CALL mp_sum(pot_nonbond_local, logger%para_env%group) 293 294 IF (iw > 0) THEN 295 ! print the energy that is minimized (this is electrostatic + induction) 296 WRITE (iw, FMT='(T5,"POL_SCF|",5X,"Final",T61,F20.10,/)') & 297 vg_coulomb + pot_nonbond_local + thermo%e_induction 298 END IF 299 300 ! Deallocate working arrays 301 DEALLOCATE (efield1) 302 DEALLOCATE (efield2) 303 CALL cp_print_key_finished_output(iw2, logger, mm_section, & 304 "PRINT%EWALD_INFO") 305 CALL cp_print_key_finished_output(iw, logger, mm_section, & 306 "PRINT%ITER_INFO") 307 308 CALL timestop(handle) 309 END SUBROUTINE fist_pol_evaluate_sc 310 311! ************************************************************************************************** 312!> \brief Conjugate-gradient solver for a polarizable force-field 313!> \param atomic_kind_set ... 314!> \param multipoles ... 315!> \param ewald_env ... 316!> \param ewald_pw ... 317!> \param fist_nonbond_env ... 318!> \param cell ... 319!> \param particle_set ... 320!> \param local_particles ... 321!> \param thermo ... 322!> \param vg_coulomb ... 323!> \param pot_nonbond ... 324!> \param f_nonbond ... 325!> \param fg_coulomb ... 326!> \param use_virial ... 327!> \param pv_g ... 328!> \param pv_nonbond ... 329!> \param mm_section ... 330!> \author Toon.Verstraelen@gmail.com (2010-03-01) 331!> \note 332!> Method: The dipoles are found by minimizing the sum of the electrostatic 333!> and the induction energy directly using a conjugate gradient method. This 334!> routine assumes that the energy is a quadratic function of the dipoles. 335!> Finding the minimum is then done by solving a linear system. This will 336!> not work for polarizable force fields that include hyperpolarizability. 337!> 338!> The implementation of the conjugate gradient solver for linear systems 339!> is described in chapter 2.7 Sparse Linear Systems, under the section 340!> "Conjugate Gradient Method for a Sparse System". Although the inducible 341!> dipoles are the solution of a dense linear system, the same algorithm is 342!> still recommended for this situation. One does not have access to the 343!> entire hardness kernel to compute the solution with conventional linear 344!> algebra routines, but one only has a function that computes the dot 345!> product of the hardness kernel and a vector. (This is the routine that 346!> computes the electrostatic field at the atoms for a given vector of 347!> inducible dipoles.) Given such function, the conjugate gradient method 348!> is an efficient way to compute the solution of a linear system, and it 349!> scales well towards many degrees of freedom in terms of efficiency and 350!> memory usage. 351! ************************************************************************************************** 352 SUBROUTINE fist_pol_evaluate_cg(atomic_kind_set, multipoles, ewald_env, ewald_pw, & 353 fist_nonbond_env, cell, particle_set, local_particles, thermo, vg_coulomb, & 354 pot_nonbond, f_nonbond, fg_coulomb, use_virial, pv_g, pv_nonbond, mm_section) 355 356 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 357 TYPE(multipole_type), POINTER :: multipoles 358 TYPE(ewald_environment_type), POINTER :: ewald_env 359 TYPE(ewald_pw_type), POINTER :: ewald_pw 360 TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env 361 TYPE(cell_type), POINTER :: cell 362 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 363 TYPE(distribution_1d_type), POINTER :: local_particles 364 TYPE(fist_energy_type), POINTER :: thermo 365 REAL(KIND=dp) :: vg_coulomb, pot_nonbond 366 REAL(KIND=dp), DIMENSION(:, :) :: f_nonbond, fg_coulomb 367 LOGICAL, INTENT(IN) :: use_virial 368 REAL(KIND=dp), DIMENSION(3, 3) :: pv_g, pv_nonbond 369 TYPE(section_vals_type), POINTER :: mm_section 370 371 CHARACTER(len=*), PARAMETER :: routineN = 'fist_pol_evaluate_cg', & 372 routineP = moduleN//':'//routineN 373 374 INTEGER :: ewald_type, handle, i, iatom, ii, ikind, & 375 iter, iw, iw2, max_ipol_iter, & 376 natom_of_kind, natoms, nkind, ntot 377 INTEGER, DIMENSION(:), POINTER :: atom_list 378 LOGICAL :: iwarn 379 REAL(KIND=dp) :: alpha, apol, beta, denom, eps_pol, & 380 pot_nonbond_local, rmsd 381 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: conjugate, conjugate_applied, efield1, & 382 efield1_ext, residual, tmp_dipoles 383 TYPE(atomic_kind_type), POINTER :: atomic_kind 384 TYPE(cp_logger_type), POINTER :: logger 385 386 CALL timeset(routineN, handle) 387 NULLIFY (logger, atomic_kind) 388 logger => cp_get_default_logger() 389 390 iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%ITER_INFO", & 391 extension=".mmLog") 392 iw2 = cp_print_key_unit_nr(logger, mm_section, "PRINT%EWALD_INFO", & 393 extension=".mmLog") 394 395 CALL ewald_env_get(ewald_env, max_ipol_iter=max_ipol_iter, eps_pol=eps_pol, & 396 ewald_type=ewald_type) 397 398 ! allocate work arrays 399 natoms = SIZE(particle_set) 400 ALLOCATE (efield1(3, natoms)) 401 ALLOCATE (tmp_dipoles(3, natoms)) 402 ALLOCATE (residual(3, natoms)) 403 ALLOCATE (conjugate(3, natoms)) 404 ALLOCATE (conjugate_applied(3, natoms)) 405 ALLOCATE (efield1_ext(3, natoms)) 406 407 ! Compute the 'external' electrostatic field (all inducible dipoles 408 ! equal to zero). this is required for the conjugate gradient solver. 409 ! We assume that all dipoles are inducible dipoles. 410 tmp_dipoles(:, :) = multipoles%dipoles ! backup of dipoles 411 multipoles%dipoles = 0.0_dp 412 CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, & 413 particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, & 414 multipoles, do_correction_bonded=.TRUE., do_forces=.FALSE., & 415 do_stress=.FALSE., do_efield=.TRUE., iw2=iw2, do_debug=.FALSE., & 416 atomic_kind_set=atomic_kind_set, mm_section=mm_section, & 417 efield1=efield1_ext) 418 multipoles%dipoles = tmp_dipoles ! restore backup 419 420 ! Compute the electric field with the initial guess of the dipoles. 421 CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, & 422 particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, & 423 multipoles, do_correction_bonded=.TRUE., do_forces=.FALSE., & 424 do_stress=.FALSE., do_efield=.TRUE., iw2=iw2, do_debug=.FALSE., & 425 atomic_kind_set=atomic_kind_set, mm_section=mm_section, & 426 efield1=efield1) 427 428 ! Compute the first residual explicitly. 429 nkind = SIZE(atomic_kind_set) 430 ntot = 0 431 residual = 0.0_dp 432 DO ikind = 1, nkind 433 atomic_kind => atomic_kind_set(ikind) 434 CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list) 435 ! ignore atoms with polarizability zero 436 IF (apol == 0) CYCLE 437 ntot = ntot + natom_of_kind 438 DO iatom = 1, natom_of_kind 439 ii = atom_list(iatom) 440 DO i = 1, 3 441 ! residual = b - A x 442 residual(i, ii) = efield1(i, ii) - multipoles%dipoles(i, ii)/apol 443 END DO 444 END DO 445 END DO 446 ! The first conjugate residual is equal to the residual. 447 conjugate(:, :) = residual 448 449 IF (iw > 0) WRITE (iw, FMT='(/,T5,"POL_SCF|","Method: conjugate-gradient")') 450 IF (iw > 0) WRITE (iw, FMT='(T5,"POL_SCF|"," Iteration",7X,"Conv.",T49,"Electrostatic & Induction Energy")') 451 pol_scf: DO iter = 1, max_ipol_iter 452 IF (debug_this_module) THEN 453 ! In principle the residual must not be computed explicitly any more. It 454 ! is obtained in an indirect way below. When the DEBUG flag is set, the 455 ! explicit residual is computed and compared with the implicitly derived 456 ! residual as a double check. 457 CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, & 458 particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, & 459 multipoles, do_correction_bonded=.TRUE., do_forces=.FALSE., & 460 do_stress=.FALSE., do_efield=.TRUE., iw2=iw2, do_debug=.FALSE., & 461 atomic_kind_set=atomic_kind_set, mm_section=mm_section, & 462 efield1=efield1) 463 ! inapropriate use of denom to check the error on the residual 464 denom = 0.0_dp 465 END IF 466 rmsd = 0.0_dp 467 ! Compute the rmsd of the residual. 468 DO ikind = 1, nkind 469 atomic_kind => atomic_kind_set(ikind) 470 CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list) 471 ! ignore atoms with polarizability zero 472 IF (apol == 0) CYCLE 473 DO iatom = 1, natom_of_kind 474 ii = atom_list(iatom) 475 DO i = 1, 3 476 ! residual = b - A x 477 rmsd = rmsd + residual(i, ii)**2 478 IF (debug_this_module) THEN 479 denom = denom + (residual(i, ii) - (efield1(i, ii) - & 480 multipoles%dipoles(i, ii)/apol))**2 481 END IF 482 END DO 483 END DO 484 END DO 485 rmsd = SQRT(rmsd/ntot) 486 IF (iw > 0) THEN 487 WRITE (iw, FMT='(T5,"POL_SCF|",5X,I5,5X,E12.6,T67,"(not computed)")') iter, rmsd 488 IF (debug_this_module) THEN 489 denom = SQRT(denom/ntot) 490 WRITE (iw, FMT='(T5,"POL_SCF|",5X,"Error on implicit residual:",5X,E12.6)') denom 491 END IF 492 END IF 493 494 ! Apply the hardness kernel to the conjugate residual. 495 ! We assume that all non-zero dipoles are inducible dipoles. 496 tmp_dipoles(:, :) = multipoles%dipoles ! backup of dipoles 497 multipoles%dipoles = conjugate 498 CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, & 499 particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, & 500 multipoles, do_correction_bonded=.TRUE., do_forces=.FALSE., & 501 do_stress=.FALSE., do_efield=.TRUE., iw2=iw2, do_debug=.FALSE., & 502 atomic_kind_set=atomic_kind_set, mm_section=mm_section, & 503 efield1=conjugate_applied) 504 multipoles%dipoles = tmp_dipoles ! restore backup 505 conjugate_applied(:, :) = efield1_ext - conjugate_applied 506 507 ! Finish conjugate_applied and compute alpha from the conjugate gradient algorithm. 508 alpha = 0.0_dp 509 denom = 0.0_dp 510 DO ikind = 1, nkind 511 atomic_kind => atomic_kind_set(ikind) 512 CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list) 513 ! ignore atoms with polarizability zero 514 IF (apol == 0) CYCLE 515 DO iatom = 1, natom_of_kind 516 ii = atom_list(iatom) 517 DO i = 1, 3 518 conjugate_applied(i, ii) = conjugate_applied(i, ii) + conjugate(i, ii)/apol 519 END DO 520 alpha = alpha + DOT_PRODUCT(residual(:, ii), residual(:, ii)) 521 denom = denom + DOT_PRODUCT(conjugate(:, ii), conjugate_applied(:, ii)) 522 END DO 523 END DO 524 alpha = alpha/denom 525 526 ! Compute the new residual and beta from the conjugate gradient method. 527 beta = 0.0_dp 528 denom = 0.0_dp 529 DO ikind = 1, nkind 530 atomic_kind => atomic_kind_set(ikind) 531 CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list) 532 IF (apol == 0) CYCLE 533 DO iatom = 1, natom_of_kind 534 ii = atom_list(iatom) 535 denom = denom + DOT_PRODUCT(residual(:, ii), residual(:, ii)) 536 DO i = 1, 3 537 residual(i, ii) = residual(i, ii) - alpha*conjugate_applied(i, ii) 538 END DO 539 beta = beta + DOT_PRODUCT(residual(:, ii), residual(:, ii)) 540 END DO 541 END DO 542 beta = beta/denom 543 544 ! Compute the new dipoles, the new conjugate residual, and the induction 545 ! energy. 546 thermo%e_induction = 0.0_dp 547 DO ikind = 1, nkind 548 atomic_kind => atomic_kind_set(ikind) 549 CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list) 550 ! ignore atoms with polarizability zero 551 IF (apol == 0) CYCLE 552 DO iatom = 1, natom_of_kind 553 ii = atom_list(iatom) 554 DO i = 1, 3 555 multipoles%dipoles(i, ii) = multipoles%dipoles(i, ii) + alpha*conjugate(i, ii) 556 conjugate(i, ii) = residual(i, ii) + beta*conjugate(i, ii) 557 thermo%e_induction = thermo%e_induction + multipoles%dipoles(i, ii)**2/apol/2.0_dp 558 END DO 559 END DO 560 END DO 561 562 ! Quit if rmsd is low enough. 563 IF (rmsd <= eps_pol) THEN 564 IF (iw > 0) WRITE (iw, FMT='(T5,"POL_SCF|",1X,"Self-consistent Polarization achieved.")') 565 EXIT pol_scf 566 END IF 567 568 ! Print warning when not converged. 569 iwarn = ((rmsd > eps_pol) .AND. (iter == max_ipol_iter)) 570 IF (iwarn .AND. iw > 0) WRITE (iw, FMT='(T5,"POL_SCF|",1X,"Self-consistent Polarization not converged!")') 571 IF (iwarn) & 572 CPWARN("Self-consistent Polarization not converged! ") 573 END DO pol_scf 574 575 IF (debug_this_module) THEN 576 ! Now evaluate after convergence to obtain forces and converged energies 577 CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, & 578 particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, & 579 multipoles, do_correction_bonded=.TRUE., do_forces=.TRUE., & 580 do_stress=use_virial, do_efield=.TRUE., iw2=iw2, do_debug=.FALSE., & 581 atomic_kind_set=atomic_kind_set, mm_section=mm_section, & 582 forces_local=fg_coulomb, forces_glob=f_nonbond, & 583 pv_local=pv_g, pv_glob=pv_nonbond, efield1=efield1) 584 585 ! Do a final check on the convergence: compute the residual explicitely 586 rmsd = 0.0_dp 587 DO ikind = 1, nkind 588 atomic_kind => atomic_kind_set(ikind) 589 CALL get_atomic_kind(atomic_kind, apol=apol, natom=natom_of_kind, atom_list=atom_list) 590 ! ignore atoms with polarizability zero 591 IF (apol == 0) CYCLE 592 DO iatom = 1, natom_of_kind 593 ii = atom_list(iatom) 594 DO i = 1, 3 595 ! residual = b - A x 596 rmsd = rmsd + (efield1(i, ii) - multipoles%dipoles(i, ii)/apol)**2 597 END DO 598 END DO 599 END DO 600 rmsd = SQRT(rmsd/ntot) 601 IF (iw > 0) WRITE (iw, FMT='(T5,"POL_SCF|",1X,"Final RMSD of residual:",5X,E12.6)') rmsd 602 ! Stop program when congergence is not reached after all 603 IF (rmsd > eps_pol) THEN 604 CPWARN("Error in the conjugate gradient method for self-consistent polarization! ") 605 END IF 606 ELSE 607 ! Now evaluate after convergence to obtain forces and converged energies 608 CALL eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, cell, & 609 particle_set, local_particles, vg_coulomb, pot_nonbond_local, thermo, & 610 multipoles, do_correction_bonded=.TRUE., do_forces=.TRUE., & 611 do_stress=use_virial, do_efield=.FALSE., iw2=iw2, do_debug=.FALSE., & 612 atomic_kind_set=atomic_kind_set, mm_section=mm_section, & 613 forces_local=fg_coulomb, forces_glob=f_nonbond, & 614 pv_local=pv_g, pv_glob=pv_nonbond) 615 ENDIF 616 pot_nonbond = pot_nonbond + pot_nonbond_local 617 CALL mp_sum(pot_nonbond_local, logger%para_env%group) 618 619 IF (iw > 0) THEN 620 WRITE (iw, FMT='(T5,"POL_SCF|",5X,"Final",T61,F20.10,/)') & 621 vg_coulomb + pot_nonbond_local + thermo%e_induction 622 END IF 623 624 ! Deallocate working arrays 625 DEALLOCATE (efield1) 626 DEALLOCATE (tmp_dipoles) 627 DEALLOCATE (residual) 628 DEALLOCATE (conjugate) 629 DEALLOCATE (conjugate_applied) 630 DEALLOCATE (efield1_ext) 631 CALL cp_print_key_finished_output(iw2, logger, mm_section, & 632 "PRINT%EWALD_INFO") 633 CALL cp_print_key_finished_output(iw, logger, mm_section, & 634 "PRINT%ITER_INFO") 635 636 CALL timestop(handle) 637 END SUBROUTINE fist_pol_evaluate_cg 638 639! ************************************************************************************************** 640!> \brief Main driver for evaluating electrostatic in polarible forcefields 641!> All the dependence on the Ewald method should go here! 642!> \param ewald_type ... 643!> \param ewald_env ... 644!> \param ewald_pw ... 645!> \param fist_nonbond_env ... 646!> \param cell ... 647!> \param particle_set ... 648!> \param local_particles ... 649!> \param vg_coulomb ... 650!> \param pot_nonbond ... 651!> \param thermo ... 652!> \param multipoles ... 653!> \param do_correction_bonded ... 654!> \param do_forces ... 655!> \param do_stress ... 656!> \param do_efield ... 657!> \param iw2 ... 658!> \param do_debug ... 659!> \param atomic_kind_set ... 660!> \param mm_section ... 661!> \param efield0 ... 662!> \param efield1 ... 663!> \param efield2 ... 664!> \param forces_local ... 665!> \param forces_glob ... 666!> \param pv_local ... 667!> \param pv_glob ... 668!> \author Teodoro Laino [tlaino] 05.2009 669! ************************************************************************************************** 670 SUBROUTINE eval_pol_ewald(ewald_type, ewald_env, ewald_pw, fist_nonbond_env, & 671 cell, particle_set, local_particles, vg_coulomb, pot_nonbond, thermo, & 672 multipoles, do_correction_bonded, do_forces, do_stress, do_efield, iw2, & 673 do_debug, atomic_kind_set, mm_section, efield0, efield1, efield2, forces_local, & 674 forces_glob, pv_local, pv_glob) 675 676 INTEGER, INTENT(IN) :: ewald_type 677 TYPE(ewald_environment_type), POINTER :: ewald_env 678 TYPE(ewald_pw_type), POINTER :: ewald_pw 679 TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env 680 TYPE(cell_type), POINTER :: cell 681 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 682 TYPE(distribution_1d_type), POINTER :: local_particles 683 REAL(KIND=dp), INTENT(OUT) :: vg_coulomb, pot_nonbond 684 TYPE(fist_energy_type), POINTER :: thermo 685 TYPE(multipole_type), POINTER :: multipoles 686 LOGICAL, INTENT(IN) :: do_correction_bonded, do_forces, & 687 do_stress, do_efield 688 INTEGER, INTENT(IN) :: iw2 689 LOGICAL, INTENT(IN) :: do_debug 690 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 691 TYPE(section_vals_type), POINTER :: mm_section 692 REAL(KIND=dp), DIMENSION(:), OPTIONAL :: efield0 693 REAL(KIND=dp), DIMENSION(:, :), OPTIONAL :: efield1, efield2 694 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), & 695 OPTIONAL :: forces_local, forces_glob, pv_local, & 696 pv_glob 697 698 CHARACTER(len=*), PARAMETER :: routineN = 'eval_pol_ewald', routineP = moduleN//':'//routineN 699 700 INTEGER :: handle 701 702 CALL timeset(routineN, handle) 703 704 pot_nonbond = 0.0_dp ! Initialization.. 705 vg_coulomb = 0.0_dp ! Initialization.. 706 SELECT CASE (ewald_type) 707 CASE (do_ewald_ewald) 708 CALL ewald_multipole_evaluate(ewald_env, ewald_pw, fist_nonbond_env, cell, & 709 particle_set, local_particles, vg_coulomb, pot_nonbond, thermo%e_neut, & 710 thermo%e_self, multipoles%task, do_correction_bonded=do_correction_bonded, & 711 do_forces=do_forces, do_stress=do_stress, do_efield=do_efield, & 712 radii=multipoles%radii, charges=multipoles%charges, & 713 dipoles=multipoles%dipoles, quadrupoles=multipoles%quadrupoles, & 714 forces_local=forces_local, forces_glob=forces_glob, pv_local=pv_local, & 715 pv_glob=pv_glob, iw=iw2, do_debug=do_debug, atomic_kind_set=atomic_kind_set, & 716 mm_section=mm_section, efield0=efield0, efield1=efield1, efield2=efield2) 717 CASE (do_ewald_pme) 718 CPABORT("Multipole Ewald not yet implemented within a PME scheme!") 719 CASE (do_ewald_spme) 720 CPABORT("Multipole Ewald not yet implemented within a SPME scheme!") 721 END SELECT 722 CALL timestop(handle) 723 END SUBROUTINE eval_pol_ewald 724 725END MODULE fist_pol_scf 726