cc_ade-implementation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0147315f",
   "metadata": {},
   "source": [
    "### Algebraic Evaluation of Coupled-Cluster Diagrams - Application"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96bdf3a6",
   "metadata": {},
   "source": [
    "The code to evaluate the coupled-cluster diagrams starts with the exponential expansion of a sum of amplitude operators, for example, $e^{T_1 + T_2}$. To accomplish this we use Sympy. The program we are using is called __cc-ade__ and we start with, to get eg the doubles amplitude at single & doubles level"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9e4521b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('../cc-ade')\n",
    "from cc_ade import expansion\n",
    "\n",
    "SD = expansion('SD', 'D')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db997fbd",
   "metadata": {},
   "source": [
    "The first thing we can get back is a Latex representation of the expansion from __.do_latex()__ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "377f131f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SD doubles =  $\\frac{1}{24}T_1^{4} + \\frac{1}{6}T_1^{3} + \\frac{1}{2}T_1^{2}T_2 + \\frac{1}{2}T_1^{2} + T_1T_2 + T_1 + \\frac{1}{2}T_2^{2} + T_2 + 1 $\n"
     ]
    }
   ],
   "source": [
    "latex_SD = SD.do_latex()\n",
    "print('SD doubles = ', latex_SD)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "135c4816",
   "metadata": {},
   "source": [
    "SD doubles =  $\\frac{1}{24}T_1^{4} + \\frac{1}{6}T_1^{3} + \\frac{1}{2}T_1^{2}T_2 + \\frac{1}{2}T_1^{2} + T_1T_2 + T_1 + \\frac{1}{2}T_2^{2} + T_2 + 1 $\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "602ae66c",
   "metadata": {},
   "source": [
    "We can get information about the terms from __.do_term(id)__. Each term is identified by a string of four characters where the entry in the ith position indicates the power of the ith amplitude in the term (0 indicates amplitude does not appear in term), so $T_1^2 T_3$ would be identified with '2010'. We are only coding for upto $T_4$ hence the four characters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "4155521c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Level is [ SD ]  Order is [ 4 ]\n",
      "Term id  [ 2100 ]   Multiplier is [  1 / 2  ]\n",
      "Latex string is [  \\frac{1}{2} T_1^{2} T_2   ]\n",
      "Vertex string is [  +- | +- | ++--  ]\n",
      "Interaction number is [  - 2  ]\n",
      "\n",
      "Contractor               Contractions\n",
      "hhpp               ['+|+|--', '+|+-|-', '+-|-|+', '+|-|+-', '-|-|++']\n"
     ]
    }
   ],
   "source": [
    "SD.do_term('2100')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26281a9e",
   "metadata": {},
   "source": [
    "In our approach each amplitude is represented by a sign-string (s-string), for example, $T_3$ is __+++---__. A term ,which is product of amplitudes, is written as __+-|++--__ with a __|__ symbol seperating the amplitudes. We see above above for '2100' the term s-string is __+-|+-|++--__. The _interaction number_ of the term is the target excitation level (for doubles this is 2) minus the actual excitation level of the term (number of __+__ signs, for '2100' this is 4). _Contractor_ are the Hamiltonian elements that can contract with the term, that is those with an _interaction number_ equal to the _interaction number_ of the term. The _contractions_ are all the valid ways in which the Hamiltonian can contract with the term. It is possible to get the Hamiltonian definitions from __.vertex__ and __.interaction__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "fa7f365b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2-body Hamiltonian hhpp  ['+', '+', '-', '-'] [-2]\n"
     ]
    }
   ],
   "source": [
    "print('2-body Hamiltonian hhpp ', SD.vertex['hhpp'], SD.interaction['hhpp'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "818ca875",
   "metadata": {},
   "source": [
    "There is a __.get(type, source, id='')__ method which can be used as __.get('t', term_id)__ to get the internal representation of a term. This returns a list as \\[term_id, factor, amplitude powers, s-strings, interaction number \\]. It can also be used as __.get('i', interaction_number)__ which returns all Hamiltonian elements with the requested _interaction number_ as a list. If the list obtained by get('i') is passed to __.get('v', hamiltonian_list)__ then a list of the Hamiltonian s-strings is returned."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f49ce868",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['2100', [1, 2], [2, 1, 0, 0], [['+', '-'], ['+', '-'], ['+', '+', '-', '-']], -2]\n",
      "['hp', 'phpp', 'hhhp']\n",
      "[['+', '-'], ['+', '+', '-'], ['+', '-', '-']]\n"
     ]
    }
   ],
   "source": [
    "print(SD.get('t', '2100'))\n",
    "print(SD.get('i', -1))\n",
    "print(SD.get('v', SD.get('i', -1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "065d1015",
   "metadata": {},
   "source": [
    "For each term in the expansion, each possible Hamiltonian contraction corresponds to a diagram. The overall sign and multiplying factor of the diagram can determined from the s-string of the Hamiltonian operator. To determine the the permutation of external lines the diagram must be labelled. The labelling algorithms are coded in the _lablib_ module. The labelling of diagrams is not unique so _lablib_ is just an example of possible schemes. The _lablib_ module is not optimized in any way, for example Hamiltonian elements +-- and -++ are just mirrors of each other and could be parameterised in a single routine. The routines not computationally intensive so this has not been done incidentally allowing for clearer code. Each diagram is stored as a dictionary in a list in the class variable __.ade__. Although this demonstration program is coded for S, D, T and Q the _lablib_ routines are designed to be applicable for any order of cluster expansion. They have been extensively tested upto Q level (see the section on validation)\n",
    "\n",
    "__.get('r', reference)__ will return the __.ade__ list item matching the Shavitt & Bartlett reference supplied. __.get('c', s-string)__ will return the __.ade__ list item matching the s-string supplied. The s-string is of the form Hamiltonian contraction s-string __:__ term amplitudes s-string, so for the '2100' example in \\[3\\] one Hamitonian contraction s-string could be '+|+|--' and the term s-string is '+-|+-|++--' so we would pass '+|+|--:+-|+-|++--'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "7a9c43a1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.25 * np.einsum('klcd,ci,dj,ak,bl->abij', gs[o,o,v,v], ts, ts, ts, ts)\n",
      "\n",
      "{'id': '2100(hhpp)', 'labels': ('klcd', 'ci|dj|abkl'), 'sign': (1, [4, 2]), 'factor': 2, 'sigma': 'dklc', 'permutation': '\\\\hat{P}(ij) ', 'latex': '\\\\frac{1}{4}\\\\hat{P}(ij)  \\\\displaystyle \\\\sum_{dklc} \\\\langle kl\\\\Vertcd \\\\rangle t^{c}_{i} t^{d}_{j} t^{ab}_{kl}', 'code': \"0.25 * np.einsum('klcd,ci,dj,abkl->abij', gs[o,o,v,v], ts, ts, td)\", 'reference': 'D_{7a}', 'contraction': '+|+|--'}\n"
     ]
    }
   ],
   "source": [
    "print(SD.get('r', 'D_{9}')['code'])\n",
    "print()\n",
    "print(SD.get('c', '+|+|--:+-|+-|++--'))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f099f4d3",
   "metadata": {},
   "source": [
    "The diagrams are computed in a list of dictionaries in __.ade__. For each diagram there is a list entry with keys \n",
    "+ _id_ - '4000(hhpp)' - the term id and the Hamiltonian contraction\n",
    "+ _labels_ -  ('klcd', 'ci|dj|ak|bl') - tuple (operator label, amplitude labels)\n",
    "+ _sign_ - (1, \\[4, 2\\]) - tuple (sign, \\[hole lines, loops\\])\n",
    "+ _factor_ - 2 - number of $\\frac{1}{2}$ multipliers\n",
    "+ _sigma_ - lkdc - summation over these indices\n",
    "+ _permutation_ - $\\hat{P}(ab) \\hat{P}(ij)$ - latex representation of permutations\n",
    "+ _latex_ - $\\frac{1}{4}\\hat{P}(ab) \\hat{P}(ij)  \\displaystyle \\sum_{lkdc} \\langle kl\\Vert cd \\rangle t^{c}_{i} t^{d}_{j} t^{a}_{k} t^{b}_{l}$ - latex representation of diagram algebra\n",
    "+ _code_ - \"0.25 * np.einsum('klcd,ci,dj,ak,bl->ijab', gs[o,o,v,v], ts, ts, ts, ts)\" - Python einsum expression.\n",
    "+ _reference_ - $D_{9}$ - Shavitt & Bartlett reference for diagram\n",
    "+ _contraction_ - '+|+|-|-' - s-string of Hamiltonian\n",
    "\n",
    "There is no obvious order for the diagram so they are listed in arbitrary order. There is a __.order_by_reference()__ method which will order the __.ade()__ list by Shavitt and Bartlett reference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b36381e9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'id': '0000(pphh)',\n",
       " 'labels': ('abij', ''),\n",
       " 'sign': (1, [2, 2]),\n",
       " 'factor': 0,\n",
       " 'sigma': '',\n",
       " 'permutation': '',\n",
       " 'latex': ' \\\\langle ab\\\\Vertij \\\\rangle',\n",
       " 'code': \"np.einsum('abij->abij', gs[v,v,o,o])\",\n",
       " 'reference': 'D_{1}',\n",
       " 'contraction': '0'}"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SD.order_by_reference()\n",
    "SD.ade[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68ba9f02",
   "metadata": {},
   "source": [
    "### Validation\n",
    "In order to verify the diagrams computed in _ade_ there is a _validation_ suite. This consists of a _reference.txt_ file which contains the Python einsum expressions for the diagrams as determined from Shavitt & Bartlett. The validation process is initiated with a call to __validate(ade_object, test_set)__. The _test\\_set_ is either \n",
    "+ 'c' - singles, doubles, triples and quadruples for single carbon atom in 3-21g basis. __unconverged amplitudes for testing only!__\n",
    "+ 'hf'- hydrogen fluoride in singles, doubles and triples in 6-31g basis (r=1.6 bohr) \n",
    "\n",
    "The _validate_ module has a __.compare(type)__ method. _Type_ can be _'all'_ in which case every diagram in the _.ade_ list will be compared with the Shavitt & Bartlett definition. If _type_ is a Shavitt & Bartlett reference then that diagram will be compared. If _type_ is a Hamiltonian id then all contractions of that Hamiltonian will be compared. The output is, for example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f9a2573d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "D_{3a}    ++|--        0.018753199846591   0.018753199846591  True True\n",
      "D_{3b}    +-|+-        0.012007054940035   0.012007054940035  True True\n",
      "D_{3c}    ++-|-        0.030511793264890   0.030511793264890  True True\n",
      "D_{3d}    +|+--        0.015795670617637   0.015795670617637  True True\n",
      "D_{7a}    +|+|--       0.000000811094776   0.000000811094776  True True\n",
      "D_{7b}    -|-|++       0.000001579422770   0.000001579422770  True True\n",
      "D_{7c}    +|-|+-       0.000000609290951   0.000000609290951  True True\n",
      "D_{7d}    +|+-|-       0.000001370808451   0.000001370808451  True True\n",
      "D_{7e}    +-|-|+       0.000000339310896   0.000000339310896  True True\n",
      "D_{9}     +|+|-|-      0.000000000065099   0.000000000065099  True True\n"
     ]
    }
   ],
   "source": [
    "from validation.validate import validate\n",
    "test = validate(SD, 'c')\n",
    "test.compare('hhpp')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95aa72cc",
   "metadata": {},
   "source": [
    "The output shows the Shavitt and Bartlett reference, the Hamiltonian contraction, np.linalg.norm of Shavitt & Bartlett and _ade_ computations, booleans are comparisons of norms followed by np.allclose on Shavitt & Bartlett and _ade_ computed amplitudes. For triples and quadruples these comparisons can take some time to run - quadruples are an 8 deep loop! The results of the carbon and hydrogen fluoride validations are given in the files 'c-3-21g-sdtq.txt' and 'hf-6-31g-sdt.txt'.\n",
    "\n",
    "The _validate_ module contains code to do the permutations where you can pass eg a string like 'ijklabcd->a/b/cd' to a subroutine to obtain the requested permutation. Although the idea of this exercise is just to generate the algebraic form of the diagrams the _validate_ module provides an __.evaluate(reference, data)__ method which will evaluate a Shavitt & Bartlett type diagram reference given spin-orbital data. As an example we could calculate the CCD energy correction for HF using the saved .npz file in the _validation_ directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "c2d1ae36",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5048370301     0.5048370301   \n",
      "0.5350832715     0.0302462414   \n",
      "0.5749348951     0.0398516236   \n",
      "0.5918903905     0.0169554954   \n",
      "0.6005635406     0.0086731501   \n",
      "0.6047685304     0.0042049898   \n",
      "0.6068430272     0.0020744968   \n",
      "0.6078606123     0.0010175852   \n",
      "0.6083614319     0.0005008196   \n",
      "0.6086077304     0.0002462985   \n",
      "0.6087289979     0.0001212675   \n",
      "0.6087887081     0.0000597102   \n",
      "0.6088181245     0.0000294164   \n",
      "0.6088326194     0.0000144948   \n",
      "0.6088397639     0.0000071445   \n",
      "0.6088432861     0.0000035222   \n",
      "0.6088450229     0.0000017368   \n",
      "HF CCD energy correction is -0.169138  Hartree\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from validation.validate import validate\n",
    "\n",
    "#load saved data \n",
    "data = np.load('../cc-ade/validation/hf.npz')\n",
    "nvir, nocc = data['ts'].shape                               #get the occupied and virtual orbital numbers\n",
    "o, v, n = slice(None, nocc), slice(nocc, None), np.newaxis  #make slices for occupied and virtual orbitals\n",
    "\n",
    "#user_data are variable and values dictionary passed to eval\n",
    "user_data = {}\n",
    "user_data['gs'], user_data['fs'] = data['g'], data['f']     #spin eri and fock\n",
    "user_data['td'] = np.zeros_like(data['td'])                 #only need doubles amplitudes, initialise to zero\n",
    "user_data['o'], user_data['v'] = o, v                       #orbital slices\n",
    "\n",
    "#fock denominator\n",
    "eps = np.diag(user_data['fs'])\n",
    "d   = -eps[v,n,n,n] - eps[n,v,n,n] + eps[n,n,o,n] + eps[n,n,n,o]\n",
    "\n",
    "#run parameters\n",
    "max_cycle, tol, norm, converged = 50, 1e-6, 0, False\n",
    "\n",
    "#get diagrams for doubles at doubles level\n",
    "D = expansion('D', 'D')\n",
    "V = validate(D, 'hf')                                        #need instance of validate for evaluate routine\n",
    "\n",
    "#iteration\n",
    "for cycle in range(0, max_cycle):\n",
    "    \n",
    "    #compute double amplitudes for this cycle\n",
    "    t = np.zeros_like(user_data['td'])\n",
    "    for diagram in D.ade:\n",
    "        t += V.evaluate(diagram['reference'], user_data)     #validate gets einsum code from ade list and handles\n",
    "                                                             #permuations automatically\n",
    "    #update user_data\n",
    "    user_data['td'] = t * np.reciprocal(d) + user_data['td'] #generated amplitudes seed for next evolution\n",
    "\n",
    "    #test convergence\n",
    "    cycle_norm = np.linalg.norm(user_data['td'])\n",
    "    if abs(cycle_norm - norm) < tol:\n",
    "        converged = True\n",
    "        break\n",
    "    \n",
    "    print('{:<15.10f}  {:<15.10f}'.format(cycle_norm, abs(cycle_norm - norm)))\n",
    "    norm = cycle_norm\n",
    "    \n",
    "if converged:\n",
    "    #get the energy\n",
    "    E = expansion('D', 'E')\n",
    "    V = validate(E, 'hf')                                     #new instance for different ade list\n",
    "    energy_correction = V.evaluate('E_{1}', user_data)\n",
    "    print('HF CCD energy correction is {:<10.6f} Hartree'.format(energy_correction))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c421756a",
   "metadata": {},
   "source": [
    "Which agrees with reference sources. We are using eri and Fock values from the saved hydrogen fluoride data together with _cc-ade_ generated doubles amplitudes which are passed to an eval expression for evaluation. The eval expression evaluated is obtained from the _code_ and _permutation_ keys of the diagram list generated by _expansion_. This is not an efficient way of getting the amplitudes but has the advantage that the contribution of each diagram could be easily monitored. Additionally for higher orders where evaluation of individual diagrams can take time, it would be possible to detect diagrams which are only contributing values below some threshold and remove them from the active list. This 'selected' approach may be worthwhile as the logic tests needed to implement it are negligible compared to evaluation of quadruple expansion terms involving loops to a depth of eight. \n",
    "\n",
    "Also provided is a way of using the __expansion__ class instantisation to deal with a single contraction by invoking the class as __expansion(contraction=s_strings)__. To evaluate, for example, the triples addition to the singles (diagram $S_7$) we would use _expansion(contraction='++--:+++---')_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "5ac22015",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'id': '0010(hhpp)', 'labels': ('klcd', 'adcikl'), 'sign': (-1, [3, 2]), 'factor': 2, 'sigma': 'dklc', 'permutation': '', 'latex': '-\\\\frac{1}{4} \\\\displaystyle \\\\sum_{dklc} \\\\langle kl\\\\Vertcd \\\\rangle t^{adc}_{ikl}', 'code': \"-0.25 * np.einsum('klcd,adcikl->ai', gs[o,o,v,v], tt)\", 'reference': 'S_{7}', 'contraction': '++--'}\n",
      "\n",
      "S_{7}     ++--         0.001552802417424   0.001552802417424  True True\n"
     ]
    }
   ],
   "source": [
    "s = expansion(contraction='++--:+++---')\n",
    "print(s.get('r', 'S_{7}'))\n",
    "V = validate(s, 'c')\n",
    "print()\n",
    "V.compare('S_{7}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4291d262",
   "metadata": {},
   "source": [
    "Using the class invocation __expansion(order=n, series='latex series')__ it is possible to explore custom series or individual terms). The *latex_series* should have each term seperated by ' + ' and powers entered explicitly Eg '1 + T_1 + T_1 T_1 + \\\\\\\\frac{1}{2} T_2 T_2 T_1'\n",
    "\n",
    "As an example Linear CCD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "d586bb0a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "D_{1}     np.einsum('abij->abij', gs[v,v,o,o])            \n",
      "D_{2a}    np.einsum('bc,acij->abij', fs[v,v], td)  \\hat{P}(ab) \n",
      "D_{2b}    -np.einsum('kj,abik->abij', fs[o,o], td)  \\hat{P}(ij) \n",
      "D_{2c}    0.5 * np.einsum('abcd,cdij->abij', gs[v,v,v,v], td)            \n",
      "D_{2d}    0.5 * np.einsum('klij,abkl->abij', gs[o,o,o,o], td)            \n",
      "D_{2e}    np.einsum('akic,cbkj->abij', gs[v,o,o,v], td)  \\hat{P}(ab) \\hat{P}(ij) \n"
     ]
    }
   ],
   "source": [
    "ld = expansion(order='D', series='1 + T_2')\n",
    "ld.order_by_reference()\n",
    "for i in ld.ade:\n",
    "    print('{:6}    {:20}  {:10}'.format(i['reference'], i['code'], i['permutation']))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00673f4b",
   "metadata": {},
   "source": [
    "Try\n",
    "\n",
    "```expansion(order='D', series='1 + T_2')``` for LCCD\n",
    "\n",
    "```expansion(order='S', series='1 + T_1 + T_2')``` and ```expansion(order='D', series='1 + T_1 + T_2')```\n",
    "for LCCSD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03bd0238",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}