1060 pw_pool, weights, xc_section, &
1061 do_triplet, calc_virial, virial_xc, deriv_set)
1063 TYPE(
pw_r3d_rs_type),
DIMENSION(:),
INTENT(IN),
POINTER :: v_xc, v_tau
1065 TYPE(
pw_r3d_rs_type),
DIMENSION(:),
INTENT(IN),
POINTER :: rho1_r, tau1_r
1066 TYPE(
pw_c1d_gs_type),
DIMENSION(:),
INTENT(IN),
POINTER :: rho1_g
1070 LOGICAL,
INTENT(IN) :: do_triplet
1071 LOGICAL,
INTENT(IN),
OPTIONAL :: calc_virial
1072 REAL(kind=
dp),
DIMENSION(3, 3),
INTENT(INOUT), &
1073 OPTIONAL :: virial_xc
1076 CHARACTER(len=*),
PARAMETER :: routinen =
'xc_calc_2nd_deriv_numerical'
1077 REAL(kind=
dp),
DIMENSION(-4:4, 4),
PARAMETER :: &
1078 rweights = reshape([0.0_dp, 0.0_dp, 0.0_dp, -0.5_dp, 0.0_dp, 0.5_dp, 0.0_dp, 0.0_dp, 0.0_dp, &
1079 0.0_dp, 0.0_dp, 1.0_dp/12.0_dp, -2.0_dp/3.0_dp, 0.0_dp, 2.0_dp/3.0_dp, -1.0_dp/12.0_dp, 0.0_dp, 0.0_dp, &
1080 0.0_dp, -1.0_dp/60.0_dp, 0.15_dp, -0.75_dp, 0.0_dp, 0.75_dp, -0.15_dp, 1.0_dp/60.0_dp, 0.0_dp, &
1081 1.0_dp/280.0_dp, -4.0_dp/105.0_dp, 0.2_dp, -0.8_dp, 0.0_dp, 0.8_dp, -0.2_dp, 4.0_dp/105.0_dp, -1.0_dp/280.0_dp], [9, 4])
1083 INTEGER :: handle, idir, ispin, nspins, istep, nsteps
1084 INTEGER,
DIMENSION(2, 3) :: bo
1085 LOGICAL :: gradient_f, lsd, my_calc_virial, tau_f, laplace_f, rho_f
1086 REAL(kind=
dp) :: exc, gradient_cut, h, rweight, step, rho_cutoff
1087 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: dr1dr, dra1dra, drb1drb
1088 REAL(kind=
dp),
DIMENSION(3, 3) :: virial_dummy
1089 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: norm_drho, norm_drho2, norm_drho2a, &
1090 norm_drho2b, norm_drhoa, norm_drhob, &
1091 rho, rho1, rho1a, rho1b, rhoa, rhob, &
1092 tau_a, tau_b, tau, tau1, tau1a, tau1b, laplace, laplace1, &
1093 laplacea, laplaceb, laplace1a, laplace1b, &
1094 laplace2, laplace2a, laplace2b, deriv_data
1095 TYPE(
cp_3d_r_cp_type),
DIMENSION(3) :: drho, drho1, drho1a, drho1b, drhoa, drhob
1100 TYPE(
pw_r3d_rs_type) :: virial_pw, v_laplace, v_laplacea, v_laplaceb
1106 CALL timeset(routinen, handle)
1108 my_calc_virial = .false.
1109 IF (
PRESENT(calc_virial) .AND.
PRESENT(virial_xc)) my_calc_virial = calc_virial
1113 NULLIFY (tau, tau_r, tau_a, tau_b)
1117 IF (nsteps < lbound(rweights, 2) .OR. nspins > ubound(rweights, 2))
THEN
1118 cpabort(
"The number of steps must be a value from 1 to 4.")
1121 IF (nspins == 2)
THEN
1122 NULLIFY (vxc_rho, rho_g, vxc_tau)
1124 DO ispin = 1, nspins
1125 CALL pw_pool%create_pw(rho_r(ispin))
1127 IF (
ASSOCIATED(tau1_r) .AND.
ASSOCIATED(v_tau))
THEN
1129 DO ispin = 1, nspins
1130 CALL pw_pool%create_pw(tau_r(ispin))
1133 CALL xc_rho_set_get(rho_set, can_return_null=.true., rhoa=rhoa, rhob=rhob, tau_a=tau_a, tau_b=tau_b)
1134 DO istep = -nsteps, nsteps
1135 IF (istep == 0) cycle
1136 rweight = rweights(istep, nsteps)/h
1137 step = real(istep,
dp)*h
1138 CALL calc_resp_potential_numer_ab(rho_r, rho_g, rho1_r, rhoa, rhob, vxc_rho, &
1139 tau_r, tau1_r, tau_a, tau_b, vxc_tau, xc_section, &
1140 weights, pw_pool, step)
1141 DO ispin = 1, nspins
1142 CALL pw_axpy(vxc_rho(ispin), v_xc(ispin), rweight)
1143 IF (
ASSOCIATED(vxc_tau) .AND.
ASSOCIATED(v_tau))
THEN
1144 CALL pw_axpy(vxc_tau(ispin), v_tau(ispin), rweight)
1147 DO ispin = 1, nspins
1148 CALL vxc_rho(ispin)%release()
1150 DEALLOCATE (vxc_rho)
1151 IF (
ASSOCIATED(vxc_tau))
THEN
1152 DO ispin = 1, nspins
1153 CALL vxc_tau(ispin)%release()
1155 DEALLOCATE (vxc_tau)
1158 ELSE IF (nspins == 1 .AND. do_triplet)
THEN
1159 NULLIFY (vxc_rho, vxc_tau, rho_g)
1162 CALL pw_pool%create_pw(rho_r(ispin))
1164 IF (
ASSOCIATED(tau1_r) .AND.
ASSOCIATED(v_tau))
THEN
1166 DO ispin = 1, nspins
1167 CALL pw_pool%create_pw(tau_r(ispin))
1170 CALL xc_rho_set_get(rho_set, can_return_null=.true., rhoa=rhoa, rhob=rhob, tau_a=tau_a, tau_b=tau_b)
1171 DO istep = -nsteps, nsteps
1172 IF (istep == 0) cycle
1173 rweight = rweights(istep, nsteps)/h
1174 step = real(istep,
dp)*h
1178 rho_r(1)%array(:, :, :) = rhoa(:, :, :) + step*rho1_r(1)%array(:, :, :)
1181 rho_r(2)%array(:, :, :) = rhob(:, :, :)
1183 IF (
ASSOCIATED(tau1_r))
THEN
1185 tau_r(1)%array(:, :, :) = tau_a(:, :, :) + step*tau1_r(1)%array(:, :, :)
1188 tau_r(2)%array(:, :, :) = tau_b(:, :, :)
1192 CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
1193 weights, pw_pool, .false., virial_dummy)
1194 CALL pw_axpy(vxc_rho(1), v_xc(1), rweight)
1195 IF (
ASSOCIATED(vxc_tau) .AND.
ASSOCIATED(v_tau))
THEN
1196 CALL pw_axpy(vxc_tau(1), v_tau(1), rweight)
1199 CALL vxc_rho(ispin)%release()
1201 DEALLOCATE (vxc_rho)
1202 IF (
ASSOCIATED(vxc_tau))
THEN
1204 CALL vxc_tau(ispin)%release()
1206 DEALLOCATE (vxc_tau)
1211 rho_r(1)%array(:, :, :) = rhoa(:, :, :)
1214 rho_r(2)%array(:, :, :) = rhob(:, :, :) + step*rho1_r(1)%array(:, :, :)
1216 IF (
ASSOCIATED(tau1_r))
THEN
1218 tau_r(1)%array(:, :, :) = tau_a(:, :, :)
1221 tau_r(2)%array(:, :, :) = tau_b(:, :, :) + step*tau1_r(1)%array(:, :, :)
1225 CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
1226 weights, pw_pool, .false., virial_dummy)
1227 CALL pw_axpy(vxc_rho(1), v_xc(1), rweight)
1228 IF (
ASSOCIATED(vxc_tau) .AND.
ASSOCIATED(v_tau))
THEN
1229 CALL pw_axpy(vxc_tau(1), v_tau(1), rweight)
1232 CALL vxc_rho(ispin)%release()
1234 DEALLOCATE (vxc_rho)
1235 IF (
ASSOCIATED(vxc_tau))
THEN
1237 CALL vxc_tau(ispin)%release()
1239 DEALLOCATE (vxc_tau)
1243 NULLIFY (vxc_rho, rho_r, rho_g, vxc_tau, tau_r, tau)
1245 CALL pw_pool%create_pw(rho_r(1))
1246 IF (
ASSOCIATED(tau1_r) .AND.
ASSOCIATED(v_tau))
THEN
1248 CALL pw_pool%create_pw(tau_r(1))
1250 CALL xc_rho_set_get(rho_set, can_return_null=.true., rho=rho, tau=tau)
1251 DO istep = -nsteps, nsteps
1252 IF (istep == 0) cycle
1253 rweight = rweights(istep, nsteps)/h
1254 step = real(istep,
dp)*h
1257 rho_r(1)%array(:, :, :) = rho(:, :, :) + step*rho1_r(1)%array(:, :, :)
1259 IF (
ASSOCIATED(tau1_r) .AND.
ASSOCIATED(tau) .AND.
ASSOCIATED(tau1_r))
THEN
1261 tau_r(1)%array(:, :, :) = tau(:, :, :) + step*tau1_r(1)%array(:, :, :)
1265 CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
1266 weights, pw_pool, .false., virial_dummy)
1267 CALL pw_axpy(vxc_rho(1), v_xc(1), rweight)
1268 IF (
ASSOCIATED(vxc_tau) .AND.
ASSOCIATED(v_tau))
THEN
1269 CALL pw_axpy(vxc_tau(1), v_tau(1), rweight)
1271 CALL vxc_rho(1)%release()
1272 DEALLOCATE (vxc_rho)
1273 IF (
ASSOCIATED(vxc_tau))
THEN
1274 CALL vxc_tau(1)%release()
1275 DEALLOCATE (vxc_tau)
1280 IF (my_calc_virial)
THEN
1282 IF (nspins == 1 .AND. do_triplet)
THEN
1286 CALL check_for_derivatives(deriv_set, (nspins == 2), rho_f, gradient_f, tau_f, laplace_f)
1292 IF (gradient_f)
THEN
1293 bo = rho_set%local_bounds
1296 CALL allocate_pw(virial_pw, pw_pool, bo)
1303 drho_cutoff=gradient_cut, &
1316 CALL xc_rho_set_get(rho_set, drhoa=drhoa, drhob=drhob, norm_drho=norm_drho, &
1317 norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, tau_a=tau_a, tau_b=tau_b, &
1318 laplace_rhoa=laplacea, laplace_rhob=laplaceb, can_return_null=.true.)
1319 CALL xc_rho_set_get(rho1_set, rhoa=rho1a, rhob=rho1b, drhoa=drho1a, drhob=drho1b, laplace_rhoa=laplace1a, &
1320 laplace_rhob=laplace1b, can_return_null=.true.)
1322 CALL calc_drho_from_ab(drho, drhoa, drhob)
1323 CALL calc_drho_from_ab(drho1, drho1a, drho1b)
1325 CALL xc_rho_set_get(rho_set, drho=drho, norm_drho=norm_drho, tau=tau, laplace_rho=laplace, can_return_null=.true.)
1326 CALL xc_rho_set_get(rho1_set, rho=rho1, drho=drho1, laplace_rho=laplace1, can_return_null=.true.)
1329 CALL prepare_dr1dr(dr1dr, drho, drho1)
1332 CALL prepare_dr1dr(dra1dra, drhoa, drho1a)
1333 CALL prepare_dr1dr(drb1drb, drhob, drho1b)
1335 CALL allocate_pw(v_drho, pw_pool, bo)
1336 CALL allocate_pw(v_drhoa, pw_pool, bo)
1337 CALL allocate_pw(v_drhob, pw_pool, bo)
1339 IF (
ASSOCIATED(norm_drhoa))
CALL apply_drho(deriv_set, [
deriv_norm_drhoa], virial_pw, &
1340 drhoa, drho1a, virial_xc, &
1341 norm_drhoa, gradient_cut, dra1dra, v_drhoa%array)
1342 IF (
ASSOCIATED(norm_drhob))
CALL apply_drho(deriv_set, [
deriv_norm_drhob], virial_pw, &
1343 drhob, drho1b, virial_xc, &
1344 norm_drhob, gradient_cut, drb1drb, v_drhob%array)
1345 IF (
ASSOCIATED(norm_drho))
CALL apply_drho(deriv_set, [
deriv_norm_drho], virial_pw, &
1346 drho, drho1, virial_xc, &
1347 norm_drho, gradient_cut, dr1dr, v_drho%array)
1350 cpassert(
ASSOCIATED(deriv_data))
1351 virial_pw%array(:, :, :) = -rho1a(:, :, :)
1352 CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
1354 CALL allocate_pw(v_laplacea, pw_pool, bo)
1357 cpassert(
ASSOCIATED(deriv_data))
1358 virial_pw%array(:, :, :) = -rho1b(:, :, :)
1359 CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
1361 CALL allocate_pw(v_laplaceb, pw_pool, bo)
1367 CALL allocate_pw(v_drho, pw_pool, bo)
1369 CALL apply_drho(deriv_set, [
deriv_norm_drho], virial_pw, drho, drho1, virial_xc, &
1370 norm_drho, gradient_cut, dr1dr, v_drho%array)
1373 cpassert(
ASSOCIATED(deriv_data))
1374 virial_pw%array(:, :, :) = -rho1(:, :, :)
1375 CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
1377 CALL allocate_pw(v_laplace, pw_pool, bo)
1383 rho_r(1)%array = rhoa
1384 rho_r(2)%array = rhob
1386 rho_r(1)%array = rho
1388 IF (
ASSOCIATED(tau1_r))
THEN
1390 tau_r(1)%array = tau_a
1391 tau_r(2)%array = tau_b
1393 tau_r(1)%array = tau
1404 rho_cutoff=rho_cutoff, &
1415 CALL xc_rho_set_get(rho1_set, rhoa=rho1a, rhob=rho1b, tau_a=tau1a, tau_b=tau1b, &
1416 laplace_rhoa=laplace1a, laplace_rhob=laplace1b, can_return_null=.true.)
1417 CALL xc_rho_set_get(rho2_set, norm_drhoa=norm_drho2a, norm_drhob=norm_drho2b, &
1418 norm_drho=norm_drho2, laplace_rhoa=laplace2a, laplace_rhob=laplace2b, can_return_null=.true.)
1420 DO istep = -nsteps, nsteps
1421 IF (istep == 0) cycle
1422 rweight = rweights(istep, nsteps)/h
1423 step = real(istep,
dp)*h
1424 IF (
ASSOCIATED(norm_drhoa))
THEN
1425 CALL get_derivs_rho(norm_drho2a, norm_drhoa, step, xc_fun_section, lsd, rho2_set, deriv_set1)
1426 CALL update_deriv_rho(deriv_set1, [
deriv_rhoa], bo, &
1427 norm_drhoa, gradient_cut, rweight, rho1a, v_drhoa%array)
1428 CALL update_deriv_rho(deriv_set1, [
deriv_rhob], bo, &
1429 norm_drhoa, gradient_cut, rweight, rho1b, v_drhoa%array)
1431 norm_drhoa, gradient_cut, rweight, dra1dra, v_drhoa%array)
1433 norm_drhoa, gradient_cut, rweight, dra1dra, drb1drb, v_drhoa%array, v_drhob%array)
1435 norm_drhoa, gradient_cut, rweight, dra1dra, dr1dr, v_drhoa%array, v_drho%array)
1437 CALL update_deriv_rho(deriv_set1, [
deriv_tau_a], bo, &
1438 norm_drhoa, gradient_cut, rweight, tau1a, v_drhoa%array)
1439 CALL update_deriv_rho(deriv_set1, [
deriv_tau_b], bo, &
1440 norm_drhoa, gradient_cut, rweight, tau1b, v_drhoa%array)
1444 norm_drhoa, gradient_cut, rweight, laplace1a, v_drhoa%array)
1446 norm_drhoa, gradient_cut, rweight, laplace1b, v_drhoa%array)
1450 IF (
ASSOCIATED(norm_drhob))
THEN
1451 CALL get_derivs_rho(norm_drho2b, norm_drhob, step, xc_fun_section, lsd, rho2_set, deriv_set1)
1452 CALL update_deriv_rho(deriv_set1, [
deriv_rhoa], bo, &
1453 norm_drhob, gradient_cut, rweight, rho1a, v_drhob%array)
1454 CALL update_deriv_rho(deriv_set1, [
deriv_rhob], bo, &
1455 norm_drhob, gradient_cut, rweight, rho1b, v_drhob%array)
1457 norm_drhob, gradient_cut, rweight, drb1drb, v_drhob%array)
1459 norm_drhob, gradient_cut, rweight, drb1drb, dra1dra, v_drhob%array, v_drhoa%array)
1461 norm_drhob, gradient_cut, rweight, drb1drb, dr1dr, v_drhob%array, v_drho%array)
1463 CALL update_deriv_rho(deriv_set1, [
deriv_tau_a], bo, &
1464 norm_drhob, gradient_cut, rweight, tau1a, v_drhob%array)
1465 CALL update_deriv_rho(deriv_set1, [
deriv_tau_b], bo, &
1466 norm_drhob, gradient_cut, rweight, tau1b, v_drhob%array)
1470 norm_drhob, gradient_cut, rweight, laplace1a, v_drhob%array)
1472 norm_drhob, gradient_cut, rweight, laplace1b, v_drhob%array)
1476 IF (
ASSOCIATED(norm_drho))
THEN
1477 CALL get_derivs_rho(norm_drho2, norm_drho, step, xc_fun_section, lsd, rho2_set, deriv_set1)
1478 CALL update_deriv_rho(deriv_set1, [
deriv_rhoa], bo, &
1479 norm_drho, gradient_cut, rweight, rho1a, v_drho%array)
1480 CALL update_deriv_rho(deriv_set1, [
deriv_rhob], bo, &
1481 norm_drho, gradient_cut, rweight, rho1b, v_drho%array)
1483 norm_drho, gradient_cut, rweight, dr1dr, v_drho%array)
1485 norm_drho, gradient_cut, rweight, dr1dr, dra1dra, v_drho%array, v_drhoa%array)
1487 norm_drho, gradient_cut, rweight, dr1dr, drb1drb, v_drho%array, v_drhob%array)
1489 CALL update_deriv_rho(deriv_set1, [
deriv_tau_a], bo, &
1490 norm_drho, gradient_cut, rweight, tau1a, v_drho%array)
1491 CALL update_deriv_rho(deriv_set1, [
deriv_tau_b], bo, &
1492 norm_drho, gradient_cut, rweight, tau1b, v_drho%array)
1496 norm_drho, gradient_cut, rweight, laplace1a, v_drho%array)
1498 norm_drho, gradient_cut, rweight, laplace1b, v_drho%array)
1504 CALL get_derivs_rho(laplace2a, laplacea, step, xc_fun_section, lsd, rho2_set, deriv_set1)
1507 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_rhoa], bo, &
1508 rweight, rho1a, v_laplacea%array)
1509 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_rhob], bo, &
1510 rweight, rho1b, v_laplacea%array)
1511 IF (
ASSOCIATED(norm_drho))
THEN
1512 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_norm_drho], bo, &
1513 rweight, dr1dr, v_laplacea%array)
1515 IF (
ASSOCIATED(norm_drhoa))
THEN
1516 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_norm_drhoa], bo, &
1517 rweight, dra1dra, v_laplacea%array)
1519 IF (
ASSOCIATED(norm_drhob))
THEN
1520 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_norm_drhob], bo, &
1521 rweight, drb1drb, v_laplacea%array)
1524 IF (
ASSOCIATED(tau1a))
THEN
1525 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_tau_a], bo, &
1526 rweight, tau1a, v_laplacea%array)
1528 IF (
ASSOCIATED(tau1b))
THEN
1529 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_tau_b], bo, &
1530 rweight, tau1b, v_laplacea%array)
1534 rweight, laplace1a, v_laplacea%array)
1537 rweight, laplace1b, v_laplacea%array)
1540 CALL get_derivs_rho(laplace2b, laplaceb, step, xc_fun_section, lsd, rho2_set, deriv_set1)
1543 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_rhoa], bo, &
1544 rweight, rho1a, v_laplaceb%array)
1545 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_rhob], bo, &
1546 rweight, rho1b, v_laplaceb%array)
1547 IF (
ASSOCIATED(norm_drho))
THEN
1548 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_norm_drho], bo, &
1549 rweight, dr1dr, v_laplaceb%array)
1551 IF (
ASSOCIATED(norm_drhoa))
THEN
1552 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_norm_drhoa], bo, &
1553 rweight, dra1dra, v_laplaceb%array)
1555 IF (
ASSOCIATED(norm_drhob))
THEN
1556 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_norm_drhob], bo, &
1557 rweight, drb1drb, v_laplaceb%array)
1561 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_tau_a], bo, &
1562 rweight, tau1a, v_laplaceb%array)
1563 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_tau_b], bo, &
1564 rweight, tau1b, v_laplaceb%array)
1568 rweight, laplace1a, v_laplaceb%array)
1571 rweight, laplace1b, v_laplaceb%array)
1575 CALL virial_drho_drho(virial_pw, drhoa, v_drhoa, virial_xc)
1576 CALL virial_drho_drho(virial_pw, drhob, v_drhob, virial_xc)
1577 CALL virial_drho_drho(virial_pw, drho, v_drho, virial_xc)
1579 CALL deallocate_pw(v_drho, pw_pool)
1580 CALL deallocate_pw(v_drhoa, pw_pool)
1581 CALL deallocate_pw(v_drhob, pw_pool)
1584 virial_pw%array(:, :, :) = -rhoa(:, :, :)
1585 CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplacea%array)
1586 CALL deallocate_pw(v_laplacea, pw_pool)
1588 virial_pw%array(:, :, :) = -rhob(:, :, :)
1589 CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplaceb%array)
1590 CALL deallocate_pw(v_laplaceb, pw_pool)
1593 CALL deallocate_pw(virial_pw, pw_pool)
1596 DEALLOCATE (drho(idir)%array)
1597 DEALLOCATE (drho1(idir)%array)
1599 DEALLOCATE (dra1dra, drb1drb)
1602 CALL xc_rho_set_get(rho1_set, rho=rho1, tau=tau1, laplace_rho=laplace1, can_return_null=.true.)
1603 CALL xc_rho_set_get(rho2_set, norm_drho=norm_drho2, laplace_rho=laplace2, can_return_null=.true.)
1605 DO istep = -nsteps, nsteps
1606 IF (istep == 0) cycle
1607 rweight = rweights(istep, nsteps)/h
1608 step = real(istep,
dp)*h
1609 CALL get_derivs_rho(norm_drho2, norm_drho, step, xc_fun_section, lsd, rho2_set, deriv_set1)
1612 CALL update_deriv_rho(deriv_set1, [
deriv_rho], bo, &
1613 norm_drho, gradient_cut, rweight, rho1, v_drho%array)
1615 norm_drho, gradient_cut, rweight, dr1dr, v_drho%array)
1618 CALL update_deriv_rho(deriv_set1, [
deriv_tau], bo, &
1619 norm_drho, gradient_cut, rweight, tau1, v_drho%array)
1623 norm_drho, gradient_cut, rweight, laplace1, v_drho%array)
1625 CALL get_derivs_rho(laplace2, laplace, step, xc_fun_section, lsd, rho2_set, deriv_set1)
1628 CALL update_deriv(deriv_set1, laplace, rho_cutoff, [
deriv_rho], bo, &
1629 rweight, rho1, v_laplace%array)
1630 CALL update_deriv(deriv_set1, laplace, rho_cutoff, [
deriv_norm_drho], bo, &
1631 rweight, dr1dr, v_laplace%array)
1634 CALL update_deriv(deriv_set1, laplace, rho_cutoff, [
deriv_tau], bo, &
1635 rweight, tau1, v_laplace%array)
1639 rweight, laplace1, v_laplace%array)
1644 CALL virial_drho_drho(virial_pw, drho, v_drho, virial_xc)
1646 CALL deallocate_pw(v_drho, pw_pool)
1649 virial_pw%array(:, :, :) = -rho(:, :, :)
1650 CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace%array)
1651 CALL deallocate_pw(v_laplace, pw_pool)
1654 CALL deallocate_pw(virial_pw, pw_pool)
1667 DO ispin = 1,
SIZE(rho_r)
1668 CALL pw_pool%give_back_pw(rho_r(ispin))
1672 IF (
ASSOCIATED(tau_r))
THEN
1673 DO ispin = 1,
SIZE(tau_r)
1674 CALL pw_pool%give_back_pw(tau_r(ispin))
1679 CALL timestop(handle)
2054 LOGICAL,
INTENT(IN),
OPTIONAL :: gapw
2055 REAL(kind=
dp),
DIMENSION(:, :, :, :),
OPTIONAL, &
2057 REAL(kind=
dp),
INTENT(in),
OPTIONAL :: tddfpt_fac
2058 LOGICAL,
INTENT(IN),
OPTIONAL :: compute_virial
2059 REAL(kind=
dp),
DIMENSION(3, 3),
INTENT(INOUT), &
2060 OPTIONAL :: virial_xc
2061 LOGICAL,
INTENT(in),
OPTIONAL :: spinflip
2063 CHARACTER(len=*),
PARAMETER :: routinen =
'xc_calc_2nd_deriv_analytical'
2065 INTEGER :: handle, i, ia, idir, ir, ispin, j, jdir, &
2066 k, nspins, xc_deriv_method_id
2067 INTEGER,
DIMENSION(2, 3) :: bo
2068 LOGICAL :: gradient_f, lsd, my_compute_virial, alda0, &
2069 my_gapw, tau_f, laplace_f, rho_f, do_spinflip
2070 REAL(kind=
dp) ::
fac, gradient_cut, tmp, factor2, s, s_thresh
2071 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: dr1dr, dra1dra, drb1drb
2072 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: deriv_data, deriv_data2, &
2073 e_drhoa, e_drhob, e_drho, norm_drho, norm_drhoa, &
2074 norm_drhob, rho1, rho1a, rho1b, &
2075 tau1, tau1a, tau1b, laplace1, laplace1a, laplace1b, &
2077 TYPE(
cp_3d_r_cp_type),
DIMENSION(3) :: drho, drho1, drho1a, drho1b, drhoa, drhob
2078 TYPE(
pw_r3d_rs_type),
DIMENSION(:),
ALLOCATABLE :: v_drhoa, v_drhob, v_drho, v_laplace
2084 CALL timeset(routinen, handle)
2086 NULLIFY (e_drhoa, e_drhob, e_drho)
2089 IF (
PRESENT(gapw)) my_gapw = gapw
2091 my_compute_virial = .false.
2092 IF (
PRESENT(compute_virial)) my_compute_virial = compute_virial
2094 cpassert(
ASSOCIATED(v_xc))
2095 cpassert(
ASSOCIATED(xc_section))
2097 cpassert(
PRESENT(vxg))
2099 IF (my_compute_virial)
THEN
2100 cpassert(
PRESENT(virial_xc))
2104 i_val=xc_deriv_method_id)
2107 lsd =
ASSOCIATED(rho_set%rhoa)
2110 IF (
PRESENT(tddfpt_fac))
fac = tddfpt_fac
2111 IF (
PRESENT(tddfpt_fac)) factor2 = tddfpt_fac
2112 do_spinflip = .false.
2113 IF (
PRESENT(spinflip)) do_spinflip = spinflip
2115 bo = rho_set%local_bounds
2117 CALL check_for_derivatives(deriv_set, lsd, rho_f, gradient_f, tau_f, laplace_f)
2120 IF (gradient_f)
THEN
2127 cpassert(
ASSOCIATED(v_xc_tau))
2130 IF (gradient_f)
THEN
2131 ALLOCATE (v_drho_r(3, nspins), v_drho(nspins))
2132 DO ispin = 1, nspins
2134 CALL allocate_pw(v_drho_r(idir, ispin), pw_pool, bo)
2136 CALL allocate_pw(v_drho(ispin), pw_pool, bo)
2140 IF (
ASSOCIATED(pw_pool))
THEN
2141 CALL pw_pool%create_pw(tmp_g)
2142 CALL pw_pool%create_pw(vxc_g)
2145 cpabort(
"XC_DERIV method is not implemented in GAPW")
2150 DO ispin = 1, nspins
2151 v_xc(ispin)%array = 0.0_dp
2155 DO ispin = 1, nspins
2156 v_xc_tau(ispin)%array = 0.0_dp
2160 IF (laplace_f .AND. my_gapw)
THEN
2161 cpabort(
"Laplace-dependent functional not implemented with GAPW!")
2164 IF (my_compute_virial .AND. (gradient_f .OR. laplace_f))
CALL allocate_pw(virial_pw, pw_pool, bo)
2172 IF (do_spinflip)
THEN
2179 IF (gradient_f)
THEN
2181 norm_drho=norm_drho, norm_drhoa=norm_drhoa, norm_drhob=norm_drhob)
2182 IF (do_spinflip)
THEN
2184 CALL calc_drho_from_a(drho1, drho1a)
2187 CALL calc_drho_from_ab(drho1, drho1a, drho1b)
2190 CALL calc_drho_from_ab(drho, drhoa, drhob)
2192 CALL prepare_dr1dr(dra1dra, drhoa, drho1a)
2193 IF (do_spinflip)
THEN
2194 CALL prepare_dr1dr(drb1drb, drhob, drho1a)
2195 CALL prepare_dr1dr(dr1dr, drho, drho1a)
2196 ELSE IF (nspins /= 1)
THEN
2197 CALL prepare_dr1dr(drb1drb, drhob, drho1b)
2198 CALL prepare_dr1dr(dr1dr, drho, drho1)
2200 CALL prepare_dr1dr(drb1drb, drhob, drho1b)
2201 CALL prepare_dr1dr_ab(dr1dr, drhoa, drhob, drho1a, drho1b,
fac)
2204 ALLOCATE (v_drhoa(nspins), v_drhob(nspins))
2205 DO ispin = 1, nspins
2206 CALL allocate_pw(v_drhoa(ispin), pw_pool, bo)
2207 CALL allocate_pw(v_drhob(ispin), pw_pool, bo)
2213 CALL xc_rho_set_get(rho1_set, laplace_rhoa=laplace1a, laplace_rhob=laplace1b)
2215 ALLOCATE (v_laplace(nspins))
2216 DO ispin = 1, nspins
2217 CALL allocate_pw(v_laplace(ispin), pw_pool, bo)
2220 IF (my_compute_virial)
CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob)
2227 IF (do_spinflip)
THEN
2236 IF (
ASSOCIATED(deriv_att))
THEN
2239 IF (
ASSOCIATED(deriv_att))
THEN
2243 DO k = bo(1, 3), bo(2, 3)
2244 DO j = bo(1, 2), bo(2, 2)
2245 DO i = bo(1, 1), bo(2, 1)
2246 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
2247 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2248 (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)/s
2260 IF (.NOT. alda0)
THEN
2262 IF (
ASSOCIATED(deriv_att))
THEN
2265 IF (
ASSOCIATED(deriv_att))
THEN
2269 DO k = bo(1, 3), bo(2, 3)
2270 DO j = bo(1, 2), bo(2, 2)
2271 DO i = bo(1, 1), bo(2, 1)
2272 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
2273 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2274 (deriv_data(i, j, k)*dra1dra(i, j, k) - &
2275 deriv_data2(i, j, k)*drb1drb(i, j, k))/s
2284 ELSE IF (nspins /= 1)
THEN
2288 IF (
ASSOCIATED(deriv_att))
THEN
2292 DO k = bo(1, 3), bo(2, 3)
2293 DO j = bo(1, 2), bo(2, 2)
2294 DO i = bo(1, 1), bo(2, 1)
2295 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2296 deriv_data(i, j, k)*rho1a(i, j, k)
2302 IF (
ASSOCIATED(deriv_att))
THEN
2306 DO k = bo(1, 3), bo(2, 3)
2307 DO j = bo(1, 2), bo(2, 2)
2308 DO i = bo(1, 1), bo(2, 1)
2309 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2310 deriv_data(i, j, k)*rho1b(i, j, k)
2316 IF (
ASSOCIATED(deriv_att))
THEN
2320 DO k = bo(1, 3), bo(2, 3)
2321 DO j = bo(1, 2), bo(2, 2)
2322 DO i = bo(1, 1), bo(2, 1)
2323 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2324 deriv_data(i, j, k)*dr1dr(i, j, k)
2330 IF (
ASSOCIATED(deriv_att))
THEN
2334 DO k = bo(1, 3), bo(2, 3)
2335 DO j = bo(1, 2), bo(2, 2)
2336 DO i = bo(1, 1), bo(2, 1)
2337 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2338 deriv_data(i, j, k)*dra1dra(i, j, k)
2344 IF (
ASSOCIATED(deriv_att))
THEN
2348 DO k = bo(1, 3), bo(2, 3)
2349 DO j = bo(1, 2), bo(2, 2)
2350 DO i = bo(1, 1), bo(2, 1)
2351 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2352 deriv_data(i, j, k)*drb1drb(i, j, k)
2358 IF (
ASSOCIATED(deriv_att))
THEN
2362 DO k = bo(1, 3), bo(2, 3)
2363 DO j = bo(1, 2), bo(2, 2)
2364 DO i = bo(1, 1), bo(2, 1)
2365 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2366 deriv_data(i, j, k)*tau1a(i, j, k)
2372 IF (
ASSOCIATED(deriv_att))
THEN
2376 DO k = bo(1, 3), bo(2, 3)
2377 DO j = bo(1, 2), bo(2, 2)
2378 DO i = bo(1, 1), bo(2, 1)
2379 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2380 deriv_data(i, j, k)*tau1b(i, j, k)
2386 IF (
ASSOCIATED(deriv_att))
THEN
2390 DO k = bo(1, 3), bo(2, 3)
2391 DO j = bo(1, 2), bo(2, 2)
2392 DO i = bo(1, 1), bo(2, 1)
2393 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2394 deriv_data(i, j, k)*laplace1a(i, j, k)
2400 IF (
ASSOCIATED(deriv_att))
THEN
2404 DO k = bo(1, 3), bo(2, 3)
2405 DO j = bo(1, 2), bo(2, 2)
2406 DO i = bo(1, 1), bo(2, 1)
2407 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2408 deriv_data(i, j, k)*laplace1b(i, j, k)
2416 IF (
ASSOCIATED(deriv_att))
THEN
2420 DO k = bo(1, 3), bo(2, 3)
2421 DO j = bo(1, 2), bo(2, 2)
2422 DO i = bo(1, 1), bo(2, 1)
2423 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2424 deriv_data(i, j, k)*rho1a(i, j, k)
2430 IF (
ASSOCIATED(deriv_att))
THEN
2434 DO k = bo(1, 3), bo(2, 3)
2435 DO j = bo(1, 2), bo(2, 2)
2436 DO i = bo(1, 1), bo(2, 1)
2437 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2438 deriv_data(i, j, k)*rho1b(i, j, k)
2444 IF (
ASSOCIATED(deriv_att))
THEN
2448 DO k = bo(1, 3), bo(2, 3)
2449 DO j = bo(1, 2), bo(2, 2)
2450 DO i = bo(1, 1), bo(2, 1)
2451 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2452 deriv_data(i, j, k)*dr1dr(i, j, k)
2458 IF (
ASSOCIATED(deriv_att))
THEN
2462 DO k = bo(1, 3), bo(2, 3)
2463 DO j = bo(1, 2), bo(2, 2)
2464 DO i = bo(1, 1), bo(2, 1)
2465 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2466 deriv_data(i, j, k)*dra1dra(i, j, k)
2472 IF (
ASSOCIATED(deriv_att))
THEN
2476 DO k = bo(1, 3), bo(2, 3)
2477 DO j = bo(1, 2), bo(2, 2)
2478 DO i = bo(1, 1), bo(2, 1)
2479 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2480 deriv_data(i, j, k)*drb1drb(i, j, k)
2486 IF (
ASSOCIATED(deriv_att))
THEN
2490 DO k = bo(1, 3), bo(2, 3)
2491 DO j = bo(1, 2), bo(2, 2)
2492 DO i = bo(1, 1), bo(2, 1)
2493 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2494 deriv_data(i, j, k)*tau1a(i, j, k)
2500 IF (
ASSOCIATED(deriv_att))
THEN
2504 DO k = bo(1, 3), bo(2, 3)
2505 DO j = bo(1, 2), bo(2, 2)
2506 DO i = bo(1, 1), bo(2, 1)
2507 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2508 deriv_data(i, j, k)*tau1b(i, j, k)
2514 IF (
ASSOCIATED(deriv_att))
THEN
2518 DO k = bo(1, 3), bo(2, 3)
2519 DO j = bo(1, 2), bo(2, 2)
2520 DO i = bo(1, 1), bo(2, 1)
2521 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2522 deriv_data(i, j, k)*laplace1a(i, j, k)
2528 IF (
ASSOCIATED(deriv_att))
THEN
2532 DO k = bo(1, 3), bo(2, 3)
2533 DO j = bo(1, 2), bo(2, 2)
2534 DO i = bo(1, 1), bo(2, 1)
2535 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2536 deriv_data(i, j, k)*laplace1b(i, j, k)
2544 IF (
ASSOCIATED(deriv_att))
THEN
2548 DO k = bo(1, 3), bo(2, 3)
2549 DO j = bo(1, 2), bo(2, 2)
2550 DO i = bo(1, 1), bo(2, 1)
2551 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2552 deriv_data(i, j, k)*rho1a(i, j, k)
2553 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2554 deriv_data(i, j, k)*rho1a(i, j, k)
2560 IF (
ASSOCIATED(deriv_att))
THEN
2564 DO k = bo(1, 3), bo(2, 3)
2565 DO j = bo(1, 2), bo(2, 2)
2566 DO i = bo(1, 1), bo(2, 1)
2567 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2568 deriv_data(i, j, k)*rho1b(i, j, k)
2569 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2570 deriv_data(i, j, k)*rho1b(i, j, k)
2576 IF (
ASSOCIATED(deriv_att))
THEN
2580 DO k = bo(1, 3), bo(2, 3)
2581 DO j = bo(1, 2), bo(2, 2)
2582 DO i = bo(1, 1), bo(2, 1)
2583 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2584 deriv_data(i, j, k)*dr1dr(i, j, k)
2585 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2586 deriv_data(i, j, k)*dr1dr(i, j, k)
2592 IF (
ASSOCIATED(deriv_att))
THEN
2596 DO k = bo(1, 3), bo(2, 3)
2597 DO j = bo(1, 2), bo(2, 2)
2598 DO i = bo(1, 1), bo(2, 1)
2599 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2600 deriv_data(i, j, k)*dra1dra(i, j, k)
2601 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2602 deriv_data(i, j, k)*dra1dra(i, j, k)
2608 IF (
ASSOCIATED(deriv_att))
THEN
2612 DO k = bo(1, 3), bo(2, 3)
2613 DO j = bo(1, 2), bo(2, 2)
2614 DO i = bo(1, 1), bo(2, 1)
2615 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2616 deriv_data(i, j, k)*drb1drb(i, j, k)
2617 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2618 deriv_data(i, j, k)*drb1drb(i, j, k)
2624 IF (
ASSOCIATED(deriv_att))
THEN
2628 DO k = bo(1, 3), bo(2, 3)
2629 DO j = bo(1, 2), bo(2, 2)
2630 DO i = bo(1, 1), bo(2, 1)
2631 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2632 deriv_data(i, j, k)*tau1a(i, j, k)
2633 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2634 deriv_data(i, j, k)*tau1a(i, j, k)
2640 IF (
ASSOCIATED(deriv_att))
THEN
2644 DO k = bo(1, 3), bo(2, 3)
2645 DO j = bo(1, 2), bo(2, 2)
2646 DO i = bo(1, 1), bo(2, 1)
2647 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2648 deriv_data(i, j, k)*tau1b(i, j, k)
2649 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2650 deriv_data(i, j, k)*tau1b(i, j, k)
2656 IF (
ASSOCIATED(deriv_att))
THEN
2660 DO k = bo(1, 3), bo(2, 3)
2661 DO j = bo(1, 2), bo(2, 2)
2662 DO i = bo(1, 1), bo(2, 1)
2663 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2664 deriv_data(i, j, k)*laplace1a(i, j, k)
2665 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2666 deriv_data(i, j, k)*laplace1a(i, j, k)
2672 IF (
ASSOCIATED(deriv_att))
THEN
2676 DO k = bo(1, 3), bo(2, 3)
2677 DO j = bo(1, 2), bo(2, 2)
2678 DO i = bo(1, 1), bo(2, 1)
2679 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2680 deriv_data(i, j, k)*laplace1b(i, j, k)
2681 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2682 deriv_data(i, j, k)*laplace1b(i, j, k)
2689 IF (
ASSOCIATED(deriv_att))
THEN
2693 IF (my_compute_virial)
THEN
2694 CALL virial_drho_drho1(virial_pw, drho, drho1, deriv_data, virial_xc)
2698 v_drho(1)%array(:, :, :) = v_drho(1)%array(:, :, :) + &
2699 deriv_data(:, :, :)*dr1dr(:, :, :)/max(gradient_cut, norm_drho(:, :, :))**2
2700 v_drho(2)%array(:, :, :) = v_drho(2)%array(:, :, :) + &
2701 deriv_data(:, :, :)*dr1dr(:, :, :)/max(gradient_cut, norm_drho(:, :, :))**2
2706 IF (
ASSOCIATED(deriv_att))
THEN
2710 DO k = bo(1, 3), bo(2, 3)
2711 DO j = bo(1, 2), bo(2, 2)
2712 DO i = bo(1, 1), bo(2, 1)
2713 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2714 deriv_data(i, j, k)*rho1a(i, j, k)
2720 IF (
ASSOCIATED(deriv_att))
THEN
2724 DO k = bo(1, 3), bo(2, 3)
2725 DO j = bo(1, 2), bo(2, 2)
2726 DO i = bo(1, 1), bo(2, 1)
2727 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2728 deriv_data(i, j, k)*rho1b(i, j, k)
2734 IF (
ASSOCIATED(deriv_att))
THEN
2738 DO k = bo(1, 3), bo(2, 3)
2739 DO j = bo(1, 2), bo(2, 2)
2740 DO i = bo(1, 1), bo(2, 1)
2741 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2742 deriv_data(i, j, k)*dr1dr(i, j, k)
2748 IF (
ASSOCIATED(deriv_att))
THEN
2752 DO k = bo(1, 3), bo(2, 3)
2753 DO j = bo(1, 2), bo(2, 2)
2754 DO i = bo(1, 1), bo(2, 1)
2755 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2756 deriv_data(i, j, k)*dra1dra(i, j, k)
2762 IF (
ASSOCIATED(deriv_att))
THEN
2766 DO k = bo(1, 3), bo(2, 3)
2767 DO j = bo(1, 2), bo(2, 2)
2768 DO i = bo(1, 1), bo(2, 1)
2769 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2770 deriv_data(i, j, k)*drb1drb(i, j, k)
2776 IF (
ASSOCIATED(deriv_att))
THEN
2780 DO k = bo(1, 3), bo(2, 3)
2781 DO j = bo(1, 2), bo(2, 2)
2782 DO i = bo(1, 1), bo(2, 1)
2783 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2784 deriv_data(i, j, k)*tau1a(i, j, k)
2790 IF (
ASSOCIATED(deriv_att))
THEN
2794 DO k = bo(1, 3), bo(2, 3)
2795 DO j = bo(1, 2), bo(2, 2)
2796 DO i = bo(1, 1), bo(2, 1)
2797 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2798 deriv_data(i, j, k)*tau1b(i, j, k)
2804 IF (
ASSOCIATED(deriv_att))
THEN
2808 DO k = bo(1, 3), bo(2, 3)
2809 DO j = bo(1, 2), bo(2, 2)
2810 DO i = bo(1, 1), bo(2, 1)
2811 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2812 deriv_data(i, j, k)*laplace1a(i, j, k)
2818 IF (
ASSOCIATED(deriv_att))
THEN
2822 DO k = bo(1, 3), bo(2, 3)
2823 DO j = bo(1, 2), bo(2, 2)
2824 DO i = bo(1, 1), bo(2, 1)
2825 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2826 deriv_data(i, j, k)*laplace1b(i, j, k)
2833 IF (
ASSOCIATED(deriv_att))
THEN
2837 IF (my_compute_virial)
THEN
2838 CALL virial_drho_drho1(virial_pw, drhoa, drho1a, deriv_data, virial_xc)
2842 v_drhoa(1)%array(:, :, :) = v_drhoa(1)%array(:, :, :) + &
2843 deriv_data(:, :, :)*dra1dra(:, :, :)/max(gradient_cut, norm_drhoa(:, :, :))**2
2848 IF (
ASSOCIATED(deriv_att))
THEN
2852 DO k = bo(1, 3), bo(2, 3)
2853 DO j = bo(1, 2), bo(2, 2)
2854 DO i = bo(1, 1), bo(2, 1)
2855 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2856 deriv_data(i, j, k)*rho1a(i, j, k)
2862 IF (
ASSOCIATED(deriv_att))
THEN
2866 DO k = bo(1, 3), bo(2, 3)
2867 DO j = bo(1, 2), bo(2, 2)
2868 DO i = bo(1, 1), bo(2, 1)
2869 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2870 deriv_data(i, j, k)*rho1b(i, j, k)
2876 IF (
ASSOCIATED(deriv_att))
THEN
2880 DO k = bo(1, 3), bo(2, 3)
2881 DO j = bo(1, 2), bo(2, 2)
2882 DO i = bo(1, 1), bo(2, 1)
2883 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2884 deriv_data(i, j, k)*dr1dr(i, j, k)
2890 IF (
ASSOCIATED(deriv_att))
THEN
2894 DO k = bo(1, 3), bo(2, 3)
2895 DO j = bo(1, 2), bo(2, 2)
2896 DO i = bo(1, 1), bo(2, 1)
2897 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2898 deriv_data(i, j, k)*dra1dra(i, j, k)
2904 IF (
ASSOCIATED(deriv_att))
THEN
2908 DO k = bo(1, 3), bo(2, 3)
2909 DO j = bo(1, 2), bo(2, 2)
2910 DO i = bo(1, 1), bo(2, 1)
2911 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2912 deriv_data(i, j, k)*drb1drb(i, j, k)
2918 IF (
ASSOCIATED(deriv_att))
THEN
2922 DO k = bo(1, 3), bo(2, 3)
2923 DO j = bo(1, 2), bo(2, 2)
2924 DO i = bo(1, 1), bo(2, 1)
2925 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2926 deriv_data(i, j, k)*tau1a(i, j, k)
2932 IF (
ASSOCIATED(deriv_att))
THEN
2936 DO k = bo(1, 3), bo(2, 3)
2937 DO j = bo(1, 2), bo(2, 2)
2938 DO i = bo(1, 1), bo(2, 1)
2939 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2940 deriv_data(i, j, k)*tau1b(i, j, k)
2946 IF (
ASSOCIATED(deriv_att))
THEN
2950 DO k = bo(1, 3), bo(2, 3)
2951 DO j = bo(1, 2), bo(2, 2)
2952 DO i = bo(1, 1), bo(2, 1)
2953 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2954 deriv_data(i, j, k)*laplace1a(i, j, k)
2960 IF (
ASSOCIATED(deriv_att))
THEN
2964 DO k = bo(1, 3), bo(2, 3)
2965 DO j = bo(1, 2), bo(2, 2)
2966 DO i = bo(1, 1), bo(2, 1)
2967 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2968 deriv_data(i, j, k)*laplace1b(i, j, k)
2975 IF (
ASSOCIATED(deriv_att))
THEN
2979 IF (my_compute_virial)
THEN
2980 CALL virial_drho_drho1(virial_pw, drhob, drho1b, deriv_data, virial_xc)
2984 v_drhob(2)%array(:, :, :) = v_drhob(2)%array(:, :, :) + &
2985 deriv_data(:, :, :)*drb1drb(:, :, :)/max(gradient_cut, norm_drhob(:, :, :))**2
2990 IF (
ASSOCIATED(deriv_att))
THEN
2994 DO k = bo(1, 3), bo(2, 3)
2995 DO j = bo(1, 2), bo(2, 2)
2996 DO i = bo(1, 1), bo(2, 1)
2997 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
2998 deriv_data(i, j, k)*rho1a(i, j, k)
3004 IF (
ASSOCIATED(deriv_att))
THEN
3008 DO k = bo(1, 3), bo(2, 3)
3009 DO j = bo(1, 2), bo(2, 2)
3010 DO i = bo(1, 1), bo(2, 1)
3011 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3012 deriv_data(i, j, k)*rho1b(i, j, k)
3018 IF (
ASSOCIATED(deriv_att))
THEN
3022 DO k = bo(1, 3), bo(2, 3)
3023 DO j = bo(1, 2), bo(2, 2)
3024 DO i = bo(1, 1), bo(2, 1)
3025 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3026 deriv_data(i, j, k)*dr1dr(i, j, k)
3032 IF (
ASSOCIATED(deriv_att))
THEN
3036 DO k = bo(1, 3), bo(2, 3)
3037 DO j = bo(1, 2), bo(2, 2)
3038 DO i = bo(1, 1), bo(2, 1)
3039 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3040 deriv_data(i, j, k)*dra1dra(i, j, k)
3046 IF (
ASSOCIATED(deriv_att))
THEN
3050 DO k = bo(1, 3), bo(2, 3)
3051 DO j = bo(1, 2), bo(2, 2)
3052 DO i = bo(1, 1), bo(2, 1)
3053 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3054 deriv_data(i, j, k)*drb1drb(i, j, k)
3060 IF (
ASSOCIATED(deriv_att))
THEN
3064 DO k = bo(1, 3), bo(2, 3)
3065 DO j = bo(1, 2), bo(2, 2)
3066 DO i = bo(1, 1), bo(2, 1)
3067 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3068 deriv_data(i, j, k)*tau1a(i, j, k)
3074 IF (
ASSOCIATED(deriv_att))
THEN
3078 DO k = bo(1, 3), bo(2, 3)
3079 DO j = bo(1, 2), bo(2, 2)
3080 DO i = bo(1, 1), bo(2, 1)
3081 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3082 deriv_data(i, j, k)*tau1b(i, j, k)
3088 IF (
ASSOCIATED(deriv_att))
THEN
3092 DO k = bo(1, 3), bo(2, 3)
3093 DO j = bo(1, 2), bo(2, 2)
3094 DO i = bo(1, 1), bo(2, 1)
3095 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3096 deriv_data(i, j, k)*laplace1a(i, j, k)
3102 IF (
ASSOCIATED(deriv_att))
THEN
3106 DO k = bo(1, 3), bo(2, 3)
3107 DO j = bo(1, 2), bo(2, 2)
3108 DO i = bo(1, 1), bo(2, 1)
3109 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3110 deriv_data(i, j, k)*laplace1b(i, j, k)
3118 IF (
ASSOCIATED(deriv_att))
THEN
3122 DO k = bo(1, 3), bo(2, 3)
3123 DO j = bo(1, 2), bo(2, 2)
3124 DO i = bo(1, 1), bo(2, 1)
3125 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3126 deriv_data(i, j, k)*rho1a(i, j, k)
3132 IF (
ASSOCIATED(deriv_att))
THEN
3136 DO k = bo(1, 3), bo(2, 3)
3137 DO j = bo(1, 2), bo(2, 2)
3138 DO i = bo(1, 1), bo(2, 1)
3139 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3140 deriv_data(i, j, k)*rho1b(i, j, k)
3146 IF (
ASSOCIATED(deriv_att))
THEN
3150 DO k = bo(1, 3), bo(2, 3)
3151 DO j = bo(1, 2), bo(2, 2)
3152 DO i = bo(1, 1), bo(2, 1)
3153 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3154 deriv_data(i, j, k)*dr1dr(i, j, k)
3160 IF (
ASSOCIATED(deriv_att))
THEN
3164 DO k = bo(1, 3), bo(2, 3)
3165 DO j = bo(1, 2), bo(2, 2)
3166 DO i = bo(1, 1), bo(2, 1)
3167 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3168 deriv_data(i, j, k)*dra1dra(i, j, k)
3174 IF (
ASSOCIATED(deriv_att))
THEN
3178 DO k = bo(1, 3), bo(2, 3)
3179 DO j = bo(1, 2), bo(2, 2)
3180 DO i = bo(1, 1), bo(2, 1)
3181 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3182 deriv_data(i, j, k)*drb1drb(i, j, k)
3188 IF (
ASSOCIATED(deriv_att))
THEN
3192 DO k = bo(1, 3), bo(2, 3)
3193 DO j = bo(1, 2), bo(2, 2)
3194 DO i = bo(1, 1), bo(2, 1)
3195 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3196 deriv_data(i, j, k)*tau1a(i, j, k)
3202 IF (
ASSOCIATED(deriv_att))
THEN
3206 DO k = bo(1, 3), bo(2, 3)
3207 DO j = bo(1, 2), bo(2, 2)
3208 DO i = bo(1, 1), bo(2, 1)
3209 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3210 deriv_data(i, j, k)*tau1b(i, j, k)
3216 IF (
ASSOCIATED(deriv_att))
THEN
3220 DO k = bo(1, 3), bo(2, 3)
3221 DO j = bo(1, 2), bo(2, 2)
3222 DO i = bo(1, 1), bo(2, 1)
3223 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3224 deriv_data(i, j, k)*laplace1a(i, j, k)
3230 IF (
ASSOCIATED(deriv_att))
THEN
3234 DO k = bo(1, 3), bo(2, 3)
3235 DO j = bo(1, 2), bo(2, 2)
3236 DO i = bo(1, 1), bo(2, 1)
3237 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3238 deriv_data(i, j, k)*laplace1b(i, j, k)
3246 IF (
ASSOCIATED(deriv_att))
THEN
3250 DO k = bo(1, 3), bo(2, 3)
3251 DO j = bo(1, 2), bo(2, 2)
3252 DO i = bo(1, 1), bo(2, 1)
3253 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3254 deriv_data(i, j, k)*rho1a(i, j, k)
3260 IF (
ASSOCIATED(deriv_att))
THEN
3264 DO k = bo(1, 3), bo(2, 3)
3265 DO j = bo(1, 2), bo(2, 2)
3266 DO i = bo(1, 1), bo(2, 1)
3267 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3268 deriv_data(i, j, k)*rho1b(i, j, k)
3274 IF (
ASSOCIATED(deriv_att))
THEN
3278 DO k = bo(1, 3), bo(2, 3)
3279 DO j = bo(1, 2), bo(2, 2)
3280 DO i = bo(1, 1), bo(2, 1)
3281 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3282 deriv_data(i, j, k)*dr1dr(i, j, k)
3288 IF (
ASSOCIATED(deriv_att))
THEN
3292 DO k = bo(1, 3), bo(2, 3)
3293 DO j = bo(1, 2), bo(2, 2)
3294 DO i = bo(1, 1), bo(2, 1)
3295 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3296 deriv_data(i, j, k)*dra1dra(i, j, k)
3302 IF (
ASSOCIATED(deriv_att))
THEN
3306 DO k = bo(1, 3), bo(2, 3)
3307 DO j = bo(1, 2), bo(2, 2)
3308 DO i = bo(1, 1), bo(2, 1)
3309 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3310 deriv_data(i, j, k)*drb1drb(i, j, k)
3316 IF (
ASSOCIATED(deriv_att))
THEN
3320 DO k = bo(1, 3), bo(2, 3)
3321 DO j = bo(1, 2), bo(2, 2)
3322 DO i = bo(1, 1), bo(2, 1)
3323 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3324 deriv_data(i, j, k)*tau1a(i, j, k)
3330 IF (
ASSOCIATED(deriv_att))
THEN
3334 DO k = bo(1, 3), bo(2, 3)
3335 DO j = bo(1, 2), bo(2, 2)
3336 DO i = bo(1, 1), bo(2, 1)
3337 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3338 deriv_data(i, j, k)*tau1b(i, j, k)
3344 IF (
ASSOCIATED(deriv_att))
THEN
3348 DO k = bo(1, 3), bo(2, 3)
3349 DO j = bo(1, 2), bo(2, 2)
3350 DO i = bo(1, 1), bo(2, 1)
3351 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3352 deriv_data(i, j, k)*laplace1a(i, j, k)
3358 IF (
ASSOCIATED(deriv_att))
THEN
3362 DO k = bo(1, 3), bo(2, 3)
3363 DO j = bo(1, 2), bo(2, 2)
3364 DO i = bo(1, 1), bo(2, 1)
3365 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3366 deriv_data(i, j, k)*laplace1b(i, j, k)
3373 IF (my_compute_virial)
THEN
3375 IF (
ASSOCIATED(deriv_att))
THEN
3378 virial_pw%array(:, :, :) = -rho1a(:, :, :)
3379 CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
3383 IF (
ASSOCIATED(deriv_att))
THEN
3387 DO k = bo(1, 3), bo(2, 3)
3388 DO j = bo(1, 2), bo(2, 2)
3389 DO i = bo(1, 1), bo(2, 1)
3390 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3391 deriv_data(i, j, k)*rho1a(i, j, k)
3397 IF (
ASSOCIATED(deriv_att))
THEN
3401 DO k = bo(1, 3), bo(2, 3)
3402 DO j = bo(1, 2), bo(2, 2)
3403 DO i = bo(1, 1), bo(2, 1)
3404 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3405 deriv_data(i, j, k)*rho1b(i, j, k)
3411 IF (
ASSOCIATED(deriv_att))
THEN
3415 DO k = bo(1, 3), bo(2, 3)
3416 DO j = bo(1, 2), bo(2, 2)
3417 DO i = bo(1, 1), bo(2, 1)
3418 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3419 deriv_data(i, j, k)*dr1dr(i, j, k)
3425 IF (
ASSOCIATED(deriv_att))
THEN
3429 DO k = bo(1, 3), bo(2, 3)
3430 DO j = bo(1, 2), bo(2, 2)
3431 DO i = bo(1, 1), bo(2, 1)
3432 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3433 deriv_data(i, j, k)*dra1dra(i, j, k)
3439 IF (
ASSOCIATED(deriv_att))
THEN
3443 DO k = bo(1, 3), bo(2, 3)
3444 DO j = bo(1, 2), bo(2, 2)
3445 DO i = bo(1, 1), bo(2, 1)
3446 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3447 deriv_data(i, j, k)*drb1drb(i, j, k)
3453 IF (
ASSOCIATED(deriv_att))
THEN
3457 DO k = bo(1, 3), bo(2, 3)
3458 DO j = bo(1, 2), bo(2, 2)
3459 DO i = bo(1, 1), bo(2, 1)
3460 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3461 deriv_data(i, j, k)*tau1a(i, j, k)
3467 IF (
ASSOCIATED(deriv_att))
THEN
3471 DO k = bo(1, 3), bo(2, 3)
3472 DO j = bo(1, 2), bo(2, 2)
3473 DO i = bo(1, 1), bo(2, 1)
3474 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3475 deriv_data(i, j, k)*tau1b(i, j, k)
3481 IF (
ASSOCIATED(deriv_att))
THEN
3485 DO k = bo(1, 3), bo(2, 3)
3486 DO j = bo(1, 2), bo(2, 2)
3487 DO i = bo(1, 1), bo(2, 1)
3488 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3489 deriv_data(i, j, k)*laplace1a(i, j, k)
3495 IF (
ASSOCIATED(deriv_att))
THEN
3499 DO k = bo(1, 3), bo(2, 3)
3500 DO j = bo(1, 2), bo(2, 2)
3501 DO i = bo(1, 1), bo(2, 1)
3502 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3503 deriv_data(i, j, k)*laplace1b(i, j, k)
3510 IF (my_compute_virial)
THEN
3512 IF (
ASSOCIATED(deriv_att))
THEN
3515 virial_pw%array(:, :, :) = -rho1b(:, :, :)
3516 CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
3525 IF (
ASSOCIATED(deriv_att))
THEN
3529 DO k = bo(1, 3), bo(2, 3)
3530 DO j = bo(1, 2), bo(2, 2)
3531 DO i = bo(1, 1), bo(2, 1)
3532 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3533 deriv_data(i, j, k)*rho1a(i, j, k)
3539 IF (
ASSOCIATED(deriv_att))
THEN
3543 DO k = bo(1, 3), bo(2, 3)
3544 DO j = bo(1, 2), bo(2, 2)
3545 DO i = bo(1, 1), bo(2, 1)
3546 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3547 deriv_data(i, j, k)*dr1dr(i, j, k)
3553 IF (
ASSOCIATED(deriv_att))
THEN
3557 DO k = bo(1, 3), bo(2, 3)
3558 DO j = bo(1, 2), bo(2, 2)
3559 DO i = bo(1, 1), bo(2, 1)
3560 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3561 deriv_data(i, j, k)*dra1dra(i, j, k)
3567 IF (
ASSOCIATED(deriv_att))
THEN
3571 DO k = bo(1, 3), bo(2, 3)
3572 DO j = bo(1, 2), bo(2, 2)
3573 DO i = bo(1, 1), bo(2, 1)
3574 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3575 deriv_data(i, j, k)*tau1a(i, j, k)
3581 IF (
ASSOCIATED(deriv_att))
THEN
3585 DO k = bo(1, 3), bo(2, 3)
3586 DO j = bo(1, 2), bo(2, 2)
3587 DO i = bo(1, 1), bo(2, 1)
3588 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3589 deriv_data(i, j, k)*laplace1a(i, j, k)
3595 IF (
ASSOCIATED(deriv_att))
THEN
3599 DO k = bo(1, 3), bo(2, 3)
3600 DO j = bo(1, 2), bo(2, 2)
3601 DO i = bo(1, 1), bo(2, 1)
3602 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3603 fac*deriv_data(i, j, k)*rho1b(i, j, k)
3609 IF (
ASSOCIATED(deriv_att))
THEN
3613 DO k = bo(1, 3), bo(2, 3)
3614 DO j = bo(1, 2), bo(2, 2)
3615 DO i = bo(1, 1), bo(2, 1)
3616 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3617 fac*deriv_data(i, j, k)*drb1drb(i, j, k)
3623 IF (
ASSOCIATED(deriv_att))
THEN
3627 DO k = bo(1, 3), bo(2, 3)
3628 DO j = bo(1, 2), bo(2, 2)
3629 DO i = bo(1, 1), bo(2, 1)
3630 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3631 fac*deriv_data(i, j, k)*tau1b(i, j, k)
3637 IF (
ASSOCIATED(deriv_att))
THEN
3641 DO k = bo(1, 3), bo(2, 3)
3642 DO j = bo(1, 2), bo(2, 2)
3643 DO i = bo(1, 1), bo(2, 1)
3644 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3645 fac*deriv_data(i, j, k)*laplace1b(i, j, k)
3653 IF (
ASSOCIATED(deriv_att))
THEN
3657 DO k = bo(1, 3), bo(2, 3)
3658 DO j = bo(1, 2), bo(2, 2)
3659 DO i = bo(1, 1), bo(2, 1)
3660 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3661 deriv_data(i, j, k)*rho1a(i, j, k)
3667 IF (
ASSOCIATED(deriv_att))
THEN
3671 DO k = bo(1, 3), bo(2, 3)
3672 DO j = bo(1, 2), bo(2, 2)
3673 DO i = bo(1, 1), bo(2, 1)
3674 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3675 deriv_data(i, j, k)*dr1dr(i, j, k)
3681 IF (
ASSOCIATED(deriv_att))
THEN
3685 DO k = bo(1, 3), bo(2, 3)
3686 DO j = bo(1, 2), bo(2, 2)
3687 DO i = bo(1, 1), bo(2, 1)
3688 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3689 deriv_data(i, j, k)*dra1dra(i, j, k)
3695 IF (
ASSOCIATED(deriv_att))
THEN
3699 DO k = bo(1, 3), bo(2, 3)
3700 DO j = bo(1, 2), bo(2, 2)
3701 DO i = bo(1, 1), bo(2, 1)
3702 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3703 deriv_data(i, j, k)*tau1a(i, j, k)
3709 IF (
ASSOCIATED(deriv_att))
THEN
3713 DO k = bo(1, 3), bo(2, 3)
3714 DO j = bo(1, 2), bo(2, 2)
3715 DO i = bo(1, 1), bo(2, 1)
3716 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3717 deriv_data(i, j, k)*laplace1a(i, j, k)
3723 IF (
ASSOCIATED(deriv_att))
THEN
3727 DO k = bo(1, 3), bo(2, 3)
3728 DO j = bo(1, 2), bo(2, 2)
3729 DO i = bo(1, 1), bo(2, 1)
3730 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3731 fac*deriv_data(i, j, k)*rho1b(i, j, k)
3737 IF (
ASSOCIATED(deriv_att))
THEN
3741 DO k = bo(1, 3), bo(2, 3)
3742 DO j = bo(1, 2), bo(2, 2)
3743 DO i = bo(1, 1), bo(2, 1)
3744 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3745 fac*deriv_data(i, j, k)*drb1drb(i, j, k)
3751 IF (
ASSOCIATED(deriv_att))
THEN
3755 DO k = bo(1, 3), bo(2, 3)
3756 DO j = bo(1, 2), bo(2, 2)
3757 DO i = bo(1, 1), bo(2, 1)
3758 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3759 fac*deriv_data(i, j, k)*tau1b(i, j, k)
3765 IF (
ASSOCIATED(deriv_att))
THEN
3769 DO k = bo(1, 3), bo(2, 3)
3770 DO j = bo(1, 2), bo(2, 2)
3771 DO i = bo(1, 1), bo(2, 1)
3772 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3773 fac*deriv_data(i, j, k)*laplace1b(i, j, k)
3780 IF (
ASSOCIATED(deriv_att))
THEN
3786 v_drho(1)%array(:, :, :) = v_drho(1)%array(:, :, :) + &
3787 deriv_data(:, :, :)*dr1dr(:, :, :)/max(gradient_cut, norm_drho(:, :, :))**2
3792 IF (
ASSOCIATED(deriv_att))
THEN
3796 DO k = bo(1, 3), bo(2, 3)
3797 DO j = bo(1, 2), bo(2, 2)
3798 DO i = bo(1, 1), bo(2, 1)
3799 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3800 deriv_data(i, j, k)*rho1a(i, j, k)
3806 IF (
ASSOCIATED(deriv_att))
THEN
3810 DO k = bo(1, 3), bo(2, 3)
3811 DO j = bo(1, 2), bo(2, 2)
3812 DO i = bo(1, 1), bo(2, 1)
3813 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3814 deriv_data(i, j, k)*dr1dr(i, j, k)
3820 IF (
ASSOCIATED(deriv_att))
THEN
3824 DO k = bo(1, 3), bo(2, 3)
3825 DO j = bo(1, 2), bo(2, 2)
3826 DO i = bo(1, 1), bo(2, 1)
3827 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3828 deriv_data(i, j, k)*dra1dra(i, j, k)
3834 IF (
ASSOCIATED(deriv_att))
THEN
3838 DO k = bo(1, 3), bo(2, 3)
3839 DO j = bo(1, 2), bo(2, 2)
3840 DO i = bo(1, 1), bo(2, 1)
3841 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3842 deriv_data(i, j, k)*tau1a(i, j, k)
3848 IF (
ASSOCIATED(deriv_att))
THEN
3852 DO k = bo(1, 3), bo(2, 3)
3853 DO j = bo(1, 2), bo(2, 2)
3854 DO i = bo(1, 1), bo(2, 1)
3855 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3856 deriv_data(i, j, k)*laplace1a(i, j, k)
3862 IF (
ASSOCIATED(deriv_att))
THEN
3866 DO k = bo(1, 3), bo(2, 3)
3867 DO j = bo(1, 2), bo(2, 2)
3868 DO i = bo(1, 1), bo(2, 1)
3869 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3870 fac*deriv_data(i, j, k)*rho1b(i, j, k)
3876 IF (
ASSOCIATED(deriv_att))
THEN
3880 DO k = bo(1, 3), bo(2, 3)
3881 DO j = bo(1, 2), bo(2, 2)
3882 DO i = bo(1, 1), bo(2, 1)
3883 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3884 fac*deriv_data(i, j, k)*drb1drb(i, j, k)
3890 IF (
ASSOCIATED(deriv_att))
THEN
3894 DO k = bo(1, 3), bo(2, 3)
3895 DO j = bo(1, 2), bo(2, 2)
3896 DO i = bo(1, 1), bo(2, 1)
3897 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3898 fac*deriv_data(i, j, k)*tau1b(i, j, k)
3904 IF (
ASSOCIATED(deriv_att))
THEN
3908 DO k = bo(1, 3), bo(2, 3)
3909 DO j = bo(1, 2), bo(2, 2)
3910 DO i = bo(1, 1), bo(2, 1)
3911 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3912 fac*deriv_data(i, j, k)*laplace1b(i, j, k)
3919 IF (
ASSOCIATED(deriv_att))
THEN
3925 v_drhoa(1)%array(:, :, :) = v_drhoa(1)%array(:, :, :) + &
3926 deriv_data(:, :, :)*dra1dra(:, :, :)/max(gradient_cut, norm_drhoa(:, :, :))**2
3931 IF (
ASSOCIATED(deriv_att))
THEN
3935 DO k = bo(1, 3), bo(2, 3)
3936 DO j = bo(1, 2), bo(2, 2)
3937 DO i = bo(1, 1), bo(2, 1)
3938 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3939 deriv_data(i, j, k)*rho1a(i, j, k)
3945 IF (
ASSOCIATED(deriv_att))
THEN
3949 DO k = bo(1, 3), bo(2, 3)
3950 DO j = bo(1, 2), bo(2, 2)
3951 DO i = bo(1, 1), bo(2, 1)
3952 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3953 deriv_data(i, j, k)*dr1dr(i, j, k)
3959 IF (
ASSOCIATED(deriv_att))
THEN
3963 DO k = bo(1, 3), bo(2, 3)
3964 DO j = bo(1, 2), bo(2, 2)
3965 DO i = bo(1, 1), bo(2, 1)
3966 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3967 deriv_data(i, j, k)*dra1dra(i, j, k)
3973 IF (
ASSOCIATED(deriv_att))
THEN
3977 DO k = bo(1, 3), bo(2, 3)
3978 DO j = bo(1, 2), bo(2, 2)
3979 DO i = bo(1, 1), bo(2, 1)
3980 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3981 deriv_data(i, j, k)*tau1a(i, j, k)
3987 IF (
ASSOCIATED(deriv_att))
THEN
3991 DO k = bo(1, 3), bo(2, 3)
3992 DO j = bo(1, 2), bo(2, 2)
3993 DO i = bo(1, 1), bo(2, 1)
3994 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3995 deriv_data(i, j, k)*laplace1a(i, j, k)
4001 IF (
ASSOCIATED(deriv_att))
THEN
4005 DO k = bo(1, 3), bo(2, 3)
4006 DO j = bo(1, 2), bo(2, 2)
4007 DO i = bo(1, 1), bo(2, 1)
4008 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4009 fac*deriv_data(i, j, k)*rho1b(i, j, k)
4015 IF (
ASSOCIATED(deriv_att))
THEN
4019 DO k = bo(1, 3), bo(2, 3)
4020 DO j = bo(1, 2), bo(2, 2)
4021 DO i = bo(1, 1), bo(2, 1)
4022 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4023 fac*deriv_data(i, j, k)*drb1drb(i, j, k)
4029 IF (
ASSOCIATED(deriv_att))
THEN
4033 DO k = bo(1, 3), bo(2, 3)
4034 DO j = bo(1, 2), bo(2, 2)
4035 DO i = bo(1, 1), bo(2, 1)
4036 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4037 fac*deriv_data(i, j, k)*tau1b(i, j, k)
4043 IF (
ASSOCIATED(deriv_att))
THEN
4047 DO k = bo(1, 3), bo(2, 3)
4048 DO j = bo(1, 2), bo(2, 2)
4049 DO i = bo(1, 1), bo(2, 1)
4050 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4051 fac*deriv_data(i, j, k)*laplace1b(i, j, k)
4059 IF (
ASSOCIATED(deriv_att))
THEN
4063 DO k = bo(1, 3), bo(2, 3)
4064 DO j = bo(1, 2), bo(2, 2)
4065 DO i = bo(1, 1), bo(2, 1)
4066 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4067 deriv_data(i, j, k)*rho1a(i, j, k)
4073 IF (
ASSOCIATED(deriv_att))
THEN
4077 DO k = bo(1, 3), bo(2, 3)
4078 DO j = bo(1, 2), bo(2, 2)
4079 DO i = bo(1, 1), bo(2, 1)
4080 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4081 deriv_data(i, j, k)*dr1dr(i, j, k)
4087 IF (
ASSOCIATED(deriv_att))
THEN
4091 DO k = bo(1, 3), bo(2, 3)
4092 DO j = bo(1, 2), bo(2, 2)
4093 DO i = bo(1, 1), bo(2, 1)
4094 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4095 deriv_data(i, j, k)*dra1dra(i, j, k)
4101 IF (
ASSOCIATED(deriv_att))
THEN
4105 DO k = bo(1, 3), bo(2, 3)
4106 DO j = bo(1, 2), bo(2, 2)
4107 DO i = bo(1, 1), bo(2, 1)
4108 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4109 deriv_data(i, j, k)*tau1a(i, j, k)
4115 IF (
ASSOCIATED(deriv_att))
THEN
4119 DO k = bo(1, 3), bo(2, 3)
4120 DO j = bo(1, 2), bo(2, 2)
4121 DO i = bo(1, 1), bo(2, 1)
4122 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4123 deriv_data(i, j, k)*laplace1a(i, j, k)
4129 IF (
ASSOCIATED(deriv_att))
THEN
4133 DO k = bo(1, 3), bo(2, 3)
4134 DO j = bo(1, 2), bo(2, 2)
4135 DO i = bo(1, 1), bo(2, 1)
4136 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4137 fac*deriv_data(i, j, k)*rho1b(i, j, k)
4143 IF (
ASSOCIATED(deriv_att))
THEN
4147 DO k = bo(1, 3), bo(2, 3)
4148 DO j = bo(1, 2), bo(2, 2)
4149 DO i = bo(1, 1), bo(2, 1)
4150 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4151 fac*deriv_data(i, j, k)*drb1drb(i, j, k)
4157 IF (
ASSOCIATED(deriv_att))
THEN
4161 DO k = bo(1, 3), bo(2, 3)
4162 DO j = bo(1, 2), bo(2, 2)
4163 DO i = bo(1, 1), bo(2, 1)
4164 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4165 fac*deriv_data(i, j, k)*tau1b(i, j, k)
4171 IF (
ASSOCIATED(deriv_att))
THEN
4175 DO k = bo(1, 3), bo(2, 3)
4176 DO j = bo(1, 2), bo(2, 2)
4177 DO i = bo(1, 1), bo(2, 1)
4178 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4179 fac*deriv_data(i, j, k)*laplace1b(i, j, k)
4190 IF (gradient_f)
THEN
4191 IF (.NOT. do_spinflip)
THEN
4193 IF (my_compute_virial)
THEN
4194 CALL virial_drho_drho(virial_pw, drhoa, v_drhoa(1), virial_xc)
4195 CALL virial_drho_drho(virial_pw, drhob, v_drhob(2), virial_xc)
4198 virial_pw%array(:, :, :) = drho(idir)%array(:, :, :)*(v_drho(1)%array(:, :, :) + v_drho(2)%array(:, :, :))
4202 drho(jdir)%array(:, :, :))
4203 virial_xc(jdir, idir) = virial_xc(jdir, idir) + tmp
4204 virial_xc(idir, jdir) = virial_xc(jdir, idir)
4214 DO ir = bo(1, 2), bo(2, 2)
4215 DO ia = bo(1, 1), bo(2, 1)
4217 DO ispin = 1, nspins
4218 vxg(idir, ia, ir, ispin) = &
4219 -(v_drhoa(ispin)%array(ia, ir, 1)*drhoa(idir)%array(ia, ir, 1) + &
4220 v_drhob(ispin)%array(ia, ir, 1)*drhob(idir)%array(ia, ir, 1) + &
4221 v_drho(ispin)%array(ia, ir, 1)*drho(idir)%array(ia, ir, 1))
4223 IF (
ASSOCIATED(e_drhoa))
THEN
4224 vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
4225 e_drhoa(ia, ir, 1)*drho1a(idir)%array(ia, ir, 1)
4227 IF (nspins /= 1 .AND.
ASSOCIATED(e_drhob))
THEN
4228 vxg(idir, ia, ir, 2) = vxg(idir, ia, ir, 2) + &
4229 e_drhob(ia, ir, 1)*drho1b(idir)%array(ia, ir, 1)
4231 IF (
ASSOCIATED(e_drho))
THEN
4232 IF (nspins /= 1)
THEN
4233 vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
4234 e_drho(ia, ir, 1)*drho1(idir)%array(ia, ir, 1)
4235 vxg(idir, ia, ir, 2) = vxg(idir, ia, ir, 2) + &
4236 e_drho(ia, ir, 1)*drho1(idir)%array(ia, ir, 1)
4238 vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
4239 e_drho(ia, ir, 1)*(drho1a(idir)%array(ia, ir, 1) + &
4240 fac*drho1b(idir)%array(ia, ir, 1))
4252 DO ispin = 1, nspins
4255 v_drho_r(idir, ispin)%array(:, :, :) = &
4256 v_drhoa(ispin)%array(:, :, :)*drhoa(idir)%array(:, :, :) + &
4257 v_drhob(ispin)%array(:, :, :)*drhob(idir)%array(:, :, :) + &
4258 v_drho(ispin)%array(:, :, :)*drho(idir)%array(:, :, :)
4261 IF (
ASSOCIATED(e_drhoa))
THEN
4264 v_drho_r(idir, 1)%array(:, :, :) = v_drho_r(idir, 1)%array(:, :, :) - &
4265 e_drhoa(:, :, :)*drho1a(idir)%array(:, :, :)
4268 IF (nspins /= 1 .AND.
ASSOCIATED(e_drhob))
THEN
4271 v_drho_r(idir, 2)%array(:, :, :) = v_drho_r(idir, 2)%array(:, :, :) - &
4272 e_drhob(:, :, :)*drho1b(idir)%array(:, :, :)
4275 IF (
ASSOCIATED(e_drho))
THEN
4278 DO k = bo(1, 3), bo(2, 3)
4279 DO j = bo(1, 2), bo(2, 2)
4280 DO i = bo(1, 1), bo(2, 1)
4281 IF (nspins /= 1)
THEN
4282 v_drho_r(idir, 1)%array(i, j, k) = v_drho_r(idir, 1)%array(i, j, k) - &
4283 e_drho(i, j, k)*drho1(idir)%array(i, j, k)
4284 v_drho_r(idir, 2)%array(i, j, k) = v_drho_r(idir, 2)%array(i, j, k) - &
4285 e_drho(i, j, k)*drho1(idir)%array(i, j, k)
4287 v_drho_r(idir, 1)%array(i, j, k) = v_drho_r(idir, 1)%array(i, j, k) - &
4288 e_drho(i, j, k)*(drho1a(idir)%array(i, j, k) + &
4289 fac*drho1b(idir)%array(i, j, k))
4299 DO ispin = 1, nspins
4300 CALL xc_pw_divergence(xc_deriv_method_id, v_drho_r(:, ispin), tmp_g, vxc_g, v_xc(ispin))
4308 DEALLOCATE (drho(idir)%array)
4309 DEALLOCATE (drho1(idir)%array)
4312 DO ispin = 1, nspins
4313 CALL deallocate_pw(v_drhoa(ispin), pw_pool)
4314 CALL deallocate_pw(v_drhob(ispin), pw_pool)
4317 DEALLOCATE (v_drhoa, v_drhob)
4321 IF (laplace_f .AND. my_compute_virial)
THEN
4322 virial_pw%array(:, :, :) = -rhoa(:, :, :)
4323 CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(1)%array)
4324 virial_pw%array(:, :, :) = -rhob(:, :, :)
4325 CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(2)%array)
4336 IF (gradient_f)
THEN
4339 CALL prepare_dr1dr(dr1dr, drho, drho1)
4345 ALLOCATE (v_laplace(nspins))
4346 DO ispin = 1, nspins
4347 CALL allocate_pw(v_laplace(ispin), pw_pool, bo)
4358 IF (
ASSOCIATED(deriv_att))
THEN
4362 DO k = bo(1, 3), bo(2, 3)
4363 DO j = bo(1, 2), bo(2, 2)
4364 DO i = bo(1, 1), bo(2, 1)
4365 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4366 deriv_data(i, j, k)*rho1(i, j, k)
4372 IF (
ASSOCIATED(deriv_att))
THEN
4376 DO k = bo(1, 3), bo(2, 3)
4377 DO j = bo(1, 2), bo(2, 2)
4378 DO i = bo(1, 1), bo(2, 1)
4379 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4380 deriv_data(i, j, k)*dr1dr(i, j, k)
4386 IF (
ASSOCIATED(deriv_att))
THEN
4390 DO k = bo(1, 3), bo(2, 3)
4391 DO j = bo(1, 2), bo(2, 2)
4392 DO i = bo(1, 1), bo(2, 1)
4393 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4394 deriv_data(i, j, k)*tau1(i, j, k)
4400 IF (
ASSOCIATED(deriv_att))
THEN
4404 DO k = bo(1, 3), bo(2, 3)
4405 DO j = bo(1, 2), bo(2, 2)
4406 DO i = bo(1, 1), bo(2, 1)
4407 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4408 deriv_data(i, j, k)*laplace1(i, j, k)
4416 IF (
ASSOCIATED(deriv_att))
THEN
4420 DO k = bo(1, 3), bo(2, 3)
4421 DO j = bo(1, 2), bo(2, 2)
4422 DO i = bo(1, 1), bo(2, 1)
4423 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
4424 deriv_data(i, j, k)*rho1(i, j, k)
4430 IF (
ASSOCIATED(deriv_att))
THEN
4434 DO k = bo(1, 3), bo(2, 3)
4435 DO j = bo(1, 2), bo(2, 2)
4436 DO i = bo(1, 1), bo(2, 1)
4437 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
4438 deriv_data(i, j, k)*dr1dr(i, j, k)
4444 IF (
ASSOCIATED(deriv_att))
THEN
4448 DO k = bo(1, 3), bo(2, 3)
4449 DO j = bo(1, 2), bo(2, 2)
4450 DO i = bo(1, 1), bo(2, 1)
4451 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
4452 deriv_data(i, j, k)*tau1(i, j, k)
4458 IF (
ASSOCIATED(deriv_att))
THEN
4462 DO k = bo(1, 3), bo(2, 3)
4463 DO j = bo(1, 2), bo(2, 2)
4464 DO i = bo(1, 1), bo(2, 1)
4465 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
4466 deriv_data(i, j, k)*laplace1(i, j, k)
4473 IF (
ASSOCIATED(deriv_att))
THEN
4477 IF (my_compute_virial)
THEN
4478 CALL virial_drho_drho1(virial_pw, drho, drho1, deriv_data, virial_xc)
4482 v_drho(1)%array(:, :, :) = v_drho(1)%array(:, :, :) + &
4483 deriv_data(:, :, :)*dr1dr(:, :, :)/max(gradient_cut, norm_drho(:, :, :))**2
4488 IF (
ASSOCIATED(deriv_att))
THEN
4492 DO k = bo(1, 3), bo(2, 3)
4493 DO j = bo(1, 2), bo(2, 2)
4494 DO i = bo(1, 1), bo(2, 1)
4495 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4496 deriv_data(i, j, k)*rho1(i, j, k)
4502 IF (
ASSOCIATED(deriv_att))
THEN
4506 DO k = bo(1, 3), bo(2, 3)
4507 DO j = bo(1, 2), bo(2, 2)
4508 DO i = bo(1, 1), bo(2, 1)
4509 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4510 deriv_data(i, j, k)*dr1dr(i, j, k)
4516 IF (
ASSOCIATED(deriv_att))
THEN
4520 DO k = bo(1, 3), bo(2, 3)
4521 DO j = bo(1, 2), bo(2, 2)
4522 DO i = bo(1, 1), bo(2, 1)
4523 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4524 deriv_data(i, j, k)*tau1(i, j, k)
4530 IF (
ASSOCIATED(deriv_att))
THEN
4534 DO k = bo(1, 3), bo(2, 3)
4535 DO j = bo(1, 2), bo(2, 2)
4536 DO i = bo(1, 1), bo(2, 1)
4537 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4538 deriv_data(i, j, k)*laplace1(i, j, k)
4546 IF (
ASSOCIATED(deriv_att))
THEN
4550 DO k = bo(1, 3), bo(2, 3)
4551 DO j = bo(1, 2), bo(2, 2)
4552 DO i = bo(1, 1), bo(2, 1)
4553 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
4554 deriv_data(i, j, k)*rho1(i, j, k)
4560 IF (
ASSOCIATED(deriv_att))
THEN
4564 DO k = bo(1, 3), bo(2, 3)
4565 DO j = bo(1, 2), bo(2, 2)
4566 DO i = bo(1, 1), bo(2, 1)
4567 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
4568 deriv_data(i, j, k)*dr1dr(i, j, k)
4574 IF (
ASSOCIATED(deriv_att))
THEN
4578 DO k = bo(1, 3), bo(2, 3)
4579 DO j = bo(1, 2), bo(2, 2)
4580 DO i = bo(1, 1), bo(2, 1)
4581 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
4582 deriv_data(i, j, k)*tau1(i, j, k)
4588 IF (
ASSOCIATED(deriv_att))
THEN
4592 DO k = bo(1, 3), bo(2, 3)
4593 DO j = bo(1, 2), bo(2, 2)
4594 DO i = bo(1, 1), bo(2, 1)
4595 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
4596 deriv_data(i, j, k)*laplace1(i, j, k)
4603 IF (my_compute_virial)
THEN
4605 IF (
ASSOCIATED(deriv_att))
THEN
4608 virial_pw%array(:, :, :) = -rho1(:, :, :)
4609 CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
4614 IF (gradient_f)
THEN
4616 IF (my_compute_virial)
THEN
4617 CALL virial_drho_drho(virial_pw, drho, v_drho(1), virial_xc)
4627 DO ia = bo(1, 1), bo(2, 1)
4628 DO ir = bo(1, 2), bo(2, 2)
4629 vxg(idir, ia, ir, 1) = -drho(idir)%array(ia, ir, 1)*v_drho(1)%array(ia, ir, 1)
4630 IF (
ASSOCIATED(e_drho))
THEN
4631 vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + factor2*drho1(idir)%array(ia, ir, 1)*e_drho(ia, ir, 1)
4642 v_drho_r(idir, 1)%array(:, :, :) = drho(idir)%array(:, :, :)*v_drho(1)%array(:, :, :) - &
4643 drho1(idir)%array(:, :, :)*e_drho(:, :, :)
4647 CALL xc_pw_divergence(xc_deriv_method_id, v_drho_r(:, 1), tmp_g, vxc_g, v_xc(1))
4652 IF (laplace_f .AND. my_compute_virial)
THEN
4653 virial_pw%array(:, :, :) = -rho(:, :, :)
4654 CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(1)%array)
4660 DO ispin = 1, nspins
4661 CALL xc_pw_laplace(v_laplace(ispin), pw_pool, xc_deriv_method_id)
4662 CALL pw_axpy(v_laplace(ispin), v_xc(ispin))
4666 IF (gradient_f)
THEN
4668 DO ispin = 1, nspins
4669 CALL deallocate_pw(v_drho(ispin), pw_pool)
4671 CALL deallocate_pw(v_drho_r(idir, ispin), pw_pool)
4674 DEALLOCATE (v_drho, v_drho_r)
4679 DO ispin = 1, nspins
4680 CALL deallocate_pw(v_laplace(ispin), pw_pool)
4682 DEALLOCATE (v_laplace)
4685 IF (
ASSOCIATED(tmp_g%pw_grid) .AND.
ASSOCIATED(pw_pool))
THEN
4686 CALL pw_pool%give_back_pw(tmp_g)
4689 IF (
ASSOCIATED(vxc_g%pw_grid) .AND.
ASSOCIATED(pw_pool))
THEN
4690 CALL pw_pool%give_back_pw(vxc_g)
4693 IF (my_compute_virial .AND. (gradient_f .OR. laplace_f))
THEN
4694 CALL deallocate_pw(virial_pw, pw_pool)
4697 CALL timestop(handle)
4733 LOGICAL,
INTENT(in),
OPTIONAL :: spinflip
4735 CHARACTER(len=*),
PARAMETER :: routinen =
'xc_calc_3rd_deriv_analytical'
4737 INTEGER :: handle, i, idir, ispin, j, &
4738 k, nspins, xc_deriv_method_id
4739 INTEGER,
DIMENSION(2, 3) :: bo
4740 LOGICAL :: lsd, do_spinflip, alda0, &
4741 rho_f, gradient_f, tau_f, laplace_f
4742 REAL(kind=
dp) :: s, s_thresh, s_thresh2, gradient_cut
4743 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: dr1dr, dra1dra, drb1drb
4744 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: deriv_data, deriv_data2, e_drhoa, e_drhob, &
4745 e_drho, norm_drho, norm_drhoa, &
4746 norm_drhob, rho1a, rho1b, &
4748 TYPE(
cp_3d_r_cp_type),
DIMENSION(3) :: drho, drho1, drho1a, drho1b, drhoa, drhob
4749 TYPE(
pw_r3d_rs_type),
DIMENSION(:),
ALLOCATABLE :: v_drhoa, v_drhob, v_drho
4754 CALL timeset(routinen, handle)
4756 NULLIFY (e_drhoa, e_drhob, e_drho)
4758 cpassert(
ASSOCIATED(v_xc))
4759 cpassert(
ASSOCIATED(xc_section))
4763 i_val=xc_deriv_method_id)
4766 lsd =
ASSOCIATED(rho_set%rhoa)
4768 do_spinflip = .false.
4769 IF (
PRESENT(spinflip)) do_spinflip = spinflip
4771 bo = rho_set%local_bounds
4773 CALL check_for_derivatives(deriv_set, lsd, rho_f, gradient_f, tau_f, laplace_f)
4783 DO ispin = 1, nspins
4785 v_xc(ispin)%array = 0.0_dp
4789 IF (gradient_f)
THEN
4790 ALLOCATE (v_drho_r(3, nspins), v_drho(nspins))
4791 DO ispin = 1, nspins
4793 CALL allocate_pw(v_drho_r(idir, ispin), pw_pool, bo)
4795 CALL allocate_pw(v_drho(ispin), pw_pool, bo)
4799 IF (
ASSOCIATED(pw_pool))
THEN
4800 CALL pw_pool%create_pw(tmp_g)
4801 CALL pw_pool%create_pw(vxc_g)
4804 cpabort(
"XC_DERIV method is not implemented in GAPW")
4812 cpassert(
ASSOCIATED(v_xc_tau))
4813 DO ispin = 1, nspins
4814 v_xc_tau(ispin)%array = 0.0_dp
4824 IF (do_spinflip)
THEN
4831 IF (gradient_f)
THEN
4833 norm_drho=norm_drho, norm_drhoa=norm_drhoa, norm_drhob=norm_drhob)
4834 IF (do_spinflip)
THEN
4836 CALL calc_drho_from_a(drho1, drho1a)
4839 CALL calc_drho_from_ab(drho1, drho1a, drho1b)
4842 CALL calc_drho_from_ab(drho, drhoa, drhob)
4844 CALL prepare_dr1dr(dra1dra, drhoa, drho1a)
4845 IF (do_spinflip)
THEN
4846 CALL prepare_dr1dr(drb1drb, drhob, drho1a)
4847 CALL prepare_dr1dr(dr1dr, drho, drho1a)
4848 ELSE IF (nspins /= 1)
THEN
4849 CALL prepare_dr1dr(drb1drb, drhob, drho1b)
4850 CALL prepare_dr1dr(dr1dr, drho, drho1)
4852 cpabort(
"Exchange-correlation's third derivative for closed-shell not yet implemented")
4856 ALLOCATE (v_drhoa(nspins), v_drhob(nspins))
4857 DO ispin = 1, nspins
4858 CALL allocate_pw(v_drhoa(ispin), pw_pool, bo)
4859 CALL allocate_pw(v_drhob(ispin), pw_pool, bo)
4865 cpabort(
"Exchange-correlation's laplace analytic third derivative not implemented")
4869 cpabort(
"Exchange-correlation's mGGA analytic third derivative not implemented")
4872 IF (nspins /= 1)
THEN
4874 IF (.NOT. do_spinflip)
THEN
4876 cpabort(
"Exchange-correlation's analytic third derivative not implemented")
4888 IF (
ASSOCIATED(deriv_att))
THEN
4891 IF (
ASSOCIATED(deriv_att))
THEN
4895 DO k = bo(1, 3), bo(2, 3)
4896 DO j = bo(1, 2), bo(2, 2)
4897 DO i = bo(1, 1), bo(2, 1)
4898 s = rhoa(i, j, k) - rhob(i, j, k)
4899 s = -sign(max(s**2, s_thresh2), s)
4900 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4901 (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)**2/s
4913 IF (.NOT. alda0)
THEN
4915 IF (
ASSOCIATED(deriv_att))
THEN
4918 IF (
ASSOCIATED(deriv_att))
THEN
4922 DO k = bo(1, 3), bo(2, 3)
4923 DO j = bo(1, 2), bo(2, 2)
4924 DO i = bo(1, 1), bo(2, 1)
4925 s = rhoa(i, j, k) - rhob(i, j, k)
4926 s = -sign(max(s**2, s_thresh2), s)
4927 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4928 (deriv_data(i, j, k)*dra1dra(i, j, k) - deriv_data2(i, j, k)*drb1drb(i, j, k)) &
4940 DO k = bo(1, 3), bo(2, 3)
4941 DO j = bo(1, 2), bo(2, 2)
4942 DO i = bo(1, 1), bo(2, 1)
4943 v_xc(2)%array(i, j, k) = -v_xc(1)%array(i, j, k)
4956 IF (
ASSOCIATED(deriv_att))
THEN
4959 IF (
ASSOCIATED(deriv_att))
THEN
4963 DO k = bo(1, 3), bo(2, 3)
4964 DO j = bo(1, 2), bo(2, 2)
4965 DO i = bo(1, 1), bo(2, 1)
4966 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
4967 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4968 rho1a(i, j, k)**2*(deriv_data(i, j, k) - deriv_data2(i, j, k))/s
4980 IF (
ASSOCIATED(deriv_att))
THEN
4983 IF (
ASSOCIATED(deriv_att))
THEN
4987 DO k = bo(1, 3), bo(2, 3)
4988 DO j = bo(1, 2), bo(2, 2)
4989 DO i = bo(1, 1), bo(2, 1)
4990 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
4991 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
4992 rho1a(i, j, k)**2*(deriv_data(i, j, k) - deriv_data2(i, j, k))/s
5000 IF (.NOT. alda0)
THEN
5005 IF (
ASSOCIATED(deriv_att))
THEN
5008 IF (
ASSOCIATED(deriv_att))
THEN
5012 DO k = bo(1, 3), bo(2, 3)
5013 DO j = bo(1, 2), bo(2, 2)
5014 DO i = bo(1, 1), bo(2, 1)
5015 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
5016 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
5017 (deriv_data(i, j, k)*dra1dra(i, j, k) - deriv_data2(i, j, k)*drb1drb(i, j, k))* &
5030 IF (
ASSOCIATED(deriv_att))
THEN
5033 IF (
ASSOCIATED(deriv_att))
THEN
5037 DO k = bo(1, 3), bo(2, 3)
5038 DO j = bo(1, 2), bo(2, 2)
5039 DO i = bo(1, 1), bo(2, 1)
5040 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
5041 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
5042 (deriv_data(i, j, k)*dra1dra(i, j, k) - deriv_data2(i, j, k)*drb1drb(i, j, k))* &
5058 IF (
ASSOCIATED(deriv_att))
THEN
5061 IF (
ASSOCIATED(deriv_att))
THEN
5065 DO k = bo(1, 3), bo(2, 3)
5066 DO j = bo(1, 2), bo(2, 2)
5067 DO i = bo(1, 1), bo(2, 1)
5068 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
5069 (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)
5081 IF (
ASSOCIATED(deriv_att))
THEN
5084 IF (
ASSOCIATED(deriv_att))
THEN
5088 DO k = bo(1, 3), bo(2, 3)
5089 DO j = bo(1, 2), bo(2, 2)
5090 DO i = bo(1, 1), bo(2, 1)
5091 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) + &
5092 (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)
5108 IF (
ASSOCIATED(deriv_att))
THEN
5111 IF (
ASSOCIATED(deriv_att))
THEN
5115 DO k = bo(1, 3), bo(2, 3)
5116 DO j = bo(1, 2), bo(2, 2)
5117 DO i = bo(1, 1), bo(2, 1)
5118 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
5119 (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)
5120 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
5121 (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)
5130 IF (
ASSOCIATED(deriv_att))
THEN
5136 IF (
ASSOCIATED(deriv_att))
THEN
5140 DO k = bo(1, 3), bo(2, 3)
5141 DO j = bo(1, 2), bo(2, 2)
5142 DO i = bo(1, 1), bo(2, 1)
5143 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
5144 deriv_data(i, j, k)*dra1dra(i, j, k) + deriv_data2(i, j, k)*drb1drb(i, j, k)
5154 IF (
ASSOCIATED(deriv_att))
THEN
5158 DO k = bo(1, 3), bo(2, 3)
5159 DO j = bo(1, 2), bo(2, 2)
5160 DO i = bo(1, 1), bo(2, 1)
5161 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
5162 deriv_data(i, j, k)*dra1dra(i, j, k) + deriv_data2(i, j, k)*drb1drb(i, j, k)
5177 IF (
ASSOCIATED(deriv_att))
THEN
5180 IF (
ASSOCIATED(deriv_att))
THEN
5184 DO k = bo(1, 3), bo(2, 3)
5185 DO j = bo(1, 2), bo(2, 2)
5186 DO i = bo(1, 1), bo(2, 1)
5187 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
5188 deriv_data(i, j, k)*dra1dra(i, j, k) + &
5189 deriv_data2(i, j, k)*drb1drb(i, j, k)
5190 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
5191 deriv_data(i, j, k)*dra1dra(i, j, k) + &
5192 deriv_data2(i, j, k)*drb1drb(i, j, k)
5207 IF (
ASSOCIATED(deriv_att))
THEN
5212 v_drhoa(1)%array(:, :, :) = v_drhoa(1)%array(:, :, :) + &
5213 deriv_data(:, :, :)*dra1dra(:, :, :)/max(gradient_cut, norm_drhoa(:, :, :))**2
5221 IF (
ASSOCIATED(deriv_att))
THEN
5226 v_drhob(2)%array(:, :, :) = v_drhob(2)%array(:, :, :) - &
5227 deriv_data(:, :, :)*drb1drb(:, :, :)/max(gradient_cut, norm_drhob(:, :, :))**2
5236 cpabort(
"Exchange-correlation's analytic third derivative not implemented")
5240 IF (gradient_f)
THEN
5241 IF (.NOT. alda0)
THEN
5255 IF (do_spinflip)
THEN
5258 DO k = bo(1, 3), bo(2, 3)
5259 DO j = bo(1, 2), bo(2, 2)
5260 DO i = bo(1, 1), bo(2, 1)
5261 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
5263 v_drho_r(idir, ispin)%array(i, j, k) = v_drho_r(idir, ispin)%array(i, j, k) + &
5264 v_drhoa(ispin)%array(i, j, k)*drhoa(idir)%array(i, j, k)*rho1a(i, j, k)/s + &
5265 v_drhob(ispin)%array(i, j, k)*drhob(idir)%array(i, j, k)*rho1a(i, j, k)/s + &
5266 v_drho(ispin)%array(i, j, k)*drho(idir)%array(i, j, k)*rho1a(i, j, k)/s
5278 IF (
ASSOCIATED(e_drhoa))
THEN
5281 DO k = bo(1, 3), bo(2, 3)
5282 DO j = bo(1, 2), bo(2, 2)
5283 DO i = bo(1, 1), bo(2, 1)
5284 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
5285 v_drho_r(idir, 1)%array(i, j, k) = v_drho_r(idir, 1)%array(i, j, k) - &
5286 e_drhoa(i, j, k)*drho1a(idir)%array(i, j, k)*rho1a(i, j, k)/s
5296 IF (
ASSOCIATED(e_drhob))
THEN
5299 DO k = bo(1, 3), bo(2, 3)
5300 DO j = bo(1, 2), bo(2, 2)
5301 DO i = bo(1, 1), bo(2, 1)
5302 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
5303 v_drho_r(idir, 2)%array(i, j, k) = v_drho_r(idir, 2)%array(i, j, k) + &
5304 e_drhob(i, j, k)*drho1a(idir)%array(i, j, k)*rho1a(i, j, k)/s
5313 DO ispin = 1, nspins
5314 CALL xc_pw_divergence(xc_deriv_method_id, v_drho_r(:, ispin), tmp_g, vxc_g, v_xc(ispin))
5319 DEALLOCATE (drho(idir)%array)
5320 DEALLOCATE (drho1(idir)%array)
5323 DO ispin = 1, nspins
5324 CALL deallocate_pw(v_drhoa(ispin), pw_pool)
5325 CALL deallocate_pw(v_drhob(ispin), pw_pool)
5328 DEALLOCATE (v_drhoa, v_drhob)
5337 cpabort(
"Exchange-correlation's analytic third derivative not implemented")
5341 IF (gradient_f)
THEN
5343 DO ispin = 1, nspins
5344 CALL deallocate_pw(v_drho(ispin), pw_pool)
5346 CALL deallocate_pw(v_drho_r(idir, ispin), pw_pool)
5349 DEALLOCATE (v_drho, v_drho_r)
5353 IF (
ASSOCIATED(tmp_g%pw_grid) .AND.
ASSOCIATED(pw_pool))
THEN
5354 CALL pw_pool%give_back_pw(tmp_g)
5357 IF (
ASSOCIATED(vxc_g%pw_grid) .AND.
ASSOCIATED(pw_pool))
THEN
5358 CALL pw_pool%give_back_pw(vxc_g)
5361 CALL timestop(handle)