{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "59e02e86-0f75-481b-9501-31d21e43eb84",
   "metadata": {},
   "source": [
    "## Initialization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b190bf82-0b95-4eb6-b7fb-e78973d5a0b3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PLAMS working folder: /path/plams/examples/ADFFrag/plams_workdir\n"
     ]
    }
   ],
   "source": [
    "from scm.plams import Settings, Molecule, init, AMSJob, Units\n",
    "from scm.plams.recipes.adffragment import ADFFragmentJob\n",
    "\n",
    "# this line is not required in AMS2025+\n",
    "init()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53117bc2-7848-43f3-9ffa-816450e36cd2",
   "metadata": {},
   "source": [
    "## Define the molecules\n",
    "For convenience we define here two molecules, normally you would read them from xyz files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "5c91b68f-e152-4128-b037-18de1c6f99ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_molecule(input_string):\n",
    "    job = AMSJob.from_input(input_string)\n",
    "    return job.molecule[\"\"]\n",
    "\n",
    "\n",
    "mol1 = get_molecule(\n",
    "    \"\"\"\n",
    "System\n",
    "    Atoms\n",
    "        C      -0.75086900       1.37782400      -2.43303700\n",
    "        C      -0.05392100       2.51281000      -2.41769100\n",
    "        H      -1.78964800       1.33942600      -2.09651100\n",
    "        H      -0.30849400       0.43896500      -2.76734700\n",
    "        H      -0.49177100       3.45043100      -2.06789100\n",
    "        H       0.98633900       2.54913500      -2.74329400\n",
    "    End\n",
    "End\n",
    "\"\"\"\n",
    ")\n",
    "\n",
    "\n",
    "mol2 = get_molecule(\n",
    "    \"\"\"\n",
    "System\n",
    "    Atoms\n",
    "        C       0.14667300      -0.21503500       0.40053800\n",
    "        C       1.45297400      -0.07836900       0.12424400\n",
    "        C       2.23119700       1.15868100       0.12912100\n",
    "        C       1.78331500       2.39701500       0.38779700\n",
    "        H      -0.48348000       0.63110600       0.67664100\n",
    "        H      -0.33261900      -1.19332100       0.35411600\n",
    "        H       2.01546300      -0.97840100      -0.14506700\n",
    "        H       3.29046200       1.03872500      -0.12139700\n",
    "        H       2.45728900       3.25301000       0.35150400\n",
    "        H       0.74193400       2.60120700       0.64028800\n",
    "    End\n",
    "End\n",
    "\"\"\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "439d8f35-16df-4d72-b84b-eb8f9e0057a7",
   "metadata": {},
   "source": [
    "## Setup and run the job"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "fd9005ee-f8d6-4eb9-8c5d-b36585748474",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[25.03|17:15:00] JOB plamsjob STARTED\n",
      "[25.03|17:15:00] JOB plamsjob RUNNING\n",
      "[25.03|17:15:00] JOB plamsjob/frag1 STARTED\n",
      "[25.03|17:15:00] JOB plamsjob/frag1 RUNNING\n",
      "[25.03|17:15:11] JOB plamsjob/frag1 FINISHED\n",
      "[25.03|17:15:12] JOB plamsjob/frag1 SUCCESSFUL\n",
      "[25.03|17:15:12] JOB plamsjob/frag2 STARTED\n",
      "[25.03|17:15:12] JOB plamsjob/frag2 RUNNING\n",
      "[25.03|17:15:23] JOB plamsjob/frag2 FINISHED\n",
      "[25.03|17:15:23] JOB plamsjob/frag2 SUCCESSFUL\n",
      "[25.03|17:15:23] JOB plamsjob/full STARTED\n",
      "[25.03|17:15:23] JOB plamsjob/full RUNNING\n",
      "[25.03|17:15:47] JOB plamsjob/full FINISHED\n",
      "[25.03|17:15:47] JOB plamsjob/full SUCCESSFUL\n",
      "[25.03|17:15:47] JOB plamsjob FINISHED\n",
      "[25.03|17:15:47] JOB plamsjob SUCCESSFUL\n"
     ]
    }
   ],
   "source": [
    "common = Settings()  # common settings for all 3 jobs\n",
    "common.input.ams.Task = \"SinglePoint\"\n",
    "common.input.adf.basis.type = \"DZP\"\n",
    "common.input.adf.xc.gga = \"PBE\"\n",
    "common.input.adf.symmetry = \"NOSYM\"\n",
    "\n",
    "full = Settings()  # additional settings for full system calculation\n",
    "full.input.adf.etsnocv  # empty block\n",
    "full.input.adf.print = \"etslowdin\"\n",
    "\n",
    "# normally you would read here the two molecules from xyz files.\n",
    "# mol1 = Molecule(\"ethene.xyz\")\n",
    "# mol2 = Molecule(\"butadiene.xyz\")\n",
    "\n",
    "j = ADFFragmentJob(fragment1=mol1, fragment2=mol2, settings=common, full_settings=full)\n",
    "r = j.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f165ff0-56e0-4c36-a3d7-10415e78e48d",
   "metadata": {},
   "source": [
    "## Print the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c4e090fe-603f-4e5e-a338-a9000261f98c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_eterm(energy_term, energy):\n",
    "    print(\n",
    "        f'{energy_term:>30s} {energy:16.4f} {Units.convert(energy, \"au\", \"eV\"):16.3f} {Units.convert(energy, \"au\", \"kcal/mol\"):16.2f} {Units.convert(energy, \"au\", \"kJ/mol\"):16.2f}'\n",
    "    )\n",
    "\n",
    "\n",
    "def print_bonding_energy_terms(r):\n",
    "    print(\"Energy terms contributing to the bond energy (with respect to the fragments):\")\n",
    "\n",
    "    bond_energy = r.get_energy()\n",
    "    decom = r.get_energy_decomposition()\n",
    "    print(f'\\n{\"term\":>30s} {\"Hartree\":>16s} {\"eV\":>16s} {\"kcal/mol\":>16s} {\"kJ/mol\":>16s}')\n",
    "    for energy_term, energy in decom.items():\n",
    "        print_eterm(energy_term, energy)\n",
    "\n",
    "    print_eterm(\"total bond energy\", bond_energy)\n",
    "    print(\"\")\n",
    "\n",
    "\n",
    "def print_eda_terms(job):\n",
    "    bond_energy = job.full.results.readrkf(\"Energy\", \"Bond Energy\", \"adf\")\n",
    "    steric_interaction = job.full.results.readrkf(\"Energy\", \"Steric Total\", \"adf\")\n",
    "    orbital_interaction = job.full.results.readrkf(\"Energy\", \"Orb.Int. Total\", \"adf\")\n",
    "    print(\"\\nFragment based energy decomposition analysis of the bond energy:\")\n",
    "    print(f'\\n{\"term\":>30s} {\"Hartree\":>16s} {\"eV\":>16s} {\"kcal/mol\":>16s} {\"kJ/mol\":>16s}')\n",
    "    print_eterm(\"Steric interaction\", steric_interaction)\n",
    "    print_eterm(\"Orbital interaction\", orbital_interaction)\n",
    "    print_eterm(\"total bond energy\", bond_energy)\n",
    "    print(\"\")\n",
    "\n",
    "\n",
    "def print_nocv_decomposition():\n",
    "    print(\"NOCV decomposition of the orbital interaction term\\n\")\n",
    "\n",
    "    print(\"The NOCV eigenvalues are occupation numbers, they should come in pairs,\")\n",
    "    print(\"with one negative value mirrored by a positive value.\")\n",
    "    print(\"The orbital interaction energy contribution is calculated for each NOCV pair.\")\n",
    "    print(\"\")\n",
    "\n",
    "    nocv_eigenvalues = j.full.results.readrkf(\"NOCV\", \"NOCV_eigenvalues_restricted\", \"engine\")\n",
    "    nocv_orbitalinteraction = j.full.results.readrkf(\"NOCV\", \"NOCV_oi_restricted\", \"engine\")\n",
    "\n",
    "    n_pairs = int(len(nocv_eigenvalues) / 2)\n",
    "    threshold = 0.001\n",
    "\n",
    "    print(f'{\"index\":>9s} {\"neg\":>9s} {\"pos\":>9s} {\"kcal/mol\":>10s}')\n",
    "    for index in range(n_pairs):\n",
    "        pop1 = nocv_eigenvalues[index]\n",
    "        pop2 = nocv_eigenvalues[len(nocv_eigenvalues) - index - 1]\n",
    "\n",
    "        if (abs(pop1) + abs(pop2)) < threshold:\n",
    "            continue\n",
    "\n",
    "        orbitalinteraction = (\n",
    "            nocv_orbitalinteraction[index] + nocv_orbitalinteraction[len(nocv_orbitalinteraction) - index - 1]\n",
    "        )\n",
    "        print(f\"{index:9d} {pop1:9.3f} {pop2:9.3f} {orbitalinteraction:10.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b5914a0c-1a36-42ac-a776-a04965c28c4d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Energy terms contributing to the bond energy (with respect to the fragments):\n",
      "\n",
      "                          term          Hartree               eV         kcal/mol           kJ/mol\n",
      "          Electrostatic Energy          -0.0059           -0.159            -3.67           -15.38\n",
      "                Kinetic Energy          -0.0109           -0.296            -6.82           -28.53\n",
      "            Elstat Interaction           0.0275            0.749            17.27            72.24\n",
      "                     XC Energy          -0.0131           -0.356            -8.21           -34.37\n",
      "             total bond energy          -0.0023           -0.062            -1.44            -6.02\n",
      "\n",
      "\n",
      "Fragment based energy decomposition analysis of the bond energy:\n",
      "\n",
      "                          term          Hartree               eV         kcal/mol           kJ/mol\n",
      "            Steric interaction           0.0010            0.028             0.64             2.68\n",
      "           Orbital interaction          -0.0033           -0.090            -2.08            -8.70\n",
      "             total bond energy          -0.0023           -0.062            -1.44            -6.02\n",
      "\n",
      "NOCV decomposition of the orbital interaction term\n",
      "\n",
      "The NOCV eigenvalues are occupation numbers, they should come in pairs,\n",
      "with one negative value mirrored by a positive value.\n",
      "The orbital interaction energy contribution is calculated for each NOCV pair.\n",
      "\n",
      "    index       neg       pos   kcal/mol\n",
      "        0    -0.098     0.098      -0.65\n",
      "        1    -0.084     0.084      -0.76\n",
      "        2    -0.045     0.045      -0.38\n",
      "        3    -0.014     0.014      -0.06\n",
      "        4    -0.012     0.012      -0.04\n",
      "        5    -0.012     0.012      -0.04\n",
      "        6    -0.010     0.010      -0.03\n",
      "        7    -0.008     0.008      -0.02\n",
      "        8    -0.008     0.008      -0.02\n",
      "        9    -0.006     0.006      -0.01\n",
      "       10    -0.006     0.006      -0.01\n",
      "       11    -0.006     0.006      -0.01\n",
      "       12    -0.005     0.005      -0.01\n",
      "       13    -0.004     0.004      -0.01\n",
      "       14    -0.003     0.003      -0.00\n",
      "       15    -0.003     0.003      -0.00\n",
      "       16    -0.002     0.002      -0.00\n"
     ]
    }
   ],
   "source": [
    "print_bonding_energy_terms(r)\n",
    "\n",
    "print_eda_terms(j)\n",
    "\n",
    "print_nocv_decomposition()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd5970dd-9910-494d-9b24-d7d638d3924c",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}