PsiCon_2023.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import qcelemental as qcel\n",
    "import molsym\n",
    "np.set_printoptions(suppress=True, precision=3, linewidth=600)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MolSym Molecule:\n",
      "   N       0.00000000     0.00000000     0.22175369\n",
      "   H       0.00000000     1.76245118    -0.51742402\n",
      "   H       1.52632801    -0.88122654    -0.51742402\n",
      "   H      -1.52632801    -0.88122654    -0.51742402\n",
      "\n",
      "Point Group is C3v\n"
     ]
    }
   ],
   "source": [
    "# Create MolSym Molecule from xyz file\n",
    "\n",
    "fn = \"test/xyz/ammonia.xyz\"\n",
    "#fn = \"test/xyz/methane.xyz\"\n",
    "#fn = \"test/xyz/benzene.xyz\"\n",
    "\n",
    "mol = molsym.Molecule.from_file(fn)\n",
    "print(mol)\n",
    "\n",
    "pg, (paxis, saxis) = molsym.find_point_group(mol)\n",
    "print(f\"Point Group is {pg}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Point Group is C3v\n",
      "Primary axis for orienting molecule: [0. 0. 1.]\n",
      "Secondary axis: [-0.866 -0.5    0.   ]\n"
     ]
    }
   ],
   "source": [
    "# Create Molecule from Schema\n",
    "\n",
    "with open(fn, \"r\") as lfn:\n",
    "    file_str = lfn.read()\n",
    "\n",
    "schema = qcel.models.Molecule.from_data(file_str).dict()\n",
    "mol = molsym.Molecule.from_schema(schema)\n",
    "\n",
    "pg, (paxis, saxis) = molsym.find_point_group(mol)\n",
    "\n",
    "print(f\"Point Group is {pg}\")\n",
    "print(f\"Primary axis for orienting molecule: {paxis}\")\n",
    "print(f\"Secondary axis: {saxis}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The default tolerance for detecting symmetry is: 1e-05\n",
      "The point group for a tolerance of 1e-12 is Cs\n",
      "The point group for a tolerance of 1e-12 is C3v after symmetrization\n"
     ]
    }
   ],
   "source": [
    "# The coordinates in test/xyz/ammonia.xyz are not symmetric to a tight tolerance\n",
    "\n",
    "print(f\"The default tolerance for detecting symmetry is: {mol.tol}\")\n",
    "\n",
    "# Increase tolerance\n",
    "mol.tol = 1e-12\n",
    "pg, (paxis, saxis) = molsym.find_point_group(mol)\n",
    "print(f\"The point group for a tolerance of {mol.tol} is {pg}\")\n",
    "mol.tol = 1e-5\n",
    "\n",
    "# molsym.symmetrize uses a loose tolerance to detect a potential point group, \n",
    "# then forces the Molecule to that symmetry to near machine precision\n",
    "mol = molsym.symmetrize(mol)\n",
    "pg, (paxis, saxis) = molsym.find_point_group(mol)\n",
    "print(f\"The point group for a tolerance of {mol.tol} is {pg} after symmetrization\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[\n",
       " Symbol:          E: [[ 1.00000  0.00000  0.00000],[ 0.00000  1.00000  0.00000],[ 0.00000  0.00000  1.00000]],\n",
       " \n",
       " Symbol:      C_3^1: [[-0.50000 -0.86603  0.00000],[ 0.86603 -0.50000  0.00000],[ 0.00000  0.00000  1.00000]],\n",
       " \n",
       " Symbol:      C_3^2: [[-0.50000  0.86603  0.00000],[-0.86603 -0.50000  0.00000],[ 0.00000  0.00000  1.00000]],\n",
       " \n",
       " Symbol: sigma_v(1): [[ 1.00000  0.00000 -0.00000],[ 0.00000 -1.00000  0.00000],[-0.00000  0.00000  1.00000]],\n",
       " \n",
       " Symbol: sigma_v(2): [[-0.50000 -0.86603  0.00000],[-0.86603  0.50000  0.00000],[ 0.00000  0.00000  1.00000]],\n",
       " \n",
       " Symbol: sigma_v(3): [[-0.50000  0.86603  0.00000],[ 0.86603  0.50000 -0.00000],[ 0.00000 -0.00000  1.00000]]]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Given a valid Schönflies symbol as a string, the symmetry elements (symels) can be generated\n",
    "molsym.symtext.main.pg_to_symels(pg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Character Table for C3v\n",
       "Irreps: ['A1', 'A2', 'E']\n",
       "Classes: ['E', '2C_3', '3sigma_v']\n",
       "Characters:\n",
       "[[ 1.  1.  1.]\n",
       " [ 1.  1. -1.]\n",
       " [ 2. -1.  0.]]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Similarly, the real valued character table of any point group can be generated\n",
    "molsym.symtext.main.pg_to_chartab(pg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "MolSym Molecule:\n",
      "   N      -0.00000000     0.00000000     0.13125886\n",
      "   H      -0.88122565    -1.52632759    -0.60791885\n",
      "   H      -0.88122565     1.52632759    -0.60791885\n",
      "   H       1.76245129    -0.00000000    -0.60791885\n",
      "\n",
      "Character Table for C3v\n",
      "Irreps: ['A1', 'A2', 'E']\n",
      "Classes: ['E', '2C_3', '3sigma_v']\n",
      "Characters:\n",
      "[[ 1.  1.  1.]\n",
      " [ 1.  1. -1.]\n",
      " [ 2. -1.  0.]]\n",
      "\n",
      "[\n",
      "Symbol:          E: [[ 1.00000  0.00000  0.00000],[ 0.00000  1.00000  0.00000],[ 0.00000  0.00000  1.00000]], \n",
      "Symbol:      C_3^1: [[-0.50000 -0.86603  0.00000],[ 0.86603 -0.50000  0.00000],[ 0.00000  0.00000  1.00000]], \n",
      "Symbol:      C_3^2: [[-0.50000  0.86603  0.00000],[-0.86603 -0.50000  0.00000],[ 0.00000  0.00000  1.00000]], \n",
      "Symbol: sigma_v(1): [[ 1.00000  0.00000 -0.00000],[ 0.00000 -1.00000  0.00000],[-0.00000  0.00000  1.00000]], \n",
      "Symbol: sigma_v(2): [[-0.50000 -0.86603  0.00000],[-0.86603  0.50000  0.00000],[ 0.00000  0.00000  1.00000]], \n",
      "Symbol: sigma_v(3): [[-0.50000  0.86603  0.00000],[ 0.86603  0.50000 -0.00000],[ 0.00000 -0.00000  1.00000]]]\n",
      "Class map:\n",
      "[0 1 1 2 2 2]\n",
      "Atom map:\n",
      "[[0 0 0 0 0 0]\n",
      " [1 3 2 2 3 1]\n",
      " [2 1 3 1 2 3]\n",
      " [3 2 1 3 1 2]]\n",
      "Multiplication Table\n",
      "[[0 1 2 3 4 5]\n",
      " [1 2 0 5 3 4]\n",
      " [2 0 1 4 5 3]\n",
      " [3 4 5 0 1 2]\n",
      " [4 5 3 2 0 1]\n",
      " [5 3 4 1 2 0]]\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "Most of the necessary symmetry information is stored in the Symtext object, such as:\n",
    "    - A molecule rotated to the symmetry element reference frame\n",
    "    - Matrices for rotating Cartesian vectors to and from the symmetry element ref. frame\n",
    "    - List of symmetry elements\n",
    "    - Character table object\n",
    "    - Map of all atoms in the Molecule under each symmetry operation\n",
    "    - Multiplication table of the symmetry operations\n",
    "\"\"\"\n",
    "\n",
    "symtext = molsym.Symtext.from_molecule(mol)\n",
    "print(symtext)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SALC from P^A1_00 (2) gamma= 2.373\n",
      "[-0.     0.     0.421  0.    -0.    -0.524  0.    -0.    -0.524  0.    -0.    -0.524]\n",
      "SALC from P^A1_00 (3) gamma= 3.464\n",
      "[-0.    -0.     0.     0.289  0.5    0.     0.289 -0.5   -0.    -0.577 -0.    -0.   ]\n",
      "SALC from P^E_00 (0) gamma= 2.695\n",
      "[ 0.371 -0.    -0.    -0.461  0.    -0.193 -0.461  0.    -0.193 -0.461  0.     0.387]\n",
      "SALC from P^E_10 (0) gamma= 2.695\n",
      "[ 0.     0.371  0.    -0.    -0.461 -0.335 -0.    -0.461  0.335 -0.    -0.461  0.   ]\n",
      "SALC from P^E_00 (3) gamma= 1.838\n",
      "[-0.315  0.     0.     0.544 -0.265  0.164  0.544  0.265  0.164  0.084  0.    -0.328]\n",
      "SALC from P^E_10 (3) gamma= 1.838\n",
      "[-0.    -0.315 -0.    -0.265  0.238  0.284  0.265  0.238 -0.284  0.     0.697 -0.   ]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "The main goal of MolSym is to return Symmetry Adapted Linear Combinations of arbitrary functions.\n",
    "Built in are SALCs of Cartesian coordinates, internal coordinates, and atomic basis functions.\n",
    "\"\"\"\n",
    "\n",
    "cart_coords = molsym.salcs.CartesianCoordinates(symtext)\n",
    "salcs = molsym.salcs.ProjectionOp(symtext, cart_coords)\n",
    "print(salcs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "B = salcs.basis_transformation_matrix\n",
    "print(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SALC from P^A1_00 (0) gamma= 1.732\n",
      "[0.577 0.577 0.577 0.    0.    0.   ]\n",
      "SALC from P^A1_00 (3) gamma= 1.732\n",
      "[0.    0.    0.    0.577 0.577 0.577]\n",
      "SALC from P^E_00 (0) gamma= 2.449\n",
      "[ 0.408  0.408 -0.816  0.     0.     0.   ]\n",
      "SALC from P^E_10 (0) gamma= 2.449\n",
      "[ 0.707 -0.707  0.     0.     0.     0.   ]\n",
      "SALC from P^E_00 (3) gamma= 1.225\n",
      "[ 0.     0.     0.     0.816 -0.408 -0.408]\n",
      "SALC from P^E_10 (3) gamma= 1.225\n",
      "[ 0.     0.     0.     0.     0.707 -0.707]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Ignore how we define internal coordinates, this is subject to change\n",
    "ic_coords = [[0,1],[0,2],[0,3],[1,0,2],[1,0,3],[2,0,3]]\n",
    "ic_names = [\"R1\", \"R2\", \"R3\", \"A1\", \"A2\", \"A3\"]\n",
    "ic_list = [[ic_coords[i], ic_names[i]] for i in range(len(ic_coords))]\n",
    "\n",
    "# SALCs of internal coordinates\n",
    "ics = molsym.salcs.InternalCoordinates(ic_list, symtext)\n",
    "salcs = molsym.salcs.ProjectionOp(symtext, ics)\n",
    "print(salcs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Symmetry",
   "language": "python",
   "name": "symmetry"
  },
  "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.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}