try1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import random\n",
    "from math import pi, sqrt\n",
    "import sys\n",
    "import os.path\n",
    "import matplotlib.pyplot as plt\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generates single particle states\n",
    "def wf_gen(N,N_pt,BC,t):\n",
    "    hop= np.diag(np.ones(N-1)+0j,1)\n",
    "    hop[N-1,0]= BC\n",
    "    H_t= -t*(hop+ np.matrix(hop).H)/2\n",
    "    energies, evecs= np.linalg.eigh(H_t)\n",
    "    exact_gs_energy=  np.sum(energies[:N_pt])\n",
    "    print('exact gs energy=', exact_gs_energy)\n",
    "#     plt.figure()\n",
    "#     plt.plot(energies,'o')\n",
    "#     plt.show()\n",
    "    return evecs[:,:N_pt]\n",
    "\n",
    "# VMC function \n",
    "def VMC_main(numconfig,t,V1,wf_1,N,N_pt,BC,n_occ1,n_pos1):\n",
    "\n",
    "    move_attempted=0.0\n",
    "    move_accepted=0.0\n",
    "\n",
    "    wf_inv=np.linalg.inv(wf_1)\n",
    "    \n",
    "    count=0 # counter for energy\n",
    "    howoften=10 # calculate energy every 10 steps\n",
    "    min_step=500\n",
    "    ratio_0=1.0+0.0j    \n",
    "    \n",
    "    ratio=np.zeros((numconfig,),dtype=np.complex64)\n",
    "    E_t=np.zeros(int(numconfig/howoften),dtype=np.complex64) # total energy\n",
    "    \n",
    "    for step in range(numconfig+min_step):\n",
    "        for moved_elec in range(N_pt):\n",
    "            move_attempted=move_attempted+1\n",
    "\n",
    "            # random walk of one step left or right\n",
    "            att_num=random.randint(1,2)\n",
    "            if att_num==1:\n",
    "                stepx=1\n",
    "            else:\n",
    "                stepx=-1\n",
    "\n",
    "            ptcls_x= np.mod( n_pos1[moved_elec]+stepx, N) # new configuration\n",
    "\n",
    "            if n_occ1[ptcls_x]==1:\n",
    "                continue\n",
    "\n",
    "            pt_wf_1=np.transpose(V1[ptcls_x,:])\n",
    "\n",
    "            rel=np.dot(wf_inv[moved_elec,:],pt_wf_1)\n",
    "\n",
    "            alpha=min(1, np.abs(rel)**2)\n",
    "            \n",
    "            random_num=random.random()\n",
    "\n",
    "            if random_num <= alpha:\n",
    "                # SMW method to compute inverse of slater det\n",
    "                u_1=pt_wf_1 - wf_1[:,moved_elec]\n",
    "                v=np.zeros((N_pt,1))\n",
    "                v[moved_elec]=1\n",
    "                wf_inv=wf_inv - np.dot(np.dot(np.dot(wf_inv,u_1),v.T),wf_inv) \\\n",
    "                                 /(1+np.dot(v.T,np.dot(wf_inv,u_1))) \n",
    "                wf_1[:,moved_elec]=pt_wf_1\n",
    "#                 wf_inv=np.linalg.inv(wf_1)\n",
    "\n",
    "                move_accepted=move_accepted+1\n",
    "                n_occ1[n_pos1[moved_elec]]= n_occ1[n_pos1[moved_elec]]-1\n",
    "                n_pos1[moved_elec]=ptcls_x\n",
    "                n_occ1[ptcls_x]= n_occ1[ptcls_x]+1\n",
    "                if step> min_step:\n",
    "                    ratio_0=ratio_0*rel\n",
    "\n",
    "# ##############################################################\n",
    "        assert n_occ1.sum()==N_pt, \"total particle number changes!!\"      \n",
    "            \n",
    "        if step> (min_step-1):\n",
    "\n",
    "            if ((step-min_step+1)%howoften)==0:\n",
    "                E_t[count] =energy_calc(t,V1,N,N_pt,BC,wf_inv,n_occ1,n_pos1)\n",
    "                count+=1\n",
    "        \n",
    "        if (step%500) ==0:\n",
    "            wf_inv=np.linalg.inv(wf_1)\n",
    "\n",
    "    acc_ratio=move_accepted/move_attempted\n",
    "    i_c=4\n",
    "    print(\"<H> for iteration %d = \" % (i_c), E_t[i_c])\n",
    "    \n",
    "    plt.figure()\n",
    "    plt.plot(np.real(E_t),'o')\n",
    "    plt.xlabel('iteration No.',fontsize=20)\n",
    "    plt.ylabel(r'$\\frac{H\\Psi}{\\Psi}$',fontsize=20)\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "    print(\"the acceptance rate=\", acc_ratio)\n",
    "\n",
    "\n",
    "# calculates energy \n",
    "def energy_calc(t,V1,N,N_pt,BC,wf_inv,n_occ1,n_pos1):\n",
    "    \n",
    "    E_t=0.0+0.0j\n",
    "    \n",
    "    for moved_elec in range(N_pt):\n",
    "    #### t1 term\n",
    "        ptcls_x= np.mod( (n_pos1[moved_elec])+1, N)\n",
    "        pt_wf_1=np.transpose(V1[ptcls_x,:])\n",
    "        Ecalc=np.dot(wf_inv[moved_elec,:],pt_wf_1)\n",
    "        if n_pos1[moved_elec]== N-1:\n",
    "            E_t += (BC)*Ecalc\n",
    "        else:\n",
    "            E_t += Ecalc\n",
    "\n",
    "        ptcls_x= np.mod( (n_pos1[moved_elec])-1, N)\n",
    "        pt_wf_1=np.transpose(V1[ptcls_x,:])\n",
    "        Ecalc=np.dot(wf_inv[moved_elec,:],pt_wf_1)\n",
    "        if n_pos1[moved_elec]== 0:\n",
    "            E_t += np.conj(BC)*Ecalc\n",
    "        else:\n",
    "            E_t += Ecalc\n",
    "\n",
    "    return -t*E_t/2\n",
    "\n",
    "# initialize wavefunction \n",
    "def initialize_wf(N,N_pt):\n",
    "    n_occ=np.zeros((N,))        \n",
    "    n_pos=random.sample(range(N),N_pt)\n",
    "    n_occ[n_pos]=np.ones((N_pt,))        \n",
    "    \n",
    "    return n_occ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "exact gs energy= -12.745494843182374\n",
      "<H> for iteration 4 =  (-12.745495+9.126979e-16j)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAEPCAYAAADiVdsmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYvUlEQVR4nO3de7RkZX3m8e8jbcR7q4go0DSKEgSNaIniXeQ2GZWLEqNOxJikRUO8JHjBXqzRNRcvZMY442SgRSNeMhq5ZmyxBQ3oci3R0zTYtC1eERswgBlQAVHgN3/sfaQoqs7pPvQ5tU+d72etWnXqffe761fV3efp/e63aqeqkCSpC+4z7gIkSZpmKEmSOsNQkiR1hqEkSeoMQ0mS1BnLxl3AYrbTTjvVypUrx12GJC0q69evv6GqHjmsz1C6F1auXMnU1NS4y5CkRSXJT0b1OX0nSeoMQ0mS1BmGkiSpMwwlSVJnGEqSpM4wlCRJnWEoSZI6w1CSJHWGoSRJ6gxDSZLUGYaSJKkzDCVJUmcYSpKkzjCUJEmdYShJkjrDUJIkdYahJEnqDENJktQZhpIkqTMmJpSSHJNkU5I7k/T62g9Jsj7Jxvb+oBHjT07y3STfTnJ2kuULVrwkCZigUAIuB44GvjrQfgPwkqp6EnAs8MkR488H9quqJwPfA06cr0IlScMtG3cB20tVbQZIMti+oe/hJmDHJPerqtsGtvtS38NvAC+fp1IlSSNM0pHS1ngZsGEwkIZ4HXDesI4kq5JMJZm6/vrrt3uBkrSULaojpSQXALsM6VpdVefOMnZf4P3AobNstxq4Hfj0sP6qWgOsAej1erUVZUuSttKiCqWqOngu45LsBpwNvKaqfjjDdscCLwZeVFUGjiQtsEUVSnPRrqJbC5xYVV+fYbvDgXcAz6+qWxaoPElSn4k5p5TkqCRbgAOBtUnWtV3HA3sBJyW5tL3t3I45rW/5+IeBBwPnt9ucstCvQZKWujhLNXe9Xq+mpqbGXYYkLSpJ1ldVb1jfxBwpSZIWP0NJktQZhpIkqTMMJUlSZxhKkqTOMJQkSZ1hKEmSOsNQkiR1hqEkSeoMQ0mS1BmGkiSpMwwlSVJnGEqSpM4wlCRJnWEoSZI6w1CSJHWGoSRJ6gxDSZLUGYaSJKkzDCVJUmcYSpKkzjCUJEmdMTGhlOSYJJuS3Jmk19d+SJL1STa29wfNsp8TklSSnea/aklSv2XjLmA7uhw4Gjh1oP0G4CVVdU2S/YB1wK7DdpBkd+AQ4Kr5LFSSNNzEHClV1eaqumJI+4aquqZ9uAnYMcn9Ruzmg8DbgZqnMiVJM5iYUNpKLwM2VNVtgx1JXgpcXVWXzbSDJKuSTCWZuv766+erTklakhbV9F2SC4BdhnStrqpzZxm7L/B+4NAhfQ8AVg/rG1RVa4A1AL1ezyMqSdqOFlUoVdXBcxmXZDfgbOA1VfXDIZs8DtgTuCwJwG7AJUkOqKqfzbVeSdK2WVShNBdJlgNrgROr6uvDtqmqjcDOfWOuBHpVdcNC1ChJakzMOaUkRyXZAhwIrE2yru06HtgLOCnJpe1t53bMaf3LxyVJ45UqT4vMVa/Xq6mpqXGXIUmLSpL1VTX0gGBijpQkSYufoSRJ6gxDSZLUGYaSJKkzDCVJUmcYSpKkzjCUJEmdYShJkjrDUJIkdYahJEnqDENJktQZhpIkqTMMJUlSZxhKkqTOMJQkSZ1hKEmSOsNQkiR1hqEkSeoMQ0mS1BmGkiSpMwwlSVJnGEqSpM6YmFBKckySTUnuTNLraz8kyfokG9v7g2bYx18luaLdzwcWpnJJ0rRl4y5gO7ocOBo4daD9BuAlVXVNkv2AdcCug4OTvBA4AnhyVd2WZOf5LliSdHcTE0pVtRkgyWD7hr6Hm4Adk9yvqm4b2MUbgPdNt1fVdfNYriRpiImZvttKLwM2DAkkgCcAz01ycZKLkjx9gWuTpCVvUR0pJbkA2GVI1+qqOneWsfsC7wcOHbHJMuBhwDOBpwP/lOSxVVUD+1kFrAJYsWLFtr0ASdKMFlUoVdXBcxmXZDfgbOA1VfXDEZttAc5qQ+ibSe4EdgKuH6hhDbAGoNfr1T32Ikmas4mfvkuyHFgLnFhVX59h03OAg9oxTwB+j2aRhCRpgUxMKCU5KskW4EBgbZJ1bdfxwF7ASUkubW87t2NO61s+/jHgsUkuBz4DHDs4dSdJml/x9+7c9Xq9mpqaGncZkrSoJFlfVb1hfRNzpCRJWvwMJUlSZxhKkqTOMJQkSZ1hKEmSOsNQkiR1hqEkSeqM7R5KSe7Y3vuUJC0N83GklNk3kSTpnmb8QtYk01+DXcAvq+rGrdinXxEhSZqT2Y6UrgR+3N7/PMmWJP85yf3nuzBJ0tIz26UrXtjeB3gM8DTgzcChSQ6qql/NZ3GSpKVlxlCqqosGmv4xyQeBrwHvBk6Yp7okSUvQ1p5TGnQK8IYkH6L5pvGrtntlkqQlZ7bpuyuZeeHCT9r+HYZ1znGhxEQ7Z8PVnLzuCq658VYes/z+vO2wvTly/12twzqso4M1WMfC1zFbKD19RPtLgT8DjhjS178k/Er6Qi3JtcDHgf9SVbdudZUT4pwNV3PiWRu59bfNR7muvvFWTjxrI8CC/uWyDuvoeh1dqME6xlPHjKvvqmr94A24CngF8Pm+trsN6/v5hTSXGH8R8CfAZ2kWSlyU5EHb5RUsIievu+J3f5jTbv3tHZy87grrsA7r6FgN1jGeOmY7p/S86R+BRwFPBV4P3ACsnm3nLpS4u2tuHH5wOKrdOqxjqdbRhRqsYzx1zPY5pQuBf2lvnwHeCnwBeHZV/Xy2nSdZMXhrn/MU4OVJdp9hMcXEeczy4R/vGtVuHdaxVOvoQg3WMZ46ZgulI2mm315A8xmlh1XVq6vquhnGDJ5T+vGQ238FdqdZKPHjOdS9KL3tsL25/33vvibk/vfdgbcdtrd1WId1dKwG6xhPHbMtdDgb+D6wvr09JMklVfXLGcb0n1Oay0KJiTV9InDcK2iswzq6XkcXarCO8dSRqtErvpO8guYIqQfsDzyEJnR+wF1Btb7/3FGSO6pq6BLxtv+RNOeULqyq47bHixiXXq9XU1NT4y5DkhaVJOurqjesb7ZvdPgszYo5kgS4AzgJeCjNood3AQ9j9OeU7tVCiW2R5BiaxRP7AAdU1VTbfgjwPuD3gN8Ab6uqrwwZ/xSac107ArcDb6yqb27PGiVJM5tt+u53qqqaXOLzVfXt6fYhCxX6zyldyF3TeQF+C5wBvHVrFkpso8uBo4FTB9pvAF5SVdck2Q9YBww71vwA8J6qOi/JH7aPX7Cda5QkzWCrQ2mUIV8xNPg5pem2XwJXVNUt9/Y5R9SxGaANzv72DX0PNwE7JrlfVd02uAua6UlojgSvmY86JUmj3etQmsmQzymN28uADUMCCeAtwLokf0uzKvFZC1mYJGn2D89+EpgCvgVc2jaP7SJ+SS4AdhnStbqqzp1l7L7A+4FDR2zyBpppxTOT/BHwUeDgIftZBawCWLFiyXzESpIWxGyr774CPAVYTrPIYQfg88BFNCvvLqmqXwyMmXH13XxLciFwwvRCh7ZtN+ArwJ9W1ddHjLsJWN6eOwtwU1U9ZNi201x9J0nb7t6svjuo3cFjaZaGT99OBB4OVJIfVtUTtm/J20+S5cBa4MRRgdS6Bng+zeKMg2g+nyVJWkBbdU6pqn4E/Aj43HRbkpU0n1966sDmYQySHAX8T+CRwNokl1bVYcDxwF7ASUlOajc/tKquS3IacEp7VPUXwIeSLAN+TTtFJ0laODNO32lmTt9J0rabafputu++kyRpwRhKkqTOMJQkSZ1hKEmSOsNQkiR1hqEkSeoMQ0mS1BmGkiSpMwwlSVJnGEqSpM4wlCRJnWEoSZI6w1CSJHWGoSRJ6gxDSZLUGYaSJKkzDCVJUmcYSpKkzjCUJEmdYShJkjrDUJIkdYahJEnqjIkJpSTHJNmU5M4kvb72A5Jc2t4uS3LUiPEPT3J+ku+39w9buOolSTBBoQRcDhwNfHVIe6+qngIcDpyaZNmQ8e8EvlxVjwe+3D6WJC2giQmlqtpcVVcMab+lqm5vH+4I1IhdHAGc3v58OnDkdi9SkjSjiQmlmSR5RpJNwEbguL6Q6veoqroWoL3feSFrlCTBsGmszkpyAbDLkK7VVXXuqHFVdTGwb5J9gNOTnFdVv55jDauAVQArVqyYyy4kSSMsqlCqqoPv5fjNSW4G9gOmBrr/Ncmjq+raJI8GrhuxjzXAGoBerzdqKlCSNAcTP32XZM/phQ1J9gD2Bq4csuk/A8e2Px8LjDzykiTNj4kJpSRHJdkCHAisTbKu7XoOcFmSS4GzgTdW1Q3tmNP6lo+/DzgkyfeBQ9rHkqQFlCpnoOaq1+vV1NTgLKAkaSZJ1ldVb1jfxBwpSZIWP0NJktQZhpIkqTMMJUlSZxhKkqTOMJQkSZ1hKEmSOsNQkiR1hqEkSeoMQ0mS1BmGkiSpMwwlSVJnGEqSpM4wlCRJnWEoSZI6w1CSJHWGoSRJ6gxDSZLUGYaSJKkzDCVJUmcYSpKkzjCUJEmdMTGhlOSYJJuS3Jmk19d+QJJL29tlSY4aMf7kJN9N8u0kZydZvmDFS5KACQol4HLgaOCrQ9p7VfUU4HDg1CTLhow/H9ivqp4MfA84cR5rlSQNMTGhVFWbq+qKIe23VNXt7cMdgRox/kt9230D2G1+KpUkjTIxoTSTJM9IsgnYCBzXFz6jvA44b/4rkyT1GzaN1VlJLgB2GdK1uqrOHTWuqi4G9k2yD3B6kvOq6tcjnmM1cDvw6RH9q4BVACtWrNjGVyBJmsmiCqWqOvhejt+c5GZgP2BqsD/JscCLgRdV1ahpvjXAGoBerzd0G0nS3Ez89F2SPacXNiTZA9gbuHLIdocD7wBeWlW3LGiRkiRggkIpyVFJtgAHAmuTrGu7ngNcluRS4GzgjVV1QzvmtL7l4x8GHgyc3y4fP2VhX4EkKSNmqbQVer1eTU3dYxZQkjSDJOurqjesb2KOlCRJi5+hJEnqDENJktQZhpIkqTMMJUlSZxhKkqTOMJQkSZ1hKEmSOsNQkiR1hqEkSeoMQ0mS1BmGkiSpMwwlSVJnGEqSpM4wlCRJnWEoSZI6w1CSJHWGoSRJ6gxDSZLUGYaSJKkzDCVJUmcYSpKkzpiYUEpyTJJNSe5M0utrPyDJpe3tsiRHzbKfE5JUkp3mv2pJUr9l4y5gO7ocOBo4dUh7r6puT/Jo4LIk/7eqbh/cQZLdgUOAq+a9WknSPUzMkVJVba6qK4a039IXQDsCNcNuPgi8fZZtJEnzZGJCaSZJnpFkE7AROG7EUdJLgaur6rIFL1CSBCyy6bskFwC7DOlaXVXnjhpXVRcD+ybZBzg9yXlV9eu+/T4AWA0cuhU1rAJWAaxYsWIbX4EkaSaLKpSq6uB7OX5zkpuB/YCpvq7HAXvSnG8C2A24JMkBVfWzgX2sAdYA9Ho9p/kkaTtaVKE0F0n2BH7aLnTYA9gbuLJ/m6raCOzcN+ZKmsURNyxgqZK05E3MOaUkRyXZAhwIrE2yru16Ds0R0KXA2cAbp8MmyWn9y8clSeOVKmeg5qrX69XU1NTsG0qSfifJ+qoaekAwMUdKkqTFz1CSJHWG03f3QpLrgZ/McfhOgAsp7uL7cXe+H3fxvbi7SXg/9qiqRw7rMJTGJMnUqDnVpcj34+58P+7ie3F3k/5+OH0nSeoMQ0mS1BmG0visGXcBHeP7cXe+H3fxvbi7iX4/PKckSeoMj5QkSZ1hKEmSOsNQGoMkhye5IskPkrxz3PWMU5Ldk/xLks3t5ezfPO6axi3JDkk2JPn8uGsZtyTLk5yR5Lvt35EDx13TuCR5a/tv5PIk/yfJjuOuaT4YSgssyQ7A/wL+HfBE4JVJnjjeqsbqduBvqmof4JnAXy7x9wPgzcDmcRfRER8CvlhVvw/8AUv0fUmyK/AmmqsX7AfsAPzxeKuaH4bSwjsA+EFV/aiqfgN8BjhizDWNTVVdW1WXtD//kuaXzq7jrWp8kuwG/HvgtHHXMm5JHgI8D/goQFX9pqpuHGtR47UMuH+SZcADgGvGXM+8MJQW3q7AT/seb2EJ/xLul2QlsD9w8ZhLGae/A94O3DnmOrrgscD1wD+005mnJXnguIsah6q6Gvhb4CrgWuCmqvrSeKuaH4bSwsuQtiW/Lj/Jg4AzgbdU1S/GXc84JHkxcF1VrR93LR2xDHgq8L+ran/gZmBJnoNN8jCaGZU9gccAD0zyH8Zb1fwwlBbeFmD3vse7MaGH4VsryX1pAunTVXXWuOsZo2cDL22vfPwZ4KAknxpvSWO1BdhSVdNHzmfQhNRSdDDw46q6vqp+C5wFPGvMNc0LQ2nhfQt4fJI9k/wezcnKfx5zTWOTJDTnDDZX1X8fdz3jVFUnVtVuVbWS5u/FV6pqIv83vDWq6mfAT5Ps3Ta9CPjOGEsap6uAZyZ5QPtv5kVM6KKPZeMuYKmpqtuTHA+so1lB87Gq2jTmssbp2cCfABvbS9YDvKuqvjC+ktQhfwV8uv0P3I+APx1zPWNRVRcnOQO4hGbF6gYm9OuG/JohSVJnOH0nSeoMQ0mS1BmGkiSpMwwlSVJnGEqSpM4wlDTxkqxMUkk+Pu5atrckFyZxCa0mhqGkJSvJx9uwWjnuWkZZDDUOmg7KJHckedKIbaZf18ELXZ+6zQ/Paim4GtgHuGnchcyD19B8Y3QX3Qc4GTh83IVo8TCUNPHa7wr77rjrmA9VddW4a5jBD4DDkhxSVeePuxgtDk7faeINO6fUnoc5tn3447a/2i9D7R/78CTvba96emuSm5J8OcmhQ57nte0+XtteXfjCdvvq2+bIJJ9K8r0kNyf5VZL1Sd6U5D4D+5u1xlHnlJLcJ8lxSb7VPsfN7c9vGHye6edq97VTkjVJrk1yW3ul07l+tc+7aL4B/+RhzzlKkqclOTPJdW0NP0ny90kePcc6tIh4pKSl6j3AkTRXM/0QcGPbPn1Pkj2AC4GVwNeALwIPBF4MfDHJ66vqI0P2/XKaKavzgFPa8dPeR3OtpItpphUfChzU1vB0mu8B3OoaZ/BJ4FU01+46jSYcjgL+HngO8OohY5YDXwd+Q/ON3Du2r+VjSe6sqtO34nn7bQA+RfOajgX+YbYB7eU7zqS5xMsZwE+ApwFvAI5I8uyqunIb69BiUlXevE30jSYUCvj4QPvH2/aVI8ZdSBMgfzzQvhy4FLgVeFRf+2vb/d0JHD5in48b0nYf4PR27DPmUGMNtL2yHXMJ8KC+9gcCU23fqwbGVHs7Ddihr/2JNF8A+p1teL8vbPe1F81lWm6luQzFA4a8roP72h4E3ADcATx3YJ/vaLf/0rj/Pnmb35vTd9IQSf4AeD5wZlV9pr+vmkty/0eaI4mXDRl+blV9cdh+q+qHQ9rupDkSAjjsXpQ97XXt/Tur6ld9z3MzzS93gD8fMu4W4K+r6o6+Md+hOXraJ8mDt7WQqvopzdV0dwX+epbNjwAeAXy2qr420PffgCuBQ5Ks2NY6tHg4fScNd2B7/9Ak7x7S/8j2fp8hfd8ctdMkjwDeBvwhzeW+By/vveu2lTnUU2mO1i4c0ncRzZHI/kP6vl/Dr/r70/Z+OfDLOdTzXuDPgLcn+UhV/euI7aYv4PeVwY5qLvnyVZqj3v1pri+kCWQoScM9or0/pL2N8qAhbT8btmGS5TQXedyTJrg+AfwbzfTYcuDNwP3mVO3dPRT4t6r6zWBH+8v9BmDnIeNuHLG/29v7HeZSTFX9Isl7gA8D76Y5PzTMQ9v7a0f0T7cvn0sdWhwMJWm46c80vbmq/sc2jh31DQt/ThNI76mqd/d3JDmQJpS2h5uAhye5bzXL4fufZxmwEzDsiGg+nQq8CfiLJKPez+n3fJcR/Y8e2E4TyHNKWsqmz50MOwL4Rnv/3O34fHu192cO6Xv+iDEz1TjKBpp/288b0ve8dl+XbMP+7rWqup3mfNYOwAdGbLahvX/BYEcbps9pHy5o7VpYhpKWsp+39/c4cV5VUzTLwI9O8rrBfoAkT0oybBpslCvb+xcM7Gd/4MRtrXEGH2vv35vkd9/20P78vvbhR7dhf9tFVZ1D856+GHj2kE3OoZnOfGWSZw70vYXmHNwF1feB4SQrkvx+/+vU4ub0nZayL9MsOvhIkjOAXwE3VtWH2/5X0Zx0/2iSN9F8tuhGYDfgycB+NAsirtvK5/tE+3x/l+SFwPeBx9P8kj4LeMUcaryHqvrHJEcAfwRsSnIOzZTikTTTh/9UVZ/eypq3txNojkL3Guyoql+1/wH4HHBRks/RLGh4GnAozbm61w8M+wTNUeYLGb6wQ4uMR0pasqpqHfA3wG+BtwL/ieaX5nT/FppfiKtpptFeTXNe5Fk0vyxfD2zchue7hmY6cC3NVNTxwB7AG4F3zqXGGbwS+EuaI63XA8cB/699zldubc3bW1V9E/jsDP3n0hxFfYFmefwJNCscTwGeVlU/Wog6NT6p8lvvJUnd4JGSJKkzDCVJUmcYSpKkzjCUJEmdYShJkjrDUJIkdYahJEnqDENJktQZhpIkqTP+P7g9xe71yOY4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the acceptance rate= 0.4544166666666667\n",
      "VMC finished, elapsed time = 1.2084949016571045 sec\n"
     ]
    }
   ],
   "source": [
    "\n",
    "scratch=\"output_files/\"\n",
    "\n",
    "# system size\n",
    "N=40\n",
    "# number of electrons\n",
    "N_pt=int(N/2)\n",
    "\n",
    "\n",
    "t=1.0 # hopping amplitude\n",
    "BC=np.exp(1j*pi) # boundary condition on a chain, you can put BC=0 for open chain\n",
    "# BC=1 periodic boundary condition and BC=-1 is anti-periodic\n",
    "# do not put BC=1 since the gs is not unique in that case\n",
    "numconfig=100\n",
    "\n",
    "t_timer=time.time()\n",
    "random.seed(time.time())\n",
    "# #################################################################\n",
    "# ##### Execution Part\n",
    "# #################################################################\n",
    "\n",
    "# ##### initialization\n",
    "n_occ1=initialize_wf(N,N_pt) # occupation numbers \n",
    "n_pos1,=n_occ1.nonzero() # position of the particles\n",
    "# #################################################################\n",
    "# ##### Slater determinant\n",
    "# #################################################################\n",
    "V1=wf_gen(N,N_pt,BC,t) # occupied single particle states\n",
    "wf_ref=np.transpose(V1[n_pos1,:]) # initial slater determinant\n",
    "\n",
    "\n",
    "##### make sure the initial config of slater determinant is not singular\n",
    "counter=0\n",
    "epsilon= 1e-8\n",
    "while np.linalg.cond(wf_ref) > 1/epsilon:\n",
    "    counter+=1\n",
    "    assert counter<=40, \"wave function cannot be constructed!\"\n",
    "    n_occ1=initialize_wf(N,N_pt) \n",
    "    n_pos1,=n_occ1.nonzero() \n",
    "    wf_ref=np.transpose(V1[n_pos1,:]) # initial slater determinant\n",
    "\n",
    "\n",
    "\n",
    "VMC_main(numconfig,t,V1,wf_ref,N,N_pt,BC,n_occ1,n_pos1)\n",
    "\n",
    "elapsed = time.time() - t_timer\n",
    "print(\"VMC finished, elapsed time =\", elapsed, \"sec\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}