1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7MODULE qmmm_per_elpot 8 9! ************************************************************************************************** 10!> \brief Setting up the potential for QM/MM periodic boundary conditions calculations 11!> \par History 12!> 07.2005 created [tlaino] 13!> \author Teodoro Laino 14 USE ao_util, ONLY: exp_radius 15 USE cell_types, ONLY: cell_type 16 USE cp_log_handling, ONLY: cp_get_default_logger,& 17 cp_logger_get_default_io_unit,& 18 cp_logger_type 19 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 20 cp_print_key_unit_nr 21 USE cp_para_types, ONLY: cp_para_env_type 22 USE ewald_environment_types, ONLY: ewald_env_create,& 23 ewald_env_get,& 24 ewald_env_set,& 25 ewald_environment_type,& 26 read_ewald_section 27 USE ewald_pw_types, ONLY: ewald_pw_create,& 28 ewald_pw_type 29 USE ewald_spline_util, ONLY: Setup_Ewald_Spline 30 USE input_constants, ONLY: do_qmmm_coulomb,& 31 do_qmmm_gauss,& 32 do_qmmm_pcharge,& 33 do_qmmm_swave 34 USE input_section_types, ONLY: section_vals_get_subs_vals,& 35 section_vals_type,& 36 section_vals_val_get 37 USE kinds, ONLY: dp 38 USE mathconstants, ONLY: pi 39 USE pw_poisson_types, ONLY: do_ewald_ewald,& 40 do_ewald_none,& 41 do_ewald_pme,& 42 do_ewald_spme 43 USE qmmm_gaussian_types, ONLY: qmmm_gaussian_p_type,& 44 qmmm_gaussian_type 45 USE qmmm_types_low, ONLY: qmmm_per_pot_p_type,& 46 qmmm_per_pot_type,& 47 qmmm_pot_p_type 48#include "./base/base_uses.f90" 49 50 IMPLICIT NONE 51 PRIVATE 52 53 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE. 54 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_per_elpot' 55 PUBLIC :: qmmm_per_potential_init, qmmm_ewald_potential_init 56 57CONTAINS 58 59! ************************************************************************************************** 60!> \brief Initialize the QMMM potential stored on vector, 61!> according the qmmm_coupl_type 62!> \param qmmm_coupl_type ... 63!> \param per_potentials ... 64!> \param potentials ... 65!> \param pgfs ... 66!> \param qm_cell_small ... 67!> \param mm_cell ... 68!> \param para_env ... 69!> \param compatibility ... 70!> \param qmmm_periodic ... 71!> \param print_section ... 72!> \param eps_mm_rspace ... 73!> \param maxchrg ... 74!> \param ncp ... 75!> \param ncpl ... 76!> \par History 77!> 09.2004 created [tlaino] 78!> \author Teodoro Laino 79! ************************************************************************************************** 80 SUBROUTINE qmmm_per_potential_init(qmmm_coupl_type, per_potentials, potentials, & 81 pgfs, qm_cell_small, mm_cell, para_env, compatibility, qmmm_periodic, print_section, & 82 eps_mm_rspace, maxchrg, ncp, ncpl) 83 INTEGER, INTENT(IN) :: qmmm_coupl_type 84 TYPE(qmmm_per_pot_p_type), DIMENSION(:), POINTER :: per_potentials 85 TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: potentials 86 TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs 87 TYPE(cell_type), POINTER :: qm_cell_small, mm_cell 88 TYPE(cp_para_env_type), POINTER :: para_env 89 LOGICAL, INTENT(IN) :: compatibility 90 TYPE(section_vals_type), POINTER :: qmmm_periodic, print_section 91 REAL(KIND=dp), INTENT(IN) :: eps_mm_rspace, maxchrg 92 INTEGER, INTENT(IN) :: ncp(3), ncpl(3) 93 94 CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_per_potential_init', & 95 routineP = moduleN//':'//routineN 96 97 INTEGER :: I, idim, ig, ig_start, iw, ix, iy, iz, & 98 K, Kmax(3), n_rep_real(3), & 99 n_rep_real_val, ncoarsel, ncoarset, & 100 Ndim, output_unit 101 INTEGER, DIMENSION(:), POINTER :: mm_atom_index 102 REAL(KIND=dp) :: Ak, alpha, box(3), Fac(3), fs, g, g2, & 103 Gk, Gmax, mymaxradius, npl, npt, & 104 Prefactor, rc, rc2, Rmax, tmp, vec(3), & 105 vol 106 REAL(KIND=dp), DIMENSION(:), POINTER :: gx, gy, gz, Lg 107 TYPE(cp_logger_type), POINTER :: logger 108 TYPE(qmmm_gaussian_type), POINTER :: pgf 109 110 NULLIFY (Lg, gx, gy, gz) 111 ncoarset = PRODUCT(ncp) 112 ncoarsel = PRODUCT(ncpl) 113 logger => cp_get_default_logger() 114 output_unit = cp_logger_get_default_io_unit(logger) 115 Rmax = SQRT(mm_cell%hmat(1, 1)**2 + & 116 mm_cell%hmat(2, 2)**2 + & 117 mm_cell%hmat(3, 3)**2) 118 CALL section_vals_val_get(qmmm_periodic, "GMAX", r_val=Gmax) 119 CALL section_vals_val_get(qmmm_periodic, "REPLICA", i_val=n_rep_real_val) 120 fac = 2.0e0_dp*Pi/(/mm_cell%hmat(1, 1), mm_cell%hmat(2, 2), mm_cell%hmat(3, 3)/) 121 Kmax = CEILING(Gmax/Fac) 122 Vol = mm_cell%hmat(1, 1)* & 123 mm_cell%hmat(2, 2)* & 124 mm_cell%hmat(3, 3) 125 Ndim = (Kmax(1) + 1)*(2*Kmax(2) + 1)*(2*Kmax(3) + 1) 126 ig_start = 1 127 n_rep_real = n_rep_real_val 128 IF (compatibility .AND. (qmmm_coupl_type == do_qmmm_gauss)) ig_start = 2 129 130 CPASSERT(.NOT. ASSOCIATED(per_potentials)) 131 ALLOCATE (per_potentials(SIZE(pgfs))) 132 CPASSERT(SIZE(pgfs) == SIZE(potentials)) 133 Potential_Type: DO K = 1, SIZE(pgfs) 134 135 rc = pgfs(K)%pgf%Elp_Radius 136 ALLOCATE (per_potentials(K)%Pot) 137 SELECT CASE (qmmm_coupl_type) 138 CASE (do_qmmm_coulomb, do_qmmm_pcharge) 139 ! Not yet implemented for this case 140 CPABORT("") 141 CASE (do_qmmm_gauss, do_qmmm_swave) 142 ALLOCATE (Lg(Ndim)) 143 ALLOCATE (gx(Ndim)) 144 ALLOCATE (gy(Ndim)) 145 ALLOCATE (gz(Ndim)) 146 END SELECT 147 148 LG = 0.0_dp 149 gx = 0.0_dp 150 gy = 0.0_dp 151 gz = 0.0_dp 152 153 SELECT CASE (qmmm_coupl_type) 154 CASE (do_qmmm_coulomb, do_qmmm_pcharge) 155 ! Not yet implemented for this case 156 CPABORT("") 157 CASE (do_qmmm_gauss, do_qmmm_swave) 158 pgf => pgfs(K)%pgf 159 idim = 0 160 DO ix = 0, kmax(1) 161 DO iy = -kmax(2), kmax(2) 162 DO iz = -kmax(3), kmax(3) 163 idim = idim + 1 164 IF (ix == 0 .AND. iy == 0 .AND. iz == 0) THEN 165 DO Ig = ig_start, pgf%number_of_gaussians 166 Gk = pgf%Gk(Ig) 167 Ak = pgf%Ak(Ig)*Pi**(3.0_dp/2.0_dp)*Gk**3.0_dp 168 LG(idim) = LG(idim) - Ak 169 END DO 170 ELSE 171 fs = 2.0_dp; IF (ix == 0) fs = 1.0_dp 172 vec = fac*(/REAL(ix, KIND=dp), REAL(iy, KIND=dp), REAL(iz, KIND=dp)/) 173 g2 = DOT_PRODUCT(vec, vec) 174 rc2 = rc*rc 175 g = SQRT(g2) 176 IF (qmmm_coupl_type == do_qmmm_gauss) THEN 177 LG(idim) = 4.0_dp*Pi/g2*EXP(-(g2*rc2)/4.0_dp) 178 ELSEIF (qmmm_coupl_type == do_qmmm_swave) THEN 179 tmp = 4.0_dp/rc2 180 LG(idim) = 4.0_dp*Pi*tmp**2/(g2*(g2 + tmp)**2) 181 END IF 182 DO Ig = ig_start, pgf%number_of_gaussians 183 Gk = pgf%Gk(Ig) 184 Ak = pgf%Ak(Ig)*Pi**(3.0_dp/2.0_dp)*Gk**3.0_dp 185 LG(idim) = LG(idim) - Ak*EXP(-(g*Gk)**2.0_dp/4.0_dp) 186 END DO 187 ENDIF 188 LG(idim) = fs*LG(idim)*1.0_dp/Vol 189 gx(idim) = fac(1)*REAL(ix, KIND=dp) 190 gy(idim) = fac(2)*REAL(iy, KIND=dp) 191 gz(idim) = fac(3)*REAL(iz, KIND=dp) 192 END DO 193 END DO 194 END DO 195 196 IF (ALL(n_rep_real == -1)) THEN 197 mymaxradius = 0.0_dp 198 DO I = 1, pgf%number_of_gaussians 199 IF (pgf%Gk(I) /= 0.0_dp) THEN 200 alpha = 1.0_dp/pgf%Gk(I) 201 alpha = alpha*alpha 202 Prefactor = pgf%Ak(I)*maxchrg 203 mymaxradius = MAX(mymaxradius, exp_radius(0, alpha, eps_mm_rspace, Prefactor)) 204 END IF 205 END DO 206 box(1) = (qm_cell_small%hmat(1, 1) - mm_cell%hmat(1, 1))/2.0_dp 207 box(2) = (qm_cell_small%hmat(2, 2) - mm_cell%hmat(2, 2))/2.0_dp 208 box(3) = (qm_cell_small%hmat(3, 3) - mm_cell%hmat(3, 3))/2.0_dp 209 IF (ANY(box > 0.0_dp)) THEN 210 CPABORT("") 211 END IF 212 n_rep_real(1) = CEILING((box(1) + mymaxradius)/mm_cell%hmat(1, 1)) 213 n_rep_real(2) = CEILING((box(2) + mymaxradius)/mm_cell%hmat(2, 2)) 214 n_rep_real(3) = CEILING((box(3) + mymaxradius)/mm_cell%hmat(3, 3)) 215 END IF 216 217 CASE DEFAULT 218 DEALLOCATE (per_potentials(K)%Pot) 219 IF (output_unit > 0) WRITE (output_unit, '(A)') " QMMM Periodic Potential - not Initialized!" 220 CYCLE Potential_Type 221 END SELECT 222 223 NULLIFY (mm_atom_index) 224 ALLOCATE (mm_atom_index(SIZE(potentials(K)%pot%mm_atom_index))) 225 mm_atom_index = potentials(K)%pot%mm_atom_index 226 227 NULLIFY (per_potentials(K)%Pot%LG, per_potentials(K)%Pot%mm_atom_index, & 228 per_potentials(K)%Pot%gx, per_potentials(K)%Pot%gy, per_potentials(K)%Pot%gz) 229 CALL qmmm_per_pot_type_create(per_potentials(K)%Pot, LG=LG, gx=gx, gy=gy, gz=gz, & 230 Gmax=Gmax, Kmax=Kmax, n_rep_real=n_rep_real, & 231 Fac=Fac, mm_atom_index=mm_atom_index, & 232 mm_cell=mm_cell, para_env=para_env, & 233 qmmm_per_section=qmmm_periodic, print_section=print_section) 234 235 iw = cp_print_key_unit_nr(logger, print_section, "PERIODIC_INFO", & 236 extension=".log") 237 IF (iw > 0) THEN 238 npt = REAL(ncoarset, KIND=dp)*REAL(ndim, KIND=dp)*REAL(SIZE(mm_atom_index), KIND=dp) 239 npl = REAL(ncoarsel, KIND=dp)*REAL(ndim, KIND=dp)*REAL(SIZE(mm_atom_index), KIND=dp) 240 WRITE (UNIT=iw, FMT="(/,T2,A)") REPEAT("-", 79) 241 WRITE (UNIT=iw, FMT="(T2,A,T20,A,T80,A)") "-", "QMMM PERIODIC BOUNDARY CONDITION INFO", "-" 242 WRITE (UNIT=iw, FMT="(T2,A)") REPEAT("-", 79) 243 WRITE (UNIT=iw, FMT="(T2,A,T10,A,F15.6,T50,A,3I5,T80,A)") "-", "RADIUS =", rc, "REPLICA =", n_rep_real, "-" 244 WRITE (UNIT=iw, FMT="(T2,A,T10,A,F15.6,T50,A,I15,T80,A)") "-", "MINGVAL =", MINVAL(ABS(Lg)), & 245 "GPOINTS =", ndim, "-" 246 WRITE (UNIT=iw, FMT="(T2,A,T10,A,3I5,T50,A,3I5,T80,A)") "-", "NCOARSL =", ncpl, & 247 "NCOARST =", ncp, "-" 248 WRITE (UNIT=iw, FMT="(T2,A,T10,A,F15.0,T50,A,F15.0,T80,A)") "-", "NFLOP-L ~", npl, & 249 "NFLOP-T ~", npt, "-" 250 WRITE (UNIT=iw, FMT="(T2,A)") REPEAT("-", 79) 251 END IF 252 CALL cp_print_key_finished_output(iw, logger, print_section, & 253 "PERIODIC_INFO") 254 255 END DO Potential_Type 256 257 END SUBROUTINE qmmm_per_potential_init 258 259! ************************************************************************************************** 260!> \brief Creates the qmmm_pot_type structure 261!> \param Pot ... 262!> \param LG ... 263!> \param gx ... 264!> \param gy ... 265!> \param gz ... 266!> \param GMax ... 267!> \param Kmax ... 268!> \param n_rep_real ... 269!> \param Fac ... 270!> \param mm_atom_index ... 271!> \param mm_cell ... 272!> \param para_env ... 273!> \param qmmm_per_section ... 274!> \param print_section ... 275!> \par History 276!> 09.2004 created [tlaino] 277!> \author Teodoro Laino 278! ************************************************************************************************** 279 SUBROUTINE qmmm_per_pot_type_create(Pot, LG, gx, gy, gz, GMax, Kmax, n_rep_real, & 280 Fac, mm_atom_index, mm_cell, para_env, qmmm_per_section, print_section) 281 TYPE(qmmm_per_pot_type), POINTER :: Pot 282 REAL(KIND=dp), DIMENSION(:), POINTER :: LG, gx, gy, gz 283 REAL(KIND=dp), INTENT(IN) :: Gmax 284 INTEGER, INTENT(IN) :: Kmax(3), n_rep_real(3) 285 REAL(KIND=dp), INTENT(IN) :: Fac(3) 286 INTEGER, DIMENSION(:), POINTER :: mm_atom_index 287 TYPE(cell_type), POINTER :: mm_cell 288 TYPE(cp_para_env_type), POINTER :: para_env 289 TYPE(section_vals_type), POINTER :: qmmm_per_section, print_section 290 291 CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_per_pot_type_create', & 292 routineP = moduleN//':'//routineN 293 294 INTEGER :: npts(3) 295 INTEGER, DIMENSION(:), POINTER :: ngrids 296 REAL(KIND=dp) :: hmat(3, 3) 297 TYPE(section_vals_type), POINTER :: grid_print_section 298 299 Pot%LG => LG 300 Pot%gx => gx 301 Pot%gy => gy 302 Pot%gz => gz 303 Pot%mm_atom_index => mm_atom_index 304 Pot%Gmax = Gmax 305 Pot%Kmax = Kmax 306 Pot%n_rep_real = n_rep_real 307 Pot%Fac = Fac 308 ! 309 ! Setting Up Fit Procedure 310 ! 311 NULLIFY (Pot%pw_grid) 312 NULLIFY (Pot%pw_pool) 313 NULLIFY (Pot%TabLR, ngrids) 314 CALL section_vals_val_get(qmmm_per_section, "ngrids", i_vals=ngrids) 315 npts = ngrids 316 hmat = mm_cell%hmat 317 318 grid_print_section => section_vals_get_subs_vals(print_section, "GRID_INFORMATION") 319 CALL Setup_Ewald_Spline(pw_grid=Pot%pw_grid, pw_pool=Pot%pw_pool, coeff=Pot%TabLR, & 320 LG=LG, gx=gx, gy=gy, gz=gz, hmat=hmat, npts=npts, param_section=qmmm_per_section, & 321 tag="qmmm", para_env=para_env, print_section=grid_print_section) 322 323 END SUBROUTINE qmmm_per_pot_type_create 324 325! ************************************************************************************************** 326!> \brief Initialize the QMMM Ewald potential needed for QM-QM Coupling using 327!> point charges 328!> \param ewald_env ... 329!> \param ewald_pw ... 330!> \param qmmm_coupl_type ... 331!> \param mm_cell ... 332!> \param para_env ... 333!> \param qmmm_periodic ... 334!> \param print_section ... 335!> \par History 336!> 10.2014 created [JGH] 337!> \author JGH 338! ************************************************************************************************** 339 SUBROUTINE qmmm_ewald_potential_init(ewald_env, ewald_pw, qmmm_coupl_type, mm_cell, para_env, & 340 qmmm_periodic, print_section) 341 TYPE(ewald_environment_type), POINTER :: ewald_env 342 TYPE(ewald_pw_type), POINTER :: ewald_pw 343 INTEGER, INTENT(IN) :: qmmm_coupl_type 344 TYPE(cell_type), POINTER :: mm_cell 345 TYPE(cp_para_env_type), POINTER :: para_env 346 TYPE(section_vals_type), POINTER :: qmmm_periodic, print_section 347 348 CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_ewald_potential_init', & 349 routineP = moduleN//':'//routineN 350 351 INTEGER :: ewald_type, gmax(3), iw, o_spline, ounit 352 LOGICAL :: do_multipoles 353 REAL(KIND=dp) :: alpha, rcut 354 TYPE(cp_logger_type), POINTER :: logger 355 TYPE(section_vals_type), POINTER :: ewald_print_section, ewald_section, & 356 poisson_section 357 358 logger => cp_get_default_logger() 359 ounit = cp_logger_get_default_io_unit(logger) 360 CPASSERT(.NOT. ASSOCIATED(ewald_env)) 361 CPASSERT(.NOT. ASSOCIATED(ewald_pw)) 362 363 ! Create Ewald environments 364 poisson_section => section_vals_get_subs_vals(qmmm_periodic, "POISSON") 365 CALL ewald_env_create(ewald_env, para_env) 366 CALL ewald_env_set(ewald_env, poisson_section=poisson_section) 367 ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD") 368 CALL read_ewald_section(ewald_env, ewald_section) 369 ewald_print_section => section_vals_get_subs_vals(print_section, "GRID_INFORMATION") 370 CALL ewald_pw_create(ewald_pw, ewald_env, mm_cell, mm_cell, print_section=ewald_print_section) 371 372 CALL ewald_env_get(ewald_env, ewald_type=ewald_type, do_multipoles=do_multipoles, & 373 gmax=gmax, o_spline=o_spline, alpha=alpha, rcut=rcut) 374 IF (do_multipoles) & 375 CPABORT("No multipole force fields allowed in QM-QM Ewald long range correction") 376 377 SELECT CASE (qmmm_coupl_type) 378 CASE (do_qmmm_coulomb) 379 CPABORT("QM-QM long range correction not possible with COULOMB coupling") 380 CASE (do_qmmm_pcharge) 381 ! OK 382 CASE (do_qmmm_gauss, do_qmmm_swave) 383 CPABORT("QM-QM long range correction not possible with GAUSS/SWAVE coupling") 384 CASE DEFAULT 385 ! We should never get to this point 386 CPABORT("") 387 END SELECT 388 389 iw = cp_print_key_unit_nr(logger, print_section, "PERIODIC_INFO", extension=".log") 390 IF (iw > 0) THEN 391 WRITE (UNIT=iw, FMT="(/,T2,A)") REPEAT("-", 79) 392 WRITE (UNIT=iw, FMT="(T2,A,T20,A,T80,A)") "-", "QMMM PERIODIC BOUNDARY CONDITION INFO", "-" 393 WRITE (UNIT=iw, FMT="(T2,A,T25,A,T80,A)") "-", "QM-QM Long Range Correction", "-" 394 WRITE (UNIT=iw, FMT="(T2,A)") REPEAT("-", 79) 395 SELECT CASE (ewald_type) 396 CASE (do_ewald_none) 397 CPABORT("QM-QM long range correction not compatible with Ewald=NONE") 398 CASE (do_ewald_pme) 399 CPABORT("QM-QM long range correction not possible with Ewald=PME") 400 CASE (do_ewald_ewald) 401 CPABORT("QM-QM long range correction not possible with Ewald method") 402 CASE (do_ewald_spme) 403 WRITE (UNIT=iw, FMT="(T2,A,T35,A,T75,A,T80,A)") "-", "Ewald type", "SPME", "-" 404 WRITE (UNIT=iw, FMT="(T2,A,T35,A,T61,3I6,T80,A)") "-", "GMAX values", gmax, "-" 405 WRITE (UNIT=iw, FMT="(T2,A,T35,A,T73,I6,T80,A)") "-", "Spline order", o_spline, "-" 406 WRITE (UNIT=iw, FMT="(T2,A,T35,A,T67,F12.4,T80,A)") "-", "Alpha value", alpha, "-" 407 WRITE (UNIT=iw, FMT="(T2,A,T35,A,T67,F12.4,T80,A)") "-", "Real space cutoff value", rcut, "-" 408 END SELECT 409 WRITE (UNIT=iw, FMT="(T2,A)") REPEAT("-", 79) 410 END IF 411 CALL cp_print_key_finished_output(iw, logger, print_section, "PERIODIC_INFO") 412 413 END SUBROUTINE qmmm_ewald_potential_init 414 415! ************************************************************************************************** 416 417END MODULE qmmm_per_elpot 418 419