1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5! ************************************************************************************************** 6!> \brief Calculation of the nuclear attraction contribution to the core Hamiltonian 7!> <a|erfc|b> :we only calculate the non-screened part 8!> \par History 9!> - core_ppnl refactored from qs_core_hamiltonian [Joost VandeVondele, 2008-11-01] 10!> - adapted for nuclear attraction [jhu, 2009-02-24] 11! ************************************************************************************************** 12MODULE core_ae 13 USE ai_verfc, ONLY: verfc 14 USE atomic_kind_types, ONLY: atomic_kind_type,& 15 get_atomic_kind_set 16 USE basis_set_types, ONLY: gto_basis_set_p_type,& 17 gto_basis_set_type 18 USE dbcsr_api, ONLY: dbcsr_add,& 19 dbcsr_get_block_p,& 20 dbcsr_p_type 21 USE external_potential_types, ONLY: all_potential_type,& 22 get_potential,& 23 sgp_potential_type 24 USE kinds, ONLY: dp 25 USE orbital_pointers, ONLY: coset,& 26 indco,& 27 init_orbital_pointers,& 28 ncoset 29 USE particle_types, ONLY: particle_type 30 USE qs_force_types, ONLY: qs_force_type 31 USE qs_kind_types, ONLY: get_qs_kind,& 32 get_qs_kind_set,& 33 qs_kind_type 34 USE qs_neighbor_list_types, ONLY: get_iterator_info,& 35 neighbor_list_iterate,& 36 neighbor_list_iterator_create,& 37 neighbor_list_iterator_p_type,& 38 neighbor_list_iterator_release,& 39 neighbor_list_set_p_type,& 40 nl_set_sub_iterator,& 41 nl_sub_iterate 42 USE virial_methods, ONLY: virial_pair_force 43 USE virial_types, ONLY: virial_type 44#include "./base/base_uses.f90" 45 46 IMPLICIT NONE 47 48 PRIVATE 49 50 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'core_ae' 51 52 PUBLIC :: build_core_ae 53 54CONTAINS 55 56! ************************************************************************************************** 57!> \brief ... 58!> \param matrix_h ... 59!> \param matrix_p ... 60!> \param force ... 61!> \param virial ... 62!> \param calculate_forces ... 63!> \param use_virial ... 64!> \param nder ... 65!> \param qs_kind_set ... 66!> \param atomic_kind_set ... 67!> \param particle_set ... 68!> \param sab_orb ... 69!> \param sac_ae ... 70!> \param nimages ... 71!> \param cell_to_index ... 72! ************************************************************************************************** 73 SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, & 74 qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, nimages, cell_to_index) 75 76 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p 77 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 78 TYPE(virial_type), POINTER :: virial 79 LOGICAL, INTENT(IN) :: calculate_forces 80 LOGICAL :: use_virial 81 INTEGER :: nder 82 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 83 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 84 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 85 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 86 POINTER :: sab_orb, sac_ae 87 INTEGER, INTENT(IN) :: nimages 88 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index 89 90 CHARACTER(LEN=*), PARAMETER :: routineN = 'build_core_ae', routineP = moduleN//':'//routineN 91 92 INTEGER :: atom_a, atom_b, atom_c, handle, iatom, icol, ikind, img, inode, irow, iset, & 93 jatom, jkind, jset, katom, kkind, ldai, ldsab, maxco, maxder, maxl, maxnset, maxsgf, & 94 na_plus, natom, nb_plus, ncoa, ncob, nij, nkind, nseta, nsetb, sgfa, sgfb 95 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind 96 INTEGER, DIMENSION(3) :: cellind 97 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, & 98 npgfb, nsgfa, nsgfb 99 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 100 LOGICAL :: dokp, found 101 REAL(KIND=dp) :: alpha_c, core_charge, core_radius, dab, & 102 dac, dbc, f0, rab2, rac2, rbc2, zeta_c 103 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ff 104 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: habd, work 105 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: hab, pab, verf, vnuc 106 REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, rab, rac, rbc 107 REAL(KIND=dp), DIMENSION(3, 3) :: pv_loc 108 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b 109 REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block, p_block, rpgfa, rpgfb, sphi_a, & 110 sphi_b, zeta, zetb 111 TYPE(all_potential_type), POINTER :: all_potential 112 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list 113 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b 114 TYPE(neighbor_list_iterator_p_type), & 115 DIMENSION(:), POINTER :: ap_iterator, nl_iterator 116 TYPE(sgp_potential_type), POINTER :: sgp_potential 117 118 IF (calculate_forces) THEN 119 CALL timeset(routineN//"_forces", handle) 120 ELSE 121 CALL timeset(routineN, handle) 122 ENDIF 123 124 nkind = SIZE(atomic_kind_set) 125 natom = SIZE(particle_set) 126 127 dokp = (nimages > 1) 128 129 ALLOCATE (atom_of_kind(natom)) 130 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind) 131 132 IF (calculate_forces) THEN 133 IF (SIZE(matrix_p, 1) == 2) THEN 134 DO img = 1, nimages 135 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, & 136 alpha_scalar=1.0_dp, beta_scalar=1.0_dp) 137 CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, & 138 alpha_scalar=-2.0_dp, beta_scalar=1.0_dp) 139 END DO 140 END IF 141 END IF 142 143 IF (use_virial) THEN 144 pv_loc = virial%pv_virial 145 END IF 146 147 maxder = ncoset(nder) 148 149 CALL get_qs_kind_set(qs_kind_set, & 150 maxco=maxco, maxlgto=maxl, maxsgf=maxsgf, maxnset=maxnset) 151 152 CALL init_orbital_pointers(maxl + nder + 1) 153 154 ldsab = MAX(maxco, maxsgf) 155 ldai = ncoset(maxl + nder + 1) 156 ALLOCATE (hab(ldsab, ldsab, maxnset*maxnset), work(ldsab, ldsab)) 157 ALLOCATE (verf(ldai, ldai, 2*maxl + nder + 1), vnuc(ldai, ldai, 2*maxl + nder + 1), ff(0:2*maxl + nder)) 158 IF (calculate_forces) THEN 159 ALLOCATE (pab(maxco, maxco, maxnset*maxnset)) 160 END IF 161 162 ! iterator for basis/potential list 163 CALL neighbor_list_iterator_create(ap_iterator, sac_ae, search=.TRUE.) 164 165 ALLOCATE (basis_set_list(nkind)) 166 DO ikind = 1, nkind 167 CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a) 168 IF (ASSOCIATED(basis_set_a)) THEN 169 basis_set_list(ikind)%gto_basis_set => basis_set_a 170 ELSE 171 NULLIFY (basis_set_list(ikind)%gto_basis_set) 172 END IF 173 END DO 174 CALL neighbor_list_iterator_create(nl_iterator, sab_orb) 175 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 176 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, & 177 iatom=iatom, jatom=jatom, r=rab, cell=cellind) 178 basis_set_a => basis_set_list(ikind)%gto_basis_set 179 IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE 180 basis_set_b => basis_set_list(jkind)%gto_basis_set 181 IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE 182 atom_a = atom_of_kind(iatom) 183 atom_b = atom_of_kind(jatom) 184 ! basis ikind 185 first_sgfa => basis_set_a%first_sgf 186 la_max => basis_set_a%lmax 187 la_min => basis_set_a%lmin 188 npgfa => basis_set_a%npgf 189 nseta = basis_set_a%nset 190 nsgfa => basis_set_a%nsgf_set 191 rpgfa => basis_set_a%pgf_radius 192 set_radius_a => basis_set_a%set_radius 193 sphi_a => basis_set_a%sphi 194 zeta => basis_set_a%zet 195 ! basis jkind 196 first_sgfb => basis_set_b%first_sgf 197 lb_max => basis_set_b%lmax 198 lb_min => basis_set_b%lmin 199 npgfb => basis_set_b%npgf 200 nsetb = basis_set_b%nset 201 nsgfb => basis_set_b%nsgf_set 202 rpgfb => basis_set_b%pgf_radius 203 set_radius_b => basis_set_b%set_radius 204 sphi_b => basis_set_b%sphi 205 zetb => basis_set_b%zet 206 207 dab = SQRT(SUM(rab*rab)) 208 209 IF (dokp) THEN 210 img = cell_to_index(cellind(1), cellind(2), cellind(3)) 211 ELSE 212 img = 1 213 END IF 214 215 ! *** Use the symmetry of the first derivatives *** 216 IF (iatom == jatom) THEN 217 f0 = 1.0_dp 218 ELSE 219 f0 = 2.0_dp 220 END IF 221 222 ! *** Create matrix blocks for a new matrix block column *** 223 IF (iatom <= jatom) THEN 224 irow = iatom 225 icol = jatom 226 ELSE 227 irow = jatom 228 icol = iatom 229 END IF 230 NULLIFY (h_block) 231 CALL dbcsr_get_block_p(matrix=matrix_h(1, img)%matrix, & 232 row=irow, col=icol, BLOCK=h_block, found=found) 233 IF (calculate_forces) THEN 234 NULLIFY (p_block) 235 CALL dbcsr_get_block_p(matrix=matrix_p(1, img)%matrix, & 236 row=irow, col=icol, BLOCK=p_block, found=found) 237 CPASSERT(ASSOCIATED(p_block)) 238 ! *** Decontract density matrix block *** 239 DO iset = 1, nseta 240 ncoa = npgfa(iset)*ncoset(la_max(iset)) 241 sgfa = first_sgfa(1, iset) 242 DO jset = 1, nsetb 243 ncob = npgfb(jset)*ncoset(lb_max(jset)) 244 sgfb = first_sgfb(1, jset) 245 nij = jset + (iset - 1)*maxnset 246 IF (iatom <= jatom) THEN 247 CALL dgemm("N", "N", ncoa, nsgfb(jset), nsgfa(iset), & 248 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), & 249 p_block(sgfa, sgfb), SIZE(p_block, 1), & 250 0.0_dp, work(1, 1), SIZE(work, 1)) 251 ELSE 252 CALL dgemm("N", "T", ncoa, nsgfb(jset), nsgfa(iset), & 253 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), & 254 p_block(sgfb, sgfa), SIZE(p_block, 1), & 255 0.0_dp, work(1, 1), SIZE(work, 1)) 256 END IF 257 CALL dgemm("N", "T", ncoa, ncob, nsgfb(jset), & 258 1.0_dp, work(1, 1), SIZE(work, 1), & 259 sphi_b(1, sgfb), SIZE(sphi_b, 1), & 260 0.0_dp, pab(1, 1, nij), SIZE(pab, 1)) 261 END DO 262 END DO 263 END IF 264 265 ! loop over all kinds for pseudopotential atoms 266 hab = 0._dp 267 DO kkind = 1, nkind 268 CALL get_qs_kind(qs_kind_set(kkind), all_potential=all_potential, & 269 sgp_potential=sgp_potential) 270 IF (ASSOCIATED(all_potential)) THEN 271 CALL get_potential(potential=all_potential, & 272 alpha_core_charge=alpha_c, zeff=zeta_c, & 273 ccore_charge=core_charge, core_charge_radius=core_radius) 274 ELSE IF (ASSOCIATED(sgp_potential)) THEN 275 CALL get_potential(potential=sgp_potential, & 276 alpha_core_charge=alpha_c, zeff=zeta_c, & 277 ccore_charge=core_charge, core_charge_radius=core_radius) 278 ELSE 279 CYCLE 280 END IF 281 282 CALL nl_set_sub_iterator(ap_iterator, ikind, kkind, iatom) 283 DO WHILE (nl_sub_iterate(ap_iterator) == 0) 284 CALL get_iterator_info(ap_iterator, jatom=katom, r=rac) 285 dac = SQRT(SUM(rac*rac)) 286 rbc(:) = rac(:) - rab(:) 287 dbc = SQRT(SUM(rbc*rbc)) 288 IF ((MAXVAL(set_radius_a(:)) + core_radius < dac) .OR. & 289 (MAXVAL(set_radius_b(:)) + core_radius < dbc)) THEN 290 CYCLE 291 END IF 292 293 DO iset = 1, nseta 294 IF (set_radius_a(iset) + core_radius < dac) CYCLE 295 ncoa = npgfa(iset)*ncoset(la_max(iset)) 296 sgfa = first_sgfa(1, iset) 297 DO jset = 1, nsetb 298 IF (set_radius_b(jset) + core_radius < dbc) CYCLE 299 ncob = npgfb(jset)*ncoset(lb_max(jset)) 300 sgfb = first_sgfb(1, jset) 301 IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE 302 rab2 = dab*dab 303 rac2 = dac*dac 304 rbc2 = dbc*dbc 305 nij = jset + (iset - 1)*maxnset 306 ! *** Calculate the GTH pseudo potential forces *** 307 IF (calculate_forces) THEN 308 na_plus = npgfa(iset)*ncoset(la_max(iset) + nder) 309 nb_plus = npgfb(jset)*ncoset(lb_max(jset)) 310 ALLOCATE (habd(na_plus, nb_plus)) 311 habd = 0._dp 312 CALL verfc( & 313 la_max(iset) + nder, npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), & 314 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), & 315 alpha_c, core_radius, zeta_c, core_charge, & 316 rab, rab2, rac, rac2, rbc2, hab(:, :, nij), verf, vnuc, ff(0:), & 317 nder, habd) 318 319 ! *** The derivatives w.r.t. atomic center c are *** 320 ! *** calculated using the translational invariance *** 321 ! *** of the first derivatives *** 322 CALL verfc_force(habd, pab(:, :, nij), force_a, force_b, nder, & 323 la_max(iset), la_min(iset), npgfa(iset), zeta(:, iset), & 324 lb_max(jset), lb_min(jset), npgfb(jset), zetb(:, jset), rab) 325 326 DEALLOCATE (habd) 327 328 atom_c = atom_of_kind(katom) 329 force(ikind)%all_potential(1, atom_a) = force(ikind)%all_potential(1, atom_a) + f0*force_a(1) 330 force(ikind)%all_potential(2, atom_a) = force(ikind)%all_potential(2, atom_a) + f0*force_a(2) 331 force(ikind)%all_potential(3, atom_a) = force(ikind)%all_potential(3, atom_a) + f0*force_a(3) 332 333 force(jkind)%all_potential(1, atom_b) = force(jkind)%all_potential(1, atom_b) + f0*force_b(1) 334 force(jkind)%all_potential(2, atom_b) = force(jkind)%all_potential(2, atom_b) + f0*force_b(2) 335 force(jkind)%all_potential(3, atom_b) = force(jkind)%all_potential(3, atom_b) + f0*force_b(3) 336 337 force(kkind)%all_potential(1, atom_c) = force(kkind)%all_potential(1, atom_c) & 338 - f0*force_a(1) - f0*force_b(1) 339 force(kkind)%all_potential(2, atom_c) = force(kkind)%all_potential(2, atom_c) & 340 - f0*force_a(2) - f0*force_b(2) 341 force(kkind)%all_potential(3, atom_c) = force(kkind)%all_potential(3, atom_c) & 342 - f0*force_a(3) - f0*force_b(3) 343 344 IF (use_virial) THEN 345 CALL virial_pair_force(virial%pv_virial, f0, force_a, rac) 346 CALL virial_pair_force(virial%pv_virial, f0, force_b, rbc) 347 END IF 348 ELSE 349 CALL verfc( & 350 la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), & 351 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), & 352 alpha_c, core_radius, zeta_c, core_charge, & 353 rab, rab2, rac, rac2, rbc2, hab(:, :, nij), verf, vnuc, ff(0:)) 354 END IF 355 END DO 356 END DO 357 END DO 358 END DO 359 ! *** Contract nuclear attraction integrals 360 DO iset = 1, nseta 361 ncoa = npgfa(iset)*ncoset(la_max(iset)) 362 sgfa = first_sgfa(1, iset) 363 DO jset = 1, nsetb 364 ncob = npgfb(jset)*ncoset(lb_max(jset)) 365 sgfb = first_sgfb(1, jset) 366 nij = jset + (iset - 1)*maxnset 367 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, & 368 1.0_dp, hab(1, 1, nij), SIZE(hab, 1), & 369 sphi_b(1, sgfb), SIZE(sphi_b, 1), & 370 0.0_dp, work(1, 1), SIZE(work, 1)) 371 IF (iatom <= jatom) THEN 372 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, & 373 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), & 374 work(1, 1), SIZE(work, 1), & 375 1.0_dp, h_block(sgfa, sgfb), SIZE(h_block, 1)) 376 ELSE 377 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, & 378 1.0_dp, work(1, 1), SIZE(work, 1), & 379 sphi_a(1, sgfa), SIZE(sphi_a, 1), & 380 1.0_dp, h_block(sgfb, sgfa), SIZE(h_block, 1)) 381 END IF 382 END DO 383 END DO 384 385 END DO 386 CALL neighbor_list_iterator_release(nl_iterator) 387 388 CALL neighbor_list_iterator_release(ap_iterator) 389 390 DEALLOCATE (atom_of_kind, basis_set_list) 391 DEALLOCATE (hab, work, verf, vnuc, ff) 392 393 IF (calculate_forces) THEN 394 DEALLOCATE (pab) 395 END IF 396 IF (calculate_forces) THEN 397 ! *** If LSD, then recover alpha density and beta density *** 398 ! *** from the total density (1) and the spin density (2) *** 399 IF (SIZE(matrix_p, 1) == 2) THEN 400 DO img = 1, nimages 401 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, & 402 alpha_scalar=0.5_dp, beta_scalar=0.5_dp) 403 CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, & 404 alpha_scalar=-1.0_dp, beta_scalar=1.0_dp) 405 END DO 406 END IF 407 END IF 408 IF (calculate_forces .AND. use_virial) THEN 409 virial%pv_ppl = virial%pv_virial - pv_loc 410 END IF 411 412 CALL timestop(handle) 413 414 END SUBROUTINE build_core_ae 415 416! ************************************************************************************************** 417!> \brief ... 418!> \param habd ... 419!> \param pab ... 420!> \param fa ... 421!> \param fb ... 422!> \param nder ... 423!> \param la_max ... 424!> \param la_min ... 425!> \param npgfa ... 426!> \param zeta ... 427!> \param lb_max ... 428!> \param lb_min ... 429!> \param npgfb ... 430!> \param zetb ... 431!> \param rab ... 432! ************************************************************************************************** 433 SUBROUTINE verfc_force(habd, pab, fa, fb, nder, la_max, la_min, npgfa, zeta, lb_max, lb_min, npgfb, zetb, rab) 434 435 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: habd, pab 436 REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: fa, fb 437 INTEGER, INTENT(IN) :: nder, la_max, la_min, npgfa 438 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta 439 INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb 440 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb 441 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab 442 443 CHARACTER(LEN=*), PARAMETER :: routineN = 'verfc_force', routineP = moduleN//':'//routineN 444 445 INTEGER :: ic_a, ic_b, icam1, icam2, icam3, icap1, & 446 icap2, icap3, icax, icbm1, icbm2, & 447 icbm3, icbx, icoa, icob, ipgfa, ipgfb, & 448 na, nap, nb 449 INTEGER, DIMENSION(3) :: la, lb 450 REAL(KIND=dp) :: zax2, zbx2 451 452 fa = 0.0_dp 453 fb = 0.0_dp 454 455 na = ncoset(la_max) 456 nap = ncoset(la_max + nder) 457 nb = ncoset(lb_max) 458 DO ipgfa = 1, npgfa 459 zax2 = zeta(ipgfa)*2.0_dp 460 DO ipgfb = 1, npgfb 461 zbx2 = zetb(ipgfb)*2.0_dp 462 DO ic_a = ncoset(la_min - 1) + 1, ncoset(la_max) 463 la(1:3) = indco(1:3, ic_a) 464 icap1 = coset(la(1) + 1, la(2), la(3)) 465 icap2 = coset(la(1), la(2) + 1, la(3)) 466 icap3 = coset(la(1), la(2), la(3) + 1) 467 icam1 = coset(la(1) - 1, la(2), la(3)) 468 icam2 = coset(la(1), la(2) - 1, la(3)) 469 icam3 = coset(la(1), la(2), la(3) - 1) 470 icoa = ic_a + (ipgfa - 1)*na 471 icax = (ipgfa - 1)*nap 472 473 DO ic_b = ncoset(lb_min - 1) + 1, ncoset(lb_max) 474 lb(1:3) = indco(1:3, ic_b) 475 icbm1 = coset(lb(1) - 1, lb(2), lb(3)) 476 icbm2 = coset(lb(1), lb(2) - 1, lb(3)) 477 icbm3 = coset(lb(1), lb(2), lb(3) - 1) 478 icob = ic_b + (ipgfb - 1)*nb 479 icbx = (ipgfb - 1)*nb 480 481 fa(1) = fa(1) - pab(icoa, icob)*(-zax2*habd(icap1 + icax, icob) + & 482 REAL(la(1), KIND=dp)*habd(icam1 + icax, icob)) 483 fa(2) = fa(2) - pab(icoa, icob)*(-zax2*habd(icap2 + icax, icob) + & 484 REAL(la(2), KIND=dp)*habd(icam2 + icax, icob)) 485 fa(3) = fa(3) - pab(icoa, icob)*(-zax2*habd(icap3 + icax, icob) + & 486 REAL(la(3), KIND=dp)*habd(icam3 + icax, icob)) 487 488 fb(1) = fb(1) - pab(icoa, icob)*( & 489 -zbx2*(habd(icap1 + icax, icob) - rab(1)*habd(ic_a + icax, icob)) + & 490 REAL(lb(1), KIND=dp)*habd(ic_a + icax, icbm1 + icbx)) 491 fb(2) = fb(2) - pab(icoa, icob)*( & 492 -zbx2*(habd(icap2 + icax, icob) - rab(2)*habd(ic_a + icax, icob)) + & 493 REAL(lb(2), KIND=dp)*habd(ic_a + icax, icbm2 + icbx)) 494 fb(3) = fb(3) - pab(icoa, icob)*( & 495 -zbx2*(habd(icap3 + icax, icob) - rab(3)*habd(ic_a + icax, icob)) + & 496 REAL(lb(3), KIND=dp)*habd(ic_a + icax, icbm3 + icbx)) 497 498 END DO ! ic_b 499 END DO ! ic_a 500 END DO ! ipgfb 501 END DO ! ipgfa 502 503 END SUBROUTINE verfc_force 504 505! ************************************************************************************************** 506 507END MODULE core_ae 508