molecular_polarizability.pynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5dea8300",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import os,sys\n",
    "from analyzetool import process, prmedit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c729f59c",
   "metadata": {},
   "outputs": [],
   "source": [
    "header = \"\"\"memory 30 GB\n",
    "set {\n",
    "    basis         aug-cc-pVTZ\n",
    "    basis_guess   False\n",
    "    scf_type      df\n",
    "    cc_type       df\n",
    "    opt_type      min\n",
    "    dynamic_level 1\n",
    "    maxiter       200\n",
    "    guess sad\n",
    "    \n",
    "    d_convergence   10\n",
    "}\n",
    "\n",
    "set_num_threads(6)\n",
    "method = 'MP2'\n",
    "\"\"\"\n",
    "\n",
    "tail_ref = f\"\"\"\n",
    "optimize(method)\n",
    "\n",
    "pert = 0.001\n",
    "lambdas = [pert, -pert, 2.0*pert, -2.0*pert]\n",
    "\n",
    "perturbed_energies = np.zeros((len(lambdas),3))\n",
    "x_perturbed_dipoles = np.zeros((len(lambdas),3))\n",
    "y_perturbed_dipoles = np.zeros((len(lambdas),3))\n",
    "z_perturbed_dipoles = np.zeros((len(lambdas),3))\n",
    "\n",
    "# start with a reference dipole calculation\n",
    "properties(method, properties=['dipole'])\n",
    "mu_x, mu_y, mu_z = psi4.core.variable(method + ' DIPOLE') / psi_dipmom_au2debye\n",
    "analytic_dipole = np.array([mu_x, mu_y, mu_z])\n",
    "\n",
    "# now compute with different applied fields\n",
    "for step, l in enumerate(lambdas):\n",
    "    set perturb_h true\n",
    "    set perturb_with dipole\n",
    "    # x pertubation\n",
    "    set perturb_dipole [$l, 0, 0]\n",
    "    perturbed_energies[step,0] = properties(method, properties=['dipole'])\n",
    "    x_perturbed_dipoles[step] = psi4.core.variable(method + ' DIPOLE')\n",
    "\n",
    "    # y pertubation\n",
    "    set perturb_dipole [0, $l, 0]\n",
    "    perturbed_energies[step,1] = properties(method, properties=['dipole'])\n",
    "    y_perturbed_dipoles[step] = psi4.core.variable(method + ' DIPOLE')\n",
    "    # z pertubation\n",
    "    set perturb_dipole [0, 0, $l]\n",
    "    perturbed_energies[step,2] = properties(method, properties=['dipole'])\n",
    "    z_perturbed_dipoles[step] = psi4.core.variable(method + ' DIPOLE')\n",
    "\n",
    "a_x_5pt = (8.0*x_perturbed_dipoles[0] - 8.0*x_perturbed_dipoles[1] - x_perturbed_dipoles[2] + x_perturbed_dipoles[3]) / (12.0*pert)\n",
    "a_y_5pt = (8.0*y_perturbed_dipoles[0] - 8.0*y_perturbed_dipoles[1] - y_perturbed_dipoles[2] + y_perturbed_dipoles[3]) / (12.0*pert)\n",
    "a_z_5pt = (8.0*z_perturbed_dipoles[0] - 8.0*z_perturbed_dipoles[1] - z_perturbed_dipoles[2] + z_perturbed_dipoles[3]) / (12.0*pert)\n",
    "alpha_5pt = np.vstack((a_x_5pt, a_y_5pt, a_z_5pt))\n",
    "\n",
    "eigen = np.linalg.eigvals(alpha_5pt)\n",
    "\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "08a52b9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "basedir = \"/home/roseane/Dropbox/WashU/Projects/HIPPO/small_molecules/org_liq_database\"\n",
    "qmdir = f\"{basedir}/qm_calc/water-mol\"\n",
    "\n",
    "workdir = \"/home/roseane/HIPPO/small_molecules/org_liq_database\"\n",
    "elecdir = f\"{workdir}/elec-2\"\n",
    "poldir = f\"{workdir}/mol-polarize\"\n",
    "atomdef = f\"{workdir}/update_params/atomdef\"\n",
    "\n",
    "for n in range(162,164):\n",
    "    pdir = f\"{poldir}/{n}\"\n",
    "    if not os.path.isdir(pdir):\n",
    "        os.system(f\"mkdir -p {pdir}\")\n",
    "    \n",
    "    f1 = f\"{atomdef}/{n}/monomer.xyz\"\n",
    "    \n",
    "    if not os.path.isfile(f1):\n",
    "        continue\n",
    "        \n",
    "    f3 = f\"{pdir}/{n}-poltensor.npy\"\n",
    "    f2 = f\"{pdir}/input\" \n",
    "\n",
    "#         os.system(f\"rm -f {f3}\")\n",
    "    if os.path.isfile(f3):\n",
    "        continue\n",
    "\n",
    "    mol = \"\"\"\n",
    "molecule mol {\n",
    "0 1\n",
    "\"\"\"\n",
    "    coords,atommap = prmedit.get_coords_from_xyz(f1)\n",
    "    for s,e in zip(coords,atommap):\n",
    "        mol += f\"{e:3s}  {s[0]:9.6f} {s[1]:9.6f} {s[2]:9.6f}\\n\"\n",
    "\n",
    "    mol += \"\"\"\n",
    "no_com\n",
    "no_reorient\n",
    "units angstrom\n",
    "}\n",
    "\"\"\"\n",
    "    tail = tail_ref + f\"\"\"np.save('{n}-poltensor.npy',alpha_5pt)\n",
    "np.save('{n}-poleigen.npy',eigen)\n",
    "\n",
    "mol.save_xyz_file('{n}-opt.xyz',1)\n",
    "\"\"\"\n",
    "    inpfile = header + mol + tail\n",
    "\n",
    "    with open(f2, 'w') as thefile:\n",
    "        thefile.write(inpfile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "e174e3fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "#### Move opt-xyz coordinates to tinker xyz\n",
    "poldir = f\"/home/roseane/HIPPO/small_molecules/org_liq_database/mol-polarize\"\n",
    "n = 43\n",
    "for n in range(162,164):\n",
    "    pdir = f\"{poldir}/{n}\"\n",
    "\n",
    "    f1 = f\"{atomdef}/{n}/monomer.xyz\"\n",
    "    f2 = f\"{pdir}/{n}-opt.xyz\"\n",
    "    f3 = f\"{pdir}/monomer.xyz\"\n",
    "\n",
    "    newcoords,atommap = prmedit.get_coords_from_xyz(f2)\n",
    "    prmedit.update_xyztypes(f1,f3,newcoords)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "755d9eff",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[ 1.15369952e-02,  1.77297659e-03, -4.91779458e-03],\n",
       "        [-3.31150681e-02,  1.08191710e-03,  1.50261090e+00],\n",
       "        [-1.43952828e+00,  1.83047880e-03,  2.04725774e+00],\n",
       "        [ 1.02087804e+00,  1.70941944e+00,  2.24782659e+00],\n",
       "        [ 1.03237685e+00, -5.16106092e-03, -3.75034541e-01],\n",
       "        [-5.01204412e-01, -8.85466937e-01, -3.77761371e-01],\n",
       "        [-4.91074697e-01,  8.81389807e-01, -4.00642788e-01],\n",
       "        [ 5.30441164e-01, -8.36840955e-01,  1.90111827e+00],\n",
       "        [-1.44816762e+00, -4.92463928e-03,  3.13308697e+00],\n",
       "        [-1.96195483e+00, -8.85480062e-01,  1.68828483e+00],\n",
       "        [-1.98018814e+00,  8.81379037e-01,  1.70517120e+00]]),\n",
       " array(['C', 'C', 'C', 'I', 'H', 'H', 'H', 'H', 'H', 'H', 'H'], dtype='<U1'))"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newcoords"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e197f0f",
   "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}