1059 pw_pool, weights, xc_section, &
1060 do_triplet, calc_virial, virial_xc, deriv_set)
1062 TYPE(
pw_r3d_rs_type),
DIMENSION(:),
INTENT(IN),
POINTER :: v_xc, v_tau
1064 TYPE(
pw_r3d_rs_type),
DIMENSION(:),
INTENT(IN),
POINTER :: rho1_r, tau1_r
1065 TYPE(
pw_c1d_gs_type),
DIMENSION(:),
INTENT(IN),
POINTER :: rho1_g
1069 LOGICAL,
INTENT(IN) :: do_triplet
1070 LOGICAL,
INTENT(IN),
OPTIONAL :: calc_virial
1071 REAL(kind=
dp),
DIMENSION(3, 3),
INTENT(INOUT), &
1072 OPTIONAL :: virial_xc
1075 CHARACTER(len=*),
PARAMETER :: routinen =
'xc_calc_2nd_deriv_numerical'
1076 REAL(kind=
dp),
DIMENSION(-4:4, 4),
PARAMETER :: &
1077 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, &
1078 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, &
1079 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, &
1080 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])
1082 INTEGER :: handle, idir, ispin, nspins, istep, nsteps
1083 INTEGER,
DIMENSION(2, 3) :: bo
1084 LOGICAL :: gradient_f, lsd, my_calc_virial, tau_f, laplace_f, rho_f
1085 REAL(kind=
dp) :: exc, gradient_cut, h, rweight, step, rho_cutoff
1086 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: dr1dr, dra1dra, drb1drb
1087 REAL(kind=
dp),
DIMENSION(3, 3) :: virial_dummy
1088 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: norm_drho, norm_drho2, norm_drho2a, &
1089 norm_drho2b, norm_drhoa, norm_drhob, &
1090 rho, rho1, rho1a, rho1b, rhoa, rhob, &
1091 tau_a, tau_b, tau, tau1, tau1a, tau1b, laplace, laplace1, &
1092 laplacea, laplaceb, laplace1a, laplace1b, &
1093 laplace2, laplace2a, laplace2b, deriv_data
1094 TYPE(
cp_3d_r_cp_type),
DIMENSION(3) :: drho, drho1, drho1a, drho1b, drhoa, drhob
1099 TYPE(
pw_r3d_rs_type) :: virial_pw, v_laplace, v_laplacea, v_laplaceb
1105 CALL timeset(routinen, handle)
1107 my_calc_virial = .false.
1108 IF (
PRESENT(calc_virial) .AND.
PRESENT(virial_xc)) my_calc_virial = calc_virial
1112 NULLIFY (tau, tau_r, tau_a, tau_b)
1116 IF (nsteps < lbound(rweights, 2) .OR. nspins > ubound(rweights, 2))
THEN
1117 cpabort(
"The number of steps must be a value from 1 to 4.")
1120 IF (nspins == 2)
THEN
1121 NULLIFY (vxc_rho, rho_g, vxc_tau)
1123 DO ispin = 1, nspins
1124 CALL pw_pool%create_pw(rho_r(ispin))
1126 IF (
ASSOCIATED(tau1_r) .AND.
ASSOCIATED(v_tau))
THEN
1128 DO ispin = 1, nspins
1129 CALL pw_pool%create_pw(tau_r(ispin))
1132 CALL xc_rho_set_get(rho_set, can_return_null=.true., rhoa=rhoa, rhob=rhob, tau_a=tau_a, tau_b=tau_b)
1133 DO istep = -nsteps, nsteps
1134 IF (istep == 0) cycle
1135 rweight = rweights(istep, nsteps)/h
1136 step = real(istep,
dp)*h
1137 CALL calc_resp_potential_numer_ab(rho_r, rho_g, rho1_r, rhoa, rhob, vxc_rho, &
1138 tau_r, tau1_r, tau_a, tau_b, vxc_tau, xc_section, &
1139 weights, pw_pool, step)
1140 DO ispin = 1, nspins
1141 CALL pw_axpy(vxc_rho(ispin), v_xc(ispin), rweight)
1142 IF (
ASSOCIATED(vxc_tau) .AND.
ASSOCIATED(v_tau))
THEN
1143 CALL pw_axpy(vxc_tau(ispin), v_tau(ispin), rweight)
1146 DO ispin = 1, nspins
1147 CALL vxc_rho(ispin)%release()
1149 DEALLOCATE (vxc_rho)
1150 IF (
ASSOCIATED(vxc_tau))
THEN
1151 DO ispin = 1, nspins
1152 CALL vxc_tau(ispin)%release()
1154 DEALLOCATE (vxc_tau)
1157 ELSE IF (nspins == 1 .AND. do_triplet)
THEN
1158 NULLIFY (vxc_rho, vxc_tau, rho_g)
1161 CALL pw_pool%create_pw(rho_r(ispin))
1163 IF (
ASSOCIATED(tau1_r) .AND.
ASSOCIATED(v_tau))
THEN
1165 DO ispin = 1, nspins
1166 CALL pw_pool%create_pw(tau_r(ispin))
1169 CALL xc_rho_set_get(rho_set, can_return_null=.true., rhoa=rhoa, rhob=rhob, tau_a=tau_a, tau_b=tau_b)
1170 DO istep = -nsteps, nsteps
1171 IF (istep == 0) cycle
1172 rweight = rweights(istep, nsteps)/h
1173 step = real(istep,
dp)*h
1177 rho_r(1)%array(:, :, :) = rhoa(:, :, :) + step*rho1_r(1)%array(:, :, :)
1180 rho_r(2)%array(:, :, :) = rhob(:, :, :)
1182 IF (
ASSOCIATED(tau1_r))
THEN
1184 tau_r(1)%array(:, :, :) = tau_a(:, :, :) + step*tau1_r(1)%array(:, :, :)
1187 tau_r(2)%array(:, :, :) = tau_b(:, :, :)
1191 CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
1192 weights, pw_pool, .false., virial_dummy)
1193 CALL pw_axpy(vxc_rho(1), v_xc(1), rweight)
1194 IF (
ASSOCIATED(vxc_tau) .AND.
ASSOCIATED(v_tau))
THEN
1195 CALL pw_axpy(vxc_tau(1), v_tau(1), rweight)
1198 CALL vxc_rho(ispin)%release()
1200 DEALLOCATE (vxc_rho)
1201 IF (
ASSOCIATED(vxc_tau))
THEN
1203 CALL vxc_tau(ispin)%release()
1205 DEALLOCATE (vxc_tau)
1210 rho_r(1)%array(:, :, :) = rhoa(:, :, :)
1213 rho_r(2)%array(:, :, :) = rhob(:, :, :) + step*rho1_r(1)%array(:, :, :)
1215 IF (
ASSOCIATED(tau1_r))
THEN
1217 tau_r(1)%array(:, :, :) = tau_a(:, :, :)
1220 tau_r(2)%array(:, :, :) = tau_b(:, :, :) + step*tau1_r(1)%array(:, :, :)
1224 CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
1225 weights, pw_pool, .false., virial_dummy)
1226 CALL pw_axpy(vxc_rho(1), v_xc(1), rweight)
1227 IF (
ASSOCIATED(vxc_tau) .AND.
ASSOCIATED(v_tau))
THEN
1228 CALL pw_axpy(vxc_tau(1), v_tau(1), rweight)
1231 CALL vxc_rho(ispin)%release()
1233 DEALLOCATE (vxc_rho)
1234 IF (
ASSOCIATED(vxc_tau))
THEN
1236 CALL vxc_tau(ispin)%release()
1238 DEALLOCATE (vxc_tau)
1242 NULLIFY (vxc_rho, rho_r, rho_g, vxc_tau, tau_r, tau)
1244 CALL pw_pool%create_pw(rho_r(1))
1245 IF (
ASSOCIATED(tau1_r) .AND.
ASSOCIATED(v_tau))
THEN
1247 CALL pw_pool%create_pw(tau_r(1))
1249 CALL xc_rho_set_get(rho_set, can_return_null=.true., rho=rho, tau=tau)
1250 DO istep = -nsteps, nsteps
1251 IF (istep == 0) cycle
1252 rweight = rweights(istep, nsteps)/h
1253 step = real(istep,
dp)*h
1256 rho_r(1)%array(:, :, :) = rho(:, :, :) + step*rho1_r(1)%array(:, :, :)
1258 IF (
ASSOCIATED(tau1_r) .AND.
ASSOCIATED(tau) .AND.
ASSOCIATED(tau1_r))
THEN
1260 tau_r(1)%array(:, :, :) = tau(:, :, :) + step*tau1_r(1)%array(:, :, :)
1264 CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
1265 weights, pw_pool, .false., virial_dummy)
1266 CALL pw_axpy(vxc_rho(1), v_xc(1), rweight)
1267 IF (
ASSOCIATED(vxc_tau) .AND.
ASSOCIATED(v_tau))
THEN
1268 CALL pw_axpy(vxc_tau(1), v_tau(1), rweight)
1270 CALL vxc_rho(1)%release()
1271 DEALLOCATE (vxc_rho)
1272 IF (
ASSOCIATED(vxc_tau))
THEN
1273 CALL vxc_tau(1)%release()
1274 DEALLOCATE (vxc_tau)
1279 IF (my_calc_virial)
THEN
1281 IF (nspins == 1 .AND. do_triplet)
THEN
1285 CALL check_for_derivatives(deriv_set, (nspins == 2), rho_f, gradient_f, tau_f, laplace_f)
1291 IF (gradient_f)
THEN
1292 bo = rho_set%local_bounds
1295 CALL allocate_pw(virial_pw, pw_pool, bo)
1302 drho_cutoff=gradient_cut, &
1315 CALL xc_rho_set_get(rho_set, drhoa=drhoa, drhob=drhob, norm_drho=norm_drho, &
1316 norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, tau_a=tau_a, tau_b=tau_b, &
1317 laplace_rhoa=laplacea, laplace_rhob=laplaceb, can_return_null=.true.)
1318 CALL xc_rho_set_get(rho1_set, rhoa=rho1a, rhob=rho1b, drhoa=drho1a, drhob=drho1b, laplace_rhoa=laplace1a, &
1319 laplace_rhob=laplace1b, can_return_null=.true.)
1321 CALL calc_drho_from_ab(drho, drhoa, drhob)
1322 CALL calc_drho_from_ab(drho1, drho1a, drho1b)
1324 CALL xc_rho_set_get(rho_set, drho=drho, norm_drho=norm_drho, tau=tau, laplace_rho=laplace, can_return_null=.true.)
1325 CALL xc_rho_set_get(rho1_set, rho=rho1, drho=drho1, laplace_rho=laplace1, can_return_null=.true.)
1328 CALL prepare_dr1dr(dr1dr, drho, drho1)
1331 CALL prepare_dr1dr(dra1dra, drhoa, drho1a)
1332 CALL prepare_dr1dr(drb1drb, drhob, drho1b)
1334 CALL allocate_pw(v_drho, pw_pool, bo)
1335 CALL allocate_pw(v_drhoa, pw_pool, bo)
1336 CALL allocate_pw(v_drhob, pw_pool, bo)
1338 IF (
ASSOCIATED(norm_drhoa))
CALL apply_drho(deriv_set, [
deriv_norm_drhoa], virial_pw, &
1339 drhoa, drho1a, virial_xc, &
1340 norm_drhoa, gradient_cut, dra1dra, v_drhoa%array)
1341 IF (
ASSOCIATED(norm_drhob))
CALL apply_drho(deriv_set, [
deriv_norm_drhob], virial_pw, &
1342 drhob, drho1b, virial_xc, &
1343 norm_drhob, gradient_cut, drb1drb, v_drhob%array)
1344 IF (
ASSOCIATED(norm_drho))
CALL apply_drho(deriv_set, [
deriv_norm_drho], virial_pw, &
1345 drho, drho1, virial_xc, &
1346 norm_drho, gradient_cut, dr1dr, v_drho%array)
1349 cpassert(
ASSOCIATED(deriv_data))
1350 virial_pw%array(:, :, :) = -rho1a(:, :, :)
1351 CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
1353 CALL allocate_pw(v_laplacea, pw_pool, bo)
1356 cpassert(
ASSOCIATED(deriv_data))
1357 virial_pw%array(:, :, :) = -rho1b(:, :, :)
1358 CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
1360 CALL allocate_pw(v_laplaceb, pw_pool, bo)
1366 CALL allocate_pw(v_drho, pw_pool, bo)
1368 CALL apply_drho(deriv_set, [
deriv_norm_drho], virial_pw, drho, drho1, virial_xc, &
1369 norm_drho, gradient_cut, dr1dr, v_drho%array)
1372 cpassert(
ASSOCIATED(deriv_data))
1373 virial_pw%array(:, :, :) = -rho1(:, :, :)
1374 CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
1376 CALL allocate_pw(v_laplace, pw_pool, bo)
1382 rho_r(1)%array = rhoa
1383 rho_r(2)%array = rhob
1385 rho_r(1)%array = rho
1387 IF (
ASSOCIATED(tau1_r))
THEN
1389 tau_r(1)%array = tau_a
1390 tau_r(2)%array = tau_b
1392 tau_r(1)%array = tau
1403 rho_cutoff=rho_cutoff, &
1414 CALL xc_rho_set_get(rho1_set, rhoa=rho1a, rhob=rho1b, tau_a=tau1a, tau_b=tau1b, &
1415 laplace_rhoa=laplace1a, laplace_rhob=laplace1b, can_return_null=.true.)
1416 CALL xc_rho_set_get(rho2_set, norm_drhoa=norm_drho2a, norm_drhob=norm_drho2b, &
1417 norm_drho=norm_drho2, laplace_rhoa=laplace2a, laplace_rhob=laplace2b, can_return_null=.true.)
1419 DO istep = -nsteps, nsteps
1420 IF (istep == 0) cycle
1421 rweight = rweights(istep, nsteps)/h
1422 step = real(istep,
dp)*h
1423 IF (
ASSOCIATED(norm_drhoa))
THEN
1424 CALL get_derivs_rho(norm_drho2a, norm_drhoa, step, xc_fun_section, lsd, rho2_set, deriv_set1)
1425 CALL update_deriv_rho(deriv_set1, [
deriv_rhoa], bo, &
1426 norm_drhoa, gradient_cut, rweight, rho1a, v_drhoa%array)
1427 CALL update_deriv_rho(deriv_set1, [
deriv_rhob], bo, &
1428 norm_drhoa, gradient_cut, rweight, rho1b, v_drhoa%array)
1430 norm_drhoa, gradient_cut, rweight, dra1dra, v_drhoa%array)
1432 norm_drhoa, gradient_cut, rweight, dra1dra, drb1drb, v_drhoa%array, v_drhob%array)
1434 norm_drhoa, gradient_cut, rweight, dra1dra, dr1dr, v_drhoa%array, v_drho%array)
1436 CALL update_deriv_rho(deriv_set1, [
deriv_tau_a], bo, &
1437 norm_drhoa, gradient_cut, rweight, tau1a, v_drhoa%array)
1438 CALL update_deriv_rho(deriv_set1, [
deriv_tau_b], bo, &
1439 norm_drhoa, gradient_cut, rweight, tau1b, v_drhoa%array)
1443 norm_drhoa, gradient_cut, rweight, laplace1a, v_drhoa%array)
1445 norm_drhoa, gradient_cut, rweight, laplace1b, v_drhoa%array)
1449 IF (
ASSOCIATED(norm_drhob))
THEN
1450 CALL get_derivs_rho(norm_drho2b, norm_drhob, step, xc_fun_section, lsd, rho2_set, deriv_set1)
1451 CALL update_deriv_rho(deriv_set1, [
deriv_rhoa], bo, &
1452 norm_drhob, gradient_cut, rweight, rho1a, v_drhob%array)
1453 CALL update_deriv_rho(deriv_set1, [
deriv_rhob], bo, &
1454 norm_drhob, gradient_cut, rweight, rho1b, v_drhob%array)
1456 norm_drhob, gradient_cut, rweight, drb1drb, v_drhob%array)
1458 norm_drhob, gradient_cut, rweight, drb1drb, dra1dra, v_drhob%array, v_drhoa%array)
1460 norm_drhob, gradient_cut, rweight, drb1drb, dr1dr, v_drhob%array, v_drho%array)
1462 CALL update_deriv_rho(deriv_set1, [
deriv_tau_a], bo, &
1463 norm_drhob, gradient_cut, rweight, tau1a, v_drhob%array)
1464 CALL update_deriv_rho(deriv_set1, [
deriv_tau_b], bo, &
1465 norm_drhob, gradient_cut, rweight, tau1b, v_drhob%array)
1469 norm_drhob, gradient_cut, rweight, laplace1a, v_drhob%array)
1471 norm_drhob, gradient_cut, rweight, laplace1b, v_drhob%array)
1475 IF (
ASSOCIATED(norm_drho))
THEN
1476 CALL get_derivs_rho(norm_drho2, norm_drho, step, xc_fun_section, lsd, rho2_set, deriv_set1)
1477 CALL update_deriv_rho(deriv_set1, [
deriv_rhoa], bo, &
1478 norm_drho, gradient_cut, rweight, rho1a, v_drho%array)
1479 CALL update_deriv_rho(deriv_set1, [
deriv_rhob], bo, &
1480 norm_drho, gradient_cut, rweight, rho1b, v_drho%array)
1482 norm_drho, gradient_cut, rweight, dr1dr, v_drho%array)
1484 norm_drho, gradient_cut, rweight, dr1dr, dra1dra, v_drho%array, v_drhoa%array)
1486 norm_drho, gradient_cut, rweight, dr1dr, drb1drb, v_drho%array, v_drhob%array)
1488 CALL update_deriv_rho(deriv_set1, [
deriv_tau_a], bo, &
1489 norm_drho, gradient_cut, rweight, tau1a, v_drho%array)
1490 CALL update_deriv_rho(deriv_set1, [
deriv_tau_b], bo, &
1491 norm_drho, gradient_cut, rweight, tau1b, v_drho%array)
1495 norm_drho, gradient_cut, rweight, laplace1a, v_drho%array)
1497 norm_drho, gradient_cut, rweight, laplace1b, v_drho%array)
1503 CALL get_derivs_rho(laplace2a, laplacea, step, xc_fun_section, lsd, rho2_set, deriv_set1)
1506 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_rhoa], bo, &
1507 rweight, rho1a, v_laplacea%array)
1508 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_rhob], bo, &
1509 rweight, rho1b, v_laplacea%array)
1510 IF (
ASSOCIATED(norm_drho))
THEN
1511 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_norm_drho], bo, &
1512 rweight, dr1dr, v_laplacea%array)
1514 IF (
ASSOCIATED(norm_drhoa))
THEN
1515 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_norm_drhoa], bo, &
1516 rweight, dra1dra, v_laplacea%array)
1518 IF (
ASSOCIATED(norm_drhob))
THEN
1519 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_norm_drhob], bo, &
1520 rweight, drb1drb, v_laplacea%array)
1523 IF (
ASSOCIATED(tau1a))
THEN
1524 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_tau_a], bo, &
1525 rweight, tau1a, v_laplacea%array)
1527 IF (
ASSOCIATED(tau1b))
THEN
1528 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_tau_b], bo, &
1529 rweight, tau1b, v_laplacea%array)
1533 rweight, laplace1a, v_laplacea%array)
1536 rweight, laplace1b, v_laplacea%array)
1539 CALL get_derivs_rho(laplace2b, laplaceb, step, xc_fun_section, lsd, rho2_set, deriv_set1)
1542 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_rhoa], bo, &
1543 rweight, rho1a, v_laplaceb%array)
1544 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_rhob], bo, &
1545 rweight, rho1b, v_laplaceb%array)
1546 IF (
ASSOCIATED(norm_drho))
THEN
1547 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_norm_drho], bo, &
1548 rweight, dr1dr, v_laplaceb%array)
1550 IF (
ASSOCIATED(norm_drhoa))
THEN
1551 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_norm_drhoa], bo, &
1552 rweight, dra1dra, v_laplaceb%array)
1554 IF (
ASSOCIATED(norm_drhob))
THEN
1555 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_norm_drhob], bo, &
1556 rweight, drb1drb, v_laplaceb%array)
1560 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_tau_a], bo, &
1561 rweight, tau1a, v_laplaceb%array)
1562 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_tau_b], bo, &
1563 rweight, tau1b, v_laplaceb%array)
1567 rweight, laplace1a, v_laplaceb%array)
1570 rweight, laplace1b, v_laplaceb%array)
1574 CALL virial_drho_drho(virial_pw, drhoa, v_drhoa, virial_xc)
1575 CALL virial_drho_drho(virial_pw, drhob, v_drhob, virial_xc)
1576 CALL virial_drho_drho(virial_pw, drho, v_drho, virial_xc)
1578 CALL deallocate_pw(v_drho, pw_pool)
1579 CALL deallocate_pw(v_drhoa, pw_pool)
1580 CALL deallocate_pw(v_drhob, pw_pool)
1583 virial_pw%array(:, :, :) = -rhoa(:, :, :)
1584 CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplacea%array)
1585 CALL deallocate_pw(v_laplacea, pw_pool)
1587 virial_pw%array(:, :, :) = -rhob(:, :, :)
1588 CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplaceb%array)
1589 CALL deallocate_pw(v_laplaceb, pw_pool)
1592 CALL deallocate_pw(virial_pw, pw_pool)
1595 DEALLOCATE (drho(idir)%array)
1596 DEALLOCATE (drho1(idir)%array)
1598 DEALLOCATE (dra1dra, drb1drb)
1601 CALL xc_rho_set_get(rho1_set, rho=rho1, tau=tau1, laplace_rho=laplace1, can_return_null=.true.)
1602 CALL xc_rho_set_get(rho2_set, norm_drho=norm_drho2, laplace_rho=laplace2, can_return_null=.true.)
1604 DO istep = -nsteps, nsteps
1605 IF (istep == 0) cycle
1606 rweight = rweights(istep, nsteps)/h
1607 step = real(istep,
dp)*h
1608 CALL get_derivs_rho(norm_drho2, norm_drho, step, xc_fun_section, lsd, rho2_set, deriv_set1)
1611 CALL update_deriv_rho(deriv_set1, [
deriv_rho], bo, &
1612 norm_drho, gradient_cut, rweight, rho1, v_drho%array)
1614 norm_drho, gradient_cut, rweight, dr1dr, v_drho%array)
1617 CALL update_deriv_rho(deriv_set1, [
deriv_tau], bo, &
1618 norm_drho, gradient_cut, rweight, tau1, v_drho%array)
1622 norm_drho, gradient_cut, rweight, laplace1, v_drho%array)
1624 CALL get_derivs_rho(laplace2, laplace, step, xc_fun_section, lsd, rho2_set, deriv_set1)
1627 CALL update_deriv(deriv_set1, laplace, rho_cutoff, [
deriv_rho], bo, &
1628 rweight, rho1, v_laplace%array)
1629 CALL update_deriv(deriv_set1, laplace, rho_cutoff, [
deriv_norm_drho], bo, &
1630 rweight, dr1dr, v_laplace%array)
1633 CALL update_deriv(deriv_set1, laplace, rho_cutoff, [
deriv_tau], bo, &
1634 rweight, tau1, v_laplace%array)
1638 rweight, laplace1, v_laplace%array)
1643 CALL virial_drho_drho(virial_pw, drho, v_drho, virial_xc)
1645 CALL deallocate_pw(v_drho, pw_pool)
1648 virial_pw%array(:, :, :) = -rho(:, :, :)
1649 CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace%array)
1650 CALL deallocate_pw(v_laplace, pw_pool)
1653 CALL deallocate_pw(virial_pw, pw_pool)
1666 DO ispin = 1,
SIZE(rho_r)
1667 CALL pw_pool%give_back_pw(rho_r(ispin))
1671 IF (
ASSOCIATED(tau_r))
THEN
1672 DO ispin = 1,
SIZE(tau_r)
1673 CALL pw_pool%give_back_pw(tau_r(ispin))
1678 CALL timestop(handle)
2053 LOGICAL,
INTENT(IN),
OPTIONAL :: gapw
2054 REAL(kind=
dp),
DIMENSION(:, :, :, :),
OPTIONAL, &
2056 REAL(kind=
dp),
INTENT(in),
OPTIONAL :: tddfpt_fac
2057 LOGICAL,
INTENT(IN),
OPTIONAL :: compute_virial
2058 REAL(kind=
dp),
DIMENSION(3, 3),
INTENT(INOUT), &
2059 OPTIONAL :: virial_xc
2060 LOGICAL,
INTENT(in),
OPTIONAL :: spinflip
2062 CHARACTER(len=*),
PARAMETER :: routinen =
'xc_calc_2nd_deriv_analytical'
2064 INTEGER :: handle, i, ia, idir, ir, ispin, j, jdir, &
2065 k, nspins, xc_deriv_method_id
2066 INTEGER,
DIMENSION(2, 3) :: bo
2067 LOGICAL :: gradient_f, lsd, my_compute_virial, alda0, &
2068 my_gapw, tau_f, laplace_f, rho_f, do_spinflip
2069 REAL(kind=
dp) ::
fac, gradient_cut, tmp, factor2, s, s_thresh
2070 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: dr1dr, dra1dra, drb1drb
2071 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: deriv_data, deriv_data2, &
2072 e_drhoa, e_drhob, e_drho, norm_drho, norm_drhoa, &
2073 norm_drhob, rho1, rho1a, rho1b, &
2074 tau1, tau1a, tau1b, laplace1, laplace1a, laplace1b, &
2076 TYPE(
cp_3d_r_cp_type),
DIMENSION(3) :: drho, drho1, drho1a, drho1b, drhoa, drhob
2077 TYPE(
pw_r3d_rs_type),
DIMENSION(:),
ALLOCATABLE :: v_drhoa, v_drhob, v_drho, v_laplace
2083 CALL timeset(routinen, handle)
2085 NULLIFY (e_drhoa, e_drhob, e_drho)
2088 IF (
PRESENT(gapw)) my_gapw = gapw
2090 my_compute_virial = .false.
2091 IF (
PRESENT(compute_virial)) my_compute_virial = compute_virial
2093 cpassert(
ASSOCIATED(v_xc))
2094 cpassert(
ASSOCIATED(xc_section))
2096 cpassert(
PRESENT(vxg))
2098 IF (my_compute_virial)
THEN
2099 cpassert(
PRESENT(virial_xc))
2103 i_val=xc_deriv_method_id)
2106 lsd =
ASSOCIATED(rho_set%rhoa)
2109 IF (
PRESENT(tddfpt_fac))
fac = tddfpt_fac
2110 IF (
PRESENT(tddfpt_fac)) factor2 = tddfpt_fac
2111 do_spinflip = .false.
2112 IF (
PRESENT(spinflip)) do_spinflip = spinflip
2114 bo = rho_set%local_bounds
2116 CALL check_for_derivatives(deriv_set, lsd, rho_f, gradient_f, tau_f, laplace_f)
2119 IF (gradient_f)
THEN
2126 cpassert(
ASSOCIATED(v_xc_tau))
2129 IF (gradient_f)
THEN
2130 ALLOCATE (v_drho_r(3, nspins), v_drho(nspins))
2131 DO ispin = 1, nspins
2133 CALL allocate_pw(v_drho_r(idir, ispin), pw_pool, bo)
2135 CALL allocate_pw(v_drho(ispin), pw_pool, bo)
2139 IF (
ASSOCIATED(pw_pool))
THEN
2140 CALL pw_pool%create_pw(tmp_g)
2141 CALL pw_pool%create_pw(vxc_g)
2144 cpabort(
"XC_DERIV method is not implemented in GAPW")
2149 DO ispin = 1, nspins
2150 v_xc(ispin)%array = 0.0_dp
2154 DO ispin = 1, nspins
2155 v_xc_tau(ispin)%array = 0.0_dp
2159 IF (laplace_f .AND. my_gapw) &
2160 cpabort(
"Laplace-dependent functional not implemented with GAPW!")
2162 IF (my_compute_virial .AND. (gradient_f .OR. laplace_f))
CALL allocate_pw(virial_pw, pw_pool, bo)
2170 IF (do_spinflip)
THEN
2177 IF (gradient_f)
THEN
2179 norm_drho=norm_drho, norm_drhoa=norm_drhoa, norm_drhob=norm_drhob)
2180 IF (do_spinflip)
THEN
2182 CALL calc_drho_from_a(drho1, drho1a)
2185 CALL calc_drho_from_ab(drho1, drho1a, drho1b)
2188 CALL calc_drho_from_ab(drho, drhoa, drhob)
2190 CALL prepare_dr1dr(dra1dra, drhoa, drho1a)
2191 IF (do_spinflip)
THEN
2192 CALL prepare_dr1dr(drb1drb, drhob, drho1a)
2193 CALL prepare_dr1dr(dr1dr, drho, drho1a)
2194 ELSE IF (nspins /= 1)
THEN
2195 CALL prepare_dr1dr(drb1drb, drhob, drho1b)
2196 CALL prepare_dr1dr(dr1dr, drho, drho1)
2198 CALL prepare_dr1dr(drb1drb, drhob, drho1b)
2199 CALL prepare_dr1dr_ab(dr1dr, drhoa, drhob, drho1a, drho1b,
fac)
2202 ALLOCATE (v_drhoa(nspins), v_drhob(nspins))
2203 DO ispin = 1, nspins
2204 CALL allocate_pw(v_drhoa(ispin), pw_pool, bo)
2205 CALL allocate_pw(v_drhob(ispin), pw_pool, bo)
2211 CALL xc_rho_set_get(rho1_set, laplace_rhoa=laplace1a, laplace_rhob=laplace1b)
2213 ALLOCATE (v_laplace(nspins))
2214 DO ispin = 1, nspins
2215 CALL allocate_pw(v_laplace(ispin), pw_pool, bo)
2218 IF (my_compute_virial)
CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob)
2225 IF (do_spinflip)
THEN
2234 IF (
ASSOCIATED(deriv_att))
THEN
2237 IF (
ASSOCIATED(deriv_att))
THEN
2241 DO k = bo(1, 3), bo(2, 3)
2242 DO j = bo(1, 2), bo(2, 2)
2243 DO i = bo(1, 1), bo(2, 1)
2244 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
2245 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2246 (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)/s
2258 IF (.NOT. alda0)
THEN
2260 IF (
ASSOCIATED(deriv_att))
THEN
2263 IF (
ASSOCIATED(deriv_att))
THEN
2267 DO k = bo(1, 3), bo(2, 3)
2268 DO j = bo(1, 2), bo(2, 2)
2269 DO i = bo(1, 1), bo(2, 1)
2270 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
2271 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2272 (deriv_data(i, j, k)*dra1dra(i, j, k) - &
2273 deriv_data2(i, j, k)*drb1drb(i, j, k))/s
2282 ELSE IF (nspins /= 1)
THEN
2286 IF (
ASSOCIATED(deriv_att))
THEN
2290 DO k = bo(1, 3), bo(2, 3)
2291 DO j = bo(1, 2), bo(2, 2)
2292 DO i = bo(1, 1), bo(2, 1)
2293 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2294 deriv_data(i, j, k)*rho1a(i, j, k)
2300 IF (
ASSOCIATED(deriv_att))
THEN
2304 DO k = bo(1, 3), bo(2, 3)
2305 DO j = bo(1, 2), bo(2, 2)
2306 DO i = bo(1, 1), bo(2, 1)
2307 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2308 deriv_data(i, j, k)*rho1b(i, j, k)
2314 IF (
ASSOCIATED(deriv_att))
THEN
2318 DO k = bo(1, 3), bo(2, 3)
2319 DO j = bo(1, 2), bo(2, 2)
2320 DO i = bo(1, 1), bo(2, 1)
2321 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2322 deriv_data(i, j, k)*dr1dr(i, j, k)
2328 IF (
ASSOCIATED(deriv_att))
THEN
2332 DO k = bo(1, 3), bo(2, 3)
2333 DO j = bo(1, 2), bo(2, 2)
2334 DO i = bo(1, 1), bo(2, 1)
2335 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2336 deriv_data(i, j, k)*dra1dra(i, j, k)
2342 IF (
ASSOCIATED(deriv_att))
THEN
2346 DO k = bo(1, 3), bo(2, 3)
2347 DO j = bo(1, 2), bo(2, 2)
2348 DO i = bo(1, 1), bo(2, 1)
2349 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2350 deriv_data(i, j, k)*drb1drb(i, j, k)
2356 IF (
ASSOCIATED(deriv_att))
THEN
2360 DO k = bo(1, 3), bo(2, 3)
2361 DO j = bo(1, 2), bo(2, 2)
2362 DO i = bo(1, 1), bo(2, 1)
2363 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2364 deriv_data(i, j, k)*tau1a(i, j, k)
2370 IF (
ASSOCIATED(deriv_att))
THEN
2374 DO k = bo(1, 3), bo(2, 3)
2375 DO j = bo(1, 2), bo(2, 2)
2376 DO i = bo(1, 1), bo(2, 1)
2377 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2378 deriv_data(i, j, k)*tau1b(i, j, k)
2384 IF (
ASSOCIATED(deriv_att))
THEN
2388 DO k = bo(1, 3), bo(2, 3)
2389 DO j = bo(1, 2), bo(2, 2)
2390 DO i = bo(1, 1), bo(2, 1)
2391 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2392 deriv_data(i, j, k)*laplace1a(i, j, k)
2398 IF (
ASSOCIATED(deriv_att))
THEN
2402 DO k = bo(1, 3), bo(2, 3)
2403 DO j = bo(1, 2), bo(2, 2)
2404 DO i = bo(1, 1), bo(2, 1)
2405 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2406 deriv_data(i, j, k)*laplace1b(i, j, k)
2414 IF (
ASSOCIATED(deriv_att))
THEN
2418 DO k = bo(1, 3), bo(2, 3)
2419 DO j = bo(1, 2), bo(2, 2)
2420 DO i = bo(1, 1), bo(2, 1)
2421 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2422 deriv_data(i, j, k)*rho1a(i, j, k)
2428 IF (
ASSOCIATED(deriv_att))
THEN
2432 DO k = bo(1, 3), bo(2, 3)
2433 DO j = bo(1, 2), bo(2, 2)
2434 DO i = bo(1, 1), bo(2, 1)
2435 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2436 deriv_data(i, j, k)*rho1b(i, j, k)
2442 IF (
ASSOCIATED(deriv_att))
THEN
2446 DO k = bo(1, 3), bo(2, 3)
2447 DO j = bo(1, 2), bo(2, 2)
2448 DO i = bo(1, 1), bo(2, 1)
2449 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2450 deriv_data(i, j, k)*dr1dr(i, j, k)
2456 IF (
ASSOCIATED(deriv_att))
THEN
2460 DO k = bo(1, 3), bo(2, 3)
2461 DO j = bo(1, 2), bo(2, 2)
2462 DO i = bo(1, 1), bo(2, 1)
2463 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2464 deriv_data(i, j, k)*dra1dra(i, j, k)
2470 IF (
ASSOCIATED(deriv_att))
THEN
2474 DO k = bo(1, 3), bo(2, 3)
2475 DO j = bo(1, 2), bo(2, 2)
2476 DO i = bo(1, 1), bo(2, 1)
2477 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2478 deriv_data(i, j, k)*drb1drb(i, j, k)
2484 IF (
ASSOCIATED(deriv_att))
THEN
2488 DO k = bo(1, 3), bo(2, 3)
2489 DO j = bo(1, 2), bo(2, 2)
2490 DO i = bo(1, 1), bo(2, 1)
2491 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2492 deriv_data(i, j, k)*tau1a(i, j, k)
2498 IF (
ASSOCIATED(deriv_att))
THEN
2502 DO k = bo(1, 3), bo(2, 3)
2503 DO j = bo(1, 2), bo(2, 2)
2504 DO i = bo(1, 1), bo(2, 1)
2505 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2506 deriv_data(i, j, k)*tau1b(i, j, k)
2512 IF (
ASSOCIATED(deriv_att))
THEN
2516 DO k = bo(1, 3), bo(2, 3)
2517 DO j = bo(1, 2), bo(2, 2)
2518 DO i = bo(1, 1), bo(2, 1)
2519 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2520 deriv_data(i, j, k)*laplace1a(i, j, k)
2526 IF (
ASSOCIATED(deriv_att))
THEN
2530 DO k = bo(1, 3), bo(2, 3)
2531 DO j = bo(1, 2), bo(2, 2)
2532 DO i = bo(1, 1), bo(2, 1)
2533 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2534 deriv_data(i, j, k)*laplace1b(i, j, k)
2542 IF (
ASSOCIATED(deriv_att))
THEN
2546 DO k = bo(1, 3), bo(2, 3)
2547 DO j = bo(1, 2), bo(2, 2)
2548 DO i = bo(1, 1), bo(2, 1)
2549 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2550 deriv_data(i, j, k)*rho1a(i, j, k)
2551 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2552 deriv_data(i, j, k)*rho1a(i, j, k)
2558 IF (
ASSOCIATED(deriv_att))
THEN
2562 DO k = bo(1, 3), bo(2, 3)
2563 DO j = bo(1, 2), bo(2, 2)
2564 DO i = bo(1, 1), bo(2, 1)
2565 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2566 deriv_data(i, j, k)*rho1b(i, j, k)
2567 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2568 deriv_data(i, j, k)*rho1b(i, j, k)
2574 IF (
ASSOCIATED(deriv_att))
THEN
2578 DO k = bo(1, 3), bo(2, 3)
2579 DO j = bo(1, 2), bo(2, 2)
2580 DO i = bo(1, 1), bo(2, 1)
2581 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2582 deriv_data(i, j, k)*dr1dr(i, j, k)
2583 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2584 deriv_data(i, j, k)*dr1dr(i, j, k)
2590 IF (
ASSOCIATED(deriv_att))
THEN
2594 DO k = bo(1, 3), bo(2, 3)
2595 DO j = bo(1, 2), bo(2, 2)
2596 DO i = bo(1, 1), bo(2, 1)
2597 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2598 deriv_data(i, j, k)*dra1dra(i, j, k)
2599 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2600 deriv_data(i, j, k)*dra1dra(i, j, k)
2606 IF (
ASSOCIATED(deriv_att))
THEN
2610 DO k = bo(1, 3), bo(2, 3)
2611 DO j = bo(1, 2), bo(2, 2)
2612 DO i = bo(1, 1), bo(2, 1)
2613 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2614 deriv_data(i, j, k)*drb1drb(i, j, k)
2615 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2616 deriv_data(i, j, k)*drb1drb(i, j, k)
2622 IF (
ASSOCIATED(deriv_att))
THEN
2626 DO k = bo(1, 3), bo(2, 3)
2627 DO j = bo(1, 2), bo(2, 2)
2628 DO i = bo(1, 1), bo(2, 1)
2629 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2630 deriv_data(i, j, k)*tau1a(i, j, k)
2631 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2632 deriv_data(i, j, k)*tau1a(i, j, k)
2638 IF (
ASSOCIATED(deriv_att))
THEN
2642 DO k = bo(1, 3), bo(2, 3)
2643 DO j = bo(1, 2), bo(2, 2)
2644 DO i = bo(1, 1), bo(2, 1)
2645 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2646 deriv_data(i, j, k)*tau1b(i, j, k)
2647 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2648 deriv_data(i, j, k)*tau1b(i, j, k)
2654 IF (
ASSOCIATED(deriv_att))
THEN
2658 DO k = bo(1, 3), bo(2, 3)
2659 DO j = bo(1, 2), bo(2, 2)
2660 DO i = bo(1, 1), bo(2, 1)
2661 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2662 deriv_data(i, j, k)*laplace1a(i, j, k)
2663 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2664 deriv_data(i, j, k)*laplace1a(i, j, k)
2670 IF (
ASSOCIATED(deriv_att))
THEN
2674 DO k = bo(1, 3), bo(2, 3)
2675 DO j = bo(1, 2), bo(2, 2)
2676 DO i = bo(1, 1), bo(2, 1)
2677 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2678 deriv_data(i, j, k)*laplace1b(i, j, k)
2679 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2680 deriv_data(i, j, k)*laplace1b(i, j, k)
2687 IF (
ASSOCIATED(deriv_att))
THEN
2691 IF (my_compute_virial)
THEN
2692 CALL virial_drho_drho1(virial_pw, drho, drho1, deriv_data, virial_xc)
2696 v_drho(1)%array(:, :, :) = v_drho(1)%array(:, :, :) + &
2697 deriv_data(:, :, :)*dr1dr(:, :, :)/max(gradient_cut, norm_drho(:, :, :))**2
2698 v_drho(2)%array(:, :, :) = v_drho(2)%array(:, :, :) + &
2699 deriv_data(:, :, :)*dr1dr(:, :, :)/max(gradient_cut, norm_drho(:, :, :))**2
2704 IF (
ASSOCIATED(deriv_att))
THEN
2708 DO k = bo(1, 3), bo(2, 3)
2709 DO j = bo(1, 2), bo(2, 2)
2710 DO i = bo(1, 1), bo(2, 1)
2711 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2712 deriv_data(i, j, k)*rho1a(i, j, k)
2718 IF (
ASSOCIATED(deriv_att))
THEN
2722 DO k = bo(1, 3), bo(2, 3)
2723 DO j = bo(1, 2), bo(2, 2)
2724 DO i = bo(1, 1), bo(2, 1)
2725 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2726 deriv_data(i, j, k)*rho1b(i, j, k)
2732 IF (
ASSOCIATED(deriv_att))
THEN
2736 DO k = bo(1, 3), bo(2, 3)
2737 DO j = bo(1, 2), bo(2, 2)
2738 DO i = bo(1, 1), bo(2, 1)
2739 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2740 deriv_data(i, j, k)*dr1dr(i, j, k)
2746 IF (
ASSOCIATED(deriv_att))
THEN
2750 DO k = bo(1, 3), bo(2, 3)
2751 DO j = bo(1, 2), bo(2, 2)
2752 DO i = bo(1, 1), bo(2, 1)
2753 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2754 deriv_data(i, j, k)*dra1dra(i, j, k)
2760 IF (
ASSOCIATED(deriv_att))
THEN
2764 DO k = bo(1, 3), bo(2, 3)
2765 DO j = bo(1, 2), bo(2, 2)
2766 DO i = bo(1, 1), bo(2, 1)
2767 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2768 deriv_data(i, j, k)*drb1drb(i, j, k)
2774 IF (
ASSOCIATED(deriv_att))
THEN
2778 DO k = bo(1, 3), bo(2, 3)
2779 DO j = bo(1, 2), bo(2, 2)
2780 DO i = bo(1, 1), bo(2, 1)
2781 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2782 deriv_data(i, j, k)*tau1a(i, j, k)
2788 IF (
ASSOCIATED(deriv_att))
THEN
2792 DO k = bo(1, 3), bo(2, 3)
2793 DO j = bo(1, 2), bo(2, 2)
2794 DO i = bo(1, 1), bo(2, 1)
2795 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2796 deriv_data(i, j, k)*tau1b(i, j, k)
2802 IF (
ASSOCIATED(deriv_att))
THEN
2806 DO k = bo(1, 3), bo(2, 3)
2807 DO j = bo(1, 2), bo(2, 2)
2808 DO i = bo(1, 1), bo(2, 1)
2809 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2810 deriv_data(i, j, k)*laplace1a(i, j, k)
2816 IF (
ASSOCIATED(deriv_att))
THEN
2820 DO k = bo(1, 3), bo(2, 3)
2821 DO j = bo(1, 2), bo(2, 2)
2822 DO i = bo(1, 1), bo(2, 1)
2823 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2824 deriv_data(i, j, k)*laplace1b(i, j, k)
2831 IF (
ASSOCIATED(deriv_att))
THEN
2835 IF (my_compute_virial)
THEN
2836 CALL virial_drho_drho1(virial_pw, drhoa, drho1a, deriv_data, virial_xc)
2840 v_drhoa(1)%array(:, :, :) = v_drhoa(1)%array(:, :, :) + &
2841 deriv_data(:, :, :)*dra1dra(:, :, :)/max(gradient_cut, norm_drhoa(:, :, :))**2
2846 IF (
ASSOCIATED(deriv_att))
THEN
2850 DO k = bo(1, 3), bo(2, 3)
2851 DO j = bo(1, 2), bo(2, 2)
2852 DO i = bo(1, 1), bo(2, 1)
2853 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2854 deriv_data(i, j, k)*rho1a(i, j, k)
2860 IF (
ASSOCIATED(deriv_att))
THEN
2864 DO k = bo(1, 3), bo(2, 3)
2865 DO j = bo(1, 2), bo(2, 2)
2866 DO i = bo(1, 1), bo(2, 1)
2867 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2868 deriv_data(i, j, k)*rho1b(i, j, k)
2874 IF (
ASSOCIATED(deriv_att))
THEN
2878 DO k = bo(1, 3), bo(2, 3)
2879 DO j = bo(1, 2), bo(2, 2)
2880 DO i = bo(1, 1), bo(2, 1)
2881 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2882 deriv_data(i, j, k)*dr1dr(i, j, k)
2888 IF (
ASSOCIATED(deriv_att))
THEN
2892 DO k = bo(1, 3), bo(2, 3)
2893 DO j = bo(1, 2), bo(2, 2)
2894 DO i = bo(1, 1), bo(2, 1)
2895 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2896 deriv_data(i, j, k)*dra1dra(i, j, k)
2902 IF (
ASSOCIATED(deriv_att))
THEN
2906 DO k = bo(1, 3), bo(2, 3)
2907 DO j = bo(1, 2), bo(2, 2)
2908 DO i = bo(1, 1), bo(2, 1)
2909 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2910 deriv_data(i, j, k)*drb1drb(i, j, k)
2916 IF (
ASSOCIATED(deriv_att))
THEN
2920 DO k = bo(1, 3), bo(2, 3)
2921 DO j = bo(1, 2), bo(2, 2)
2922 DO i = bo(1, 1), bo(2, 1)
2923 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2924 deriv_data(i, j, k)*tau1a(i, j, k)
2930 IF (
ASSOCIATED(deriv_att))
THEN
2934 DO k = bo(1, 3), bo(2, 3)
2935 DO j = bo(1, 2), bo(2, 2)
2936 DO i = bo(1, 1), bo(2, 1)
2937 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2938 deriv_data(i, j, k)*tau1b(i, j, k)
2944 IF (
ASSOCIATED(deriv_att))
THEN
2948 DO k = bo(1, 3), bo(2, 3)
2949 DO j = bo(1, 2), bo(2, 2)
2950 DO i = bo(1, 1), bo(2, 1)
2951 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2952 deriv_data(i, j, k)*laplace1a(i, j, k)
2958 IF (
ASSOCIATED(deriv_att))
THEN
2962 DO k = bo(1, 3), bo(2, 3)
2963 DO j = bo(1, 2), bo(2, 2)
2964 DO i = bo(1, 1), bo(2, 1)
2965 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2966 deriv_data(i, j, k)*laplace1b(i, j, k)
2973 IF (
ASSOCIATED(deriv_att))
THEN
2977 IF (my_compute_virial)
THEN
2978 CALL virial_drho_drho1(virial_pw, drhob, drho1b, deriv_data, virial_xc)
2982 v_drhob(2)%array(:, :, :) = v_drhob(2)%array(:, :, :) + &
2983 deriv_data(:, :, :)*drb1drb(:, :, :)/max(gradient_cut, norm_drhob(:, :, :))**2
2988 IF (
ASSOCIATED(deriv_att))
THEN
2992 DO k = bo(1, 3), bo(2, 3)
2993 DO j = bo(1, 2), bo(2, 2)
2994 DO i = bo(1, 1), bo(2, 1)
2995 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
2996 deriv_data(i, j, k)*rho1a(i, j, k)
3002 IF (
ASSOCIATED(deriv_att))
THEN
3006 DO k = bo(1, 3), bo(2, 3)
3007 DO j = bo(1, 2), bo(2, 2)
3008 DO i = bo(1, 1), bo(2, 1)
3009 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3010 deriv_data(i, j, k)*rho1b(i, j, k)
3016 IF (
ASSOCIATED(deriv_att))
THEN
3020 DO k = bo(1, 3), bo(2, 3)
3021 DO j = bo(1, 2), bo(2, 2)
3022 DO i = bo(1, 1), bo(2, 1)
3023 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3024 deriv_data(i, j, k)*dr1dr(i, j, k)
3030 IF (
ASSOCIATED(deriv_att))
THEN
3034 DO k = bo(1, 3), bo(2, 3)
3035 DO j = bo(1, 2), bo(2, 2)
3036 DO i = bo(1, 1), bo(2, 1)
3037 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3038 deriv_data(i, j, k)*dra1dra(i, j, k)
3044 IF (
ASSOCIATED(deriv_att))
THEN
3048 DO k = bo(1, 3), bo(2, 3)
3049 DO j = bo(1, 2), bo(2, 2)
3050 DO i = bo(1, 1), bo(2, 1)
3051 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3052 deriv_data(i, j, k)*drb1drb(i, j, k)
3058 IF (
ASSOCIATED(deriv_att))
THEN
3062 DO k = bo(1, 3), bo(2, 3)
3063 DO j = bo(1, 2), bo(2, 2)
3064 DO i = bo(1, 1), bo(2, 1)
3065 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3066 deriv_data(i, j, k)*tau1a(i, j, k)
3072 IF (
ASSOCIATED(deriv_att))
THEN
3076 DO k = bo(1, 3), bo(2, 3)
3077 DO j = bo(1, 2), bo(2, 2)
3078 DO i = bo(1, 1), bo(2, 1)
3079 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3080 deriv_data(i, j, k)*tau1b(i, j, k)
3086 IF (
ASSOCIATED(deriv_att))
THEN
3090 DO k = bo(1, 3), bo(2, 3)
3091 DO j = bo(1, 2), bo(2, 2)
3092 DO i = bo(1, 1), bo(2, 1)
3093 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3094 deriv_data(i, j, k)*laplace1a(i, j, k)
3100 IF (
ASSOCIATED(deriv_att))
THEN
3104 DO k = bo(1, 3), bo(2, 3)
3105 DO j = bo(1, 2), bo(2, 2)
3106 DO i = bo(1, 1), bo(2, 1)
3107 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3108 deriv_data(i, j, k)*laplace1b(i, j, k)
3116 IF (
ASSOCIATED(deriv_att))
THEN
3120 DO k = bo(1, 3), bo(2, 3)
3121 DO j = bo(1, 2), bo(2, 2)
3122 DO i = bo(1, 1), bo(2, 1)
3123 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3124 deriv_data(i, j, k)*rho1a(i, j, k)
3130 IF (
ASSOCIATED(deriv_att))
THEN
3134 DO k = bo(1, 3), bo(2, 3)
3135 DO j = bo(1, 2), bo(2, 2)
3136 DO i = bo(1, 1), bo(2, 1)
3137 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3138 deriv_data(i, j, k)*rho1b(i, j, k)
3144 IF (
ASSOCIATED(deriv_att))
THEN
3148 DO k = bo(1, 3), bo(2, 3)
3149 DO j = bo(1, 2), bo(2, 2)
3150 DO i = bo(1, 1), bo(2, 1)
3151 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3152 deriv_data(i, j, k)*dr1dr(i, j, k)
3158 IF (
ASSOCIATED(deriv_att))
THEN
3162 DO k = bo(1, 3), bo(2, 3)
3163 DO j = bo(1, 2), bo(2, 2)
3164 DO i = bo(1, 1), bo(2, 1)
3165 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3166 deriv_data(i, j, k)*dra1dra(i, j, k)
3172 IF (
ASSOCIATED(deriv_att))
THEN
3176 DO k = bo(1, 3), bo(2, 3)
3177 DO j = bo(1, 2), bo(2, 2)
3178 DO i = bo(1, 1), bo(2, 1)
3179 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3180 deriv_data(i, j, k)*drb1drb(i, j, k)
3186 IF (
ASSOCIATED(deriv_att))
THEN
3190 DO k = bo(1, 3), bo(2, 3)
3191 DO j = bo(1, 2), bo(2, 2)
3192 DO i = bo(1, 1), bo(2, 1)
3193 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3194 deriv_data(i, j, k)*tau1a(i, j, k)
3200 IF (
ASSOCIATED(deriv_att))
THEN
3204 DO k = bo(1, 3), bo(2, 3)
3205 DO j = bo(1, 2), bo(2, 2)
3206 DO i = bo(1, 1), bo(2, 1)
3207 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3208 deriv_data(i, j, k)*tau1b(i, j, k)
3214 IF (
ASSOCIATED(deriv_att))
THEN
3218 DO k = bo(1, 3), bo(2, 3)
3219 DO j = bo(1, 2), bo(2, 2)
3220 DO i = bo(1, 1), bo(2, 1)
3221 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3222 deriv_data(i, j, k)*laplace1a(i, j, k)
3228 IF (
ASSOCIATED(deriv_att))
THEN
3232 DO k = bo(1, 3), bo(2, 3)
3233 DO j = bo(1, 2), bo(2, 2)
3234 DO i = bo(1, 1), bo(2, 1)
3235 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3236 deriv_data(i, j, k)*laplace1b(i, j, k)
3244 IF (
ASSOCIATED(deriv_att))
THEN
3248 DO k = bo(1, 3), bo(2, 3)
3249 DO j = bo(1, 2), bo(2, 2)
3250 DO i = bo(1, 1), bo(2, 1)
3251 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3252 deriv_data(i, j, k)*rho1a(i, j, k)
3258 IF (
ASSOCIATED(deriv_att))
THEN
3262 DO k = bo(1, 3), bo(2, 3)
3263 DO j = bo(1, 2), bo(2, 2)
3264 DO i = bo(1, 1), bo(2, 1)
3265 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3266 deriv_data(i, j, k)*rho1b(i, j, k)
3272 IF (
ASSOCIATED(deriv_att))
THEN
3276 DO k = bo(1, 3), bo(2, 3)
3277 DO j = bo(1, 2), bo(2, 2)
3278 DO i = bo(1, 1), bo(2, 1)
3279 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3280 deriv_data(i, j, k)*dr1dr(i, j, k)
3286 IF (
ASSOCIATED(deriv_att))
THEN
3290 DO k = bo(1, 3), bo(2, 3)
3291 DO j = bo(1, 2), bo(2, 2)
3292 DO i = bo(1, 1), bo(2, 1)
3293 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3294 deriv_data(i, j, k)*dra1dra(i, j, k)
3300 IF (
ASSOCIATED(deriv_att))
THEN
3304 DO k = bo(1, 3), bo(2, 3)
3305 DO j = bo(1, 2), bo(2, 2)
3306 DO i = bo(1, 1), bo(2, 1)
3307 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3308 deriv_data(i, j, k)*drb1drb(i, j, k)
3314 IF (
ASSOCIATED(deriv_att))
THEN
3318 DO k = bo(1, 3), bo(2, 3)
3319 DO j = bo(1, 2), bo(2, 2)
3320 DO i = bo(1, 1), bo(2, 1)
3321 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3322 deriv_data(i, j, k)*tau1a(i, j, k)
3328 IF (
ASSOCIATED(deriv_att))
THEN
3332 DO k = bo(1, 3), bo(2, 3)
3333 DO j = bo(1, 2), bo(2, 2)
3334 DO i = bo(1, 1), bo(2, 1)
3335 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3336 deriv_data(i, j, k)*tau1b(i, j, k)
3342 IF (
ASSOCIATED(deriv_att))
THEN
3346 DO k = bo(1, 3), bo(2, 3)
3347 DO j = bo(1, 2), bo(2, 2)
3348 DO i = bo(1, 1), bo(2, 1)
3349 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3350 deriv_data(i, j, k)*laplace1a(i, j, k)
3356 IF (
ASSOCIATED(deriv_att))
THEN
3360 DO k = bo(1, 3), bo(2, 3)
3361 DO j = bo(1, 2), bo(2, 2)
3362 DO i = bo(1, 1), bo(2, 1)
3363 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3364 deriv_data(i, j, k)*laplace1b(i, j, k)
3371 IF (my_compute_virial)
THEN
3373 IF (
ASSOCIATED(deriv_att))
THEN
3376 virial_pw%array(:, :, :) = -rho1a(:, :, :)
3377 CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
3381 IF (
ASSOCIATED(deriv_att))
THEN
3385 DO k = bo(1, 3), bo(2, 3)
3386 DO j = bo(1, 2), bo(2, 2)
3387 DO i = bo(1, 1), bo(2, 1)
3388 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3389 deriv_data(i, j, k)*rho1a(i, j, k)
3395 IF (
ASSOCIATED(deriv_att))
THEN
3399 DO k = bo(1, 3), bo(2, 3)
3400 DO j = bo(1, 2), bo(2, 2)
3401 DO i = bo(1, 1), bo(2, 1)
3402 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3403 deriv_data(i, j, k)*rho1b(i, j, k)
3409 IF (
ASSOCIATED(deriv_att))
THEN
3413 DO k = bo(1, 3), bo(2, 3)
3414 DO j = bo(1, 2), bo(2, 2)
3415 DO i = bo(1, 1), bo(2, 1)
3416 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3417 deriv_data(i, j, k)*dr1dr(i, j, k)
3423 IF (
ASSOCIATED(deriv_att))
THEN
3427 DO k = bo(1, 3), bo(2, 3)
3428 DO j = bo(1, 2), bo(2, 2)
3429 DO i = bo(1, 1), bo(2, 1)
3430 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3431 deriv_data(i, j, k)*dra1dra(i, j, k)
3437 IF (
ASSOCIATED(deriv_att))
THEN
3441 DO k = bo(1, 3), bo(2, 3)
3442 DO j = bo(1, 2), bo(2, 2)
3443 DO i = bo(1, 1), bo(2, 1)
3444 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3445 deriv_data(i, j, k)*drb1drb(i, j, k)
3451 IF (
ASSOCIATED(deriv_att))
THEN
3455 DO k = bo(1, 3), bo(2, 3)
3456 DO j = bo(1, 2), bo(2, 2)
3457 DO i = bo(1, 1), bo(2, 1)
3458 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3459 deriv_data(i, j, k)*tau1a(i, j, k)
3465 IF (
ASSOCIATED(deriv_att))
THEN
3469 DO k = bo(1, 3), bo(2, 3)
3470 DO j = bo(1, 2), bo(2, 2)
3471 DO i = bo(1, 1), bo(2, 1)
3472 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3473 deriv_data(i, j, k)*tau1b(i, j, k)
3479 IF (
ASSOCIATED(deriv_att))
THEN
3483 DO k = bo(1, 3), bo(2, 3)
3484 DO j = bo(1, 2), bo(2, 2)
3485 DO i = bo(1, 1), bo(2, 1)
3486 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3487 deriv_data(i, j, k)*laplace1a(i, j, k)
3493 IF (
ASSOCIATED(deriv_att))
THEN
3497 DO k = bo(1, 3), bo(2, 3)
3498 DO j = bo(1, 2), bo(2, 2)
3499 DO i = bo(1, 1), bo(2, 1)
3500 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3501 deriv_data(i, j, k)*laplace1b(i, j, k)
3508 IF (my_compute_virial)
THEN
3510 IF (
ASSOCIATED(deriv_att))
THEN
3513 virial_pw%array(:, :, :) = -rho1b(:, :, :)
3514 CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
3523 IF (
ASSOCIATED(deriv_att))
THEN
3527 DO k = bo(1, 3), bo(2, 3)
3528 DO j = bo(1, 2), bo(2, 2)
3529 DO i = bo(1, 1), bo(2, 1)
3530 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3531 deriv_data(i, j, k)*rho1a(i, j, k)
3537 IF (
ASSOCIATED(deriv_att))
THEN
3541 DO k = bo(1, 3), bo(2, 3)
3542 DO j = bo(1, 2), bo(2, 2)
3543 DO i = bo(1, 1), bo(2, 1)
3544 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3545 deriv_data(i, j, k)*dr1dr(i, j, k)
3551 IF (
ASSOCIATED(deriv_att))
THEN
3555 DO k = bo(1, 3), bo(2, 3)
3556 DO j = bo(1, 2), bo(2, 2)
3557 DO i = bo(1, 1), bo(2, 1)
3558 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3559 deriv_data(i, j, k)*dra1dra(i, j, k)
3565 IF (
ASSOCIATED(deriv_att))
THEN
3569 DO k = bo(1, 3), bo(2, 3)
3570 DO j = bo(1, 2), bo(2, 2)
3571 DO i = bo(1, 1), bo(2, 1)
3572 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3573 deriv_data(i, j, k)*tau1a(i, j, k)
3579 IF (
ASSOCIATED(deriv_att))
THEN
3583 DO k = bo(1, 3), bo(2, 3)
3584 DO j = bo(1, 2), bo(2, 2)
3585 DO i = bo(1, 1), bo(2, 1)
3586 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3587 deriv_data(i, j, k)*laplace1a(i, j, k)
3593 IF (
ASSOCIATED(deriv_att))
THEN
3597 DO k = bo(1, 3), bo(2, 3)
3598 DO j = bo(1, 2), bo(2, 2)
3599 DO i = bo(1, 1), bo(2, 1)
3600 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3601 fac*deriv_data(i, j, k)*rho1b(i, j, k)
3607 IF (
ASSOCIATED(deriv_att))
THEN
3611 DO k = bo(1, 3), bo(2, 3)
3612 DO j = bo(1, 2), bo(2, 2)
3613 DO i = bo(1, 1), bo(2, 1)
3614 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3615 fac*deriv_data(i, j, k)*drb1drb(i, j, k)
3621 IF (
ASSOCIATED(deriv_att))
THEN
3625 DO k = bo(1, 3), bo(2, 3)
3626 DO j = bo(1, 2), bo(2, 2)
3627 DO i = bo(1, 1), bo(2, 1)
3628 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3629 fac*deriv_data(i, j, k)*tau1b(i, j, k)
3635 IF (
ASSOCIATED(deriv_att))
THEN
3639 DO k = bo(1, 3), bo(2, 3)
3640 DO j = bo(1, 2), bo(2, 2)
3641 DO i = bo(1, 1), bo(2, 1)
3642 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3643 fac*deriv_data(i, j, k)*laplace1b(i, j, k)
3651 IF (
ASSOCIATED(deriv_att))
THEN
3655 DO k = bo(1, 3), bo(2, 3)
3656 DO j = bo(1, 2), bo(2, 2)
3657 DO i = bo(1, 1), bo(2, 1)
3658 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3659 deriv_data(i, j, k)*rho1a(i, j, k)
3665 IF (
ASSOCIATED(deriv_att))
THEN
3669 DO k = bo(1, 3), bo(2, 3)
3670 DO j = bo(1, 2), bo(2, 2)
3671 DO i = bo(1, 1), bo(2, 1)
3672 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3673 deriv_data(i, j, k)*dr1dr(i, j, k)
3679 IF (
ASSOCIATED(deriv_att))
THEN
3683 DO k = bo(1, 3), bo(2, 3)
3684 DO j = bo(1, 2), bo(2, 2)
3685 DO i = bo(1, 1), bo(2, 1)
3686 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3687 deriv_data(i, j, k)*dra1dra(i, j, k)
3693 IF (
ASSOCIATED(deriv_att))
THEN
3697 DO k = bo(1, 3), bo(2, 3)
3698 DO j = bo(1, 2), bo(2, 2)
3699 DO i = bo(1, 1), bo(2, 1)
3700 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3701 deriv_data(i, j, k)*tau1a(i, j, k)
3707 IF (
ASSOCIATED(deriv_att))
THEN
3711 DO k = bo(1, 3), bo(2, 3)
3712 DO j = bo(1, 2), bo(2, 2)
3713 DO i = bo(1, 1), bo(2, 1)
3714 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3715 deriv_data(i, j, k)*laplace1a(i, j, k)
3721 IF (
ASSOCIATED(deriv_att))
THEN
3725 DO k = bo(1, 3), bo(2, 3)
3726 DO j = bo(1, 2), bo(2, 2)
3727 DO i = bo(1, 1), bo(2, 1)
3728 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3729 fac*deriv_data(i, j, k)*rho1b(i, j, k)
3735 IF (
ASSOCIATED(deriv_att))
THEN
3739 DO k = bo(1, 3), bo(2, 3)
3740 DO j = bo(1, 2), bo(2, 2)
3741 DO i = bo(1, 1), bo(2, 1)
3742 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3743 fac*deriv_data(i, j, k)*drb1drb(i, j, k)
3749 IF (
ASSOCIATED(deriv_att))
THEN
3753 DO k = bo(1, 3), bo(2, 3)
3754 DO j = bo(1, 2), bo(2, 2)
3755 DO i = bo(1, 1), bo(2, 1)
3756 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3757 fac*deriv_data(i, j, k)*tau1b(i, j, k)
3763 IF (
ASSOCIATED(deriv_att))
THEN
3767 DO k = bo(1, 3), bo(2, 3)
3768 DO j = bo(1, 2), bo(2, 2)
3769 DO i = bo(1, 1), bo(2, 1)
3770 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3771 fac*deriv_data(i, j, k)*laplace1b(i, j, k)
3778 IF (
ASSOCIATED(deriv_att))
THEN
3784 v_drho(1)%array(:, :, :) = v_drho(1)%array(:, :, :) + &
3785 deriv_data(:, :, :)*dr1dr(:, :, :)/max(gradient_cut, norm_drho(:, :, :))**2
3790 IF (
ASSOCIATED(deriv_att))
THEN
3794 DO k = bo(1, 3), bo(2, 3)
3795 DO j = bo(1, 2), bo(2, 2)
3796 DO i = bo(1, 1), bo(2, 1)
3797 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3798 deriv_data(i, j, k)*rho1a(i, j, k)
3804 IF (
ASSOCIATED(deriv_att))
THEN
3808 DO k = bo(1, 3), bo(2, 3)
3809 DO j = bo(1, 2), bo(2, 2)
3810 DO i = bo(1, 1), bo(2, 1)
3811 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3812 deriv_data(i, j, k)*dr1dr(i, j, k)
3818 IF (
ASSOCIATED(deriv_att))
THEN
3822 DO k = bo(1, 3), bo(2, 3)
3823 DO j = bo(1, 2), bo(2, 2)
3824 DO i = bo(1, 1), bo(2, 1)
3825 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3826 deriv_data(i, j, k)*dra1dra(i, j, k)
3832 IF (
ASSOCIATED(deriv_att))
THEN
3836 DO k = bo(1, 3), bo(2, 3)
3837 DO j = bo(1, 2), bo(2, 2)
3838 DO i = bo(1, 1), bo(2, 1)
3839 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3840 deriv_data(i, j, k)*tau1a(i, j, k)
3846 IF (
ASSOCIATED(deriv_att))
THEN
3850 DO k = bo(1, 3), bo(2, 3)
3851 DO j = bo(1, 2), bo(2, 2)
3852 DO i = bo(1, 1), bo(2, 1)
3853 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3854 deriv_data(i, j, k)*laplace1a(i, j, k)
3860 IF (
ASSOCIATED(deriv_att))
THEN
3864 DO k = bo(1, 3), bo(2, 3)
3865 DO j = bo(1, 2), bo(2, 2)
3866 DO i = bo(1, 1), bo(2, 1)
3867 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3868 fac*deriv_data(i, j, k)*rho1b(i, j, k)
3874 IF (
ASSOCIATED(deriv_att))
THEN
3878 DO k = bo(1, 3), bo(2, 3)
3879 DO j = bo(1, 2), bo(2, 2)
3880 DO i = bo(1, 1), bo(2, 1)
3881 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3882 fac*deriv_data(i, j, k)*drb1drb(i, j, k)
3888 IF (
ASSOCIATED(deriv_att))
THEN
3892 DO k = bo(1, 3), bo(2, 3)
3893 DO j = bo(1, 2), bo(2, 2)
3894 DO i = bo(1, 1), bo(2, 1)
3895 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3896 fac*deriv_data(i, j, k)*tau1b(i, j, k)
3902 IF (
ASSOCIATED(deriv_att))
THEN
3906 DO k = bo(1, 3), bo(2, 3)
3907 DO j = bo(1, 2), bo(2, 2)
3908 DO i = bo(1, 1), bo(2, 1)
3909 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3910 fac*deriv_data(i, j, k)*laplace1b(i, j, k)
3917 IF (
ASSOCIATED(deriv_att))
THEN
3923 v_drhoa(1)%array(:, :, :) = v_drhoa(1)%array(:, :, :) + &
3924 deriv_data(:, :, :)*dra1dra(:, :, :)/max(gradient_cut, norm_drhoa(:, :, :))**2
3929 IF (
ASSOCIATED(deriv_att))
THEN
3933 DO k = bo(1, 3), bo(2, 3)
3934 DO j = bo(1, 2), bo(2, 2)
3935 DO i = bo(1, 1), bo(2, 1)
3936 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3937 deriv_data(i, j, k)*rho1a(i, j, k)
3943 IF (
ASSOCIATED(deriv_att))
THEN
3947 DO k = bo(1, 3), bo(2, 3)
3948 DO j = bo(1, 2), bo(2, 2)
3949 DO i = bo(1, 1), bo(2, 1)
3950 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3951 deriv_data(i, j, k)*dr1dr(i, j, k)
3957 IF (
ASSOCIATED(deriv_att))
THEN
3961 DO k = bo(1, 3), bo(2, 3)
3962 DO j = bo(1, 2), bo(2, 2)
3963 DO i = bo(1, 1), bo(2, 1)
3964 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3965 deriv_data(i, j, k)*dra1dra(i, j, k)
3971 IF (
ASSOCIATED(deriv_att))
THEN
3975 DO k = bo(1, 3), bo(2, 3)
3976 DO j = bo(1, 2), bo(2, 2)
3977 DO i = bo(1, 1), bo(2, 1)
3978 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3979 deriv_data(i, j, k)*tau1a(i, j, k)
3985 IF (
ASSOCIATED(deriv_att))
THEN
3989 DO k = bo(1, 3), bo(2, 3)
3990 DO j = bo(1, 2), bo(2, 2)
3991 DO i = bo(1, 1), bo(2, 1)
3992 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3993 deriv_data(i, j, k)*laplace1a(i, j, k)
3999 IF (
ASSOCIATED(deriv_att))
THEN
4003 DO k = bo(1, 3), bo(2, 3)
4004 DO j = bo(1, 2), bo(2, 2)
4005 DO i = bo(1, 1), bo(2, 1)
4006 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4007 fac*deriv_data(i, j, k)*rho1b(i, j, k)
4013 IF (
ASSOCIATED(deriv_att))
THEN
4017 DO k = bo(1, 3), bo(2, 3)
4018 DO j = bo(1, 2), bo(2, 2)
4019 DO i = bo(1, 1), bo(2, 1)
4020 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4021 fac*deriv_data(i, j, k)*drb1drb(i, j, k)
4027 IF (
ASSOCIATED(deriv_att))
THEN
4031 DO k = bo(1, 3), bo(2, 3)
4032 DO j = bo(1, 2), bo(2, 2)
4033 DO i = bo(1, 1), bo(2, 1)
4034 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4035 fac*deriv_data(i, j, k)*tau1b(i, j, k)
4041 IF (
ASSOCIATED(deriv_att))
THEN
4045 DO k = bo(1, 3), bo(2, 3)
4046 DO j = bo(1, 2), bo(2, 2)
4047 DO i = bo(1, 1), bo(2, 1)
4048 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4049 fac*deriv_data(i, j, k)*laplace1b(i, j, k)
4057 IF (
ASSOCIATED(deriv_att))
THEN
4061 DO k = bo(1, 3), bo(2, 3)
4062 DO j = bo(1, 2), bo(2, 2)
4063 DO i = bo(1, 1), bo(2, 1)
4064 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4065 deriv_data(i, j, k)*rho1a(i, j, k)
4071 IF (
ASSOCIATED(deriv_att))
THEN
4075 DO k = bo(1, 3), bo(2, 3)
4076 DO j = bo(1, 2), bo(2, 2)
4077 DO i = bo(1, 1), bo(2, 1)
4078 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4079 deriv_data(i, j, k)*dr1dr(i, j, k)
4085 IF (
ASSOCIATED(deriv_att))
THEN
4089 DO k = bo(1, 3), bo(2, 3)
4090 DO j = bo(1, 2), bo(2, 2)
4091 DO i = bo(1, 1), bo(2, 1)
4092 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4093 deriv_data(i, j, k)*dra1dra(i, j, k)
4099 IF (
ASSOCIATED(deriv_att))
THEN
4103 DO k = bo(1, 3), bo(2, 3)
4104 DO j = bo(1, 2), bo(2, 2)
4105 DO i = bo(1, 1), bo(2, 1)
4106 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4107 deriv_data(i, j, k)*tau1a(i, j, k)
4113 IF (
ASSOCIATED(deriv_att))
THEN
4117 DO k = bo(1, 3), bo(2, 3)
4118 DO j = bo(1, 2), bo(2, 2)
4119 DO i = bo(1, 1), bo(2, 1)
4120 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4121 deriv_data(i, j, k)*laplace1a(i, j, k)
4127 IF (
ASSOCIATED(deriv_att))
THEN
4131 DO k = bo(1, 3), bo(2, 3)
4132 DO j = bo(1, 2), bo(2, 2)
4133 DO i = bo(1, 1), bo(2, 1)
4134 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4135 fac*deriv_data(i, j, k)*rho1b(i, j, k)
4141 IF (
ASSOCIATED(deriv_att))
THEN
4145 DO k = bo(1, 3), bo(2, 3)
4146 DO j = bo(1, 2), bo(2, 2)
4147 DO i = bo(1, 1), bo(2, 1)
4148 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4149 fac*deriv_data(i, j, k)*drb1drb(i, j, k)
4155 IF (
ASSOCIATED(deriv_att))
THEN
4159 DO k = bo(1, 3), bo(2, 3)
4160 DO j = bo(1, 2), bo(2, 2)
4161 DO i = bo(1, 1), bo(2, 1)
4162 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4163 fac*deriv_data(i, j, k)*tau1b(i, j, k)
4169 IF (
ASSOCIATED(deriv_att))
THEN
4173 DO k = bo(1, 3), bo(2, 3)
4174 DO j = bo(1, 2), bo(2, 2)
4175 DO i = bo(1, 1), bo(2, 1)
4176 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4177 fac*deriv_data(i, j, k)*laplace1b(i, j, k)
4188 IF (gradient_f)
THEN
4189 IF (.NOT. do_spinflip)
THEN
4191 IF (my_compute_virial)
THEN
4192 CALL virial_drho_drho(virial_pw, drhoa, v_drhoa(1), virial_xc)
4193 CALL virial_drho_drho(virial_pw, drhob, v_drhob(2), virial_xc)
4196 virial_pw%array(:, :, :) = drho(idir)%array(:, :, :)*(v_drho(1)%array(:, :, :) + v_drho(2)%array(:, :, :))
4200 drho(jdir)%array(:, :, :))
4201 virial_xc(jdir, idir) = virial_xc(jdir, idir) + tmp
4202 virial_xc(idir, jdir) = virial_xc(jdir, idir)
4212 DO ir = bo(1, 2), bo(2, 2)
4213 DO ia = bo(1, 1), bo(2, 1)
4215 DO ispin = 1, nspins
4216 vxg(idir, ia, ir, ispin) = &
4217 -(v_drhoa(ispin)%array(ia, ir, 1)*drhoa(idir)%array(ia, ir, 1) + &
4218 v_drhob(ispin)%array(ia, ir, 1)*drhob(idir)%array(ia, ir, 1) + &
4219 v_drho(ispin)%array(ia, ir, 1)*drho(idir)%array(ia, ir, 1))
4221 IF (
ASSOCIATED(e_drhoa))
THEN
4222 vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
4223 e_drhoa(ia, ir, 1)*drho1a(idir)%array(ia, ir, 1)
4225 IF (nspins /= 1 .AND.
ASSOCIATED(e_drhob))
THEN
4226 vxg(idir, ia, ir, 2) = vxg(idir, ia, ir, 2) + &
4227 e_drhob(ia, ir, 1)*drho1b(idir)%array(ia, ir, 1)
4229 IF (
ASSOCIATED(e_drho))
THEN
4230 IF (nspins /= 1)
THEN
4231 vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
4232 e_drho(ia, ir, 1)*drho1(idir)%array(ia, ir, 1)
4233 vxg(idir, ia, ir, 2) = vxg(idir, ia, ir, 2) + &
4234 e_drho(ia, ir, 1)*drho1(idir)%array(ia, ir, 1)
4236 vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
4237 e_drho(ia, ir, 1)*(drho1a(idir)%array(ia, ir, 1) + &
4238 fac*drho1b(idir)%array(ia, ir, 1))
4250 DO ispin = 1, nspins
4253 v_drho_r(idir, ispin)%array(:, :, :) = &
4254 v_drhoa(ispin)%array(:, :, :)*drhoa(idir)%array(:, :, :) + &
4255 v_drhob(ispin)%array(:, :, :)*drhob(idir)%array(:, :, :) + &
4256 v_drho(ispin)%array(:, :, :)*drho(idir)%array(:, :, :)
4259 IF (
ASSOCIATED(e_drhoa))
THEN
4262 v_drho_r(idir, 1)%array(:, :, :) = v_drho_r(idir, 1)%array(:, :, :) - &
4263 e_drhoa(:, :, :)*drho1a(idir)%array(:, :, :)
4266 IF (nspins /= 1 .AND.
ASSOCIATED(e_drhob))
THEN
4269 v_drho_r(idir, 2)%array(:, :, :) = v_drho_r(idir, 2)%array(:, :, :) - &
4270 e_drhob(:, :, :)*drho1b(idir)%array(:, :, :)
4273 IF (
ASSOCIATED(e_drho))
THEN
4276 DO k = bo(1, 3), bo(2, 3)
4277 DO j = bo(1, 2), bo(2, 2)
4278 DO i = bo(1, 1), bo(2, 1)
4279 IF (nspins /= 1)
THEN
4280 v_drho_r(idir, 1)%array(i, j, k) = v_drho_r(idir, 1)%array(i, j, k) - &
4281 e_drho(i, j, k)*drho1(idir)%array(i, j, k)
4282 v_drho_r(idir, 2)%array(i, j, k) = v_drho_r(idir, 2)%array(i, j, k) - &
4283 e_drho(i, j, k)*drho1(idir)%array(i, j, k)
4285 v_drho_r(idir, 1)%array(i, j, k) = v_drho_r(idir, 1)%array(i, j, k) - &
4286 e_drho(i, j, k)*(drho1a(idir)%array(i, j, k) + &
4287 fac*drho1b(idir)%array(i, j, k))
4297 DO ispin = 1, nspins
4298 CALL xc_pw_divergence(xc_deriv_method_id, v_drho_r(:, ispin), tmp_g, vxc_g, v_xc(ispin))
4306 DEALLOCATE (drho(idir)%array)
4307 DEALLOCATE (drho1(idir)%array)
4310 DO ispin = 1, nspins
4311 CALL deallocate_pw(v_drhoa(ispin), pw_pool)
4312 CALL deallocate_pw(v_drhob(ispin), pw_pool)
4315 DEALLOCATE (v_drhoa, v_drhob)
4319 IF (laplace_f .AND. my_compute_virial)
THEN
4320 virial_pw%array(:, :, :) = -rhoa(:, :, :)
4321 CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(1)%array)
4322 virial_pw%array(:, :, :) = -rhob(:, :, :)
4323 CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(2)%array)
4334 IF (gradient_f)
THEN
4337 CALL prepare_dr1dr(dr1dr, drho, drho1)
4343 ALLOCATE (v_laplace(nspins))
4344 DO ispin = 1, nspins
4345 CALL allocate_pw(v_laplace(ispin), pw_pool, bo)
4356 IF (
ASSOCIATED(deriv_att))
THEN
4360 DO k = bo(1, 3), bo(2, 3)
4361 DO j = bo(1, 2), bo(2, 2)
4362 DO i = bo(1, 1), bo(2, 1)
4363 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4364 deriv_data(i, j, k)*rho1(i, j, k)
4370 IF (
ASSOCIATED(deriv_att))
THEN
4374 DO k = bo(1, 3), bo(2, 3)
4375 DO j = bo(1, 2), bo(2, 2)
4376 DO i = bo(1, 1), bo(2, 1)
4377 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4378 deriv_data(i, j, k)*dr1dr(i, j, k)
4384 IF (
ASSOCIATED(deriv_att))
THEN
4388 DO k = bo(1, 3), bo(2, 3)
4389 DO j = bo(1, 2), bo(2, 2)
4390 DO i = bo(1, 1), bo(2, 1)
4391 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4392 deriv_data(i, j, k)*tau1(i, j, k)
4398 IF (
ASSOCIATED(deriv_att))
THEN
4402 DO k = bo(1, 3), bo(2, 3)
4403 DO j = bo(1, 2), bo(2, 2)
4404 DO i = bo(1, 1), bo(2, 1)
4405 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4406 deriv_data(i, j, k)*laplace1(i, j, k)
4414 IF (
ASSOCIATED(deriv_att))
THEN
4418 DO k = bo(1, 3), bo(2, 3)
4419 DO j = bo(1, 2), bo(2, 2)
4420 DO i = bo(1, 1), bo(2, 1)
4421 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
4422 deriv_data(i, j, k)*rho1(i, j, k)
4428 IF (
ASSOCIATED(deriv_att))
THEN
4432 DO k = bo(1, 3), bo(2, 3)
4433 DO j = bo(1, 2), bo(2, 2)
4434 DO i = bo(1, 1), bo(2, 1)
4435 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
4436 deriv_data(i, j, k)*dr1dr(i, j, k)
4442 IF (
ASSOCIATED(deriv_att))
THEN
4446 DO k = bo(1, 3), bo(2, 3)
4447 DO j = bo(1, 2), bo(2, 2)
4448 DO i = bo(1, 1), bo(2, 1)
4449 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
4450 deriv_data(i, j, k)*tau1(i, j, k)
4456 IF (
ASSOCIATED(deriv_att))
THEN
4460 DO k = bo(1, 3), bo(2, 3)
4461 DO j = bo(1, 2), bo(2, 2)
4462 DO i = bo(1, 1), bo(2, 1)
4463 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
4464 deriv_data(i, j, k)*laplace1(i, j, k)
4471 IF (
ASSOCIATED(deriv_att))
THEN
4475 IF (my_compute_virial)
THEN
4476 CALL virial_drho_drho1(virial_pw, drho, drho1, deriv_data, virial_xc)
4480 v_drho(1)%array(:, :, :) = v_drho(1)%array(:, :, :) + &
4481 deriv_data(:, :, :)*dr1dr(:, :, :)/max(gradient_cut, norm_drho(:, :, :))**2
4486 IF (
ASSOCIATED(deriv_att))
THEN
4490 DO k = bo(1, 3), bo(2, 3)
4491 DO j = bo(1, 2), bo(2, 2)
4492 DO i = bo(1, 1), bo(2, 1)
4493 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4494 deriv_data(i, j, k)*rho1(i, j, k)
4500 IF (
ASSOCIATED(deriv_att))
THEN
4504 DO k = bo(1, 3), bo(2, 3)
4505 DO j = bo(1, 2), bo(2, 2)
4506 DO i = bo(1, 1), bo(2, 1)
4507 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4508 deriv_data(i, j, k)*dr1dr(i, j, k)
4514 IF (
ASSOCIATED(deriv_att))
THEN
4518 DO k = bo(1, 3), bo(2, 3)
4519 DO j = bo(1, 2), bo(2, 2)
4520 DO i = bo(1, 1), bo(2, 1)
4521 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4522 deriv_data(i, j, k)*tau1(i, j, k)
4528 IF (
ASSOCIATED(deriv_att))
THEN
4532 DO k = bo(1, 3), bo(2, 3)
4533 DO j = bo(1, 2), bo(2, 2)
4534 DO i = bo(1, 1), bo(2, 1)
4535 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4536 deriv_data(i, j, k)*laplace1(i, j, k)
4544 IF (
ASSOCIATED(deriv_att))
THEN
4548 DO k = bo(1, 3), bo(2, 3)
4549 DO j = bo(1, 2), bo(2, 2)
4550 DO i = bo(1, 1), bo(2, 1)
4551 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
4552 deriv_data(i, j, k)*rho1(i, j, k)
4558 IF (
ASSOCIATED(deriv_att))
THEN
4562 DO k = bo(1, 3), bo(2, 3)
4563 DO j = bo(1, 2), bo(2, 2)
4564 DO i = bo(1, 1), bo(2, 1)
4565 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
4566 deriv_data(i, j, k)*dr1dr(i, j, k)
4572 IF (
ASSOCIATED(deriv_att))
THEN
4576 DO k = bo(1, 3), bo(2, 3)
4577 DO j = bo(1, 2), bo(2, 2)
4578 DO i = bo(1, 1), bo(2, 1)
4579 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
4580 deriv_data(i, j, k)*tau1(i, j, k)
4586 IF (
ASSOCIATED(deriv_att))
THEN
4590 DO k = bo(1, 3), bo(2, 3)
4591 DO j = bo(1, 2), bo(2, 2)
4592 DO i = bo(1, 1), bo(2, 1)
4593 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
4594 deriv_data(i, j, k)*laplace1(i, j, k)
4601 IF (my_compute_virial)
THEN
4603 IF (
ASSOCIATED(deriv_att))
THEN
4606 virial_pw%array(:, :, :) = -rho1(:, :, :)
4607 CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
4612 IF (gradient_f)
THEN
4614 IF (my_compute_virial)
THEN
4615 CALL virial_drho_drho(virial_pw, drho, v_drho(1), virial_xc)
4625 DO ia = bo(1, 1), bo(2, 1)
4626 DO ir = bo(1, 2), bo(2, 2)
4627 vxg(idir, ia, ir, 1) = -drho(idir)%array(ia, ir, 1)*v_drho(1)%array(ia, ir, 1)
4628 IF (
ASSOCIATED(e_drho))
THEN
4629 vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + factor2*drho1(idir)%array(ia, ir, 1)*e_drho(ia, ir, 1)
4640 v_drho_r(idir, 1)%array(:, :, :) = drho(idir)%array(:, :, :)*v_drho(1)%array(:, :, :) - &
4641 drho1(idir)%array(:, :, :)*e_drho(:, :, :)
4645 CALL xc_pw_divergence(xc_deriv_method_id, v_drho_r(:, 1), tmp_g, vxc_g, v_xc(1))
4650 IF (laplace_f .AND. my_compute_virial)
THEN
4651 virial_pw%array(:, :, :) = -rho(:, :, :)
4652 CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(1)%array)
4658 DO ispin = 1, nspins
4659 CALL xc_pw_laplace(v_laplace(ispin), pw_pool, xc_deriv_method_id)
4660 CALL pw_axpy(v_laplace(ispin), v_xc(ispin))
4664 IF (gradient_f)
THEN
4666 DO ispin = 1, nspins
4667 CALL deallocate_pw(v_drho(ispin), pw_pool)
4669 CALL deallocate_pw(v_drho_r(idir, ispin), pw_pool)
4672 DEALLOCATE (v_drho, v_drho_r)
4677 DO ispin = 1, nspins
4678 CALL deallocate_pw(v_laplace(ispin), pw_pool)
4680 DEALLOCATE (v_laplace)
4683 IF (
ASSOCIATED(tmp_g%pw_grid) .AND.
ASSOCIATED(pw_pool))
THEN
4684 CALL pw_pool%give_back_pw(tmp_g)
4687 IF (
ASSOCIATED(vxc_g%pw_grid) .AND.
ASSOCIATED(pw_pool))
THEN
4688 CALL pw_pool%give_back_pw(vxc_g)
4691 IF (my_compute_virial .AND. (gradient_f .OR. laplace_f))
THEN
4692 CALL deallocate_pw(virial_pw, pw_pool)
4695 CALL timestop(handle)
4731 LOGICAL,
INTENT(in),
OPTIONAL :: spinflip
4733 CHARACTER(len=*),
PARAMETER :: routinen =
'xc_calc_3rd_deriv_analytical'
4735 INTEGER :: handle, i, idir, ispin, j, &
4736 k, nspins, xc_deriv_method_id
4737 INTEGER,
DIMENSION(2, 3) :: bo
4738 LOGICAL :: lsd, do_spinflip, alda0, &
4739 rho_f, gradient_f, tau_f, laplace_f
4740 REAL(kind=
dp) :: s, s_thresh, s_thresh2, gradient_cut
4741 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: dr1dr, dra1dra, drb1drb
4742 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: deriv_data, deriv_data2, e_drhoa, e_drhob, &
4743 e_drho, norm_drho, norm_drhoa, &
4744 norm_drhob, rho1a, rho1b, &
4746 TYPE(
cp_3d_r_cp_type),
DIMENSION(3) :: drho, drho1, drho1a, drho1b, drhoa, drhob
4747 TYPE(
pw_r3d_rs_type),
DIMENSION(:),
ALLOCATABLE :: v_drhoa, v_drhob, v_drho
4752 CALL timeset(routinen, handle)
4754 NULLIFY (e_drhoa, e_drhob, e_drho)
4756 cpassert(
ASSOCIATED(v_xc))
4757 cpassert(
ASSOCIATED(xc_section))
4761 i_val=xc_deriv_method_id)
4764 lsd =
ASSOCIATED(rho_set%rhoa)
4766 do_spinflip = .false.
4767 IF (
PRESENT(spinflip)) do_spinflip = spinflip
4769 bo = rho_set%local_bounds
4771 CALL check_for_derivatives(deriv_set, lsd, rho_f, gradient_f, tau_f, laplace_f)
4781 DO ispin = 1, nspins
4783 v_xc(ispin)%array = 0.0_dp
4787 IF (gradient_f)
THEN
4788 ALLOCATE (v_drho_r(3, nspins), v_drho(nspins))
4789 DO ispin = 1, nspins
4791 CALL allocate_pw(v_drho_r(idir, ispin), pw_pool, bo)
4793 CALL allocate_pw(v_drho(ispin), pw_pool, bo)
4797 IF (
ASSOCIATED(pw_pool))
THEN
4798 CALL pw_pool%create_pw(tmp_g)
4799 CALL pw_pool%create_pw(vxc_g)
4802 cpabort(
"XC_DERIV method is not implemented in GAPW")
4810 cpassert(
ASSOCIATED(v_xc_tau))
4811 DO ispin = 1, nspins
4812 v_xc_tau(ispin)%array = 0.0_dp
4822 IF (do_spinflip)
THEN
4829 IF (gradient_f)
THEN
4831 norm_drho=norm_drho, norm_drhoa=norm_drhoa, norm_drhob=norm_drhob)
4832 IF (do_spinflip)
THEN
4834 CALL calc_drho_from_a(drho1, drho1a)
4837 CALL calc_drho_from_ab(drho1, drho1a, drho1b)
4840 CALL calc_drho_from_ab(drho, drhoa, drhob)
4842 CALL prepare_dr1dr(dra1dra, drhoa, drho1a)
4843 IF (do_spinflip)
THEN
4844 CALL prepare_dr1dr(drb1drb, drhob, drho1a)
4845 CALL prepare_dr1dr(dr1dr, drho, drho1a)
4846 ELSE IF (nspins /= 1)
THEN
4847 CALL prepare_dr1dr(drb1drb, drhob, drho1b)
4848 CALL prepare_dr1dr(dr1dr, drho, drho1)
4850 cpabort(
"Exchange-correlation's third derivative for closed-shell not yet implemented")
4854 ALLOCATE (v_drhoa(nspins), v_drhob(nspins))
4855 DO ispin = 1, nspins
4856 CALL allocate_pw(v_drhoa(ispin), pw_pool, bo)
4857 CALL allocate_pw(v_drhob(ispin), pw_pool, bo)
4863 cpabort(
"Exchange-correlation's laplace analytic third derivative not implemented")
4867 cpabort(
"Exchange-correlation's mGGA analytic third derivative not implemented")
4870 IF (nspins /= 1)
THEN
4872 IF (.NOT. do_spinflip)
THEN
4874 cpabort(
"Exchange-correlation's analytic third derivative not implemented")
4886 IF (
ASSOCIATED(deriv_att))
THEN
4889 IF (
ASSOCIATED(deriv_att))
THEN
4893 DO k = bo(1, 3), bo(2, 3)
4894 DO j = bo(1, 2), bo(2, 2)
4895 DO i = bo(1, 1), bo(2, 1)
4896 s = rhoa(i, j, k) - rhob(i, j, k)
4897 s = -sign(max(s**2, s_thresh2), s)
4898 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4899 (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)**2/s
4911 IF (.NOT. alda0)
THEN
4913 IF (
ASSOCIATED(deriv_att))
THEN
4916 IF (
ASSOCIATED(deriv_att))
THEN
4920 DO k = bo(1, 3), bo(2, 3)
4921 DO j = bo(1, 2), bo(2, 2)
4922 DO i = bo(1, 1), bo(2, 1)
4923 s = rhoa(i, j, k) - rhob(i, j, k)
4924 s = -sign(max(s**2, s_thresh2), s)
4925 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4926 (deriv_data(i, j, k)*dra1dra(i, j, k) - deriv_data2(i, j, k)*drb1drb(i, j, k)) &
4938 DO k = bo(1, 3), bo(2, 3)
4939 DO j = bo(1, 2), bo(2, 2)
4940 DO i = bo(1, 1), bo(2, 1)
4941 v_xc(2)%array(i, j, k) = -v_xc(1)%array(i, j, k)
4954 IF (
ASSOCIATED(deriv_att))
THEN
4957 IF (
ASSOCIATED(deriv_att))
THEN
4961 DO k = bo(1, 3), bo(2, 3)
4962 DO j = bo(1, 2), bo(2, 2)
4963 DO i = bo(1, 1), bo(2, 1)
4964 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
4965 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4966 rho1a(i, j, k)**2*(deriv_data(i, j, k) - deriv_data2(i, j, k))/s
4978 IF (
ASSOCIATED(deriv_att))
THEN
4981 IF (
ASSOCIATED(deriv_att))
THEN
4985 DO k = bo(1, 3), bo(2, 3)
4986 DO j = bo(1, 2), bo(2, 2)
4987 DO i = bo(1, 1), bo(2, 1)
4988 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
4989 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
4990 rho1a(i, j, k)**2*(deriv_data(i, j, k) - deriv_data2(i, j, k))/s
4998 IF (.NOT. alda0)
THEN
5003 IF (
ASSOCIATED(deriv_att))
THEN
5006 IF (
ASSOCIATED(deriv_att))
THEN
5010 DO k = bo(1, 3), bo(2, 3)
5011 DO j = bo(1, 2), bo(2, 2)
5012 DO i = bo(1, 1), bo(2, 1)
5013 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
5014 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
5015 (deriv_data(i, j, k)*dra1dra(i, j, k) - deriv_data2(i, j, k)*drb1drb(i, j, k))* &
5028 IF (
ASSOCIATED(deriv_att))
THEN
5031 IF (
ASSOCIATED(deriv_att))
THEN
5035 DO k = bo(1, 3), bo(2, 3)
5036 DO j = bo(1, 2), bo(2, 2)
5037 DO i = bo(1, 1), bo(2, 1)
5038 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
5039 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
5040 (deriv_data(i, j, k)*dra1dra(i, j, k) - deriv_data2(i, j, k)*drb1drb(i, j, k))* &
5056 IF (
ASSOCIATED(deriv_att))
THEN
5059 IF (
ASSOCIATED(deriv_att))
THEN
5063 DO k = bo(1, 3), bo(2, 3)
5064 DO j = bo(1, 2), bo(2, 2)
5065 DO i = bo(1, 1), bo(2, 1)
5066 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
5067 (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)
5079 IF (
ASSOCIATED(deriv_att))
THEN
5082 IF (
ASSOCIATED(deriv_att))
THEN
5086 DO k = bo(1, 3), bo(2, 3)
5087 DO j = bo(1, 2), bo(2, 2)
5088 DO i = bo(1, 1), bo(2, 1)
5089 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) + &
5090 (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)
5106 IF (
ASSOCIATED(deriv_att))
THEN
5109 IF (
ASSOCIATED(deriv_att))
THEN
5113 DO k = bo(1, 3), bo(2, 3)
5114 DO j = bo(1, 2), bo(2, 2)
5115 DO i = bo(1, 1), bo(2, 1)
5116 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
5117 (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)
5118 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
5119 (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)
5128 IF (
ASSOCIATED(deriv_att))
THEN
5134 IF (
ASSOCIATED(deriv_att))
THEN
5138 DO k = bo(1, 3), bo(2, 3)
5139 DO j = bo(1, 2), bo(2, 2)
5140 DO i = bo(1, 1), bo(2, 1)
5141 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
5142 deriv_data(i, j, k)*dra1dra(i, j, k) + deriv_data2(i, j, k)*drb1drb(i, j, k)
5152 IF (
ASSOCIATED(deriv_att))
THEN
5156 DO k = bo(1, 3), bo(2, 3)
5157 DO j = bo(1, 2), bo(2, 2)
5158 DO i = bo(1, 1), bo(2, 1)
5159 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
5160 deriv_data(i, j, k)*dra1dra(i, j, k) + deriv_data2(i, j, k)*drb1drb(i, j, k)
5175 IF (
ASSOCIATED(deriv_att))
THEN
5178 IF (
ASSOCIATED(deriv_att))
THEN
5182 DO k = bo(1, 3), bo(2, 3)
5183 DO j = bo(1, 2), bo(2, 2)
5184 DO i = bo(1, 1), bo(2, 1)
5185 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
5186 deriv_data(i, j, k)*dra1dra(i, j, k) + &
5187 deriv_data2(i, j, k)*drb1drb(i, j, k)
5188 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
5189 deriv_data(i, j, k)*dra1dra(i, j, k) + &
5190 deriv_data2(i, j, k)*drb1drb(i, j, k)
5205 IF (
ASSOCIATED(deriv_att))
THEN
5210 v_drhoa(1)%array(:, :, :) = v_drhoa(1)%array(:, :, :) + &
5211 deriv_data(:, :, :)*dra1dra(:, :, :)/max(gradient_cut, norm_drhoa(:, :, :))**2
5219 IF (
ASSOCIATED(deriv_att))
THEN
5224 v_drhob(2)%array(:, :, :) = v_drhob(2)%array(:, :, :) - &
5225 deriv_data(:, :, :)*drb1drb(:, :, :)/max(gradient_cut, norm_drhob(:, :, :))**2
5234 cpabort(
"Exchange-correlation's analytic third derivative not implemented")
5238 IF (gradient_f)
THEN
5239 IF (.NOT. alda0)
THEN
5253 IF (do_spinflip)
THEN
5256 DO k = bo(1, 3), bo(2, 3)
5257 DO j = bo(1, 2), bo(2, 2)
5258 DO i = bo(1, 1), bo(2, 1)
5259 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
5261 v_drho_r(idir, ispin)%array(i, j, k) = v_drho_r(idir, ispin)%array(i, j, k) + &
5262 v_drhoa(ispin)%array(i, j, k)*drhoa(idir)%array(i, j, k)*rho1a(i, j, k)/s + &
5263 v_drhob(ispin)%array(i, j, k)*drhob(idir)%array(i, j, k)*rho1a(i, j, k)/s + &
5264 v_drho(ispin)%array(i, j, k)*drho(idir)%array(i, j, k)*rho1a(i, j, k)/s
5276 IF (
ASSOCIATED(e_drhoa))
THEN
5279 DO k = bo(1, 3), bo(2, 3)
5280 DO j = bo(1, 2), bo(2, 2)
5281 DO i = bo(1, 1), bo(2, 1)
5282 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
5283 v_drho_r(idir, 1)%array(i, j, k) = v_drho_r(idir, 1)%array(i, j, k) - &
5284 e_drhoa(i, j, k)*drho1a(idir)%array(i, j, k)*rho1a(i, j, k)/s
5294 IF (
ASSOCIATED(e_drhob))
THEN
5297 DO k = bo(1, 3), bo(2, 3)
5298 DO j = bo(1, 2), bo(2, 2)
5299 DO i = bo(1, 1), bo(2, 1)
5300 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
5301 v_drho_r(idir, 2)%array(i, j, k) = v_drho_r(idir, 2)%array(i, j, k) + &
5302 e_drhob(i, j, k)*drho1a(idir)%array(i, j, k)*rho1a(i, j, k)/s
5311 DO ispin = 1, nspins
5312 CALL xc_pw_divergence(xc_deriv_method_id, v_drho_r(:, ispin), tmp_g, vxc_g, v_xc(ispin))
5317 DEALLOCATE (drho(idir)%array)
5318 DEALLOCATE (drho1(idir)%array)
5321 DO ispin = 1, nspins
5322 CALL deallocate_pw(v_drhoa(ispin), pw_pool)
5323 CALL deallocate_pw(v_drhob(ispin), pw_pool)
5326 DEALLOCATE (v_drhoa, v_drhob)
5335 cpabort(
"Exchange-correlation's analytic third derivative not implemented")
5339 IF (gradient_f)
THEN
5341 DO ispin = 1, nspins
5342 CALL deallocate_pw(v_drho(ispin), pw_pool)
5344 CALL deallocate_pw(v_drho_r(idir, ispin), pw_pool)
5347 DEALLOCATE (v_drho, v_drho_r)
5351 IF (
ASSOCIATED(tmp_g%pw_grid) .AND.
ASSOCIATED(pw_pool))
THEN
5352 CALL pw_pool%give_back_pw(tmp_g)
5355 IF (
ASSOCIATED(vxc_g%pw_grid) .AND.
ASSOCIATED(pw_pool))
THEN
5356 CALL pw_pool%give_back_pw(vxc_g)
5359 CALL timestop(handle)