1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Set of procedures to modify the density matrix, in order to collocate 8!> on the real space grid specific functions of the primitives 9!> \author MI 06.2006 10! ************************************************************************************************** 11MODULE qs_modify_pab_block 12 13 USE kinds, ONLY: dp 14 USE orbital_pointers, ONLY: coset 15#include "./base/base_uses.f90" 16 17 IMPLICIT NONE 18 19 PRIVATE 20 21 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_modify_pab_block' 22 23! *** Public subroutines *** 24 25 PUBLIC :: prepare_dadb, prepare_adb_m_dab, prepare_dab_p_adb, prepare_ardb_m_darb, prepare_arb 26 PUBLIC :: prepare_diadib, prepare_dijadijb, prepare_diiadiib 27 PUBLIC :: FUNC_AB, FUNC_ARB, FUNC_DADB, FUNC_ADBmDAB, FUNC_ARDBmDARB, FUNC_DABpADB 28 PUBLIC :: FUNC_DER, FUNC_DXDY, FUNC_DYDZ, FUNC_DZDX, FUNC_DX, FUNC_DY, FUNC_DZ 29 PUBLIC :: FUNC_DXDX, FUNC_DYDY, FUNC_DZDZ 30 31 ! Flags for the type of function to be collocated 32 INTEGER, PARAMETER :: FUNC_AB = 401, FUNC_DADB = 402 33 INTEGER, PARAMETER :: FUNC_ADBmDAB = 403, FUNC_ARDBmDARB = 404 34 INTEGER, PARAMETER :: FUNC_DABpADB = 405, FUNC_ARB = 406 35 INTEGER, DIMENSION(0:9), PARAMETER :: FUNC_DER = (/401, 501, 502, 503, 601, 602, 603, 604, 605, 606/) 36 INTEGER, PARAMETER :: FUNC_DXDY = 601, FUNC_DYDZ = 602, FUNC_DZDX = 603 37 INTEGER, PARAMETER :: FUNC_DXDX = 604, FUNC_DYDY = 605, FUNC_DZDZ = 606 38 INTEGER, PARAMETER :: FUNC_DX = 501, FUNC_DY = 502, FUNC_DZ = 503 39 40CONTAINS 41 42! 43 44! ************************************************************************************************** 45!> \brief create a new pab_local so that mapping pab_local with pgf_a pgf_b 46!> is equivalent to mapping pab with 0.5 * (nabla pgf_a) . (nabla pgf_b) 47!> (ddx pgf_a ) (ddx pgf_b) = (lax pgf_{a-1x} - 2*zeta*pgf_{a+1x})*(lbx pgf_{b-1x} - 2*zetb*pgf_{b+1x}) 48!> \param pab_local ... 49!> \param pab ... 50!> \param lxa ... 51!> \param lya ... 52!> \param lza ... 53!> \param lxb ... 54!> \param lyb ... 55!> \param lzb ... 56!> \param o1 ... 57!> \param o2 ... 58!> \param zeta ... 59!> \param zetb ... 60! ************************************************************************************************** 61 SUBROUTINE prepare_dadb(pab_local, pab, lxa, lya, lza, lxb, lyb, lzb, o1, o2, zeta, zetb) 62 63 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: pab_local 64 REAL(dp), DIMENSION(:, :), INTENT(IN) :: pab 65 INTEGER, INTENT(IN) :: lxa, lya, lza, lxb, lyb, lzb, o1, o2 66 REAL(dp), INTENT(IN) :: zeta, zetb 67 68 INTEGER :: ico, ico_l, jco, jco_l 69 70! this element of pab results in 12 elements of pab_local 71 72 ico = coset(lxa, lya, lza) 73 jco = coset(lxb, lyb, lzb) 74 ! x (all safe if lxa = 0, as the spurious added terms have zero prefactor) 75 76 ico_l = coset(MAX(lxa - 1, 0), lya, lza) 77 jco_l = coset(MAX(lxb - 1, 0), lyb, lzb) 78 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lxa*lxb*pab(o1 + ico, o2 + jco) 79 ico_l = coset(MAX(lxa - 1, 0), lya, lza) 80 jco_l = coset((lxb + 1), lyb, lzb) 81 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*lxa*zetb*pab(o1 + ico, o2 + jco) 82 ico_l = coset((lxa + 1), lya, lza) 83 jco_l = coset(MAX(lxb - 1, 0), lyb, lzb) 84 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zeta*lxb*pab(o1 + ico, o2 + jco) 85 ico_l = coset((lxa + 1), lya, lza) 86 jco_l = coset((lxb + 1), lyb, lzb) 87 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 4.0_dp*zeta*zetb*pab(o1 + ico, o2 + jco) 88 89 ! y 90 91 ico_l = coset(lxa, MAX(lya - 1, 0), lza) 92 jco_l = coset(lxb, MAX(lyb - 1, 0), lzb) 93 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lya*lyb*pab(o1 + ico, o2 + jco) 94 ico_l = coset(lxa, MAX(lya - 1, 0), lza) 95 jco_l = coset(lxb, (lyb + 1), lzb) 96 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*lya*zetb*pab(o1 + ico, o2 + jco) 97 ico_l = coset(lxa, (lya + 1), lza) 98 jco_l = coset(lxb, MAX(lyb - 1, 0), lzb) 99 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zeta*lyb*pab(o1 + ico, o2 + jco) 100 ico_l = coset(lxa, (lya + 1), lza) 101 jco_l = coset(lxb, (lyb + 1), lzb) 102 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 4.0_dp*zeta*zetb*pab(o1 + ico, o2 + jco) 103 104 ! z 105 106 ico_l = coset(lxa, lya, MAX(lza - 1, 0)) 107 jco_l = coset(lxb, lyb, MAX(lzb - 1, 0)) 108 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lza*lzb*pab(o1 + ico, o2 + jco) 109 ico_l = coset(lxa, lya, MAX(lza - 1, 0)) 110 jco_l = coset(lxb, lyb, (lzb + 1)) 111 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*lza*zetb*pab(o1 + ico, o2 + jco) 112 ico_l = coset(lxa, lya, (lza + 1)) 113 jco_l = coset(lxb, lyb, MAX(lzb - 1, 0)) 114 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zeta*lzb*pab(o1 + ico, o2 + jco) 115 ico_l = coset(lxa, lya, (lza + 1)) 116 jco_l = coset(lxb, lyb, (lzb + 1)) 117 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 4.0_dp*zeta*zetb*pab(o1 + ico, o2 + jco) 118 119 END SUBROUTINE prepare_dadb 120 121! 122! ************************************************************************************************** 123!> \brief create a new pab_local so that mapping pab_local with pgf_a pgf_b 124!> is equivalent to mapping pab with (ddi pgf_a) . (ddi pgf_b) 125!> (ddx pgf_a ) (ddx pgf_b) = (lax pgf_{a-1x} - 2*zeta*pgf_{a+1x})*(lbx pgf_{b-1x} - 2*zetb*pgf_{b+1x}) 126!> \param pab_local ... 127!> \param pab ... 128!> \param ider ... 129!> \param lxa ... 130!> \param lya ... 131!> \param lza ... 132!> \param lxb ... 133!> \param lyb ... 134!> \param lzb ... 135!> \param o1 ... 136!> \param o2 ... 137!> \param zeta ... 138!> \param zetb ... 139! ************************************************************************************************** 140 SUBROUTINE prepare_diadib(pab_local, pab, ider, lxa, lya, lza, lxb, lyb, lzb, o1, o2, zeta, zetb) 141 142 REAL(dp), DIMENSION(:, :), INTENT(inout) :: pab_local 143 REAL(dp), DIMENSION(:, :), INTENT(in) :: pab 144 INTEGER, INTENT(in) :: ider, lxa, lya, lza, lxb, lyb, lzb, o1, & 145 o2 146 REAL(dp), INTENT(in) :: zeta, zetb 147 148 INTEGER :: ico, ico_l, jco, jco_l 149 150! this element of pab results in 12 elements of pab_local 151 152 ico = coset(lxa, lya, lza) 153 jco = coset(lxb, lyb, lzb) 154 IF (ider == 1) THEN 155 ! x (all safe if lxa = 0, as the spurious added terms have zero prefactor) 156 ico_l = coset(MAX(lxa - 1, 0), lya, lza) 157 jco_l = coset(MAX(lxb - 1, 0), lyb, lzb) 158 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lxa*lxb*pab(o1 + ico, o2 + jco) 159 ico_l = coset(MAX(lxa - 1, 0), lya, lza) 160 jco_l = coset((lxb + 1), lyb, lzb) 161 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*lxa*zetb*pab(o1 + ico, o2 + jco) 162 ico_l = coset((lxa + 1), lya, lza) 163 jco_l = coset(MAX(lxb - 1, 0), lyb, lzb) 164 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zeta*lxb*pab(o1 + ico, o2 + jco) 165 ico_l = coset((lxa + 1), lya, lza) 166 jco_l = coset((lxb + 1), lyb, lzb) 167 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 4.0_dp*zeta*zetb*pab(o1 + ico, o2 + jco) 168 169 ELSEIF (ider == 2) THEN 170 ! y 171 172 ico_l = coset(lxa, MAX(lya - 1, 0), lza) 173 jco_l = coset(lxb, MAX(lyb - 1, 0), lzb) 174 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lya*lyb*pab(o1 + ico, o2 + jco) 175 ico_l = coset(lxa, MAX(lya - 1, 0), lza) 176 jco_l = coset(lxb, (lyb + 1), lzb) 177 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*lya*zetb*pab(o1 + ico, o2 + jco) 178 ico_l = coset(lxa, (lya + 1), lza) 179 jco_l = coset(lxb, MAX(lyb - 1, 0), lzb) 180 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zeta*lyb*pab(o1 + ico, o2 + jco) 181 ico_l = coset(lxa, (lya + 1), lza) 182 jco_l = coset(lxb, (lyb + 1), lzb) 183 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 4.0_dp*zeta*zetb*pab(o1 + ico, o2 + jco) 184 ELSEIF (ider == 3) THEN 185 ! z 186 ico_l = coset(lxa, lya, MAX(lza - 1, 0)) 187 jco_l = coset(lxb, lyb, MAX(lzb - 1, 0)) 188 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lza*lzb*pab(o1 + ico, o2 + jco) 189 ico_l = coset(lxa, lya, MAX(lza - 1, 0)) 190 jco_l = coset(lxb, lyb, (lzb + 1)) 191 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*lza*zetb*pab(o1 + ico, o2 + jco) 192 ico_l = coset(lxa, lya, (lza + 1)) 193 jco_l = coset(lxb, lyb, MAX(lzb - 1, 0)) 194 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zeta*lzb*pab(o1 + ico, o2 + jco) 195 ico_l = coset(lxa, lya, (lza + 1)) 196 jco_l = coset(lxb, lyb, (lzb + 1)) 197 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 4.0_dp*zeta*zetb*pab(o1 + ico, o2 + jco) 198 END IF 199 200 END SUBROUTINE prepare_diadib 201! 202! ************************************************************************************************** 203!> \brief create a new pab_local so that mapping pab_local with pgf_a pgf_b 204!> is equivalent to mapping pab with (ddidj pgf_a) . (ddidj pgf_b) 205!> (ddxdy pgf_a ) (ddxdy pgf_b) = 206!> \param pab_local ... 207!> \param pab ... 208!> \param ider1 ... 209!> \param ider2 ... 210!> \param lxa ... 211!> \param lya ... 212!> \param lza ... 213!> \param lxb ... 214!> \param lyb ... 215!> \param lzb ... 216!> \param o1 ... 217!> \param o2 ... 218!> \param zeta ... 219!> \param zetb ... 220! ************************************************************************************************** 221 SUBROUTINE prepare_dijadijb(pab_local, pab, ider1, ider2, lxa, lya, lza, lxb, lyb, lzb, o1, o2, zeta, zetb) 222 223 REAL(dp), DIMENSION(:, :), INTENT(inout) :: pab_local 224 REAL(dp), DIMENSION(:, :), INTENT(in) :: pab 225 INTEGER, INTENT(in) :: ider1, ider2, lxa, lya, lza, lxb, lyb, & 226 lzb, o1, o2 227 REAL(dp), INTENT(in) :: zeta, zetb 228 229 INTEGER :: ico, ico_l, jco 230 REAL(dp) :: func_a 231 232! this element of pab results in 12 elements of pab_local 233 234 ico = coset(lxa, lya, lza) 235 jco = coset(lxb, lyb, lzb) 236 237 IF ((ider1 == 1 .AND. ider2 == 2) .OR. (ider1 == 2 .AND. ider2 == 1)) THEN ! xy 238 ico_l = coset(MAX(lxa - 1, 0), MAX(lya - 1, 0), lza) 239 func_a = lxa*lya*pab(o1 + ico, o2 + jco) 240 CALL oneterm_dijdij(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=1) 241 ico_l = coset(lxa + 1, MAX(lya - 1, 0), lza) 242 func_a = -2.0_dp*zeta*lya*pab(o1 + ico, o2 + jco) 243 CALL oneterm_dijdij(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=1) 244 ico_l = coset(MAX(lxa - 1, 0), lya + 1, lza) 245 func_a = -2.0_dp*zeta*lxa*pab(o1 + ico, o2 + jco) 246 CALL oneterm_dijdij(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=1) 247 ico_l = coset(lxa + 1, lya + 1, lza) 248 func_a = 4.0_dp*zeta*zeta*pab(o1 + ico, o2 + jco) 249 CALL oneterm_dijdij(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=1) 250 ELSEIF ((ider1 == 2 .AND. ider2 == 3) .OR. (ider1 == 3 .AND. ider2 == 2)) THEN ! yz 251 ico_l = coset(lxa, MAX(lya - 1, 0), MAX(lza - 1, 0)) 252 func_a = lya*lza*pab(o1 + ico, o2 + jco) 253 CALL oneterm_dijdij(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=2) 254 ico_l = coset(lxa, lya + 1, MAX(lza - 1, 0)) 255 func_a = -2.0_dp*zeta*lza*pab(o1 + ico, o2 + jco) 256 CALL oneterm_dijdij(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=2) 257 ico_l = coset(lxa, MAX(lya - 1, 0), lza + 1) 258 func_a = -2.0_dp*zeta*lya*pab(o1 + ico, o2 + jco) 259 CALL oneterm_dijdij(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=2) 260 ico_l = coset(lxa, lya + 1, lza + 1) 261 func_a = 4.0_dp*zeta*zeta*pab(o1 + ico, o2 + jco) 262 CALL oneterm_dijdij(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=2) 263 ELSEIF ((ider1 == 3 .AND. ider2 == 1) .OR. (ider1 == 1 .AND. ider2 == 3)) THEN ! zx 264 ico_l = coset(MAX(lxa - 1, 0), lya, MAX(lza - 1, 0)) 265 func_a = lza*lxa*pab(o1 + ico, o2 + jco) 266 CALL oneterm_dijdij(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=3) 267 ico_l = coset(MAX(lxa - 1, 0), lya, lza + 1) 268 func_a = -2.0_dp*zeta*lxa*pab(o1 + ico, o2 + jco) 269 CALL oneterm_dijdij(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=3) 270 ico_l = coset(lxa + 1, lya, MAX(lza - 1, 0)) 271 func_a = -2.0_dp*zeta*lza*pab(o1 + ico, o2 + jco) 272 CALL oneterm_dijdij(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=3) 273 ico_l = coset(lxa + 1, lya, lza + 1) 274 func_a = 4.0_dp*zeta*zeta*pab(o1 + ico, o2 + jco) 275 CALL oneterm_dijdij(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=3) 276 END IF 277 278 END SUBROUTINE prepare_dijadijb 279 280! ************************************************************************************************** 281!> \brief ... 282!> \param pab_local ... 283!> \param func_a ... 284!> \param ico_l ... 285!> \param lx ... 286!> \param ly ... 287!> \param lz ... 288!> \param zet ... 289!> \param idir ... 290! ************************************************************************************************** 291 SUBROUTINE oneterm_dijdij(pab_local, func_a, ico_l, lx, ly, lz, zet, idir) 292 293 REAL(dp), DIMENSION(:, :), INTENT(inout) :: pab_local 294 REAL(dp), INTENT(in) :: func_a 295 INTEGER, INTENT(in) :: ico_l, lx, ly, lz 296 REAL(dp), INTENT(in) :: zet 297 INTEGER, INTENT(in) :: idir 298 299 INTEGER :: jco_l, l1, l2 300 301 IF (idir == 1) THEN 302 l1 = lx 303 l2 = ly 304 jco_l = coset(MAX(lx - 1, 0), MAX(ly - 1, 0), lz) 305 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + l1*l2*func_a 306 jco_l = coset(lx + 1, MAX(ly - 1, 0), lz) 307 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zet*l2*func_a 308 jco_l = coset(MAX(lx - 1, 0), ly + 1, lz) 309 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zet*l1*func_a 310 jco_l = coset(lx + 1, ly + 1, lz) 311 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 4.0_dp*zet*zet*func_a 312 ELSEIF (idir == 2) THEN 313 l1 = ly 314 l2 = lz 315 jco_l = coset(lx, MAX(ly - 1, 0), MAX(lz - 1, 0)) 316 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + l1*l2*func_a 317 jco_l = coset(lx, ly + 1, MAX(lz - 1, 0)) 318 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zet*l2*func_a 319 jco_l = coset(lx, MAX(ly - 1, 0), lz + 1) 320 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zet*l1*func_a 321 jco_l = coset(lx, ly + 1, lz + 1) 322 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 4.0_dp*zet*zet*func_a 323 ELSEIF (idir == 3) THEN 324 l1 = lz 325 l2 = lx 326 jco_l = coset(MAX(lx - 1, 0), ly, MAX(lz - 1, 0)) 327 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + l1*l2*func_a 328 jco_l = coset(MAX(lx - 1, 0), ly, lz + 1) 329 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zet*l2*func_a 330 jco_l = coset(lx + 1, ly, MAX(lz - 1, 0)) 331 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zet*l1*func_a 332 jco_l = coset(lx + 1, ly, lz + 1) 333 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 4.0_dp*zet*zet*func_a 334 335 END IF 336 END SUBROUTINE oneterm_dijdij 337 338! 339! ************************************************************************************************** 340!> \brief create a new pab_local so that mapping pab_local with pgf_a pgf_b 341!> is equivalent to mapping pab with (ddidj pgf_a) . (ddidj pgf_b) 342!> (ddxdx pgf_a ) (ddxdx pgf_b) = 343!> \param pab_local ... 344!> \param pab ... 345!> \param ider ... 346!> \param lxa ... 347!> \param lya ... 348!> \param lza ... 349!> \param lxb ... 350!> \param lyb ... 351!> \param lzb ... 352!> \param o1 ... 353!> \param o2 ... 354!> \param zeta ... 355!> \param zetb ... 356! ************************************************************************************************** 357 SUBROUTINE prepare_diiadiib(pab_local, pab, ider, lxa, lya, lza, lxb, lyb, lzb, o1, o2, zeta, zetb) 358 359 REAL(dp), DIMENSION(:, :), INTENT(inout) :: pab_local 360 REAL(dp), DIMENSION(:, :), INTENT(in) :: pab 361 INTEGER, INTENT(in) :: ider, lxa, lya, lza, lxb, lyb, lzb, o1, & 362 o2 363 REAL(dp), INTENT(in) :: zeta, zetb 364 365 INTEGER :: ico, ico_l, jco 366 REAL(dp) :: func_a 367 368! this element of pab results in 9 elements of pab_local 369 370 ico = coset(lxa, lya, lza) 371 jco = coset(lxb, lyb, lzb) 372 373 IF (ider == 1) THEN ! x 374 ico_l = coset(MAX(lxa - 2, 0), lya, lza) 375 func_a = lxa*(lxa - 1)*pab(o1 + ico, o2 + jco) 376 CALL oneterm_diidii(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=1) 377 ico_l = coset(lxa, lya, lza) 378 func_a = -2.0_dp*zeta*(2*lxa + 1)*pab(o1 + ico, o2 + jco) 379 CALL oneterm_diidii(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=1) 380 ico_l = coset(lxa + 2, lya, lza) 381 func_a = 4.0_dp*zeta*zeta*pab(o1 + ico, o2 + jco) 382 CALL oneterm_diidii(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=1) 383 ELSEIF (ider == 2) THEN ! y 384 ico_l = coset(lxa, MAX(lya - 2, 0), lza) 385 func_a = lya*(lya - 1)*pab(o1 + ico, o2 + jco) 386 CALL oneterm_diidii(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=2) 387 ico_l = coset(lxa, lya, lza) 388 func_a = -2.0_dp*zeta*(2*lya + 1)*pab(o1 + ico, o2 + jco) 389 CALL oneterm_diidii(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=2) 390 ico_l = coset(lxa, lya + 2, lza) 391 func_a = 4.0_dp*zeta*zeta*pab(o1 + ico, o2 + jco) 392 CALL oneterm_diidii(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=2) 393 ELSEIF (ider == 3) THEN ! z 394 ico_l = coset(lxa, lya, MAX(lza - 2, 0)) 395 func_a = lza*(lza - 1)*pab(o1 + ico, o2 + jco) 396 CALL oneterm_diidii(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=3) 397 ico_l = coset(lxa, lya, lza) 398 func_a = -2.0_dp*zeta*(2*lza + 1)*pab(o1 + ico, o2 + jco) 399 CALL oneterm_diidii(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=3) 400 ico_l = coset(lxa, lya, lza + 2) 401 func_a = 4.0_dp*zeta*zeta*pab(o1 + ico, o2 + jco) 402 CALL oneterm_diidii(pab_local, func_a, ico_l, lxb, lyb, lzb, zetb, idir=3) 403 END IF 404 405 END SUBROUTINE prepare_diiadiib 406 407! ************************************************************************************************** 408!> \brief ... 409!> \param pab_local ... 410!> \param func_a ... 411!> \param ico_l ... 412!> \param lx ... 413!> \param ly ... 414!> \param lz ... 415!> \param zet ... 416!> \param idir ... 417! ************************************************************************************************** 418 SUBROUTINE oneterm_diidii(pab_local, func_a, ico_l, lx, ly, lz, zet, idir) 419 REAL(dp), DIMENSION(:, :), INTENT(inout) :: pab_local 420 REAL(dp), INTENT(in) :: func_a 421 INTEGER, INTENT(in) :: ico_l, lx, ly, lz 422 REAL(dp), INTENT(in) :: zet 423 INTEGER, INTENT(in) :: idir 424 425 INTEGER :: jco_l, l1 426 427 IF (idir == 1) THEN 428 l1 = lx 429 jco_l = coset(MAX(lx - 2, 0), ly, lz) 430 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + l1*(l1 - 1)*func_a 431 jco_l = coset(lx, ly, lz) 432 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zet*(2*l1 + 1)*func_a 433 jco_l = coset(lx + 2, ly, lz) 434 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 4.0_dp*zet*zet*func_a 435 ELSEIF (idir == 2) THEN 436 l1 = ly 437 jco_l = coset(lx, MAX(ly - 2, 0), lz) 438 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + l1*(l1 - 1)*func_a 439 jco_l = coset(lx, ly, lz) 440 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zet*(2*l1 + 1)*func_a 441 jco_l = coset(lx, ly + 2, lz) 442 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 4.0_dp*zet*zet*func_a 443 ELSEIF (idir == 3) THEN 444 l1 = lz 445 jco_l = coset(lx, ly, MAX(lz - 2, 0)) 446 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + l1*(l1 - 1)*func_a 447 jco_l = coset(lx, ly, lz) 448 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zet*(2*l1 + 1)*func_a 449 jco_l = coset(lx, ly, lz + 2) 450 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 4.0_dp*zet*zet*func_a 451 END IF 452 453 END SUBROUTINE oneterm_diidii 454 455! ************************************************************************************************** 456!> \brief create a new pab_local so that mapping pab_local with pgf_a pgf_b 457!> is equivalent to mapping pab with pgf_a (nabla_{idir} pgf_b) - (nabla_{idir} pgf_a) pgf_b 458!> ( pgf_a ) (ddx pgf_b) - (ddx pgf_a)( pgf_b ) = 459!> pgf_a *(lbx pgf_{b-1x} - 2*zetb*pgf_{b+1x}) - (lax pgf_{a-1x} - 2*zeta*pgf_{a+1x}) pgf_b 460!> \param pab_local ... 461!> \param pab ... 462!> \param idir ... 463!> \param lxa ... 464!> \param lya ... 465!> \param lza ... 466!> \param lxb ... 467!> \param lyb ... 468!> \param lzb ... 469!> \param o1 ... 470!> \param o2 ... 471!> \param zeta ... 472!> \param zetb ... 473! ************************************************************************************************** 474 SUBROUTINE prepare_adb_m_dab(pab_local, pab, idir, lxa, lya, lza, lxb, lyb, lzb, o1, o2, zeta, zetb) 475 476 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: pab_local 477 REAL(dp), DIMENSION(:, :), INTENT(IN) :: pab 478 INTEGER, INTENT(IN) :: idir, lxa, lya, lza, lxb, lyb, lzb, o1, & 479 o2 480 REAL(dp), INTENT(IN) :: zeta, zetb 481 482 INTEGER :: ico, ico_l, jco, jco_l 483 484! this element of pab results in 4 elements of pab_local 485 486 ico = coset(lxa, lya, lza) 487 jco = coset(lxb, lyb, lzb) 488 489 IF (idir == 1) THEN ! x 490 ico_l = coset(lxa, lya, lza) 491 jco_l = coset(MAX(lxb - 1, 0), lyb, lzb) 492 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lxb*pab(o1 + ico, o2 + jco) 493 ico_l = coset(lxa, lya, lza) 494 jco_l = coset((lxb + 1), lyb, lzb) 495 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 496 ico_l = coset(MAX(lxa - 1, 0), lya, lza) 497 jco_l = coset(lxb, lyb, lzb) 498 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - lxa*pab(o1 + ico, o2 + jco) 499 ico_l = coset((lxa + 1), lya, lza) 500 jco_l = coset(lxb, lyb, lzb) 501 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 502 ELSEIF (idir == 2) THEN ! y 503 ico_l = coset(lxa, lya, lza) 504 jco_l = coset(lxb, MAX(lyb - 1, 0), lzb) 505 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lyb*pab(o1 + ico, o2 + jco) 506 ico_l = coset(lxa, lya, lza) 507 jco_l = coset(lxb, (lyb + 1), lzb) 508 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 509 ico_l = coset(lxa, MAX(lya - 1, 0), lza) 510 jco_l = coset(lxb, lyb, lzb) 511 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - lya*pab(o1 + ico, o2 + jco) 512 ico_l = coset(lxa, (lya + 1), lza) 513 jco_l = coset(lxb, lyb, lzb) 514 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 515 ELSE ! z 516 ico_l = coset(lxa, lya, lza) 517 jco_l = coset(lxb, lyb, MAX(lzb - 1, 0)) 518 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lzb*pab(o1 + ico, o2 + jco) 519 ico_l = coset(lxa, lya, lza) 520 jco_l = coset(lxb, lyb, (lzb + 1)) 521 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 522 ico_l = coset(lxa, lya, MAX(lza - 1, 0)) 523 jco_l = coset(lxb, lyb, lzb) 524 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - lza*pab(o1 + ico, o2 + jco) 525 ico_l = coset(lxa, lya, (lza + 1)) 526 jco_l = coset(lxb, lyb, lzb) 527 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 528 END IF 529 530 END SUBROUTINE prepare_adb_m_dab 531 532! 533! ************************************************************************************************** 534!> \brief create a new pab_local so that mapping pab_local with pgf_a pgf_b 535!> is equivalent to mapping pab with pgf_a (nabla_{idir} pgf_b) + (nabla_{idir} pgf_a) pgf_b 536!> ( pgf_a ) (ddx pgf_b) + (ddx pgf_a)( pgf_b ) = 537!> pgf_a *(lbx pgf_{b-1x} - 2*zetb*pgf_{b+1x}) + (lax pgf_{a-1x} - 2*zeta*pgf_{a+1x}) pgf_b 538!> \param pab_local ... 539!> \param pab ... 540!> \param idir ... 541!> \param lxa ... 542!> \param lya ... 543!> \param lza ... 544!> \param lxb ... 545!> \param lyb ... 546!> \param lzb ... 547!> \param o1 ... 548!> \param o2 ... 549!> \param zeta ... 550!> \param zetb ... 551! ************************************************************************************************** 552 SUBROUTINE prepare_dab_p_adb(pab_local, pab, idir, lxa, lya, lza, lxb, lyb, lzb, o1, o2, zeta, zetb) 553 554 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: pab_local 555 REAL(dp), DIMENSION(:, :), INTENT(IN) :: pab 556 INTEGER, INTENT(IN) :: idir, lxa, lya, lza, lxb, lyb, lzb, o1, & 557 o2 558 REAL(dp), INTENT(IN) :: zeta, zetb 559 560 INTEGER :: ico, ico_l, jco, jco_l 561 562! this element of pab results in 4 elements of pab_local 563 564 ico = coset(lxa, lya, lza) 565 jco = coset(lxb, lyb, lzb) 566 567 IF (idir == 1) THEN ! x 568 ico_l = coset(lxa, lya, lza) 569 jco_l = coset(MAX(lxb - 1, 0), lyb, lzb) 570 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lxb*pab(o1 + ico, o2 + jco) 571 ico_l = coset(lxa, lya, lza) 572 jco_l = coset((lxb + 1), lyb, lzb) 573 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 574 ico_l = coset(MAX(lxa - 1, 0), lya, lza) 575 jco_l = coset(lxb, lyb, lzb) 576 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lxa*pab(o1 + ico, o2 + jco) 577 ico_l = coset((lxa + 1), lya, lza) 578 jco_l = coset(lxb, lyb, lzb) 579 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 580 ELSEIF (idir == 2) THEN ! y 581 ico_l = coset(lxa, lya, lza) 582 jco_l = coset(lxb, MAX(lyb - 1, 0), lzb) 583 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lyb*pab(o1 + ico, o2 + jco) 584 ico_l = coset(lxa, lya, lza) 585 jco_l = coset(lxb, (lyb + 1), lzb) 586 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 587 ico_l = coset(lxa, MAX(lya - 1, 0), lza) 588 jco_l = coset(lxb, lyb, lzb) 589 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lya*pab(o1 + ico, o2 + jco) 590 ico_l = coset(lxa, (lya + 1), lza) 591 jco_l = coset(lxb, lyb, lzb) 592 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 593 ELSE ! z 594 ico_l = coset(lxa, lya, lza) 595 jco_l = coset(lxb, lyb, MAX(lzb - 1, 0)) 596 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lzb*pab(o1 + ico, o2 + jco) 597 ico_l = coset(lxa, lya, lza) 598 jco_l = coset(lxb, lyb, (lzb + 1)) 599 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 600 ico_l = coset(lxa, lya, MAX(lza - 1, 0)) 601 jco_l = coset(lxb, lyb, lzb) 602 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lza*pab(o1 + ico, o2 + jco) 603 ico_l = coset(lxa, lya, (lza + 1)) 604 jco_l = coset(lxb, lyb, lzb) 605 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 606 END IF 607 608 END SUBROUTINE prepare_dab_p_adb 609 610! 611 612! ************************************************************************************************** 613!> \brief create a new pab_local so that mapping pab_local with pgf_a pgf_b 614!> pgf_a (r-Rb)_{ir} (nabla_{idir} pgf_b) - (nabla_{idir} pgf_a) (r-Rb)_{ir} pgf_b 615!> ( pgf_a )(r-Rb)_{ir} (ddx pgf_b) - (ddx pgf_a) (r-Rb)_{ir} ( pgf_b ) = 616!> pgf_a *(lbx pgf_{b-1x+1ir} - 2*zetb*pgf_{b+1x+1ir}) - 617!> (lax pgf_{a-1x} - 2*zeta*pgf_{a+1x}) pgf_{b+1ir} 618!> \param pab_local ... 619!> \param pab ... 620!> \param idir ... 621!> \param ir ... 622!> \param lxa ... 623!> \param lya ... 624!> \param lza ... 625!> \param lxb ... 626!> \param lyb ... 627!> \param lzb ... 628!> \param o1 ... 629!> \param o2 ... 630!> \param zeta ... 631!> \param zetb ... 632! ************************************************************************************************** 633 SUBROUTINE prepare_ardb_m_darb(pab_local, pab, idir, ir, lxa, lya, lza, lxb, lyb, lzb, o1, o2, zeta, zetb) 634 635 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: pab_local 636 REAL(dp), DIMENSION(:, :), INTENT(IN) :: pab 637 INTEGER, INTENT(IN) :: idir, ir, lxa, lya, lza, lxb, lyb, lzb, & 638 o1, o2 639 REAL(dp), INTENT(IN) :: zeta, zetb 640 641 INTEGER :: ico, ico_l, jco, jco_l 642 643! this element of pab results in 4 elements of pab_local 644 645 ico = coset(lxa, lya, lza) 646 jco = coset(lxb, lyb, lzb) 647 648 IF (idir == 1 .AND. ir == 1) THEN 649 650 ico_l = coset(lxa, lya, lza) 651 jco_l = coset(lxb, lyb, lzb) 652 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lxb*pab(o1 + ico, o2 + jco) 653 654 ico_l = coset(lxa, lya, lza) 655 jco_l = coset((lxb + 2), lyb, lzb) 656 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 657 658 ico_l = coset(MAX(lxa - 1, 0), lya, lza) 659 jco_l = coset((lxb + 1), lyb, lzb) 660 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - lxa*pab(o1 + ico, o2 + jco) 661 662 ico_l = coset((lxa + 1), lya, lza) 663 jco_l = coset((lxb + 1), lyb, lzb) 664 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 665 666 ELSEIF (idir == 1 .AND. ir == 2) THEN 667 668 ico_l = coset(lxa, lya, lza) 669 jco_l = coset(MAX(lxb - 1, 0), (lyb + 1), lzb) 670 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lxb*pab(o1 + ico, o2 + jco) 671 672 ico_l = coset(lxa, lya, lza) 673 jco_l = coset((lxb + 1), (lyb + 1), lzb) 674 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 675 676 ico_l = coset(MAX(lxa - 1, 0), lya, lza) 677 jco_l = coset(lxb, (lyb + 1), lzb) 678 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - lxa*pab(o1 + ico, o2 + jco) 679 680 ico_l = coset((lxa + 1), lya, lza) 681 jco_l = coset(lxb, (lyb + 1), lzb) 682 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 683 684 ELSEIF (idir == 1 .AND. ir == 3) THEN 685 686 ico_l = coset(lxa, lya, lza) 687 jco_l = coset(MAX(lxb - 1, 0), lyb, (lzb + 1)) 688 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lxb*pab(o1 + ico, o2 + jco) 689 690 ico_l = coset(lxa, lya, lza) 691 jco_l = coset((lxb + 1), lyb, (lzb + 1)) 692 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 693 694 ico_l = coset(MAX(lxa - 1, 0), lya, lza) 695 jco_l = coset(lxb, lyb, (lzb + 1)) 696 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - lxa*pab(o1 + ico, o2 + jco) 697 698 ico_l = coset((lxa + 1), lya, lza) 699 jco_l = coset(lxb, lyb, (lzb + 1)) 700 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 701 702 ELSEIF (idir == 2 .AND. ir == 1) THEN 703 704 ico_l = coset(lxa, lya, lza) 705 jco_l = coset((lxb + 1), MAX(lyb - 1, 0), lzb) 706 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lyb*pab(o1 + ico, o2 + jco) 707 708 ico_l = coset(lxa, lya, lza) 709 jco_l = coset((lxb + 1), (lyb + 1), lzb) 710 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 711 712 ico_l = coset(lxa, MAX(lya - 1, 0), lza) 713 jco_l = coset((lxb + 1), lyb, lzb) 714 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - lya*pab(o1 + ico, o2 + jco) 715 716 ico_l = coset(lxa, (lya + 1), lza) 717 jco_l = coset((lxb + 1), lyb, lzb) 718 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 719 720 ELSEIF (idir == 2 .AND. ir == 2) THEN 721 722 ico_l = coset(lxa, lya, lza) 723 jco_l = coset(lxb, lyb, lzb) 724 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lyb*pab(o1 + ico, o2 + jco) 725 726 ico_l = coset(lxa, lya, lza) 727 jco_l = coset(lxb, (lyb + 2), lzb) 728 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 729 730 ico_l = coset(lxa, MAX(lya - 1, 0), lza) 731 jco_l = coset(lxb, (lyb + 1), lzb) 732 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - lya*pab(o1 + ico, o2 + jco) 733 734 ico_l = coset(lxa, (lya + 1), lza) 735 jco_l = coset(lxb, (lyb + 1), lzb) 736 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 737 738 ELSEIF (idir == 2 .AND. ir == 3) THEN 739 740 ico_l = coset(lxa, lya, lza) 741 jco_l = coset(lxb, MAX(lyb - 1, 0), (lzb + 1)) 742 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lyb*pab(o1 + ico, o2 + jco) 743 744 ico_l = coset(lxa, lya, lza) 745 jco_l = coset(lxb, (lyb + 1), (lzb + 1)) 746 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 747 748 ico_l = coset(lxa, MAX(lya - 1, 0), lza) 749 jco_l = coset(lxb, lyb, (lzb + 1)) 750 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - lya*pab(o1 + ico, o2 + jco) 751 752 ico_l = coset(lxa, (lya + 1), lza) 753 jco_l = coset(lxb, lyb, (lzb + 1)) 754 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 755 756 ELSEIF (idir == 3 .AND. ir == 1) THEN 757 758 ico_l = coset(lxa, lya, lza) 759 jco_l = coset((lxb + 1), lyb, MAX(lzb - 1, 0)) 760 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lzb*pab(o1 + ico, o2 + jco) 761 762 ico_l = coset(lxa, lya, lza) 763 jco_l = coset((lxb + 1), lyb, (lzb + 1)) 764 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 765 766 ico_l = coset(lxa, lya, MAX(lza - 1, 0)) 767 jco_l = coset((lxb + 1), lyb, lzb) 768 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - lza*pab(o1 + ico, o2 + jco) 769 770 ico_l = coset(lxa, lya, (lza + 1)) 771 jco_l = coset((lxb + 1), lyb, lzb) 772 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 773 774 ELSEIF (idir == 3 .AND. ir == 2) THEN 775 776 ico_l = coset(lxa, lya, lza) 777 jco_l = coset(lxb, (lyb + 1), MAX(lzb - 1, 0)) 778 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lzb*pab(o1 + ico, o2 + jco) 779 780 ico_l = coset(lxa, lya, lza) 781 jco_l = coset(lxb, (lyb + 1), (lzb + 1)) 782 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 783 784 ico_l = coset(lxa, lya, MAX(lza - 1, 0)) 785 jco_l = coset(lxb, (lyb + 1), lzb) 786 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - lza*pab(o1 + ico, o2 + jco) 787 788 ico_l = coset(lxa, lya, (lza + 1)) 789 jco_l = coset(lxb, (lyb + 1), lzb) 790 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 791 792 ELSEIF (idir == 3 .AND. ir == 3) THEN 793 794 ico_l = coset(lxa, lya, lza) 795 jco_l = coset(lxb, lyb, lzb) 796 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + lzb*pab(o1 + ico, o2 + jco) 797 798 ico_l = coset(lxa, lya, lza) 799 jco_l = coset(lxb, lyb, (lzb + 2)) 800 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - 2.0_dp*zetb*pab(o1 + ico, o2 + jco) 801 802 ico_l = coset(lxa, lya, MAX(lza - 1, 0)) 803 jco_l = coset(lxb, lyb, (lzb + 1)) 804 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) - lza*pab(o1 + ico, o2 + jco) 805 806 ico_l = coset(lxa, lya, (lza + 1)) 807 jco_l = coset(lxb, lyb, (lzb + 1)) 808 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + 2.0_dp*zeta*pab(o1 + ico, o2 + jco) 809 810 END IF 811 812 END SUBROUTINE prepare_ardb_m_darb 813 814! ************************************************************************************************** 815!> \brief create a new pab_local so that mapping pab_local with pgf_a pgf_b 816!> pgf_a (r-Rb)_{ir} pgf_b = pgf_a * pgf_b{+1ir} 817!> \param pab_local ... 818!> \param pab ... 819!> \param ir ... 820!> \param lxa ... 821!> \param lya ... 822!> \param lza ... 823!> \param lxb ... 824!> \param lyb ... 825!> \param lzb ... 826!> \param o1 ... 827!> \param o2 ... 828! ************************************************************************************************** 829 SUBROUTINE prepare_arb(pab_local, pab, ir, lxa, lya, lza, lxb, lyb, lzb, o1, o2) 830 831 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: pab_local 832 REAL(dp), DIMENSION(:, :), INTENT(IN) :: pab 833 INTEGER, INTENT(IN) :: ir, lxa, lya, lza, lxb, lyb, lzb, o1, o2 834 835 INTEGER :: ico, ico_l, jco, jco_l 836 837 ico = coset(lxa, lya, lza) 838 jco = coset(lxb, lyb, lzb) 839 840 IF (ir == 1) THEN 841 842 ico_l = coset(lxa, lya, lza) 843 jco_l = coset((lxb + 1), lyb, lzb) 844 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + pab(o1 + ico, o2 + jco) 845 846 ELSEIF (ir == 2) THEN 847 848 ico_l = coset(lxa, lya, lza) 849 jco_l = coset(lxb, (lyb + 1), lzb) 850 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + pab(o1 + ico, o2 + jco) 851 852 ELSEIF (ir == 3) THEN 853 854 ico_l = coset(lxa, lya, lza) 855 jco_l = coset(lxb, lyb, (lzb + 1)) 856 pab_local(ico_l, jco_l) = pab_local(ico_l, jco_l) + pab(o1 + ico, o2 + jco) 857 858 END IF 859 860 END SUBROUTINE prepare_arb 861 862END MODULE qs_modify_pab_block 863 864