QeqEwald.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "from math import erfc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\"\"\"This section reads a geometry\"\"\"\n",
    "\n",
    "def vol(alpha,beta,gamma):\n",
    "\n",
    "    \"\"\" Calculate volume \"\"\"\n",
    "\n",
    "    v = math.sqrt(1-math.cos(alpha)**2-math.cos(beta)**2-math.cos(gamma)**2+2*math.cos(alpha)*math.cos(beta)*math.cos(gamma))\n",
    "    return v\n",
    "    \n",
    "def red2cart(xr,yr,zr,aa,bb,cc,alpha,beta,gamma):\n",
    "\n",
    "    \"\"\" Convert fractional coordinates to cartesian\"\"\"\n",
    "\n",
    "    vv = vol(alpha,beta,gamma)\n",
    "    x = aa*xr+bb*math.cos(gamma)*yr+cc*math.cos(beta)*zr\n",
    "    y = bb*math.sin(gamma)*yr+zr*cc*(math.cos(alpha)-math.cos(beta)*math.cos(gamma))/math.sin(gamma)\n",
    "    z = zr*cc*vv/math.sin(gamma)\n",
    "\n",
    "    return x, y, z\n",
    "\n",
    "def ang(x, y):\n",
    "    \n",
    "    \"\"\" Computes angle between 2 vectors \"\"\"\n",
    "    \n",
    "    return math.acos(np.dot(x,y)/np.dot(x,x))\n",
    "    \n",
    "def read_poscar(filename):\n",
    "    \n",
    "    \"\"\" Read POSCAR format \"\"\"\n",
    "    \n",
    "    fin = open(filename,\"r\")\n",
    "    r = fin.readlines()\n",
    "    scale_pos = float(r[1])\n",
    "    cell0 = []\n",
    "    for i in range(3):\n",
    "        for j in range(3):\n",
    "            cell0.append(float(r[2+i].split()[j])*scale_pos)\n",
    "    ntypes = int(len(r[5].split()))\n",
    "\n",
    "    t,nat = [],[]\n",
    "    for i in range(ntypes):\n",
    "        t.append(r[5].split()[i])\n",
    "    for i in range(len(r[6].split())):\n",
    "        nat.append(int(r[6].split()[i]))\n",
    "\n",
    "    a = [cell0[0],cell0[1],cell0[2]]\n",
    "    b = [cell0[3],cell0[4],cell0[5]]\n",
    "    c = [cell0[6],cell0[7],cell0[8]]\n",
    "    aa = math.sqrt(np.dot(a,a))\n",
    "    bb = math.sqrt(np.dot(b,b))\n",
    "    cc = math.sqrt(np.dot(c,c))\n",
    "    alpha = ang(b,c)\n",
    "    beta = ang(c,a)\n",
    "    gamma = ang(a,b)\n",
    "    cell = [aa,bb,cc,alpha,beta,gamma]\n",
    "\n",
    "    if r[7].split()[0] == \"Selective\":\n",
    "        start = 9\n",
    "    else:\n",
    "        start = 8\n",
    "\n",
    "    ctype = r[start-1].split()[0]\n",
    "\n",
    "    if ctype == \"Direct\":\n",
    "        xr,yr,zr,x,y,z=[],[],[],[],[],[]\n",
    "        for i in range(start,start+sum(nat)):\n",
    "            xr.append(float(r[i].split()[0]))\n",
    "            yr.append(float(r[i].split()[1]))\n",
    "            zr.append(float(r[i].split()[2]))\n",
    "        for i in range(len(xr)):\n",
    "            xx,yy,zz = red2cart(xr[i],yr[i],zr[i],aa,bb,cc,alpha,beta,gamma)\n",
    "            x.append(xx)\n",
    "            y.append(yy)\n",
    "            z.append(zz)\n",
    "\n",
    "    if ctype == \"Cartesian\":\n",
    "        x,y,z=[],[],[]\n",
    "        for i in range(start,len(r)):\n",
    "            x.append(float(r[i].split()[0]))\n",
    "            y.append(float(r[i].split()[1]))\n",
    "            z.append(float(r[i].split()[2]))\n",
    "\n",
    "        fin.close()\n",
    "    return nat,t,x,y,z,cell0,cell,ctype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nat_list,t,x,y,z,cell0,cell,ctype = read_poscar(\"SiC_Nico_cart\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def EEqeq(x,y,z,t,nat,a1,a2,a3):\n",
    "    \n",
    "    \"\"\" This compute EQEq according to Ref. [] \"\"\"\n",
    "\n",
    "    Nmax =2\n",
    "    K = 14.4\n",
    "    zeta = 0.27                           # denotes the Ewald splitting parameter, choice of it is made based on your system\n",
    "    lmda = 1.2\n",
    "    chi = {\"Si\":1.8038,\"C\":5.7254}       # denotes the electronegativity in eV\n",
    "    eta = {\"Si\":14.7704,\"C\":13.847}   # Eta denotes the self-Coulomb potential in eV\n",
    "    gamma = {\"Si\":0.8218,\"C\":0.8712}     # denotes shielded Coulomb constant in distance units. \n",
    "    \n",
    "    \"\"\"Calculation of the K-space vectors\"\"\"\n",
    "    vf = 2*np.pi/np.dot(a1,np.cross(a2,a3))\n",
    "    b1 = 2*np.pi*np.cross(a2,a3)/np.dot(a1,np.cross(a2,a3))\n",
    "    b2 = 2*np.pi*np.cross(a3,a1)/np.dot(a2,np.cross(a3,a1))\n",
    "    b3 = 2*np.pi*np.cross(a1,a2)/np.dot(a3,np.cross(a1,a2))\n",
    "    print(b1,b2,b3)\n",
    "\n",
    "    def JEEQeq(i,j,x,y,z,t,nat,a1,a2,a3): \n",
    "\n",
    "        if i==j: # Looking at one atom interactin with itslef\n",
    "            Astar = 0\n",
    "            Bstar = 0\n",
    "            for u in range(-Nmax,Nmax+1):\n",
    "                for v in range(-Nmax,Nmax+1):\n",
    "                    for w in range(-Nmax,Nmax+1):\n",
    "                        if not (u==0 and v==0 and w==0): # Interacting with its periodic images\n",
    "                            dx = u*a1[0]+v*a2[0]+w*a3[0]\n",
    "                            dy = u*a1[1]+v*a2[1]+w*a3[1]\n",
    "                            dz = u*a1[2]+v*a2[2]+w*a3[2]\n",
    "                            dh1 = u*b1[0]+v*b2[0]+w*b3[0]   \n",
    "                            dh2 = u*b1[1]+v*b2[1]+w*b3[1]\n",
    "                            dh3 = u*b1[2]+v*b2[2]+w*b3[2]\n",
    "                            dr = np.sqrt(dx**2+dy**2+dz**2)\n",
    "                            dh = np.sqrt(dh1**2+dh2**2+dh3**2)\n",
    "                            Astar += erfc(dr*np.sqrt(zeta))*1./dr\n",
    "                            Bstar += (2*vf/dh**2)*np.exp(-(dh*0.5)**2/zeta)\n",
    "            return eta[t[i]]+0.5*lmda*K*(Astar+Bstar-2*np.sqrt(zeta/np.pi))\n",
    "        \n",
    "        else: # Looking at one atom interacting with others\n",
    "            Alpha = 0\n",
    "            Beta = 0\n",
    "            for u in range(-Nmax,Nmax+1):\n",
    "                for v in range(-Nmax,Nmax+1):\n",
    "                    for w in range(-Nmax,Nmax+1): # Interacting with any other atom in any other cell\n",
    "                        xi,yi,zi = x[i],y[i],z[i]\n",
    "                        xj,yj,zj = x[j],y[j],z[j]\n",
    "                        dx = xj-xi+u*a1[0]+v*a2[0]+w*a3[0]\n",
    "                        dy = yj-yi+u*a1[1]+v*a2[1]+w*a3[1]\n",
    "                        dz = zj-zi+u*a1[2]+v*a2[2]+w*a3[2]\n",
    "                        dr = np.sqrt(dx**2+dy**2+dz**2)\n",
    "                        Alpha += erfc(dr*np.sqrt(zeta))*1/dr\n",
    "                        if not (u==0 and v==0 and w==0): # Interacting with other atoms in other cell than itself\n",
    "                            dh1 = u*b1[0]+v*b2[0]+w*b3[0]   \n",
    "                            dh2 = u*b1[1]+v*b2[1]+w*b3[1]\n",
    "                            dh3 = u*b1[2]+v*b2[2]+w*b3[2]\n",
    "                            xi,yi,zi = x[i],y[i],z[i]\n",
    "                            xj,yj,zj = x[j],y[j],z[j]\n",
    "                            dx = xj-xi \n",
    "                            dy = yj-yi \n",
    "                            dz = zj-zi\n",
    "                            dh = np.sqrt(dh1**2+dh2**2+dh3**2)\n",
    "                            Beta += (2*vf/dh**2)*np.exp(-(dh*0.5)**2/zeta)*np.cos(dx*dh1+dy*dh2+dz*dh3)\n",
    "            return 0.5*lmda*K*(Alpha+Beta)\n",
    "    Qtot = 0.0    \n",
    "    b = np.array([0.0]*nat)\n",
    "    b[0] = Qtot\n",
    "    A = np.matrix([[0.0]*nat]*nat)\n",
    "    A[0] = [1]*nat  # First row = 1\n",
    "    for i in range(1,nat):\n",
    "        b[i] = chi[t[i]]-chi[t[0]] \n",
    "        for j in range(nat):\n",
    "            A[i,j] = JEEQeq(i,j,x,y,z,t,nat,a1,a2,a3)-JEEQeq(0,j,x,y,z,t,nat,a1,a2,a3)\n",
    "        print(i,nat)\n",
    "            \n",
    "    #print(\"This is the A matrix\")\n",
    "    #print(A)\n",
    "    #print(\"This is the b vector\")\n",
    "    #print(b)\n",
    "    Q = np.linalg.solve(A, -b)\n",
    "    \n",
    "    print(\"This is the partial charge of each atom type\")\n",
    "    for i in range(nat):\n",
    "        print(t[i],Q[i])\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "type_list = []\n",
    "for i in range(len(t)):\n",
    "    type_list += [t[i]]*nat_list[i]\n",
    "\n",
    "#print(type_list)\n",
    "a1 = cell0[:3]\n",
    "a2 = cell0[3:6]\n",
    "a3 = cell0[6:]\n",
    "nat = sum(nat_list)\n",
    "\n",
    "#print(nat,type_list,a1,a2,a3)\n",
    "EEqeq(x,y,z,type_list,nat,a1,a2,a3)"
   ]
  }
 ],
 "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.8"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}