1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief implementation of dipole and three-body part of Siepmann-Sprik potential 8!> dipole term: 3rd term in Eq. (1) in J. Chem. Phys., Vol. 102, p.511 9!> three-body term: Eq. (4) in J. Chem. Phys., Vol. 102, p. 511 10!> remaining terms of Siepmann-Sprik potential can be given via the GENPOT section 11!> \par History 12!> 12.2012 created 13!> \author Dorothea Golze 14! ************************************************************************************************** 15MODULE manybody_siepmann 16 17 USE atomic_kind_types, ONLY: get_atomic_kind 18 USE cell_types, ONLY: cell_type,& 19 pbc 20 USE cp_log_handling, ONLY: cp_get_default_logger,& 21 cp_logger_type 22 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 23 cp_print_key_unit_nr 24 USE cp_para_types, ONLY: cp_para_env_type 25 USE fist_neighbor_list_types, ONLY: fist_neighbor_type,& 26 neighbor_kind_pairs_type 27 USE fist_nonbond_env_types, ONLY: pos_type 28 USE input_section_types, ONLY: section_vals_type 29 USE kinds, ONLY: dp 30 USE mathlib, ONLY: matvec_3x3 31 USE message_passing, ONLY: mp_sum 32 USE pair_potential_types, ONLY: pair_potential_pp_type,& 33 pair_potential_single_type,& 34 siepmann_pot_type,& 35 siepmann_type 36 USE particle_types, ONLY: particle_type 37 USE util, ONLY: sort 38#include "./base/base_uses.f90" 39 40 IMPLICIT NONE 41 42 PRIVATE 43 PUBLIC :: setup_siepmann_arrays, destroy_siepmann_arrays, & 44 siepmann_energy, siepmann_forces_v2, siepmann_forces_v3, & 45 print_nr_ions_siepmann 46 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'manybody_siepmann' 47 48CONTAINS 49 50! ************************************************************************************************** 51!> \brief energy of two-body dipole term and three-body term 52!> \param pot_loc ... 53!> \param siepmann ... 54!> \param r_last_update_pbc ... 55!> \param atom_a ... 56!> \param atom_b ... 57!> \param nloc_size ... 58!> \param full_loc_list ... 59!> \param cell_v ... 60!> \param cell ... 61!> \param drij ... 62!> \param particle_set ... 63!> \param nr_oh number of OH- ions near surface 64!> \param nr_h3o number of hydronium ions near surface 65!> \param nr_o number of O^(2-) ions near surface 66!> \author Dorothea Golze - 11.2012 - University of Zurich 67! ************************************************************************************************** 68 SUBROUTINE siepmann_energy(pot_loc, siepmann, r_last_update_pbc, atom_a, atom_b, & 69 nloc_size, full_loc_list, cell_v, cell, drij, particle_set, & 70 nr_oh, nr_h3o, nr_o) 71 72 REAL(KIND=dp), INTENT(OUT) :: pot_loc 73 TYPE(siepmann_pot_type), POINTER :: siepmann 74 TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc 75 INTEGER, INTENT(IN) :: atom_a, atom_b, nloc_size 76 INTEGER, DIMENSION(2, 1:nloc_size) :: full_loc_list 77 REAL(KIND=dp), DIMENSION(3) :: cell_v 78 TYPE(cell_type), POINTER :: cell 79 REAL(KIND=dp) :: drij 80 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 81 INTEGER, INTENT(INOUT) :: nr_oh, nr_h3o, nr_o 82 83 CHARACTER(LEN=*), PARAMETER :: routineN = 'siepmann_energy', & 84 routineP = moduleN//':'//routineN 85 86 REAL(KIND=dp) :: a_ij, D, E, f2, Phi_ij, pot_loc_v2, & 87 pot_loc_v3 88 89 a_ij = siep_a_ij(siepmann, r_last_update_pbc, atom_a, atom_b, nloc_size, & 90 full_loc_list, cell_v, siepmann%rcutsq, particle_set, & 91 cell) 92 Phi_ij = siep_Phi_ij(siepmann, r_last_update_pbc, atom_a, atom_b, & 93 cell_v, cell, siepmann%rcutsq, particle_set, nr_oh, nr_h3o, nr_o) 94 f2 = siep_f2(siepmann, drij) 95 D = siepmann%D 96 E = siepmann%E 97 98 !two-body part --> dipole term 99 pot_loc_v2 = -D*f2*drij**(-3)*Phi_ij 100 101 !three-body part 102 pot_loc_v3 = E*f2*drij**(-siepmann%beta)*a_ij 103 104 pot_loc = pot_loc_v2 + pot_loc_v3 105 106 END SUBROUTINE siepmann_energy 107 108! ************************************************************************************************** 109!> \brief f2(r) corresponds to Equation (2) in J. Chem. Phys., Vol. 102, p.511 110!> \param siepmann ... 111!> \param r distance between oxygen and metal atom 112!> \return ... 113! ************************************************************************************************** 114 FUNCTION siep_f2(siepmann, r) 115 TYPE(siepmann_pot_type), POINTER :: siepmann 116 REAL(KIND=dp), INTENT(IN) :: r 117 REAL(KIND=dp) :: siep_f2 118 119 CHARACTER(LEN=*), PARAMETER :: routineN = 'siep_f2', routineP = moduleN//':'//routineN 120 121 REAL(KIND=dp) :: rcut 122 123 rcut = SQRT(siepmann%rcutsq) 124 siep_f2 = 0.0_dp 125 IF (r < rcut) THEN 126 siep_f2 = EXP(siepmann%B/(r - rcut)) 127 END IF 128 END FUNCTION siep_f2 129 130! ************************************************************************************************** 131!> \brief f2(r)_d derivative of f2 132!> \param siepmann ... 133!> \param r distance between oxygen and metal atom 134!> \return ... 135! ************************************************************************************************** 136 FUNCTION siep_f2_d(siepmann, r) 137 TYPE(siepmann_pot_type), POINTER :: siepmann 138 REAL(KIND=dp), INTENT(IN) :: r 139 REAL(KIND=dp) :: siep_f2_d 140 141 CHARACTER(LEN=*), PARAMETER :: routineN = 'siep_f2_d', routineP = moduleN//':'//routineN 142 143 REAL(KIND=dp) :: B, rcut 144 145 rcut = SQRT(siepmann%rcutsq) 146 B = siepmann%B 147 siep_f2_d = 0.0_dp 148 IF (r < rcut) THEN 149 siep_f2_d = -B*EXP(B/(r - rcut))/(r - rcut)**2 150 END IF 151 152 END FUNCTION siep_f2_d 153 154! ************************************************************************************************** 155!> \brief exponential part of three-body term, see Equation (4) in J. Chem. 156!> Phys., Vol. 102, p.511 157!> \param siepmann ... 158!> \param r_last_update_pbc ... 159!> \param iparticle ... 160!> \param jparticle ... 161!> \param n_loc_size ... 162!> \param full_loc_list ... 163!> \param cell_v ... 164!> \param rcutsq ... 165!> \param particle_set ... 166!> \param cell ... 167!> \return ... 168!> \par History 169!> Using a local list of neighbors 170! ************************************************************************************************** 171 FUNCTION siep_a_ij(siepmann, r_last_update_pbc, iparticle, jparticle, n_loc_size, & 172 full_loc_list, cell_v, rcutsq, particle_set, cell) 173 TYPE(siepmann_pot_type), POINTER :: siepmann 174 TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc 175 INTEGER, INTENT(IN) :: iparticle, jparticle, n_loc_size 176 INTEGER, DIMENSION(2, 1:n_loc_size) :: full_loc_list 177 REAL(KIND=dp), DIMENSION(3) :: cell_v 178 REAL(KIND=dp), INTENT(IN) :: rcutsq 179 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 180 TYPE(cell_type), POINTER :: cell 181 REAL(KIND=dp) :: siep_a_ij 182 183 CHARACTER(LEN=*), PARAMETER :: routineN = 'siep_a_ij', routineP = moduleN//':'//routineN 184 185 CHARACTER(LEN=2) :: element_symbol 186 INTEGER :: ilist, kparticle 187 REAL(KIND=dp) :: costheta, drji, drjk, F, rab2_max, & 188 rji(3), rjk(3), theta 189 190 siep_a_ij = 0.0_dp 191 rab2_max = rcutsq 192 F = siepmann%F 193 CALL get_atomic_kind(atomic_kind=particle_set(iparticle)%atomic_kind, & 194 element_symbol=element_symbol) 195 IF (element_symbol /= "O") RETURN 196 rji(:) = -1.0_dp*(r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v) 197 drji = SQRT(DOT_PRODUCT(rji, rji)) 198 DO ilist = 1, n_loc_size 199 kparticle = full_loc_list(2, ilist) 200 IF (kparticle == jparticle) CYCLE 201 rjk(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(kparticle)%r(:), cell) 202 drjk = DOT_PRODUCT(rjk, rjk) 203 IF (drjk > rab2_max) CYCLE 204 drjk = SQRT(drjk) 205 costheta = DOT_PRODUCT(rji, rjk)/(drji*drjk) 206 IF (costheta < -1.0_dp) costheta = -1.0_dp 207 IF (costheta > +1.0_dp) costheta = +1.0_dp 208 theta = ACOS(costheta) 209 siep_a_ij = siep_a_ij + EXP(F*(COS(theta/2.0_dp))**2) 210 END DO 211 END FUNCTION siep_a_ij 212 213! ************************************************************************************************** 214!> \brief derivative of a_ij 215!> \param siepmann ... 216!> \param r_last_update_pbc ... 217!> \param iparticle ... 218!> \param jparticle ... 219!> \param f_nonbond ... 220!> \param prefactor ... 221!> \param n_loc_size ... 222!> \param full_loc_list ... 223!> \param cell_v ... 224!> \param cell ... 225!> \param rcutsq ... 226!> \param use_virial ... 227!> \par History 228!> Using a local list of neighbors 229! ************************************************************************************************** 230 SUBROUTINE siep_a_ij_d(siepmann, r_last_update_pbc, iparticle, jparticle, f_nonbond, & 231 prefactor, n_loc_size, full_loc_list, & 232 cell_v, cell, rcutsq, use_virial) 233 TYPE(siepmann_pot_type), POINTER :: siepmann 234 TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc 235 INTEGER, INTENT(IN) :: iparticle, jparticle 236 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: f_nonbond 237 REAL(KIND=dp), INTENT(IN) :: prefactor 238 INTEGER, INTENT(IN) :: n_loc_size 239 INTEGER, DIMENSION(2, 1:n_loc_size) :: full_loc_list 240 REAL(KIND=dp), DIMENSION(3) :: cell_v 241 TYPE(cell_type), POINTER :: cell 242 REAL(KIND=dp), INTENT(IN) :: rcutsq 243 LOGICAL, INTENT(IN) :: use_virial 244 245 CHARACTER(LEN=*), PARAMETER :: routineN = 'siep_a_ij_d', routineP = moduleN//':'//routineN 246 247 INTEGER :: ilist, kparticle, nparticle 248 REAL(KIND=dp) :: costheta, d_expterm, dcos_thetahalf, & 249 drji, drjk, F, rab2_max, theta 250 REAL(KIND=dp), DIMENSION(3) :: dcosdri, dcosdrj, dcosdrk, dri, drj, & 251 drk, rji, rji_hat, rjk, rjk_hat 252 253 rab2_max = rcutsq 254 F = siepmann%F 255 256 rji(:) = -1.0_dp*(r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v) 257 drji = SQRT(DOT_PRODUCT(rji, rji)) 258 rji_hat(:) = rji(:)/drji 259 260 nparticle = SIZE(r_last_update_pbc) 261 DO ilist = 1, n_loc_size 262 kparticle = full_loc_list(2, ilist) 263 IF (kparticle == jparticle) CYCLE 264 rjk(:) = pbc(r_last_update_pbc(jparticle)%r(:), r_last_update_pbc(kparticle)%r(:), cell) 265 drjk = DOT_PRODUCT(rjk, rjk) 266 IF (drjk > rab2_max) CYCLE 267 drjk = SQRT(drjk) 268 rjk_hat(:) = rjk(:)/drjk 269 costheta = DOT_PRODUCT(rji, rjk)/(drji*drjk) 270 IF (costheta < -1.0_dp) costheta = -1.0_dp 271 IF (costheta > +1.0_dp) costheta = +1.0_dp 272 273 dcosdri(:) = (1.0_dp/(drji))*(rjk_hat(:) - costheta*rji_hat(:)) 274 dcosdrk(:) = (1.0_dp/(drjk))*(rji_hat(:) - costheta*rjk_hat(:)) 275 dcosdrj(:) = -(dcosdri(:) + dcosdrk(:)) 276 277 theta = ACOS(costheta) 278 dcos_thetahalf = -1.0_dp/(2.0_dp*SQRT(1 - costheta**2)) 279 d_expterm = -2.0_dp*F*COS(theta/2.0_dp)*SIN(theta/2.0_dp) & 280 *EXP(F*(COS(theta/2.0_dp))**2) 281 282 dri = d_expterm*dcos_thetahalf*dcosdri 283 284 drj = d_expterm*dcos_thetahalf*dcosdrj 285 286 drk = d_expterm*dcos_thetahalf*dcosdrk 287 288 f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + prefactor*dri(1) 289 f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + prefactor*dri(2) 290 f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + prefactor*dri(3) 291 292 f_nonbond(1, jparticle) = f_nonbond(1, jparticle) + prefactor*drj(1) 293 f_nonbond(2, jparticle) = f_nonbond(2, jparticle) + prefactor*drj(2) 294 f_nonbond(3, jparticle) = f_nonbond(3, jparticle) + prefactor*drj(3) 295 296 f_nonbond(1, kparticle) = f_nonbond(1, kparticle) + prefactor*drk(1) 297 f_nonbond(2, kparticle) = f_nonbond(2, kparticle) + prefactor*drk(2) 298 f_nonbond(3, kparticle) = f_nonbond(3, kparticle) + prefactor*drk(3) 299 300 IF (use_virial) THEN 301 CALL cp_abort(__LOCATION__, & 302 "using virial with Siepmann-Sprik"// & 303 " not implemented") 304 END IF 305 END DO 306 END SUBROUTINE siep_a_ij_d 307! ************************************************************************************************** 308!> \brief Phi_ij corresponds to Equation (3) in J. Chem. Phys., Vol. 102, p.511 309!> \param siepmann ... 310!> \param r_last_update_pbc ... 311!> \param iparticle ... 312!> \param jparticle ... 313!> \param cell_v ... 314!> \param cell ... 315!> \param rcutsq ... 316!> \param particle_set ... 317!> \param nr_oh ... 318!> \param nr_h3o ... 319!> \param nr_o ... 320!> \return ... 321!> \par History 322!> Using a local list of neighbors 323! ************************************************************************************************** 324 FUNCTION siep_Phi_ij(siepmann, r_last_update_pbc, iparticle, jparticle, & 325 cell_v, cell, rcutsq, particle_set, nr_oh, nr_h3o, nr_o) 326 TYPE(siepmann_pot_type), POINTER :: siepmann 327 TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc 328 INTEGER, INTENT(IN) :: iparticle, jparticle 329 REAL(KIND=dp), DIMENSION(3) :: cell_v 330 TYPE(cell_type), POINTER :: cell 331 REAL(KIND=dp), INTENT(IN) :: rcutsq 332 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 333 INTEGER, INTENT(INOUT), OPTIONAL :: nr_oh, nr_h3o, nr_o 334 REAL(KIND=dp) :: siep_Phi_ij 335 336 CHARACTER(LEN=*), PARAMETER :: routineN = 'siep_Phi_ij', routineP = moduleN//':'//routineN 337 338 CHARACTER(LEN=2) :: element_symbol 339 INTEGER :: count_h, iatom, index_h1, index_h2, natom 340 REAL(KIND=dp) :: cosphi, drih, drix, drji, h_max_dist, & 341 rab2_max, rih(3), rih1(3), rih2(3), & 342 rix(3), rji(3) 343 344 siep_Phi_ij = 0.0_dp 345 count_h = 0 346 index_h1 = 0 347 index_h2 = 0 348 rab2_max = rcutsq 349 h_max_dist = 2.27_dp ! 1.2 angstrom 350 natom = SIZE(particle_set) 351 CALL get_atomic_kind(atomic_kind=particle_set(iparticle)%atomic_kind, & 352 element_symbol=element_symbol) 353 IF (element_symbol /= "O") RETURN 354 rji(:) = -1.0_dp*(r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v) 355 drji = SQRT(DOT_PRODUCT(rji, rji)) 356 357 DO iatom = 1, natom 358 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, & 359 element_symbol=element_symbol) 360 IF (element_symbol /= "H") CYCLE 361 rih(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(iatom)%r(:), cell) 362 drih = SQRT(DOT_PRODUCT(rih, rih)) 363 IF (drih >= h_max_dist) CYCLE 364 count_h = count_h + 1 365 IF (count_h == 1) THEN 366 index_h1 = iatom 367 ELSEIF (count_h == 2) THEN 368 index_h2 = iatom 369 ENDIF 370 ENDDO 371 372 IF (count_h == 0) THEN 373 IF (siepmann%allow_o_formation) THEN 374 IF (PRESENT(nr_o)) nr_o = nr_o + 1 375 siep_Phi_ij = 0.0_dp 376 ELSE 377 CPABORT("No H atoms for O found") 378 ENDIF 379 ELSEIF (count_h == 1) THEN 380 IF (siepmann%allow_oh_formation) THEN 381 IF (PRESENT(nr_oh)) nr_oh = nr_oh + 1 382 siep_Phi_ij = 0.0_dp 383 ELSE 384 CPABORT("Only one H atom of O atom found") 385 ENDIF 386 ELSEIF (count_h == 3) THEN 387 IF (siepmann%allow_h3o_formation) THEN 388 IF (PRESENT(nr_h3o)) nr_h3o = nr_h3o + 1 389 siep_Phi_ij = 0.0_dp 390 ELSE 391 CPABORT("Three H atoms for O atom found") 392 ENDIF 393 ELSEIF (count_h > 3) THEN 394 CPABORT("Error in Siepmann-Sprik part: too many H atoms for O") 395 ENDIF 396 397 IF (count_h == 2) THEN 398 !dipole vector rix of the H2O molecule 399 rih1(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(index_h1)%r(:), cell) 400 rih2(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(index_h2)%r(:), cell) 401 rix(:) = rih1(:) + rih2(:) 402 drix = SQRT(DOT_PRODUCT(rix, rix)) 403 cosphi = DOT_PRODUCT(rji, rix)/(drji*drix) 404 IF (cosphi < -1.0_dp) cosphi = -1.0_dp 405 IF (cosphi > +1.0_dp) cosphi = +1.0_dp 406 siep_Phi_ij = EXP(-8.0_dp*((cosphi - 1)/4.0_dp)**4) 407 ENDIF 408 END FUNCTION siep_Phi_ij 409 410! ************************************************************************************************** 411!> \brief derivative of Phi_ij 412!> \param siepmann ... 413!> \param r_last_update_pbc ... 414!> \param iparticle ... 415!> \param jparticle ... 416!> \param f_nonbond ... 417!> \param prefactor ... 418!> \param cell_v ... 419!> \param cell ... 420!> \param rcutsq ... 421!> \param use_virial ... 422!> \param particle_set ... 423!> \par History 424!> Using a local list of neighbors 425! ************************************************************************************************** 426 SUBROUTINE siep_Phi_ij_d(siepmann, r_last_update_pbc, iparticle, jparticle, f_nonbond, & 427 prefactor, cell_v, cell, rcutsq, use_virial, particle_set) 428 TYPE(siepmann_pot_type), POINTER :: siepmann 429 TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc 430 INTEGER, INTENT(IN) :: iparticle, jparticle 431 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: f_nonbond 432 REAL(KIND=dp), INTENT(IN) :: prefactor 433 REAL(KIND=dp), DIMENSION(3) :: cell_v 434 TYPE(cell_type), POINTER :: cell 435 REAL(KIND=dp), INTENT(IN) :: rcutsq 436 LOGICAL, INTENT(IN) :: use_virial 437 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 438 439 CHARACTER(LEN=*), PARAMETER :: routineN = 'siep_Phi_ij_d', routineP = moduleN//':'//routineN 440 441 CHARACTER(LEN=2) :: element_symbol 442 INTEGER :: count_h, iatom, index_h1, index_h2, natom 443 REAL(KIND=dp) :: cosphi, dphi, drih, drix, drji, & 444 h_max_dist, Phi_ij, rab2_max 445 REAL(KIND=dp), DIMENSION(3) :: dcosdrh, dcosdri, dcosdrj, drh, dri, & 446 drj, rih, rih1, rih2, rix, rix_hat, & 447 rji, rji_hat 448 449 count_h = 0 450 index_h1 = 0 451 index_h2 = 0 452 rab2_max = rcutsq 453 h_max_dist = 2.27_dp ! 1.2 angstrom 454 natom = SIZE(particle_set) 455 Phi_ij = siep_Phi_ij(siepmann, r_last_update_pbc, iparticle, jparticle, & 456 cell_v, cell, rcutsq, & 457 particle_set) 458 rji(:) = -1.0_dp*(r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v) 459 drji = SQRT(DOT_PRODUCT(rji, rji)) 460 rji_hat(:) = rji(:)/drji 461 462 DO iatom = 1, natom 463 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, & 464 element_symbol=element_symbol) 465 IF (element_symbol /= "H") CYCLE 466 rih(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(iatom)%r(:), cell) 467 drih = SQRT(DOT_PRODUCT(rih, rih)) 468 IF (drih >= h_max_dist) CYCLE 469 count_h = count_h + 1 470 IF (count_h == 1) THEN 471 index_h1 = iatom 472 ELSEIF (count_h == 2) THEN 473 index_h2 = iatom 474 ENDIF 475 ENDDO 476 477 IF (count_h == 0 .AND. .NOT. siepmann%allow_o_formation) THEN 478 CPABORT("No H atoms for O found") 479 ELSEIF (count_h == 1 .AND. .NOT. siepmann%allow_oh_formation) THEN 480 CPABORT("Only one H atom for O atom found") 481 ELSEIF (count_h == 3 .AND. .NOT. siepmann%allow_h3o_formation) THEN 482 CPABORT("Three H atoms for O atom found") 483 ELSEIF (count_h > 3) THEN 484 CPABORT("Error in Siepmann-Sprik part: too many H atoms for O") 485 ENDIF 486 487 IF (count_h == 2) THEN 488 !dipole vector rix of the H2O molecule 489 rih1(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(index_h1)%r(:), cell) 490 rih2(:) = pbc(r_last_update_pbc(iparticle)%r(:), r_last_update_pbc(index_h2)%r(:), cell) 491 rix(:) = rih1(:) + rih2(:) 492 drix = SQRT(DOT_PRODUCT(rix, rix)) 493 rix_hat(:) = rix(:)/drix 494 cosphi = DOT_PRODUCT(rji, rix)/(drji*drix) 495 IF (cosphi < -1.0_dp) cosphi = -1.0_dp 496 IF (cosphi > +1.0_dp) cosphi = +1.0_dp 497 498 dcosdrj(:) = (1.0_dp/(drji))*(-rix_hat(:) + cosphi*rji_hat(:)) 499 ! for H atoms: 500 dcosdrh(:) = (1.0_dp/(drix))*(rji_hat(:) - cosphi*rix_hat(:)) 501 dcosdri(:) = -dcosdrj - 2.0_dp*dcosdrh 502 503 dphi = Phi_ij*(-8.0_dp)*((cosphi - 1)/4.0_dp)**3 504 505 dri = dphi*dcosdri 506 drj = dphi*dcosdrj 507 drh = dphi*dcosdrh 508 509 f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + prefactor*dri(1) 510 f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + prefactor*dri(2) 511 f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + prefactor*dri(3) 512 513 f_nonbond(1, jparticle) = f_nonbond(1, jparticle) + prefactor*drj(1) 514 f_nonbond(2, jparticle) = f_nonbond(2, jparticle) + prefactor*drj(2) 515 f_nonbond(3, jparticle) = f_nonbond(3, jparticle) + prefactor*drj(3) 516 517 f_nonbond(1, index_h1) = f_nonbond(1, index_h1) + prefactor*drh(1) 518 f_nonbond(2, index_h1) = f_nonbond(2, index_h1) + prefactor*drh(2) 519 f_nonbond(3, index_h1) = f_nonbond(3, index_h1) + prefactor*drh(3) 520 521 f_nonbond(1, index_h2) = f_nonbond(1, index_h2) + prefactor*drh(1) 522 f_nonbond(2, index_h2) = f_nonbond(2, index_h2) + prefactor*drh(2) 523 f_nonbond(3, index_h2) = f_nonbond(3, index_h2) + prefactor*drh(3) 524 525 IF (use_virial) THEN 526 CALL cp_abort(__LOCATION__, & 527 "using virial with Siepmann-Sprik"// & 528 " not implemented") 529 END IF 530 531 ENDIF 532 END SUBROUTINE siep_Phi_ij_d 533 534! ************************************************************************************************** 535!> \brief forces generated by the three-body term 536!> \param siepmann ... 537!> \param r_last_update_pbc ... 538!> \param cell_v ... 539!> \param n_loc_size ... 540!> \param full_loc_list ... 541!> \param iparticle ... 542!> \param jparticle ... 543!> \param f_nonbond ... 544!> \param use_virial ... 545!> \param rcutsq ... 546!> \param cell ... 547!> \param particle_set ... 548!> \par History 549!> Using a local list of neighbors 550! ************************************************************************************************** 551 SUBROUTINE siepmann_forces_v3(siepmann, r_last_update_pbc, cell_v, n_loc_size, & 552 full_loc_list, iparticle, jparticle, f_nonbond, & 553 use_virial, rcutsq, cell, particle_set) 554 TYPE(siepmann_pot_type), POINTER :: siepmann 555 TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc 556 REAL(KIND=dp), DIMENSION(3) :: cell_v 557 INTEGER, INTENT(IN) :: n_loc_size 558 INTEGER, DIMENSION(2, 1:n_loc_size) :: full_loc_list 559 INTEGER, INTENT(IN) :: iparticle, jparticle 560 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: f_nonbond 561 LOGICAL, INTENT(IN) :: use_virial 562 REAL(KIND=dp), INTENT(IN) :: rcutsq 563 TYPE(cell_type), POINTER :: cell 564 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 565 566 CHARACTER(LEN=*), PARAMETER :: routineN = 'siepmann_forces_v3', & 567 routineP = moduleN//':'//routineN 568 569 CHARACTER(LEN=2) :: element_symbol 570 REAL(KIND=dp) :: a_ij, beta, drji, E, f2, f2_d, f_A1, & 571 f_A2, fac, prefactor, rji(3), & 572 rji_hat(3) 573 574 beta = siepmann%beta 575 E = siepmann%E 576 577 CALL get_atomic_kind(atomic_kind=particle_set(iparticle)%atomic_kind, & 578 element_symbol=element_symbol) 579 IF (element_symbol /= "O") RETURN 580 581 rji(:) = -1.0_dp*(r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v) 582 drji = SQRT(DOT_PRODUCT(rji, rji)) 583 rji_hat(:) = rji(:)/drji 584 585 fac = -1.0_dp !gradient to force 586 a_ij = siep_a_ij(siepmann, r_last_update_pbc, iparticle, jparticle, n_loc_size, & 587 full_loc_list, cell_v, rcutsq, particle_set, cell) 588 f2 = siep_f2(siepmann, drji) 589 f2_d = siep_f2_d(siepmann, drji) 590 591 ! Lets do the f_A1 piece derivative of f2 592 f_A1 = E*f2_d*drji**(-beta)*a_ij*fac*(1.0_dp/drji) 593 f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_A1*rji(1) 594 f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_A1*rji(2) 595 f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_A1*rji(3) 596 f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_A1*rji(1) 597 f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_A1*rji(2) 598 f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_A1*rji(3) 599 600 IF (use_virial) THEN 601 CALL cp_abort(__LOCATION__, & 602 "using virial with Siepmann-Sprik"// & 603 " not implemented") 604 END IF 605 606 ! Lets do the f_A2 piece derivative of rji**(-beta) 607 f_A2 = E*f2*(-beta)*drji**(-beta - 1)*a_ij*fac*(1.0_dp/drji) 608 f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_A2*rji(1) 609 f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_A2*rji(2) 610 f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_A2*rji(3) 611 f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_A2*rji(1) 612 f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_A2*rji(2) 613 f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_A2*rji(3) 614 615 IF (use_virial) THEN 616 CALL cp_abort(__LOCATION__, & 617 "using virial with Siepmann-Sprik"// & 618 " not implemented") 619 END IF 620 621 ! Lets do the f_A3 piece derivative: of a_ij 622 prefactor = E*f2*drji**(-beta)*fac 623 CALL siep_a_ij_d(siepmann, r_last_update_pbc, iparticle, jparticle, f_nonbond, & 624 prefactor, n_loc_size, full_loc_list, cell_v, & 625 cell, rcutsq, use_virial) 626 627 END SUBROUTINE siepmann_forces_v3 628 629! ************************************************************************************************** 630!> \brief forces generated by the dipole term 631!> \param siepmann ... 632!> \param r_last_update_pbc ... 633!> \param cell_v ... 634!> \param cell ... 635!> \param iparticle ... 636!> \param jparticle ... 637!> \param f_nonbond ... 638!> \param use_virial ... 639!> \param rcutsq ... 640!> \param particle_set ... 641!> \par History 642!> Using a local list of neighbors 643! ************************************************************************************************** 644 SUBROUTINE siepmann_forces_v2(siepmann, r_last_update_pbc, cell_v, cell, & 645 iparticle, jparticle, f_nonbond, use_virial, rcutsq, particle_set) 646 TYPE(siepmann_pot_type), POINTER :: siepmann 647 TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc 648 REAL(KIND=dp), DIMENSION(3) :: cell_v 649 TYPE(cell_type), POINTER :: cell 650 INTEGER, INTENT(IN) :: iparticle, jparticle 651 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: f_nonbond 652 LOGICAL, INTENT(IN) :: use_virial 653 REAL(KIND=dp), INTENT(IN) :: rcutsq 654 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 655 656 CHARACTER(LEN=*), PARAMETER :: routineN = 'siepmann_forces_v2', & 657 routineP = moduleN//':'//routineN 658 659 CHARACTER(LEN=2) :: element_symbol 660 REAL(KIND=dp) :: D, drji, f2, f2_d, f_A1, f_A2, fac, & 661 Phi_ij, prefactor, rji(3) 662 663 D = siepmann%D 664 665 CALL get_atomic_kind(atomic_kind=particle_set(iparticle)%atomic_kind, & 666 element_symbol=element_symbol) 667 IF (element_symbol /= "O") RETURN 668 669 rji(:) = -1.0_dp*(r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v) 670 drji = SQRT(DOT_PRODUCT(rji, rji)) 671 672 fac = -1.0_dp 673 Phi_ij = siep_Phi_ij(siepmann, r_last_update_pbc, iparticle, jparticle, & 674 cell_v, cell, rcutsq, particle_set) 675 f2 = siep_f2(siepmann, drji) 676 f2_d = siep_f2_d(siepmann, drji) 677 678 ! Lets do the f_A1 piece derivative of f2 679 f_A1 = -D*f2_d*drji**(-3)*Phi_ij*fac*(1.0_dp/drji) 680 f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_A1*rji(1) 681 f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_A1*rji(2) 682 f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_A1*rji(3) 683 f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_A1*rji(1) 684 f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_A1*rji(2) 685 f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_A1*rji(3) 686 687 IF (use_virial) THEN 688 CALL cp_abort(__LOCATION__, & 689 "using virial with Siepmann-Sprik"// & 690 " not implemented") 691 END IF 692 693! ! Lets do the f_A2 piece derivative of rji**(-3) 694 f_A2 = -D*f2*(-3.0_dp)*drji**(-4)*Phi_ij*fac*(1.0_dp/drji) 695 f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_A2*rji(1) 696 f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_A2*rji(2) 697 f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_A2*rji(3) 698 f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_A2*rji(1) 699 f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_A2*rji(2) 700 f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_A2*rji(3) 701 702 IF (use_virial) THEN 703 CALL cp_abort(__LOCATION__, & 704 "using virial with Siepmann-Sprik"// & 705 " not implemented") 706 END IF 707 708 ! Lets do the f_A3 piece derivative: of Phi_ij 709 prefactor = -D*f2*drji**(-3)*fac 710 CALL siep_Phi_ij_d(siepmann, r_last_update_pbc, iparticle, jparticle, f_nonbond, & 711 prefactor, cell_v, cell, & 712 rcutsq, use_virial, particle_set) 713 714 END SUBROUTINE siepmann_forces_v2 715 716! ************************************************************************************************** 717!> \brief ... 718!> \param nonbonded ... 719!> \param potparm ... 720!> \param glob_loc_list ... 721!> \param glob_cell_v ... 722!> \param glob_loc_list_a ... 723!> \param cell ... 724!> \par History 725! ************************************************************************************************** 726 SUBROUTINE setup_siepmann_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, & 727 glob_loc_list_a, cell) 728 TYPE(fist_neighbor_type), POINTER :: nonbonded 729 TYPE(pair_potential_pp_type), POINTER :: potparm 730 INTEGER, DIMENSION(:, :), POINTER :: glob_loc_list 731 REAL(KIND=dp), DIMENSION(:, :), POINTER :: glob_cell_v 732 INTEGER, DIMENSION(:), POINTER :: glob_loc_list_a 733 TYPE(cell_type), POINTER :: cell 734 735 CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_siepmann_arrays', & 736 routineP = moduleN//':'//routineN 737 738 INTEGER :: handle, i, iend, igrp, ikind, ilist, & 739 ipair, istart, jkind, nkinds, npairs, & 740 npairs_tot 741 INTEGER, DIMENSION(:), POINTER :: work_list, work_list2 742 INTEGER, DIMENSION(:, :), POINTER :: list 743 REAL(KIND=dp), DIMENSION(3) :: cell_v, cvi 744 REAL(KIND=dp), DIMENSION(:, :), POINTER :: rwork_list 745 TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair 746 TYPE(pair_potential_single_type), POINTER :: pot 747 748 CPASSERT(.NOT. ASSOCIATED(glob_loc_list)) 749 CPASSERT(.NOT. ASSOCIATED(glob_loc_list_a)) 750 CPASSERT(.NOT. ASSOCIATED(glob_cell_v)) 751 CALL timeset(routineN, handle) 752 npairs_tot = 0 753 nkinds = SIZE(potparm%pot, 1) 754 DO ilist = 1, nonbonded%nlists 755 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist) 756 npairs = neighbor_kind_pair%npairs 757 IF (npairs == 0) CYCLE 758 Kind_Group_Loop1: DO igrp = 1, neighbor_kind_pair%ngrp_kind 759 istart = neighbor_kind_pair%grp_kind_start(igrp) 760 iend = neighbor_kind_pair%grp_kind_end(igrp) 761 ikind = neighbor_kind_pair%ij_kind(1, igrp) 762 jkind = neighbor_kind_pair%ij_kind(2, igrp) 763 pot => potparm%pot(ikind, jkind)%pot 764 npairs = iend - istart + 1 765 IF (pot%no_mb) CYCLE 766 DO i = 1, SIZE(pot%type) 767 IF (pot%type(i) == siepmann_type) npairs_tot = npairs_tot + npairs 768 END DO 769 END DO Kind_Group_Loop1 770 END DO 771 ALLOCATE (work_list(npairs_tot)) 772 ALLOCATE (work_list2(npairs_tot)) 773 ALLOCATE (glob_loc_list(2, npairs_tot)) 774 ALLOCATE (glob_cell_v(3, npairs_tot)) 775 ! Fill arrays with data 776 npairs_tot = 0 777 DO ilist = 1, nonbonded%nlists 778 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist) 779 npairs = neighbor_kind_pair%npairs 780 IF (npairs == 0) CYCLE 781 Kind_Group_Loop2: DO igrp = 1, neighbor_kind_pair%ngrp_kind 782 istart = neighbor_kind_pair%grp_kind_start(igrp) 783 iend = neighbor_kind_pair%grp_kind_end(igrp) 784 ikind = neighbor_kind_pair%ij_kind(1, igrp) 785 jkind = neighbor_kind_pair%ij_kind(2, igrp) 786 list => neighbor_kind_pair%list 787 cvi = neighbor_kind_pair%cell_vector 788 pot => potparm%pot(ikind, jkind)%pot 789 npairs = iend - istart + 1 790 IF (pot%no_mb) CYCLE 791 CALL matvec_3x3(cell_v, cell%hmat, cvi) 792 DO i = 1, SIZE(pot%type) 793 ! SIEPMANN 794 IF (pot%type(i) == siepmann_type) THEN 795 DO ipair = 1, npairs 796 glob_loc_list(:, npairs_tot + ipair) = list(:, istart - 1 + ipair) 797 glob_cell_v(1:3, npairs_tot + ipair) = cell_v(1:3) 798 END DO 799 npairs_tot = npairs_tot + npairs 800 END IF 801 END DO 802 END DO Kind_Group_Loop2 803 END DO 804 ! Order the arrays w.r.t. the first index of glob_loc_list 805 CALL sort(glob_loc_list(1, :), npairs_tot, work_list) 806 DO ipair = 1, npairs_tot 807 work_list2(ipair) = glob_loc_list(2, work_list(ipair)) 808 END DO 809 glob_loc_list(2, :) = work_list2 810 DEALLOCATE (work_list2) 811 ALLOCATE (rwork_list(3, npairs_tot)) 812 DO ipair = 1, npairs_tot 813 rwork_list(:, ipair) = glob_cell_v(:, work_list(ipair)) 814 END DO 815 glob_cell_v = rwork_list 816 DEALLOCATE (rwork_list) 817 DEALLOCATE (work_list) 818 ALLOCATE (glob_loc_list_a(npairs_tot)) 819 glob_loc_list_a = glob_loc_list(1, :) 820 CALL timestop(handle) 821 END SUBROUTINE setup_siepmann_arrays 822 823! ************************************************************************************************** 824!> \brief ... 825!> \param glob_loc_list ... 826!> \param glob_cell_v ... 827!> \param glob_loc_list_a ... 828! ************************************************************************************************** 829 SUBROUTINE destroy_siepmann_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a) 830 INTEGER, DIMENSION(:, :), POINTER :: glob_loc_list 831 REAL(KIND=dp), DIMENSION(:, :), POINTER :: glob_cell_v 832 INTEGER, DIMENSION(:), POINTER :: glob_loc_list_a 833 834 CHARACTER(LEN=*), PARAMETER :: routineN = 'destroy_siepmann_arrays', & 835 routineP = moduleN//':'//routineN 836 837 IF (ASSOCIATED(glob_loc_list)) THEN 838 DEALLOCATE (glob_loc_list) 839 END IF 840 IF (ASSOCIATED(glob_loc_list_a)) THEN 841 DEALLOCATE (glob_loc_list_a) 842 END IF 843 IF (ASSOCIATED(glob_cell_v)) THEN 844 DEALLOCATE (glob_cell_v) 845 END IF 846 847 END SUBROUTINE destroy_siepmann_arrays 848 849! ************************************************************************************************** 850!> \brief prints the number of OH- ions or H3O+ ions near surface 851!> \param nr_ions number of ions 852!> \param mm_section ... 853!> \param para_env ... 854!> \param print_oh flag indicating if number OH- is printed 855!> \param print_h3o flag indicating if number H3O+ is printed 856!> \param print_o flag indicating if number O^(2-) is printed 857! ************************************************************************************************** 858 SUBROUTINE print_nr_ions_siepmann(nr_ions, mm_section, para_env, print_oh, & 859 print_h3o, print_o) 860 INTEGER, INTENT(INOUT) :: nr_ions 861 TYPE(section_vals_type), POINTER :: mm_section 862 TYPE(cp_para_env_type), OPTIONAL, POINTER :: para_env 863 LOGICAL, INTENT(IN) :: print_oh, print_h3o, print_o 864 865 CHARACTER(LEN=*), PARAMETER :: routineN = 'print_nr_ions_siepmann', & 866 routineP = moduleN//':'//routineN 867 868 INTEGER :: iw 869 TYPE(cp_logger_type), POINTER :: logger 870 871 NULLIFY (logger) 872 873 CALL mp_sum(nr_ions, para_env%group) 874 logger => cp_get_default_logger() 875 876 iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%PROGRAM_RUN_INFO", & 877 extension=".mmLog") 878 879 IF (iw > 0 .AND. nr_ions > 0 .AND. print_oh) THEN 880 WRITE (iw, '(/,A,T71,I10,/)') " SIEPMANN: number of OH- ions at surface", nr_ions 881 ENDIF 882 IF (iw > 0 .AND. nr_ions > 0 .AND. print_h3o) THEN 883 WRITE (iw, '(/,A,T71,I10,/)') " SIEPMANN: number of H3O+ ions at surface", nr_ions 884 ENDIF 885 IF (iw > 0 .AND. nr_ions > 0 .AND. print_o) THEN 886 WRITE (iw, '(/,A,T71,I10,/)') " SIEPMANN: number of O^2- ions at surface", nr_ions 887 ENDIF 888 889 CALL cp_print_key_finished_output(iw, logger, mm_section, "PRINT%PROGRAM_RUN_INFO") 890 891 END SUBROUTINE print_nr_ions_siepmann 892 893END MODULE manybody_siepmann 894 895