1055 pw_pool, weights, xc_section, &
1056 do_triplet, calc_virial, virial_xc, deriv_set)
1058 TYPE(
pw_r3d_rs_type),
DIMENSION(:),
INTENT(IN),
POINTER :: v_xc, v_tau
1060 TYPE(
pw_r3d_rs_type),
DIMENSION(:),
INTENT(IN),
POINTER :: rho1_r, tau1_r
1061 TYPE(
pw_c1d_gs_type),
DIMENSION(:),
INTENT(IN),
POINTER :: rho1_g
1065 LOGICAL,
INTENT(IN) :: do_triplet
1066 LOGICAL,
INTENT(IN),
OPTIONAL :: calc_virial
1067 REAL(kind=
dp),
DIMENSION(3, 3),
INTENT(INOUT), &
1068 OPTIONAL :: virial_xc
1071 CHARACTER(len=*),
PARAMETER :: routinen =
'xc_calc_2nd_deriv_numerical'
1072 REAL(kind=
dp),
DIMENSION(-4:4, 4),
PARAMETER :: &
1073 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, &
1074 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, &
1075 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, &
1076 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])
1078 INTEGER :: handle, idir, ispin, nspins, istep, nsteps
1079 INTEGER,
DIMENSION(2, 3) :: bo
1080 LOGICAL :: gradient_f, lsd, my_calc_virial, tau_f, laplace_f, rho_f
1081 REAL(kind=
dp) :: exc, gradient_cut, h, rweight, step, rho_cutoff
1082 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: dr1dr, dra1dra, drb1drb
1083 REAL(kind=
dp),
DIMENSION(3, 3) :: virial_dummy
1084 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: norm_drho, norm_drho2, norm_drho2a, &
1085 norm_drho2b, norm_drhoa, norm_drhob, &
1086 rho, rho1, rho1a, rho1b, rhoa, rhob, &
1087 tau_a, tau_b, tau, tau1, tau1a, tau1b, laplace, laplace1, &
1088 laplacea, laplaceb, laplace1a, laplace1b, &
1089 laplace2, laplace2a, laplace2b, deriv_data
1090 TYPE(
cp_3d_r_cp_type),
DIMENSION(3) :: drho, drho1, drho1a, drho1b, drhoa, drhob
1095 TYPE(
pw_r3d_rs_type) :: virial_pw, v_laplace, v_laplacea, v_laplaceb
1101 CALL timeset(routinen, handle)
1103 my_calc_virial = .false.
1104 IF (
PRESENT(calc_virial) .AND.
PRESENT(virial_xc)) my_calc_virial = calc_virial
1108 NULLIFY (tau, tau_r, tau_a, tau_b)
1112 IF (nsteps < lbound(rweights, 2) .OR. nspins > ubound(rweights, 2))
THEN
1113 cpabort(
"The number of steps must be a value from 1 to 4.")
1116 IF (nspins == 2)
THEN
1117 NULLIFY (vxc_rho, rho_g, vxc_tau)
1119 DO ispin = 1, nspins
1120 CALL pw_pool%create_pw(rho_r(ispin))
1122 IF (
ASSOCIATED(tau1_r) .AND.
ASSOCIATED(v_tau))
THEN
1124 DO ispin = 1, nspins
1125 CALL pw_pool%create_pw(tau_r(ispin))
1128 CALL xc_rho_set_get(rho_set, can_return_null=.true., rhoa=rhoa, rhob=rhob, tau_a=tau_a, tau_b=tau_b)
1129 DO istep = -nsteps, nsteps
1130 IF (istep == 0) cycle
1131 rweight = rweights(istep, nsteps)/h
1132 step = real(istep,
dp)*h
1133 CALL calc_resp_potential_numer_ab(rho_r, rho_g, rho1_r, rhoa, rhob, vxc_rho, &
1134 tau_r, tau1_r, tau_a, tau_b, vxc_tau, xc_section, &
1135 weights, pw_pool, step)
1136 DO ispin = 1, nspins
1137 CALL pw_axpy(vxc_rho(ispin), v_xc(ispin), rweight)
1138 IF (
ASSOCIATED(vxc_tau) .AND.
ASSOCIATED(v_tau))
THEN
1139 CALL pw_axpy(vxc_tau(ispin), v_tau(ispin), rweight)
1142 DO ispin = 1, nspins
1143 CALL vxc_rho(ispin)%release()
1145 DEALLOCATE (vxc_rho)
1146 IF (
ASSOCIATED(vxc_tau))
THEN
1147 DO ispin = 1, nspins
1148 CALL vxc_tau(ispin)%release()
1150 DEALLOCATE (vxc_tau)
1153 ELSE IF (nspins == 1 .AND. do_triplet)
THEN
1154 NULLIFY (vxc_rho, vxc_tau, rho_g)
1157 CALL pw_pool%create_pw(rho_r(ispin))
1159 IF (
ASSOCIATED(tau1_r) .AND.
ASSOCIATED(v_tau))
THEN
1161 DO ispin = 1, nspins
1162 CALL pw_pool%create_pw(tau_r(ispin))
1165 CALL xc_rho_set_get(rho_set, can_return_null=.true., rhoa=rhoa, rhob=rhob, tau_a=tau_a, tau_b=tau_b)
1166 DO istep = -nsteps, nsteps
1167 IF (istep == 0) cycle
1168 rweight = rweights(istep, nsteps)/h
1169 step = real(istep,
dp)*h
1173 rho_r(1)%array(:, :, :) = rhoa(:, :, :) + step*rho1_r(1)%array(:, :, :)
1176 rho_r(2)%array(:, :, :) = rhob(:, :, :)
1178 IF (
ASSOCIATED(tau1_r))
THEN
1180 tau_r(1)%array(:, :, :) = tau_a(:, :, :) + step*tau1_r(1)%array(:, :, :)
1183 tau_r(2)%array(:, :, :) = tau_b(:, :, :)
1187 CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
1188 weights, pw_pool, .false., virial_dummy)
1189 CALL pw_axpy(vxc_rho(1), v_xc(1), rweight)
1190 IF (
ASSOCIATED(vxc_tau) .AND.
ASSOCIATED(v_tau))
THEN
1191 CALL pw_axpy(vxc_tau(1), v_tau(1), rweight)
1194 CALL vxc_rho(ispin)%release()
1196 DEALLOCATE (vxc_rho)
1197 IF (
ASSOCIATED(vxc_tau))
THEN
1199 CALL vxc_tau(ispin)%release()
1201 DEALLOCATE (vxc_tau)
1206 rho_r(1)%array(:, :, :) = rhoa(:, :, :)
1209 rho_r(2)%array(:, :, :) = rhob(:, :, :) + step*rho1_r(1)%array(:, :, :)
1211 IF (
ASSOCIATED(tau1_r))
THEN
1213 tau_r(1)%array(:, :, :) = tau_a(:, :, :)
1216 tau_r(2)%array(:, :, :) = tau_b(:, :, :) + step*tau1_r(1)%array(:, :, :)
1220 CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
1221 weights, pw_pool, .false., virial_dummy)
1222 CALL pw_axpy(vxc_rho(1), v_xc(1), rweight)
1223 IF (
ASSOCIATED(vxc_tau) .AND.
ASSOCIATED(v_tau))
THEN
1224 CALL pw_axpy(vxc_tau(1), v_tau(1), rweight)
1227 CALL vxc_rho(ispin)%release()
1229 DEALLOCATE (vxc_rho)
1230 IF (
ASSOCIATED(vxc_tau))
THEN
1232 CALL vxc_tau(ispin)%release()
1234 DEALLOCATE (vxc_tau)
1238 NULLIFY (vxc_rho, rho_r, rho_g, vxc_tau, tau_r, tau)
1240 CALL pw_pool%create_pw(rho_r(1))
1241 IF (
ASSOCIATED(tau1_r) .AND.
ASSOCIATED(v_tau))
THEN
1243 CALL pw_pool%create_pw(tau_r(1))
1245 CALL xc_rho_set_get(rho_set, can_return_null=.true., rho=rho, tau=tau)
1246 DO istep = -nsteps, nsteps
1247 IF (istep == 0) cycle
1248 rweight = rweights(istep, nsteps)/h
1249 step = real(istep,
dp)*h
1252 rho_r(1)%array(:, :, :) = rho(:, :, :) + step*rho1_r(1)%array(:, :, :)
1254 IF (
ASSOCIATED(tau1_r) .AND.
ASSOCIATED(tau) .AND.
ASSOCIATED(tau1_r))
THEN
1256 tau_r(1)%array(:, :, :) = tau(:, :, :) + step*tau1_r(1)%array(:, :, :)
1260 CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
1261 weights, pw_pool, .false., virial_dummy)
1262 CALL pw_axpy(vxc_rho(1), v_xc(1), rweight)
1263 IF (
ASSOCIATED(vxc_tau) .AND.
ASSOCIATED(v_tau))
THEN
1264 CALL pw_axpy(vxc_tau(1), v_tau(1), rweight)
1266 CALL vxc_rho(1)%release()
1267 DEALLOCATE (vxc_rho)
1268 IF (
ASSOCIATED(vxc_tau))
THEN
1269 CALL vxc_tau(1)%release()
1270 DEALLOCATE (vxc_tau)
1275 IF (my_calc_virial)
THEN
1277 IF (nspins == 1 .AND. do_triplet)
THEN
1281 CALL check_for_derivatives(deriv_set, (nspins == 2), rho_f, gradient_f, tau_f, laplace_f)
1287 IF (gradient_f)
THEN
1288 bo = rho_set%local_bounds
1291 CALL allocate_pw(virial_pw, pw_pool, bo)
1298 drho_cutoff=gradient_cut, &
1311 CALL xc_rho_set_get(rho_set, drhoa=drhoa, drhob=drhob, norm_drho=norm_drho, &
1312 norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, tau_a=tau_a, tau_b=tau_b, &
1313 laplace_rhoa=laplacea, laplace_rhob=laplaceb, can_return_null=.true.)
1314 CALL xc_rho_set_get(rho1_set, rhoa=rho1a, rhob=rho1b, drhoa=drho1a, drhob=drho1b, laplace_rhoa=laplace1a, &
1315 laplace_rhob=laplace1b, can_return_null=.true.)
1317 CALL calc_drho_from_ab(drho, drhoa, drhob)
1318 CALL calc_drho_from_ab(drho1, drho1a, drho1b)
1320 CALL xc_rho_set_get(rho_set, drho=drho, norm_drho=norm_drho, tau=tau, laplace_rho=laplace, can_return_null=.true.)
1321 CALL xc_rho_set_get(rho1_set, rho=rho1, drho=drho1, laplace_rho=laplace1, can_return_null=.true.)
1324 CALL prepare_dr1dr(dr1dr, drho, drho1)
1327 CALL prepare_dr1dr(dra1dra, drhoa, drho1a)
1328 CALL prepare_dr1dr(drb1drb, drhob, drho1b)
1330 CALL allocate_pw(v_drho, pw_pool, bo)
1331 CALL allocate_pw(v_drhoa, pw_pool, bo)
1332 CALL allocate_pw(v_drhob, pw_pool, bo)
1334 IF (
ASSOCIATED(norm_drhoa))
CALL apply_drho(deriv_set, [
deriv_norm_drhoa], virial_pw, &
1335 drhoa, drho1a, virial_xc, &
1336 norm_drhoa, gradient_cut, dra1dra, v_drhoa%array)
1337 IF (
ASSOCIATED(norm_drhob))
CALL apply_drho(deriv_set, [
deriv_norm_drhob], virial_pw, &
1338 drhob, drho1b, virial_xc, &
1339 norm_drhob, gradient_cut, drb1drb, v_drhob%array)
1340 IF (
ASSOCIATED(norm_drho))
CALL apply_drho(deriv_set, [
deriv_norm_drho], virial_pw, &
1341 drho, drho1, virial_xc, &
1342 norm_drho, gradient_cut, dr1dr, v_drho%array)
1345 cpassert(
ASSOCIATED(deriv_data))
1346 virial_pw%array(:, :, :) = -rho1a(:, :, :)
1347 CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
1349 CALL allocate_pw(v_laplacea, pw_pool, bo)
1352 cpassert(
ASSOCIATED(deriv_data))
1353 virial_pw%array(:, :, :) = -rho1b(:, :, :)
1354 CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
1356 CALL allocate_pw(v_laplaceb, pw_pool, bo)
1362 CALL allocate_pw(v_drho, pw_pool, bo)
1364 CALL apply_drho(deriv_set, [
deriv_norm_drho], virial_pw, drho, drho1, virial_xc, &
1365 norm_drho, gradient_cut, dr1dr, v_drho%array)
1368 cpassert(
ASSOCIATED(deriv_data))
1369 virial_pw%array(:, :, :) = -rho1(:, :, :)
1370 CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
1372 CALL allocate_pw(v_laplace, pw_pool, bo)
1378 rho_r(1)%array = rhoa
1379 rho_r(2)%array = rhob
1381 rho_r(1)%array = rho
1383 IF (
ASSOCIATED(tau1_r))
THEN
1385 tau_r(1)%array = tau_a
1386 tau_r(2)%array = tau_b
1388 tau_r(1)%array = tau
1399 rho_cutoff=rho_cutoff, &
1410 CALL xc_rho_set_get(rho1_set, rhoa=rho1a, rhob=rho1b, tau_a=tau1a, tau_b=tau1b, &
1411 laplace_rhoa=laplace1a, laplace_rhob=laplace1b, can_return_null=.true.)
1412 CALL xc_rho_set_get(rho2_set, norm_drhoa=norm_drho2a, norm_drhob=norm_drho2b, &
1413 norm_drho=norm_drho2, laplace_rhoa=laplace2a, laplace_rhob=laplace2b, can_return_null=.true.)
1415 DO istep = -nsteps, nsteps
1416 IF (istep == 0) cycle
1417 rweight = rweights(istep, nsteps)/h
1418 step = real(istep,
dp)*h
1419 IF (
ASSOCIATED(norm_drhoa))
THEN
1420 CALL get_derivs_rho(norm_drho2a, norm_drhoa, step, xc_fun_section, lsd, rho2_set, deriv_set1)
1421 CALL update_deriv_rho(deriv_set1, [
deriv_rhoa], bo, &
1422 norm_drhoa, gradient_cut, rweight, rho1a, v_drhoa%array)
1423 CALL update_deriv_rho(deriv_set1, [
deriv_rhob], bo, &
1424 norm_drhoa, gradient_cut, rweight, rho1b, v_drhoa%array)
1426 norm_drhoa, gradient_cut, rweight, dra1dra, v_drhoa%array)
1428 norm_drhoa, gradient_cut, rweight, dra1dra, drb1drb, v_drhoa%array, v_drhob%array)
1430 norm_drhoa, gradient_cut, rweight, dra1dra, dr1dr, v_drhoa%array, v_drho%array)
1432 CALL update_deriv_rho(deriv_set1, [
deriv_tau_a], bo, &
1433 norm_drhoa, gradient_cut, rweight, tau1a, v_drhoa%array)
1434 CALL update_deriv_rho(deriv_set1, [
deriv_tau_b], bo, &
1435 norm_drhoa, gradient_cut, rweight, tau1b, v_drhoa%array)
1439 norm_drhoa, gradient_cut, rweight, laplace1a, v_drhoa%array)
1441 norm_drhoa, gradient_cut, rweight, laplace1b, v_drhoa%array)
1445 IF (
ASSOCIATED(norm_drhob))
THEN
1446 CALL get_derivs_rho(norm_drho2b, norm_drhob, step, xc_fun_section, lsd, rho2_set, deriv_set1)
1447 CALL update_deriv_rho(deriv_set1, [
deriv_rhoa], bo, &
1448 norm_drhob, gradient_cut, rweight, rho1a, v_drhob%array)
1449 CALL update_deriv_rho(deriv_set1, [
deriv_rhob], bo, &
1450 norm_drhob, gradient_cut, rweight, rho1b, v_drhob%array)
1452 norm_drhob, gradient_cut, rweight, drb1drb, v_drhob%array)
1454 norm_drhob, gradient_cut, rweight, drb1drb, dra1dra, v_drhob%array, v_drhoa%array)
1456 norm_drhob, gradient_cut, rweight, drb1drb, dr1dr, v_drhob%array, v_drho%array)
1458 CALL update_deriv_rho(deriv_set1, [
deriv_tau_a], bo, &
1459 norm_drhob, gradient_cut, rweight, tau1a, v_drhob%array)
1460 CALL update_deriv_rho(deriv_set1, [
deriv_tau_b], bo, &
1461 norm_drhob, gradient_cut, rweight, tau1b, v_drhob%array)
1465 norm_drhob, gradient_cut, rweight, laplace1a, v_drhob%array)
1467 norm_drhob, gradient_cut, rweight, laplace1b, v_drhob%array)
1471 IF (
ASSOCIATED(norm_drho))
THEN
1472 CALL get_derivs_rho(norm_drho2, norm_drho, step, xc_fun_section, lsd, rho2_set, deriv_set1)
1473 CALL update_deriv_rho(deriv_set1, [
deriv_rhoa], bo, &
1474 norm_drho, gradient_cut, rweight, rho1a, v_drho%array)
1475 CALL update_deriv_rho(deriv_set1, [
deriv_rhob], bo, &
1476 norm_drho, gradient_cut, rweight, rho1b, v_drho%array)
1478 norm_drho, gradient_cut, rweight, dr1dr, v_drho%array)
1480 norm_drho, gradient_cut, rweight, dr1dr, dra1dra, v_drho%array, v_drhoa%array)
1482 norm_drho, gradient_cut, rweight, dr1dr, drb1drb, v_drho%array, v_drhob%array)
1484 CALL update_deriv_rho(deriv_set1, [
deriv_tau_a], bo, &
1485 norm_drho, gradient_cut, rweight, tau1a, v_drho%array)
1486 CALL update_deriv_rho(deriv_set1, [
deriv_tau_b], bo, &
1487 norm_drho, gradient_cut, rweight, tau1b, v_drho%array)
1491 norm_drho, gradient_cut, rweight, laplace1a, v_drho%array)
1493 norm_drho, gradient_cut, rweight, laplace1b, v_drho%array)
1499 CALL get_derivs_rho(laplace2a, laplacea, step, xc_fun_section, lsd, rho2_set, deriv_set1)
1502 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_rhoa], bo, &
1503 rweight, rho1a, v_laplacea%array)
1504 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_rhob], bo, &
1505 rweight, rho1b, v_laplacea%array)
1506 IF (
ASSOCIATED(norm_drho))
THEN
1507 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_norm_drho], bo, &
1508 rweight, dr1dr, v_laplacea%array)
1510 IF (
ASSOCIATED(norm_drhoa))
THEN
1511 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_norm_drhoa], bo, &
1512 rweight, dra1dra, v_laplacea%array)
1514 IF (
ASSOCIATED(norm_drhob))
THEN
1515 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_norm_drhob], bo, &
1516 rweight, drb1drb, v_laplacea%array)
1519 IF (
ASSOCIATED(tau1a))
THEN
1520 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_tau_a], bo, &
1521 rweight, tau1a, v_laplacea%array)
1523 IF (
ASSOCIATED(tau1b))
THEN
1524 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_tau_b], bo, &
1525 rweight, tau1b, v_laplacea%array)
1529 rweight, laplace1a, v_laplacea%array)
1532 rweight, laplace1b, v_laplacea%array)
1535 CALL get_derivs_rho(laplace2b, laplaceb, step, xc_fun_section, lsd, rho2_set, deriv_set1)
1538 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_rhoa], bo, &
1539 rweight, rho1a, v_laplaceb%array)
1540 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_rhob], bo, &
1541 rweight, rho1b, v_laplaceb%array)
1542 IF (
ASSOCIATED(norm_drho))
THEN
1543 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_norm_drho], bo, &
1544 rweight, dr1dr, v_laplaceb%array)
1546 IF (
ASSOCIATED(norm_drhoa))
THEN
1547 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_norm_drhoa], bo, &
1548 rweight, dra1dra, v_laplaceb%array)
1550 IF (
ASSOCIATED(norm_drhob))
THEN
1551 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_norm_drhob], bo, &
1552 rweight, drb1drb, v_laplaceb%array)
1556 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_tau_a], bo, &
1557 rweight, tau1a, v_laplaceb%array)
1558 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_tau_b], bo, &
1559 rweight, tau1b, v_laplaceb%array)
1563 rweight, laplace1a, v_laplaceb%array)
1566 rweight, laplace1b, v_laplaceb%array)
1570 CALL virial_drho_drho(virial_pw, drhoa, v_drhoa, virial_xc)
1571 CALL virial_drho_drho(virial_pw, drhob, v_drhob, virial_xc)
1572 CALL virial_drho_drho(virial_pw, drho, v_drho, virial_xc)
1574 CALL deallocate_pw(v_drho, pw_pool)
1575 CALL deallocate_pw(v_drhoa, pw_pool)
1576 CALL deallocate_pw(v_drhob, pw_pool)
1579 virial_pw%array(:, :, :) = -rhoa(:, :, :)
1580 CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplacea%array)
1581 CALL deallocate_pw(v_laplacea, pw_pool)
1583 virial_pw%array(:, :, :) = -rhob(:, :, :)
1584 CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplaceb%array)
1585 CALL deallocate_pw(v_laplaceb, pw_pool)
1588 CALL deallocate_pw(virial_pw, pw_pool)
1591 DEALLOCATE (drho(idir)%array)
1592 DEALLOCATE (drho1(idir)%array)
1594 DEALLOCATE (dra1dra, drb1drb)
1597 CALL xc_rho_set_get(rho1_set, rho=rho1, tau=tau1, laplace_rho=laplace1, can_return_null=.true.)
1598 CALL xc_rho_set_get(rho2_set, norm_drho=norm_drho2, laplace_rho=laplace2, can_return_null=.true.)
1600 DO istep = -nsteps, nsteps
1601 IF (istep == 0) cycle
1602 rweight = rweights(istep, nsteps)/h
1603 step = real(istep,
dp)*h
1604 CALL get_derivs_rho(norm_drho2, norm_drho, step, xc_fun_section, lsd, rho2_set, deriv_set1)
1607 CALL update_deriv_rho(deriv_set1, [
deriv_rho], bo, &
1608 norm_drho, gradient_cut, rweight, rho1, v_drho%array)
1610 norm_drho, gradient_cut, rweight, dr1dr, v_drho%array)
1613 CALL update_deriv_rho(deriv_set1, [
deriv_tau], bo, &
1614 norm_drho, gradient_cut, rweight, tau1, v_drho%array)
1618 norm_drho, gradient_cut, rweight, laplace1, v_drho%array)
1620 CALL get_derivs_rho(laplace2, laplace, step, xc_fun_section, lsd, rho2_set, deriv_set1)
1623 CALL update_deriv(deriv_set1, laplace, rho_cutoff, [
deriv_rho], bo, &
1624 rweight, rho1, v_laplace%array)
1625 CALL update_deriv(deriv_set1, laplace, rho_cutoff, [
deriv_norm_drho], bo, &
1626 rweight, dr1dr, v_laplace%array)
1629 CALL update_deriv(deriv_set1, laplace, rho_cutoff, [
deriv_tau], bo, &
1630 rweight, tau1, v_laplace%array)
1634 rweight, laplace1, v_laplace%array)
1639 CALL virial_drho_drho(virial_pw, drho, v_drho, virial_xc)
1641 CALL deallocate_pw(v_drho, pw_pool)
1644 virial_pw%array(:, :, :) = -rho(:, :, :)
1645 CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace%array)
1646 CALL deallocate_pw(v_laplace, pw_pool)
1649 CALL deallocate_pw(virial_pw, pw_pool)
1662 DO ispin = 1,
SIZE(rho_r)
1663 CALL pw_pool%give_back_pw(rho_r(ispin))
1667 IF (
ASSOCIATED(tau_r))
THEN
1668 DO ispin = 1,
SIZE(tau_r)
1669 CALL pw_pool%give_back_pw(tau_r(ispin))
1674 CALL timestop(handle)
2049 LOGICAL,
INTENT(IN),
OPTIONAL :: gapw
2050 REAL(kind=
dp),
DIMENSION(:, :, :, :),
OPTIONAL, &
2052 REAL(kind=
dp),
INTENT(in),
OPTIONAL :: tddfpt_fac
2053 LOGICAL,
INTENT(IN),
OPTIONAL :: compute_virial
2054 REAL(kind=
dp),
DIMENSION(3, 3),
INTENT(INOUT), &
2055 OPTIONAL :: virial_xc
2056 LOGICAL,
INTENT(in),
OPTIONAL :: spinflip
2058 CHARACTER(len=*),
PARAMETER :: routinen =
'xc_calc_2nd_deriv_analytical'
2060 INTEGER :: handle, i, ia, idir, ir, ispin, j, jdir, &
2061 k, nspins, xc_deriv_method_id
2062 INTEGER,
DIMENSION(2, 3) :: bo
2063 LOGICAL :: gradient_f, lsd, my_compute_virial, alda0, &
2064 my_gapw, tau_f, laplace_f, rho_f, do_spinflip
2065 REAL(kind=
dp) ::
fac, gradient_cut, tmp, factor2, s, s_thresh
2066 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: dr1dr, dra1dra, drb1drb
2067 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: deriv_data, deriv_data2, &
2068 e_drhoa, e_drhob, e_drho, norm_drho, norm_drhoa, &
2069 norm_drhob, rho1, rho1a, rho1b, &
2070 tau1, tau1a, tau1b, laplace1, laplace1a, laplace1b, &
2072 TYPE(
cp_3d_r_cp_type),
DIMENSION(3) :: drho, drho1, drho1a, drho1b, drhoa, drhob
2073 TYPE(
pw_r3d_rs_type),
DIMENSION(:),
ALLOCATABLE :: v_drhoa, v_drhob, v_drho, v_laplace
2079 CALL timeset(routinen, handle)
2081 NULLIFY (e_drhoa, e_drhob, e_drho)
2084 IF (
PRESENT(gapw)) my_gapw = gapw
2086 my_compute_virial = .false.
2087 IF (
PRESENT(compute_virial)) my_compute_virial = compute_virial
2089 cpassert(
ASSOCIATED(v_xc))
2090 cpassert(
ASSOCIATED(xc_section))
2092 cpassert(
PRESENT(vxg))
2094 IF (my_compute_virial)
THEN
2095 cpassert(
PRESENT(virial_xc))
2099 i_val=xc_deriv_method_id)
2102 lsd =
ASSOCIATED(rho_set%rhoa)
2105 IF (
PRESENT(tddfpt_fac))
fac = tddfpt_fac
2106 IF (
PRESENT(tddfpt_fac)) factor2 = tddfpt_fac
2107 do_spinflip = .false.
2108 IF (
PRESENT(spinflip)) do_spinflip = spinflip
2110 bo = rho_set%local_bounds
2112 CALL check_for_derivatives(deriv_set, lsd, rho_f, gradient_f, tau_f, laplace_f)
2115 IF (gradient_f)
THEN
2122 cpassert(
ASSOCIATED(v_xc_tau))
2125 IF (gradient_f)
THEN
2126 ALLOCATE (v_drho_r(3, nspins), v_drho(nspins))
2127 DO ispin = 1, nspins
2129 CALL allocate_pw(v_drho_r(idir, ispin), pw_pool, bo)
2131 CALL allocate_pw(v_drho(ispin), pw_pool, bo)
2135 IF (
ASSOCIATED(pw_pool))
THEN
2136 CALL pw_pool%create_pw(tmp_g)
2137 CALL pw_pool%create_pw(vxc_g)
2140 cpabort(
"XC_DERIV method is not implemented in GAPW")
2145 DO ispin = 1, nspins
2146 v_xc(ispin)%array = 0.0_dp
2150 DO ispin = 1, nspins
2151 v_xc_tau(ispin)%array = 0.0_dp
2155 IF (laplace_f .AND. my_gapw) &
2156 cpabort(
"Laplace-dependent functional not implemented with GAPW!")
2158 IF (my_compute_virial .AND. (gradient_f .OR. laplace_f))
CALL allocate_pw(virial_pw, pw_pool, bo)
2166 IF (do_spinflip)
THEN
2173 IF (gradient_f)
THEN
2175 norm_drho=norm_drho, norm_drhoa=norm_drhoa, norm_drhob=norm_drhob)
2176 IF (do_spinflip)
THEN
2178 CALL calc_drho_from_a(drho1, drho1a)
2181 CALL calc_drho_from_ab(drho1, drho1a, drho1b)
2184 CALL calc_drho_from_ab(drho, drhoa, drhob)
2186 CALL prepare_dr1dr(dra1dra, drhoa, drho1a)
2187 IF (do_spinflip)
THEN
2188 CALL prepare_dr1dr(drb1drb, drhob, drho1a)
2189 CALL prepare_dr1dr(dr1dr, drho, drho1a)
2190 ELSE IF (nspins /= 1)
THEN
2191 CALL prepare_dr1dr(drb1drb, drhob, drho1b)
2192 CALL prepare_dr1dr(dr1dr, drho, drho1)
2194 CALL prepare_dr1dr(drb1drb, drhob, drho1b)
2195 CALL prepare_dr1dr_ab(dr1dr, drhoa, drhob, drho1a, drho1b,
fac)
2198 ALLOCATE (v_drhoa(nspins), v_drhob(nspins))
2199 DO ispin = 1, nspins
2200 CALL allocate_pw(v_drhoa(ispin), pw_pool, bo)
2201 CALL allocate_pw(v_drhob(ispin), pw_pool, bo)
2207 CALL xc_rho_set_get(rho1_set, laplace_rhoa=laplace1a, laplace_rhob=laplace1b)
2209 ALLOCATE (v_laplace(nspins))
2210 DO ispin = 1, nspins
2211 CALL allocate_pw(v_laplace(ispin), pw_pool, bo)
2214 IF (my_compute_virial)
CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob)
2221 IF (do_spinflip)
THEN
2230 IF (
ASSOCIATED(deriv_att))
THEN
2233 IF (
ASSOCIATED(deriv_att))
THEN
2237 DO k = bo(1, 3), bo(2, 3)
2238 DO j = bo(1, 2), bo(2, 2)
2239 DO i = bo(1, 1), bo(2, 1)
2240 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
2241 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2242 (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)/s
2254 IF (.NOT. alda0)
THEN
2256 IF (
ASSOCIATED(deriv_att))
THEN
2259 IF (
ASSOCIATED(deriv_att))
THEN
2263 DO k = bo(1, 3), bo(2, 3)
2264 DO j = bo(1, 2), bo(2, 2)
2265 DO i = bo(1, 1), bo(2, 1)
2266 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
2267 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2268 (deriv_data(i, j, k)*dra1dra(i, j, k) - &
2269 deriv_data2(i, j, k)*drb1drb(i, j, k))/s
2278 ELSE IF (nspins /= 1)
THEN
2282 IF (
ASSOCIATED(deriv_att))
THEN
2286 DO k = bo(1, 3), bo(2, 3)
2287 DO j = bo(1, 2), bo(2, 2)
2288 DO i = bo(1, 1), bo(2, 1)
2289 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2290 deriv_data(i, j, k)*rho1a(i, j, k)
2296 IF (
ASSOCIATED(deriv_att))
THEN
2300 DO k = bo(1, 3), bo(2, 3)
2301 DO j = bo(1, 2), bo(2, 2)
2302 DO i = bo(1, 1), bo(2, 1)
2303 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2304 deriv_data(i, j, k)*rho1b(i, j, k)
2310 IF (
ASSOCIATED(deriv_att))
THEN
2314 DO k = bo(1, 3), bo(2, 3)
2315 DO j = bo(1, 2), bo(2, 2)
2316 DO i = bo(1, 1), bo(2, 1)
2317 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2318 deriv_data(i, j, k)*dr1dr(i, j, k)
2324 IF (
ASSOCIATED(deriv_att))
THEN
2328 DO k = bo(1, 3), bo(2, 3)
2329 DO j = bo(1, 2), bo(2, 2)
2330 DO i = bo(1, 1), bo(2, 1)
2331 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2332 deriv_data(i, j, k)*dra1dra(i, j, k)
2338 IF (
ASSOCIATED(deriv_att))
THEN
2342 DO k = bo(1, 3), bo(2, 3)
2343 DO j = bo(1, 2), bo(2, 2)
2344 DO i = bo(1, 1), bo(2, 1)
2345 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2346 deriv_data(i, j, k)*drb1drb(i, j, k)
2352 IF (
ASSOCIATED(deriv_att))
THEN
2356 DO k = bo(1, 3), bo(2, 3)
2357 DO j = bo(1, 2), bo(2, 2)
2358 DO i = bo(1, 1), bo(2, 1)
2359 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2360 deriv_data(i, j, k)*tau1a(i, j, k)
2366 IF (
ASSOCIATED(deriv_att))
THEN
2370 DO k = bo(1, 3), bo(2, 3)
2371 DO j = bo(1, 2), bo(2, 2)
2372 DO i = bo(1, 1), bo(2, 1)
2373 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2374 deriv_data(i, j, k)*tau1b(i, j, k)
2380 IF (
ASSOCIATED(deriv_att))
THEN
2384 DO k = bo(1, 3), bo(2, 3)
2385 DO j = bo(1, 2), bo(2, 2)
2386 DO i = bo(1, 1), bo(2, 1)
2387 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2388 deriv_data(i, j, k)*laplace1a(i, j, k)
2394 IF (
ASSOCIATED(deriv_att))
THEN
2398 DO k = bo(1, 3), bo(2, 3)
2399 DO j = bo(1, 2), bo(2, 2)
2400 DO i = bo(1, 1), bo(2, 1)
2401 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2402 deriv_data(i, j, k)*laplace1b(i, j, k)
2410 IF (
ASSOCIATED(deriv_att))
THEN
2414 DO k = bo(1, 3), bo(2, 3)
2415 DO j = bo(1, 2), bo(2, 2)
2416 DO i = bo(1, 1), bo(2, 1)
2417 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2418 deriv_data(i, j, k)*rho1a(i, j, k)
2424 IF (
ASSOCIATED(deriv_att))
THEN
2428 DO k = bo(1, 3), bo(2, 3)
2429 DO j = bo(1, 2), bo(2, 2)
2430 DO i = bo(1, 1), bo(2, 1)
2431 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2432 deriv_data(i, j, k)*rho1b(i, j, k)
2438 IF (
ASSOCIATED(deriv_att))
THEN
2442 DO k = bo(1, 3), bo(2, 3)
2443 DO j = bo(1, 2), bo(2, 2)
2444 DO i = bo(1, 1), bo(2, 1)
2445 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2446 deriv_data(i, j, k)*dr1dr(i, j, k)
2452 IF (
ASSOCIATED(deriv_att))
THEN
2456 DO k = bo(1, 3), bo(2, 3)
2457 DO j = bo(1, 2), bo(2, 2)
2458 DO i = bo(1, 1), bo(2, 1)
2459 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2460 deriv_data(i, j, k)*dra1dra(i, j, k)
2466 IF (
ASSOCIATED(deriv_att))
THEN
2470 DO k = bo(1, 3), bo(2, 3)
2471 DO j = bo(1, 2), bo(2, 2)
2472 DO i = bo(1, 1), bo(2, 1)
2473 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2474 deriv_data(i, j, k)*drb1drb(i, j, k)
2480 IF (
ASSOCIATED(deriv_att))
THEN
2484 DO k = bo(1, 3), bo(2, 3)
2485 DO j = bo(1, 2), bo(2, 2)
2486 DO i = bo(1, 1), bo(2, 1)
2487 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2488 deriv_data(i, j, k)*tau1a(i, j, k)
2494 IF (
ASSOCIATED(deriv_att))
THEN
2498 DO k = bo(1, 3), bo(2, 3)
2499 DO j = bo(1, 2), bo(2, 2)
2500 DO i = bo(1, 1), bo(2, 1)
2501 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2502 deriv_data(i, j, k)*tau1b(i, j, k)
2508 IF (
ASSOCIATED(deriv_att))
THEN
2512 DO k = bo(1, 3), bo(2, 3)
2513 DO j = bo(1, 2), bo(2, 2)
2514 DO i = bo(1, 1), bo(2, 1)
2515 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2516 deriv_data(i, j, k)*laplace1a(i, j, k)
2522 IF (
ASSOCIATED(deriv_att))
THEN
2526 DO k = bo(1, 3), bo(2, 3)
2527 DO j = bo(1, 2), bo(2, 2)
2528 DO i = bo(1, 1), bo(2, 1)
2529 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
2530 deriv_data(i, j, k)*laplace1b(i, j, k)
2538 IF (
ASSOCIATED(deriv_att))
THEN
2542 DO k = bo(1, 3), bo(2, 3)
2543 DO j = bo(1, 2), bo(2, 2)
2544 DO i = bo(1, 1), bo(2, 1)
2545 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2546 deriv_data(i, j, k)*rho1a(i, j, k)
2547 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2548 deriv_data(i, j, k)*rho1a(i, j, k)
2554 IF (
ASSOCIATED(deriv_att))
THEN
2558 DO k = bo(1, 3), bo(2, 3)
2559 DO j = bo(1, 2), bo(2, 2)
2560 DO i = bo(1, 1), bo(2, 1)
2561 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2562 deriv_data(i, j, k)*rho1b(i, j, k)
2563 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2564 deriv_data(i, j, k)*rho1b(i, j, k)
2570 IF (
ASSOCIATED(deriv_att))
THEN
2574 DO k = bo(1, 3), bo(2, 3)
2575 DO j = bo(1, 2), bo(2, 2)
2576 DO i = bo(1, 1), bo(2, 1)
2577 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2578 deriv_data(i, j, k)*dr1dr(i, j, k)
2579 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2580 deriv_data(i, j, k)*dr1dr(i, j, k)
2586 IF (
ASSOCIATED(deriv_att))
THEN
2590 DO k = bo(1, 3), bo(2, 3)
2591 DO j = bo(1, 2), bo(2, 2)
2592 DO i = bo(1, 1), bo(2, 1)
2593 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2594 deriv_data(i, j, k)*dra1dra(i, j, k)
2595 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2596 deriv_data(i, j, k)*dra1dra(i, j, k)
2602 IF (
ASSOCIATED(deriv_att))
THEN
2606 DO k = bo(1, 3), bo(2, 3)
2607 DO j = bo(1, 2), bo(2, 2)
2608 DO i = bo(1, 1), bo(2, 1)
2609 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2610 deriv_data(i, j, k)*drb1drb(i, j, k)
2611 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2612 deriv_data(i, j, k)*drb1drb(i, j, k)
2618 IF (
ASSOCIATED(deriv_att))
THEN
2622 DO k = bo(1, 3), bo(2, 3)
2623 DO j = bo(1, 2), bo(2, 2)
2624 DO i = bo(1, 1), bo(2, 1)
2625 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2626 deriv_data(i, j, k)*tau1a(i, j, k)
2627 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2628 deriv_data(i, j, k)*tau1a(i, j, k)
2634 IF (
ASSOCIATED(deriv_att))
THEN
2638 DO k = bo(1, 3), bo(2, 3)
2639 DO j = bo(1, 2), bo(2, 2)
2640 DO i = bo(1, 1), bo(2, 1)
2641 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2642 deriv_data(i, j, k)*tau1b(i, j, k)
2643 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2644 deriv_data(i, j, k)*tau1b(i, j, k)
2650 IF (
ASSOCIATED(deriv_att))
THEN
2654 DO k = bo(1, 3), bo(2, 3)
2655 DO j = bo(1, 2), bo(2, 2)
2656 DO i = bo(1, 1), bo(2, 1)
2657 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2658 deriv_data(i, j, k)*laplace1a(i, j, k)
2659 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2660 deriv_data(i, j, k)*laplace1a(i, j, k)
2666 IF (
ASSOCIATED(deriv_att))
THEN
2670 DO k = bo(1, 3), bo(2, 3)
2671 DO j = bo(1, 2), bo(2, 2)
2672 DO i = bo(1, 1), bo(2, 1)
2673 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
2674 deriv_data(i, j, k)*laplace1b(i, j, k)
2675 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
2676 deriv_data(i, j, k)*laplace1b(i, j, k)
2683 IF (
ASSOCIATED(deriv_att))
THEN
2687 IF (my_compute_virial)
THEN
2688 CALL virial_drho_drho1(virial_pw, drho, drho1, deriv_data, virial_xc)
2692 v_drho(1)%array(:, :, :) = v_drho(1)%array(:, :, :) + &
2693 deriv_data(:, :, :)*dr1dr(:, :, :)/max(gradient_cut, norm_drho(:, :, :))**2
2694 v_drho(2)%array(:, :, :) = v_drho(2)%array(:, :, :) + &
2695 deriv_data(:, :, :)*dr1dr(:, :, :)/max(gradient_cut, norm_drho(:, :, :))**2
2700 IF (
ASSOCIATED(deriv_att))
THEN
2704 DO k = bo(1, 3), bo(2, 3)
2705 DO j = bo(1, 2), bo(2, 2)
2706 DO i = bo(1, 1), bo(2, 1)
2707 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2708 deriv_data(i, j, k)*rho1a(i, j, k)
2714 IF (
ASSOCIATED(deriv_att))
THEN
2718 DO k = bo(1, 3), bo(2, 3)
2719 DO j = bo(1, 2), bo(2, 2)
2720 DO i = bo(1, 1), bo(2, 1)
2721 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2722 deriv_data(i, j, k)*rho1b(i, j, k)
2728 IF (
ASSOCIATED(deriv_att))
THEN
2732 DO k = bo(1, 3), bo(2, 3)
2733 DO j = bo(1, 2), bo(2, 2)
2734 DO i = bo(1, 1), bo(2, 1)
2735 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2736 deriv_data(i, j, k)*dr1dr(i, j, k)
2742 IF (
ASSOCIATED(deriv_att))
THEN
2746 DO k = bo(1, 3), bo(2, 3)
2747 DO j = bo(1, 2), bo(2, 2)
2748 DO i = bo(1, 1), bo(2, 1)
2749 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2750 deriv_data(i, j, k)*dra1dra(i, j, k)
2756 IF (
ASSOCIATED(deriv_att))
THEN
2760 DO k = bo(1, 3), bo(2, 3)
2761 DO j = bo(1, 2), bo(2, 2)
2762 DO i = bo(1, 1), bo(2, 1)
2763 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2764 deriv_data(i, j, k)*drb1drb(i, j, k)
2770 IF (
ASSOCIATED(deriv_att))
THEN
2774 DO k = bo(1, 3), bo(2, 3)
2775 DO j = bo(1, 2), bo(2, 2)
2776 DO i = bo(1, 1), bo(2, 1)
2777 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2778 deriv_data(i, j, k)*tau1a(i, j, k)
2784 IF (
ASSOCIATED(deriv_att))
THEN
2788 DO k = bo(1, 3), bo(2, 3)
2789 DO j = bo(1, 2), bo(2, 2)
2790 DO i = bo(1, 1), bo(2, 1)
2791 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2792 deriv_data(i, j, k)*tau1b(i, j, k)
2798 IF (
ASSOCIATED(deriv_att))
THEN
2802 DO k = bo(1, 3), bo(2, 3)
2803 DO j = bo(1, 2), bo(2, 2)
2804 DO i = bo(1, 1), bo(2, 1)
2805 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2806 deriv_data(i, j, k)*laplace1a(i, j, k)
2812 IF (
ASSOCIATED(deriv_att))
THEN
2816 DO k = bo(1, 3), bo(2, 3)
2817 DO j = bo(1, 2), bo(2, 2)
2818 DO i = bo(1, 1), bo(2, 1)
2819 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
2820 deriv_data(i, j, k)*laplace1b(i, j, k)
2827 IF (
ASSOCIATED(deriv_att))
THEN
2831 IF (my_compute_virial)
THEN
2832 CALL virial_drho_drho1(virial_pw, drhoa, drho1a, deriv_data, virial_xc)
2836 v_drhoa(1)%array(:, :, :) = v_drhoa(1)%array(:, :, :) + &
2837 deriv_data(:, :, :)*dra1dra(:, :, :)/max(gradient_cut, norm_drhoa(:, :, :))**2
2842 IF (
ASSOCIATED(deriv_att))
THEN
2846 DO k = bo(1, 3), bo(2, 3)
2847 DO j = bo(1, 2), bo(2, 2)
2848 DO i = bo(1, 1), bo(2, 1)
2849 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2850 deriv_data(i, j, k)*rho1a(i, j, k)
2856 IF (
ASSOCIATED(deriv_att))
THEN
2860 DO k = bo(1, 3), bo(2, 3)
2861 DO j = bo(1, 2), bo(2, 2)
2862 DO i = bo(1, 1), bo(2, 1)
2863 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2864 deriv_data(i, j, k)*rho1b(i, j, k)
2870 IF (
ASSOCIATED(deriv_att))
THEN
2874 DO k = bo(1, 3), bo(2, 3)
2875 DO j = bo(1, 2), bo(2, 2)
2876 DO i = bo(1, 1), bo(2, 1)
2877 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2878 deriv_data(i, j, k)*dr1dr(i, j, k)
2884 IF (
ASSOCIATED(deriv_att))
THEN
2888 DO k = bo(1, 3), bo(2, 3)
2889 DO j = bo(1, 2), bo(2, 2)
2890 DO i = bo(1, 1), bo(2, 1)
2891 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2892 deriv_data(i, j, k)*dra1dra(i, j, k)
2898 IF (
ASSOCIATED(deriv_att))
THEN
2902 DO k = bo(1, 3), bo(2, 3)
2903 DO j = bo(1, 2), bo(2, 2)
2904 DO i = bo(1, 1), bo(2, 1)
2905 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2906 deriv_data(i, j, k)*drb1drb(i, j, k)
2912 IF (
ASSOCIATED(deriv_att))
THEN
2916 DO k = bo(1, 3), bo(2, 3)
2917 DO j = bo(1, 2), bo(2, 2)
2918 DO i = bo(1, 1), bo(2, 1)
2919 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2920 deriv_data(i, j, k)*tau1a(i, j, k)
2926 IF (
ASSOCIATED(deriv_att))
THEN
2930 DO k = bo(1, 3), bo(2, 3)
2931 DO j = bo(1, 2), bo(2, 2)
2932 DO i = bo(1, 1), bo(2, 1)
2933 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2934 deriv_data(i, j, k)*tau1b(i, j, k)
2940 IF (
ASSOCIATED(deriv_att))
THEN
2944 DO k = bo(1, 3), bo(2, 3)
2945 DO j = bo(1, 2), bo(2, 2)
2946 DO i = bo(1, 1), bo(2, 1)
2947 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2948 deriv_data(i, j, k)*laplace1a(i, j, k)
2954 IF (
ASSOCIATED(deriv_att))
THEN
2958 DO k = bo(1, 3), bo(2, 3)
2959 DO j = bo(1, 2), bo(2, 2)
2960 DO i = bo(1, 1), bo(2, 1)
2961 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
2962 deriv_data(i, j, k)*laplace1b(i, j, k)
2969 IF (
ASSOCIATED(deriv_att))
THEN
2973 IF (my_compute_virial)
THEN
2974 CALL virial_drho_drho1(virial_pw, drhob, drho1b, deriv_data, virial_xc)
2978 v_drhob(2)%array(:, :, :) = v_drhob(2)%array(:, :, :) + &
2979 deriv_data(:, :, :)*drb1drb(:, :, :)/max(gradient_cut, norm_drhob(:, :, :))**2
2984 IF (
ASSOCIATED(deriv_att))
THEN
2988 DO k = bo(1, 3), bo(2, 3)
2989 DO j = bo(1, 2), bo(2, 2)
2990 DO i = bo(1, 1), bo(2, 1)
2991 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
2992 deriv_data(i, j, k)*rho1a(i, j, k)
2998 IF (
ASSOCIATED(deriv_att))
THEN
3002 DO k = bo(1, 3), bo(2, 3)
3003 DO j = bo(1, 2), bo(2, 2)
3004 DO i = bo(1, 1), bo(2, 1)
3005 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3006 deriv_data(i, j, k)*rho1b(i, j, k)
3012 IF (
ASSOCIATED(deriv_att))
THEN
3016 DO k = bo(1, 3), bo(2, 3)
3017 DO j = bo(1, 2), bo(2, 2)
3018 DO i = bo(1, 1), bo(2, 1)
3019 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3020 deriv_data(i, j, k)*dr1dr(i, j, k)
3026 IF (
ASSOCIATED(deriv_att))
THEN
3030 DO k = bo(1, 3), bo(2, 3)
3031 DO j = bo(1, 2), bo(2, 2)
3032 DO i = bo(1, 1), bo(2, 1)
3033 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3034 deriv_data(i, j, k)*dra1dra(i, j, k)
3040 IF (
ASSOCIATED(deriv_att))
THEN
3044 DO k = bo(1, 3), bo(2, 3)
3045 DO j = bo(1, 2), bo(2, 2)
3046 DO i = bo(1, 1), bo(2, 1)
3047 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3048 deriv_data(i, j, k)*drb1drb(i, j, k)
3054 IF (
ASSOCIATED(deriv_att))
THEN
3058 DO k = bo(1, 3), bo(2, 3)
3059 DO j = bo(1, 2), bo(2, 2)
3060 DO i = bo(1, 1), bo(2, 1)
3061 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3062 deriv_data(i, j, k)*tau1a(i, j, k)
3068 IF (
ASSOCIATED(deriv_att))
THEN
3072 DO k = bo(1, 3), bo(2, 3)
3073 DO j = bo(1, 2), bo(2, 2)
3074 DO i = bo(1, 1), bo(2, 1)
3075 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3076 deriv_data(i, j, k)*tau1b(i, j, k)
3082 IF (
ASSOCIATED(deriv_att))
THEN
3086 DO k = bo(1, 3), bo(2, 3)
3087 DO j = bo(1, 2), bo(2, 2)
3088 DO i = bo(1, 1), bo(2, 1)
3089 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3090 deriv_data(i, j, k)*laplace1a(i, j, k)
3096 IF (
ASSOCIATED(deriv_att))
THEN
3100 DO k = bo(1, 3), bo(2, 3)
3101 DO j = bo(1, 2), bo(2, 2)
3102 DO i = bo(1, 1), bo(2, 1)
3103 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3104 deriv_data(i, j, k)*laplace1b(i, j, k)
3112 IF (
ASSOCIATED(deriv_att))
THEN
3116 DO k = bo(1, 3), bo(2, 3)
3117 DO j = bo(1, 2), bo(2, 2)
3118 DO i = bo(1, 1), bo(2, 1)
3119 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3120 deriv_data(i, j, k)*rho1a(i, j, k)
3126 IF (
ASSOCIATED(deriv_att))
THEN
3130 DO k = bo(1, 3), bo(2, 3)
3131 DO j = bo(1, 2), bo(2, 2)
3132 DO i = bo(1, 1), bo(2, 1)
3133 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3134 deriv_data(i, j, k)*rho1b(i, j, k)
3140 IF (
ASSOCIATED(deriv_att))
THEN
3144 DO k = bo(1, 3), bo(2, 3)
3145 DO j = bo(1, 2), bo(2, 2)
3146 DO i = bo(1, 1), bo(2, 1)
3147 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3148 deriv_data(i, j, k)*dr1dr(i, j, k)
3154 IF (
ASSOCIATED(deriv_att))
THEN
3158 DO k = bo(1, 3), bo(2, 3)
3159 DO j = bo(1, 2), bo(2, 2)
3160 DO i = bo(1, 1), bo(2, 1)
3161 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3162 deriv_data(i, j, k)*dra1dra(i, j, k)
3168 IF (
ASSOCIATED(deriv_att))
THEN
3172 DO k = bo(1, 3), bo(2, 3)
3173 DO j = bo(1, 2), bo(2, 2)
3174 DO i = bo(1, 1), bo(2, 1)
3175 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3176 deriv_data(i, j, k)*drb1drb(i, j, k)
3182 IF (
ASSOCIATED(deriv_att))
THEN
3186 DO k = bo(1, 3), bo(2, 3)
3187 DO j = bo(1, 2), bo(2, 2)
3188 DO i = bo(1, 1), bo(2, 1)
3189 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3190 deriv_data(i, j, k)*tau1a(i, j, k)
3196 IF (
ASSOCIATED(deriv_att))
THEN
3200 DO k = bo(1, 3), bo(2, 3)
3201 DO j = bo(1, 2), bo(2, 2)
3202 DO i = bo(1, 1), bo(2, 1)
3203 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3204 deriv_data(i, j, k)*tau1b(i, j, k)
3210 IF (
ASSOCIATED(deriv_att))
THEN
3214 DO k = bo(1, 3), bo(2, 3)
3215 DO j = bo(1, 2), bo(2, 2)
3216 DO i = bo(1, 1), bo(2, 1)
3217 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3218 deriv_data(i, j, k)*laplace1a(i, j, k)
3224 IF (
ASSOCIATED(deriv_att))
THEN
3228 DO k = bo(1, 3), bo(2, 3)
3229 DO j = bo(1, 2), bo(2, 2)
3230 DO i = bo(1, 1), bo(2, 1)
3231 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3232 deriv_data(i, j, k)*laplace1b(i, j, k)
3240 IF (
ASSOCIATED(deriv_att))
THEN
3244 DO k = bo(1, 3), bo(2, 3)
3245 DO j = bo(1, 2), bo(2, 2)
3246 DO i = bo(1, 1), bo(2, 1)
3247 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3248 deriv_data(i, j, k)*rho1a(i, j, k)
3254 IF (
ASSOCIATED(deriv_att))
THEN
3258 DO k = bo(1, 3), bo(2, 3)
3259 DO j = bo(1, 2), bo(2, 2)
3260 DO i = bo(1, 1), bo(2, 1)
3261 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3262 deriv_data(i, j, k)*rho1b(i, j, k)
3268 IF (
ASSOCIATED(deriv_att))
THEN
3272 DO k = bo(1, 3), bo(2, 3)
3273 DO j = bo(1, 2), bo(2, 2)
3274 DO i = bo(1, 1), bo(2, 1)
3275 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3276 deriv_data(i, j, k)*dr1dr(i, j, k)
3282 IF (
ASSOCIATED(deriv_att))
THEN
3286 DO k = bo(1, 3), bo(2, 3)
3287 DO j = bo(1, 2), bo(2, 2)
3288 DO i = bo(1, 1), bo(2, 1)
3289 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3290 deriv_data(i, j, k)*dra1dra(i, j, k)
3296 IF (
ASSOCIATED(deriv_att))
THEN
3300 DO k = bo(1, 3), bo(2, 3)
3301 DO j = bo(1, 2), bo(2, 2)
3302 DO i = bo(1, 1), bo(2, 1)
3303 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3304 deriv_data(i, j, k)*drb1drb(i, j, k)
3310 IF (
ASSOCIATED(deriv_att))
THEN
3314 DO k = bo(1, 3), bo(2, 3)
3315 DO j = bo(1, 2), bo(2, 2)
3316 DO i = bo(1, 1), bo(2, 1)
3317 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3318 deriv_data(i, j, k)*tau1a(i, j, k)
3324 IF (
ASSOCIATED(deriv_att))
THEN
3328 DO k = bo(1, 3), bo(2, 3)
3329 DO j = bo(1, 2), bo(2, 2)
3330 DO i = bo(1, 1), bo(2, 1)
3331 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3332 deriv_data(i, j, k)*tau1b(i, j, k)
3338 IF (
ASSOCIATED(deriv_att))
THEN
3342 DO k = bo(1, 3), bo(2, 3)
3343 DO j = bo(1, 2), bo(2, 2)
3344 DO i = bo(1, 1), bo(2, 1)
3345 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3346 deriv_data(i, j, k)*laplace1a(i, j, k)
3352 IF (
ASSOCIATED(deriv_att))
THEN
3356 DO k = bo(1, 3), bo(2, 3)
3357 DO j = bo(1, 2), bo(2, 2)
3358 DO i = bo(1, 1), bo(2, 1)
3359 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3360 deriv_data(i, j, k)*laplace1b(i, j, k)
3367 IF (my_compute_virial)
THEN
3369 IF (
ASSOCIATED(deriv_att))
THEN
3372 virial_pw%array(:, :, :) = -rho1a(:, :, :)
3373 CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
3377 IF (
ASSOCIATED(deriv_att))
THEN
3381 DO k = bo(1, 3), bo(2, 3)
3382 DO j = bo(1, 2), bo(2, 2)
3383 DO i = bo(1, 1), bo(2, 1)
3384 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3385 deriv_data(i, j, k)*rho1a(i, j, k)
3391 IF (
ASSOCIATED(deriv_att))
THEN
3395 DO k = bo(1, 3), bo(2, 3)
3396 DO j = bo(1, 2), bo(2, 2)
3397 DO i = bo(1, 1), bo(2, 1)
3398 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3399 deriv_data(i, j, k)*rho1b(i, j, k)
3405 IF (
ASSOCIATED(deriv_att))
THEN
3409 DO k = bo(1, 3), bo(2, 3)
3410 DO j = bo(1, 2), bo(2, 2)
3411 DO i = bo(1, 1), bo(2, 1)
3412 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3413 deriv_data(i, j, k)*dr1dr(i, j, k)
3419 IF (
ASSOCIATED(deriv_att))
THEN
3423 DO k = bo(1, 3), bo(2, 3)
3424 DO j = bo(1, 2), bo(2, 2)
3425 DO i = bo(1, 1), bo(2, 1)
3426 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3427 deriv_data(i, j, k)*dra1dra(i, j, k)
3433 IF (
ASSOCIATED(deriv_att))
THEN
3437 DO k = bo(1, 3), bo(2, 3)
3438 DO j = bo(1, 2), bo(2, 2)
3439 DO i = bo(1, 1), bo(2, 1)
3440 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3441 deriv_data(i, j, k)*drb1drb(i, j, k)
3447 IF (
ASSOCIATED(deriv_att))
THEN
3451 DO k = bo(1, 3), bo(2, 3)
3452 DO j = bo(1, 2), bo(2, 2)
3453 DO i = bo(1, 1), bo(2, 1)
3454 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3455 deriv_data(i, j, k)*tau1a(i, j, k)
3461 IF (
ASSOCIATED(deriv_att))
THEN
3465 DO k = bo(1, 3), bo(2, 3)
3466 DO j = bo(1, 2), bo(2, 2)
3467 DO i = bo(1, 1), bo(2, 1)
3468 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3469 deriv_data(i, j, k)*tau1b(i, j, k)
3475 IF (
ASSOCIATED(deriv_att))
THEN
3479 DO k = bo(1, 3), bo(2, 3)
3480 DO j = bo(1, 2), bo(2, 2)
3481 DO i = bo(1, 1), bo(2, 1)
3482 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3483 deriv_data(i, j, k)*laplace1a(i, j, k)
3489 IF (
ASSOCIATED(deriv_att))
THEN
3493 DO k = bo(1, 3), bo(2, 3)
3494 DO j = bo(1, 2), bo(2, 2)
3495 DO i = bo(1, 1), bo(2, 1)
3496 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3497 deriv_data(i, j, k)*laplace1b(i, j, k)
3504 IF (my_compute_virial)
THEN
3506 IF (
ASSOCIATED(deriv_att))
THEN
3509 virial_pw%array(:, :, :) = -rho1b(:, :, :)
3510 CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
3519 IF (
ASSOCIATED(deriv_att))
THEN
3523 DO k = bo(1, 3), bo(2, 3)
3524 DO j = bo(1, 2), bo(2, 2)
3525 DO i = bo(1, 1), bo(2, 1)
3526 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3527 deriv_data(i, j, k)*rho1a(i, j, k)
3533 IF (
ASSOCIATED(deriv_att))
THEN
3537 DO k = bo(1, 3), bo(2, 3)
3538 DO j = bo(1, 2), bo(2, 2)
3539 DO i = bo(1, 1), bo(2, 1)
3540 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3541 deriv_data(i, j, k)*dr1dr(i, j, k)
3547 IF (
ASSOCIATED(deriv_att))
THEN
3551 DO k = bo(1, 3), bo(2, 3)
3552 DO j = bo(1, 2), bo(2, 2)
3553 DO i = bo(1, 1), bo(2, 1)
3554 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3555 deriv_data(i, j, k)*dra1dra(i, j, k)
3561 IF (
ASSOCIATED(deriv_att))
THEN
3565 DO k = bo(1, 3), bo(2, 3)
3566 DO j = bo(1, 2), bo(2, 2)
3567 DO i = bo(1, 1), bo(2, 1)
3568 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3569 deriv_data(i, j, k)*tau1a(i, j, k)
3575 IF (
ASSOCIATED(deriv_att))
THEN
3579 DO k = bo(1, 3), bo(2, 3)
3580 DO j = bo(1, 2), bo(2, 2)
3581 DO i = bo(1, 1), bo(2, 1)
3582 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3583 deriv_data(i, j, k)*laplace1a(i, j, k)
3589 IF (
ASSOCIATED(deriv_att))
THEN
3593 DO k = bo(1, 3), bo(2, 3)
3594 DO j = bo(1, 2), bo(2, 2)
3595 DO i = bo(1, 1), bo(2, 1)
3596 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3597 fac*deriv_data(i, j, k)*rho1b(i, j, k)
3603 IF (
ASSOCIATED(deriv_att))
THEN
3607 DO k = bo(1, 3), bo(2, 3)
3608 DO j = bo(1, 2), bo(2, 2)
3609 DO i = bo(1, 1), bo(2, 1)
3610 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3611 fac*deriv_data(i, j, k)*drb1drb(i, j, k)
3617 IF (
ASSOCIATED(deriv_att))
THEN
3621 DO k = bo(1, 3), bo(2, 3)
3622 DO j = bo(1, 2), bo(2, 2)
3623 DO i = bo(1, 1), bo(2, 1)
3624 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3625 fac*deriv_data(i, j, k)*tau1b(i, j, k)
3631 IF (
ASSOCIATED(deriv_att))
THEN
3635 DO k = bo(1, 3), bo(2, 3)
3636 DO j = bo(1, 2), bo(2, 2)
3637 DO i = bo(1, 1), bo(2, 1)
3638 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3639 fac*deriv_data(i, j, k)*laplace1b(i, j, k)
3647 IF (
ASSOCIATED(deriv_att))
THEN
3651 DO k = bo(1, 3), bo(2, 3)
3652 DO j = bo(1, 2), bo(2, 2)
3653 DO i = bo(1, 1), bo(2, 1)
3654 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3655 deriv_data(i, j, k)*rho1a(i, j, k)
3661 IF (
ASSOCIATED(deriv_att))
THEN
3665 DO k = bo(1, 3), bo(2, 3)
3666 DO j = bo(1, 2), bo(2, 2)
3667 DO i = bo(1, 1), bo(2, 1)
3668 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3669 deriv_data(i, j, k)*dr1dr(i, j, k)
3675 IF (
ASSOCIATED(deriv_att))
THEN
3679 DO k = bo(1, 3), bo(2, 3)
3680 DO j = bo(1, 2), bo(2, 2)
3681 DO i = bo(1, 1), bo(2, 1)
3682 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3683 deriv_data(i, j, k)*dra1dra(i, j, k)
3689 IF (
ASSOCIATED(deriv_att))
THEN
3693 DO k = bo(1, 3), bo(2, 3)
3694 DO j = bo(1, 2), bo(2, 2)
3695 DO i = bo(1, 1), bo(2, 1)
3696 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3697 deriv_data(i, j, k)*tau1a(i, j, k)
3703 IF (
ASSOCIATED(deriv_att))
THEN
3707 DO k = bo(1, 3), bo(2, 3)
3708 DO j = bo(1, 2), bo(2, 2)
3709 DO i = bo(1, 1), bo(2, 1)
3710 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3711 deriv_data(i, j, k)*laplace1a(i, j, k)
3717 IF (
ASSOCIATED(deriv_att))
THEN
3721 DO k = bo(1, 3), bo(2, 3)
3722 DO j = bo(1, 2), bo(2, 2)
3723 DO i = bo(1, 1), bo(2, 1)
3724 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3725 fac*deriv_data(i, j, k)*rho1b(i, j, k)
3731 IF (
ASSOCIATED(deriv_att))
THEN
3735 DO k = bo(1, 3), bo(2, 3)
3736 DO j = bo(1, 2), bo(2, 2)
3737 DO i = bo(1, 1), bo(2, 1)
3738 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3739 fac*deriv_data(i, j, k)*drb1drb(i, j, k)
3745 IF (
ASSOCIATED(deriv_att))
THEN
3749 DO k = bo(1, 3), bo(2, 3)
3750 DO j = bo(1, 2), bo(2, 2)
3751 DO i = bo(1, 1), bo(2, 1)
3752 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3753 fac*deriv_data(i, j, k)*tau1b(i, j, k)
3759 IF (
ASSOCIATED(deriv_att))
THEN
3763 DO k = bo(1, 3), bo(2, 3)
3764 DO j = bo(1, 2), bo(2, 2)
3765 DO i = bo(1, 1), bo(2, 1)
3766 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3767 fac*deriv_data(i, j, k)*laplace1b(i, j, k)
3774 IF (
ASSOCIATED(deriv_att))
THEN
3780 v_drho(1)%array(:, :, :) = v_drho(1)%array(:, :, :) + &
3781 deriv_data(:, :, :)*dr1dr(:, :, :)/max(gradient_cut, norm_drho(:, :, :))**2
3786 IF (
ASSOCIATED(deriv_att))
THEN
3790 DO k = bo(1, 3), bo(2, 3)
3791 DO j = bo(1, 2), bo(2, 2)
3792 DO i = bo(1, 1), bo(2, 1)
3793 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3794 deriv_data(i, j, k)*rho1a(i, j, k)
3800 IF (
ASSOCIATED(deriv_att))
THEN
3804 DO k = bo(1, 3), bo(2, 3)
3805 DO j = bo(1, 2), bo(2, 2)
3806 DO i = bo(1, 1), bo(2, 1)
3807 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3808 deriv_data(i, j, k)*dr1dr(i, j, k)
3814 IF (
ASSOCIATED(deriv_att))
THEN
3818 DO k = bo(1, 3), bo(2, 3)
3819 DO j = bo(1, 2), bo(2, 2)
3820 DO i = bo(1, 1), bo(2, 1)
3821 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3822 deriv_data(i, j, k)*dra1dra(i, j, k)
3828 IF (
ASSOCIATED(deriv_att))
THEN
3832 DO k = bo(1, 3), bo(2, 3)
3833 DO j = bo(1, 2), bo(2, 2)
3834 DO i = bo(1, 1), bo(2, 1)
3835 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3836 deriv_data(i, j, k)*tau1a(i, j, k)
3842 IF (
ASSOCIATED(deriv_att))
THEN
3846 DO k = bo(1, 3), bo(2, 3)
3847 DO j = bo(1, 2), bo(2, 2)
3848 DO i = bo(1, 1), bo(2, 1)
3849 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3850 deriv_data(i, j, k)*laplace1a(i, j, k)
3856 IF (
ASSOCIATED(deriv_att))
THEN
3860 DO k = bo(1, 3), bo(2, 3)
3861 DO j = bo(1, 2), bo(2, 2)
3862 DO i = bo(1, 1), bo(2, 1)
3863 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3864 fac*deriv_data(i, j, k)*rho1b(i, j, k)
3870 IF (
ASSOCIATED(deriv_att))
THEN
3874 DO k = bo(1, 3), bo(2, 3)
3875 DO j = bo(1, 2), bo(2, 2)
3876 DO i = bo(1, 1), bo(2, 1)
3877 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3878 fac*deriv_data(i, j, k)*drb1drb(i, j, k)
3884 IF (
ASSOCIATED(deriv_att))
THEN
3888 DO k = bo(1, 3), bo(2, 3)
3889 DO j = bo(1, 2), bo(2, 2)
3890 DO i = bo(1, 1), bo(2, 1)
3891 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3892 fac*deriv_data(i, j, k)*tau1b(i, j, k)
3898 IF (
ASSOCIATED(deriv_att))
THEN
3902 DO k = bo(1, 3), bo(2, 3)
3903 DO j = bo(1, 2), bo(2, 2)
3904 DO i = bo(1, 1), bo(2, 1)
3905 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3906 fac*deriv_data(i, j, k)*laplace1b(i, j, k)
3913 IF (
ASSOCIATED(deriv_att))
THEN
3919 v_drhoa(1)%array(:, :, :) = v_drhoa(1)%array(:, :, :) + &
3920 deriv_data(:, :, :)*dra1dra(:, :, :)/max(gradient_cut, norm_drhoa(:, :, :))**2
3925 IF (
ASSOCIATED(deriv_att))
THEN
3929 DO k = bo(1, 3), bo(2, 3)
3930 DO j = bo(1, 2), bo(2, 2)
3931 DO i = bo(1, 1), bo(2, 1)
3932 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3933 deriv_data(i, j, k)*rho1a(i, j, k)
3939 IF (
ASSOCIATED(deriv_att))
THEN
3943 DO k = bo(1, 3), bo(2, 3)
3944 DO j = bo(1, 2), bo(2, 2)
3945 DO i = bo(1, 1), bo(2, 1)
3946 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3947 deriv_data(i, j, k)*dr1dr(i, j, k)
3953 IF (
ASSOCIATED(deriv_att))
THEN
3957 DO k = bo(1, 3), bo(2, 3)
3958 DO j = bo(1, 2), bo(2, 2)
3959 DO i = bo(1, 1), bo(2, 1)
3960 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3961 deriv_data(i, j, k)*dra1dra(i, j, k)
3967 IF (
ASSOCIATED(deriv_att))
THEN
3971 DO k = bo(1, 3), bo(2, 3)
3972 DO j = bo(1, 2), bo(2, 2)
3973 DO i = bo(1, 1), bo(2, 1)
3974 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3975 deriv_data(i, j, k)*tau1a(i, j, k)
3981 IF (
ASSOCIATED(deriv_att))
THEN
3985 DO k = bo(1, 3), bo(2, 3)
3986 DO j = bo(1, 2), bo(2, 2)
3987 DO i = bo(1, 1), bo(2, 1)
3988 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3989 deriv_data(i, j, k)*laplace1a(i, j, k)
3995 IF (
ASSOCIATED(deriv_att))
THEN
3999 DO k = bo(1, 3), bo(2, 3)
4000 DO j = bo(1, 2), bo(2, 2)
4001 DO i = bo(1, 1), bo(2, 1)
4002 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4003 fac*deriv_data(i, j, k)*rho1b(i, j, k)
4009 IF (
ASSOCIATED(deriv_att))
THEN
4013 DO k = bo(1, 3), bo(2, 3)
4014 DO j = bo(1, 2), bo(2, 2)
4015 DO i = bo(1, 1), bo(2, 1)
4016 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4017 fac*deriv_data(i, j, k)*drb1drb(i, j, k)
4023 IF (
ASSOCIATED(deriv_att))
THEN
4027 DO k = bo(1, 3), bo(2, 3)
4028 DO j = bo(1, 2), bo(2, 2)
4029 DO i = bo(1, 1), bo(2, 1)
4030 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4031 fac*deriv_data(i, j, k)*tau1b(i, j, k)
4037 IF (
ASSOCIATED(deriv_att))
THEN
4041 DO k = bo(1, 3), bo(2, 3)
4042 DO j = bo(1, 2), bo(2, 2)
4043 DO i = bo(1, 1), bo(2, 1)
4044 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4045 fac*deriv_data(i, j, k)*laplace1b(i, j, k)
4053 IF (
ASSOCIATED(deriv_att))
THEN
4057 DO k = bo(1, 3), bo(2, 3)
4058 DO j = bo(1, 2), bo(2, 2)
4059 DO i = bo(1, 1), bo(2, 1)
4060 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4061 deriv_data(i, j, k)*rho1a(i, j, k)
4067 IF (
ASSOCIATED(deriv_att))
THEN
4071 DO k = bo(1, 3), bo(2, 3)
4072 DO j = bo(1, 2), bo(2, 2)
4073 DO i = bo(1, 1), bo(2, 1)
4074 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4075 deriv_data(i, j, k)*dr1dr(i, j, k)
4081 IF (
ASSOCIATED(deriv_att))
THEN
4085 DO k = bo(1, 3), bo(2, 3)
4086 DO j = bo(1, 2), bo(2, 2)
4087 DO i = bo(1, 1), bo(2, 1)
4088 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4089 deriv_data(i, j, k)*dra1dra(i, j, k)
4095 IF (
ASSOCIATED(deriv_att))
THEN
4099 DO k = bo(1, 3), bo(2, 3)
4100 DO j = bo(1, 2), bo(2, 2)
4101 DO i = bo(1, 1), bo(2, 1)
4102 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4103 deriv_data(i, j, k)*tau1a(i, j, k)
4109 IF (
ASSOCIATED(deriv_att))
THEN
4113 DO k = bo(1, 3), bo(2, 3)
4114 DO j = bo(1, 2), bo(2, 2)
4115 DO i = bo(1, 1), bo(2, 1)
4116 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4117 deriv_data(i, j, k)*laplace1a(i, j, k)
4123 IF (
ASSOCIATED(deriv_att))
THEN
4127 DO k = bo(1, 3), bo(2, 3)
4128 DO j = bo(1, 2), bo(2, 2)
4129 DO i = bo(1, 1), bo(2, 1)
4130 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4131 fac*deriv_data(i, j, k)*rho1b(i, j, k)
4137 IF (
ASSOCIATED(deriv_att))
THEN
4141 DO k = bo(1, 3), bo(2, 3)
4142 DO j = bo(1, 2), bo(2, 2)
4143 DO i = bo(1, 1), bo(2, 1)
4144 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4145 fac*deriv_data(i, j, k)*drb1drb(i, j, k)
4151 IF (
ASSOCIATED(deriv_att))
THEN
4155 DO k = bo(1, 3), bo(2, 3)
4156 DO j = bo(1, 2), bo(2, 2)
4157 DO i = bo(1, 1), bo(2, 1)
4158 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4159 fac*deriv_data(i, j, k)*tau1b(i, j, k)
4165 IF (
ASSOCIATED(deriv_att))
THEN
4169 DO k = bo(1, 3), bo(2, 3)
4170 DO j = bo(1, 2), bo(2, 2)
4171 DO i = bo(1, 1), bo(2, 1)
4172 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4173 fac*deriv_data(i, j, k)*laplace1b(i, j, k)
4184 IF (gradient_f)
THEN
4185 IF (.NOT. do_spinflip)
THEN
4187 IF (my_compute_virial)
THEN
4188 CALL virial_drho_drho(virial_pw, drhoa, v_drhoa(1), virial_xc)
4189 CALL virial_drho_drho(virial_pw, drhob, v_drhob(2), virial_xc)
4192 virial_pw%array(:, :, :) = drho(idir)%array(:, :, :)*(v_drho(1)%array(:, :, :) + v_drho(2)%array(:, :, :))
4196 drho(jdir)%array(:, :, :))
4197 virial_xc(jdir, idir) = virial_xc(jdir, idir) + tmp
4198 virial_xc(idir, jdir) = virial_xc(jdir, idir)
4208 DO ir = bo(1, 2), bo(2, 2)
4209 DO ia = bo(1, 1), bo(2, 1)
4211 DO ispin = 1, nspins
4212 vxg(idir, ia, ir, ispin) = &
4213 -(v_drhoa(ispin)%array(ia, ir, 1)*drhoa(idir)%array(ia, ir, 1) + &
4214 v_drhob(ispin)%array(ia, ir, 1)*drhob(idir)%array(ia, ir, 1) + &
4215 v_drho(ispin)%array(ia, ir, 1)*drho(idir)%array(ia, ir, 1))
4217 IF (
ASSOCIATED(e_drhoa))
THEN
4218 vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
4219 e_drhoa(ia, ir, 1)*drho1a(idir)%array(ia, ir, 1)
4221 IF (nspins /= 1 .AND.
ASSOCIATED(e_drhob))
THEN
4222 vxg(idir, ia, ir, 2) = vxg(idir, ia, ir, 2) + &
4223 e_drhob(ia, ir, 1)*drho1b(idir)%array(ia, ir, 1)
4225 IF (
ASSOCIATED(e_drho))
THEN
4226 IF (nspins /= 1)
THEN
4227 vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
4228 e_drho(ia, ir, 1)*drho1(idir)%array(ia, ir, 1)
4229 vxg(idir, ia, ir, 2) = vxg(idir, ia, ir, 2) + &
4230 e_drho(ia, ir, 1)*drho1(idir)%array(ia, ir, 1)
4232 vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
4233 e_drho(ia, ir, 1)*(drho1a(idir)%array(ia, ir, 1) + &
4234 fac*drho1b(idir)%array(ia, ir, 1))
4246 DO ispin = 1, nspins
4249 v_drho_r(idir, ispin)%array(:, :, :) = &
4250 v_drhoa(ispin)%array(:, :, :)*drhoa(idir)%array(:, :, :) + &
4251 v_drhob(ispin)%array(:, :, :)*drhob(idir)%array(:, :, :) + &
4252 v_drho(ispin)%array(:, :, :)*drho(idir)%array(:, :, :)
4255 IF (
ASSOCIATED(e_drhoa))
THEN
4258 v_drho_r(idir, 1)%array(:, :, :) = v_drho_r(idir, 1)%array(:, :, :) - &
4259 e_drhoa(:, :, :)*drho1a(idir)%array(:, :, :)
4262 IF (nspins /= 1 .AND.
ASSOCIATED(e_drhob))
THEN
4265 v_drho_r(idir, 2)%array(:, :, :) = v_drho_r(idir, 2)%array(:, :, :) - &
4266 e_drhob(:, :, :)*drho1b(idir)%array(:, :, :)
4269 IF (
ASSOCIATED(e_drho))
THEN
4272 DO k = bo(1, 3), bo(2, 3)
4273 DO j = bo(1, 2), bo(2, 2)
4274 DO i = bo(1, 1), bo(2, 1)
4275 IF (nspins /= 1)
THEN
4276 v_drho_r(idir, 1)%array(i, j, k) = v_drho_r(idir, 1)%array(i, j, k) - &
4277 e_drho(i, j, k)*drho1(idir)%array(i, j, k)
4278 v_drho_r(idir, 2)%array(i, j, k) = v_drho_r(idir, 2)%array(i, j, k) - &
4279 e_drho(i, j, k)*drho1(idir)%array(i, j, k)
4281 v_drho_r(idir, 1)%array(i, j, k) = v_drho_r(idir, 1)%array(i, j, k) - &
4282 e_drho(i, j, k)*(drho1a(idir)%array(i, j, k) + &
4283 fac*drho1b(idir)%array(i, j, k))
4293 DO ispin = 1, nspins
4294 CALL xc_pw_divergence(xc_deriv_method_id, v_drho_r(:, ispin), tmp_g, vxc_g, v_xc(ispin))
4302 DEALLOCATE (drho(idir)%array)
4303 DEALLOCATE (drho1(idir)%array)
4306 DO ispin = 1, nspins
4307 CALL deallocate_pw(v_drhoa(ispin), pw_pool)
4308 CALL deallocate_pw(v_drhob(ispin), pw_pool)
4311 DEALLOCATE (v_drhoa, v_drhob)
4315 IF (laplace_f .AND. my_compute_virial)
THEN
4316 virial_pw%array(:, :, :) = -rhoa(:, :, :)
4317 CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(1)%array)
4318 virial_pw%array(:, :, :) = -rhob(:, :, :)
4319 CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(2)%array)
4330 IF (gradient_f)
THEN
4333 CALL prepare_dr1dr(dr1dr, drho, drho1)
4339 ALLOCATE (v_laplace(nspins))
4340 DO ispin = 1, nspins
4341 CALL allocate_pw(v_laplace(ispin), pw_pool, bo)
4352 IF (
ASSOCIATED(deriv_att))
THEN
4356 DO k = bo(1, 3), bo(2, 3)
4357 DO j = bo(1, 2), bo(2, 2)
4358 DO i = bo(1, 1), bo(2, 1)
4359 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4360 deriv_data(i, j, k)*rho1(i, j, k)
4366 IF (
ASSOCIATED(deriv_att))
THEN
4370 DO k = bo(1, 3), bo(2, 3)
4371 DO j = bo(1, 2), bo(2, 2)
4372 DO i = bo(1, 1), bo(2, 1)
4373 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4374 deriv_data(i, j, k)*dr1dr(i, j, k)
4380 IF (
ASSOCIATED(deriv_att))
THEN
4384 DO k = bo(1, 3), bo(2, 3)
4385 DO j = bo(1, 2), bo(2, 2)
4386 DO i = bo(1, 1), bo(2, 1)
4387 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4388 deriv_data(i, j, k)*tau1(i, j, k)
4394 IF (
ASSOCIATED(deriv_att))
THEN
4398 DO k = bo(1, 3), bo(2, 3)
4399 DO j = bo(1, 2), bo(2, 2)
4400 DO i = bo(1, 1), bo(2, 1)
4401 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4402 deriv_data(i, j, k)*laplace1(i, j, k)
4410 IF (
ASSOCIATED(deriv_att))
THEN
4414 DO k = bo(1, 3), bo(2, 3)
4415 DO j = bo(1, 2), bo(2, 2)
4416 DO i = bo(1, 1), bo(2, 1)
4417 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
4418 deriv_data(i, j, k)*rho1(i, j, k)
4424 IF (
ASSOCIATED(deriv_att))
THEN
4428 DO k = bo(1, 3), bo(2, 3)
4429 DO j = bo(1, 2), bo(2, 2)
4430 DO i = bo(1, 1), bo(2, 1)
4431 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
4432 deriv_data(i, j, k)*dr1dr(i, j, k)
4438 IF (
ASSOCIATED(deriv_att))
THEN
4442 DO k = bo(1, 3), bo(2, 3)
4443 DO j = bo(1, 2), bo(2, 2)
4444 DO i = bo(1, 1), bo(2, 1)
4445 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
4446 deriv_data(i, j, k)*tau1(i, j, k)
4452 IF (
ASSOCIATED(deriv_att))
THEN
4456 DO k = bo(1, 3), bo(2, 3)
4457 DO j = bo(1, 2), bo(2, 2)
4458 DO i = bo(1, 1), bo(2, 1)
4459 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
4460 deriv_data(i, j, k)*laplace1(i, j, k)
4467 IF (
ASSOCIATED(deriv_att))
THEN
4471 IF (my_compute_virial)
THEN
4472 CALL virial_drho_drho1(virial_pw, drho, drho1, deriv_data, virial_xc)
4476 v_drho(1)%array(:, :, :) = v_drho(1)%array(:, :, :) + &
4477 deriv_data(:, :, :)*dr1dr(:, :, :)/max(gradient_cut, norm_drho(:, :, :))**2
4482 IF (
ASSOCIATED(deriv_att))
THEN
4486 DO k = bo(1, 3), bo(2, 3)
4487 DO j = bo(1, 2), bo(2, 2)
4488 DO i = bo(1, 1), bo(2, 1)
4489 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4490 deriv_data(i, j, k)*rho1(i, j, k)
4496 IF (
ASSOCIATED(deriv_att))
THEN
4500 DO k = bo(1, 3), bo(2, 3)
4501 DO j = bo(1, 2), bo(2, 2)
4502 DO i = bo(1, 1), bo(2, 1)
4503 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4504 deriv_data(i, j, k)*dr1dr(i, j, k)
4510 IF (
ASSOCIATED(deriv_att))
THEN
4514 DO k = bo(1, 3), bo(2, 3)
4515 DO j = bo(1, 2), bo(2, 2)
4516 DO i = bo(1, 1), bo(2, 1)
4517 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4518 deriv_data(i, j, k)*tau1(i, j, k)
4524 IF (
ASSOCIATED(deriv_att))
THEN
4528 DO k = bo(1, 3), bo(2, 3)
4529 DO j = bo(1, 2), bo(2, 2)
4530 DO i = bo(1, 1), bo(2, 1)
4531 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4532 deriv_data(i, j, k)*laplace1(i, j, k)
4540 IF (
ASSOCIATED(deriv_att))
THEN
4544 DO k = bo(1, 3), bo(2, 3)
4545 DO j = bo(1, 2), bo(2, 2)
4546 DO i = bo(1, 1), bo(2, 1)
4547 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
4548 deriv_data(i, j, k)*rho1(i, j, k)
4554 IF (
ASSOCIATED(deriv_att))
THEN
4558 DO k = bo(1, 3), bo(2, 3)
4559 DO j = bo(1, 2), bo(2, 2)
4560 DO i = bo(1, 1), bo(2, 1)
4561 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
4562 deriv_data(i, j, k)*dr1dr(i, j, k)
4568 IF (
ASSOCIATED(deriv_att))
THEN
4572 DO k = bo(1, 3), bo(2, 3)
4573 DO j = bo(1, 2), bo(2, 2)
4574 DO i = bo(1, 1), bo(2, 1)
4575 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
4576 deriv_data(i, j, k)*tau1(i, j, k)
4582 IF (
ASSOCIATED(deriv_att))
THEN
4586 DO k = bo(1, 3), bo(2, 3)
4587 DO j = bo(1, 2), bo(2, 2)
4588 DO i = bo(1, 1), bo(2, 1)
4589 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
4590 deriv_data(i, j, k)*laplace1(i, j, k)
4597 IF (my_compute_virial)
THEN
4599 IF (
ASSOCIATED(deriv_att))
THEN
4602 virial_pw%array(:, :, :) = -rho1(:, :, :)
4603 CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
4608 IF (gradient_f)
THEN
4610 IF (my_compute_virial)
THEN
4611 CALL virial_drho_drho(virial_pw, drho, v_drho(1), virial_xc)
4621 DO ia = bo(1, 1), bo(2, 1)
4622 DO ir = bo(1, 2), bo(2, 2)
4623 vxg(idir, ia, ir, 1) = -drho(idir)%array(ia, ir, 1)*v_drho(1)%array(ia, ir, 1)
4624 IF (
ASSOCIATED(e_drho))
THEN
4625 vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + factor2*drho1(idir)%array(ia, ir, 1)*e_drho(ia, ir, 1)
4636 v_drho_r(idir, 1)%array(:, :, :) = drho(idir)%array(:, :, :)*v_drho(1)%array(:, :, :) - &
4637 drho1(idir)%array(:, :, :)*e_drho(:, :, :)
4641 CALL xc_pw_divergence(xc_deriv_method_id, v_drho_r(:, 1), tmp_g, vxc_g, v_xc(1))
4646 IF (laplace_f .AND. my_compute_virial)
THEN
4647 virial_pw%array(:, :, :) = -rho(:, :, :)
4648 CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(1)%array)
4654 DO ispin = 1, nspins
4655 CALL xc_pw_laplace(v_laplace(ispin), pw_pool, xc_deriv_method_id)
4656 CALL pw_axpy(v_laplace(ispin), v_xc(ispin))
4660 IF (gradient_f)
THEN
4662 DO ispin = 1, nspins
4663 CALL deallocate_pw(v_drho(ispin), pw_pool)
4665 CALL deallocate_pw(v_drho_r(idir, ispin), pw_pool)
4668 DEALLOCATE (v_drho, v_drho_r)
4673 DO ispin = 1, nspins
4674 CALL deallocate_pw(v_laplace(ispin), pw_pool)
4676 DEALLOCATE (v_laplace)
4679 IF (
ASSOCIATED(tmp_g%pw_grid) .AND.
ASSOCIATED(pw_pool))
THEN
4680 CALL pw_pool%give_back_pw(tmp_g)
4683 IF (
ASSOCIATED(vxc_g%pw_grid) .AND.
ASSOCIATED(pw_pool))
THEN
4684 CALL pw_pool%give_back_pw(vxc_g)
4687 IF (my_compute_virial .AND. (gradient_f .OR. laplace_f))
THEN
4688 CALL deallocate_pw(virial_pw, pw_pool)
4691 CALL timestop(handle)
4727 LOGICAL,
INTENT(in),
OPTIONAL :: spinflip
4729 CHARACTER(len=*),
PARAMETER :: routinen =
'xc_calc_3rd_deriv_analytical'
4731 INTEGER :: handle, i, idir, ispin, j, &
4732 k, nspins, xc_deriv_method_id
4733 INTEGER,
DIMENSION(2, 3) :: bo
4734 LOGICAL :: lsd, do_spinflip, alda0, &
4735 rho_f, gradient_f, tau_f, laplace_f
4736 REAL(kind=
dp) :: s, s_thresh, s_thresh2, gradient_cut
4737 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: dr1dr, dra1dra, drb1drb
4738 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: deriv_data, deriv_data2, e_drhoa, e_drhob, &
4739 e_drho, norm_drho, norm_drhoa, &
4740 norm_drhob, rho1a, rho1b, &
4742 TYPE(
cp_3d_r_cp_type),
DIMENSION(3) :: drho, drho1, drho1a, drho1b, drhoa, drhob
4743 TYPE(
pw_r3d_rs_type),
DIMENSION(:),
ALLOCATABLE :: v_drhoa, v_drhob, v_drho
4748 CALL timeset(routinen, handle)
4750 NULLIFY (e_drhoa, e_drhob, e_drho)
4752 cpassert(
ASSOCIATED(v_xc))
4753 cpassert(
ASSOCIATED(xc_section))
4757 i_val=xc_deriv_method_id)
4760 lsd =
ASSOCIATED(rho_set%rhoa)
4762 do_spinflip = .false.
4763 IF (
PRESENT(spinflip)) do_spinflip = spinflip
4765 bo = rho_set%local_bounds
4767 CALL check_for_derivatives(deriv_set, lsd, rho_f, gradient_f, tau_f, laplace_f)
4777 DO ispin = 1, nspins
4779 v_xc(ispin)%array = 0.0_dp
4783 IF (gradient_f)
THEN
4784 ALLOCATE (v_drho_r(3, nspins), v_drho(nspins))
4785 DO ispin = 1, nspins
4787 CALL allocate_pw(v_drho_r(idir, ispin), pw_pool, bo)
4789 CALL allocate_pw(v_drho(ispin), pw_pool, bo)
4793 IF (
ASSOCIATED(pw_pool))
THEN
4794 CALL pw_pool%create_pw(tmp_g)
4795 CALL pw_pool%create_pw(vxc_g)
4798 cpabort(
"XC_DERIV method is not implemented in GAPW")
4806 cpassert(
ASSOCIATED(v_xc_tau))
4807 DO ispin = 1, nspins
4808 v_xc_tau(ispin)%array = 0.0_dp
4818 IF (do_spinflip)
THEN
4825 IF (gradient_f)
THEN
4827 norm_drho=norm_drho, norm_drhoa=norm_drhoa, norm_drhob=norm_drhob)
4828 IF (do_spinflip)
THEN
4830 CALL calc_drho_from_a(drho1, drho1a)
4833 CALL calc_drho_from_ab(drho1, drho1a, drho1b)
4836 CALL calc_drho_from_ab(drho, drhoa, drhob)
4838 CALL prepare_dr1dr(dra1dra, drhoa, drho1a)
4839 IF (do_spinflip)
THEN
4840 CALL prepare_dr1dr(drb1drb, drhob, drho1a)
4841 CALL prepare_dr1dr(dr1dr, drho, drho1a)
4842 ELSE IF (nspins /= 1)
THEN
4843 CALL prepare_dr1dr(drb1drb, drhob, drho1b)
4844 CALL prepare_dr1dr(dr1dr, drho, drho1)
4846 cpabort(
"Exchange-correlation's third derivative for closed-shell not yet implemented")
4850 ALLOCATE (v_drhoa(nspins), v_drhob(nspins))
4851 DO ispin = 1, nspins
4852 CALL allocate_pw(v_drhoa(ispin), pw_pool, bo)
4853 CALL allocate_pw(v_drhob(ispin), pw_pool, bo)
4859 cpabort(
"Exchange-correlation's laplace analytic third derivative not implemented")
4863 cpabort(
"Exchange-correlation's mGGA analytic third derivative not implemented")
4866 IF (nspins /= 1)
THEN
4868 IF (.NOT. do_spinflip)
THEN
4870 cpabort(
"Exchange-correlation's analytic third derivative not implemented")
4882 IF (
ASSOCIATED(deriv_att))
THEN
4885 IF (
ASSOCIATED(deriv_att))
THEN
4889 DO k = bo(1, 3), bo(2, 3)
4890 DO j = bo(1, 2), bo(2, 2)
4891 DO i = bo(1, 1), bo(2, 1)
4892 s = rhoa(i, j, k) - rhob(i, j, k)
4893 s = -sign(max(s**2, s_thresh2), s)
4894 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4895 (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)**2/s
4907 IF (.NOT. alda0)
THEN
4909 IF (
ASSOCIATED(deriv_att))
THEN
4912 IF (
ASSOCIATED(deriv_att))
THEN
4916 DO k = bo(1, 3), bo(2, 3)
4917 DO j = bo(1, 2), bo(2, 2)
4918 DO i = bo(1, 1), bo(2, 1)
4919 s = rhoa(i, j, k) - rhob(i, j, k)
4920 s = -sign(max(s**2, s_thresh2), s)
4921 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4922 (deriv_data(i, j, k)*dra1dra(i, j, k) - deriv_data2(i, j, k)*drb1drb(i, j, k)) &
4934 DO k = bo(1, 3), bo(2, 3)
4935 DO j = bo(1, 2), bo(2, 2)
4936 DO i = bo(1, 1), bo(2, 1)
4937 v_xc(2)%array(i, j, k) = -v_xc(1)%array(i, j, k)
4950 IF (
ASSOCIATED(deriv_att))
THEN
4953 IF (
ASSOCIATED(deriv_att))
THEN
4957 DO k = bo(1, 3), bo(2, 3)
4958 DO j = bo(1, 2), bo(2, 2)
4959 DO i = bo(1, 1), bo(2, 1)
4960 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
4961 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4962 rho1a(i, j, k)**2*(deriv_data(i, j, k) - deriv_data2(i, j, k))/s
4974 IF (
ASSOCIATED(deriv_att))
THEN
4977 IF (
ASSOCIATED(deriv_att))
THEN
4981 DO k = bo(1, 3), bo(2, 3)
4982 DO j = bo(1, 2), bo(2, 2)
4983 DO i = bo(1, 1), bo(2, 1)
4984 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
4985 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
4986 rho1a(i, j, k)**2*(deriv_data(i, j, k) - deriv_data2(i, j, k))/s
4994 IF (.NOT. alda0)
THEN
4999 IF (
ASSOCIATED(deriv_att))
THEN
5002 IF (
ASSOCIATED(deriv_att))
THEN
5006 DO k = bo(1, 3), bo(2, 3)
5007 DO j = bo(1, 2), bo(2, 2)
5008 DO i = bo(1, 1), bo(2, 1)
5009 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
5010 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
5011 (deriv_data(i, j, k)*dra1dra(i, j, k) - deriv_data2(i, j, k)*drb1drb(i, j, k))* &
5024 IF (
ASSOCIATED(deriv_att))
THEN
5027 IF (
ASSOCIATED(deriv_att))
THEN
5031 DO k = bo(1, 3), bo(2, 3)
5032 DO j = bo(1, 2), bo(2, 2)
5033 DO i = bo(1, 1), bo(2, 1)
5034 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
5035 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
5036 (deriv_data(i, j, k)*dra1dra(i, j, k) - deriv_data2(i, j, k)*drb1drb(i, j, k))* &
5052 IF (
ASSOCIATED(deriv_att))
THEN
5055 IF (
ASSOCIATED(deriv_att))
THEN
5059 DO k = bo(1, 3), bo(2, 3)
5060 DO j = bo(1, 2), bo(2, 2)
5061 DO i = bo(1, 1), bo(2, 1)
5062 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
5063 (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)
5075 IF (
ASSOCIATED(deriv_att))
THEN
5078 IF (
ASSOCIATED(deriv_att))
THEN
5082 DO k = bo(1, 3), bo(2, 3)
5083 DO j = bo(1, 2), bo(2, 2)
5084 DO i = bo(1, 1), bo(2, 1)
5085 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) + &
5086 (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)
5102 IF (
ASSOCIATED(deriv_att))
THEN
5105 IF (
ASSOCIATED(deriv_att))
THEN
5109 DO k = bo(1, 3), bo(2, 3)
5110 DO j = bo(1, 2), bo(2, 2)
5111 DO i = bo(1, 1), bo(2, 1)
5112 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
5113 (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)
5114 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
5115 (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)
5124 IF (
ASSOCIATED(deriv_att))
THEN
5130 IF (
ASSOCIATED(deriv_att))
THEN
5134 DO k = bo(1, 3), bo(2, 3)
5135 DO j = bo(1, 2), bo(2, 2)
5136 DO i = bo(1, 1), bo(2, 1)
5137 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
5138 deriv_data(i, j, k)*dra1dra(i, j, k) + deriv_data2(i, j, k)*drb1drb(i, j, k)
5148 IF (
ASSOCIATED(deriv_att))
THEN
5152 DO k = bo(1, 3), bo(2, 3)
5153 DO j = bo(1, 2), bo(2, 2)
5154 DO i = bo(1, 1), bo(2, 1)
5155 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
5156 deriv_data(i, j, k)*dra1dra(i, j, k) + deriv_data2(i, j, k)*drb1drb(i, j, k)
5171 IF (
ASSOCIATED(deriv_att))
THEN
5174 IF (
ASSOCIATED(deriv_att))
THEN
5178 DO k = bo(1, 3), bo(2, 3)
5179 DO j = bo(1, 2), bo(2, 2)
5180 DO i = bo(1, 1), bo(2, 1)
5181 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
5182 deriv_data(i, j, k)*dra1dra(i, j, k) + &
5183 deriv_data2(i, j, k)*drb1drb(i, j, k)
5184 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
5185 deriv_data(i, j, k)*dra1dra(i, j, k) + &
5186 deriv_data2(i, j, k)*drb1drb(i, j, k)
5201 IF (
ASSOCIATED(deriv_att))
THEN
5206 v_drhoa(1)%array(:, :, :) = v_drhoa(1)%array(:, :, :) + &
5207 deriv_data(:, :, :)*dra1dra(:, :, :)/max(gradient_cut, norm_drhoa(:, :, :))**2
5215 IF (
ASSOCIATED(deriv_att))
THEN
5220 v_drhob(2)%array(:, :, :) = v_drhob(2)%array(:, :, :) - &
5221 deriv_data(:, :, :)*drb1drb(:, :, :)/max(gradient_cut, norm_drhob(:, :, :))**2
5230 cpabort(
"Exchange-correlation's analytic third derivative not implemented")
5234 IF (gradient_f)
THEN
5235 IF (.NOT. alda0)
THEN
5249 IF (do_spinflip)
THEN
5252 DO k = bo(1, 3), bo(2, 3)
5253 DO j = bo(1, 2), bo(2, 2)
5254 DO i = bo(1, 1), bo(2, 1)
5255 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
5257 v_drho_r(idir, ispin)%array(i, j, k) = v_drho_r(idir, ispin)%array(i, j, k) + &
5258 v_drhoa(ispin)%array(i, j, k)*drhoa(idir)%array(i, j, k)*rho1a(i, j, k)/s + &
5259 v_drhob(ispin)%array(i, j, k)*drhob(idir)%array(i, j, k)*rho1a(i, j, k)/s + &
5260 v_drho(ispin)%array(i, j, k)*drho(idir)%array(i, j, k)*rho1a(i, j, k)/s
5272 IF (
ASSOCIATED(e_drhoa))
THEN
5275 DO k = bo(1, 3), bo(2, 3)
5276 DO j = bo(1, 2), bo(2, 2)
5277 DO i = bo(1, 1), bo(2, 1)
5278 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
5279 v_drho_r(idir, 1)%array(i, j, k) = v_drho_r(idir, 1)%array(i, j, k) - &
5280 e_drhoa(i, j, k)*drho1a(idir)%array(i, j, k)*rho1a(i, j, k)/s
5290 IF (
ASSOCIATED(e_drhob))
THEN
5293 DO k = bo(1, 3), bo(2, 3)
5294 DO j = bo(1, 2), bo(2, 2)
5295 DO i = bo(1, 1), bo(2, 1)
5296 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
5297 v_drho_r(idir, 2)%array(i, j, k) = v_drho_r(idir, 2)%array(i, j, k) + &
5298 e_drhob(i, j, k)*drho1a(idir)%array(i, j, k)*rho1a(i, j, k)/s
5307 DO ispin = 1, nspins
5308 CALL xc_pw_divergence(xc_deriv_method_id, v_drho_r(:, ispin), tmp_g, vxc_g, v_xc(ispin))
5313 DEALLOCATE (drho(idir)%array)
5314 DEALLOCATE (drho1(idir)%array)
5317 DO ispin = 1, nspins
5318 CALL deallocate_pw(v_drhoa(ispin), pw_pool)
5319 CALL deallocate_pw(v_drhob(ispin), pw_pool)
5322 DEALLOCATE (v_drhoa, v_drhob)
5331 cpabort(
"Exchange-correlation's analytic third derivative not implemented")
5335 IF (gradient_f)
THEN
5337 DO ispin = 1, nspins
5338 CALL deallocate_pw(v_drho(ispin), pw_pool)
5340 CALL deallocate_pw(v_drho_r(idir, ispin), pw_pool)
5343 DEALLOCATE (v_drho, v_drho_r)
5347 IF (
ASSOCIATED(tmp_g%pw_grid) .AND.
ASSOCIATED(pw_pool))
THEN
5348 CALL pw_pool%give_back_pw(tmp_g)
5351 IF (
ASSOCIATED(vxc_g%pw_grid) .AND.
ASSOCIATED(pw_pool))
THEN
5352 CALL pw_pool%give_back_pw(vxc_g)
5355 CALL timestop(handle)