1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \par History 8!> none 9! ************************************************************************************************** 10MODULE dg_rho0_types 11 12 USE kinds, ONLY: dp 13 USE pw_grid_types, ONLY: pw_grid_type 14 USE pw_methods, ONLY: pw_zero 15 USE pw_poisson_types, ONLY: do_ewald_ewald,& 16 do_ewald_none,& 17 do_ewald_pme,& 18 do_ewald_spme 19 USE pw_types, ONLY: REALDATA3D,& 20 pw_create,& 21 pw_p_type,& 22 pw_release 23#include "../base/base_uses.f90" 24 25 IMPLICIT NONE 26 27 PRIVATE 28 PUBLIC:: dg_rho0_type, dg_rho0_init, dg_rho0_set, dg_rho0_get, & 29 dg_rho0_create, dg_rho0_retain, dg_rho0_release 30 31 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dg_rho0_types' 32 INTEGER, PRIVATE, SAVE :: last_dg_rho0_id_nr = 0 33 34! ************************************************************************************************** 35!> \brief Type for Gaussian Densities 36!> type = type of gaussian (PME) 37!> grid = grid number 38!> gcc = Gaussian contraction coefficient 39!> zet = Gaussian exponent 40! ************************************************************************************************** 41 TYPE dg_rho0_type 42 INTEGER :: ref_count, id_nr 43 INTEGER :: TYPE 44 INTEGER :: grid 45 INTEGER :: kind 46 REAL(KIND=dp) :: cutoff_radius 47 REAL(KIND=dp), DIMENSION(:), POINTER :: gcc 48 REAL(KIND=dp), DIMENSION(:), POINTER :: zet 49 TYPE(pw_p_type) :: density 50 END TYPE dg_rho0_type 51 52CONTAINS 53 54! ************************************************************************************************** 55!> \brief Set the dg_rho0_type 56!> \param dg_rho0 ... 57!> \param TYPE ... 58!> \param grid ... 59!> \param kind ... 60!> \param cutoff_radius ... 61!> \param gcc ... 62!> \param zet ... 63!> \param density ... 64!> \version 1.0 65! ************************************************************************************************** 66 SUBROUTINE dg_rho0_set(dg_rho0, TYPE, grid, kind, cutoff_radius, & 67 gcc, zet, density) 68 INTEGER, OPTIONAL :: TYPE 69 TYPE(dg_rho0_type), POINTER :: dg_rho0 70 INTEGER, OPTIONAL :: grid, kind 71 REAL(KIND=dp), OPTIONAL :: cutoff_radius 72 REAL(KIND=dp), OPTIONAL, POINTER :: gcc(:), zet(:) 73 TYPE(pw_p_type), OPTIONAL :: density 74 75 CHARACTER(LEN=*), PARAMETER :: routineN = 'dg_rho0_set', routineP = moduleN//':'//routineN 76 77 IF (PRESENT(grid)) dg_rho0%grid = grid 78 IF (PRESENT(kind)) dg_rho0%kind = kind 79 IF (PRESENT(density)) dg_rho0%density = density 80 IF (PRESENT(gcc)) dg_rho0%gcc => gcc 81 IF (PRESENT(zet)) dg_rho0%zet => zet 82 IF (PRESENT(TYPE)) dg_rho0%type = TYPE 83 IF (PRESENT(cutoff_radius)) dg_rho0%cutoff_radius = cutoff_radius 84 85 END SUBROUTINE dg_rho0_set 86 87! ************************************************************************************************** 88!> \brief Get the dg_rho0_type 89!> \param dg_rho0 ... 90!> \param cutoff_radius ... 91!> \param TYPE ... 92!> \param grid ... 93!> \param kind ... 94!> \param gcc ... 95!> \param zet ... 96!> \param density ... 97!> \version 1.0 98! ************************************************************************************************** 99 SUBROUTINE dg_rho0_get(dg_rho0, cutoff_radius, TYPE, grid, kind, gcc, zet, density) 100 INTEGER, OPTIONAL :: TYPE 101 REAL(KIND=dp), OPTIONAL :: cutoff_radius 102 TYPE(dg_rho0_type), POINTER :: dg_rho0 103 INTEGER, OPTIONAL :: grid, kind 104 REAL(KIND=dp), OPTIONAL, POINTER :: gcc(:), zet(:) 105 TYPE(pw_p_type), OPTIONAL, POINTER :: density 106 107 CHARACTER(LEN=*), PARAMETER :: routineN = 'dg_rho0_get', routineP = moduleN//':'//routineN 108 109 IF (PRESENT(grid)) grid = dg_rho0%grid 110 IF (PRESENT(kind)) kind = dg_rho0%kind 111 IF (PRESENT(density)) density = dg_rho0%density 112 IF (PRESENT(gcc)) gcc => dg_rho0%gcc 113 IF (PRESENT(zet)) zet => dg_rho0%zet 114 IF (PRESENT(TYPE)) TYPE = dg_rho0%type 115 IF (PRESENT(cutoff_radius)) cutoff_radius = dg_rho0%cutoff_radius 116 117 END SUBROUTINE dg_rho0_get 118 119! ************************************************************************************************** 120!> \brief create the dg_rho0 structure 121!> \param dg_rho0 ... 122!> \version 1.0 123! ************************************************************************************************** 124 SUBROUTINE dg_rho0_create(dg_rho0) 125 TYPE(dg_rho0_type), POINTER :: dg_rho0 126 127 CHARACTER(LEN=*), PARAMETER :: routineN = 'dg_rho0_create', routineP = moduleN//':'//routineN 128 129 ALLOCATE (dg_rho0) 130 NULLIFY (dg_rho0%gcc) 131 NULLIFY (dg_rho0%zet) 132 dg_rho0%cutoff_radius = 0.0_dp 133 dg_rho0%grid = 0 134 dg_rho0%kind = 0 135 dg_rho0%type = do_ewald_none 136 last_dg_rho0_id_nr = last_dg_rho0_id_nr + 1 137 dg_rho0%id_nr = last_dg_rho0_id_nr 138 dg_rho0%ref_count = 1 139 NULLIFY (dg_rho0%density%pw) 140 141 END SUBROUTINE dg_rho0_create 142 143! ************************************************************************************************** 144!> \brief retains the given dg_rho0_type 145!> \param dg_rho0 the dg_rho0_type to retain 146!> \par History 147!> 04.2003 created [fawzi] 148!> \author fawzi 149!> \note 150!> see doc/ReferenceCounting.html 151! ************************************************************************************************** 152 SUBROUTINE dg_rho0_retain(dg_rho0) 153 TYPE(dg_rho0_type), POINTER :: dg_rho0 154 155 CHARACTER(len=*), PARAMETER :: routineN = 'dg_rho0_retain', routineP = moduleN//':'//routineN 156 157 CPASSERT(ASSOCIATED(dg_rho0)) 158 CPASSERT(dg_rho0%ref_count > 0) 159 dg_rho0%ref_count = dg_rho0%ref_count + 1 160 END SUBROUTINE dg_rho0_retain 161 162! ************************************************************************************************** 163!> \brief releases the given dg_rho0_type 164!> \param dg_rho0 the dg_rho0_type to release 165!> \par History 166!> 04.2003 created [fawzi] 167!> \author fawzi 168!> \note 169!> see doc/ReferenceCounting.html 170! ************************************************************************************************** 171 SUBROUTINE dg_rho0_release(dg_rho0) 172 TYPE(dg_rho0_type), POINTER :: dg_rho0 173 174 CHARACTER(len=*), PARAMETER :: routineN = 'dg_rho0_release', & 175 routineP = moduleN//':'//routineN 176 177 IF (ASSOCIATED(dg_rho0)) THEN 178 CPASSERT(dg_rho0%ref_count > 0) 179 dg_rho0%ref_count = dg_rho0%ref_count - 1 180 IF (dg_rho0%ref_count == 0) THEN 181 IF (ASSOCIATED(dg_rho0%gcc)) THEN 182 DEALLOCATE (dg_rho0%gcc) 183 END IF 184 IF (ASSOCIATED(dg_rho0%zet)) THEN 185 DEALLOCATE (dg_rho0%zet) 186 END IF 187 CALL pw_release(dg_rho0%density%pw) 188 NULLIFY (dg_rho0%gcc) 189 NULLIFY (dg_rho0%zet) 190 DEALLOCATE (dg_rho0) 191 END IF 192 END IF 193 NULLIFY (dg_rho0) 194 END SUBROUTINE dg_rho0_release 195 196! ************************************************************************************************** 197!> \brief ... 198!> \param dg_rho0 ... 199!> \param pw_grid ... 200! ************************************************************************************************** 201 SUBROUTINE dg_rho0_init(dg_rho0, pw_grid) 202 TYPE(dg_rho0_type), POINTER :: dg_rho0 203 TYPE(pw_grid_type), POINTER :: pw_grid 204 205 CHARACTER(LEN=*), PARAMETER :: routineN = 'dg_rho0_init', routineP = moduleN//':'//routineN 206 207 CALL pw_release(dg_rho0%density%pw) 208 SELECT CASE (dg_rho0%type) 209 CASE (do_ewald_ewald) 210 CALL pw_create(dg_rho0%density%pw, pw_grid, REALDATA3D) 211 CALL dg_rho0_pme_gauss(dg_rho0%density, dg_rho0%zet(1)) 212 CASE (do_ewald_pme) 213 CALL pw_create(dg_rho0%density%pw, pw_grid, REALDATA3D) 214 CALL dg_rho0_pme_gauss(dg_rho0%density, dg_rho0%zet(1)) 215 CASE (do_ewald_spme) 216 CPABORT('SPME type not implemented') 217 END SELECT 218 219 END SUBROUTINE dg_rho0_init 220 221! ************************************************************************************************** 222!> \brief ... 223!> \param dg_rho0 ... 224!> \param alpha ... 225! ************************************************************************************************** 226 SUBROUTINE dg_rho0_pme_gauss(dg_rho0, alpha) 227 228 TYPE(pw_p_type), INTENT(INOUT) :: dg_rho0 229 REAL(KIND=dp), INTENT(IN) :: alpha 230 231 CHARACTER(LEN=*), PARAMETER :: routineN = 'dg_rho0_pme_gauss', & 232 routineP = moduleN//':'//routineN 233 INTEGER, PARAMETER :: IMPOSSIBLE = 10000 234 235 INTEGER :: gpt, l0, ln, lp, m0, mn, mp, n0, nn, np 236 INTEGER, DIMENSION(:), POINTER :: ghat 237 INTEGER, DIMENSION(:, :), POINTER :: bds 238 REAL(KIND=dp) :: const, e_gsq 239 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: rho0 240 TYPE(pw_grid_type), POINTER :: pw_grid 241 242 const = 1.0_dp/(8.0_dp*alpha**2) 243 244 pw_grid => dg_rho0%pw%pw_grid 245 bds => pw_grid%bounds 246 247 IF (-bds(1, 1) == bds(2, 1)) THEN 248 l0 = IMPOSSIBLE 249 ELSE 250 l0 = bds(1, 1) 251 END IF 252 253 IF (-bds(1, 2) == bds(2, 2)) THEN 254 m0 = IMPOSSIBLE 255 ELSE 256 m0 = bds(1, 2) 257 END IF 258 259 IF (-bds(1, 3) == bds(2, 3)) THEN 260 n0 = IMPOSSIBLE 261 ELSE 262 n0 = bds(1, 3) 263 END IF 264 265 CALL pw_zero(dg_rho0%pw) 266 267 rho0 => dg_rho0%pw%cr3d 268 269 DO gpt = 1, pw_grid%ngpts_cut_local 270 ghat => pw_grid%g_hat(:, gpt) 271 272 lp = pw_grid%mapl%pos(ghat(1)) 273 ln = pw_grid%mapl%neg(ghat(1)) 274 mp = pw_grid%mapm%pos(ghat(2)) 275 mn = pw_grid%mapm%neg(ghat(2)) 276 np = pw_grid%mapn%pos(ghat(3)) 277 nn = pw_grid%mapn%neg(ghat(3)) 278 279 e_gsq = EXP(-const*pw_grid%gsq(gpt))/pw_grid%vol 280 281 !*apsi 282 lp = lp + bds(1, 1) 283 mp = mp + bds(1, 2) 284 np = np + bds(1, 3) 285 ln = ln + bds(1, 1) 286 mn = mn + bds(1, 2) 287 nn = nn + bds(1, 3) 288 289 rho0(lp, mp, np) = e_gsq 290 rho0(ln, mn, nn) = e_gsq 291 292 IF (ghat(1) == l0 .OR. ghat(2) == m0 .OR. ghat(3) == n0) THEN 293 rho0(lp, mp, np) = 0.0_dp 294 rho0(ln, mn, nn) = 0.0_dp 295 END IF 296 297 END DO 298 299 END SUBROUTINE dg_rho0_pme_gauss 300 301END MODULE dg_rho0_types 302