diff --git a/notebooks/fanci_tutorial.ipynb b/notebooks/fanci_tutorial.ipynb index 26a1d52..54da886 100644 --- a/notebooks/fanci_tutorial.ipynb +++ b/notebooks/fanci_tutorial.ipynb @@ -51,7 +51,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 68, "metadata": {}, "outputs": [], "source": [ @@ -63,9 +63,12 @@ "from pyci.test import datafile\n", "\n", "# System information\n", - "filename = datafile(\"lih_sto6g.fcidump\")\n", + "# filename = datafile(\"lih_sto6g.fcidump\")\n", + "filename = datafile(\"/workspaces/PyCI/pyci/test/data/h4_2_0.0_STO-6G.fcidump\")\n", + "\n", "ham = pyci.hamiltonian(filename)\n", - "e_dict = {}" + "e_dict = {}\n", + "f_dict = {}" ] }, { @@ -79,7 +82,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 59, "metadata": {}, "outputs": [], "source": [ @@ -101,7 +104,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 60, "metadata": {}, "outputs": [], "source": [ @@ -124,13 +127,28 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 61, "metadata": {}, "outputs": [], "source": [ "# Optimize wavefunction\n", "ap1rog_results = ap1rog.optimize(ap1_params, use_jac=True)\n", - "e_dict[\"AP1roG\"] = ap1rog_results.x[-1]" + "e_dict[\"AP1roG\"] = (ap1rog_results.x[-1])\n", + "\n", + "# f_dict['AP1roG'] = (np.std(ap1rog_results.fun)) #Calculating the std of the residuals \n", + "\n", + "#Calculation of root mean square deviation (RMSD) \n", + "\"\"\" \n", + "RMSD = \\sqrt ( \\sum_i^n (residuals)^2 / (n - 2)) \n", + "n = size of sample \n", + "resilduals = Vector of residuals at the solution.\n", + "\"\"\"\n", + "f_dict['AP1roG'] = (np.sqrt(np.sum((ap1rog_results.fun)**2)/(len(ap1rog_results.fun)-2)))\n", + "\n", + "# print(f_dict['AP1roG'])\n", + "\n", + "# print(ap1rog_results) \n", + "\n" ] }, { @@ -146,7 +164,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 62, "metadata": {}, "outputs": [], "source": [ @@ -160,7 +178,22 @@ "\n", "# Optimize wavefunction\n", "pccds_results = pccds.optimize(pccds_params, use_jac=False)\n", - "e_dict[\"pCCD+S\"] = pccds_results.x[-1]" + "e_dict[\"pCCD+S\"] = pccds_results.x[-1]\n", + "# f_dict[\"pCCD+S\"] = (np.std(ap1rog_results.fun))\n", + "\n", + "#Calculation of root mean square deviation (RMSD) \n", + "\"\"\" \n", + "RMSD = \\sqrt ( \\sum_i^n (residuals)^2 / (n - 2)) \n", + "n = size of sample \n", + "resilduals = Vector of residuals at the solution.\n", + "\"\"\"\n", + "\n", + "f_dict['pCCD+S'] = (np.sqrt(np.sum((pccds_results.fun)**2)/(len(pccds_results.fun)-2)))\n", + "\n", + "\n", + "# print(pccds_results)\n", + "# print(f_dict['pCCD+S'])\n", + "\n" ] }, { @@ -174,7 +207,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 63, "metadata": {}, "outputs": [ { @@ -182,9 +215,9 @@ "output_type": "stream", "text": [ " METHOD, ENERGY\n", - " HF, -8.947289175e+00\n", - " AP1roG, -8.963531034e+00\n", - " pCCD+S, -8.963613544e+00\n" + " HF, -3.875079433e+00\n", + " AP1roG, -3.955402403e+00\n", + " pCCD+S, -3.816787845e+00\n" ] } ], @@ -195,6 +228,54 @@ " print(f\"{name:>8s}, {energy:>16.9e}\")" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Calculate the standard deviation " + ] + }, + { + "cell_type": "code", + "execution_count": 66, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2.4683556322835916e-14\n", + "1.9667674453325034e-09\n" + ] + } + ], + "source": [ + "print(f_dict[\"AP1roG\"])\n", + "print(f_dict[\"pCCD+S\"])\n" + ] + }, + { + "cell_type": "code", + "execution_count": 65, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " METHOD, RMSD\n", + " AP1roG, 2.468355632e-14\n", + " pCCD+S, 1.966767445e-09\n" + ] + } + ], + "source": [ + "# Print RMSD from various methods\n", + "print(f\"{'METHOD':>8s}, {'RMSD':>16s}\")\n", + "for name, fun in f_dict.items():\n", + " print(f\"{name:>8s}, {fun:>16.9e}\")" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -242,7 +323,7 @@ "output_type": "stream", "text": [ "Number of electrons = 4.1\n", - "Number of pairs = 6.1\n" + "Number of pairs = 6.2\n" ] } ], @@ -269,7 +350,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.6" + "version": "3.10.13" } }, "nbformat": 4, diff --git a/pyci/rdm/algorithms.py b/pyci/rdm/algorithms.py index e69de29..15b261c 100644 --- a/pyci/rdm/algorithms.py +++ b/pyci/rdm/algorithms.py @@ -0,0 +1,65 @@ +import numpy as np +from abc import ABCMeta, abstractmethod +class Projection(metaclass=ABCMeta): + + @abstractmethod + def __init__(self, initial_guess,constraints): + self.initial_guess = initial_guess + self.constraints = constraints + + @abstractmethod + def optimize(self): + pass + +class Dykstra(metaclass=ABCMeta): + def __init__(self, initial_guess:float, constraints:list, alpha:int =1, max_iterations:int =100, eps:int =1e-6): + """ + Dykstra's algorithm for projection onto convex sets. + + Parameters + ---------- + initial_guess : float + Initializing with guessed density matrix \Gamma_0. + constraints : list + A list of projection operators onto all $J$ of the constraints as a single set, $\{\mathcal{P}_j\}_{j=0}^{J-1}$. + alpha : int, optional + Tunning parameter, by default 1 + max_iterations : int, optional + Number of maximum iterations, by default 100 + eps : int, optional + Tolerance, by default 1e-6 + """ + super().__init__(initial_guess, constraints) + self.alpha = alpha + self.max_iterations = max_iterations + self.eps = eps + def optimize(self): + X = [np.zeros(self.initial_guess) for i in range(len(self.constraints))] + D = np.copy(self.initial_guess) + norm = [] + for i in range (self.max_iterations): + for j, projection in enumerate(self.constraints): + C = D - X[j] + D = projection(C) + X[j] = D - C + norm.append(np.linalg.norm(D - self.initial_guess)) + is_stop = self.alpha * abs(norm[i] - norm[i - 1]) + (1 - self.alpha) * norm[i] < self.eps + if is_stop: + break +class Neumann(metaclass=ABCMeta): + def __init__(self, initial_guess, constraints, alpha, max_iterations=100, eps=1e-6): + super().__init__(initial_guess, constraints) + self.alpha = alpha + self.max_iterations = max_iterations + self.eps = eps + def optimize(self): + D = np.copy(self.initial_guess) + norm = [] + for i in range (self.max_iterations): + for projection in enumerate(self.constraints): + D = projection(D) + norm.append(np.linalg.norm(D - self.initial_guess)) + is_stop = self.alpha * abs(norm[i] - norm[i - 1]) + (1 - self.alpha) * norm[i] < self.eps + if is_stop: + break + diff --git a/pyci/test/data/h4_2_0.0_STO-6G.fcidump b/pyci/test/data/h4_2_0.0_STO-6G.fcidump new file mode 100644 index 0000000..c2f8a60 --- /dev/null +++ b/pyci/test/data/h4_2_0.0_STO-6G.fcidump @@ -0,0 +1,115 @@ + &FCI NORB= 4,NELEC= 4,MS2=0, + ORBSYM=1,1,1,1, + ISYM=1, + &END + 0.5068898550551018 1 1 1 1 + 0.003794527559320846 1 1 2 1 + 0.4098635263371225 1 1 2 2 + 1.025916213093225e-09 1 1 3 1 + -8.035957177465747e-10 1 1 3 2 + 0.445100515768005 1 1 3 3 + 0.09719871204645167 1 1 4 1 + 0.05096181960669316 1 1 4 2 + -2.511664404236669e-09 1 1 4 3 + 0.5183452230109675 1 1 4 4 + 0.003794527559320839 2 1 1 1 + 0.08813159971644449 2 1 2 1 + 0.07782243723669213 2 1 2 2 + -5.653522719417658e-10 2 1 3 1 + -4.283259379123106e-09 2 1 3 2 + -0.08764576192895299 2 1 3 3 + 0.02826446097635334 2 1 4 1 + -0.08001497578900814 2 1 4 2 + 7.330390461302017e-11 2 1 4 3 + 0.05589849747028398 2 1 4 4 + 0.4098635263371225 2 2 1 1 + 0.07782243723669217 2 2 2 1 + 0.521293056868916 2 2 2 2 + -5.046337625547181e-09 2 2 3 1 + 1.404773519730895e-11 2 2 3 2 + 0.3561543663891216 2 2 3 3 + -0.001427950323337701 2 2 4 1 + -0.09114783741363394 2 2 4 2 + 3.876907371525018e-09 2 2 4 3 + 0.4713483506720944 2 2 4 4 + 1.025916174929309e-09 3 1 1 1 + -5.653522754112128e-10 3 1 2 1 + -5.046337653302757e-09 3 1 2 2 + 0.1293028829330219 3 1 3 1 + -0.1072466481774506 3 1 3 2 + 4.05264671821115e-09 3 1 3 3 + -1.581443448417374e-09 3 1 4 1 + 7.217758335453794e-11 3 1 4 2 + -0.0846859118356684 3 1 4 3 + -1.893284105880522e-09 3 1 4 4 + -8.035956726437643e-10 3 2 1 1 + -4.283259382592552e-09 3 2 2 1 + 1.404779764735409e-11 3 2 2 2 + -0.1072466481774506 3 2 3 1 + 0.09458993631478847 3 2 3 2 + 1.40467162901281e-09 3 2 3 3 + 7.391295561709477e-11 3 2 4 1 + 4.802033475720391e-09 3 2 4 2 + 0.07788008022122561 3 2 4 3 + -1.244524983068374e-09 3 2 4 4 + 0.445100515768005 3 3 1 1 + -0.08764576192895296 3 3 2 1 + 0.3561543663891216 3 3 2 2 + 4.052646773722302e-09 3 3 3 1 + 1.404671601257235e-09 3 3 3 2 + 0.5359374140556439 3 3 3 3 + -0.006072330211121087 3 3 4 1 + 0.1003613149338521 3 3 4 2 + -3.081953497963141e-09 3 3 4 3 + 0.4053549128731515 3 3 4 4 + 0.09719871204645164 4 1 1 1 + 0.02826446097635335 4 1 2 1 + -0.001427950323337657 4 1 2 2 + -1.581443447550013e-09 4 1 3 1 + 7.391295756865868e-11 4 1 3 2 + -0.006072330211121068 4 1 3 3 + 0.1186180412469477 4 1 4 1 + 0.03612665777791915 4 1 4 2 + -1.119209735306326e-09 4 1 4 3 + 0.122452034039052 4 1 4 4 + 0.05096181960669319 4 2 1 1 + -0.08001497578900814 4 2 2 1 + -0.09114783741363391 4 2 2 2 + 7.217762845734832e-11 4 2 3 1 + 4.802033409800899e-09 4 2 3 2 + 0.1003613149338522 4 2 3 3 + 0.03612665777791913 4 2 4 1 + 0.1272246486516095 4 2 4 2 + -1.005679005350313e-09 4 2 4 3 + 0.004098379798186072 4 2 4 4 + -2.511664473625608e-09 4 3 1 1 + 7.330390461302017e-11 4 3 2 1 + 3.876907420097275e-09 4 3 2 2 + -0.08468591183566843 4 3 3 1 + 0.07788008022122564 4 3 3 2 + -3.081953497963141e-09 4 3 3 3 + -1.119209738775773e-09 4 3 4 1 + -1.005678960247502e-09 4 3 4 2 + 0.06715293748280152 4 3 4 3 + -8.629244367674982e-10 4 3 4 4 + 0.5183452230109674 4 4 1 1 + 0.05589849747028395 4 4 2 1 + 0.4713483506720945 4 4 2 2 + -1.89328408506384e-09 4 4 3 1 + -1.244525038579525e-09 4 4 3 2 + 0.4053549128731514 4 4 3 3 + 0.1224520340390518 4 4 4 1 + 0.004098379798186058 4 4 4 2 + -8.629244610536269e-10 4 4 4 3 + 0.5934456368866865 4 4 4 4 + -1.80985796887928 1 1 0 0 + -0.08161693543269405 2 1 0 0 + -1.373368656657632 2 2 0 0 + 4.184586471194904e-09 3 1 0 0 + 5.930225890093923e-10 3 2 0 0 + -1.366056716783876 3 3 0 0 + -0.174357790674463 4 1 0 0 + 0.01748866944921174 4 2 0 0 + 5.075449520787198e-10 4 3 0 0 + -1.115558063357834 4 4 0 0 + 2.168608539661471 0 0 0 0