diff --git a/README.md b/README.md index 6136a48..01c6df1 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,6 @@ # Py_xDH + Python approach of realization to xDH type of functional (written in Chinese) + +Published web page: https://py-xdh.readthedocs.io/zh_CN/latest/ + diff --git a/source/dft_nuc_hess.ipynb b/source/dft_nuc_hess.ipynb index 990110f..a1236e7 100644 --- a/source/dft_nuc_hess.ipynb +++ b/source/dft_nuc_hess.ipynb @@ -22,7 +22,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -50,9 +50,20 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "mol = gto.Mole()\n", "mol.atom = \"\"\"\n", @@ -67,9 +78,17 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "converged SCF energy = -151.477320518606\n" + ] + } + ], "source": [ "scf_eng = dft.RKS(mol)\n", "scf_eng.xc = \"b3lypg\" # compare that to gaussian\n", @@ -81,9 +100,23 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "--------------- RKS gradients ---------------\n", + " x y z\n", + "0 O -0.0271975153 0.0107158349 0.0176473933\n", + "1 O 0.0107158349 -0.0271975153 -0.0176473933\n", + "2 H 0.0099981076 0.0064834995 0.0257989619\n", + "3 H 0.0064834995 0.0099981076 -0.0257989619\n", + "----------------------------------------------\n" + ] + } + ], "source": [ "scf_grad = grad.rks.Gradients(scf_eng)\n", "grad_RKS = scf_grad.kernel()" @@ -91,7 +124,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ @@ -99,6 +132,13 @@ "hess_RKS = scf_hess.kernel()" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, { "cell_type": "markdown", "metadata": {}, @@ -108,7 +148,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "metadata": {}, "outputs": [], "source": [ @@ -128,7 +168,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 7, "metadata": {}, "outputs": [], "source": [ @@ -146,9 +186,19 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "True\n", + "True\n", + "True\n" + ] + } + ], "source": [ "print(np.allclose(energy_RKS, energy_Gaussian))\n", "print(np.allclose(grad_RKS.reshape(-1), grad_Gaussian, atol=1e-5))\n", @@ -165,7 +215,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 9, "metadata": {}, "outputs": [], "source": [ @@ -191,7 +241,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 10, "metadata": {}, "outputs": [], "source": [ @@ -217,7 +267,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 11, "metadata": {}, "outputs": [], "source": [ @@ -249,7 +299,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 12, "metadata": {}, "outputs": [], "source": [ @@ -390,7 +440,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 13, "metadata": {}, "outputs": [], "source": [ @@ -438,9 +488,20 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "np.allclose(grad_hf + grad_gga, scf_grad.grad_elec())" ] @@ -461,7 +522,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 15, "metadata": {}, "outputs": [], "source": [ @@ -471,7 +532,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 16, "metadata": {}, "outputs": [], "source": [ @@ -496,7 +557,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 17, "metadata": {}, "outputs": [], "source": [ @@ -508,7 +569,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 18, "metadata": {}, "outputs": [], "source": [ @@ -526,7 +587,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 19, "metadata": {}, "outputs": [], "source": [ @@ -540,7 +601,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 20, "metadata": {}, "outputs": [], "source": [ @@ -549,16 +610,27 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "execution_count": 21, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 21, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "np.allclose(vxc_diag_pyscf, vxc_diag)" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 22, "metadata": {}, "outputs": [], "source": [ @@ -577,7 +649,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 23, "metadata": {}, "outputs": [], "source": [ @@ -586,54 +658,267 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 24, "metadata": {}, "outputs": [], "source": [ - "from pyscf.dft import numint" + "from pyscf.dft import numint\n", + "import time" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 25, "metadata": {}, "outputs": [], "source": [ - "%%time\n", + "vxc_deriv2_pyscf = hessian.rks._get_vxc_deriv2(scf_hess, C, mo_occ, 16000)\n", + "hess_deriv2_pyscf = np.zeros((natm, natm, 3, 3))\n", + "for A in range(natm):\n", + " for B in range(natm):\n", + " sA, sB = mol_slice(A), mol_slice(B)\n", + " hess_deriv2_pyscf[A, B] = np.einsum(\"tsuv, uv -> ts\", vxc_deriv2_pyscf[A, :, :, sB], D[sB]) * 2" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Time" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "3.644152879714966\n", + "9.904576301574707\n", + "10.93254017829895\n", + "0.004670143127441406\n" + ] + } + ], + "source": [ + "# %%time\n", "\n", + "#---\n", + "time0 = time.time();\n", "vmat = np.zeros((mol.natm,3,3,nao,nao))\n", - "ipip = np.zeros((3,3,nao,nao))\n", + "dR_rho1 = np.zeros((3, 4, ngrid))\n", + "dR_rho1_0 = 2 * np.einsum(\"tgk , gl , li -> tgki \", grid_ao_1, grid_ao_0, Co, optimize=True)\n", + "dR_rho1_1 = 2 * np.einsum(\"tsgk, gl , li -> tsgki\", grid_ao_2, grid_ao_0, Co, optimize=True)\n", + "dR_rho1_1 += 2 * np.einsum(\"tgk , sgl, li -> tsgki\", grid_ao_1, grid_ao_1, Co, optimize=True)\n", + "print(time.time() - time0)\n", "\n", - "ni = scf_eng._numint\n", - "xctype = ni._xc_type(scf_eng.xc)\n", - "aoslices = mol.aoslice_by_atom()\n", - "shls_slice = (0, mol.nbas)\n", - "ao_loc = mol.ao_loc_nr()\n", + "#---\n", + "time0 = time.time();\n", + "wv_0 = np.zeros((3, ngrid, nao, nocc))\n", + "wv_1 = np.zeros((3, 3, ngrid, nao, nocc))\n", + "wv_0 += .5* np.einsum(\"g, g, tgki -> tgki \", grid_weight, grid_frr, dR_rho1_0, optimize=True)\n", + "wv_1 += 2 * np.einsum(\"g, g, tgki , rg -> trgki\", grid_weight, grid_frg, dR_rho1_0, grid_rho_1, optimize=True)\n", + "wv_0 += np.einsum(\"g, g, trgki, rg -> tgki \", grid_weight, grid_frg, dR_rho1_1, grid_rho_1, optimize=True)\n", + "wv_1 += 4 * np.einsum(\"g, g, tsgki, sg, rg -> trgki\", grid_weight, grid_fgg, dR_rho1_1, grid_rho_1, grid_rho_1, optimize=True)\n", + "wv_1 += 2 * np.einsum(\"g, g, trgki -> trgki\", grid_weight, grid_fg, dR_rho1_1, optimize=True)\n", + "print(time.time() - time0)\n", "\n", + "#---\n", + "time0 = time.time();\n", + "vmat = np.einsum(\"tgki, sgu, gv , vj -> tskuij\", wv_0, grid_ao_1, grid_ao_0, Co, optimize=True)\n", + "vmat += np.einsum(\"tgki, gv, sgu, vj -> tskuij\", wv_0, grid_ao_0, grid_ao_1, Co, optimize=True)\n", + "vmat += np.einsum(\"trgki, srgu, gv , vj -> tskuij\", wv_1, grid_ao_2, grid_ao_0, Co, optimize=True)\n", + "vmat += np.einsum(\"trgki, rgv, sgu, vj -> tskuij\", wv_1, grid_ao_1, grid_ao_1, Co, optimize=True)\n", + "print(time.time() - time0)\n", "\n", - "ipip += 0.5 * np.einsum(\"g, g, sgu, tgv -> tsuv\", grid_weight, grid_fr, grid_ao_1, grid_ao_1, optimize=True)\n", + "#---\n", + "time0 = time.time();\n", + "hess_vmat_tmp = np.zeros((natm, natm, 3, 3))\n", + "for A in range(natm):\n", + " for B in range(natm):\n", + " sA, sB = mol_slice(A), mol_slice(B)\n", + " hess_vmat_tmp[A, B] = 8 * np.einsum(\"tskuij, ki, uj -> ts\", vmat[:, :, sA, sB], Co[sA], Co[sB])\n", + "print(time.time() - time0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 53, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.3527677059173584\n", + "0.12494540214538574\n", + "1.4512510299682617\n" + ] + } + ], + "source": [ + "timeA = timeB = timeC = 0\n", + "hess_deriv2 = np.zeros((natm, natm, 3, 3))\n", + "\n", + "for A in range(natm):\n", + " sA = mol_slice(A)\n", + "\n", + " #---\n", + " time0 = time.time();\n", + " grid_A_rho_1 = 2 * np.einsum(\"tgk , gl , kl -> tg \", grid_ao_1[:, :, sA], grid_ao_0, D[sA], optimize=True)\n", + " grid_A_rho_2 = 2 * np.einsum(\"trgk, gl , kl -> trg\", grid_ao_2[:, :, :, sA], grid_ao_0, D[sA], optimize=True)\n", + " grid_A_rho_2 += 2 * np.einsum(\"tgk , rgl, kl -> trg\", grid_ao_1[:, :, sA], grid_ao_1, D[sA], optimize=True)\n", + " timeA += time.time() - time0\n", + "\n", + " #---\n", + " time0 = time.time();\n", + " wv_0 = np.zeros((3, ngrid))\n", + " wv_1 = np.zeros((3, 3, ngrid))\n", + " wv_0 += .5* np.einsum(\"g, g, tg -> tg \", grid_weight, grid_frr, grid_A_rho_1, optimize=True)\n", + " wv_1 += 2 * np.einsum(\"g, g, tg , rg -> trg\", grid_weight, grid_frg, grid_A_rho_1, grid_rho_1, optimize=True)\n", + " wv_0 += np.einsum(\"g, g, trg, rg -> tg \", grid_weight, grid_frg, grid_A_rho_2, grid_rho_1, optimize=True)\n", + " wv_1 += 4 * np.einsum(\"g, g, twg, rg, wg -> trg\", grid_weight, grid_fgg, grid_A_rho_2, grid_rho_1, grid_rho_1, optimize=True)\n", + " wv_1 += 2 * np.einsum(\"g, g, trg -> trg\", grid_weight, grid_fg, grid_A_rho_2, optimize=True)\n", + " timeB += time.time() - time0\n", + "\n", + " #---\n", + " time0 = time.time();\n", + " for B in range(natm):\n", + " sB = mol_slice(B)\n", + " hess_deriv2[A, B] = 4 * np.einsum(\"tg, sgu, gv , uv -> ts\", wv_0, grid_ao_1[:, :, sB], grid_ao_0, D[sB], optimize=True) \n", + " hess_deriv2[A, B] += 2 * np.einsum(\"trg, srgu, gv , uv -> ts\", wv_1, grid_ao_2[:, :, :, sB], grid_ao_0, D[sB], optimize=True)\n", + " hess_deriv2[A, B] += 2 * np.einsum(\"trg, sgu, rgv, uv -> ts\", wv_1, grid_ao_1[:, :, sB], grid_ao_1, D[sB], optimize=True)\n", + " hess_deriv2[B, A] = hess_deriv2[A, B].T\n", + " timeC += time.time() - time0\n", + "\n", + "#---\n", + " \n", + "print(timeA)\n", + "print(timeB)\n", + "print(timeC)" + ] + }, + { + "cell_type": "code", + "execution_count": 54, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 54, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ipip = 0.5 * np.einsum(\"g, g, sgu, tgv -> tsuv\", grid_weight, grid_fr, grid_ao_1, grid_ao_1, optimize=True)\n", "ipip += 2 * np.einsum(\"g, g, rg, rsgu, tgv -> tsuv\", grid_weight, grid_fg, grid_rho_1, grid_ao_2, grid_ao_1, optimize=True)\n", "\n", - "for ia in range(natm):\n", - " sA = mol_slice(ia)\n", - " \n", - " wv = np.zeros((3, 4, ngrid))\n", - " dR_rho1 = np.zeros((3, 4, ngrid))\n", - " dR_rho1[:, 0] = 2 * np.einsum(\"tgu, uv, gv -> tg\", grid_ao_1[:, :, sA], D[sA], grid_ao_0, optimize=True)\n", - " dR_rho1[:, 1:4] = 2 * np.einsum(\"tsgu, uv, gv -> tsg\", grid_ao_2[:, :, :, sA], D[sA], grid_ao_0, optimize=True)\n", - " dR_rho1[:, 1:4] += 2 * np.einsum(\"tgu, uv, sgv -> tsg\", grid_ao_1[:, :, sA], D[sA], grid_ao_1, optimize=True)\n", - " \n", - " wv[:, 0] = 0.5 * np.einsum(\"g, g, tg -> tg\", grid_weight, grid_frr, dR_rho1[:, 0], optimize=True)\n", - " wv[:, 0] += np.einsum(\"g, g, trg, rg -> tg\", grid_weight, grid_frg, dR_rho1[:, 1:4], grid_rho_1, optimize=True)\n", - " wv[:, 1:4] = 2 * np.einsum(\"g, g, tg, sg -> tsg\", grid_weight, grid_frg, dR_rho1[:, 0], grid_rho_1, optimize=True)\n", - " wv[:, 1:4] += 4 * np.einsum(\"g, g, rg, trg, sg -> tsg\", grid_weight, grid_fgg, grid_rho_1, dR_rho1[:, 1:4], grid_rho_1, optimize=True)\n", - " wv[:, 1:4] += 2 * np.einsum(\"g, g, tsg -> tsg\", grid_weight, grid_fg, dR_rho1[:, 1:4], optimize=True)\n", - " \n", - " vmat[ia] += np.einsum(\"tg, sgu, gv -> tsuv\", wv[:, 0], grid_ao_1, grid_ao_0, optimize=True)\n", - " vmat[ia] += np.einsum(\"trg, srgu, gv -> tsuv\", wv[:, 1:4], grid_ao_2, grid_ao_0, optimize=True)\n", - " vmat[ia] += np.einsum(\"trg, rgv, sgu -> tsuv\", wv, grid_ao[:4], grid_ao_1, optimize=True)\n", - " \n", - "print(np.allclose(vmat, vxc_deriv2_pyscf))" + "for A in range(natm):\n", + " for B in range(natm):\n", + " sA, sB = mol_slice(A), mol_slice(B)\n", + " hess_deriv2[A, B] += np.einsum(\"tsuv, uv -> ts\", ipip[:, :, sB, sA], D[sB, sA]) * 2\n", + " hess_deriv2[A, B] += np.einsum(\"stvu, uv -> ts\", ipip[:, :, sA, sB], D[sB, sA]) * 2\n", + "np.allclose(hess_deriv2, hess_deriv2_pyscf)" + ] + }, + { + "cell_type": "code", + "execution_count": 49, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[[[ 0. , 0. , -0. ],\n", + " [-0. , 0. , -0. ],\n", + " [-0. , 0. , 0. ]],\n", + "\n", + " [[-0. , -0.00236, 0.03204],\n", + " [ 0.00236, 0. , 0.03204],\n", + " [-0.03204, -0.03204, -0. ]],\n", + "\n", + " [[-0. , -0.00206, -0.022 ],\n", + " [ 0.00206, 0. , 0.00012],\n", + " [ 0.022 , -0.00012, 0. ]],\n", + "\n", + " [[-0. , 0.00128, 0.00187],\n", + " [-0.00128, 0. , -0.00794],\n", + " [-0.00187, 0.00794, -0. ]]],\n", + "\n", + "\n", + " [[[-0. , -0. , 0. ],\n", + " [ 0. , 0. , 0. ],\n", + " [-0. , 0. , -0. ]],\n", + "\n", + " [[ 0. , -0. , 0. ],\n", + " [ 0. , 0. , 0. ],\n", + " [-0. , -0. , 0. ]],\n", + "\n", + " [[-0. , -0.00128, 0.00794],\n", + " [ 0.00128, 0. , -0.00187],\n", + " [-0.00794, 0.00187, -0. ]],\n", + "\n", + " [[-0. , 0.00206, -0.00012],\n", + " [-0.00206, 0. , 0.022 ],\n", + " [ 0.00012, -0.022 , 0. ]]],\n", + "\n", + "\n", + " [[[-0. , 0. , 0. ],\n", + " [-0. , 0. , 0. ],\n", + " [-0. , 0. , 0. ]],\n", + "\n", + " [[-0. , 0. , 0. ],\n", + " [-0. , 0. , 0. ],\n", + " [ 0. , -0. , 0. ]],\n", + "\n", + " [[ 0. , -0. , -0. ],\n", + " [-0. , 0. , 0. ],\n", + " [-0. , 0. , 0. ]],\n", + "\n", + " [[-0. , -0.00012, -0.00007],\n", + " [ 0.00012, -0. , -0.00007],\n", + " [ 0.00007, 0.00007, -0. ]]],\n", + "\n", + "\n", + " [[[-0. , 0. , -0. ],\n", + " [ 0. , -0. , -0. ],\n", + " [ 0. , -0. , 0. ]],\n", + "\n", + " [[-0. , -0. , -0. ],\n", + " [-0. , 0. , -0. ],\n", + " [ 0. , -0. , 0. ]],\n", + "\n", + " [[-0. , 0. , -0. ],\n", + " [ 0. , -0. , 0. ],\n", + " [-0. , 0. , -0. ]],\n", + "\n", + " [[ 0. , -0. , 0. ],\n", + " [ 0. , 0. , 0. ],\n", + " [-0. , 0. , -0. ]]]])" + ] + }, + "execution_count": 49, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "hess_deriv2 - hess_deriv2_pyscf" ] }, { @@ -655,76 +940,84 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "vmat = np.zeros((mol.natm,3,3,nao,nao))\n", - "\n", - "for ia in range(natm):\n", - " sA = mol_slice(ia)\n", - " \n", - " wv = np.zeros((3, 4, ngrid))\n", - " dR_rho1 = np.zeros((3, 4, ngrid))\n", - " dR_rho1[:, 0] = 2 * np.einsum(\"tgu, uv, gv -> tg\", grid_ao_1[:, :, sA], D[sA], grid_ao_0, optimize=True)\n", - " dR_rho1[:, 1:4] = 2 * np.einsum(\"tsgu, uv, gv -> tsg\", grid_ao_2[:, :, :, sA], D[sA], grid_ao_0, optimize=True)\n", - " dR_rho1[:, 1:4] += 2 * np.einsum(\"tgu, uv, sgv -> tsg\", grid_ao_1[:, :, sA], D[sA], grid_ao_1, optimize=True)\n", - " \n", - " wv[:, 0] = 0.5 * np.einsum(\"g, g, tg -> tg\", grid_weight, grid_frr, dR_rho1[:, 0], optimize=True)\n", - " wv[:, 0] += np.einsum(\"g, g, trg, rg -> tg\", grid_weight, grid_frg, dR_rho1[:, 1:4], grid_rho_1, optimize=True)\n", - " wv[:, 1:4] = 2 * np.einsum(\"g, g, tg, sg -> tsg\", grid_weight, grid_frg, dR_rho1[:, 0], grid_rho_1, optimize=True)\n", - " wv[:, 1:4] += 4 * np.einsum(\"g, g, rg, trg, sg -> tsg\", grid_weight, grid_fgg, grid_rho_1, dR_rho1[:, 1:4], grid_rho_1, optimize=True)\n", - " wv[:, 1:4] += 2 * np.einsum(\"g, g, tsg -> tsg\", grid_weight, grid_fg, dR_rho1[:, 1:4], optimize=True)\n", - " \n", - " vmat[ia] += np.einsum(\"tg, sgu, gv -> tsuv\", wv[:, 0], grid_ao_1, grid_ao_0, optimize=True)\n", - " vmat[ia] += np.einsum(\"trg, srgu, gv -> tsuv\", wv[:, 1:4], grid_ao_2, grid_ao_0, optimize=True)\n", - " vmat[ia] += np.einsum(\"trg, rgv, sgu -> tsuv\", wv, grid_ao[:4], grid_ao_1, optimize=True)\n", - " \n", - "hess_vmat = np.zeros((natm, natm, 3, 3))\n", - "for A in range(natm):\n", - " for B in range(natm):\n", - " sA, sB = mol_slice(A), mol_slice(B)\n", - " hess_vmat[A, B] = np.einsum(\"tsuv, uv -> ts\", vmat[A, :, :, sB], D[sB]) * 2" - ] + "source": [] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 32, "metadata": {}, "outputs": [], "source": [ - "%%time\n", - "\n", - "vmat = np.zeros((mol.natm,3,3,nao,nao))\n", - "\n", - "sA = mol_slice(ia)\n", + "def get_hess_ao_noU_hcore(A, B):\n", + " ao_matrix = np.zeros((3 * 3, nao, nao))\n", + " zA, zB = mol.atom_charge(A), mol.atom_charge(B)\n", + " sA, sB = mol_slice(A), mol_slice(B)\n", + " if (A == B):\n", + " ao_matrix[:, sA] += int1e_ipipkin[:, sA]\n", + " ao_matrix[:, sA] += int1e_ipipnuc[:, sA]\n", + " with mol.with_rinv_as_nucleus(A):\n", + " ao_matrix -= zA * mol.intor('int1e_ipiprinv')\n", + " ao_matrix -= zA * mol.intor('int1e_iprinvip')\n", + " ao_matrix[:, sA, sB] += int1e_ipkinip[:, sA, sB]\n", + " ao_matrix[:, sA, sB] += int1e_ipnucip[:, sA, sB]\n", + " with mol.with_rinv_as_nucleus(B):\n", + " ao_matrix[:, sA] += zB * mol.intor('int1e_ipiprinv')[:, sA]\n", + " ao_matrix[:, sA] += zB * mol.intor('int1e_iprinvip')[:, sA]\n", + " with mol.with_rinv_as_nucleus(A):\n", + " ao_matrix[:, sB] += zA * mol.intor('int1e_ipiprinv')[:, sB]\n", + " ao_matrix[:, sB] += zA * mol.intor('int1e_iprinvip').swapaxes(1, 2)[:, sB]\n", + " ao_matrix += ao_matrix.swapaxes(1, 2)\n", + " return ao_matrix\n", "\n", - "wv = np.zeros((3, 4, ngrid))\n", - "dR_rho1 = np.zeros((3, 4, ngrid))\n", - "dR_rho1_0 = 2 * np.einsum(\"tgu, gv -> tguv\", grid_ao_1, grid_ao_0, optimize=True)\n", - "dR_rho1_1 = 2 * np.einsum(\"tsgu, gv -> tsguv\", grid_ao_2, grid_ao_0, optimize=True)\n", - "dR_rho1_1 += 2 * np.einsum(\"tgu, sgv -> tsguv\", grid_ao_1, grid_ao_1, optimize=True)\n", + "eri_mat1 = np.einsum(\"Tuvkl, kl -> Tuv\", int2e_ipip1, D) * 2 - c_x * np.einsum(\"Tukvl, kl -> Tuv\", int2e_ipip1, D)\n", + "eri_mat2 = np.einsum(\"Tuvkl, kl -> Tuv\", int2e_ipvip1, D) * 2 - c_x * np.einsum(\"Tukvl, kl -> Tuv\", int2e_ip1ip2, D)\n", + "eri_tensor1 = int2e_ip1ip2 * 4 - c_x * int2e_ipvip1.swapaxes(2, 3) - c_x * int2e_ip1ip2.swapaxes(2, 4)\n", "\n", - "wv_0 = np.zeros((3, ngrid, nao))\n", - "wv_1 = np.zeros((3, 3, ngrid, nao))\n", + "def get_hess_ao_noU_eri(A, B):\n", + " ao_matrix = np.zeros((9, nao, nao))\n", + " sA, sB = mol_slice(A), mol_slice(B)\n", + " if (A == B):\n", + " ao_matrix[:, sA] += eri_mat1[:, sA]\n", + " ao_matrix[:, sA, sB] += eri_mat2[:, sA, sB]\n", + " ao_matrix[:, sA] += np.einsum(\"Tuvkl, kl -> Tuv\", eri_tensor1[:, sA, :, sB], D[sB])\n", + " return ao_matrix\n", "\n", - "wv_0 += .5* np.einsum(\"g, g, tguv -> tguv\", grid_weight, grid_frr, dR_rho1_0, optimize=True)\n", - "wv_1 += 2 * np.einsum(\"g, g, tguv, sg -> tsguv\", grid_weight, grid_frg, dR_rho1_0, grid_rho_1, optimize=True)\n", - "wv_0 += np.einsum(\"g, g, trguv, rg -> tguv\", grid_weight, grid_frg, dR_rho1_1, grid_rho_1, optimize=True)\n", - "wv_1 += 4 * np.einsum(\"g, g, trguv, rg, sg -> tsguv\", grid_weight, grid_fgg, dR_rho1_1, grid_rho_1, grid_rho_1, optimize=True)\n", - "wv_1 += 2 * np.einsum(\"g, g, trguv -> trguv\", grid_weight, grid_fg, dR_rho1_1, optimize=True)\n", + "def get_hess_ao_noU_S(A, B):\n", + " ao_matrix = np.zeros((9, nao, nao))\n", + " sA, sB = mol_slice(A), mol_slice(B)\n", + " if (A == B):\n", + " ao_matrix[:, sA] -= int1e_ipipovlp[:, sA] * 2\n", + " ao_matrix[:, sA, sB] -= int1e_ipovlpip[:, sA, sB] * 2\n", + " return ao_matrix\n", "\n", - "for ia in range(natm):\n", - " \n", - " vmat[ia] = np.einsum(\"tgkl, sgu, gv , kl -> tsuv\", wv_0[:, :, sA], grid_ao_1, grid_ao_0, D[sA], optimize=True)\n", - " vmat[ia] += np.einsum(\"tgkl, gv, sgu , kl -> tsuv\", wv_0[:, :, sA], grid_ao_0, grid_ao_1, D[sA], optimize=True)\n", - " vmat[ia] += np.einsum(\"trgkl, srgu, gv , kl -> tsuv\", wv_1[:, :, :, sA], grid_ao_2, grid_ao_0, D[sA], optimize=True)\n", - " vmat[ia] += np.einsum(\"trgkl, rgv, sgu , kl -> tsuv\", wv_1[:, :, :, sA], grid_ao_1, grid_ao_1, D[sA], optimize=True)\n", + "def get_hess_noU(A, B):\n", + " return (np.einsum(\"Tuv, uv -> T\", get_hess_ao_noU_hcore(A, B) + get_hess_ao_noU_eri(A, B), D).reshape(3, 3)\n", + " + np.einsum(\"Tuv, uv -> T\", + get_hess_ao_noU_S(A, B), De).reshape(3, 3))\n", "\n", - " vmat[ia] += vmat_all[:, :, sA].sum(axis=2)\n", - " \n", - "hess_vmat_tmp = np.zeros((natm, natm, 3, 3))\n", - "for A in range(natm):\n", - " for B in range(natm):\n", - " sA, sB = mol_slice(A), mol_slice(B)\n", - " hess_vmat_tmp[A, B] = np.einsum(\"tsuv, uv -> ts\", vmat[A, :, :, sB], D[sB]) * 2" + "hess_noU_noGGA = np.array([ [ get_hess_noU(A, B) for B in range(natm) ] for A in range(natm) ])" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 33, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np.allclose(\n", + " hess_noU_noGGA + hess_veff_diag + hess_deriv2,\n", + " partial_hess_pyscf\n", + ")" ] }, { @@ -732,18 +1025,14 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "np.allclose(hess_vmat_tmp, hess_vmat)" - ] + "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "np.einsum(\"g, g, tgv -> tgv\", grid_weight, grid_frr, dR_rho1_0, optimize=True).shape" - ] + "source": [] }, { "cell_type": "code", @@ -764,29 +1053,14 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "hess_deriv2 = np.zeros((natm, natm, 3, 3))\n", - "for A in range(natm):\n", - " for B in range(natm):\n", - " sA, sB = mol_slice(A), mol_slice(B)\n", - " hess_deriv2[A, B] = np.einsum(\"tsuv, uv -> ts\", vmat[A, :, :, sB], D[sB]) * 2\n", - " hess_deriv2[A, B] += np.einsum(\"tsuv, uv -> ts\", ipip[:, :, sB, sA], D[sB, sA]) * 2\n", - " hess_deriv2[A, B] += np.einsum(\"stvu, uv -> ts\", ipip[:, :, sA, sB], D[sB, sA]) * 2" - ] + "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "hess_deriv2_pyscf = np.zeros((natm, natm, 3, 3))\n", - "for A in range(natm):\n", - " for B in range(natm):\n", - " sA, sB = mol_slice(A), mol_slice(B)\n", - " hess_deriv2_pyscf[A, B] = np.einsum(\"tsuv, uv -> ts\", vxc_deriv2_pyscf[A, :, :, sB], D[sB]) * 2\n", - "np.allclose(hess_deriv2, hess_deriv2_pyscf)" - ] + "source": [] }, { "cell_type": "code", @@ -821,12 +1095,7 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "np.allclose(\n", - " hess_noU_noGGA + hess_veff_diag + hess_deriv2,\n", - " partial_hess_pyscf\n", - ")" - ] + "source": [] }, { "cell_type": "code", @@ -1167,7 +1436,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.0" + "version": "3.6.6" } }, "nbformat": 4,