Heisenberg_cupy2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import time\n",
    "import cupy as cp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cross = cp.ElementwiseKernel(\n",
    "     'raw float64 a,raw float64 b, float64 c', '',\n",
    "     '''\n",
    "     int z=i%3;\n",
    "     if (z==0){\n",
    "     c=a[i+1]*b[i+2]-a[i+2]*b[i+1];\n",
    "     }\n",
    "     else if (z==1){\n",
    "     c=a[i+1]*b[i-1]-a[i-1]*b[i+1];\n",
    "     }\n",
    "     else{\n",
    "     c=a[i-2]*b[i-1]-a[i-1]*b[i-2];\n",
    "     }\n",
    "     ''', 'my_kernel')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "L=900\n",
    "s=cp.zeros((L,L,3))\n",
    "n_iter=1001 #10001\n",
    "K=2\n",
    "imp_prob=0.25 #impurity probability\n",
    "R=cp.random.choice(a=[True, False], size=(L, L), p=[imp_prob, 1-imp_prob]) #impurtiy lattice matrix\n",
    "\n",
    "f=open('corr_K2_p0.25.csv','ab')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "white = cp.ones((L, L),dtype='float32')\n",
    "\n",
    "kernel = cp.ElementwiseKernel(\n",
    "     'float32 x, int16 L', 'float32 z',\n",
    "     '''\n",
    "     int b=i%L;\n",
    "     int a=(i-b)/L;\n",
    "     if ((a+b)%2 == 0) {\n",
    "       z = 1;\n",
    "     } \n",
    "     else {\n",
    "       z = 0;\n",
    "     }''', 'my_kernel')\n",
    "\n",
    "white=kernel(white,L)\n",
    "black=1-white"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def monte(s,theta,phi,B,T,K,colour):\n",
    "    J=1\n",
    "    D=6**0.5\n",
    "    beta=1/T\n",
    "    e=0.2#1\n",
    "    \n",
    "    del_theta=e*(cp.random.rand(L,L)-0.5)\n",
    "    del_phi=e*(cp.random.rand(L,L)-0.5)\n",
    "    theta2=(theta+del_theta)%(np.pi) #theta has to be in range (0,pi)\n",
    "    phi2=(phi+del_phi)%(2*np.pi)\n",
    "    \n",
    "    s2=cp.zeros((L,L,3))\n",
    "    s2[:,:,0]=cp.sin(theta2)*cp.cos(phi2)\n",
    "    s2[:,:,1]=cp.sin(theta2)*cp.sin(phi2)\n",
    "    s2[:,:,2]=cp.cos(theta2)\n",
    "    y=cp.array([0,1,0])\n",
    "    x=cp.array([1,0,0])\n",
    "    \n",
    "    sx=cp.roll(s,-1,axis=0)\n",
    "    sx_=cp.roll(s,1,axis=0)\n",
    "    sy=cp.roll(s,-1,axis=1)\n",
    "    sy_=cp.roll(s,1,axis=1)\n",
    "    \n",
    "    E_J=-J*cp.dot(cp.multiply(s2-s,sx+sx_+sy+sy_),cp.ones(3))\n",
    "    E_B=-B*(s2-s)[:,:,2]\n",
    "    t1=cp.zeros((L,L,3))\n",
    "    cross(s2-s,sx-sx_,t1)\n",
    "    t2=cp.zeros((L,L,3))\n",
    "    cross(s2-s,sy-sy_,t2)\n",
    "    E_D=+D*cp.dot(t1,y) - D*cp.dot(t2,x)\n",
    "    E_T=-T*(cp.log(cp.sin(theta2)+0.00001)-cp.log(cp.sin(theta)+0.00001))\n",
    "    E_A=-K*R*(s2[:,:,2]*s2[:,:,2]-s[:,:,2]*s[:,:,2])\n",
    "    del_E=E_J+E_B+E_D+E_T+E_A\n",
    "    \n",
    "    trans_prob= cp.logical_or(del_E<0, cp.exp(-beta*del_E)>cp.random.rand(L,L))\n",
    "    if colour==1:\n",
    "        trans_prob=cp.multiply(trans_prob,white)\n",
    "    else:\n",
    "        trans_prob=cp.multiply(trans_prob,black)\n",
    "\n",
    "    theta=(theta+cp.multiply(trans_prob,del_theta))%(np.pi) #theta has to be in range (0,pi)\n",
    "    phi=(phi+cp.multiply(trans_prob,del_phi))%(2*np.pi)\n",
    "    p=cp.zeros((L,L,3))\n",
    "    p[:,:,0]=s[:,:,0]+cp.multiply(trans_prob,(s2-s)[:,:,0])\n",
    "    p[:,:,1]=s[:,:,1]+cp.multiply(trans_prob,(s2-s)[:,:,1])\n",
    "    p[:,:,2]=s[:,:,2]+cp.multiply(trans_prob,(s2-s)[:,:,2])\n",
    "    \n",
    "    return [p,theta,phi]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "203.8 208.7 212.1 "
     ]
    }
   ],
   "source": [
    "for B in [0,1.45,1.9,295.3,4]:\n",
    "    for n in range(1):\n",
    "        #print(n,end=' ')\n",
    "        theta=cp.random.rand(L,L)*(np.pi)\n",
    "        phi=cp.random.rand(L,L)*(2*np.pi)\n",
    "        \n",
    "        for T in range(101):\n",
    "            cp.zeros((L,L,3))\n",
    "            temp=2*np.exp(-0.04*T)\n",
    "            t = time.time()\n",
    "            for m in range(n_iter):\n",
    "                if T>95:\n",
    "                    s_temp=cp.reshape(s,(1,-1))\n",
    "                    cp.savetxt(f,s_temp,delimiter=',')\n",
    "                [s,theta,phi]=monte(s,theta,phi,B,temp,K,0)\n",
    "                [s,theta,phi]=monte(s,theta,phi,B,temp,K,1)\n",
    "            print(\"{0:.1f}\".format(time.time() - t),end=' ')\n",
    "            \n",
    "    print('\\n',\"Magnetic Field: \",B)\n",
    "    print(\"\\n\")\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}