1661 pw_pool, xc_section, &
1662 do_triplet, calc_virial, virial_xc, deriv_set)
1664 TYPE(
pw_r3d_rs_type),
DIMENSION(:),
INTENT(IN),
POINTER :: v_xc, v_tau
1666 TYPE(
pw_r3d_rs_type),
DIMENSION(:),
INTENT(IN),
POINTER :: rho1_r, tau1_r
1667 TYPE(
pw_c1d_gs_type),
DIMENSION(:),
INTENT(IN),
POINTER :: rho1_g
1670 LOGICAL,
INTENT(IN) :: do_triplet
1671 LOGICAL,
INTENT(IN),
OPTIONAL :: calc_virial
1672 REAL(kind=
dp),
DIMENSION(3, 3),
INTENT(INOUT), &
1673 OPTIONAL :: virial_xc
1676 CHARACTER(len=*),
PARAMETER :: routinen =
'xc_calc_2nd_deriv_numerical'
1677 REAL(kind=
dp),
DIMENSION(-4:4, 4),
PARAMETER :: &
1678 weights = 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, &
1679 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, &
1680 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, &
1681 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])
1683 INTEGER :: handle, idir, ispin, nspins, istep, nsteps
1684 INTEGER,
DIMENSION(2, 3) :: bo
1685 LOGICAL :: gradient_f, lsd, my_calc_virial, tau_f, laplace_f, rho_f
1686 REAL(kind=
dp) :: exc, gradient_cut, h, weight, step, rho_cutoff
1687 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: dr1dr, dra1dra, drb1drb
1688 REAL(kind=
dp),
DIMENSION(3, 3) :: virial_dummy
1689 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: norm_drho, norm_drho2, norm_drho2a, &
1690 norm_drho2b, norm_drhoa, norm_drhob, &
1691 rho, rho1, rho1a, rho1b, rhoa, rhob, &
1692 tau_a, tau_b, tau, tau1, tau1a, tau1b, laplace, laplace1, &
1693 laplacea, laplaceb, laplace1a, laplace1b, &
1694 laplace2, laplace2a, laplace2b, deriv_data
1695 TYPE(
cp_3d_r_cp_type),
DIMENSION(3) :: drho, drho1, drho1a, drho1b, drhoa, drhob
1700 TYPE(
pw_r3d_rs_type) :: virial_pw, v_laplace, v_laplacea, v_laplaceb
1706 CALL timeset(routinen, handle)
1708 my_calc_virial = .false.
1709 IF (
PRESENT(calc_virial) .AND.
PRESENT(virial_xc)) my_calc_virial = calc_virial
1713 NULLIFY (tau, tau_r, tau_a, tau_b)
1717 IF (nsteps < lbound(weights, 2) .OR. nspins > ubound(weights, 2))
THEN
1718 cpabort(
"The number of steps must be a value from 1 to 4.")
1721 IF (nspins == 2)
THEN
1722 NULLIFY (vxc_rho, rho_g, vxc_tau)
1724 DO ispin = 1, nspins
1725 CALL pw_pool%create_pw(rho_r(ispin))
1727 IF (
ASSOCIATED(tau1_r) .AND.
ASSOCIATED(v_tau))
THEN
1729 DO ispin = 1, nspins
1730 CALL pw_pool%create_pw(tau_r(ispin))
1733 CALL xc_rho_set_get(rho_set, can_return_null=.true., rhoa=rhoa, rhob=rhob, tau_a=tau_a, tau_b=tau_b)
1734 DO istep = -nsteps, nsteps
1735 IF (istep == 0) cycle
1736 weight = weights(istep, nsteps)/h
1737 step = real(istep,
dp)*h
1738 CALL calc_resp_potential_numer_ab(rho_r, rho_g, rho1_r, rhoa, rhob, vxc_rho, &
1739 tau_r, tau1_r, tau_a, tau_b, vxc_tau, xc_section, pw_pool, step)
1740 DO ispin = 1, nspins
1741 CALL pw_axpy(vxc_rho(ispin), v_xc(ispin), weight)
1742 IF (
ASSOCIATED(vxc_tau) .AND.
ASSOCIATED(v_tau))
THEN
1743 CALL pw_axpy(vxc_tau(ispin), v_tau(ispin), weight)
1746 DO ispin = 1, nspins
1747 CALL vxc_rho(ispin)%release()
1749 DEALLOCATE (vxc_rho)
1750 IF (
ASSOCIATED(vxc_tau))
THEN
1751 DO ispin = 1, nspins
1752 CALL vxc_tau(ispin)%release()
1754 DEALLOCATE (vxc_tau)
1757 ELSE IF (nspins == 1 .AND. do_triplet)
THEN
1758 NULLIFY (vxc_rho, vxc_tau, rho_g)
1761 CALL pw_pool%create_pw(rho_r(ispin))
1763 IF (
ASSOCIATED(tau1_r) .AND.
ASSOCIATED(v_tau))
THEN
1765 DO ispin = 1, nspins
1766 CALL pw_pool%create_pw(tau_r(ispin))
1769 CALL xc_rho_set_get(rho_set, can_return_null=.true., rhoa=rhoa, rhob=rhob, tau_a=tau_a, tau_b=tau_b)
1770 DO istep = -nsteps, nsteps
1771 IF (istep == 0) cycle
1772 weight = weights(istep, nsteps)/h
1773 step = real(istep,
dp)*h
1777 rho_r(1)%array(:, :, :) = rhoa(:, :, :) + step*rho1_r(1)%array(:, :, :)
1780 rho_r(2)%array(:, :, :) = rhob(:, :, :)
1782 IF (
ASSOCIATED(tau1_r))
THEN
1784 tau_r(1)%array(:, :, :) = tau_a(:, :, :) + step*tau1_r(1)%array(:, :, :)
1787 tau_r(2)%array(:, :, :) = tau_b(:, :, :)
1791 CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
1792 pw_pool, .false., virial_dummy)
1793 CALL pw_axpy(vxc_rho(1), v_xc(1), weight)
1794 IF (
ASSOCIATED(vxc_tau) .AND.
ASSOCIATED(v_tau))
THEN
1795 CALL pw_axpy(vxc_tau(1), v_tau(1), weight)
1798 CALL vxc_rho(ispin)%release()
1800 DEALLOCATE (vxc_rho)
1801 IF (
ASSOCIATED(vxc_tau))
THEN
1803 CALL vxc_tau(ispin)%release()
1805 DEALLOCATE (vxc_tau)
1810 rho_r(1)%array(:, :, :) = rhoa(:, :, :)
1813 rho_r(2)%array(:, :, :) = rhob(:, :, :) + step*rho1_r(1)%array(:, :, :)
1815 IF (
ASSOCIATED(tau1_r))
THEN
1817 tau_r(1)%array(:, :, :) = tau_a(:, :, :)
1820 tau_r(2)%array(:, :, :) = tau_b(:, :, :) + step*tau1_r(1)%array(:, :, :)
1824 CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
1825 pw_pool, .false., virial_dummy)
1826 CALL pw_axpy(vxc_rho(1), v_xc(1), weight)
1827 IF (
ASSOCIATED(vxc_tau) .AND.
ASSOCIATED(v_tau))
THEN
1828 CALL pw_axpy(vxc_tau(1), v_tau(1), weight)
1831 CALL vxc_rho(ispin)%release()
1833 DEALLOCATE (vxc_rho)
1834 IF (
ASSOCIATED(vxc_tau))
THEN
1836 CALL vxc_tau(ispin)%release()
1838 DEALLOCATE (vxc_tau)
1842 NULLIFY (vxc_rho, rho_r, rho_g, vxc_tau, tau_r, tau)
1844 CALL pw_pool%create_pw(rho_r(1))
1845 IF (
ASSOCIATED(tau1_r) .AND.
ASSOCIATED(v_tau))
THEN
1847 CALL pw_pool%create_pw(tau_r(1))
1849 CALL xc_rho_set_get(rho_set, can_return_null=.true., rho=rho, tau=tau)
1850 DO istep = -nsteps, nsteps
1851 IF (istep == 0) cycle
1852 weight = weights(istep, nsteps)/h
1853 step = real(istep,
dp)*h
1856 rho_r(1)%array(:, :, :) = rho(:, :, :) + step*rho1_r(1)%array(:, :, :)
1858 IF (
ASSOCIATED(tau1_r) .AND.
ASSOCIATED(tau) .AND.
ASSOCIATED(tau1_r))
THEN
1860 tau_r(1)%array(:, :, :) = tau(:, :, :) + step*tau1_r(1)%array(:, :, :)
1864 CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
1865 pw_pool, .false., virial_dummy)
1866 CALL pw_axpy(vxc_rho(1), v_xc(1), weight)
1867 IF (
ASSOCIATED(vxc_tau) .AND.
ASSOCIATED(v_tau))
THEN
1868 CALL pw_axpy(vxc_tau(1), v_tau(1), weight)
1870 CALL vxc_rho(1)%release()
1871 DEALLOCATE (vxc_rho)
1872 IF (
ASSOCIATED(vxc_tau))
THEN
1873 CALL vxc_tau(1)%release()
1874 DEALLOCATE (vxc_tau)
1879 IF (my_calc_virial)
THEN
1881 IF (nspins == 1 .AND. do_triplet)
THEN
1885 CALL check_for_derivatives(deriv_set, (nspins == 2), rho_f, gradient_f, tau_f, laplace_f)
1891 IF (gradient_f)
THEN
1892 bo = rho_set%local_bounds
1895 CALL allocate_pw(virial_pw, pw_pool, bo)
1902 drho_cutoff=gradient_cut, &
1915 CALL xc_rho_set_get(rho_set, drhoa=drhoa, drhob=drhob, norm_drho=norm_drho, &
1916 norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, tau_a=tau_a, tau_b=tau_b, &
1917 laplace_rhoa=laplacea, laplace_rhob=laplaceb, can_return_null=.true.)
1918 CALL xc_rho_set_get(rho1_set, rhoa=rho1a, rhob=rho1b, drhoa=drho1a, drhob=drho1b, laplace_rhoa=laplace1a, &
1919 laplace_rhob=laplace1b, can_return_null=.true.)
1921 CALL calc_drho_from_ab(drho, drhoa, drhob)
1922 CALL calc_drho_from_ab(drho1, drho1a, drho1b)
1924 CALL xc_rho_set_get(rho_set, drho=drho, norm_drho=norm_drho, tau=tau, laplace_rho=laplace, can_return_null=.true.)
1925 CALL xc_rho_set_get(rho1_set, rho=rho1, drho=drho1, laplace_rho=laplace1, can_return_null=.true.)
1928 CALL prepare_dr1dr(dr1dr, drho, drho1)
1931 CALL prepare_dr1dr(dra1dra, drhoa, drho1a)
1932 CALL prepare_dr1dr(drb1drb, drhob, drho1b)
1934 CALL allocate_pw(v_drho, pw_pool, bo)
1935 CALL allocate_pw(v_drhoa, pw_pool, bo)
1936 CALL allocate_pw(v_drhob, pw_pool, bo)
1938 IF (
ASSOCIATED(norm_drhoa))
CALL apply_drho(deriv_set, [
deriv_norm_drhoa], virial_pw, drhoa, drho1a, virial_xc, &
1939 norm_drhoa, gradient_cut, dra1dra, v_drhoa%array)
1940 IF (
ASSOCIATED(norm_drhob))
CALL apply_drho(deriv_set, [
deriv_norm_drhob], virial_pw, drhob, drho1b, virial_xc, &
1941 norm_drhob, gradient_cut, drb1drb, v_drhob%array)
1942 IF (
ASSOCIATED(norm_drho))
CALL apply_drho(deriv_set, [
deriv_norm_drho], virial_pw, drho, drho1, virial_xc, &
1943 norm_drho, gradient_cut, dr1dr, v_drho%array)
1946 cpassert(
ASSOCIATED(deriv_data))
1947 virial_pw%array(:, :, :) = -rho1a(:, :, :)
1948 CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
1950 CALL allocate_pw(v_laplacea, pw_pool, bo)
1953 cpassert(
ASSOCIATED(deriv_data))
1954 virial_pw%array(:, :, :) = -rho1b(:, :, :)
1955 CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
1957 CALL allocate_pw(v_laplaceb, pw_pool, bo)
1963 CALL allocate_pw(v_drho, pw_pool, bo)
1965 CALL apply_drho(deriv_set, [
deriv_norm_drho], virial_pw, drho, drho1, virial_xc, &
1966 norm_drho, gradient_cut, dr1dr, v_drho%array)
1969 cpassert(
ASSOCIATED(deriv_data))
1970 virial_pw%array(:, :, :) = -rho1(:, :, :)
1971 CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
1973 CALL allocate_pw(v_laplace, pw_pool, bo)
1979 rho_r(1)%array = rhoa
1980 rho_r(2)%array = rhob
1982 rho_r(1)%array = rho
1984 IF (
ASSOCIATED(tau1_r))
THEN
1986 tau_r(1)%array = tau_a
1987 tau_r(2)%array = tau_b
1989 tau_r(1)%array = tau
2000 rho_cutoff=rho_cutoff, &
2011 CALL xc_rho_set_get(rho1_set, rhoa=rho1a, rhob=rho1b, tau_a=tau1a, tau_b=tau1b, &
2012 laplace_rhoa=laplace1a, laplace_rhob=laplace1b, can_return_null=.true.)
2013 CALL xc_rho_set_get(rho2_set, norm_drhoa=norm_drho2a, norm_drhob=norm_drho2b, &
2014 norm_drho=norm_drho2, laplace_rhoa=laplace2a, laplace_rhob=laplace2b, can_return_null=.true.)
2016 DO istep = -nsteps, nsteps
2017 IF (istep == 0) cycle
2018 weight = weights(istep, nsteps)/h
2019 step = real(istep,
dp)*h
2020 IF (
ASSOCIATED(norm_drhoa))
THEN
2021 CALL get_derivs_rho(norm_drho2a, norm_drhoa, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2022 CALL update_deriv_rho(deriv_set1, [
deriv_rhoa], bo, &
2023 norm_drhoa, gradient_cut, weight, rho1a, v_drhoa%array)
2024 CALL update_deriv_rho(deriv_set1, [
deriv_rhob], bo, &
2025 norm_drhoa, gradient_cut, weight, rho1b, v_drhoa%array)
2027 norm_drhoa, gradient_cut, weight, dra1dra, v_drhoa%array)
2029 norm_drhoa, gradient_cut, weight, dra1dra, drb1drb, v_drhoa%array, v_drhob%array)
2031 norm_drhoa, gradient_cut, weight, dra1dra, dr1dr, v_drhoa%array, v_drho%array)
2033 CALL update_deriv_rho(deriv_set1, [
deriv_tau_a], bo, &
2034 norm_drhoa, gradient_cut, weight, tau1a, v_drhoa%array)
2035 CALL update_deriv_rho(deriv_set1, [
deriv_tau_b], bo, &
2036 norm_drhoa, gradient_cut, weight, tau1b, v_drhoa%array)
2040 norm_drhoa, gradient_cut, weight, laplace1a, v_drhoa%array)
2042 norm_drhoa, gradient_cut, weight, laplace1b, v_drhoa%array)
2046 IF (
ASSOCIATED(norm_drhob))
THEN
2047 CALL get_derivs_rho(norm_drho2b, norm_drhob, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2048 CALL update_deriv_rho(deriv_set1, [
deriv_rhoa], bo, &
2049 norm_drhob, gradient_cut, weight, rho1a, v_drhob%array)
2050 CALL update_deriv_rho(deriv_set1, [
deriv_rhob], bo, &
2051 norm_drhob, gradient_cut, weight, rho1b, v_drhob%array)
2053 norm_drhob, gradient_cut, weight, drb1drb, v_drhob%array)
2055 norm_drhob, gradient_cut, weight, drb1drb, dra1dra, v_drhob%array, v_drhoa%array)
2057 norm_drhob, gradient_cut, weight, drb1drb, dr1dr, v_drhob%array, v_drho%array)
2059 CALL update_deriv_rho(deriv_set1, [
deriv_tau_a], bo, &
2060 norm_drhob, gradient_cut, weight, tau1a, v_drhob%array)
2061 CALL update_deriv_rho(deriv_set1, [
deriv_tau_b], bo, &
2062 norm_drhob, gradient_cut, weight, tau1b, v_drhob%array)
2066 norm_drhob, gradient_cut, weight, laplace1a, v_drhob%array)
2068 norm_drhob, gradient_cut, weight, laplace1b, v_drhob%array)
2072 IF (
ASSOCIATED(norm_drho))
THEN
2073 CALL get_derivs_rho(norm_drho2, norm_drho, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2074 CALL update_deriv_rho(deriv_set1, [
deriv_rhoa], bo, &
2075 norm_drho, gradient_cut, weight, rho1a, v_drho%array)
2076 CALL update_deriv_rho(deriv_set1, [
deriv_rhob], bo, &
2077 norm_drho, gradient_cut, weight, rho1b, v_drho%array)
2079 norm_drho, gradient_cut, weight, dr1dr, v_drho%array)
2081 norm_drho, gradient_cut, weight, dr1dr, dra1dra, v_drho%array, v_drhoa%array)
2083 norm_drho, gradient_cut, weight, dr1dr, drb1drb, v_drho%array, v_drhob%array)
2085 CALL update_deriv_rho(deriv_set1, [
deriv_tau_a], bo, &
2086 norm_drho, gradient_cut, weight, tau1a, v_drho%array)
2087 CALL update_deriv_rho(deriv_set1, [
deriv_tau_b], bo, &
2088 norm_drho, gradient_cut, weight, tau1b, v_drho%array)
2092 norm_drho, gradient_cut, weight, laplace1a, v_drho%array)
2094 norm_drho, gradient_cut, weight, laplace1b, v_drho%array)
2100 CALL get_derivs_rho(laplace2a, laplacea, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2103 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_rhoa], bo, &
2104 weight, rho1a, v_laplacea%array)
2105 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_rhob], bo, &
2106 weight, rho1b, v_laplacea%array)
2107 IF (
ASSOCIATED(norm_drho))
THEN
2108 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_norm_drho], bo, &
2109 weight, dr1dr, v_laplacea%array)
2111 IF (
ASSOCIATED(norm_drhoa))
THEN
2112 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_norm_drhoa], bo, &
2113 weight, dra1dra, v_laplacea%array)
2115 IF (
ASSOCIATED(norm_drhob))
THEN
2116 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_norm_drhob], bo, &
2117 weight, drb1drb, v_laplacea%array)
2120 IF (
ASSOCIATED(tau1a))
THEN
2121 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_tau_a], bo, &
2122 weight, tau1a, v_laplacea%array)
2124 IF (
ASSOCIATED(tau1b))
THEN
2125 CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [
deriv_tau_b], bo, &
2126 weight, tau1b, v_laplacea%array)
2130 weight, laplace1a, v_laplacea%array)
2133 weight, laplace1b, v_laplacea%array)
2136 CALL get_derivs_rho(laplace2b, laplaceb, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2139 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_rhoa], bo, &
2140 weight, rho1a, v_laplaceb%array)
2141 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_rhob], bo, &
2142 weight, rho1b, v_laplaceb%array)
2143 IF (
ASSOCIATED(norm_drho))
THEN
2144 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_norm_drho], bo, &
2145 weight, dr1dr, v_laplaceb%array)
2147 IF (
ASSOCIATED(norm_drhoa))
THEN
2148 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_norm_drhoa], bo, &
2149 weight, dra1dra, v_laplaceb%array)
2151 IF (
ASSOCIATED(norm_drhob))
THEN
2152 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_norm_drhob], bo, &
2153 weight, drb1drb, v_laplaceb%array)
2157 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_tau_a], bo, &
2158 weight, tau1a, v_laplaceb%array)
2159 CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [
deriv_tau_b], bo, &
2160 weight, tau1b, v_laplaceb%array)
2164 weight, laplace1a, v_laplaceb%array)
2167 weight, laplace1b, v_laplaceb%array)
2171 CALL virial_drho_drho(virial_pw, drhoa, v_drhoa, virial_xc)
2172 CALL virial_drho_drho(virial_pw, drhob, v_drhob, virial_xc)
2173 CALL virial_drho_drho(virial_pw, drho, v_drho, virial_xc)
2175 CALL deallocate_pw(v_drho, pw_pool)
2176 CALL deallocate_pw(v_drhoa, pw_pool)
2177 CALL deallocate_pw(v_drhob, pw_pool)
2180 virial_pw%array(:, :, :) = -rhoa(:, :, :)
2181 CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplacea%array)
2182 CALL deallocate_pw(v_laplacea, pw_pool)
2184 virial_pw%array(:, :, :) = -rhob(:, :, :)
2185 CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplaceb%array)
2186 CALL deallocate_pw(v_laplaceb, pw_pool)
2189 CALL deallocate_pw(virial_pw, pw_pool)
2192 DEALLOCATE (drho(idir)%array)
2193 DEALLOCATE (drho1(idir)%array)
2195 DEALLOCATE (dra1dra, drb1drb)
2198 CALL xc_rho_set_get(rho1_set, rho=rho1, tau=tau1, laplace_rho=laplace1, can_return_null=.true.)
2199 CALL xc_rho_set_get(rho2_set, norm_drho=norm_drho2, laplace_rho=laplace2, can_return_null=.true.)
2201 DO istep = -nsteps, nsteps
2202 IF (istep == 0) cycle
2203 weight = weights(istep, nsteps)/h
2204 step = real(istep,
dp)*h
2205 CALL get_derivs_rho(norm_drho2, norm_drho, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2208 CALL update_deriv_rho(deriv_set1, [
deriv_rho], bo, &
2209 norm_drho, gradient_cut, weight, rho1, v_drho%array)
2211 norm_drho, gradient_cut, weight, dr1dr, v_drho%array)
2214 CALL update_deriv_rho(deriv_set1, [
deriv_tau], bo, &
2215 norm_drho, gradient_cut, weight, tau1, v_drho%array)
2219 norm_drho, gradient_cut, weight, laplace1, v_drho%array)
2221 CALL get_derivs_rho(laplace2, laplace, step, xc_fun_section, lsd, rho2_set, deriv_set1)
2224 CALL update_deriv(deriv_set1, laplace, rho_cutoff, [
deriv_rho], bo, &
2225 weight, rho1, v_laplace%array)
2226 CALL update_deriv(deriv_set1, laplace, rho_cutoff, [
deriv_norm_drho], bo, &
2227 weight, dr1dr, v_laplace%array)
2230 CALL update_deriv(deriv_set1, laplace, rho_cutoff, [
deriv_tau], bo, &
2231 weight, tau1, v_laplace%array)
2235 weight, laplace1, v_laplace%array)
2240 CALL virial_drho_drho(virial_pw, drho, v_drho, virial_xc)
2242 CALL deallocate_pw(v_drho, pw_pool)
2245 virial_pw%array(:, :, :) = -rho(:, :, :)
2246 CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace%array)
2247 CALL deallocate_pw(v_laplace, pw_pool)
2250 CALL deallocate_pw(virial_pw, pw_pool)
2263 DO ispin = 1,
SIZE(rho_r)
2264 CALL pw_pool%give_back_pw(rho_r(ispin))
2268 IF (
ASSOCIATED(tau_r))
THEN
2269 DO ispin = 1,
SIZE(tau_r)
2270 CALL pw_pool%give_back_pw(tau_r(ispin))
2275 CALL timestop(handle)
2648 LOGICAL,
INTENT(IN),
OPTIONAL :: gapw
2649 REAL(kind=
dp),
DIMENSION(:, :, :, :),
OPTIONAL, &
2651 REAL(kind=
dp),
INTENT(in),
OPTIONAL :: tddfpt_fac
2652 LOGICAL,
INTENT(IN),
OPTIONAL :: compute_virial
2653 REAL(kind=
dp),
DIMENSION(3, 3),
INTENT(INOUT), &
2654 OPTIONAL :: virial_xc
2655 LOGICAL,
INTENT(in),
OPTIONAL :: spinflip
2657 CHARACTER(len=*),
PARAMETER :: routinen =
'xc_calc_2nd_deriv_analytical'
2659 INTEGER :: handle, i, ia, idir, ir, ispin, j, jdir, &
2660 k, nspins, xc_deriv_method_id
2661 INTEGER,
DIMENSION(2, 3) :: bo
2662 LOGICAL :: gradient_f, lsd, my_compute_virial, alda0, &
2663 my_gapw, tau_f, laplace_f, rho_f, do_spinflip
2664 REAL(kind=
dp) ::
fac, gradient_cut, tmp, factor2, s, s_thresh
2665 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: dr1dr, dra1dra, drb1drb
2666 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: deriv_data, deriv_data2, &
2667 e_drhoa, e_drhob, e_drho, norm_drho, norm_drhoa, &
2668 norm_drhob, rho1, rho1a, rho1b, &
2669 tau1, tau1a, tau1b, laplace1, laplace1a, laplace1b, &
2671 TYPE(
cp_3d_r_cp_type),
DIMENSION(3) :: drho, drho1, drho1a, drho1b, drhoa, drhob
2672 TYPE(
pw_r3d_rs_type),
DIMENSION(:),
ALLOCATABLE :: v_drhoa, v_drhob, v_drho, v_laplace
2678 CALL timeset(routinen, handle)
2680 NULLIFY (e_drhoa, e_drhob, e_drho)
2683 IF (
PRESENT(gapw)) my_gapw = gapw
2685 my_compute_virial = .false.
2686 IF (
PRESENT(compute_virial)) my_compute_virial = compute_virial
2688 cpassert(
ASSOCIATED(v_xc))
2689 cpassert(
ASSOCIATED(xc_section))
2691 cpassert(
PRESENT(vxg))
2693 IF (my_compute_virial)
THEN
2694 cpassert(
PRESENT(virial_xc))
2698 i_val=xc_deriv_method_id)
2701 lsd =
ASSOCIATED(rho_set%rhoa)
2704 IF (
PRESENT(tddfpt_fac))
fac = tddfpt_fac
2705 IF (
PRESENT(tddfpt_fac)) factor2 = tddfpt_fac
2706 do_spinflip = .false.
2707 IF (
PRESENT(spinflip)) do_spinflip = spinflip
2709 bo = rho_set%local_bounds
2711 CALL check_for_derivatives(deriv_set, lsd, rho_f, gradient_f, tau_f, laplace_f)
2714 IF (gradient_f)
THEN
2721 cpassert(
ASSOCIATED(v_xc_tau))
2724 IF (gradient_f)
THEN
2725 ALLOCATE (v_drho_r(3, nspins), v_drho(nspins))
2726 DO ispin = 1, nspins
2728 CALL allocate_pw(v_drho_r(idir, ispin), pw_pool, bo)
2730 CALL allocate_pw(v_drho(ispin), pw_pool, bo)
2734 IF (
ASSOCIATED(pw_pool))
THEN
2735 CALL pw_pool%create_pw(tmp_g)
2736 CALL pw_pool%create_pw(vxc_g)
2739 cpabort(
"XC_DERIV method is not implemented in GAPW")
2744 DO ispin = 1, nspins
2745 v_xc(ispin)%array = 0.0_dp
2749 DO ispin = 1, nspins
2750 v_xc_tau(ispin)%array = 0.0_dp
2754 IF (laplace_f .AND. my_gapw) &
2755 cpabort(
"Laplace-dependent functional not implemented with GAPW!")
2757 IF (my_compute_virial .AND. (gradient_f .OR. laplace_f))
CALL allocate_pw(virial_pw, pw_pool, bo)
2765 IF (do_spinflip)
THEN
2772 IF (gradient_f)
THEN
2774 norm_drho=norm_drho, norm_drhoa=norm_drhoa, norm_drhob=norm_drhob)
2775 IF (do_spinflip)
THEN
2777 CALL calc_drho_from_a(drho1, drho1a)
2780 CALL calc_drho_from_ab(drho1, drho1a, drho1b)
2783 CALL calc_drho_from_ab(drho, drhoa, drhob)
2785 CALL prepare_dr1dr(dra1dra, drhoa, drho1a)
2786 IF (do_spinflip)
THEN
2787 CALL prepare_dr1dr(drb1drb, drhob, drho1a)
2788 CALL prepare_dr1dr(dr1dr, drho, drho1a)
2789 ELSE IF (nspins /= 1)
THEN
2790 CALL prepare_dr1dr(drb1drb, drhob, drho1b)
2791 CALL prepare_dr1dr(dr1dr, drho, drho1)
2793 CALL prepare_dr1dr(drb1drb, drhob, drho1b)
2794 CALL prepare_dr1dr_ab(dr1dr, drhoa, drhob, drho1a, drho1b,
fac)
2797 ALLOCATE (v_drhoa(nspins), v_drhob(nspins))
2798 DO ispin = 1, nspins
2799 CALL allocate_pw(v_drhoa(ispin), pw_pool, bo)
2800 CALL allocate_pw(v_drhob(ispin), pw_pool, bo)
2806 CALL xc_rho_set_get(rho1_set, laplace_rhoa=laplace1a, laplace_rhob=laplace1b)
2808 ALLOCATE (v_laplace(nspins))
2809 DO ispin = 1, nspins
2810 CALL allocate_pw(v_laplace(ispin), pw_pool, bo)
2813 IF (my_compute_virial)
CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob)
2820 IF (do_spinflip)
THEN
2829 IF (
ASSOCIATED(deriv_att))
THEN
2832 IF (
ASSOCIATED(deriv_att))
THEN
2836 DO k = bo(1, 3), bo(2, 3)
2837 DO j = bo(1, 2), bo(2, 2)
2838 DO i = bo(1, 1), bo(2, 1)
2839 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
2840 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2841 (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)/s
2853 IF (.NOT. alda0)
THEN
2855 IF (
ASSOCIATED(deriv_att))
THEN
2858 IF (
ASSOCIATED(deriv_att))
THEN
2862 DO k = bo(1, 3), bo(2, 3)
2863 DO j = bo(1, 2), bo(2, 2)
2864 DO i = bo(1, 1), bo(2, 1)
2865 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
2866 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2867 (deriv_data(i, j, k)*dra1dra(i, j, k) - &
2868 deriv_data2(i, j, k)*drb1drb(i, j, k))/s
2877 ELSE IF (nspins /= 1)
THEN
2881 IF (
ASSOCIATED(deriv_att))
THEN
2885 DO k = bo(1, 3), bo(2, 3)
2886 DO j = bo(1, 2), bo(2, 2)
2887 DO i = bo(1, 1), bo(2, 1)
2888 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2889 deriv_data(i, j, k)*rho1a(i, j, k)
2895 IF (
ASSOCIATED(deriv_att))
THEN
2899 DO k = bo(1, 3), bo(2, 3)
2900 DO j = bo(1, 2), bo(2, 2)
2901 DO i = bo(1, 1), bo(2, 1)
2902 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2903 deriv_data(i, j, k)*rho1b(i, j, k)
2909 IF (
ASSOCIATED(deriv_att))
THEN
2913 DO k = bo(1, 3), bo(2, 3)
2914 DO j = bo(1, 2), bo(2, 2)
2915 DO i = bo(1, 1), bo(2, 1)
2916 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2917 deriv_data(i, j, k)*dr1dr(i, j, k)
2923 IF (
ASSOCIATED(deriv_att))
THEN
2927 DO k = bo(1, 3), bo(2, 3)
2928 DO j = bo(1, 2), bo(2, 2)
2929 DO i = bo(1, 1), bo(2, 1)
2930 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2931 deriv_data(i, j, k)*dra1dra(i, j, k)
2937 IF (
ASSOCIATED(deriv_att))
THEN
2941 DO k = bo(1, 3), bo(2, 3)
2942 DO j = bo(1, 2), bo(2, 2)
2943 DO i = bo(1, 1), bo(2, 1)
2944 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2945 deriv_data(i, j, k)*drb1drb(i, j, k)
2951 IF (
ASSOCIATED(deriv_att))
THEN
2955 DO k = bo(1, 3), bo(2, 3)
2956 DO j = bo(1, 2), bo(2, 2)
2957 DO i = bo(1, 1), bo(2, 1)
2958 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2959 deriv_data(i, j, k)*tau1a(i, j, k)
2965 IF (
ASSOCIATED(deriv_att))
THEN
2969 DO k = bo(1, 3), bo(2, 3)
2970 DO j = bo(1, 2), bo(2, 2)
2971 DO i = bo(1, 1), bo(2, 1)
2972 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2973 deriv_data(i, j, k)*tau1b(i, j, k)
2979 IF (
ASSOCIATED(deriv_att))
THEN
2983 DO k = bo(1, 3), bo(2, 3)
2984 DO j = bo(1, 2), bo(2, 2)
2985 DO i = bo(1, 1), bo(2, 1)
2986 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
2987 deriv_data(i, j, k)*laplace1a(i, j, k)
2993 IF (
ASSOCIATED(deriv_att))
THEN
2997 DO k = bo(1, 3), bo(2, 3)
2998 DO j = bo(1, 2), bo(2, 2)
2999 DO i = bo(1, 1), bo(2, 1)
3000 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
3001 deriv_data(i, j, k)*laplace1b(i, j, k)
3009 IF (
ASSOCIATED(deriv_att))
THEN
3013 DO k = bo(1, 3), bo(2, 3)
3014 DO j = bo(1, 2), bo(2, 2)
3015 DO i = bo(1, 1), bo(2, 1)
3016 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
3017 deriv_data(i, j, k)*rho1a(i, j, k)
3023 IF (
ASSOCIATED(deriv_att))
THEN
3027 DO k = bo(1, 3), bo(2, 3)
3028 DO j = bo(1, 2), bo(2, 2)
3029 DO i = bo(1, 1), bo(2, 1)
3030 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
3031 deriv_data(i, j, k)*rho1b(i, j, k)
3037 IF (
ASSOCIATED(deriv_att))
THEN
3041 DO k = bo(1, 3), bo(2, 3)
3042 DO j = bo(1, 2), bo(2, 2)
3043 DO i = bo(1, 1), bo(2, 1)
3044 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
3045 deriv_data(i, j, k)*dr1dr(i, j, k)
3051 IF (
ASSOCIATED(deriv_att))
THEN
3055 DO k = bo(1, 3), bo(2, 3)
3056 DO j = bo(1, 2), bo(2, 2)
3057 DO i = bo(1, 1), bo(2, 1)
3058 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
3059 deriv_data(i, j, k)*dra1dra(i, j, k)
3065 IF (
ASSOCIATED(deriv_att))
THEN
3069 DO k = bo(1, 3), bo(2, 3)
3070 DO j = bo(1, 2), bo(2, 2)
3071 DO i = bo(1, 1), bo(2, 1)
3072 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
3073 deriv_data(i, j, k)*drb1drb(i, j, k)
3079 IF (
ASSOCIATED(deriv_att))
THEN
3083 DO k = bo(1, 3), bo(2, 3)
3084 DO j = bo(1, 2), bo(2, 2)
3085 DO i = bo(1, 1), bo(2, 1)
3086 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
3087 deriv_data(i, j, k)*tau1a(i, j, k)
3093 IF (
ASSOCIATED(deriv_att))
THEN
3097 DO k = bo(1, 3), bo(2, 3)
3098 DO j = bo(1, 2), bo(2, 2)
3099 DO i = bo(1, 1), bo(2, 1)
3100 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
3101 deriv_data(i, j, k)*tau1b(i, j, k)
3107 IF (
ASSOCIATED(deriv_att))
THEN
3111 DO k = bo(1, 3), bo(2, 3)
3112 DO j = bo(1, 2), bo(2, 2)
3113 DO i = bo(1, 1), bo(2, 1)
3114 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
3115 deriv_data(i, j, k)*laplace1a(i, j, k)
3121 IF (
ASSOCIATED(deriv_att))
THEN
3125 DO k = bo(1, 3), bo(2, 3)
3126 DO j = bo(1, 2), bo(2, 2)
3127 DO i = bo(1, 1), bo(2, 1)
3128 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
3129 deriv_data(i, j, k)*laplace1b(i, j, k)
3137 IF (
ASSOCIATED(deriv_att))
THEN
3141 DO k = bo(1, 3), bo(2, 3)
3142 DO j = bo(1, 2), bo(2, 2)
3143 DO i = bo(1, 1), bo(2, 1)
3144 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3145 deriv_data(i, j, k)*rho1a(i, j, k)
3146 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
3147 deriv_data(i, j, k)*rho1a(i, j, k)
3153 IF (
ASSOCIATED(deriv_att))
THEN
3157 DO k = bo(1, 3), bo(2, 3)
3158 DO j = bo(1, 2), bo(2, 2)
3159 DO i = bo(1, 1), bo(2, 1)
3160 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3161 deriv_data(i, j, k)*rho1b(i, j, k)
3162 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
3163 deriv_data(i, j, k)*rho1b(i, j, k)
3169 IF (
ASSOCIATED(deriv_att))
THEN
3173 DO k = bo(1, 3), bo(2, 3)
3174 DO j = bo(1, 2), bo(2, 2)
3175 DO i = bo(1, 1), bo(2, 1)
3176 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3177 deriv_data(i, j, k)*dr1dr(i, j, k)
3178 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
3179 deriv_data(i, j, k)*dr1dr(i, j, k)
3185 IF (
ASSOCIATED(deriv_att))
THEN
3189 DO k = bo(1, 3), bo(2, 3)
3190 DO j = bo(1, 2), bo(2, 2)
3191 DO i = bo(1, 1), bo(2, 1)
3192 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3193 deriv_data(i, j, k)*dra1dra(i, j, k)
3194 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
3195 deriv_data(i, j, k)*dra1dra(i, j, k)
3201 IF (
ASSOCIATED(deriv_att))
THEN
3205 DO k = bo(1, 3), bo(2, 3)
3206 DO j = bo(1, 2), bo(2, 2)
3207 DO i = bo(1, 1), bo(2, 1)
3208 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3209 deriv_data(i, j, k)*drb1drb(i, j, k)
3210 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
3211 deriv_data(i, j, k)*drb1drb(i, j, k)
3217 IF (
ASSOCIATED(deriv_att))
THEN
3221 DO k = bo(1, 3), bo(2, 3)
3222 DO j = bo(1, 2), bo(2, 2)
3223 DO i = bo(1, 1), bo(2, 1)
3224 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3225 deriv_data(i, j, k)*tau1a(i, j, k)
3226 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
3227 deriv_data(i, j, k)*tau1a(i, j, k)
3233 IF (
ASSOCIATED(deriv_att))
THEN
3237 DO k = bo(1, 3), bo(2, 3)
3238 DO j = bo(1, 2), bo(2, 2)
3239 DO i = bo(1, 1), bo(2, 1)
3240 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3241 deriv_data(i, j, k)*tau1b(i, j, k)
3242 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
3243 deriv_data(i, j, k)*tau1b(i, j, k)
3249 IF (
ASSOCIATED(deriv_att))
THEN
3253 DO k = bo(1, 3), bo(2, 3)
3254 DO j = bo(1, 2), bo(2, 2)
3255 DO i = bo(1, 1), bo(2, 1)
3256 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3257 deriv_data(i, j, k)*laplace1a(i, j, k)
3258 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
3259 deriv_data(i, j, k)*laplace1a(i, j, k)
3265 IF (
ASSOCIATED(deriv_att))
THEN
3269 DO k = bo(1, 3), bo(2, 3)
3270 DO j = bo(1, 2), bo(2, 2)
3271 DO i = bo(1, 1), bo(2, 1)
3272 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
3273 deriv_data(i, j, k)*laplace1b(i, j, k)
3274 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
3275 deriv_data(i, j, k)*laplace1b(i, j, k)
3282 IF (
ASSOCIATED(deriv_att))
THEN
3286 IF (my_compute_virial)
THEN
3287 CALL virial_drho_drho1(virial_pw, drho, drho1, deriv_data, virial_xc)
3291 v_drho(1)%array(:, :, :) = v_drho(1)%array(:, :, :) + &
3292 deriv_data(:, :, :)*dr1dr(:, :, :)/max(gradient_cut, norm_drho(:, :, :))**2
3293 v_drho(2)%array(:, :, :) = v_drho(2)%array(:, :, :) + &
3294 deriv_data(:, :, :)*dr1dr(:, :, :)/max(gradient_cut, norm_drho(:, :, :))**2
3299 IF (
ASSOCIATED(deriv_att))
THEN
3303 DO k = bo(1, 3), bo(2, 3)
3304 DO j = bo(1, 2), bo(2, 2)
3305 DO i = bo(1, 1), bo(2, 1)
3306 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3307 deriv_data(i, j, k)*rho1a(i, j, k)
3313 IF (
ASSOCIATED(deriv_att))
THEN
3317 DO k = bo(1, 3), bo(2, 3)
3318 DO j = bo(1, 2), bo(2, 2)
3319 DO i = bo(1, 1), bo(2, 1)
3320 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3321 deriv_data(i, j, k)*rho1b(i, j, k)
3327 IF (
ASSOCIATED(deriv_att))
THEN
3331 DO k = bo(1, 3), bo(2, 3)
3332 DO j = bo(1, 2), bo(2, 2)
3333 DO i = bo(1, 1), bo(2, 1)
3334 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3335 deriv_data(i, j, k)*dr1dr(i, j, k)
3341 IF (
ASSOCIATED(deriv_att))
THEN
3345 DO k = bo(1, 3), bo(2, 3)
3346 DO j = bo(1, 2), bo(2, 2)
3347 DO i = bo(1, 1), bo(2, 1)
3348 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3349 deriv_data(i, j, k)*dra1dra(i, j, k)
3355 IF (
ASSOCIATED(deriv_att))
THEN
3359 DO k = bo(1, 3), bo(2, 3)
3360 DO j = bo(1, 2), bo(2, 2)
3361 DO i = bo(1, 1), bo(2, 1)
3362 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3363 deriv_data(i, j, k)*drb1drb(i, j, k)
3369 IF (
ASSOCIATED(deriv_att))
THEN
3373 DO k = bo(1, 3), bo(2, 3)
3374 DO j = bo(1, 2), bo(2, 2)
3375 DO i = bo(1, 1), bo(2, 1)
3376 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3377 deriv_data(i, j, k)*tau1a(i, j, k)
3383 IF (
ASSOCIATED(deriv_att))
THEN
3387 DO k = bo(1, 3), bo(2, 3)
3388 DO j = bo(1, 2), bo(2, 2)
3389 DO i = bo(1, 1), bo(2, 1)
3390 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3391 deriv_data(i, j, k)*tau1b(i, j, k)
3397 IF (
ASSOCIATED(deriv_att))
THEN
3401 DO k = bo(1, 3), bo(2, 3)
3402 DO j = bo(1, 2), bo(2, 2)
3403 DO i = bo(1, 1), bo(2, 1)
3404 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3405 deriv_data(i, j, k)*laplace1a(i, j, k)
3411 IF (
ASSOCIATED(deriv_att))
THEN
3415 DO k = bo(1, 3), bo(2, 3)
3416 DO j = bo(1, 2), bo(2, 2)
3417 DO i = bo(1, 1), bo(2, 1)
3418 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
3419 deriv_data(i, j, k)*laplace1b(i, j, k)
3426 IF (
ASSOCIATED(deriv_att))
THEN
3430 IF (my_compute_virial)
THEN
3431 CALL virial_drho_drho1(virial_pw, drhoa, drho1a, deriv_data, virial_xc)
3435 v_drhoa(1)%array(:, :, :) = v_drhoa(1)%array(:, :, :) + &
3436 deriv_data(:, :, :)*dra1dra(:, :, :)/max(gradient_cut, norm_drhoa(:, :, :))**2
3441 IF (
ASSOCIATED(deriv_att))
THEN
3445 DO k = bo(1, 3), bo(2, 3)
3446 DO j = bo(1, 2), bo(2, 2)
3447 DO i = bo(1, 1), bo(2, 1)
3448 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
3449 deriv_data(i, j, k)*rho1a(i, j, k)
3455 IF (
ASSOCIATED(deriv_att))
THEN
3459 DO k = bo(1, 3), bo(2, 3)
3460 DO j = bo(1, 2), bo(2, 2)
3461 DO i = bo(1, 1), bo(2, 1)
3462 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
3463 deriv_data(i, j, k)*rho1b(i, j, k)
3469 IF (
ASSOCIATED(deriv_att))
THEN
3473 DO k = bo(1, 3), bo(2, 3)
3474 DO j = bo(1, 2), bo(2, 2)
3475 DO i = bo(1, 1), bo(2, 1)
3476 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
3477 deriv_data(i, j, k)*dr1dr(i, j, k)
3483 IF (
ASSOCIATED(deriv_att))
THEN
3487 DO k = bo(1, 3), bo(2, 3)
3488 DO j = bo(1, 2), bo(2, 2)
3489 DO i = bo(1, 1), bo(2, 1)
3490 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
3491 deriv_data(i, j, k)*dra1dra(i, j, k)
3497 IF (
ASSOCIATED(deriv_att))
THEN
3501 DO k = bo(1, 3), bo(2, 3)
3502 DO j = bo(1, 2), bo(2, 2)
3503 DO i = bo(1, 1), bo(2, 1)
3504 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
3505 deriv_data(i, j, k)*drb1drb(i, j, k)
3511 IF (
ASSOCIATED(deriv_att))
THEN
3515 DO k = bo(1, 3), bo(2, 3)
3516 DO j = bo(1, 2), bo(2, 2)
3517 DO i = bo(1, 1), bo(2, 1)
3518 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
3519 deriv_data(i, j, k)*tau1a(i, j, k)
3525 IF (
ASSOCIATED(deriv_att))
THEN
3529 DO k = bo(1, 3), bo(2, 3)
3530 DO j = bo(1, 2), bo(2, 2)
3531 DO i = bo(1, 1), bo(2, 1)
3532 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
3533 deriv_data(i, j, k)*tau1b(i, j, k)
3539 IF (
ASSOCIATED(deriv_att))
THEN
3543 DO k = bo(1, 3), bo(2, 3)
3544 DO j = bo(1, 2), bo(2, 2)
3545 DO i = bo(1, 1), bo(2, 1)
3546 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
3547 deriv_data(i, j, k)*laplace1a(i, j, k)
3553 IF (
ASSOCIATED(deriv_att))
THEN
3557 DO k = bo(1, 3), bo(2, 3)
3558 DO j = bo(1, 2), bo(2, 2)
3559 DO i = bo(1, 1), bo(2, 1)
3560 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
3561 deriv_data(i, j, k)*laplace1b(i, j, k)
3568 IF (
ASSOCIATED(deriv_att))
THEN
3572 IF (my_compute_virial)
THEN
3573 CALL virial_drho_drho1(virial_pw, drhob, drho1b, deriv_data, virial_xc)
3577 v_drhob(2)%array(:, :, :) = v_drhob(2)%array(:, :, :) + &
3578 deriv_data(:, :, :)*drb1drb(:, :, :)/max(gradient_cut, norm_drhob(:, :, :))**2
3583 IF (
ASSOCIATED(deriv_att))
THEN
3587 DO k = bo(1, 3), bo(2, 3)
3588 DO j = bo(1, 2), bo(2, 2)
3589 DO i = bo(1, 1), bo(2, 1)
3590 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3591 deriv_data(i, j, k)*rho1a(i, j, k)
3597 IF (
ASSOCIATED(deriv_att))
THEN
3601 DO k = bo(1, 3), bo(2, 3)
3602 DO j = bo(1, 2), bo(2, 2)
3603 DO i = bo(1, 1), bo(2, 1)
3604 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3605 deriv_data(i, j, k)*rho1b(i, j, k)
3611 IF (
ASSOCIATED(deriv_att))
THEN
3615 DO k = bo(1, 3), bo(2, 3)
3616 DO j = bo(1, 2), bo(2, 2)
3617 DO i = bo(1, 1), bo(2, 1)
3618 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3619 deriv_data(i, j, k)*dr1dr(i, j, k)
3625 IF (
ASSOCIATED(deriv_att))
THEN
3629 DO k = bo(1, 3), bo(2, 3)
3630 DO j = bo(1, 2), bo(2, 2)
3631 DO i = bo(1, 1), bo(2, 1)
3632 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3633 deriv_data(i, j, k)*dra1dra(i, j, k)
3639 IF (
ASSOCIATED(deriv_att))
THEN
3643 DO k = bo(1, 3), bo(2, 3)
3644 DO j = bo(1, 2), bo(2, 2)
3645 DO i = bo(1, 1), bo(2, 1)
3646 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3647 deriv_data(i, j, k)*drb1drb(i, j, k)
3653 IF (
ASSOCIATED(deriv_att))
THEN
3657 DO k = bo(1, 3), bo(2, 3)
3658 DO j = bo(1, 2), bo(2, 2)
3659 DO i = bo(1, 1), bo(2, 1)
3660 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3661 deriv_data(i, j, k)*tau1a(i, j, k)
3667 IF (
ASSOCIATED(deriv_att))
THEN
3671 DO k = bo(1, 3), bo(2, 3)
3672 DO j = bo(1, 2), bo(2, 2)
3673 DO i = bo(1, 1), bo(2, 1)
3674 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3675 deriv_data(i, j, k)*tau1b(i, j, k)
3681 IF (
ASSOCIATED(deriv_att))
THEN
3685 DO k = bo(1, 3), bo(2, 3)
3686 DO j = bo(1, 2), bo(2, 2)
3687 DO i = bo(1, 1), bo(2, 1)
3688 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3689 deriv_data(i, j, k)*laplace1a(i, j, k)
3695 IF (
ASSOCIATED(deriv_att))
THEN
3699 DO k = bo(1, 3), bo(2, 3)
3700 DO j = bo(1, 2), bo(2, 2)
3701 DO i = bo(1, 1), bo(2, 1)
3702 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
3703 deriv_data(i, j, k)*laplace1b(i, j, k)
3711 IF (
ASSOCIATED(deriv_att))
THEN
3715 DO k = bo(1, 3), bo(2, 3)
3716 DO j = bo(1, 2), bo(2, 2)
3717 DO i = bo(1, 1), bo(2, 1)
3718 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3719 deriv_data(i, j, k)*rho1a(i, j, k)
3725 IF (
ASSOCIATED(deriv_att))
THEN
3729 DO k = bo(1, 3), bo(2, 3)
3730 DO j = bo(1, 2), bo(2, 2)
3731 DO i = bo(1, 1), bo(2, 1)
3732 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3733 deriv_data(i, j, k)*rho1b(i, j, k)
3739 IF (
ASSOCIATED(deriv_att))
THEN
3743 DO k = bo(1, 3), bo(2, 3)
3744 DO j = bo(1, 2), bo(2, 2)
3745 DO i = bo(1, 1), bo(2, 1)
3746 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3747 deriv_data(i, j, k)*dr1dr(i, j, k)
3753 IF (
ASSOCIATED(deriv_att))
THEN
3757 DO k = bo(1, 3), bo(2, 3)
3758 DO j = bo(1, 2), bo(2, 2)
3759 DO i = bo(1, 1), bo(2, 1)
3760 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3761 deriv_data(i, j, k)*dra1dra(i, j, k)
3767 IF (
ASSOCIATED(deriv_att))
THEN
3771 DO k = bo(1, 3), bo(2, 3)
3772 DO j = bo(1, 2), bo(2, 2)
3773 DO i = bo(1, 1), bo(2, 1)
3774 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3775 deriv_data(i, j, k)*drb1drb(i, j, k)
3781 IF (
ASSOCIATED(deriv_att))
THEN
3785 DO k = bo(1, 3), bo(2, 3)
3786 DO j = bo(1, 2), bo(2, 2)
3787 DO i = bo(1, 1), bo(2, 1)
3788 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3789 deriv_data(i, j, k)*tau1a(i, j, k)
3795 IF (
ASSOCIATED(deriv_att))
THEN
3799 DO k = bo(1, 3), bo(2, 3)
3800 DO j = bo(1, 2), bo(2, 2)
3801 DO i = bo(1, 1), bo(2, 1)
3802 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3803 deriv_data(i, j, k)*tau1b(i, j, k)
3809 IF (
ASSOCIATED(deriv_att))
THEN
3813 DO k = bo(1, 3), bo(2, 3)
3814 DO j = bo(1, 2), bo(2, 2)
3815 DO i = bo(1, 1), bo(2, 1)
3816 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3817 deriv_data(i, j, k)*laplace1a(i, j, k)
3823 IF (
ASSOCIATED(deriv_att))
THEN
3827 DO k = bo(1, 3), bo(2, 3)
3828 DO j = bo(1, 2), bo(2, 2)
3829 DO i = bo(1, 1), bo(2, 1)
3830 v_xc_tau(2)%array(i, j, k) = v_xc_tau(2)%array(i, j, k) + &
3831 deriv_data(i, j, k)*laplace1b(i, j, k)
3839 IF (
ASSOCIATED(deriv_att))
THEN
3843 DO k = bo(1, 3), bo(2, 3)
3844 DO j = bo(1, 2), bo(2, 2)
3845 DO i = bo(1, 1), bo(2, 1)
3846 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3847 deriv_data(i, j, k)*rho1a(i, j, k)
3853 IF (
ASSOCIATED(deriv_att))
THEN
3857 DO k = bo(1, 3), bo(2, 3)
3858 DO j = bo(1, 2), bo(2, 2)
3859 DO i = bo(1, 1), bo(2, 1)
3860 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3861 deriv_data(i, j, k)*rho1b(i, j, k)
3867 IF (
ASSOCIATED(deriv_att))
THEN
3871 DO k = bo(1, 3), bo(2, 3)
3872 DO j = bo(1, 2), bo(2, 2)
3873 DO i = bo(1, 1), bo(2, 1)
3874 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3875 deriv_data(i, j, k)*dr1dr(i, j, k)
3881 IF (
ASSOCIATED(deriv_att))
THEN
3885 DO k = bo(1, 3), bo(2, 3)
3886 DO j = bo(1, 2), bo(2, 2)
3887 DO i = bo(1, 1), bo(2, 1)
3888 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3889 deriv_data(i, j, k)*dra1dra(i, j, k)
3895 IF (
ASSOCIATED(deriv_att))
THEN
3899 DO k = bo(1, 3), bo(2, 3)
3900 DO j = bo(1, 2), bo(2, 2)
3901 DO i = bo(1, 1), bo(2, 1)
3902 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3903 deriv_data(i, j, k)*drb1drb(i, j, k)
3909 IF (
ASSOCIATED(deriv_att))
THEN
3913 DO k = bo(1, 3), bo(2, 3)
3914 DO j = bo(1, 2), bo(2, 2)
3915 DO i = bo(1, 1), bo(2, 1)
3916 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3917 deriv_data(i, j, k)*tau1a(i, j, k)
3923 IF (
ASSOCIATED(deriv_att))
THEN
3927 DO k = bo(1, 3), bo(2, 3)
3928 DO j = bo(1, 2), bo(2, 2)
3929 DO i = bo(1, 1), bo(2, 1)
3930 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3931 deriv_data(i, j, k)*tau1b(i, j, k)
3937 IF (
ASSOCIATED(deriv_att))
THEN
3941 DO k = bo(1, 3), bo(2, 3)
3942 DO j = bo(1, 2), bo(2, 2)
3943 DO i = bo(1, 1), bo(2, 1)
3944 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3945 deriv_data(i, j, k)*laplace1a(i, j, k)
3951 IF (
ASSOCIATED(deriv_att))
THEN
3955 DO k = bo(1, 3), bo(2, 3)
3956 DO j = bo(1, 2), bo(2, 2)
3957 DO i = bo(1, 1), bo(2, 1)
3958 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
3959 deriv_data(i, j, k)*laplace1b(i, j, k)
3966 IF (my_compute_virial)
THEN
3968 IF (
ASSOCIATED(deriv_att))
THEN
3971 virial_pw%array(:, :, :) = -rho1a(:, :, :)
3972 CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
3976 IF (
ASSOCIATED(deriv_att))
THEN
3980 DO k = bo(1, 3), bo(2, 3)
3981 DO j = bo(1, 2), bo(2, 2)
3982 DO i = bo(1, 1), bo(2, 1)
3983 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3984 deriv_data(i, j, k)*rho1a(i, j, k)
3990 IF (
ASSOCIATED(deriv_att))
THEN
3994 DO k = bo(1, 3), bo(2, 3)
3995 DO j = bo(1, 2), bo(2, 2)
3996 DO i = bo(1, 1), bo(2, 1)
3997 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
3998 deriv_data(i, j, k)*rho1b(i, j, k)
4004 IF (
ASSOCIATED(deriv_att))
THEN
4008 DO k = bo(1, 3), bo(2, 3)
4009 DO j = bo(1, 2), bo(2, 2)
4010 DO i = bo(1, 1), bo(2, 1)
4011 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4012 deriv_data(i, j, k)*dr1dr(i, j, k)
4018 IF (
ASSOCIATED(deriv_att))
THEN
4022 DO k = bo(1, 3), bo(2, 3)
4023 DO j = bo(1, 2), bo(2, 2)
4024 DO i = bo(1, 1), bo(2, 1)
4025 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4026 deriv_data(i, j, k)*dra1dra(i, j, k)
4032 IF (
ASSOCIATED(deriv_att))
THEN
4036 DO k = bo(1, 3), bo(2, 3)
4037 DO j = bo(1, 2), bo(2, 2)
4038 DO i = bo(1, 1), bo(2, 1)
4039 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4040 deriv_data(i, j, k)*drb1drb(i, j, k)
4046 IF (
ASSOCIATED(deriv_att))
THEN
4050 DO k = bo(1, 3), bo(2, 3)
4051 DO j = bo(1, 2), bo(2, 2)
4052 DO i = bo(1, 1), bo(2, 1)
4053 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4054 deriv_data(i, j, k)*tau1a(i, j, k)
4060 IF (
ASSOCIATED(deriv_att))
THEN
4064 DO k = bo(1, 3), bo(2, 3)
4065 DO j = bo(1, 2), bo(2, 2)
4066 DO i = bo(1, 1), bo(2, 1)
4067 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4068 deriv_data(i, j, k)*tau1b(i, j, k)
4074 IF (
ASSOCIATED(deriv_att))
THEN
4078 DO k = bo(1, 3), bo(2, 3)
4079 DO j = bo(1, 2), bo(2, 2)
4080 DO i = bo(1, 1), bo(2, 1)
4081 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4082 deriv_data(i, j, k)*laplace1a(i, j, k)
4088 IF (
ASSOCIATED(deriv_att))
THEN
4092 DO k = bo(1, 3), bo(2, 3)
4093 DO j = bo(1, 2), bo(2, 2)
4094 DO i = bo(1, 1), bo(2, 1)
4095 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4096 deriv_data(i, j, k)*laplace1b(i, j, k)
4103 IF (my_compute_virial)
THEN
4105 IF (
ASSOCIATED(deriv_att))
THEN
4108 virial_pw%array(:, :, :) = -rho1b(:, :, :)
4109 CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
4118 IF (
ASSOCIATED(deriv_att))
THEN
4122 DO k = bo(1, 3), bo(2, 3)
4123 DO j = bo(1, 2), bo(2, 2)
4124 DO i = bo(1, 1), bo(2, 1)
4125 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4126 deriv_data(i, j, k)*rho1a(i, j, k)
4132 IF (
ASSOCIATED(deriv_att))
THEN
4136 DO k = bo(1, 3), bo(2, 3)
4137 DO j = bo(1, 2), bo(2, 2)
4138 DO i = bo(1, 1), bo(2, 1)
4139 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4140 deriv_data(i, j, k)*dr1dr(i, j, k)
4146 IF (
ASSOCIATED(deriv_att))
THEN
4150 DO k = bo(1, 3), bo(2, 3)
4151 DO j = bo(1, 2), bo(2, 2)
4152 DO i = bo(1, 1), bo(2, 1)
4153 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4154 deriv_data(i, j, k)*dra1dra(i, j, k)
4160 IF (
ASSOCIATED(deriv_att))
THEN
4164 DO k = bo(1, 3), bo(2, 3)
4165 DO j = bo(1, 2), bo(2, 2)
4166 DO i = bo(1, 1), bo(2, 1)
4167 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4168 deriv_data(i, j, k)*tau1a(i, j, k)
4174 IF (
ASSOCIATED(deriv_att))
THEN
4178 DO k = bo(1, 3), bo(2, 3)
4179 DO j = bo(1, 2), bo(2, 2)
4180 DO i = bo(1, 1), bo(2, 1)
4181 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4182 deriv_data(i, j, k)*laplace1a(i, j, k)
4188 IF (
ASSOCIATED(deriv_att))
THEN
4192 DO k = bo(1, 3), bo(2, 3)
4193 DO j = bo(1, 2), bo(2, 2)
4194 DO i = bo(1, 1), bo(2, 1)
4195 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4196 fac*deriv_data(i, j, k)*rho1b(i, j, k)
4202 IF (
ASSOCIATED(deriv_att))
THEN
4206 DO k = bo(1, 3), bo(2, 3)
4207 DO j = bo(1, 2), bo(2, 2)
4208 DO i = bo(1, 1), bo(2, 1)
4209 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4210 fac*deriv_data(i, j, k)*drb1drb(i, j, k)
4216 IF (
ASSOCIATED(deriv_att))
THEN
4220 DO k = bo(1, 3), bo(2, 3)
4221 DO j = bo(1, 2), bo(2, 2)
4222 DO i = bo(1, 1), bo(2, 1)
4223 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4224 fac*deriv_data(i, j, k)*tau1b(i, j, k)
4230 IF (
ASSOCIATED(deriv_att))
THEN
4234 DO k = bo(1, 3), bo(2, 3)
4235 DO j = bo(1, 2), bo(2, 2)
4236 DO i = bo(1, 1), bo(2, 1)
4237 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4238 fac*deriv_data(i, j, k)*laplace1b(i, j, k)
4246 IF (
ASSOCIATED(deriv_att))
THEN
4250 DO k = bo(1, 3), bo(2, 3)
4251 DO j = bo(1, 2), bo(2, 2)
4252 DO i = bo(1, 1), bo(2, 1)
4253 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
4254 deriv_data(i, j, k)*rho1a(i, j, k)
4260 IF (
ASSOCIATED(deriv_att))
THEN
4264 DO k = bo(1, 3), bo(2, 3)
4265 DO j = bo(1, 2), bo(2, 2)
4266 DO i = bo(1, 1), bo(2, 1)
4267 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
4268 deriv_data(i, j, k)*dr1dr(i, j, k)
4274 IF (
ASSOCIATED(deriv_att))
THEN
4278 DO k = bo(1, 3), bo(2, 3)
4279 DO j = bo(1, 2), bo(2, 2)
4280 DO i = bo(1, 1), bo(2, 1)
4281 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
4282 deriv_data(i, j, k)*dra1dra(i, j, k)
4288 IF (
ASSOCIATED(deriv_att))
THEN
4292 DO k = bo(1, 3), bo(2, 3)
4293 DO j = bo(1, 2), bo(2, 2)
4294 DO i = bo(1, 1), bo(2, 1)
4295 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
4296 deriv_data(i, j, k)*tau1a(i, j, k)
4302 IF (
ASSOCIATED(deriv_att))
THEN
4306 DO k = bo(1, 3), bo(2, 3)
4307 DO j = bo(1, 2), bo(2, 2)
4308 DO i = bo(1, 1), bo(2, 1)
4309 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
4310 deriv_data(i, j, k)*laplace1a(i, j, k)
4316 IF (
ASSOCIATED(deriv_att))
THEN
4320 DO k = bo(1, 3), bo(2, 3)
4321 DO j = bo(1, 2), bo(2, 2)
4322 DO i = bo(1, 1), bo(2, 1)
4323 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
4324 fac*deriv_data(i, j, k)*rho1b(i, j, k)
4330 IF (
ASSOCIATED(deriv_att))
THEN
4334 DO k = bo(1, 3), bo(2, 3)
4335 DO j = bo(1, 2), bo(2, 2)
4336 DO i = bo(1, 1), bo(2, 1)
4337 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
4338 fac*deriv_data(i, j, k)*drb1drb(i, j, k)
4344 IF (
ASSOCIATED(deriv_att))
THEN
4348 DO k = bo(1, 3), bo(2, 3)
4349 DO j = bo(1, 2), bo(2, 2)
4350 DO i = bo(1, 1), bo(2, 1)
4351 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
4352 fac*deriv_data(i, j, k)*tau1b(i, j, k)
4358 IF (
ASSOCIATED(deriv_att))
THEN
4362 DO k = bo(1, 3), bo(2, 3)
4363 DO j = bo(1, 2), bo(2, 2)
4364 DO i = bo(1, 1), bo(2, 1)
4365 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
4366 fac*deriv_data(i, j, k)*laplace1b(i, j, k)
4373 IF (
ASSOCIATED(deriv_att))
THEN
4379 v_drho(1)%array(:, :, :) = v_drho(1)%array(:, :, :) + &
4380 deriv_data(:, :, :)*dr1dr(:, :, :)/max(gradient_cut, norm_drho(:, :, :))**2
4385 IF (
ASSOCIATED(deriv_att))
THEN
4389 DO k = bo(1, 3), bo(2, 3)
4390 DO j = bo(1, 2), bo(2, 2)
4391 DO i = bo(1, 1), bo(2, 1)
4392 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
4393 deriv_data(i, j, k)*rho1a(i, j, k)
4399 IF (
ASSOCIATED(deriv_att))
THEN
4403 DO k = bo(1, 3), bo(2, 3)
4404 DO j = bo(1, 2), bo(2, 2)
4405 DO i = bo(1, 1), bo(2, 1)
4406 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
4407 deriv_data(i, j, k)*dr1dr(i, j, k)
4413 IF (
ASSOCIATED(deriv_att))
THEN
4417 DO k = bo(1, 3), bo(2, 3)
4418 DO j = bo(1, 2), bo(2, 2)
4419 DO i = bo(1, 1), bo(2, 1)
4420 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
4421 deriv_data(i, j, k)*dra1dra(i, j, k)
4427 IF (
ASSOCIATED(deriv_att))
THEN
4431 DO k = bo(1, 3), bo(2, 3)
4432 DO j = bo(1, 2), bo(2, 2)
4433 DO i = bo(1, 1), bo(2, 1)
4434 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
4435 deriv_data(i, j, k)*tau1a(i, j, k)
4441 IF (
ASSOCIATED(deriv_att))
THEN
4445 DO k = bo(1, 3), bo(2, 3)
4446 DO j = bo(1, 2), bo(2, 2)
4447 DO i = bo(1, 1), bo(2, 1)
4448 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
4449 deriv_data(i, j, k)*laplace1a(i, j, k)
4455 IF (
ASSOCIATED(deriv_att))
THEN
4459 DO k = bo(1, 3), bo(2, 3)
4460 DO j = bo(1, 2), bo(2, 2)
4461 DO i = bo(1, 1), bo(2, 1)
4462 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
4463 fac*deriv_data(i, j, k)*rho1b(i, j, k)
4469 IF (
ASSOCIATED(deriv_att))
THEN
4473 DO k = bo(1, 3), bo(2, 3)
4474 DO j = bo(1, 2), bo(2, 2)
4475 DO i = bo(1, 1), bo(2, 1)
4476 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
4477 fac*deriv_data(i, j, k)*drb1drb(i, j, k)
4483 IF (
ASSOCIATED(deriv_att))
THEN
4487 DO k = bo(1, 3), bo(2, 3)
4488 DO j = bo(1, 2), bo(2, 2)
4489 DO i = bo(1, 1), bo(2, 1)
4490 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
4491 fac*deriv_data(i, j, k)*tau1b(i, j, k)
4497 IF (
ASSOCIATED(deriv_att))
THEN
4501 DO k = bo(1, 3), bo(2, 3)
4502 DO j = bo(1, 2), bo(2, 2)
4503 DO i = bo(1, 1), bo(2, 1)
4504 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
4505 fac*deriv_data(i, j, k)*laplace1b(i, j, k)
4512 IF (
ASSOCIATED(deriv_att))
THEN
4518 v_drhoa(1)%array(:, :, :) = v_drhoa(1)%array(:, :, :) + &
4519 deriv_data(:, :, :)*dra1dra(:, :, :)/max(gradient_cut, norm_drhoa(:, :, :))**2
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)*rho1a(i, j, k)
4538 IF (
ASSOCIATED(deriv_att))
THEN
4542 DO k = bo(1, 3), bo(2, 3)
4543 DO j = bo(1, 2), bo(2, 2)
4544 DO i = bo(1, 1), bo(2, 1)
4545 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4546 deriv_data(i, j, k)*dr1dr(i, j, k)
4552 IF (
ASSOCIATED(deriv_att))
THEN
4556 DO k = bo(1, 3), bo(2, 3)
4557 DO j = bo(1, 2), bo(2, 2)
4558 DO i = bo(1, 1), bo(2, 1)
4559 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4560 deriv_data(i, j, k)*dra1dra(i, j, k)
4566 IF (
ASSOCIATED(deriv_att))
THEN
4570 DO k = bo(1, 3), bo(2, 3)
4571 DO j = bo(1, 2), bo(2, 2)
4572 DO i = bo(1, 1), bo(2, 1)
4573 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4574 deriv_data(i, j, k)*tau1a(i, j, k)
4580 IF (
ASSOCIATED(deriv_att))
THEN
4584 DO k = bo(1, 3), bo(2, 3)
4585 DO j = bo(1, 2), bo(2, 2)
4586 DO i = bo(1, 1), bo(2, 1)
4587 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4588 deriv_data(i, j, k)*laplace1a(i, j, k)
4594 IF (
ASSOCIATED(deriv_att))
THEN
4598 DO k = bo(1, 3), bo(2, 3)
4599 DO j = bo(1, 2), bo(2, 2)
4600 DO i = bo(1, 1), bo(2, 1)
4601 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4602 fac*deriv_data(i, j, k)*rho1b(i, j, k)
4608 IF (
ASSOCIATED(deriv_att))
THEN
4612 DO k = bo(1, 3), bo(2, 3)
4613 DO j = bo(1, 2), bo(2, 2)
4614 DO i = bo(1, 1), bo(2, 1)
4615 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4616 fac*deriv_data(i, j, k)*drb1drb(i, j, k)
4622 IF (
ASSOCIATED(deriv_att))
THEN
4626 DO k = bo(1, 3), bo(2, 3)
4627 DO j = bo(1, 2), bo(2, 2)
4628 DO i = bo(1, 1), bo(2, 1)
4629 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4630 fac*deriv_data(i, j, k)*tau1b(i, j, k)
4636 IF (
ASSOCIATED(deriv_att))
THEN
4640 DO k = bo(1, 3), bo(2, 3)
4641 DO j = bo(1, 2), bo(2, 2)
4642 DO i = bo(1, 1), bo(2, 1)
4643 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
4644 fac*deriv_data(i, j, k)*laplace1b(i, j, k)
4652 IF (
ASSOCIATED(deriv_att))
THEN
4656 DO k = bo(1, 3), bo(2, 3)
4657 DO j = bo(1, 2), bo(2, 2)
4658 DO i = bo(1, 1), bo(2, 1)
4659 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4660 deriv_data(i, j, k)*rho1a(i, j, k)
4666 IF (
ASSOCIATED(deriv_att))
THEN
4670 DO k = bo(1, 3), bo(2, 3)
4671 DO j = bo(1, 2), bo(2, 2)
4672 DO i = bo(1, 1), bo(2, 1)
4673 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4674 deriv_data(i, j, k)*dr1dr(i, j, k)
4680 IF (
ASSOCIATED(deriv_att))
THEN
4684 DO k = bo(1, 3), bo(2, 3)
4685 DO j = bo(1, 2), bo(2, 2)
4686 DO i = bo(1, 1), bo(2, 1)
4687 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4688 deriv_data(i, j, k)*dra1dra(i, j, k)
4694 IF (
ASSOCIATED(deriv_att))
THEN
4698 DO k = bo(1, 3), bo(2, 3)
4699 DO j = bo(1, 2), bo(2, 2)
4700 DO i = bo(1, 1), bo(2, 1)
4701 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4702 deriv_data(i, j, k)*tau1a(i, j, k)
4708 IF (
ASSOCIATED(deriv_att))
THEN
4712 DO k = bo(1, 3), bo(2, 3)
4713 DO j = bo(1, 2), bo(2, 2)
4714 DO i = bo(1, 1), bo(2, 1)
4715 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4716 deriv_data(i, j, k)*laplace1a(i, j, k)
4722 IF (
ASSOCIATED(deriv_att))
THEN
4726 DO k = bo(1, 3), bo(2, 3)
4727 DO j = bo(1, 2), bo(2, 2)
4728 DO i = bo(1, 1), bo(2, 1)
4729 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4730 fac*deriv_data(i, j, k)*rho1b(i, j, k)
4736 IF (
ASSOCIATED(deriv_att))
THEN
4740 DO k = bo(1, 3), bo(2, 3)
4741 DO j = bo(1, 2), bo(2, 2)
4742 DO i = bo(1, 1), bo(2, 1)
4743 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4744 fac*deriv_data(i, j, k)*drb1drb(i, j, k)
4750 IF (
ASSOCIATED(deriv_att))
THEN
4754 DO k = bo(1, 3), bo(2, 3)
4755 DO j = bo(1, 2), bo(2, 2)
4756 DO i = bo(1, 1), bo(2, 1)
4757 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4758 fac*deriv_data(i, j, k)*tau1b(i, j, k)
4764 IF (
ASSOCIATED(deriv_att))
THEN
4768 DO k = bo(1, 3), bo(2, 3)
4769 DO j = bo(1, 2), bo(2, 2)
4770 DO i = bo(1, 1), bo(2, 1)
4771 v_laplace(2)%array(i, j, k) = v_laplace(2)%array(i, j, k) + &
4772 fac*deriv_data(i, j, k)*laplace1b(i, j, k)
4783 IF (gradient_f)
THEN
4784 IF (.NOT. do_spinflip)
THEN
4786 IF (my_compute_virial)
THEN
4787 CALL virial_drho_drho(virial_pw, drhoa, v_drhoa(1), virial_xc)
4788 CALL virial_drho_drho(virial_pw, drhob, v_drhob(2), virial_xc)
4791 virial_pw%array(:, :, :) = drho(idir)%array(:, :, :)*(v_drho(1)%array(:, :, :) + v_drho(2)%array(:, :, :))
4795 drho(jdir)%array(:, :, :))
4796 virial_xc(jdir, idir) = virial_xc(jdir, idir) + tmp
4797 virial_xc(idir, jdir) = virial_xc(jdir, idir)
4807 DO ir = bo(1, 2), bo(2, 2)
4808 DO ia = bo(1, 1), bo(2, 1)
4810 DO ispin = 1, nspins
4811 vxg(idir, ia, ir, ispin) = &
4812 -(v_drhoa(ispin)%array(ia, ir, 1)*drhoa(idir)%array(ia, ir, 1) + &
4813 v_drhob(ispin)%array(ia, ir, 1)*drhob(idir)%array(ia, ir, 1) + &
4814 v_drho(ispin)%array(ia, ir, 1)*drho(idir)%array(ia, ir, 1))
4816 IF (
ASSOCIATED(e_drhoa))
THEN
4817 vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
4818 e_drhoa(ia, ir, 1)*drho1a(idir)%array(ia, ir, 1)
4820 IF (nspins /= 1 .AND.
ASSOCIATED(e_drhob))
THEN
4821 vxg(idir, ia, ir, 2) = vxg(idir, ia, ir, 2) + &
4822 e_drhob(ia, ir, 1)*drho1b(idir)%array(ia, ir, 1)
4824 IF (
ASSOCIATED(e_drho))
THEN
4825 IF (nspins /= 1)
THEN
4826 vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
4827 e_drho(ia, ir, 1)*drho1(idir)%array(ia, ir, 1)
4828 vxg(idir, ia, ir, 2) = vxg(idir, ia, ir, 2) + &
4829 e_drho(ia, ir, 1)*drho1(idir)%array(ia, ir, 1)
4831 vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
4832 e_drho(ia, ir, 1)*(drho1a(idir)%array(ia, ir, 1) + &
4833 fac*drho1b(idir)%array(ia, ir, 1))
4845 DO ispin = 1, nspins
4848 v_drho_r(idir, ispin)%array(:, :, :) = &
4849 v_drhoa(ispin)%array(:, :, :)*drhoa(idir)%array(:, :, :) + &
4850 v_drhob(ispin)%array(:, :, :)*drhob(idir)%array(:, :, :) + &
4851 v_drho(ispin)%array(:, :, :)*drho(idir)%array(:, :, :)
4854 IF (
ASSOCIATED(e_drhoa))
THEN
4857 v_drho_r(idir, 1)%array(:, :, :) = v_drho_r(idir, 1)%array(:, :, :) - &
4858 e_drhoa(:, :, :)*drho1a(idir)%array(:, :, :)
4861 IF (nspins /= 1 .AND.
ASSOCIATED(e_drhob))
THEN
4864 v_drho_r(idir, 2)%array(:, :, :) = v_drho_r(idir, 2)%array(:, :, :) - &
4865 e_drhob(:, :, :)*drho1b(idir)%array(:, :, :)
4868 IF (
ASSOCIATED(e_drho))
THEN
4871 DO k = bo(1, 3), bo(2, 3)
4872 DO j = bo(1, 2), bo(2, 2)
4873 DO i = bo(1, 1), bo(2, 1)
4874 IF (nspins /= 1)
THEN
4875 v_drho_r(idir, 1)%array(i, j, k) = v_drho_r(idir, 1)%array(i, j, k) - &
4876 e_drho(i, j, k)*drho1(idir)%array(i, j, k)
4877 v_drho_r(idir, 2)%array(i, j, k) = v_drho_r(idir, 2)%array(i, j, k) - &
4878 e_drho(i, j, k)*drho1(idir)%array(i, j, k)
4880 v_drho_r(idir, 1)%array(i, j, k) = v_drho_r(idir, 1)%array(i, j, k) - &
4881 e_drho(i, j, k)*(drho1a(idir)%array(i, j, k) + &
4882 fac*drho1b(idir)%array(i, j, k))
4892 DO ispin = 1, nspins
4893 CALL xc_pw_divergence(xc_deriv_method_id, v_drho_r(:, ispin), tmp_g, vxc_g, v_xc(ispin))
4901 DEALLOCATE (drho(idir)%array)
4902 DEALLOCATE (drho1(idir)%array)
4905 DO ispin = 1, nspins
4906 CALL deallocate_pw(v_drhoa(ispin), pw_pool)
4907 CALL deallocate_pw(v_drhob(ispin), pw_pool)
4910 DEALLOCATE (v_drhoa, v_drhob)
4914 IF (laplace_f .AND. my_compute_virial)
THEN
4915 virial_pw%array(:, :, :) = -rhoa(:, :, :)
4916 CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(1)%array)
4917 virial_pw%array(:, :, :) = -rhob(:, :, :)
4918 CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(2)%array)
4929 IF (gradient_f)
THEN
4932 CALL prepare_dr1dr(dr1dr, drho, drho1)
4938 ALLOCATE (v_laplace(nspins))
4939 DO ispin = 1, nspins
4940 CALL allocate_pw(v_laplace(ispin), pw_pool, bo)
4951 IF (
ASSOCIATED(deriv_att))
THEN
4955 DO k = bo(1, 3), bo(2, 3)
4956 DO j = bo(1, 2), bo(2, 2)
4957 DO i = bo(1, 1), bo(2, 1)
4958 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4959 deriv_data(i, j, k)*rho1(i, j, k)
4965 IF (
ASSOCIATED(deriv_att))
THEN
4969 DO k = bo(1, 3), bo(2, 3)
4970 DO j = bo(1, 2), bo(2, 2)
4971 DO i = bo(1, 1), bo(2, 1)
4972 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4973 deriv_data(i, j, k)*dr1dr(i, j, k)
4979 IF (
ASSOCIATED(deriv_att))
THEN
4983 DO k = bo(1, 3), bo(2, 3)
4984 DO j = bo(1, 2), bo(2, 2)
4985 DO i = bo(1, 1), bo(2, 1)
4986 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
4987 deriv_data(i, j, k)*tau1(i, j, k)
4993 IF (
ASSOCIATED(deriv_att))
THEN
4997 DO k = bo(1, 3), bo(2, 3)
4998 DO j = bo(1, 2), bo(2, 2)
4999 DO i = bo(1, 1), bo(2, 1)
5000 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
5001 deriv_data(i, j, k)*laplace1(i, j, k)
5009 IF (
ASSOCIATED(deriv_att))
THEN
5013 DO k = bo(1, 3), bo(2, 3)
5014 DO j = bo(1, 2), bo(2, 2)
5015 DO i = bo(1, 1), bo(2, 1)
5016 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
5017 deriv_data(i, j, k)*rho1(i, j, k)
5023 IF (
ASSOCIATED(deriv_att))
THEN
5027 DO k = bo(1, 3), bo(2, 3)
5028 DO j = bo(1, 2), bo(2, 2)
5029 DO i = bo(1, 1), bo(2, 1)
5030 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
5031 deriv_data(i, j, k)*dr1dr(i, j, k)
5037 IF (
ASSOCIATED(deriv_att))
THEN
5041 DO k = bo(1, 3), bo(2, 3)
5042 DO j = bo(1, 2), bo(2, 2)
5043 DO i = bo(1, 1), bo(2, 1)
5044 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
5045 deriv_data(i, j, k)*tau1(i, j, k)
5051 IF (
ASSOCIATED(deriv_att))
THEN
5055 DO k = bo(1, 3), bo(2, 3)
5056 DO j = bo(1, 2), bo(2, 2)
5057 DO i = bo(1, 1), bo(2, 1)
5058 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
5059 deriv_data(i, j, k)*laplace1(i, j, k)
5066 IF (
ASSOCIATED(deriv_att))
THEN
5070 IF (my_compute_virial)
THEN
5071 CALL virial_drho_drho1(virial_pw, drho, drho1, deriv_data, virial_xc)
5075 v_drho(1)%array(:, :, :) = v_drho(1)%array(:, :, :) + &
5076 deriv_data(:, :, :)*dr1dr(:, :, :)/max(gradient_cut, norm_drho(:, :, :))**2
5081 IF (
ASSOCIATED(deriv_att))
THEN
5085 DO k = bo(1, 3), bo(2, 3)
5086 DO j = bo(1, 2), bo(2, 2)
5087 DO i = bo(1, 1), bo(2, 1)
5088 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
5089 deriv_data(i, j, k)*rho1(i, j, k)
5095 IF (
ASSOCIATED(deriv_att))
THEN
5099 DO k = bo(1, 3), bo(2, 3)
5100 DO j = bo(1, 2), bo(2, 2)
5101 DO i = bo(1, 1), bo(2, 1)
5102 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
5103 deriv_data(i, j, k)*dr1dr(i, j, k)
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_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
5117 deriv_data(i, j, k)*tau1(i, j, k)
5123 IF (
ASSOCIATED(deriv_att))
THEN
5127 DO k = bo(1, 3), bo(2, 3)
5128 DO j = bo(1, 2), bo(2, 2)
5129 DO i = bo(1, 1), bo(2, 1)
5130 v_xc_tau(1)%array(i, j, k) = v_xc_tau(1)%array(i, j, k) + &
5131 deriv_data(i, j, k)*laplace1(i, j, k)
5139 IF (
ASSOCIATED(deriv_att))
THEN
5143 DO k = bo(1, 3), bo(2, 3)
5144 DO j = bo(1, 2), bo(2, 2)
5145 DO i = bo(1, 1), bo(2, 1)
5146 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
5147 deriv_data(i, j, k)*rho1(i, j, k)
5153 IF (
ASSOCIATED(deriv_att))
THEN
5157 DO k = bo(1, 3), bo(2, 3)
5158 DO j = bo(1, 2), bo(2, 2)
5159 DO i = bo(1, 1), bo(2, 1)
5160 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
5161 deriv_data(i, j, k)*dr1dr(i, j, k)
5167 IF (
ASSOCIATED(deriv_att))
THEN
5171 DO k = bo(1, 3), bo(2, 3)
5172 DO j = bo(1, 2), bo(2, 2)
5173 DO i = bo(1, 1), bo(2, 1)
5174 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
5175 deriv_data(i, j, k)*tau1(i, j, k)
5181 IF (
ASSOCIATED(deriv_att))
THEN
5185 DO k = bo(1, 3), bo(2, 3)
5186 DO j = bo(1, 2), bo(2, 2)
5187 DO i = bo(1, 1), bo(2, 1)
5188 v_laplace(1)%array(i, j, k) = v_laplace(1)%array(i, j, k) + &
5189 deriv_data(i, j, k)*laplace1(i, j, k)
5196 IF (my_compute_virial)
THEN
5198 IF (
ASSOCIATED(deriv_att))
THEN
5201 virial_pw%array(:, :, :) = -rho1(:, :, :)
5202 CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
5207 IF (gradient_f)
THEN
5209 IF (my_compute_virial)
THEN
5210 CALL virial_drho_drho(virial_pw, drho, v_drho(1), virial_xc)
5220 DO ia = bo(1, 1), bo(2, 1)
5221 DO ir = bo(1, 2), bo(2, 2)
5222 vxg(idir, ia, ir, 1) = -drho(idir)%array(ia, ir, 1)*v_drho(1)%array(ia, ir, 1)
5223 IF (
ASSOCIATED(e_drho))
THEN
5224 vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + factor2*drho1(idir)%array(ia, ir, 1)*e_drho(ia, ir, 1)
5235 v_drho_r(idir, 1)%array(:, :, :) = drho(idir)%array(:, :, :)*v_drho(1)%array(:, :, :) - &
5236 drho1(idir)%array(:, :, :)*e_drho(:, :, :)
5240 CALL xc_pw_divergence(xc_deriv_method_id, v_drho_r(:, 1), tmp_g, vxc_g, v_xc(1))
5245 IF (laplace_f .AND. my_compute_virial)
THEN
5246 virial_pw%array(:, :, :) = -rho(:, :, :)
5247 CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(1)%array)
5253 DO ispin = 1, nspins
5254 CALL xc_pw_laplace(v_laplace(ispin), pw_pool, xc_deriv_method_id)
5255 CALL pw_axpy(v_laplace(ispin), v_xc(ispin))
5259 IF (gradient_f)
THEN
5261 DO ispin = 1, nspins
5262 CALL deallocate_pw(v_drho(ispin), pw_pool)
5264 CALL deallocate_pw(v_drho_r(idir, ispin), pw_pool)
5267 DEALLOCATE (v_drho, v_drho_r)
5272 DO ispin = 1, nspins
5273 CALL deallocate_pw(v_laplace(ispin), pw_pool)
5275 DEALLOCATE (v_laplace)
5278 IF (
ASSOCIATED(tmp_g%pw_grid) .AND.
ASSOCIATED(pw_pool))
THEN
5279 CALL pw_pool%give_back_pw(tmp_g)
5282 IF (
ASSOCIATED(vxc_g%pw_grid) .AND.
ASSOCIATED(pw_pool))
THEN
5283 CALL pw_pool%give_back_pw(vxc_g)
5286 IF (my_compute_virial .AND. (gradient_f .OR. laplace_f))
THEN
5287 CALL deallocate_pw(virial_pw, pw_pool)
5290 CALL timestop(handle)
5326 LOGICAL,
INTENT(in),
OPTIONAL :: spinflip
5328 CHARACTER(len=*),
PARAMETER :: routinen =
'xc_calc_3rd_deriv_analytical'
5330 INTEGER :: handle, i, idir, ispin, j, &
5331 k, nspins, xc_deriv_method_id
5332 INTEGER,
DIMENSION(2, 3) :: bo
5333 LOGICAL :: lsd, do_spinflip, alda0, &
5334 rho_f, gradient_f, tau_f, laplace_f
5335 REAL(kind=
dp) :: s, s_thresh, s_thresh2, gradient_cut
5336 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: dr1dr, dra1dra, drb1drb
5337 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: deriv_data, deriv_data2, e_drhoa, e_drhob, &
5338 e_drho, norm_drho, norm_drhoa, &
5339 norm_drhob, rho1a, rho1b, &
5341 TYPE(
cp_3d_r_cp_type),
DIMENSION(3) :: drho, drho1, drho1a, drho1b, drhoa, drhob
5342 TYPE(
pw_r3d_rs_type),
DIMENSION(:),
ALLOCATABLE :: v_drhoa, v_drhob, v_drho
5347 CALL timeset(routinen, handle)
5349 NULLIFY (e_drhoa, e_drhob, e_drho)
5351 cpassert(
ASSOCIATED(v_xc))
5352 cpassert(
ASSOCIATED(xc_section))
5356 i_val=xc_deriv_method_id)
5359 lsd =
ASSOCIATED(rho_set%rhoa)
5361 do_spinflip = .false.
5362 IF (
PRESENT(spinflip)) do_spinflip = spinflip
5364 bo = rho_set%local_bounds
5366 CALL check_for_derivatives(deriv_set, lsd, rho_f, gradient_f, tau_f, laplace_f)
5376 DO ispin = 1, nspins
5378 v_xc(ispin)%array = 0.0_dp
5382 IF (gradient_f)
THEN
5383 ALLOCATE (v_drho_r(3, nspins), v_drho(nspins))
5384 DO ispin = 1, nspins
5386 CALL allocate_pw(v_drho_r(idir, ispin), pw_pool, bo)
5388 CALL allocate_pw(v_drho(ispin), pw_pool, bo)
5392 IF (
ASSOCIATED(pw_pool))
THEN
5393 CALL pw_pool%create_pw(tmp_g)
5394 CALL pw_pool%create_pw(vxc_g)
5397 cpabort(
"XC_DERIV method is not implemented in GAPW")
5405 cpassert(
ASSOCIATED(v_xc_tau))
5406 DO ispin = 1, nspins
5407 v_xc_tau(ispin)%array = 0.0_dp
5417 IF (do_spinflip)
THEN
5424 IF (gradient_f)
THEN
5426 norm_drho=norm_drho, norm_drhoa=norm_drhoa, norm_drhob=norm_drhob)
5427 IF (do_spinflip)
THEN
5429 CALL calc_drho_from_a(drho1, drho1a)
5432 CALL calc_drho_from_ab(drho1, drho1a, drho1b)
5435 CALL calc_drho_from_ab(drho, drhoa, drhob)
5437 CALL prepare_dr1dr(dra1dra, drhoa, drho1a)
5438 IF (do_spinflip)
THEN
5439 CALL prepare_dr1dr(drb1drb, drhob, drho1a)
5440 CALL prepare_dr1dr(dr1dr, drho, drho1a)
5441 ELSE IF (nspins /= 1)
THEN
5442 CALL prepare_dr1dr(drb1drb, drhob, drho1b)
5443 CALL prepare_dr1dr(dr1dr, drho, drho1)
5445 cpabort(
"Exchange-correlation's third derivative for closed-shell not yet implemented")
5449 ALLOCATE (v_drhoa(nspins), v_drhob(nspins))
5450 DO ispin = 1, nspins
5451 CALL allocate_pw(v_drhoa(ispin), pw_pool, bo)
5452 CALL allocate_pw(v_drhob(ispin), pw_pool, bo)
5458 cpabort(
"Exchange-correlation's laplace analytic third derivative not implemented")
5462 cpabort(
"Exchange-correlation's mGGA analytic third derivative not implemented")
5465 IF (nspins /= 1)
THEN
5467 IF (.NOT. do_spinflip)
THEN
5469 cpabort(
"Exchange-correlation's analytic third derivative not implemented")
5481 IF (
ASSOCIATED(deriv_att))
THEN
5484 IF (
ASSOCIATED(deriv_att))
THEN
5488 DO k = bo(1, 3), bo(2, 3)
5489 DO j = bo(1, 2), bo(2, 2)
5490 DO i = bo(1, 1), bo(2, 1)
5491 s = rhoa(i, j, k) - rhob(i, j, k)
5492 s = -sign(max(s**2, s_thresh2), s)
5493 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
5494 (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)**2/s
5506 IF (.NOT. alda0)
THEN
5508 IF (
ASSOCIATED(deriv_att))
THEN
5511 IF (
ASSOCIATED(deriv_att))
THEN
5515 DO k = bo(1, 3), bo(2, 3)
5516 DO j = bo(1, 2), bo(2, 2)
5517 DO i = bo(1, 1), bo(2, 1)
5518 s = rhoa(i, j, k) - rhob(i, j, k)
5519 s = -sign(max(s**2, s_thresh2), s)
5520 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
5521 (deriv_data(i, j, k)*dra1dra(i, j, k) - deriv_data2(i, j, k)*drb1drb(i, j, k)) &
5533 DO k = bo(1, 3), bo(2, 3)
5534 DO j = bo(1, 2), bo(2, 2)
5535 DO i = bo(1, 1), bo(2, 1)
5536 v_xc(2)%array(i, j, k) = -v_xc(1)%array(i, j, k)
5549 IF (
ASSOCIATED(deriv_att))
THEN
5552 IF (
ASSOCIATED(deriv_att))
THEN
5556 DO k = bo(1, 3), bo(2, 3)
5557 DO j = bo(1, 2), bo(2, 2)
5558 DO i = bo(1, 1), bo(2, 1)
5559 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
5560 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
5561 rho1a(i, j, k)**2*(deriv_data(i, j, k) - deriv_data2(i, j, k))/s
5573 IF (
ASSOCIATED(deriv_att))
THEN
5576 IF (
ASSOCIATED(deriv_att))
THEN
5580 DO k = bo(1, 3), bo(2, 3)
5581 DO j = bo(1, 2), bo(2, 2)
5582 DO i = bo(1, 1), bo(2, 1)
5583 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
5584 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
5585 rho1a(i, j, k)**2*(deriv_data(i, j, k) - deriv_data2(i, j, k))/s
5593 IF (.NOT. alda0)
THEN
5598 IF (
ASSOCIATED(deriv_att))
THEN
5601 IF (
ASSOCIATED(deriv_att))
THEN
5605 DO k = bo(1, 3), bo(2, 3)
5606 DO j = bo(1, 2), bo(2, 2)
5607 DO i = bo(1, 1), bo(2, 1)
5608 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
5609 v_xc(1)%array(i, j, k) = v_xc(1)%array(i, j, k) + &
5610 (deriv_data(i, j, k)*dra1dra(i, j, k) - deriv_data2(i, j, k)*drb1drb(i, j, k))* &
5623 IF (
ASSOCIATED(deriv_att))
THEN
5626 IF (
ASSOCIATED(deriv_att))
THEN
5630 DO k = bo(1, 3), bo(2, 3)
5631 DO j = bo(1, 2), bo(2, 2)
5632 DO i = bo(1, 1), bo(2, 1)
5633 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
5634 v_xc(2)%array(i, j, k) = v_xc(2)%array(i, j, k) + &
5635 (deriv_data(i, j, k)*dra1dra(i, j, k) - deriv_data2(i, j, k)*drb1drb(i, j, k))* &
5651 IF (
ASSOCIATED(deriv_att))
THEN
5654 IF (
ASSOCIATED(deriv_att))
THEN
5658 DO k = bo(1, 3), bo(2, 3)
5659 DO j = bo(1, 2), bo(2, 2)
5660 DO i = bo(1, 1), bo(2, 1)
5661 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
5662 (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)
5674 IF (
ASSOCIATED(deriv_att))
THEN
5677 IF (
ASSOCIATED(deriv_att))
THEN
5681 DO k = bo(1, 3), bo(2, 3)
5682 DO j = bo(1, 2), bo(2, 2)
5683 DO i = bo(1, 1), bo(2, 1)
5684 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) + &
5685 (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)
5701 IF (
ASSOCIATED(deriv_att))
THEN
5704 IF (
ASSOCIATED(deriv_att))
THEN
5708 DO k = bo(1, 3), bo(2, 3)
5709 DO j = bo(1, 2), bo(2, 2)
5710 DO i = bo(1, 1), bo(2, 1)
5711 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
5712 (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)
5713 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
5714 (deriv_data(i, j, k) - deriv_data2(i, j, k))*rho1a(i, j, k)
5723 IF (
ASSOCIATED(deriv_att))
THEN
5729 IF (
ASSOCIATED(deriv_att))
THEN
5733 DO k = bo(1, 3), bo(2, 3)
5734 DO j = bo(1, 2), bo(2, 2)
5735 DO i = bo(1, 1), bo(2, 1)
5736 v_drhoa(1)%array(i, j, k) = v_drhoa(1)%array(i, j, k) - &
5737 deriv_data(i, j, k)*dra1dra(i, j, k) + deriv_data2(i, j, k)*drb1drb(i, j, k)
5747 IF (
ASSOCIATED(deriv_att))
THEN
5751 DO k = bo(1, 3), bo(2, 3)
5752 DO j = bo(1, 2), bo(2, 2)
5753 DO i = bo(1, 1), bo(2, 1)
5754 v_drhob(2)%array(i, j, k) = v_drhob(2)%array(i, j, k) - &
5755 deriv_data(i, j, k)*dra1dra(i, j, k) + deriv_data2(i, j, k)*drb1drb(i, j, k)
5770 IF (
ASSOCIATED(deriv_att))
THEN
5773 IF (
ASSOCIATED(deriv_att))
THEN
5777 DO k = bo(1, 3), bo(2, 3)
5778 DO j = bo(1, 2), bo(2, 2)
5779 DO i = bo(1, 1), bo(2, 1)
5780 v_drho(1)%array(i, j, k) = v_drho(1)%array(i, j, k) - &
5781 deriv_data(i, j, k)*dra1dra(i, j, k) + &
5782 deriv_data2(i, j, k)*drb1drb(i, j, k)
5783 v_drho(2)%array(i, j, k) = v_drho(2)%array(i, j, k) - &
5784 deriv_data(i, j, k)*dra1dra(i, j, k) + &
5785 deriv_data2(i, j, k)*drb1drb(i, j, k)
5800 IF (
ASSOCIATED(deriv_att))
THEN
5805 v_drhoa(1)%array(:, :, :) = v_drhoa(1)%array(:, :, :) + &
5806 deriv_data(:, :, :)*dra1dra(:, :, :)/max(gradient_cut, norm_drhoa(:, :, :))**2
5814 IF (
ASSOCIATED(deriv_att))
THEN
5819 v_drhob(2)%array(:, :, :) = v_drhob(2)%array(:, :, :) - &
5820 deriv_data(:, :, :)*drb1drb(:, :, :)/max(gradient_cut, norm_drhob(:, :, :))**2
5829 cpabort(
"Exchange-correlation's analytic third derivative not implemented")
5833 IF (gradient_f)
THEN
5834 IF (.NOT. alda0)
THEN
5848 IF (do_spinflip)
THEN
5851 DO k = bo(1, 3), bo(2, 3)
5852 DO j = bo(1, 2), bo(2, 2)
5853 DO i = bo(1, 1), bo(2, 1)
5854 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
5856 v_drho_r(idir, ispin)%array(i, j, k) = v_drho_r(idir, ispin)%array(i, j, k) + &
5857 v_drhoa(ispin)%array(i, j, k)*drhoa(idir)%array(i, j, k)*rho1a(i, j, k)/s + &
5858 v_drhob(ispin)%array(i, j, k)*drhob(idir)%array(i, j, k)*rho1a(i, j, k)/s + &
5859 v_drho(ispin)%array(i, j, k)*drho(idir)%array(i, j, k)*rho1a(i, j, k)/s
5871 IF (
ASSOCIATED(e_drhoa))
THEN
5874 DO k = bo(1, 3), bo(2, 3)
5875 DO j = bo(1, 2), bo(2, 2)
5876 DO i = bo(1, 1), bo(2, 1)
5877 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
5878 v_drho_r(idir, 1)%array(i, j, k) = v_drho_r(idir, 1)%array(i, j, k) - &
5879 e_drhoa(i, j, k)*drho1a(idir)%array(i, j, k)*rho1a(i, j, k)/s
5889 IF (
ASSOCIATED(e_drhob))
THEN
5892 DO k = bo(1, 3), bo(2, 3)
5893 DO j = bo(1, 2), bo(2, 2)
5894 DO i = bo(1, 1), bo(2, 1)
5895 s = max(abs(rhoa(i, j, k) - rhob(i, j, k)), s_thresh)
5896 v_drho_r(idir, 2)%array(i, j, k) = v_drho_r(idir, 2)%array(i, j, k) + &
5897 e_drhob(i, j, k)*drho1a(idir)%array(i, j, k)*rho1a(i, j, k)/s
5906 DO ispin = 1, nspins
5907 CALL xc_pw_divergence(xc_deriv_method_id, v_drho_r(:, ispin), tmp_g, vxc_g, v_xc(ispin))
5912 DEALLOCATE (drho(idir)%array)
5913 DEALLOCATE (drho1(idir)%array)
5916 DO ispin = 1, nspins
5917 CALL deallocate_pw(v_drhoa(ispin), pw_pool)
5918 CALL deallocate_pw(v_drhob(ispin), pw_pool)
5921 DEALLOCATE (v_drhoa, v_drhob)
5930 cpabort(
"Exchange-correlation's analytic third derivative not implemented")
5934 IF (gradient_f)
THEN
5936 DO ispin = 1, nspins
5937 CALL deallocate_pw(v_drho(ispin), pw_pool)
5939 CALL deallocate_pw(v_drho_r(idir, ispin), pw_pool)
5942 DEALLOCATE (v_drho, v_drho_r)
5946 IF (
ASSOCIATED(tmp_g%pw_grid) .AND.
ASSOCIATED(pw_pool))
THEN
5947 CALL pw_pool%give_back_pw(tmp_g)
5950 IF (
ASSOCIATED(vxc_g%pw_grid) .AND.
ASSOCIATED(pw_pool))
THEN
5951 CALL pw_pool%give_back_pw(vxc_g)
5954 CALL timestop(handle)