test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# '''\n",
    "# @breif: function for debugging\n",
    "# @param x: input\n",
    "# @return y\n",
    "# '''\n",
    "\n",
    "# def fun3(x):\n",
    "#     y = (1.0/3.0)*x**3 - 0.5*x**2 - x - 1.0\n",
    "#     return y\n",
    "\n",
    "# def func(x):\n",
    "#     y = math.cos(x)\n",
    "#     return y\n",
    "\n",
    "# '''\n",
    "# @breif: parabolic interpolation method\n",
    "# @param x1, y1, x2, y2, x3, y3: three points\n",
    "# @return p, q: x4 = x3 + p/q               \n",
    "# '''\n",
    "\n",
    "# # x3<x2, x1\n",
    "# #x = (x1^2*(y2-y3)+x2^2*(y3-y1)+x3^2*(y1-y2))/(2*(x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2)))\n",
    "# #which is invariant under the permutation of the three points\n",
    "# def parabolic_interpolation(x1, y1, x2, y2, x3, y3):\n",
    "#     p = (x2 - x3)**2 * (y3 - y1) + (x1 - x3)**2 * (y2 - y3)\n",
    "#     q = 2.0 * ((x2 - x3) * (y3 - y1) + (x1 - x3) * (y2 - y3))\n",
    "#     return p, q\n",
    "    \n",
    "# '''\n",
    "# @breif: brent method in one iteration, in one iteration, only one new point is evaluated(fu)\n",
    "\n",
    "# @param a, b: boundary of the interval\n",
    "# @param v, fv: the old value of w\n",
    "# @param w, fw: second least value\n",
    "# @param x, fx: the least value\n",
    "# @param fun_p: the function to be minimized\n",
    "# @param dx: the last step length\n",
    "# @param dxold: the last last step length\n",
    "\n",
    "# @return fv, fw, fx: updated values\n",
    "# @return a, b, v, w, x, dx, dxold: updated values\n",
    "# @return u: last evaluated point\n",
    "# '''\n",
    "\n",
    "# #* here I updated the values of fv, fw, fx to save the computation time, only fu is evaluated in one iteration\n",
    "# def brent(a, b, v, fv, w, fw, x, fx, fun_p, dx, dxold):\n",
    "#     CGOLD = 0.381966013\n",
    "#     xm = 0.5 * (a + b)\n",
    "\n",
    "#     #perform parabolic interpolation using v, w, x (x the minimum point)\n",
    "#     p, q = parabolic_interpolation(v, fv, w, fw, x, fx)\n",
    "\n",
    "#     #Boolean variables for the safety of parabolic interpolation\n",
    "#     #* qSafe: q is not zero\n",
    "#     qSafe = True if (q != 0) else False\n",
    "\n",
    "#     tempdx = p/q if qSafe else 0.0\n",
    "#     #* paraSafe:\n",
    "#     #* 1. q is not zero\n",
    "#     #* 2. the new point is within the interval\n",
    "#     #* 3. the new step length is less than half of the last last step length\n",
    "#     paraSafe = True if (qSafe and (a < x+tempdx) and (x+tempdx < b) and (abs(tempdx) < 0.5*abs(dxold))) else False\n",
    "\n",
    "#     #update the dx and store the last last step length for next iteration\n",
    "\n",
    "#     if paraSafe:\n",
    "#         dxold = dx\n",
    "#         dx = tempdx\n",
    "#     else:\n",
    "#         if x < xm:\n",
    "#             dxold = b - x\n",
    "#         else:\n",
    "#             dxold = a - x\n",
    "#         dx = CGOLD * dxold\n",
    "\n",
    "#     #! evaluate the new point\n",
    "#     u = x + dx\n",
    "#     fu = fun_p(u)\n",
    "\n",
    "#     if (fu < fx):\n",
    "#         #* update the boundary of interval\n",
    "#         a = a if (u < x) else x\n",
    "#         b = x if (u < x) else b\n",
    "#         #* update the values of v, w, x\n",
    "#         v = w\n",
    "#         fv = fw\n",
    "#         w = x\n",
    "#         fw = fx\n",
    "#         x = u\n",
    "#         fx = fu\n",
    "#     else:\n",
    "#         #* update the boundary of interval\n",
    "#         a = u if (u < x) else a\n",
    "#         b = b if (u < x) else u\n",
    "\n",
    "#         #* update the values of v, w, x\n",
    "#         if (fu < fw) or (w == x):\n",
    "#             v = w\n",
    "#             fv = fw\n",
    "#             w = u\n",
    "#             fw = fu\n",
    "#         elif (fu < fv) or (v == x) or (v == w):\n",
    "#             v = u\n",
    "#             fv = fu\n",
    "    \n",
    "#     return a, b, v, fv, w, fw, x, fx, dx, dxold\n",
    "\n",
    "# '''\n",
    "# @breif: brent method\n",
    "# @param a0, b0: initial boundary of the interval\n",
    "# @param xguess: initial guess of the minimum point\n",
    "# @param fun_p: the function to be minimized\n",
    "# @param tol: tolerance\n",
    "\n",
    "# @return xmin: the minimum point\n",
    "# @return icount: the number of iterations\n",
    "# '''\n",
    "# def brentMethod(a0, b0, xguess, fun_p, tol):\n",
    "#     Maxloop = 100\n",
    "\n",
    "#     icount = 0\n",
    "#     #! evaluate the function at the initial guess\n",
    "#     y0 = fun_p(xguess)\n",
    "#     a, b, v, fv, w, fw, x, fx, dx, dxold = brent(a0, b0, xguess, y0, xguess, y0, xguess, y0, fun_p, 0.0, 0.0)\n",
    "\n",
    "#     print(f\"x:\\t{x:.6f}\\ty:\\t{fx:.12f}\")\n",
    "\n",
    "#     for i in range(Maxloop):\n",
    "#         a, b, v, fv, w, fw, x, fx, dx, dxold = brent(a, b, v, fv, w, fw, x, fx, fun_p, dx, dxold)\n",
    "\n",
    "#         print(f\"x:\\t{x:.6f}\\ty:\\t{fx:.12f}\")\n",
    "\n",
    "#         icount += 1\n",
    "#         if (abs(b - a) < tol):\n",
    "#             xmin = 0.5 * (a + b)\n",
    "#             break\n",
    "            \n",
    "#     return xmin, icount\n",
    "        \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# '''\n",
    "# @breif: find the spin index in the gjf file\n",
    "# by locating the second blank line index, then add 1\n",
    "\n",
    "# @param filename: default filename is \"template.gjf\"\n",
    "\n",
    "# @retrun iopIndex: the index of the iop line (start from 1)\n",
    "# @return spinIndex: the index of the spin line (start from 1)\n",
    "# '''\n",
    "# def findSpinIndex(filename = \"template.gjf\"):\n",
    "\n",
    "#     with open(filename, \"r\") as gjf:\n",
    "#         icount = 0\n",
    "#         blankCount = 0\n",
    "        \n",
    "#         for line in gjf:\n",
    "#             icount += 1\n",
    "#             if line == \"\\n\":\n",
    "#                 blankCount += 1\n",
    "#             if blankCount == 2:\n",
    "#                 spinIndex = icount + 1\n",
    "#                 break\n",
    "        \n",
    "#         iopIndex = spinIndex - 4\n",
    "\n",
    "#     return iopIndex, spinIndex\n",
    "\n",
    "# '''\n",
    "# @brief: generate gjf with N+1, N, N-1 spins in different IOp(3/107, 3/108) values (int(wpara*1000):05 + 0:05 format)\n",
    "\n",
    "# @param Nspin: \"charge spin\" of N\n",
    "# @param Naspin: \"charge spin\" of N+1\n",
    "# @param Nmspin: \"charge spin\" of N-1\n",
    "# @param wpara: the parameter for the IOp(3/107, 3/108) value\n",
    "# @param filename: the template file name, default is \"template.gjf\"\n",
    "\n",
    "\n",
    "# @return: void, three gjf files: N.gjf, N+1.gjf, N-1.gjf\n",
    "# '''\n",
    "# #! here I use f\"{int(wpara*10000):05}\" to convert the float to string, by adding 5 zeros\n",
    "# #! it's not good, I suggest to use round(x,4) functon instead\n",
    "\n",
    "# def g16input(Nspin, Naspin, Nmspin, wpara, filename = \"template.gjf\"):\n",
    "\n",
    "#     iopIndex, spinIndex = findSpinIndex()\n",
    "\n",
    "#     with open(\"N.gjf\", \"w\") as Ngjf, open(\"N+1.gjf\", \"w\") as Nagjf, open(\"N-1.gjf\", \"w\") as Nmgjf:\n",
    "#         with open(filename, \"r\") as gjf:\n",
    "#             icount = 0\n",
    "#             for line in gjf:\n",
    "#                 icount += 1\n",
    "#                 if icount == iopIndex:\n",
    "#                     tempIop = line.strip()  #*remove the newline character\n",
    "#                     wpara = round(wpara, 4)    #*round the float to 4 decimal places\n",
    "#                     str_wpara = f\"{int(wpara*10000):05}\"   #*convert the float to string, by adding 5 zeros\n",
    "#                     str_iop = str_wpara + f\"{0:05}\"\n",
    "                    \n",
    "#                     Ngjf.write(f\"{tempIop} IOp(3/107={str_iop},3/108={str_iop})\\n\")\n",
    "#                     Nagjf.write(f\"{tempIop} IOp(3/107={str_iop},3/108={str_iop})\\n\")\n",
    "#                     Nmgjf.write(f\"{tempIop} IOp(3/107={str_iop},3/108={str_iop})\\n\")\n",
    "\n",
    "#                 elif icount == spinIndex:\n",
    "#                     Ngjf.write(f\"{Nspin}\\n\")\n",
    "#                     Nagjf.write(f\"{Naspin}\\n\")\n",
    "#                     Nmgjf.write(f\"{Nmspin}\\n\")\n",
    "#                 else:\n",
    "#                     Ngjf.write(line)\n",
    "#                     Nagjf.write(line)\n",
    "#                     Nmgjf.write(line)\n",
    "\n",
    "\n",
    "# '''\n",
    "# @brief: read the SCF energy and eigenvalues of HOMO orbital\n",
    "\n",
    "# @param: filename the log file\n",
    "\n",
    "# @return: SCFenergy, eHOMO\n",
    "# '''\n",
    "\n",
    "# #* NOTICE: sometimes there will be not one SCF iterations in a log file,\n",
    "# #* such as stable=opt, and we take the newest values\n",
    "# def g16read(filename):\n",
    "#     success = False\n",
    "\n",
    "#     with open(filename, \"r\") as log:\n",
    "#         for line in log:\n",
    "#             if \"SCF Done\" in line:\n",
    "#                 templine = line.split()\n",
    "#                 SCFenergy = float(templine[4])\n",
    "#             if \"Alpha virt. eigenvalues\" in line and \"Alpha  occ. eigenvalues\" in prev_line:\n",
    "#                 templine = prev_line.split()\n",
    "#                 eHOMO = float(templine[-1])\n",
    "#             #* check if succeed\n",
    "#             if \"Normal termination\" in line:\n",
    "#                 success = True\n",
    "#             prev_line = line    #! prev_line should be updated in the end of line\n",
    "\n",
    "#         if not success:\n",
    "#             raise Exception(\"Error: Gaussian doesn't finish correctly. (Does not finished normally.)\")\n",
    "    \n",
    "#     return SCFenergy, eHOMO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# '''\n",
    "# @brief: read the energy from the Gaussian output file, as a function y(x)\n",
    "\n",
    "# @param wpara: parameter of functional LC-wPBE, as x\n",
    "# @return: J2, as y\n",
    "\n",
    "# @global: Nspin, Naspin, Nmspin\n",
    "# '''\n",
    "\n",
    "# def calIP(wpara):\n",
    "#     global Nspin, Naspin, Nmspin\n",
    "\n",
    "#     g16input(Nspin, Naspin, Nmspin, wpara)\n",
    "#     #* run g16\n",
    "#     #! remove the module load command later\n",
    "#     exitcode1 = os.system(\"module load Gaussian16-B01;g16 < N.gjf > N.log\")\n",
    "#     exitcode2 = os.system(\"module load Gaussian16-B01;g16 < N+1.gjf > N+1.log\")\n",
    "#     exitcode3 = os.system(\"module load Gaussian16-B01;g16 < N-1.gjf > N-1.log\")\n",
    "\n",
    "#     #* check the exit code (non-zero means error)\n",
    "#     if exitcode1 or exitcode2 or exitcode3:\n",
    "#         print(\"Error in Gaussian16\")\n",
    "#         raise Exception(\"Error: Gaussian doesn't finish correctly. (exitcode does not equal to 0)\")\n",
    "\n",
    "#     #* store the output files\n",
    "#     tempstr = f\"{int(round(wpara,4)*10000):05}\" \n",
    "#     os.system(f\"mkdir {tempstr}\")\n",
    "#     os.system(f\"cp *.log {tempstr}\")\n",
    "\n",
    "#     #* extract the energy\n",
    "#     EN, eN = g16read(\"N.log\")\n",
    "#     ENa, eNa = g16read(\"N+1.log\")\n",
    "#     ENm, eNm = g16read(\"N-1.log\")\n",
    "\n",
    "#     JN = abs(eN + ENm - EN)\n",
    "#     JNa = abs(eNa + EN - ENa)\n",
    "#     J = JN + JNa\n",
    "#     J2 = JN**2 + JNa**2\n",
    "\n",
    "#     return J, J2\n",
    "\n",
    "# '''\n",
    "# @brief: point to the function calIP, only return J2\n",
    "# #! Here wpara -> round(wpara, 4)\n",
    "\n",
    "# @param wpara: parameter of functional LC-wPBE\n",
    "# '''\n",
    "# def funJ2(wpara):\n",
    "#     global wparaList, J2List, JList, Nspin, Naspin, Nmspin\n",
    "#     if round(wpara, 4) in wparaList:\n",
    "#         #print(f\"Already calculated.{round(wpara, 4)} {wparaList}\")\n",
    "#         return J2List[wparaList.index(wpara)]\n",
    "#     J, J2 = calIP(wpara)\n",
    "#     wparaList.append(round(wpara, 4))\n",
    "#     JList.append(J)\n",
    "#     J2List.append(J2)\n",
    "#     return J2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "rm: cannot remove '06000': No such file or directory\n",
      "rm: cannot remove '03899': No such file or directory\n",
      "rm: cannot remove '02601': No such file or directory\n",
      "rm: cannot remove '04688': No such file or directory\n",
      "rm: cannot remove '04464': No such file or directory\n",
      "rm: cannot remove '04413': No such file or directory\n",
      "rm: cannot remove '04392': No such file or directory\n",
      "rm: cannot remove '04390': No such file or directory\n",
      "rm: cannot remove '04203': No such file or directory\n",
      "rm: cannot remove '04319': No such file or directory\n",
      "rm: cannot remove '04363': No such file or directory\n",
      "rm: cannot remove '04380': No such file or directory\n",
      "rm: cannot remove '04388': No such file or directory\n"
     ]
    }
   ],
   "source": [
    "from fileIO import fileClean\n",
    "fileClean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "wparaList = []\n",
    "J2List = []\n",
    "JList = []\n",
    "\n",
    "chargeList = []\n",
    "chargeSpinList = []\n",
    "\n",
    "UDFT = []\n",
    "\n",
    "oldchkBool = False\n",
    "restartBool = False\n",
    "cleanBool = False\n",
    "#!NOTICE: N in line -> N+1 and N-1 also satisfy the condition\n",
    "#!NOTICE: so it should be \"N:\" in line\n",
    "\n",
    "#! global variables:\n",
    "# chargeList, chargeSpinList, xmin, xmax, xguess, tol, cleanBool, restartBool, oldchkBool, wparaList, J2List, JList\n",
    "\n",
    "readCharge = False\n",
    "tmplist2 = []\n",
    "with open(\"input\", \"r\") as input:\n",
    "    for line in input:\n",
    "        if \"restart\" in line:\n",
    "            restartBool = True\n",
    "        \n",
    "        if \"oldchk\" in line:\n",
    "            oldchkBool = True\n",
    "\n",
    "        if \"clean\" in line:\n",
    "            cleanBool = True    \n",
    "\n",
    "        if \"orbit\" in line:\n",
    "            templine = line.strip().split()\n",
    "            for i in range(1,len(templine)):\n",
    "                tcharge1, tcharge2 = templine[i].split(\",\")\n",
    "                #*there should be one extra file (N-1)\n",
    "                tmplist = list(range(int(tcharge1)-1, int(tcharge2)+1))\n",
    "                chargeList.extend(tmplist)\n",
    "            #*remove the duplicate\n",
    "            chargeList = list(set(chargeList))\n",
    "            chargeList.sort()\n",
    "\n",
    "        #* turn on the readCharge flag\n",
    "        if \"charge spin start\" in line:\n",
    "            readCharge = True\n",
    "            continue\n",
    "\n",
    "        if \"charge spin end\" in line:\n",
    "            readCharge = False\n",
    "\n",
    "        #* finish reading the charge spin\n",
    "\n",
    "        if readCharge:\n",
    "            tmplist1 = []\n",
    "            templine = line.strip().split()\n",
    "            #* templine: [charge in chargeList, charge, spin]\n",
    "            for i in range(3):\n",
    "                tmplist1.append(int(templine[i]))\n",
    "            #* append the last element (may be 'U' or just duplicate of the charge spin)\n",
    "            tmplist1.append(templine[-1])\n",
    "            tmplist2.append(tmplist1)\n",
    "\n",
    "            #* sort the charge spin list, by charge\n",
    "            tmplist2 = sorted(tmplist2, key = lambda x: x[0])\n",
    "\n",
    "        if \"xmin\" in line:\n",
    "            templine = line.strip().split()\n",
    "            xmin = float(templine[-1])\n",
    "        if \"xmax\" in line:\n",
    "            templine = line.strip().split()\n",
    "            xmax = float(templine[-1])\n",
    "        if \"xguess\" in line:\n",
    "            templine = line.strip().split()\n",
    "            xguess = float(templine[-1])\n",
    "        if \"tolerance\" in line:\n",
    "            templine = line.strip().split()\n",
    "            tol = float(templine[-1])\n",
    "    \n",
    "\n",
    "#* check if the charge spin is complete\n",
    "if len(tmplist2) != len(chargeList):\n",
    "    raise Exception(\"Error: charge spin is not complete.\")\n",
    "#* convert the charge spin list to string\n",
    "for i in range(len(tmplist2)):\n",
    "    tmpstr = f\"{tmplist2[i][1]} {tmplist2[i][2]}\"\n",
    "    if tmplist2[i][-1] == \"U\":\n",
    "        UDFT.append(True)\n",
    "    else:\n",
    "        UDFT.append(False)\n",
    "    chargeSpinList.append(tmpstr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"output.txt\", \"a\") as foutput:\n",
    "    foutput.write(\"chargeList:\\n\")\n",
    "    foutput.flush()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['#guess=read', 'stable=opt', 'ULC-wPBE']\n",
      "#guess=read stable=opt ULC-wPBE/6-311+G** opt freq\n"
     ]
    }
   ],
   "source": [
    "print(tempIop1)\n",
    "print(ttt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "beta = np.zeros((3,3,3))\n",
    "temp_beta = np.zeros((3,3,3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [],
   "source": [
    "#* input is XXX XXY XYY YYY XXZ XYZ YYZ XZZ YZZ ZZZ\n",
    "ccsd = \"\"\"\n",
    "XXX=       -150.20500\n",
    "XXY=          3.12910\n",
    "XYY=         66.98390\n",
    "YYY=        -99.80810\n",
    "XXZ=          3.31752\n",
    "XYZ=          1.76751\n",
    "YYZ=          0.25054\n",
    "XZZ=         66.74340\n",
    "YZZ=         98.96630\n",
    "ZZZ=          0.00000\n",
    "\"\"\"\n",
    "cam = \"\"\"\n",
    " XXX=       -127.52800\n",
    " XXY=         -0.00293\n",
    " XYY=         63.76490\n",
    " YYY=        -90.16700\n",
    " XXZ=          0.00000\n",
    " XYZ=          0.00000\n",
    " YYZ=          0.00000\n",
    " XZZ=         63.75850\n",
    " YZZ=         90.16990\n",
    " ZZZ=          0.00000\n",
    "\"\"\"\n",
    "mp2 = \"\"\"\n",
    " XXX=       -194.58000\n",
    " XXY=         -0.00788\n",
    " XYY=         97.28350\n",
    " YYY=       -137.56900\n",
    " XXZ=          0.00000\n",
    " XYZ=          0.00000\n",
    " YYZ=          0.00000\n",
    " XZZ=         97.27890\n",
    " YZZ=        137.58000\n",
    " ZZZ=          0.00000\n",
    "\"\"\"\n",
    "\n",
    "w0pbe = \"\"\"\n",
    " XXX=        -99.28240\n",
    " XXY=         -0.00671\n",
    " XYY=         49.64040\n",
    " YYY=        -70.19600\n",
    " XXZ=          0.00000\n",
    " XYZ=          0.00000\n",
    " YYZ=          0.00000\n",
    " XZZ=         49.65120\n",
    " YZZ=         70.21150\n",
    " ZZZ=          0.00000\n",
    "\"\"\"\n",
    "\n",
    "wpbe = \"\"\"\n",
    " XXX=        -103.928\n",
    " XXY=         -0.02928\n",
    " XYY=         51.9619\n",
    " YYY=        -73.4829\n",
    " XXZ=          0.00000\n",
    " XYZ=          0.00000\n",
    " YYZ=          0.00000\n",
    " XZZ=         51.9616\n",
    " YZZ=         73.4802\n",
    " ZZZ=          0.00000\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calAvBeta(functional):\n",
    "\n",
    "    lines = functional.split(\"\\n\")\n",
    "\n",
    "    beta_static = []\n",
    "\n",
    "    icount = 0\n",
    "    for line in lines:\n",
    "        icount += 1\n",
    "        #* skip the first and last line\n",
    "        if icount == 1 or icount == len(lines):\n",
    "            continue\n",
    "        beta_static.append(float(line.split(\"=\")[1]))\n",
    "\n",
    "    temp_beta[0,0,0] = beta_static[0]\n",
    "    temp_beta[0,0,1] = beta_static[1]\n",
    "    temp_beta[0,1,1] = beta_static[2]\n",
    "    temp_beta[1,1,1] = beta_static[3]\n",
    "    temp_beta[0,0,2] = beta_static[4]\n",
    "    temp_beta[0,1,2] = beta_static[5]\n",
    "    temp_beta[1,1,2] = beta_static[6]\n",
    "    temp_beta[0,2,2] = beta_static[7]\n",
    "    temp_beta[1,2,2] = beta_static[8]\n",
    "    temp_beta[2,2,2] = beta_static[9]\n",
    "\n",
    "    #* symmetrize the beta tensor\n",
    "    for i in range(3):\n",
    "        for j in range(3):\n",
    "            for k in range(3):\n",
    "                #* index of temp_beta are sorted\n",
    "                itmp, jtmp, ktmp = sorted([i,j,k])\n",
    "\n",
    "                beta[i,j,k] = temp_beta[itmp,jtmp,ktmp]\n",
    "                beta[i,k,j] = temp_beta[itmp,jtmp,ktmp]\n",
    "                beta[j,i,k] = temp_beta[itmp,jtmp,ktmp]\n",
    "                beta[j,k,i] = temp_beta[itmp,jtmp,ktmp]\n",
    "                beta[k,i,j] = temp_beta[itmp,jtmp,ktmp]\n",
    "                beta[k,j,i] = temp_beta[itmp,jtmp,ktmp]\n",
    "\n",
    "\n",
    "    sum1 = 0.0  # sum_i beta_iii**2\n",
    "    sum2 = 0.0  # sum_i!=j beta_iii*beta_ijj\n",
    "    sum3 = 0.0  # sum_i!=j beta_ijj**2\n",
    "    sum4 = 0.0  # sum_i!=j!=k beta_ijj*beta_ikk\n",
    "    sum5 = 0.0  # sum_i!=j!=k beta_ijk**2\n",
    "\n",
    "    for i in range(3):\n",
    "        for j in range(3):\n",
    "            for k in range(3):\n",
    "                if i == j == k:\n",
    "                    sum1 += beta[i,i,i]**2\n",
    "                elif j == k and i != j:\n",
    "                    sum2 += beta[i,i,i]*beta[i,j,j]\n",
    "                    sum3 += beta[i,j,j]**2\n",
    "                elif i != j != k:\n",
    "                    sum4 += beta[i,j,j]*beta[i,k,k]\n",
    "                    sum5 += beta[i,j,k]**2\n",
    "\n",
    "    ZZZ = 1.0/7.0 * sum1 + 6.0/35.0 * sum2 + 9.0/35.0 * sum3 + 3.0/35.0 * sum4 + 2.0/35.0 * sum5\n",
    "    ZXX = 1.0/35.0 * sum1 - 2.0/105.0 * sum2 + 11.0/105.0 * sum3 - 1.0/105.0 * sum4 + 4.0/105.0 * sum5\n",
    "\n",
    "    return ZZZ, ZXX"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ZZZ: 3576.475729172228\n",
      "ZXX: 4383.44912601458\n",
      "DR 0.8159044684576847\n"
     ]
    }
   ],
   "source": [
    "ZZZ, ZXX = calAvBeta(ccsd)\n",
    "\n",
    "print(f\"ZZZ: {ZZZ}\")\n",
    "print(f\"ZXX: {ZXX}\")\n",
    "print(f\"DR {ZZZ/ZXX}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ZZZ: 3020.058478364339\n",
      "ZXX: 3639.5649401114997\n",
      "DR 0.8297855727426087\n"
     ]
    }
   ],
   "source": [
    "ZZZ, ZXX = calAvBeta(cam)\n",
    "\n",
    "print(f\"ZZZ: {ZZZ}\")\n",
    "print(f\"ZXX: {ZXX}\")\n",
    "print(f\"DR {ZZZ/ZXX}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ZZZ: 1830.956053262431\n",
      "ZXX: 2206.370943826681\n",
      "DR 0.8298496036604168\n"
     ]
    }
   ],
   "source": [
    "ZZZ, ZXX = calAvBeta(w0pbe)\n",
    "\n",
    "print(f\"ZZZ: {ZZZ}\")\n",
    "print(f\"ZXX: {ZXX}\")\n",
    "print(f\"DR {ZZZ/ZXX}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ZZZ: 2005.7940266901253\n",
      "ZXX: 2417.059156079924\n",
      "DR 0.8298489599001774\n"
     ]
    }
   ],
   "source": [
    "ZZZ, ZXX = calAvBeta(wpbe)\n",
    "\n",
    "print(f\"ZZZ: {ZZZ}\")\n",
    "print(f\"ZXX: {ZXX}\")\n",
    "print(f\"DR {ZZZ/ZXX}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ZZZ: 7030.472649487953\n",
      "ZXX: 8472.555796244344\n",
      "DR 0.8297936087484219\n"
     ]
    }
   ],
   "source": [
    "ZZZ, ZXX = calAvBeta(mp2)\n",
    "\n",
    "print(f\"ZZZ: {ZZZ}\")\n",
    "print(f\"ZXX: {ZXX}\")\n",
    "print(f\"DR {ZZZ/ZXX}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}