1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Treats the electrostatic for multipoles (up to quadrupoles) 8!> \author Teodoro Laino [tlaino] - 12.2007 - University of Zurich 9!> inclusion of optional electric field damping for the polarizable atoms 10!> Rodolphe Vuilleumier and Mathieu Salanne - 12.2009 11! ************************************************************************************************** 12MODULE ewalds_multipole 13 USE atomic_kind_types, ONLY: atomic_kind_type 14 USE bibliography, ONLY: Aguado2003,& 15 Laino2008,& 16 cite_reference 17 USE cell_types, ONLY: cell_type,& 18 pbc 19 USE cp_log_handling, ONLY: cp_get_default_logger,& 20 cp_logger_type 21 USE damping_dipole_types, ONLY: damping_type,& 22 no_damping,& 23 tang_toennies 24 USE dg_rho0_types, ONLY: dg_rho0_type 25 USE dg_types, ONLY: dg_get,& 26 dg_type 27 USE distribution_1d_types, ONLY: distribution_1d_type 28 USE ewald_environment_types, ONLY: ewald_env_get,& 29 ewald_environment_type 30 USE ewald_pw_types, ONLY: ewald_pw_get,& 31 ewald_pw_type 32 USE fist_neighbor_list_control, ONLY: list_control 33 USE fist_neighbor_list_types, ONLY: fist_neighbor_type,& 34 neighbor_kind_pairs_type 35 USE fist_nonbond_env_types, ONLY: fist_nonbond_env_get,& 36 fist_nonbond_env_type,& 37 pos_type 38 USE input_section_types, ONLY: section_vals_type 39 USE kinds, ONLY: dp 40 USE mathconstants, ONLY: fourpi,& 41 oorootpi,& 42 pi,& 43 sqrthalf 44 USE mathlib, ONLY: matvec_3x3 45 USE message_passing, ONLY: mp_sum 46 USE parallel_rng_types, ONLY: UNIFORM,& 47 create_rng_stream,& 48 delete_rng_stream,& 49 random_numbers,& 50 rng_stream_type 51 USE particle_types, ONLY: particle_type 52 USE pw_grid_types, ONLY: pw_grid_type 53 USE pw_pool_types, ONLY: pw_pool_type 54 USE structure_factor_types, ONLY: structure_factor_type 55 USE structure_factors, ONLY: structure_factor_allocate,& 56 structure_factor_deallocate,& 57 structure_factor_evaluate 58#include "./base/base_uses.f90" 59 60#:include "ewalds_multipole_sr.fypp" 61 62 IMPLICIT NONE 63 PRIVATE 64 65 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE. 66 LOGICAL, PRIVATE, PARAMETER :: debug_r_space = .FALSE. 67 LOGICAL, PRIVATE, PARAMETER :: debug_g_space = .FALSE. 68 LOGICAL, PRIVATE, PARAMETER :: debug_e_field = .FALSE. 69 LOGICAL, PRIVATE, PARAMETER :: debug_e_field_en = .FALSE. 70 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ewalds_multipole' 71 72 PUBLIC :: ewald_multipole_evaluate 73 74CONTAINS 75 76! ************************************************************************************************** 77!> \brief Computes the potential and the force for a lattice sum of multipoles (up to quadrupole) 78!> \param ewald_env ... 79!> \param ewald_pw ... 80!> \param nonbond_env ... 81!> \param cell ... 82!> \param particle_set ... 83!> \param local_particles ... 84!> \param energy_local ... 85!> \param energy_glob ... 86!> \param e_neut ... 87!> \param e_self ... 88!> \param task ... 89!> \param do_correction_bonded ... 90!> \param do_forces ... 91!> \param do_stress ... 92!> \param do_efield ... 93!> \param radii ... 94!> \param charges ... 95!> \param dipoles ... 96!> \param quadrupoles ... 97!> \param forces_local ... 98!> \param forces_glob ... 99!> \param pv_local ... 100!> \param pv_glob ... 101!> \param efield0 ... 102!> \param efield1 ... 103!> \param efield2 ... 104!> \param iw ... 105!> \param do_debug ... 106!> \param atomic_kind_set ... 107!> \param mm_section ... 108!> \par Note 109!> atomic_kind_set and mm_section are between the arguments only 110!> for debug purpose (therefore optional) and can be avoided when this 111!> function is called in other part of the program 112!> \par Note 113!> When a gaussian multipole is used instead of point multipole, i.e. 114!> when radii(i)>0, the electrostatic fields (efield0, efield1, efield2) 115!> become derivatives of the electrostatic potential energy towards 116!> these gaussian multipoles. 117!> \author Teodoro Laino [tlaino] - 12.2007 - University of Zurich 118! ************************************************************************************************** 119 RECURSIVE SUBROUTINE ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, & 120 cell, particle_set, local_particles, energy_local, energy_glob, e_neut, e_self, & 121 task, do_correction_bonded, do_forces, do_stress, & 122 do_efield, radii, charges, dipoles, & 123 quadrupoles, forces_local, forces_glob, pv_local, pv_glob, efield0, efield1, & 124 efield2, iw, do_debug, atomic_kind_set, mm_section) 125 TYPE(ewald_environment_type), POINTER :: ewald_env 126 TYPE(ewald_pw_type), POINTER :: ewald_pw 127 TYPE(fist_nonbond_env_type), POINTER :: nonbond_env 128 TYPE(cell_type), POINTER :: cell 129 TYPE(particle_type), POINTER :: particle_set(:) 130 TYPE(distribution_1d_type), POINTER :: local_particles 131 REAL(KIND=dp), INTENT(INOUT) :: energy_local, energy_glob 132 REAL(KIND=dp), INTENT(OUT) :: e_neut, e_self 133 LOGICAL, DIMENSION(3), INTENT(IN) :: task 134 LOGICAL, INTENT(IN) :: do_correction_bonded, do_forces, & 135 do_stress, do_efield 136 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: radii, charges 137 REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: dipoles 138 REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, & 139 POINTER :: quadrupoles 140 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), & 141 OPTIONAL :: forces_local, forces_glob, pv_local, & 142 pv_glob 143 REAL(KIND=dp), DIMENSION(:), INTENT(OUT), OPTIONAL :: efield0 144 REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT), & 145 OPTIONAL :: efield1, efield2 146 INTEGER, INTENT(IN) :: iw 147 LOGICAL, INTENT(IN) :: do_debug 148 TYPE(atomic_kind_type), DIMENSION(:), OPTIONAL, & 149 POINTER :: atomic_kind_set 150 TYPE(section_vals_type), OPTIONAL, POINTER :: mm_section 151 152 CHARACTER(len=*), PARAMETER :: routineN = 'ewald_multipole_evaluate', & 153 routineP = moduleN//':'//routineN 154 155 INTEGER :: group, handle, i, j, size1, size2 156 LOGICAL :: check_debug, check_efield, check_forces, & 157 do_task(3) 158 LOGICAL, DIMENSION(3, 3) :: my_task 159 REAL(KIND=dp) :: e_bonded, e_bonded_t, e_rspace, & 160 e_rspace_t, energy_glob_t 161 REAL(KIND=dp), DIMENSION(:), POINTER :: efield0_lr, efield0_sr 162 REAL(KIND=dp), DIMENSION(:, :), POINTER :: efield1_lr, efield1_sr, efield2_lr, & 163 efield2_sr 164 165 CALL cite_reference(Aguado2003) 166 CALL cite_reference(Laino2008) 167 CALL timeset(routineN, handle) 168 CPASSERT(ASSOCIATED(nonbond_env)) 169 check_debug = (debug_this_module .OR. debug_r_space .OR. debug_g_space .OR. debug_e_field .OR. debug_e_field_en) & 170 .EQV. debug_this_module 171 CPASSERT(check_debug) 172 check_forces = do_forces .EQV. (PRESENT(forces_local) .AND. PRESENT(forces_glob)) 173 CPASSERT(check_forces) 174 check_efield = do_efield .EQV. (PRESENT(efield0) .OR. PRESENT(efield1) .OR. PRESENT(efield2)) 175 CPASSERT(check_efield) 176 ! Debugging this module 177 IF (debug_this_module .AND. do_debug) THEN 178 ! Debug specifically real space part 179 IF (debug_r_space) THEN 180 CALL debug_ewald_multipoles(ewald_env, ewald_pw, nonbond_env, cell, & 181 particle_set, local_particles, iw, debug_r_space) 182 CPABORT("Debug Multipole Requested: Real Part!") 183 END IF 184 ! Debug electric fields and gradients as pure derivatives 185 IF (debug_e_field) THEN 186 CPASSERT(PRESENT(atomic_kind_set)) 187 CPASSERT(PRESENT(mm_section)) 188 CALL debug_ewald_multipoles_fields(ewald_env, ewald_pw, nonbond_env, & 189 cell, particle_set, local_particles, radii, charges, dipoles, & 190 quadrupoles, task, iw, atomic_kind_set, mm_section) 191 CPABORT("Debug Multipole Requested: POT+EFIELDS+GRAD!") 192 END IF 193 ! Debug the potential, electric fields and electric fields gradient in oder 194 ! to retrieve the correct energy 195 IF (debug_e_field_en) THEN 196 CALL debug_ewald_multipoles_fields2(ewald_env, ewald_pw, nonbond_env, & 197 cell, particle_set, local_particles, radii, charges, dipoles, & 198 quadrupoles, task, iw) 199 CPABORT("Debug Multipole Requested: POT+EFIELDS+GRAD to give the correct energy!!") 200 END IF 201 END IF 202 203 ! Setup the tasks (needed to skip useless parts in the real-space term) 204 do_task = task 205 DO i = 1, 3 206 IF (do_task(i)) THEN 207 SELECT CASE (i) 208 CASE (1) 209 do_task(1) = ANY(charges /= 0.0_dp) 210 CASE (2) 211 do_task(2) = ANY(dipoles /= 0.0_dp) 212 CASE (3) 213 do_task(3) = ANY(quadrupoles /= 0.0_dp) 214 END SELECT 215 END IF 216 END DO 217 DO i = 1, 3 218 DO j = i, 3 219 my_task(j, i) = do_task(i) .AND. do_task(j) 220 my_task(i, j) = my_task(j, i) 221 END DO 222 END DO 223 224 ! Allocate arrays for the evaluation of the potential, fields and electrostatic field gradients 225 NULLIFY (efield0_sr, efield0_lr, efield1_sr, efield1_lr, efield2_sr, efield2_lr) 226 IF (do_efield) THEN 227 IF (PRESENT(efield0)) THEN 228 size1 = SIZE(efield0) 229 ALLOCATE (efield0_sr(size1)) 230 ALLOCATE (efield0_lr(size1)) 231 efield0_sr = 0.0_dp 232 efield0_lr = 0.0_dp 233 END IF 234 IF (PRESENT(efield1)) THEN 235 size1 = SIZE(efield1, 1) 236 size2 = SIZE(efield1, 2) 237 ALLOCATE (efield1_sr(size1, size2)) 238 ALLOCATE (efield1_lr(size1, size2)) 239 efield1_sr = 0.0_dp 240 efield1_lr = 0.0_dp 241 END IF 242 IF (PRESENT(efield2)) THEN 243 size1 = SIZE(efield2, 1) 244 size2 = SIZE(efield2, 2) 245 ALLOCATE (efield2_sr(size1, size2)) 246 ALLOCATE (efield2_lr(size1, size2)) 247 efield2_sr = 0.0_dp 248 efield2_lr = 0.0_dp 249 END IF 250 END IF 251 252 e_rspace = 0.0_dp 253 e_bonded = 0.0_dp 254 IF ((.NOT. debug_g_space) .AND. (nonbond_env%do_nonbonded)) THEN 255 ! Compute the Real Space (Short-Range) part of the Ewald sum. 256 ! This contribution is only added when the nonbonded flag in the input 257 ! is set, because these contributions depend. the neighborlists. 258 CALL ewald_multipole_SR(nonbond_env, ewald_env, atomic_kind_set, & 259 particle_set, cell, e_rspace, my_task, & 260 do_forces, do_efield, do_stress, radii, charges, dipoles, quadrupoles, & 261 forces_glob, pv_glob, efield0_sr, efield1_sr, efield2_sr) 262 energy_glob = energy_glob+e_rspace 263 264 IF (do_correction_bonded) THEN 265 ! The corrections for bonded interactions are stored in the Real Space 266 ! (Short-Range) part of the fields array. 267 CALL ewald_multipole_bonded(nonbond_env, particle_set, ewald_env, & 268 cell, e_bonded, my_task, do_forces, do_efield, do_stress, & 269 charges, dipoles, quadrupoles, forces_glob, pv_glob, & 270 efield0_sr, efield1_sr, efield2_sr) 271 energy_glob = energy_glob+e_bonded 272 END IF 273 END IF 274 275 e_neut = 0.0_dp 276 e_self = 0.0_dp 277 energy_local = 0.0_dp 278 IF (.NOT. debug_r_space) THEN 279 ! Compute the Reciprocal Space (Long-Range) part of the Ewald sum 280 CALL ewald_multipole_LR(ewald_env, ewald_pw, cell, particle_set, & 281 local_particles, energy_local, my_task, do_forces, do_efield, do_stress, & 282 charges, dipoles, quadrupoles, forces_local, pv_local, efield0_lr, efield1_lr, & 283 efield2_lr) 284 285 ! Self-Interactions corrections 286 CALL ewald_multipole_self(ewald_env, cell, local_particles, e_self, & 287 e_neut, my_task, do_efield, radii, charges, dipoles, quadrupoles, & 288 efield0_lr, efield1_lr, efield2_lr) 289 END IF 290 291 ! Sumup energy contributions for possible IO 292 CALL ewald_env_get(ewald_env, group=group) 293 energy_glob_t = energy_glob 294 e_rspace_t = e_rspace 295 e_bonded_t = e_bonded 296 CALL mp_sum(energy_glob_t, group) 297 CALL mp_sum(e_rspace_t, group) 298 CALL mp_sum(e_bonded_t, group) 299 ! Print some info about energetics 300 CALL ewald_multipole_print(iw, energy_local, e_rspace_t, e_bonded_t, e_self, e_neut) 301 302 ! Gather the components of the potential, fields and electrostatic field gradients 303 IF (do_efield) THEN 304 IF (PRESENT(efield0)) THEN 305 efield0 = efield0_sr+efield0_lr 306 CALL mp_sum(efield0, group) 307 DEALLOCATE (efield0_sr) 308 DEALLOCATE (efield0_lr) 309 END IF 310 IF (PRESENT(efield1)) THEN 311 efield1 = efield1_sr+efield1_lr 312 CALL mp_sum(efield1, group) 313 DEALLOCATE (efield1_sr) 314 DEALLOCATE (efield1_lr) 315 END IF 316 IF (PRESENT(efield2)) THEN 317 efield2 = efield2_sr+efield2_lr 318 CALL mp_sum(efield2, group) 319 DEALLOCATE (efield2_sr) 320 DEALLOCATE (efield2_lr) 321 END IF 322 END IF 323 CALL timestop(handle) 324 END SUBROUTINE ewald_multipole_evaluate 325 326! ************************************************************************************************** 327!> \brief computes the potential and the force for a lattice sum of multipoles 328!> up to quadrupole - Short Range (Real Space) Term 329!> \param nonbond_env ... 330!> \param ewald_env ... 331!> \param atomic_kind_set ... 332!> \param particle_set ... 333!> \param cell ... 334!> \param energy ... 335!> \param task ... 336!> \param do_forces ... 337!> \param do_efield ... 338!> \param do_stress ... 339!> \param radii ... 340!> \param charges ... 341!> \param dipoles ... 342!> \param quadrupoles ... 343!> \param forces ... 344!> \param pv ... 345!> \param efield0 ... 346!> \param efield1 ... 347!> \param efield2 ... 348!> \author Teodoro Laino [tlaino] - 12.2007 - University of Zurich 349! ************************************************************************************************** 350 SUBROUTINE ewald_multipole_SR(nonbond_env, ewald_env, atomic_kind_set, & 351 particle_set, cell, energy, task, & 352 do_forces, do_efield, do_stress, radii, charges, dipoles, quadrupoles, & 353 forces, pv, efield0, efield1, efield2) 354 TYPE(fist_nonbond_env_type), POINTER :: nonbond_env 355 TYPE(ewald_environment_type), POINTER :: ewald_env 356 TYPE(atomic_kind_type), DIMENSION(:), OPTIONAL, & 357 POINTER :: atomic_kind_set 358 TYPE(particle_type), POINTER :: particle_set(:) 359 TYPE(cell_type), POINTER :: cell 360 REAL(KIND=dp), INTENT(INOUT) :: energy 361 LOGICAL, DIMENSION(3, 3), INTENT(IN) :: task 362 LOGICAL, INTENT(IN) :: do_forces, do_efield, do_stress 363 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: radii, charges 364 REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: dipoles 365 REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, & 366 POINTER :: quadrupoles 367 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), & 368 OPTIONAL :: forces, pv 369 REAL(KIND=dp), DIMENSION(:), POINTER :: efield0 370 REAL(KIND=dp), DIMENSION(:, :), POINTER :: efield1, efield2 371 372 CHARACTER(len=*), PARAMETER :: routineN = 'ewald_multipole_SR', & 373 routineP = moduleN//':'//routineN 374 375 INTEGER :: a, atom_a, atom_b, b, c, d, e, handle, i, iend, igrp, ikind, ilist, ipair, & 376 istart, itype_ij, itype_ji, jkind, k, kind_a, kind_b, kk, nkdamp_ij, nkdamp_ji, nkinds, & 377 npairs 378 INTEGER, DIMENSION(:, :), POINTER :: list 379 LOGICAL :: do_efield0, do_efield1, do_efield2, & 380 force_eval 381 REAL(KIND=dp) :: alpha, beta, ch_i, ch_j, dampa_ij, dampa_ji, dampaexpi, dampaexpj, & 382 dampfac_ij, dampfac_ji, dampfuncdiffi, dampfuncdiffj, dampfunci, dampfuncj, dampsumfi, & 383 dampsumfj, ef0_i, ef0_j, eloc, fac, fac_ij, factorial, ir, irab2, ptens11, ptens12, & 384 ptens13, ptens21, ptens22, ptens23, ptens31, ptens32, ptens33, r, rab2, rab2_max, radius, & 385 rcut, tij, tmp, tmp1, tmp11, tmp12, tmp13, tmp2, tmp21, tmp22, tmp23, tmp31, tmp32, & 386 tmp33, tmp_ij, tmp_ji, xf 387 REAL(KIND=dp), DIMENSION(0:5) :: f 388 REAL(KIND=dp), DIMENSION(3) :: cell_v, cvi, damptij_a, damptji_a, dp_i, & 389 dp_j, ef1_i, ef1_j, fr, rab, tij_a 390 REAL(KIND=dp), DIMENSION(3, 3) :: damptij_ab, damptji_ab, ef2_i, ef2_j, & 391 qp_i, qp_j, tij_ab 392 REAL(KIND=dp), DIMENSION(3, 3, 3) :: tij_abc 393 REAL(KIND=dp), DIMENSION(3, 3, 3, 3) :: tij_abcd 394 REAL(KIND=dp), DIMENSION(3, 3, 3, 3, 3) :: tij_abcde 395 TYPE(damping_type) :: damping_ij, damping_ji 396 TYPE(fist_neighbor_type), POINTER :: nonbonded 397 TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair 398 TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update, r_last_update_pbc 399 400 CALL timeset(routineN, handle) 401 NULLIFY (nonbonded, r_last_update, r_last_update_pbc) 402 do_efield0 = do_efield .AND. ASSOCIATED(efield0) 403 do_efield1 = do_efield .AND. ASSOCIATED(efield1) 404 do_efield2 = do_efield .AND. ASSOCIATED(efield2) 405 IF (do_stress) THEN 406 ptens11 = 0.0_dp; ptens12 = 0.0_dp; ptens13 = 0.0_dp 407 ptens21 = 0.0_dp; ptens22 = 0.0_dp; ptens23 = 0.0_dp 408 ptens31 = 0.0_dp; ptens32 = 0.0_dp; ptens33 = 0.0_dp 409 END IF 410 ! Get nonbond_env info 411 CALL fist_nonbond_env_get(nonbond_env, nonbonded=nonbonded, natom_types=nkinds, & 412 r_last_update=r_last_update, r_last_update_pbc=r_last_update_pbc) 413 CALL ewald_env_get(ewald_env, alpha=alpha, rcut=rcut) 414 rab2_max = rcut**2 415 IF (debug_r_space) THEN 416 rab2_max = HUGE(0.0_dp) 417 END IF 418 ! Starting the force loop 419 Lists: DO ilist = 1, nonbonded%nlists 420 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist) 421 npairs = neighbor_kind_pair%npairs 422 IF (npairs == 0) CYCLE 423 list => neighbor_kind_pair%list 424 cvi = neighbor_kind_pair%cell_vector 425 CALL matvec_3x3(cell_v, cell%hmat, cvi) 426 Kind_Group_Loop: DO igrp = 1, neighbor_kind_pair%ngrp_kind 427 istart = neighbor_kind_pair%grp_kind_start(igrp) 428 iend = neighbor_kind_pair%grp_kind_end(igrp) 429 ikind = neighbor_kind_pair%ij_kind(1, igrp) 430 jkind = neighbor_kind_pair%ij_kind(2, igrp) 431 432 itype_ij = no_damping 433 nkdamp_ij = 1 434 dampa_ij = 1.0_dp 435 dampfac_ij = 0.0_dp 436 437 itype_ji = no_damping 438 nkdamp_ji = 1 439 dampa_ji = 1.0_dp 440 dampfac_ji = 0.0_dp 441 IF (PRESENT(atomic_kind_set)) THEN 442 IF (ASSOCIATED(atomic_kind_set(jkind)%damping)) THEN 443 damping_ij = atomic_kind_set(jkind)%damping%damp(ikind) 444 itype_ij = damping_ij%itype 445 nkdamp_ij = damping_ij%order 446 dampa_ij = damping_ij%bij 447 dampfac_ij = damping_ij%cij 448 END IF 449 450 IF (ASSOCIATED(atomic_kind_set(ikind)%damping)) THEN 451 damping_ji = atomic_kind_set(ikind)%damping%damp(jkind) 452 itype_ji = damping_ji%itype 453 nkdamp_ji = damping_ji%order 454 dampa_ji = damping_ji%bij 455 dampfac_ji = damping_ji%cij 456 END IF 457 END IF 458 459 Pairs: DO ipair = istart, iend 460 IF (ipair <= neighbor_kind_pair%nscale) THEN 461 ! scale the electrostatic interaction if needed 462 ! (most often scaled to zero) 463 fac_ij = neighbor_kind_pair%ei_scale(ipair) 464 IF (fac_ij <= 0) CYCLE 465 ELSE 466 fac_ij = 1.0_dp 467 END IF 468 atom_a = list(1, ipair) 469 atom_b = list(2, ipair) 470 kind_a = particle_set(atom_a)%atomic_kind%kind_number 471 kind_b = particle_set(atom_b)%atomic_kind%kind_number 472 IF (atom_a == atom_b) fac_ij = 0.5_dp 473 rab = r_last_update_pbc(atom_b)%r-r_last_update_pbc(atom_a)%r 474 rab = rab+cell_v 475 rab2 = rab(1)**2+rab(2)**2+rab(3)**2 476 IF (rab2 <= rab2_max) THEN 477 IF (PRESENT(radii)) THEN 478 radius = SQRT(radii(atom_a)*radii(atom_a)+radii(atom_b)*radii(atom_b)) 479 ELSE 480 radius = 0.0_dp 481 END IF 482 IF (radius > 0.0_dp) THEN 483 beta = sqrthalf/radius 484$: ewalds_multipole_sr_macro(mode="SCREENED_COULOMB_GAUSS", damping=False, store_energy=True, store_forces=True) 485 ELSE 486$: ewalds_multipole_sr_macro(mode="SCREENED_COULOMB_ERFC", damping=True, store_energy=True, store_forces=True ) 487 END IF 488 END IF 489 END DO Pairs 490 END DO Kind_Group_Loop 491 END DO Lists 492 IF (do_stress) THEN 493 pv(1, 1) = pv(1, 1)+ptens11 494 pv(1, 2) = pv(1, 2)+(ptens12+ptens21)*0.5_dp 495 pv(1, 3) = pv(1, 3)+(ptens13+ptens31)*0.5_dp 496 pv(2, 1) = pv(1, 2) 497 pv(2, 2) = pv(2, 2)+ptens22 498 pv(2, 3) = pv(2, 3)+(ptens23+ptens32)*0.5_dp 499 pv(3, 1) = pv(1, 3) 500 pv(3, 2) = pv(2, 3) 501 pv(3, 3) = pv(3, 3)+ptens33 502 END IF 503 504 CALL timestop(handle) 505 END SUBROUTINE ewald_multipole_SR 506 507! ************************************************************************************************** 508!> \brief computes the bonded correction for the potential and the force for a 509!> lattice sum of multipoles up to quadrupole 510!> \param nonbond_env ... 511!> \param particle_set ... 512!> \param ewald_env ... 513!> \param cell ... 514!> \param energy ... 515!> \param task ... 516!> \param do_forces ... 517!> \param do_efield ... 518!> \param do_stress ... 519!> \param charges ... 520!> \param dipoles ... 521!> \param quadrupoles ... 522!> \param forces ... 523!> \param pv ... 524!> \param efield0 ... 525!> \param efield1 ... 526!> \param efield2 ... 527!> \author Teodoro Laino [tlaino] - 05.2009 528! ************************************************************************************************** 529 SUBROUTINE ewald_multipole_bonded(nonbond_env, particle_set, ewald_env, & 530 cell, energy, task, do_forces, do_efield, do_stress, charges, & 531 dipoles, quadrupoles, forces, pv, efield0, efield1, efield2) 532 533 TYPE(fist_nonbond_env_type), POINTER :: nonbond_env 534 TYPE(particle_type), POINTER :: particle_set(:) 535 TYPE(ewald_environment_type), POINTER :: ewald_env 536 TYPE(cell_type), POINTER :: cell 537 REAL(KIND=dp), INTENT(INOUT) :: energy 538 LOGICAL, DIMENSION(3, 3), INTENT(IN) :: task 539 LOGICAL, INTENT(IN) :: do_forces, do_efield, do_stress 540 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: charges 541 REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: dipoles 542 REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, & 543 POINTER :: quadrupoles 544 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), & 545 OPTIONAL :: forces, pv 546 REAL(KIND=dp), DIMENSION(:), POINTER :: efield0 547 REAL(KIND=dp), DIMENSION(:, :), POINTER :: efield1, efield2 548 549 CHARACTER(len=*), PARAMETER :: routineN = 'ewald_multipole_bonded', & 550 routineP = moduleN//':'//routineN 551 552 INTEGER :: a, atom_a, atom_b, b, c, d, e, handle, & 553 i, iend, igrp, ilist, ipair, istart, & 554 k, nscale 555 INTEGER, DIMENSION(:, :), POINTER :: list 556 LOGICAL :: do_efield0, do_efield1, do_efield2, & 557 force_eval 558 REAL(KIND=dp) :: alpha, ch_i, ch_j, ef0_i, ef0_j, eloc, fac, fac_ij, ir, irab2, ptens11, & 559 ptens12, ptens13, ptens21, ptens22, ptens23, ptens31, ptens32, ptens33, r, rab2, tij, & 560 tmp, tmp1, tmp11, tmp12, tmp13, tmp2, tmp21, tmp22, tmp23, tmp31, tmp32, tmp33, tmp_ij, & 561 tmp_ji 562 REAL(KIND=dp), DIMENSION(0:5) :: f 563 REAL(KIND=dp), DIMENSION(3) :: dp_i, dp_j, ef1_i, ef1_j, fr, rab, tij_a 564 REAL(KIND=dp), DIMENSION(3, 3) :: ef2_i, ef2_j, qp_i, qp_j, tij_ab 565 REAL(KIND=dp), DIMENSION(3, 3, 3) :: tij_abc 566 REAL(KIND=dp), DIMENSION(3, 3, 3, 3) :: tij_abcd 567 REAL(KIND=dp), DIMENSION(3, 3, 3, 3, 3) :: tij_abcde 568 TYPE(fist_neighbor_type), POINTER :: nonbonded 569 TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair 570 571 CALL timeset(routineN, handle) 572 do_efield0 = do_efield .AND. ASSOCIATED(efield0) 573 do_efield1 = do_efield .AND. ASSOCIATED(efield1) 574 do_efield2 = do_efield .AND. ASSOCIATED(efield2) 575 IF (do_stress) THEN 576 ptens11 = 0.0_dp; ptens12 = 0.0_dp; ptens13 = 0.0_dp 577 ptens21 = 0.0_dp; ptens22 = 0.0_dp; ptens23 = 0.0_dp 578 ptens31 = 0.0_dp; ptens32 = 0.0_dp; ptens33 = 0.0_dp 579 END IF 580 CALL ewald_env_get(ewald_env, alpha=alpha) 581 CALL fist_nonbond_env_get(nonbond_env, nonbonded=nonbonded) 582 583 ! Starting the force loop 584 Lists: DO ilist = 1, nonbonded%nlists 585 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist) 586 nscale = neighbor_kind_pair%nscale 587 IF (nscale == 0) CYCLE 588 list => neighbor_kind_pair%list 589 Kind_Group_Loop: DO igrp = 1, neighbor_kind_pair%ngrp_kind 590 istart = neighbor_kind_pair%grp_kind_start(igrp) 591 IF (istart > nscale) CYCLE 592 iend = MIN(neighbor_kind_pair%grp_kind_end(igrp), nscale) 593 Pairs: DO ipair = istart, iend 594 ! only use pairs that are (partially) excluded for electrostatics 595 fac_ij = -1.0_dp+neighbor_kind_pair%ei_scale(ipair) 596 IF (fac_ij >= 0) CYCLE 597 598 atom_a = list(1, ipair) 599 atom_b = list(2, ipair) 600 601 rab = particle_set(atom_b)%r-particle_set(atom_a)%r 602 rab = pbc(rab, cell) 603 rab2 = rab(1)**2+rab(2)**2+rab(3)**2 604$: ewalds_multipole_sr_macro(mode="SCREENED_COULOMB_ERF", store_energy=True, store_forces=True) 605 END DO Pairs 606 END DO Kind_Group_Loop 607 END DO Lists 608 IF (do_stress) THEN 609 pv(1, 1) = pv(1, 1)+ptens11 610 pv(1, 2) = pv(1, 2)+(ptens12+ptens21)*0.5_dp 611 pv(1, 3) = pv(1, 3)+(ptens13+ptens31)*0.5_dp 612 pv(2, 1) = pv(1, 2) 613 pv(2, 2) = pv(2, 2)+ptens22 614 pv(2, 3) = pv(2, 3)+(ptens23+ptens32)*0.5_dp 615 pv(3, 1) = pv(1, 3) 616 pv(3, 2) = pv(2, 3) 617 pv(3, 3) = pv(3, 3)+ptens33 618 END IF 619 620 CALL timestop(handle) 621 END SUBROUTINE ewald_multipole_bonded 622 623! ************************************************************************************************** 624!> \brief computes the potential and the force for a lattice sum of multipoles 625!> up to quadrupole - Long Range (Reciprocal Space) Term 626!> \param ewald_env ... 627!> \param ewald_pw ... 628!> \param cell ... 629!> \param particle_set ... 630!> \param local_particles ... 631!> \param energy ... 632!> \param task ... 633!> \param do_forces ... 634!> \param do_efield ... 635!> \param do_stress ... 636!> \param charges ... 637!> \param dipoles ... 638!> \param quadrupoles ... 639!> \param forces ... 640!> \param pv ... 641!> \param efield0 ... 642!> \param efield1 ... 643!> \param efield2 ... 644!> \author Teodoro Laino [tlaino] - 12.2007 - University of Zurich 645! ************************************************************************************************** 646 SUBROUTINE ewald_multipole_LR(ewald_env, ewald_pw, cell, particle_set, & 647 local_particles, energy, task, do_forces, do_efield, do_stress, & 648 charges, dipoles, quadrupoles, forces, pv, efield0, efield1, efield2) 649 TYPE(ewald_environment_type), POINTER :: ewald_env 650 TYPE(ewald_pw_type), POINTER :: ewald_pw 651 TYPE(cell_type), POINTER :: cell 652 TYPE(particle_type), POINTER :: particle_set(:) 653 TYPE(distribution_1d_type), POINTER :: local_particles 654 REAL(KIND=dp), INTENT(INOUT) :: energy 655 LOGICAL, DIMENSION(3, 3), INTENT(IN) :: task 656 LOGICAL, INTENT(IN) :: do_forces, do_efield, do_stress 657 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: charges 658 REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: dipoles 659 REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, & 660 POINTER :: quadrupoles 661 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), & 662 OPTIONAL :: forces, pv 663 REAL(KIND=dp), DIMENSION(:), POINTER :: efield0 664 REAL(KIND=dp), DIMENSION(:, :), POINTER :: efield1, efield2 665 666 CHARACTER(len=*), PARAMETER :: routineN = 'ewald_multipole_LR', & 667 routineP = moduleN//':'//routineN 668 669 COMPLEX(KIND=dp) :: atm_factor, atm_factor_st(3), cnjg_fac, & 670 fac, summe_tmp 671 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: summe_ef 672 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: summe_st 673 INTEGER :: gpt, group, handle, iparticle, iparticle_kind, iparticle_local, lp, mp, nnodes, & 674 node, np, nparticle_kind, nparticle_local 675 INTEGER, DIMENSION(:, :), POINTER :: bds 676 LOGICAL :: do_efield0, do_efield1, do_efield2 677 REAL(KIND=dp) :: alpha, denom, dipole_t(3), f0, factor, & 678 four_alpha_sq, gauss, pref, q_t, tmp, & 679 trq_t 680 REAL(KIND=dp), DIMENSION(3) :: tmp_v, vec 681 REAL(KIND=dp), DIMENSION(3, 3) :: pv_tmp 682 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: rho0 683 TYPE(dg_rho0_type), POINTER :: dg_rho0 684 TYPE(dg_type), POINTER :: dg 685 TYPE(pw_grid_type), POINTER :: pw_grid 686 TYPE(pw_pool_type), POINTER :: pw_pool 687 TYPE(structure_factor_type) :: exp_igr 688 689 CALL timeset(routineN, handle) 690 do_efield0 = do_efield .AND. ASSOCIATED(efield0) 691 do_efield1 = do_efield .AND. ASSOCIATED(efield1) 692 do_efield2 = do_efield .AND. ASSOCIATED(efield2) 693 694 ! Gathering data from the ewald environment 695 CALL ewald_env_get(ewald_env, alpha=alpha, group=group) 696 CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, dg=dg) 697 CALL dg_get(dg, dg_rho0=dg_rho0) 698 rho0 => dg_rho0%density%pw%cr3d 699 pw_grid => pw_pool%pw_grid 700 bds => pw_grid%bounds 701 702 ! Allocation of working arrays 703 nparticle_kind = SIZE(local_particles%n_el) 704 nnodes = 0 705 DO iparticle_kind = 1, nparticle_kind 706 nnodes = nnodes+local_particles%n_el(iparticle_kind) 707 ENDDO 708 CALL structure_factor_allocate(pw_grid%bounds, nnodes, exp_igr) 709 710 ALLOCATE (summe_ef(1:pw_grid%ngpts_cut)) 711 summe_ef = CMPLX(0.0_dp, 0.0_dp, KIND=dp) 712 ! Stress Tensor 713 IF (do_stress) THEN 714 pv_tmp = 0.0_dp 715 ALLOCATE (summe_st(3, 1:pw_grid%ngpts_cut)) 716 summe_st = CMPLX(0.0_dp, 0.0_dp, KIND=dp) 717 END IF 718 719 ! Defining four_alpha_sq 720 four_alpha_sq = 4.0_dp*alpha**2 721 dipole_t = 0.0_dp 722 q_t = 0.0_dp 723 trq_t = 0.0_dp 724 ! Zero node count 725 node = 0 726 DO iparticle_kind = 1, nparticle_kind 727 nparticle_local = local_particles%n_el(iparticle_kind) 728 DO iparticle_local = 1, nparticle_local 729 node = node+1 730 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 731 CALL matvec_3x3(vec, cell%h_inv, particle_set(iparticle)%r) 732 CALL structure_factor_evaluate(vec, exp_igr%lb, & 733 exp_igr%ex(:, node), exp_igr%ey(:, node), exp_igr%ez(:, node)) 734 735 ! Computing the total charge, dipole and quadrupole trace (if any) 736 IF (ANY(task(1, :))) THEN 737 q_t = q_t+charges(iparticle) 738 END IF 739 IF (ANY(task(2, :))) THEN 740 dipole_t = dipole_t+dipoles(:, iparticle) 741 END IF 742 IF (ANY(task(3, :))) THEN 743 trq_t = trq_t+quadrupoles(1, 1, iparticle)+ & 744 quadrupoles(2, 2, iparticle)+ & 745 quadrupoles(3, 3, iparticle) 746 END IF 747 END DO 748 END DO 749 750 ! Looping over the positive g-vectors 751 DO gpt = 1, pw_grid%ngpts_cut_local 752 lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt)) 753 mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt)) 754 np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt)) 755 756 lp = lp+bds(1, 1) 757 mp = mp+bds(1, 2) 758 np = np+bds(1, 3) 759 760 ! Initializing sum to be used in the energy and force 761 node = 0 762 DO iparticle_kind = 1, nparticle_kind 763 nparticle_local = local_particles%n_el(iparticle_kind) 764 DO iparticle_local = 1, nparticle_local 765 node = node+1 766 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 767 ! Density for energy and forces 768 CALL get_atom_factor(atm_factor, pw_grid, gpt, iparticle, task, charges, & 769 dipoles, quadrupoles) 770 summe_tmp = exp_igr%ex(lp, node)*exp_igr%ey(mp, node)*exp_igr%ez(np, node) 771 summe_ef(gpt) = summe_ef(gpt)+atm_factor*summe_tmp 772 773 ! Precompute pseudo-density for stress tensor calculation 774 IF (do_stress) THEN 775 CALL get_atom_factor_stress(atm_factor_st, pw_grid, gpt, iparticle, task, & 776 dipoles, quadrupoles) 777 summe_st(1:3, gpt) = summe_st(1:3, gpt)+atm_factor_st(1:3)*summe_tmp 778 END IF 779 END DO 780 END DO 781 END DO 782 CALL mp_sum(q_t, group) 783 CALL mp_sum(dipole_t, group) 784 CALL mp_sum(trq_t, group) 785 CALL mp_sum(summe_ef, group) 786 IF (do_stress) CALL mp_sum(summe_st, group) 787 788 ! Looping over the positive g-vectors 789 DO gpt = 1, pw_grid%ngpts_cut_local 790 ! computing the potential energy 791 lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt)) 792 mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt)) 793 np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt)) 794 795 lp = lp+bds(1, 1) 796 mp = mp+bds(1, 2) 797 np = np+bds(1, 3) 798 799 IF (pw_grid%gsq(gpt) == 0.0_dp) THEN 800 ! G=0 vector for dipole-dipole and charge-quadrupole 801 energy = energy+(1.0_dp/6.0_dp)*DOT_PRODUCT(dipole_t, dipole_t) & 802 -(1.0_dp/9.0_dp)*q_t*trq_t 803 ! Stress tensor 804 IF (do_stress) THEN 805 pv_tmp(1, 1) = pv_tmp(1, 1)+(1.0_dp/6.0_dp)*DOT_PRODUCT(dipole_t, dipole_t) 806 pv_tmp(2, 2) = pv_tmp(2, 2)+(1.0_dp/6.0_dp)*DOT_PRODUCT(dipole_t, dipole_t) 807 pv_tmp(3, 3) = pv_tmp(3, 3)+(1.0_dp/6.0_dp)*DOT_PRODUCT(dipole_t, dipole_t) 808 END IF 809 ! Corrections for G=0 to potential, field and field gradient 810 IF (do_efield .AND. (debug_e_field_en .OR. (.NOT. debug_this_module))) THEN 811 ! This term is important and may give problems if one is debugging 812 ! VS finite differences since it comes from a residual integral in 813 ! the complex plane (cannot be reproduced with finite differences) 814 node = 0 815 DO iparticle_kind = 1, nparticle_kind 816 nparticle_local = local_particles%n_el(iparticle_kind) 817 DO iparticle_local = 1, nparticle_local 818 node = node+1 819 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 820 821 ! Potential 822 IF (do_efield0) THEN 823 efield0(iparticle) = efield0(iparticle) 824 END IF 825 ! Electrostatic field 826 IF (do_efield1) THEN 827 efield1(1:3, iparticle) = efield1(1:3, iparticle)-(1.0_dp/6.0_dp)*dipole_t(1:3) 828 END IF 829 ! Electrostatic field gradients 830 IF (do_efield2) THEN 831 efield2(1, iparticle) = efield2(1, iparticle)-(1.0_dp/(18.0_dp))*q_t 832 efield2(5, iparticle) = efield2(5, iparticle)-(1.0_dp/(18.0_dp))*q_t 833 efield2(9, iparticle) = efield2(9, iparticle)-(1.0_dp/(18.0_dp))*q_t 834 END IF 835 END DO 836 END DO 837 END IF 838 CYCLE 839 END IF 840 gauss = (rho0(lp, mp, np)*pw_grid%vol)**2/pw_grid%gsq(gpt) 841 factor = gauss*REAL(summe_ef(gpt)*CONJG(summe_ef(gpt)), KIND=dp) 842 energy = energy+factor 843 844 IF (do_forces .OR. do_efield) THEN 845 node = 0 846 DO iparticle_kind = 1, nparticle_kind 847 nparticle_local = local_particles%n_el(iparticle_kind) 848 DO iparticle_local = 1, nparticle_local 849 node = node+1 850 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 851 fac = exp_igr%ex(lp, node)*exp_igr%ey(mp, node)*exp_igr%ez(np, node) 852 cnjg_fac = CONJG(fac) 853 854 ! Forces 855 IF (do_forces) THEN 856 CALL get_atom_factor(atm_factor, pw_grid, gpt, iparticle, task, charges, & 857 dipoles, quadrupoles) 858 859 tmp = gauss*AIMAG(summe_ef(gpt)*(cnjg_fac*CONJG(atm_factor))) 860 forces(1, node) = forces(1, node)+tmp*pw_grid%g(1, gpt) 861 forces(2, node) = forces(2, node)+tmp*pw_grid%g(2, gpt) 862 forces(3, node) = forces(3, node)+tmp*pw_grid%g(3, gpt) 863 END IF 864 865 ! Electric field 866 IF (do_efield) THEN 867 ! Potential 868 IF (do_efield0) THEN 869 efield0(iparticle) = efield0(iparticle)+gauss*REAL(fac*CONJG(summe_ef(gpt)), KIND=dp) 870 END IF 871 ! Electric field 872 IF (do_efield1) THEN 873 tmp = AIMAG(fac*CONJG(summe_ef(gpt)))*gauss 874 efield1(1, iparticle) = efield1(1, iparticle)-tmp*pw_grid%g(1, gpt) 875 efield1(2, iparticle) = efield1(2, iparticle)-tmp*pw_grid%g(2, gpt) 876 efield1(3, iparticle) = efield1(3, iparticle)-tmp*pw_grid%g(3, gpt) 877 END IF 878 ! Electric field gradient 879 IF (do_efield2) THEN 880 tmp_v(1) = REAL(fac*CONJG(summe_ef(gpt)), KIND=dp)*pw_grid%g(1, gpt)*gauss 881 tmp_v(2) = REAL(fac*CONJG(summe_ef(gpt)), KIND=dp)*pw_grid%g(2, gpt)*gauss 882 tmp_v(3) = REAL(fac*CONJG(summe_ef(gpt)), KIND=dp)*pw_grid%g(3, gpt)*gauss 883 884 efield2(1, iparticle) = efield2(1, iparticle)+tmp_v(1)*pw_grid%g(1, gpt) 885 efield2(2, iparticle) = efield2(2, iparticle)+tmp_v(1)*pw_grid%g(2, gpt) 886 efield2(3, iparticle) = efield2(3, iparticle)+tmp_v(1)*pw_grid%g(3, gpt) 887 efield2(4, iparticle) = efield2(4, iparticle)+tmp_v(2)*pw_grid%g(1, gpt) 888 efield2(5, iparticle) = efield2(5, iparticle)+tmp_v(2)*pw_grid%g(2, gpt) 889 efield2(6, iparticle) = efield2(6, iparticle)+tmp_v(2)*pw_grid%g(3, gpt) 890 efield2(7, iparticle) = efield2(7, iparticle)+tmp_v(3)*pw_grid%g(1, gpt) 891 efield2(8, iparticle) = efield2(8, iparticle)+tmp_v(3)*pw_grid%g(2, gpt) 892 efield2(9, iparticle) = efield2(9, iparticle)+tmp_v(3)*pw_grid%g(3, gpt) 893 END IF 894 END IF 895 END DO 896 END DO 897 END IF 898 899 ! Compute the virial P*V 900 IF (do_stress) THEN 901 ! The Stress Tensor can be decomposed in two main components. 902 ! The first one is just a normal ewald component for reciprocal space 903 denom = 1.0_dp/four_alpha_sq+1.0_dp/pw_grid%gsq(gpt) 904 pv_tmp(1, 1) = pv_tmp(1, 1)+factor*(1.0_dp-2.0_dp*pw_grid%g(1, gpt)*pw_grid%g(1, gpt)*denom) 905 pv_tmp(1, 2) = pv_tmp(1, 2)-factor*(2.0_dp*pw_grid%g(1, gpt)*pw_grid%g(2, gpt)*denom) 906 pv_tmp(1, 3) = pv_tmp(1, 3)-factor*(2.0_dp*pw_grid%g(1, gpt)*pw_grid%g(3, gpt)*denom) 907 pv_tmp(2, 1) = pv_tmp(2, 1)-factor*(2.0_dp*pw_grid%g(2, gpt)*pw_grid%g(1, gpt)*denom) 908 pv_tmp(2, 2) = pv_tmp(2, 2)+factor*(1.0_dp-2.0_dp*pw_grid%g(2, gpt)*pw_grid%g(2, gpt)*denom) 909 pv_tmp(2, 3) = pv_tmp(2, 3)-factor*(2.0_dp*pw_grid%g(2, gpt)*pw_grid%g(3, gpt)*denom) 910 pv_tmp(3, 1) = pv_tmp(3, 1)-factor*(2.0_dp*pw_grid%g(3, gpt)*pw_grid%g(1, gpt)*denom) 911 pv_tmp(3, 2) = pv_tmp(3, 2)-factor*(2.0_dp*pw_grid%g(3, gpt)*pw_grid%g(2, gpt)*denom) 912 pv_tmp(3, 3) = pv_tmp(3, 3)+factor*(1.0_dp-2.0_dp*pw_grid%g(3, gpt)*pw_grid%g(3, gpt)*denom) 913 ! The second one can be written in the following way 914 f0 = 2.0_dp*gauss 915 pv_tmp(1, 1) = pv_tmp(1, 1)+f0*pw_grid%g(1, gpt)*REAL(summe_st(1, gpt)*CONJG(summe_ef(gpt)), KIND=dp) 916 pv_tmp(1, 2) = pv_tmp(1, 2)+f0*pw_grid%g(1, gpt)*REAL(summe_st(2, gpt)*CONJG(summe_ef(gpt)), KIND=dp) 917 pv_tmp(1, 3) = pv_tmp(1, 3)+f0*pw_grid%g(1, gpt)*REAL(summe_st(3, gpt)*CONJG(summe_ef(gpt)), KIND=dp) 918 pv_tmp(2, 1) = pv_tmp(2, 1)+f0*pw_grid%g(2, gpt)*REAL(summe_st(1, gpt)*CONJG(summe_ef(gpt)), KIND=dp) 919 pv_tmp(2, 2) = pv_tmp(2, 2)+f0*pw_grid%g(2, gpt)*REAL(summe_st(2, gpt)*CONJG(summe_ef(gpt)), KIND=dp) 920 pv_tmp(2, 3) = pv_tmp(2, 3)+f0*pw_grid%g(2, gpt)*REAL(summe_st(3, gpt)*CONJG(summe_ef(gpt)), KIND=dp) 921 pv_tmp(3, 1) = pv_tmp(3, 1)+f0*pw_grid%g(3, gpt)*REAL(summe_st(1, gpt)*CONJG(summe_ef(gpt)), KIND=dp) 922 pv_tmp(3, 2) = pv_tmp(3, 2)+f0*pw_grid%g(3, gpt)*REAL(summe_st(2, gpt)*CONJG(summe_ef(gpt)), KIND=dp) 923 pv_tmp(3, 3) = pv_tmp(3, 3)+f0*pw_grid%g(3, gpt)*REAL(summe_st(3, gpt)*CONJG(summe_ef(gpt)), KIND=dp) 924 END IF 925 END DO 926 pref = fourpi/pw_grid%vol 927 energy = energy*pref 928 929 CALL structure_factor_deallocate(exp_igr) 930 DEALLOCATE (summe_ef) 931 IF (do_stress) THEN 932 pv_tmp = pv_tmp*pref 933 ! Symmetrize the tensor 934 pv(1, 1) = pv(1, 1)+pv_tmp(1, 1) 935 pv(1, 2) = pv(1, 2)+(pv_tmp(1, 2)+pv_tmp(2, 1))*0.5_dp 936 pv(1, 3) = pv(1, 3)+(pv_tmp(1, 3)+pv_tmp(3, 1))*0.5_dp 937 pv(2, 1) = pv(1, 2) 938 pv(2, 2) = pv(2, 2)+pv_tmp(2, 2) 939 pv(2, 3) = pv(2, 3)+(pv_tmp(2, 3)+pv_tmp(3, 2))*0.5_dp 940 pv(3, 1) = pv(1, 3) 941 pv(3, 2) = pv(2, 3) 942 pv(3, 3) = pv(3, 3)+pv_tmp(3, 3) 943 DEALLOCATE (summe_st) 944 END IF 945 IF (do_forces) THEN 946 forces = 2.0_dp*forces*pref 947 END IF 948 IF (do_efield0) THEN 949 efield0 = 2.0_dp*efield0*pref 950 END IF 951 IF (do_efield1) THEN 952 efield1 = 2.0_dp*efield1*pref 953 END IF 954 IF (do_efield2) THEN 955 efield2 = 2.0_dp*efield2*pref 956 END IF 957 CALL timestop(handle) 958 959 END SUBROUTINE ewald_multipole_LR 960 961! ************************************************************************************************** 962!> \brief Computes the atom factor including charge, dipole and quadrupole 963!> \param atm_factor ... 964!> \param pw_grid ... 965!> \param gpt ... 966!> \param iparticle ... 967!> \param task ... 968!> \param charges ... 969!> \param dipoles ... 970!> \param quadrupoles ... 971!> \par History 972!> none 973!> \author Teodoro Laino [tlaino] - 12.2007 - University of Zurich 974! ************************************************************************************************** 975 SUBROUTINE get_atom_factor(atm_factor, pw_grid, gpt, iparticle, task, charges, & 976 dipoles, quadrupoles) 977 COMPLEX(KIND=dp), INTENT(OUT) :: atm_factor 978 TYPE(pw_grid_type), POINTER :: pw_grid 979 INTEGER, INTENT(IN) :: gpt 980 INTEGER :: iparticle 981 LOGICAL, DIMENSION(3, 3), INTENT(IN) :: task 982 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: charges 983 REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: dipoles 984 REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, & 985 POINTER :: quadrupoles 986 987 CHARACTER(len=*), PARAMETER :: routineN = 'get_atom_factor', & 988 routineP = moduleN//':'//routineN 989 990 COMPLEX(KIND=dp) :: tmp 991 INTEGER :: i, j 992 993 atm_factor = CMPLX(0.0_dp, 0.0_dp, KIND=dp) 994 IF (task(1, 1)) THEN 995 ! Charge 996 atm_factor = atm_factor+charges(iparticle) 997 END IF 998 IF (task(2, 2)) THEN 999 ! Dipole 1000 tmp = CMPLX(0.0_dp, 0.0_dp, KIND=dp) 1001 DO i = 1, 3 1002 tmp = tmp+dipoles(i, iparticle)*pw_grid%g(i, gpt) 1003 END DO 1004 atm_factor = atm_factor+tmp*CMPLX(0.0_dp, -1.0_dp, KIND=dp) 1005 END IF 1006 IF (task(3, 3)) THEN 1007 ! Quadrupole 1008 tmp = CMPLX(0.0_dp, 0.0_dp, KIND=dp) 1009 DO i = 1, 3 1010 DO j = 1, 3 1011 tmp = tmp+quadrupoles(j, i, iparticle)*pw_grid%g(j, gpt)*pw_grid%g(i, gpt) 1012 END DO 1013 END DO 1014 atm_factor = atm_factor-1.0_dp/3.0_dp*tmp 1015 END IF 1016 1017 END SUBROUTINE get_atom_factor 1018 1019! ************************************************************************************************** 1020!> \brief Computes the atom factor including charge, dipole and quadrupole 1021!> \param atm_factor ... 1022!> \param pw_grid ... 1023!> \param gpt ... 1024!> \param iparticle ... 1025!> \param task ... 1026!> \param dipoles ... 1027!> \param quadrupoles ... 1028!> \par History 1029!> none 1030!> \author Teodoro Laino [tlaino] - 12.2007 - University of Zurich 1031! ************************************************************************************************** 1032 SUBROUTINE get_atom_factor_stress(atm_factor, pw_grid, gpt, iparticle, task, & 1033 dipoles, quadrupoles) 1034 COMPLEX(KIND=dp), INTENT(OUT) :: atm_factor(3) 1035 TYPE(pw_grid_type), POINTER :: pw_grid 1036 INTEGER, INTENT(IN) :: gpt 1037 INTEGER :: iparticle 1038 LOGICAL, DIMENSION(3, 3), INTENT(IN) :: task 1039 REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: dipoles 1040 REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, & 1041 POINTER :: quadrupoles 1042 1043 CHARACTER(len=*), PARAMETER :: routineN = 'get_atom_factor_stress', & 1044 routineP = moduleN//':'//routineN 1045 1046 INTEGER :: i 1047 1048 atm_factor = CMPLX(0.0_dp, 0.0_dp, KIND=dp) 1049 IF (ANY(task(2, :))) THEN 1050 ! Dipole 1051 atm_factor = dipoles(:, iparticle)*CMPLX(0.0_dp, -1.0_dp, KIND=dp) 1052 END IF 1053 IF (ANY(task(3, :))) THEN 1054 ! Quadrupole 1055 DO i = 1, 3 1056 atm_factor(1) = atm_factor(1)-1.0_dp/3.0_dp* & 1057 (quadrupoles(1, i, iparticle)*pw_grid%g(i, gpt)+ & 1058 quadrupoles(i, 1, iparticle)*pw_grid%g(i, gpt)) 1059 atm_factor(2) = atm_factor(2)-1.0_dp/3.0_dp* & 1060 (quadrupoles(2, i, iparticle)*pw_grid%g(i, gpt)+ & 1061 quadrupoles(i, 2, iparticle)*pw_grid%g(i, gpt)) 1062 atm_factor(3) = atm_factor(3)-1.0_dp/3.0_dp* & 1063 (quadrupoles(3, i, iparticle)*pw_grid%g(i, gpt)+ & 1064 quadrupoles(i, 3, iparticle)*pw_grid%g(i, gpt)) 1065 END DO 1066 END IF 1067 1068 END SUBROUTINE get_atom_factor_stress 1069 1070! ************************************************************************************************** 1071!> \brief Computes the self interaction from g-space and the neutralizing background 1072!> when using multipoles 1073!> \param ewald_env ... 1074!> \param cell ... 1075!> \param local_particles ... 1076!> \param e_self ... 1077!> \param e_neut ... 1078!> \param task ... 1079!> \param do_efield ... 1080!> \param radii ... 1081!> \param charges ... 1082!> \param dipoles ... 1083!> \param quadrupoles ... 1084!> \param efield0 ... 1085!> \param efield1 ... 1086!> \param efield2 ... 1087!> \author Teodoro Laino [tlaino] - University of Zurich - 12.2007 1088! ************************************************************************************************** 1089 SUBROUTINE ewald_multipole_self(ewald_env, cell, local_particles, e_self, & 1090 e_neut, task, do_efield, radii, charges, dipoles, quadrupoles, efield0, & 1091 efield1, efield2) 1092 TYPE(ewald_environment_type), POINTER :: ewald_env 1093 TYPE(cell_type), INTENT(IN) :: cell 1094 TYPE(distribution_1d_type), POINTER :: local_particles 1095 REAL(KIND=dp), INTENT(OUT) :: e_self, e_neut 1096 LOGICAL, DIMENSION(3, 3), INTENT(IN) :: task 1097 LOGICAL, INTENT(IN) :: do_efield 1098 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: radii, charges 1099 REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: dipoles 1100 REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, & 1101 POINTER :: quadrupoles 1102 REAL(KIND=dp), DIMENSION(:), POINTER :: efield0 1103 REAL(KIND=dp), DIMENSION(:, :), POINTER :: efield1, efield2 1104 1105 CHARACTER(len=*), PARAMETER :: routineN = 'ewald_multipole_self', & 1106 routineP = moduleN//':'//routineN 1107 REAL(KIND=dp), PARAMETER :: f23 = 2.0_dp/3.0_dp, & 1108 f415 = 4.0_dp/15.0_dp 1109 1110 INTEGER :: ewald_type, group, i, iparticle, & 1111 iparticle_kind, iparticle_local, j, & 1112 nparticle_local 1113 LOGICAL :: do_efield0, do_efield1, do_efield2, & 1114 lradii 1115 REAL(KIND=dp) :: alpha, ch_qu_self, ch_qu_self_tmp, & 1116 dipole_self, fac1, fac2, fac3, fac4, & 1117 q, q_neutg, q_self, q_sum, qu_qu_self, & 1118 radius 1119 1120 CALL ewald_env_get(ewald_env, ewald_type=ewald_type, alpha=alpha, & 1121 group=group) 1122 1123 do_efield0 = do_efield .AND. ASSOCIATED(efield0) 1124 do_efield1 = do_efield .AND. ASSOCIATED(efield1) 1125 do_efield2 = do_efield .AND. ASSOCIATED(efield2) 1126 q_self = 0.0_dp 1127 q_sum = 0.0_dp 1128 dipole_self = 0.0_dp 1129 ch_qu_self = 0.0_dp 1130 qu_qu_self = 0.0_dp 1131 fac1 = 2.0_dp*alpha*oorootpi 1132 fac2 = 6.0_dp*(f23**2)*(alpha**3)*oorootpi 1133 fac3 = (2.0_dp*oorootpi)*f23*alpha**3 1134 fac4 = (4.0_dp*oorootpi)*f415*alpha**5 1135 lradii = PRESENT(radii) 1136 radius = 0.0_dp 1137 q_neutg = 0.0_dp 1138 DO iparticle_kind = 1, SIZE(local_particles%n_el) 1139 nparticle_local = local_particles%n_el(iparticle_kind) 1140 DO iparticle_local = 1, nparticle_local 1141 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 1142 IF (ANY(task(1, :))) THEN 1143 ! Charge - Charge 1144 q = charges(iparticle) 1145 IF (lradii) radius = radii(iparticle) 1146 IF (radius > 0) THEN 1147 q_neutg = q_neutg+2.0_dp*q*radius**2 1148 END IF 1149 q_self = q_self+q*q 1150 q_sum = q_sum+q 1151 ! Potential 1152 IF (do_efield0) THEN 1153 efield0(iparticle) = efield0(iparticle)-q*fac1 1154 END IF 1155 1156 IF (task(1, 3)) THEN 1157 ! Charge - Quadrupole 1158 ch_qu_self_tmp = 0.0_dp 1159 DO i = 1, 3 1160 ch_qu_self_tmp = ch_qu_self_tmp+quadrupoles(i, i, iparticle)*q 1161 END DO 1162 ch_qu_self = ch_qu_self+ch_qu_self_tmp 1163 ! Electric Field Gradient 1164 IF (do_efield2) THEN 1165 efield2(1, iparticle) = efield2(1, iparticle)+fac2*q 1166 efield2(5, iparticle) = efield2(5, iparticle)+fac2*q 1167 efield2(9, iparticle) = efield2(9, iparticle)+fac2*q 1168 END IF 1169 END IF 1170 END IF 1171 IF (ANY(task(2, :))) THEN 1172 ! Dipole - Dipole 1173 DO i = 1, 3 1174 dipole_self = dipole_self+dipoles(i, iparticle)**2 1175 END DO 1176 ! Electric Field 1177 IF (do_efield1) THEN 1178 efield1(1, iparticle) = efield1(1, iparticle)+fac3*dipoles(1, iparticle) 1179 efield1(2, iparticle) = efield1(2, iparticle)+fac3*dipoles(2, iparticle) 1180 efield1(3, iparticle) = efield1(3, iparticle)+fac3*dipoles(3, iparticle) 1181 END IF 1182 END IF 1183 IF (ANY(task(3, :))) THEN 1184 ! Quadrupole - Quadrupole 1185 DO i = 1, 3 1186 DO j = 1, 3 1187 qu_qu_self = qu_qu_self+quadrupoles(j, i, iparticle)**2 1188 END DO 1189 END DO 1190 ! Electric Field Gradient 1191 IF (do_efield2) THEN 1192 efield2(1, iparticle) = efield2(1, iparticle)+fac4*quadrupoles(1, 1, iparticle) 1193 efield2(2, iparticle) = efield2(2, iparticle)+fac4*quadrupoles(2, 1, iparticle) 1194 efield2(3, iparticle) = efield2(3, iparticle)+fac4*quadrupoles(3, 1, iparticle) 1195 efield2(4, iparticle) = efield2(4, iparticle)+fac4*quadrupoles(1, 2, iparticle) 1196 efield2(5, iparticle) = efield2(5, iparticle)+fac4*quadrupoles(2, 2, iparticle) 1197 efield2(6, iparticle) = efield2(6, iparticle)+fac4*quadrupoles(3, 2, iparticle) 1198 efield2(7, iparticle) = efield2(7, iparticle)+fac4*quadrupoles(1, 3, iparticle) 1199 efield2(8, iparticle) = efield2(8, iparticle)+fac4*quadrupoles(2, 3, iparticle) 1200 efield2(9, iparticle) = efield2(9, iparticle)+fac4*quadrupoles(3, 3, iparticle) 1201 END IF 1202 END IF 1203 END DO 1204 END DO 1205 1206 CALL mp_sum(q_neutg, group) 1207 CALL mp_sum(q_self, group) 1208 CALL mp_sum(q_sum, group) 1209 CALL mp_sum(dipole_self, group) 1210 CALL mp_sum(ch_qu_self, group) 1211 CALL mp_sum(qu_qu_self, group) 1212 1213 e_self = -(q_self+f23*(dipole_self-f23*ch_qu_self+f415*qu_qu_self*alpha**2)*alpha**2)*alpha*oorootpi 1214 fac1 = pi/(2.0_dp*cell%deth) 1215 e_neut = -q_sum*fac1*(q_sum/alpha**2-q_neutg) 1216 1217 ! Correcting Potential for the neutralizing background charge 1218 DO iparticle_kind = 1, SIZE(local_particles%n_el) 1219 nparticle_local = local_particles%n_el(iparticle_kind) 1220 DO iparticle_local = 1, nparticle_local 1221 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 1222 IF (ANY(task(1, :))) THEN 1223 ! Potential energy 1224 IF (do_efield0) THEN 1225 efield0(iparticle) = efield0(iparticle)-q_sum*2.0_dp*fac1/alpha**2 1226 IF (lradii) radius = radii(iparticle) 1227 IF (radius > 0) THEN 1228 q = charges(iparticle) 1229 efield0(iparticle) = efield0(iparticle)+fac1*radius**2*(q_sum+q) 1230 END IF 1231 END IF 1232 END IF 1233 END DO 1234 END DO 1235 END SUBROUTINE ewald_multipole_self 1236 1237! ************************************************************************************************** 1238!> \brief ... 1239!> \param iw ... 1240!> \param e_gspace ... 1241!> \param e_rspace ... 1242!> \param e_bonded ... 1243!> \param e_self ... 1244!> \param e_neut ... 1245!> \author Teodoro Laino [tlaino] - University of Zurich - 12.2007 1246! ************************************************************************************************** 1247 SUBROUTINE ewald_multipole_print(iw, e_gspace, e_rspace, e_bonded, e_self, e_neut) 1248 1249 INTEGER, INTENT(IN) :: iw 1250 REAL(KIND=dp), INTENT(IN) :: e_gspace, e_rspace, e_bonded, e_self, & 1251 e_neut 1252 1253 CHARACTER(len=*), PARAMETER :: routineN = 'ewald_multipole_print', & 1254 routineP = moduleN//':'//routineN 1255 1256 IF (iw > 0) THEN 1257 WRITE (iw, '( A, A )') ' *********************************', & 1258 '**********************************************' 1259 WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' INITIAL GSPACE ENERGY', & 1260 '[hartree]', '= ', e_gspace 1261 WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' INITIAL RSPACE ENERGY', & 1262 '[hartree]', '= ', e_rspace 1263 WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' BONDED CORRECTION', & 1264 '[hartree]', '= ', e_bonded 1265 WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' SELF ENERGY CORRECTION', & 1266 '[hartree]', '= ', e_self 1267 WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' NEUTRALIZ. BCKGR. ENERGY', & 1268 '[hartree]', '= ', e_neut 1269 WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' TOTAL ELECTROSTATIC EN.', & 1270 '[hartree]', '= ', e_rspace+e_bonded+e_gspace+e_self+e_neut 1271 WRITE (iw, '( A, A )') ' *********************************', & 1272 '**********************************************' 1273 END IF 1274 END SUBROUTINE ewald_multipole_print 1275 1276! ************************************************************************************************** 1277!> \brief Debug routines for multipoles 1278!> \param ewald_env ... 1279!> \param ewald_pw ... 1280!> \param nonbond_env ... 1281!> \param cell ... 1282!> \param particle_set ... 1283!> \param local_particles ... 1284!> \param iw ... 1285!> \param debug_r_space ... 1286!> \date 05.2008 1287!> \author Teodoro Laino [tlaino] - University of Zurich - 05.2008 1288! ************************************************************************************************** 1289 SUBROUTINE debug_ewald_multipoles(ewald_env, ewald_pw, nonbond_env, cell, & 1290 particle_set, local_particles, iw, debug_r_space) 1291 TYPE charge_mono_type 1292 REAL(KIND=dp), DIMENSION(:), & 1293 POINTER :: charge 1294 REAL(KIND=dp), DIMENSION(:, :), & 1295 POINTER :: pos 1296 END TYPE charge_mono_type 1297 TYPE multi_charge_type 1298 TYPE(charge_mono_type), DIMENSION(:), & 1299 POINTER :: charge_typ 1300 END TYPE multi_charge_type 1301 TYPE(ewald_environment_type), POINTER :: ewald_env 1302 TYPE(ewald_pw_type), POINTER :: ewald_pw 1303 TYPE(fist_nonbond_env_type), POINTER :: nonbond_env 1304 TYPE(cell_type), POINTER :: cell 1305 TYPE(particle_type), DIMENSION(:), & 1306 POINTER :: particle_set 1307 TYPE(distribution_1d_type), POINTER :: local_particles 1308 INTEGER, INTENT(IN) :: iw 1309 LOGICAL, INTENT(IN) :: debug_r_space 1310 1311 CHARACTER(len=*), PARAMETER :: routineN = 'debug_ewald_multipoles', & 1312 routineP = "ewalds_multipole_debug"//':'//routineN 1313 1314 INTEGER :: nparticles 1315 LOGICAL, DIMENSION(3) :: task 1316 REAL(KIND=dp) :: e_neut, e_self, g_energy, & 1317 r_energy, debug_energy 1318 REAL(KIND=dp), POINTER, DIMENSION(:) :: charges 1319 REAL(KIND=dp), POINTER, & 1320 DIMENSION(:, :) :: dipoles, g_forces, g_pv, & 1321 r_forces, r_pv, e_field1, & 1322 e_field2 1323 REAL(KIND=dp), POINTER, & 1324 DIMENSION(:, :, :) :: quadrupoles 1325 TYPE(rng_stream_type), POINTER :: random_stream 1326 TYPE(multi_charge_type), DIMENSION(:), & 1327 POINTER :: multipoles 1328 1329 NULLIFY (random_stream, multipoles, charges, dipoles, g_forces, g_pv, & 1330 r_forces, r_pv, e_field1, e_field2) 1331 CALL create_rng_stream(random_stream, name="DEBUG_EWALD_MULTIPOLE", & 1332 distribution_type=UNIFORM) 1333 ! check: charge - charge 1334 task = .FALSE. 1335 nparticles = SIZE(particle_set) 1336 1337 ! Allocate charges, dipoles, quadrupoles 1338 ALLOCATE (charges(nparticles)) 1339 ALLOCATE (dipoles(3, nparticles)) 1340 ALLOCATE (quadrupoles(3, 3, nparticles)) 1341 1342 ! Allocate arrays for forces 1343 ALLOCATE (r_forces(3, nparticles)) 1344 ALLOCATE (g_forces(3, nparticles)) 1345 ALLOCATE (e_field1(3, nparticles)) 1346 ALLOCATE (e_field2(3, nparticles)) 1347 ALLOCATE (g_pv(3, 3)) 1348 ALLOCATE (r_pv(3, 3)) 1349 1350 ! Debug CHARGES-CHARGES 1351 task(1) = .TRUE. 1352 charges = 0.0_dp 1353 dipoles = 0.0_dp 1354 quadrupoles = 0.0_dp 1355 r_forces = 0.0_dp 1356 g_forces = 0.0_dp 1357 e_field1 = 0.0_dp 1358 e_field2 = 0.0_dp 1359 g_pv = 0.0_dp 1360 r_pv = 0.0_dp 1361 g_energy = 0.0_dp 1362 r_energy = 0.0_dp 1363 e_neut = 0.0_dp 1364 e_self = 0.0_dp 1365 1366 CALL create_multi_type(multipoles, nparticles, 1, nparticles/2, "CHARGE", echarge=-1.0_dp, & 1367 random_stream=random_stream, charges=charges) 1368 CALL create_multi_type(multipoles, nparticles, nparticles/2+1, nparticles, "CHARGE", echarge=1.0_dp, & 1369 random_stream=random_stream, charges=charges) 1370 CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, & 1371 debug_r_space) 1372 1373 WRITE (iw, *) "DEBUG ENERGY (CHARGE-CHARGE): ", debug_energy 1374 CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, & 1375 particle_set, local_particles, g_energy, r_energy, e_neut, e_self, & 1376 task, do_correction_bonded=.FALSE., do_forces=.TRUE., do_stress=.TRUE., do_efield=.FALSE., & 1377 charges=charges, dipoles=dipoles, quadrupoles=quadrupoles, forces_local=g_forces, & 1378 forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.FALSE.) 1379 CALL release_multi_type(multipoles) 1380 1381 ! Debug CHARGES-DIPOLES 1382 task(1) = .TRUE. 1383 task(2) = .TRUE. 1384 charges = 0.0_dp 1385 dipoles = 0.0_dp 1386 quadrupoles = 0.0_dp 1387 r_forces = 0.0_dp 1388 g_forces = 0.0_dp 1389 e_field1 = 0.0_dp 1390 e_field2 = 0.0_dp 1391 g_pv = 0.0_dp 1392 r_pv = 0.0_dp 1393 g_energy = 0.0_dp 1394 r_energy = 0.0_dp 1395 e_neut = 0.0_dp 1396 e_self = 0.0_dp 1397 1398 CALL create_multi_type(multipoles, nparticles, 1, nparticles/2, "CHARGE", echarge=-1.0_dp, & 1399 random_stream=random_stream, charges=charges) 1400 CALL create_multi_type(multipoles, nparticles, nparticles/2+1, nparticles, "DIPOLE", echarge=0.5_dp, & 1401 random_stream=random_stream, dipoles=dipoles) 1402 WRITE (iw, '("CHARGES",F15.9)') charges 1403 WRITE (iw, '("DIPOLES",3F15.9)') dipoles 1404 CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, & 1405 debug_r_space) 1406 1407 WRITE (iw, *) "DEBUG ENERGY (CHARGE-DIPOLE): ", debug_energy 1408 CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, & 1409 particle_set, local_particles, g_energy, r_energy, e_neut, e_self, & 1410 task, do_correction_bonded=.FALSE., do_forces=.TRUE., do_stress=.TRUE., do_efield=.FALSE., & 1411 charges=charges, dipoles=dipoles, quadrupoles=quadrupoles, forces_local=g_forces, & 1412 forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.FALSE.) 1413 CALL release_multi_type(multipoles) 1414 1415 ! Debug DIPOLES-DIPOLES 1416 task(2) = .TRUE. 1417 charges = 0.0_dp 1418 dipoles = 0.0_dp 1419 quadrupoles = 0.0_dp 1420 r_forces = 0.0_dp 1421 g_forces = 0.0_dp 1422 e_field1 = 0.0_dp 1423 e_field2 = 0.0_dp 1424 g_pv = 0.0_dp 1425 r_pv = 0.0_dp 1426 g_energy = 0.0_dp 1427 r_energy = 0.0_dp 1428 e_neut = 0.0_dp 1429 e_self = 0.0_dp 1430 1431 CALL create_multi_type(multipoles, nparticles, 1, nparticles/2, "DIPOLE", echarge=10000.0_dp, & 1432 random_stream=random_stream, dipoles=dipoles) 1433 CALL create_multi_type(multipoles, nparticles, nparticles/2+1, nparticles, "DIPOLE", echarge=20000._dp, & 1434 random_stream=random_stream, dipoles=dipoles) 1435 WRITE (iw, '("DIPOLES",3F15.9)') dipoles 1436 CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, & 1437 debug_r_space) 1438 1439 WRITE (iw, *) "DEBUG ENERGY (DIPOLE-DIPOLE): ", debug_energy 1440 CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, & 1441 particle_set, local_particles, g_energy, r_energy, e_neut, e_self, & 1442 task, do_correction_bonded=.FALSE., do_forces=.TRUE., do_stress=.TRUE., do_efield=.FALSE., & 1443 charges=charges, dipoles=dipoles, quadrupoles=quadrupoles, forces_local=g_forces, & 1444 forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.FALSE.) 1445 CALL release_multi_type(multipoles) 1446 1447 ! Debug CHARGES-QUADRUPOLES 1448 task(1) = .TRUE. 1449 task(3) = .TRUE. 1450 charges = 0.0_dp 1451 dipoles = 0.0_dp 1452 quadrupoles = 0.0_dp 1453 r_forces = 0.0_dp 1454 g_forces = 0.0_dp 1455 e_field1 = 0.0_dp 1456 e_field2 = 0.0_dp 1457 g_pv = 0.0_dp 1458 r_pv = 0.0_dp 1459 g_energy = 0.0_dp 1460 r_energy = 0.0_dp 1461 e_neut = 0.0_dp 1462 e_self = 0.0_dp 1463 1464 CALL create_multi_type(multipoles, nparticles, 1, nparticles/2, "CHARGE", echarge=-1.0_dp, & 1465 random_stream=random_stream, charges=charges) 1466 CALL create_multi_type(multipoles, nparticles, nparticles/2+1, nparticles, "QUADRUPOLE", echarge=10.0_dp, & 1467 random_stream=random_stream, quadrupoles=quadrupoles) 1468 WRITE (iw, '("CHARGES",F15.9)') charges 1469 WRITE (iw, '("QUADRUPOLES",9F15.9)') quadrupoles 1470 CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, & 1471 debug_r_space) 1472 1473 WRITE (iw, *) "DEBUG ENERGY (CHARGE-QUADRUPOLE): ", debug_energy 1474 CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, & 1475 particle_set, local_particles, g_energy, r_energy, e_neut, e_self, & 1476 task, do_correction_bonded=.FALSE., do_forces=.TRUE., do_stress=.TRUE., do_efield=.FALSE., & 1477 charges=charges, dipoles=dipoles, quadrupoles=quadrupoles, forces_local=g_forces, & 1478 forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.FALSE.) 1479 CALL release_multi_type(multipoles) 1480 1481 ! Debug DIPOLES-QUADRUPOLES 1482 task(2) = .TRUE. 1483 task(3) = .TRUE. 1484 charges = 0.0_dp 1485 dipoles = 0.0_dp 1486 quadrupoles = 0.0_dp 1487 r_forces = 0.0_dp 1488 g_forces = 0.0_dp 1489 e_field1 = 0.0_dp 1490 e_field2 = 0.0_dp 1491 g_pv = 0.0_dp 1492 r_pv = 0.0_dp 1493 g_energy = 0.0_dp 1494 r_energy = 0.0_dp 1495 e_neut = 0.0_dp 1496 e_self = 0.0_dp 1497 1498 CALL create_multi_type(multipoles, nparticles, 1, nparticles/2, "DIPOLE", echarge=10000.0_dp, & 1499 random_stream=random_stream, dipoles=dipoles) 1500 CALL create_multi_type(multipoles, nparticles, nparticles/2+1, nparticles, "QUADRUPOLE", echarge=10000.0_dp, & 1501 random_stream=random_stream, quadrupoles=quadrupoles) 1502 WRITE (iw, '("DIPOLES",3F15.9)') dipoles 1503 WRITE (iw, '("QUADRUPOLES",9F15.9)') quadrupoles 1504 CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, & 1505 debug_r_space) 1506 1507 WRITE (iw, *) "DEBUG ENERGY (DIPOLE-QUADRUPOLE): ", debug_energy 1508 CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, & 1509 particle_set, local_particles, g_energy, r_energy, e_neut, e_self, & 1510 task, do_correction_bonded=.FALSE., do_forces=.TRUE., do_stress=.TRUE., do_efield=.FALSE., & 1511 charges=charges, dipoles=dipoles, quadrupoles=quadrupoles, forces_local=g_forces, & 1512 forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.FALSE.) 1513 CALL release_multi_type(multipoles) 1514 1515 ! Debug QUADRUPOLES-QUADRUPOLES 1516 task(3) = .TRUE. 1517 charges = 0.0_dp 1518 dipoles = 0.0_dp 1519 quadrupoles = 0.0_dp 1520 r_forces = 0.0_dp 1521 g_forces = 0.0_dp 1522 e_field1 = 0.0_dp 1523 e_field2 = 0.0_dp 1524 g_pv = 0.0_dp 1525 r_pv = 0.0_dp 1526 g_energy = 0.0_dp 1527 r_energy = 0.0_dp 1528 e_neut = 0.0_dp 1529 e_self = 0.0_dp 1530 1531 CALL create_multi_type(multipoles, nparticles, 1, nparticles/2, "QUADRUPOLE", echarge=-20000.0_dp, & 1532 random_stream=random_stream, quadrupoles=quadrupoles) 1533 CALL create_multi_type(multipoles, nparticles, nparticles/2+1, nparticles, "QUADRUPOLE", echarge=10000.0_dp, & 1534 random_stream=random_stream, quadrupoles=quadrupoles) 1535 WRITE (iw, '("QUADRUPOLES",9F15.9)') quadrupoles 1536 CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, & 1537 debug_r_space) 1538 1539 WRITE (iw, *) "DEBUG ENERGY (QUADRUPOLE-QUADRUPOLE): ", debug_energy 1540 CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, & 1541 particle_set, local_particles, g_energy, r_energy, e_neut, e_self, & 1542 task, do_correction_bonded=.FALSE., do_forces=.TRUE., do_stress=.TRUE., do_efield=.FALSE., & 1543 charges=charges, dipoles=dipoles, quadrupoles=quadrupoles, forces_local=g_forces, & 1544 forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.FALSE.) 1545 CALL release_multi_type(multipoles) 1546 1547 DEALLOCATE (charges) 1548 DEALLOCATE (dipoles) 1549 DEALLOCATE (quadrupoles) 1550 DEALLOCATE (r_forces) 1551 DEALLOCATE (g_forces) 1552 DEALLOCATE (e_field1) 1553 DEALLOCATE (e_field2) 1554 DEALLOCATE (g_pv) 1555 DEALLOCATE (r_pv) 1556 CALL delete_rng_stream(random_stream) 1557 1558 CONTAINS 1559! ************************************************************************************************** 1560!> \brief Debug routines for multipoles - low level - charge interactions 1561!> \param particle_set ... 1562!> \param cell ... 1563!> \param nonbond_env ... 1564!> \param multipoles ... 1565!> \param energy ... 1566!> \param debug_r_space ... 1567!> \date 05.2008 1568!> \author Teodoro Laino [tlaino] - University of Zurich - 05.2008 1569! ************************************************************************************************** 1570 SUBROUTINE debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, & 1571 energy, debug_r_space) 1572 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1573 TYPE(cell_type), POINTER :: cell 1574 TYPE(fist_nonbond_env_type), POINTER :: nonbond_env 1575 TYPE(multi_charge_type), DIMENSION(:), POINTER :: multipoles 1576 REAL(KIND=dp), INTENT(OUT) :: energy 1577 LOGICAL, INTENT(IN) :: debug_r_space 1578 1579 CHARACTER(len=*), PARAMETER :: routineN = 'debug_ewald_multipole_low', & 1580 routineP = moduleN//':'//routineN 1581 1582 INTEGER :: atom_a, atom_b, icell, iend, igrp, & 1583 ikind, ilist, ipair, istart, jcell, & 1584 jkind, k, k1, kcell, l, l1, ncells, & 1585 nkinds, npairs 1586 INTEGER, DIMENSION(:, :), POINTER :: list 1587 REAL(KIND=dp) :: fac_ij, q, r, rab2, rab2_max 1588 REAL(KIND=dp), DIMENSION(3) :: cell_v, cvi, rab, rab0, rm 1589 TYPE(fist_neighbor_type), POINTER :: nonbonded 1590 TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair 1591 TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update, r_last_update_pbc 1592 1593 energy = 0.0_dp 1594 CALL fist_nonbond_env_get(nonbond_env, nonbonded=nonbonded, natom_types=nkinds, & 1595 r_last_update=r_last_update, r_last_update_pbc=r_last_update_pbc) 1596 rab2_max = HUGE(0.0_dp) 1597 IF (debug_r_space) THEN 1598 ! This debugs the real space part of the multipole Ewald summation scheme 1599 ! Starting the force loop 1600 Lists: DO ilist = 1, nonbonded%nlists 1601 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist) 1602 npairs = neighbor_kind_pair%npairs 1603 IF (npairs == 0) CYCLE 1604 list => neighbor_kind_pair%list 1605 cvi = neighbor_kind_pair%cell_vector 1606 CALL matvec_3x3(cell_v, cell%hmat, cvi) 1607 Kind_Group_Loop: DO igrp = 1, neighbor_kind_pair%ngrp_kind 1608 istart = neighbor_kind_pair%grp_kind_start(igrp) 1609 iend = neighbor_kind_pair%grp_kind_end(igrp) 1610 ikind = neighbor_kind_pair%ij_kind(1, igrp) 1611 jkind = neighbor_kind_pair%ij_kind(2, igrp) 1612 Pairs: DO ipair = istart, iend 1613 fac_ij = 1.0_dp 1614 atom_a = list(1, ipair) 1615 atom_b = list(2, ipair) 1616 IF (atom_a == atom_b) fac_ij = 0.5_dp 1617 rab = r_last_update_pbc(atom_b)%r-r_last_update_pbc(atom_a)%r 1618 rab = rab+cell_v 1619 rab2 = rab(1)**2+rab(2)**2+rab(3)**2 1620 IF (rab2 <= rab2_max) THEN 1621 1622 DO k = 1, SIZE(multipoles(atom_a)%charge_typ) 1623 DO k1 = 1, SIZE(multipoles(atom_a)%charge_typ(k)%charge) 1624 1625 DO l = 1, SIZE(multipoles(atom_b)%charge_typ) 1626 DO l1 = 1, SIZE(multipoles(atom_b)%charge_typ(l)%charge) 1627 1628 rm = rab+multipoles(atom_b)%charge_typ(l)%pos(:, l1)-multipoles(atom_a)%charge_typ(k)%pos(:, k1) 1629 r = SQRT(DOT_PRODUCT(rm, rm)) 1630 q = multipoles(atom_b)%charge_typ(l)%charge(l1)*multipoles(atom_a)%charge_typ(k)%charge(k1) 1631 energy = energy+q/r*fac_ij 1632 END DO 1633 END DO 1634 1635 END DO 1636 END DO 1637 1638 END IF 1639 END DO Pairs 1640 END DO Kind_Group_Loop 1641 END DO Lists 1642 ELSE 1643 ncells = 6 1644 !Debugs the sum of real + space terms.. (Charge-Charge and Charge-Dipole should be anyway wrong but 1645 !all the other terms should be correct) 1646 DO atom_a = 1, SIZE(particle_set) 1647 DO atom_b = atom_a, SIZE(particle_set) 1648 fac_ij = 1.0_dp 1649 IF (atom_a == atom_b) fac_ij = 0.5_dp 1650 rab0 = r_last_update_pbc(atom_b)%r-r_last_update_pbc(atom_a)%r 1651 ! Loop over cells 1652 DO icell = -ncells, ncells 1653 DO jcell = -ncells, ncells 1654 DO kcell = -ncells, ncells 1655 cell_v = MATMUL(cell%hmat, REAL((/icell, jcell, kcell/), KIND=dp)) 1656 IF (ALL(cell_v == 0.0_dp) .AND. (atom_a == atom_b)) CYCLE 1657 rab = rab0+cell_v 1658 rab2 = rab(1)**2+rab(2)**2+rab(3)**2 1659 IF (rab2 <= rab2_max) THEN 1660 1661 DO k = 1, SIZE(multipoles(atom_a)%charge_typ) 1662 DO k1 = 1, SIZE(multipoles(atom_a)%charge_typ(k)%charge) 1663 1664 DO l = 1, SIZE(multipoles(atom_b)%charge_typ) 1665 DO l1 = 1, SIZE(multipoles(atom_b)%charge_typ(l)%charge) 1666 1667 rm = rab+multipoles(atom_b)%charge_typ(l)%pos(:, l1)-multipoles(atom_a)%charge_typ(k)%pos(:, k1) 1668 r = SQRT(DOT_PRODUCT(rm, rm)) 1669 q = multipoles(atom_b)%charge_typ(l)%charge(l1)*multipoles(atom_a)%charge_typ(k)%charge(k1) 1670 energy = energy+q/r*fac_ij 1671 END DO 1672 END DO 1673 1674 END DO 1675 END DO 1676 1677 END IF 1678 END DO 1679 END DO 1680 END DO 1681 END DO 1682 END DO 1683 END IF 1684 END SUBROUTINE debug_ewald_multipole_low 1685 1686! ************************************************************************************************** 1687!> \brief create multi_type for multipoles 1688!> \param multipoles ... 1689!> \param idim ... 1690!> \param istart ... 1691!> \param iend ... 1692!> \param label ... 1693!> \param echarge ... 1694!> \param random_stream ... 1695!> \param charges ... 1696!> \param dipoles ... 1697!> \param quadrupoles ... 1698!> \date 05.2008 1699!> \author Teodoro Laino [tlaino] - University of Zurich - 05.2008 1700! ************************************************************************************************** 1701 SUBROUTINE create_multi_type(multipoles, idim, istart, iend, label, echarge, & 1702 random_stream, charges, dipoles, quadrupoles) 1703 TYPE(multi_charge_type), DIMENSION(:), POINTER :: multipoles 1704 INTEGER, INTENT(IN) :: idim, istart, iend 1705 CHARACTER(LEN=*), INTENT(IN) :: label 1706 REAL(KIND=dp), INTENT(IN) :: echarge 1707 TYPE(rng_stream_type), POINTER :: random_stream 1708 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: charges 1709 REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: dipoles 1710 REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, & 1711 POINTER :: quadrupoles 1712 1713 CHARACTER(len=*), PARAMETER :: routineN = 'create_multi_type', & 1714 routineP = moduleN//':'//routineN 1715 1716 INTEGER :: i, isize, k, l, m 1717 REAL(KIND=dp) :: dx, r2, rvec(3), rvec1(3), rvec2(3) 1718 1719 IF (ASSOCIATED(multipoles)) THEN 1720 CPASSERT(SIZE(multipoles) == idim) 1721 ELSE 1722 ALLOCATE (multipoles(idim)) 1723 DO i = 1, idim 1724 NULLIFY (multipoles(i)%charge_typ) 1725 END DO 1726 END IF 1727 DO i = istart, iend 1728 IF (ASSOCIATED(multipoles(i)%charge_typ)) THEN 1729 ! make a copy of the array and enlarge the array type by 1 1730 isize = SIZE(multipoles(i)%charge_typ)+1 1731 ELSE 1732 isize = 1 1733 END IF 1734 CALL reallocate_charge_type(multipoles(i)%charge_typ, 1, isize) 1735 SELECT CASE (label) 1736 CASE ("CHARGE") 1737 CPASSERT(PRESENT(charges)) 1738 CPASSERT(ASSOCIATED(charges)) 1739 ALLOCATE (multipoles(i)%charge_typ(isize)%charge(1)) 1740 ALLOCATE (multipoles(i)%charge_typ(isize)%pos(3, 1)) 1741 1742 multipoles(i)%charge_typ(isize)%charge(1) = echarge 1743 multipoles(i)%charge_typ(isize)%pos(1:3, 1) = 0.0_dp 1744 charges(i) = charges(i)+echarge 1745 CASE ("DIPOLE") 1746 dx = 1.0E-4_dp 1747 CPASSERT(PRESENT(dipoles)) 1748 CPASSERT(ASSOCIATED(dipoles)) 1749 ALLOCATE (multipoles(i)%charge_typ(isize)%charge(2)) 1750 ALLOCATE (multipoles(i)%charge_typ(isize)%pos(3, 2)) 1751 CALL random_numbers(rvec, random_stream) 1752 rvec = rvec/(2.0_dp*SQRT(DOT_PRODUCT(rvec, rvec)))*dx 1753 multipoles(i)%charge_typ(isize)%charge(1) = echarge 1754 multipoles(i)%charge_typ(isize)%pos(1:3, 1) = rvec 1755 multipoles(i)%charge_typ(isize)%charge(2) = -echarge 1756 multipoles(i)%charge_typ(isize)%pos(1:3, 2) = -rvec 1757 1758 dipoles(:, i) = dipoles(:, i)+2.0_dp*echarge*rvec 1759 CASE ("QUADRUPOLE") 1760 dx = 1.0E-2_dp 1761 CPASSERT(PRESENT(quadrupoles)) 1762 CPASSERT(ASSOCIATED(quadrupoles)) 1763 ALLOCATE (multipoles(i)%charge_typ(isize)%charge(4)) 1764 ALLOCATE (multipoles(i)%charge_typ(isize)%pos(3, 4)) 1765 CALL random_numbers(rvec1, random_stream) 1766 CALL random_numbers(rvec2, random_stream) 1767 rvec1 = rvec1/SQRT(DOT_PRODUCT(rvec1, rvec1)) 1768 rvec2 = rvec2-DOT_PRODUCT(rvec2, rvec1)*rvec1 1769 rvec2 = rvec2/SQRT(DOT_PRODUCT(rvec2, rvec2)) 1770 ! 1771 rvec1 = rvec1/2.0_dp*dx 1772 rvec2 = rvec2/2.0_dp*dx 1773 ! + (4) ^ - (1) 1774 ! |rvec2 1775 ! | 1776 ! 0------> rvec1 1777 ! 1778 ! 1779 ! - (3) + (2) 1780 multipoles(i)%charge_typ(isize)%charge(1) = -echarge 1781 multipoles(i)%charge_typ(isize)%pos(1:3, 1) = rvec1+rvec2 1782 multipoles(i)%charge_typ(isize)%charge(2) = echarge 1783 multipoles(i)%charge_typ(isize)%pos(1:3, 2) = rvec1-rvec2 1784 multipoles(i)%charge_typ(isize)%charge(3) = -echarge 1785 multipoles(i)%charge_typ(isize)%pos(1:3, 3) = -rvec1-rvec2 1786 multipoles(i)%charge_typ(isize)%charge(4) = echarge 1787 multipoles(i)%charge_typ(isize)%pos(1:3, 4) = -rvec1+rvec2 1788 1789 DO k = 1, 4 1790 r2 = DOT_PRODUCT(multipoles(i)%charge_typ(isize)%pos(:, k), multipoles(i)%charge_typ(isize)%pos(:, k)) 1791 DO l = 1, 3 1792 DO m = 1, 3 1793 quadrupoles(m, l, i) = quadrupoles(m, l, i)+3.0_dp*0.5_dp*multipoles(i)%charge_typ(isize)%charge(k)* & 1794 multipoles(i)%charge_typ(isize)%pos(l, k)* & 1795 multipoles(i)%charge_typ(isize)%pos(m, k) 1796 IF (m == l) quadrupoles(m, l, i) = quadrupoles(m, l, i)-0.5_dp*multipoles(i)%charge_typ(isize)%charge(k)*r2 1797 END DO 1798 END DO 1799 END DO 1800 1801 END SELECT 1802 END DO 1803 END SUBROUTINE create_multi_type 1804 1805! ************************************************************************************************** 1806!> \brief release multi_type for multipoles 1807!> \param multipoles ... 1808!> \date 05.2008 1809!> \author Teodoro Laino [tlaino] - University of Zurich - 05.2008 1810! ************************************************************************************************** 1811 SUBROUTINE release_multi_type(multipoles) 1812 TYPE(multi_charge_type), DIMENSION(:), POINTER :: multipoles 1813 1814 CHARACTER(len=*), PARAMETER :: routineN = 'release_multi_type', & 1815 routineP = moduleN//':'//routineN 1816 1817 INTEGER :: i, j 1818 1819 IF (ASSOCIATED(multipoles)) THEN 1820 DO i = 1, SIZE(multipoles) 1821 DO j = 1, SIZE(multipoles(i)%charge_typ) 1822 DEALLOCATE (multipoles(i)%charge_typ(j)%charge) 1823 DEALLOCATE (multipoles(i)%charge_typ(j)%pos) 1824 END DO 1825 DEALLOCATE (multipoles(i)%charge_typ) 1826 END DO 1827 END IF 1828 END SUBROUTINE release_multi_type 1829 1830! ************************************************************************************************** 1831!> \brief reallocates multi_type for multipoles 1832!> \param charge_typ ... 1833!> \param istart ... 1834!> \param iend ... 1835!> \date 05.2008 1836!> \author Teodoro Laino [tlaino] - University of Zurich - 05.2008 1837! ************************************************************************************************** 1838 SUBROUTINE reallocate_charge_type(charge_typ, istart, iend) 1839 TYPE(charge_mono_type), DIMENSION(:), POINTER :: charge_typ 1840 INTEGER, INTENT(IN) :: istart, iend 1841 1842 CHARACTER(len=*), PARAMETER :: routineN = 'reallocate_charge_type', & 1843 routineP = moduleN//':'//routineN 1844 1845 INTEGER :: i, isize, j, jsize, jsize1, jsize2 1846 TYPE(charge_mono_type), DIMENSION(:), POINTER :: charge_typ_bk 1847 1848 IF (ASSOCIATED(charge_typ)) THEN 1849 isize = SIZE(charge_typ) 1850 ALLOCATE (charge_typ_bk(1:isize)) 1851 DO j = 1, isize 1852 jsize = SIZE(charge_typ(j)%charge) 1853 ALLOCATE (charge_typ_bk(j)%charge(jsize)) 1854 jsize1 = SIZE(charge_typ(j)%pos, 1) 1855 jsize2 = SIZE(charge_typ(j)%pos, 2) 1856 ALLOCATE (charge_typ_bk(j)%pos(jsize1, jsize2)) 1857 charge_typ_bk(j)%pos = charge_typ(j)%pos 1858 charge_typ_bk(j)%charge = charge_typ(j)%charge 1859 END DO 1860 DO j = 1, SIZE(charge_typ) 1861 DEALLOCATE (charge_typ(j)%charge) 1862 DEALLOCATE (charge_typ(j)%pos) 1863 END DO 1864 DEALLOCATE (charge_typ) 1865 ! Reallocate 1866 ALLOCATE (charge_typ_bk(istart:iend)) 1867 DO i = istart, isize 1868 jsize = SIZE(charge_typ_bk(j)%charge) 1869 ALLOCATE (charge_typ(j)%charge(jsize)) 1870 jsize1 = SIZE(charge_typ_bk(j)%pos, 1) 1871 jsize2 = SIZE(charge_typ_bk(j)%pos, 2) 1872 ALLOCATE (charge_typ(j)%pos(jsize1, jsize2)) 1873 charge_typ(j)%pos = charge_typ_bk(j)%pos 1874 charge_typ(j)%charge = charge_typ_bk(j)%charge 1875 END DO 1876 DO j = 1, SIZE(charge_typ_bk) 1877 DEALLOCATE (charge_typ_bk(j)%charge) 1878 DEALLOCATE (charge_typ_bk(j)%pos) 1879 END DO 1880 DEALLOCATE (charge_typ_bk) 1881 ELSE 1882 ALLOCATE (charge_typ(istart:iend)) 1883 END IF 1884 1885 END SUBROUTINE reallocate_charge_type 1886 1887 END SUBROUTINE debug_ewald_multipoles 1888 1889! ************************************************************************************************** 1890!> \brief Routine to debug potential, field and electric field gradients 1891!> \param ewald_env ... 1892!> \param ewald_pw ... 1893!> \param nonbond_env ... 1894!> \param cell ... 1895!> \param particle_set ... 1896!> \param local_particles ... 1897!> \param radii ... 1898!> \param charges ... 1899!> \param dipoles ... 1900!> \param quadrupoles ... 1901!> \param task ... 1902!> \param iw ... 1903!> \param atomic_kind_set ... 1904!> \param mm_section ... 1905!> \date 05.2008 1906!> \author Teodoro Laino [tlaino] - University of Zurich - 05.2008 1907! ************************************************************************************************** 1908 SUBROUTINE debug_ewald_multipoles_fields(ewald_env, ewald_pw, nonbond_env, cell, & 1909 particle_set, local_particles, radii, charges, dipoles, quadrupoles, task, iw, & 1910 atomic_kind_set, mm_section) 1911 TYPE(ewald_environment_type), POINTER :: ewald_env 1912 TYPE(ewald_pw_type), POINTER :: ewald_pw 1913 TYPE(fist_nonbond_env_type), POINTER :: nonbond_env 1914 TYPE(cell_type), POINTER :: cell 1915 TYPE(particle_type), POINTER :: particle_set(:) 1916 TYPE(distribution_1d_type), POINTER :: local_particles 1917 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: radii, charges 1918 REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: dipoles 1919 REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, & 1920 POINTER :: quadrupoles 1921 LOGICAL, DIMENSION(3), INTENT(IN) :: task 1922 INTEGER, INTENT(IN) :: iw 1923 TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:) 1924 TYPE(section_vals_type), POINTER :: mm_section 1925 1926 CHARACTER(len=*), PARAMETER :: routineN = 'debug_ewald_multipoles_fields', & 1927 routineP = moduleN//':'//routineN 1928 1929 INTEGER :: i, iparticle_kind, j, k, & 1930 nparticle_local, nparticles 1931 REAL(KIND=dp) :: coord(3), dq, e_neut, e_self, efield1n(3), efield2n(3, 3), ene(2), & 1932 energy_glob, energy_local, enev(3, 2), o_tot_ene, pot, pv_glob(3, 3), pv_local(3, 3), & 1933 tot_ene 1934 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: efield1, efield2, forces_glob, & 1935 forces_local 1936 REAL(KIND=dp), DIMENSION(:), POINTER :: efield0, lcharges 1937 TYPE(cp_logger_type), POINTER :: logger 1938 TYPE(particle_type), DIMENSION(:), POINTER :: core_particle_set, shell_particle_set 1939 1940 NULLIFY (lcharges, shell_particle_set, core_particle_set) 1941 NULLIFY (logger) 1942 logger => cp_get_default_logger() 1943 1944 nparticles = SIZE(particle_set) 1945 nparticle_local = 0 1946 DO iparticle_kind = 1, SIZE(local_particles%n_el) 1947 nparticle_local = nparticle_local+local_particles%n_el(iparticle_kind) 1948 END DO 1949 ALLOCATE (lcharges(nparticles)) 1950 ALLOCATE (forces_glob(3, nparticles)) 1951 ALLOCATE (forces_local(3, nparticle_local)) 1952 ALLOCATE (efield0(nparticles)) 1953 ALLOCATE (efield1(3, nparticles)) 1954 ALLOCATE (efield2(9, nparticles)) 1955 forces_glob = 0.0_dp 1956 forces_local = 0.0_dp 1957 efield0 = 0.0_dp 1958 efield1 = 0.0_dp 1959 efield2 = 0.0_dp 1960 pv_local = 0.0_dp 1961 pv_glob = 0.0_dp 1962 energy_glob = 0.0_dp 1963 energy_local = 0.0_dp 1964 e_neut = 0.0_dp 1965 e_self = 0.0_dp 1966 CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, particle_set, & 1967 local_particles, energy_local, energy_glob, e_neut, e_self, task, .FALSE., .TRUE., .TRUE., & 1968 .TRUE., radii, charges, dipoles, quadrupoles, forces_local, forces_glob, pv_local, pv_glob, & 1969 efield0, efield1, efield2, iw, do_debug=.FALSE.) 1970 o_tot_ene = energy_local+energy_glob+e_neut+e_self 1971 WRITE (iw, *) "TOTAL ENERGY :: ========>", o_tot_ene 1972 ! Debug Potential 1973 dq = 0.001_dp 1974 tot_ene = 0.0_dp 1975 DO i = 1, nparticles 1976 DO k = 1, 2 1977 lcharges = charges 1978 lcharges(i) = charges(i)+(-1.0_dp)**k*dq 1979 forces_glob = 0.0_dp 1980 forces_local = 0.0_dp 1981 pv_local = 0.0_dp 1982 pv_glob = 0.0_dp 1983 energy_glob = 0.0_dp 1984 energy_local = 0.0_dp 1985 e_neut = 0.0_dp 1986 e_self = 0.0_dp 1987 CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, particle_set, & 1988 local_particles, energy_local, energy_glob, e_neut, e_self, & 1989 task, .FALSE., .FALSE., .FALSE., .FALSE., radii, & 1990 lcharges, dipoles, quadrupoles, iw=iw, do_debug=.FALSE.) 1991 ene(k) = energy_local+energy_glob+e_neut+e_self 1992 END DO 1993 pot = (ene(2)-ene(1))/(2.0_dp*dq) 1994 WRITE (iw, '(A,I8,3(A,F15.9))') "POTENTIAL FOR ATOM: ", i, " NUMERICAL: ", pot, " ANALYTICAL: ", efield0(i), & 1995 " ERROR: ", pot-efield0(i) 1996 tot_ene = tot_ene+0.5_dp*efield0(i)*charges(i) 1997 END DO 1998 WRITE (iw, *) "ENERGIES: ", o_tot_ene, tot_ene, o_tot_ene-tot_ene 1999 WRITE (iw, '(/,/,/)') 2000 ! Debug Field 2001 dq = 0.001_dp 2002 DO i = 1, nparticles 2003 coord = particle_set(i)%r 2004 DO j = 1, 3 2005 DO k = 1, 2 2006 particle_set(i)%r(j) = coord(j)+(-1.0_dp)**k*dq 2007 2008 ! Rebuild neighbor lists 2009 CALL list_control(atomic_kind_set, particle_set, local_particles, & 2010 cell, nonbond_env, logger%para_env, mm_section, & 2011 shell_particle_set, core_particle_set) 2012 2013 forces_glob = 0.0_dp 2014 forces_local = 0.0_dp 2015 pv_local = 0.0_dp 2016 pv_glob = 0.0_dp 2017 energy_glob = 0.0_dp 2018 energy_local = 0.0_dp 2019 e_neut = 0.0_dp 2020 e_self = 0.0_dp 2021 efield0 = 0.0_dp 2022 CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, particle_set, & 2023 local_particles, energy_local, energy_glob, e_neut, e_self, & 2024 task, .FALSE., .TRUE., .TRUE., .TRUE., radii, & 2025 charges, dipoles, quadrupoles, forces_local, forces_glob, & 2026 pv_local, pv_glob, efield0, iw=iw, do_debug=.FALSE.) 2027 ene(k) = efield0(i) 2028 particle_set(i)%r(j) = coord(j) 2029 END DO 2030 efield1n(j) = -(ene(2)-ene(1))/(2.0_dp*dq) 2031 END DO 2032 WRITE (iw, '(/,A,I8)') "FIELD FOR ATOM: ", i 2033 WRITE (iw, '(A,3F15.9)') " NUMERICAL: ", efield1n, " ANALYTICAL: ", efield1(:, i), & 2034 " ERROR: ", efield1n-efield1(:, i) 2035 IF (task(2)) THEN 2036 tot_ene = tot_ene-0.5_dp*DOT_PRODUCT(efield1(:, i), dipoles(:, i)) 2037 END IF 2038 END DO 2039 WRITE (iw, *) "ENERGIES: ", o_tot_ene, tot_ene, o_tot_ene-tot_ene 2040 2041 ! Debug Field Gradient 2042 dq = 0.0001_dp 2043 DO i = 1, nparticles 2044 coord = particle_set(i)%r 2045 DO j = 1, 3 2046 DO k = 1, 2 2047 particle_set(i)%r(j) = coord(j)+(-1.0_dp)**k*dq 2048 2049 ! Rebuild neighbor lists 2050 CALL list_control(atomic_kind_set, particle_set, local_particles, & 2051 cell, nonbond_env, logger%para_env, mm_section, & 2052 shell_particle_set, core_particle_set) 2053 2054 forces_glob = 0.0_dp 2055 forces_local = 0.0_dp 2056 pv_local = 0.0_dp 2057 pv_glob = 0.0_dp 2058 energy_glob = 0.0_dp 2059 energy_local = 0.0_dp 2060 e_neut = 0.0_dp 2061 e_self = 0.0_dp 2062 efield1 = 0.0_dp 2063 CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, particle_set, & 2064 local_particles, energy_local, energy_glob, e_neut, e_self, & 2065 task, .FALSE., .TRUE., .TRUE., .TRUE., radii, & 2066 charges, dipoles, quadrupoles, forces_local, forces_glob, & 2067 pv_local, pv_glob, efield1=efield1, iw=iw, do_debug=.FALSE.) 2068 enev(:, k) = efield1(:, i) 2069 particle_set(i)%r(j) = coord(j) 2070 END DO 2071 efield2n(:, j) = (enev(:, 2)-enev(:, 1))/(2.0_dp*dq) 2072 END DO 2073 WRITE (iw, '(/,A,I8)') "FIELD GRADIENT FOR ATOM: ", i 2074 WRITE (iw, '(A,9F15.9)') " NUMERICAL: ", efield2n, & 2075 " ANALYTICAL: ", efield2(:, i), & 2076 " ERROR: ", RESHAPE(efield2n, (/9/))-efield2(:, i) 2077 END DO 2078 END SUBROUTINE debug_ewald_multipoles_fields 2079 2080! ************************************************************************************************** 2081!> \brief Routine to debug potential, field and electric field gradients 2082!> \param ewald_env ... 2083!> \param ewald_pw ... 2084!> \param nonbond_env ... 2085!> \param cell ... 2086!> \param particle_set ... 2087!> \param local_particles ... 2088!> \param radii ... 2089!> \param charges ... 2090!> \param dipoles ... 2091!> \param quadrupoles ... 2092!> \param task ... 2093!> \param iw ... 2094!> \date 05.2008 2095!> \author Teodoro Laino [tlaino] - University of Zurich - 05.2008 2096! ************************************************************************************************** 2097 SUBROUTINE debug_ewald_multipoles_fields2(ewald_env, ewald_pw, nonbond_env, cell, & 2098 particle_set, local_particles, radii, charges, dipoles, quadrupoles, task, iw) 2099 TYPE(ewald_environment_type), POINTER :: ewald_env 2100 TYPE(ewald_pw_type), POINTER :: ewald_pw 2101 TYPE(fist_nonbond_env_type), POINTER :: nonbond_env 2102 TYPE(cell_type), POINTER :: cell 2103 TYPE(particle_type), POINTER :: particle_set(:) 2104 TYPE(distribution_1d_type), POINTER :: local_particles 2105 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: radii, charges 2106 REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: dipoles 2107 REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, & 2108 POINTER :: quadrupoles 2109 LOGICAL, DIMENSION(3), INTENT(IN) :: task 2110 INTEGER, INTENT(IN) :: iw 2111 2112 CHARACTER(len=*), PARAMETER :: routineN = 'debug_ewald_multipoles_fields2', & 2113 routineP = moduleN//':'//routineN 2114 2115 INTEGER :: i, ind, iparticle_kind, j, k, & 2116 nparticle_local, nparticles 2117 REAL(KIND=dp) :: e_neut, e_self, energy_glob, & 2118 energy_local, o_tot_ene, prod, & 2119 pv_glob(3, 3), pv_local(3, 3), tot_ene 2120 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: efield1, efield2, forces_glob, & 2121 forces_local 2122 REAL(KIND=dp), DIMENSION(:), POINTER :: efield0 2123 TYPE(cp_logger_type), POINTER :: logger 2124 2125 NULLIFY (logger) 2126 logger => cp_get_default_logger() 2127 2128 nparticles = SIZE(particle_set) 2129 nparticle_local = 0 2130 DO iparticle_kind = 1, SIZE(local_particles%n_el) 2131 nparticle_local = nparticle_local+local_particles%n_el(iparticle_kind) 2132 END DO 2133 ALLOCATE (forces_glob(3, nparticles)) 2134 ALLOCATE (forces_local(3, nparticle_local)) 2135 ALLOCATE (efield0(nparticles)) 2136 ALLOCATE (efield1(3, nparticles)) 2137 ALLOCATE (efield2(9, nparticles)) 2138 forces_glob = 0.0_dp 2139 forces_local = 0.0_dp 2140 efield0 = 0.0_dp 2141 efield1 = 0.0_dp 2142 efield2 = 0.0_dp 2143 pv_local = 0.0_dp 2144 pv_glob = 0.0_dp 2145 energy_glob = 0.0_dp 2146 energy_local = 0.0_dp 2147 e_neut = 0.0_dp 2148 e_self = 0.0_dp 2149 CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, particle_set, & 2150 local_particles, energy_local, energy_glob, e_neut, e_self, task, .FALSE., .TRUE., .TRUE., & 2151 .TRUE., radii, charges, dipoles, quadrupoles, forces_local, forces_glob, pv_local, pv_glob, & 2152 efield0, efield1, efield2, iw, do_debug=.FALSE.) 2153 o_tot_ene = energy_local+energy_glob+e_neut+e_self 2154 WRITE (iw, *) "TOTAL ENERGY :: ========>", o_tot_ene 2155 2156 ! Debug Potential 2157 tot_ene = 0.0_dp 2158 IF (task(1)) THEN 2159 DO i = 1, nparticles 2160 tot_ene = tot_ene+0.5_dp*efield0(i)*charges(i) 2161 END DO 2162 WRITE (iw, *) "ENERGIES: ", o_tot_ene, tot_ene, o_tot_ene-tot_ene 2163 WRITE (iw, '(/,/,/)') 2164 END IF 2165 2166 ! Debug Field 2167 IF (task(2)) THEN 2168 DO i = 1, nparticles 2169 tot_ene = tot_ene-0.5_dp*DOT_PRODUCT(efield1(:, i), dipoles(:, i)) 2170 END DO 2171 WRITE (iw, *) "ENERGIES: ", o_tot_ene, tot_ene, o_tot_ene-tot_ene 2172 WRITE (iw, '(/,/,/)') 2173 END IF 2174 2175 ! Debug Field Gradient 2176 IF (task(3)) THEN 2177 DO i = 1, nparticles 2178 ind = 0 2179 prod = 0.0_dp 2180 DO j = 1, 3 2181 DO k = 1, 3 2182 ind = ind+1 2183 prod = prod+efield2(ind, i)*quadrupoles(j, k, i) 2184 END DO 2185 END DO 2186 tot_ene = tot_ene-0.5_dp*(1.0_dp/3.0_dp)*prod 2187 END DO 2188 WRITE (iw, *) "ENERGIES: ", o_tot_ene, tot_ene, o_tot_ene-tot_ene 2189 WRITE (iw, '(/,/,/)') 2190 END IF 2191 2192 END SUBROUTINE debug_ewald_multipoles_fields2 2193 2194END MODULE ewalds_multipole 2195