VMC-Renyi_AS.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": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generates single particle states\n",
    "def wf_gen(N,N_pt,BC,t1,t2):\n",
    "    h1=np.ones(N-1)\n",
    "    h1[0::2]=0\n",
    "    h2=np.ones(N-1)\n",
    "    h2[1::2]=0\n",
    "    hop= np.diag(t1*h1+t2*h2+0j,1)\n",
    "    hop[N-1,0]= t1*BC\n",
    "    H_t= -(hop+ np.matrix(hop).H)/2 \n",
    "    energies, evecs= np.linalg.eigh(H_t)\n",
    "    return evecs[:,:N_pt]\n",
    "\n",
    "\n",
    "def exact_renyi_calc(r,GA,epsilon=1e-9):\n",
    "    chi0, _ =np.linalg.eigh(GA)\n",
    "    chi1=chi0[np.nonzero(np.abs(chi0)>epsilon)]\n",
    "    chi2=chi1[np.nonzero(np.abs(chi1-1)>epsilon)]\n",
    "#     return -np.sum((1-chi2)*np.log(1-chi2)+chi2*np.log(chi2))\n",
    "    return np.sum(np.log((1-chi2)**r+chi2**r))/(1-r)\n",
    "\n",
    "\n",
    "# VMC function \n",
    "def VMC_equal_number(numconfig,V1,wf_1,wf_1_anc,N,inds_A,N_pt,\\\n",
    "             n_occ1,coords1,\\\n",
    "             n_occ1_anc,coords1_anc):\n",
    "\n",
    "    move_attempted = 0\n",
    "    move_accepted = 0\n",
    "\n",
    "    inv_wf=np.linalg.inv(wf_1)\n",
    "    inv_wf_anc=np.linalg.inv(wf_1_anc)\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",
    "    fraction=np.zeros(len(inds_A)+1,dtype=np.float64) \n",
    "    \n",
    "    # possible steps\n",
    "    step_abs=np.arange(1,Ns+1)\n",
    "    step_vals=np.sort(np.concatenate((-step_abs,step_abs),axis=0)).tolist()\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 Ns steps left or right\n",
    "            stepx=random.sample(step_vals,1)[0]\n",
    "            ptcls_x= np.mod( coords1[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(inv_wf[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",
    "                u_1=np.reshape(pt_wf_1,(N_pt,1))  - np.reshape(wf_1[:,moved_elec],(N_pt,1))\n",
    "                v=np.zeros((N_pt,1))\n",
    "                v[moved_elec]=1\n",
    "                inv_wf=inv_wf - np.dot(np.dot(np.dot(inv_wf,u_1),v.T),inv_wf) \\\n",
    "                                 /(1+np.dot(v.T,np.dot(inv_wf,u_1))) \n",
    "                wf_1[:,moved_elec]=np.reshape(pt_wf_1,(N_pt,))\n",
    "#                 wf_inv=np.linalg.inv(wf_1)\n",
    "    \n",
    "             \n",
    "                move_accepted=move_accepted+1\n",
    "                n_occ1[coords1[moved_elec]]= n_occ1[coords1[moved_elec]]-1\n",
    "\n",
    "                coords1[moved_elec]=ptcls_x\n",
    "                n_occ1[ptcls_x]= n_occ1[ptcls_x]+1\n",
    "                \n",
    "        assert np.sum(n_occ1)== N_pt, 'n_occ anc ptcle is %d' % (np.sum(n_occ1))\n",
    "\n",
    "####################### ancilla ############################\n",
    "\n",
    "        for moved_elec in range(N_pt):\n",
    "            move_attempted=move_attempted+1\n",
    "\n",
    "            # random walk of Ns steps left or right\n",
    "            stepx=random.sample(step_vals,1)[0]\n",
    "            ptcls_x= np.mod( coords1_anc[moved_elec]+stepx, N) # new configuration\n",
    "\n",
    "            if n_occ1_anc[ptcls_x]==1:\n",
    "                continue\n",
    "                            \n",
    "            pt_wf_1=np.transpose(V1[ptcls_x,:])\n",
    "\n",
    "            rel=np.dot(inv_wf_anc[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",
    "                u_1=np.reshape(pt_wf_1,(N_pt,1))  - np.reshape(wf_1_anc[:,moved_elec],(N_pt,1))\n",
    "                v=np.zeros((N_pt,1))\n",
    "                v[moved_elec]=1\n",
    "                inv_wf_anc=inv_wf_anc - np.dot(np.dot(np.dot(inv_wf_anc,u_1),v.T),inv_wf_anc) \\\n",
    "                         /(1+np.dot(v.T,np.dot(inv_wf_anc,u_1))) \n",
    "                wf_1_anc[:,moved_elec]=np.reshape(pt_wf_1,(N_pt,))\n",
    "#                 wf_inv_anc=np.linalg.inv(wf_1_anc)\n",
    "\n",
    "                \n",
    "                move_accepted=move_accepted+1\n",
    "                n_occ1_anc[coords1_anc[moved_elec]]= n_occ1_anc[coords1_anc[moved_elec]]-1\n",
    "            \n",
    "                coords1_anc[moved_elec]=ptcls_x\n",
    "                n_occ1_anc[ptcls_x]= n_occ1_anc[ptcls_x]+1\n",
    "            \n",
    "\n",
    "        assert np.sum(n_occ1_anc)== N_pt, 'n_occ anc ptcle is %d' % (np.sum(n_occ1_anc))\n",
    "\n",
    "\n",
    "# ##############################################################\n",
    "            \n",
    "        if step> (min_step-1):\n",
    "            Na = np.sum(n_occ1[inds_A])\n",
    "            Na_anc = np.sum(n_occ1_anc[inds_A])\n",
    "            if Na==Na_anc:\n",
    "                fraction[int(Na)] += 1\n",
    "\n",
    "        if (step%500) ==0:\n",
    "            inv_wf=np.linalg.inv(wf_1)\n",
    "            inv_wf_anc=np.linalg.inv(wf_1_anc)\n",
    "\n",
    "\n",
    "    acc_ratio=move_accepted/move_attempted\n",
    "    fraction = fraction/numconfig\n",
    "    \n",
    "#     print(\"fractions are \", fraction)\n",
    "    print(\"fraction acceptance rate=\", acc_ratio)\n",
    "    return fraction\n",
    "\n",
    "\n",
    "# VMC function \n",
    "def VMC_amplitude_ratio(numconfig,V1,wf_1,wf_1_anc,N,inds_A,N_pt,N_pt_A,\\\n",
    "             n_occ1,n_pos1,coords1,\\\n",
    "             n_occ1_anc,n_pos1_anc,coords1_anc):\n",
    "\n",
    "    move_attempted = 0\n",
    "    move_accepted = 0\n",
    "\n",
    "    inv_wf=np.linalg.inv(wf_1)\n",
    "    inv_wf_anc=np.linalg.inv(wf_1_anc)\n",
    "\n",
    "    inside_A = 0\n",
    "    \n",
    "    pt_num_inside = np.argwhere( n_occ1[inds_A]>0 )\n",
    "    pt_num_inside = np.reshape( pt_num_inside, (N_pt_A,)).tolist()\n",
    "    wf_inds=np.sort( n_pos1[ inds_A[pt_num_inside] ] )-1\n",
    "\n",
    "    pt_num_inside_anc = np.argwhere( n_occ1_anc[inds_A]>0 )\n",
    "    pt_num_inside_anc = np.reshape( pt_num_inside_anc, (N_pt_A,)).tolist()\n",
    "    wf_inds_anc = np.sort( n_pos1_anc[inds_A[pt_num_inside_anc] ] )-1\n",
    "\n",
    "    # swapping between original wf and ancilla wf\n",
    "    wf_1_swap = np.copy(wf_1)\n",
    "    wf_1_swap[:,wf_inds] = wf_1_anc[:,wf_inds_anc]\n",
    "    wf_1_swap_anc = np.copy(wf_1_anc)  \n",
    "    wf_1_swap_anc[:,wf_inds_anc] = wf_1[:,wf_inds] \n",
    "        \n",
    "    count=0 # counter for energy\n",
    "    howoften=10 # calculate SWAP amplitude every 10 steps\n",
    "    min_step=500\n",
    "    \n",
    "    ent_ratio=np.zeros(int(numconfig/howoften),dtype=np.float64) \n",
    "    \n",
    "    # possible steps\n",
    "    step_abs=np.arange(1,Ns+1)\n",
    "    step_vals=np.sort(np.concatenate((-step_abs,step_abs),axis=0)).tolist()\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 Ns steps left or right\n",
    "            stepx=random.sample(step_vals,1)[0]\n",
    "            ptcls_x= np.mod( coords1[moved_elec]+stepx, N) # new configuration\n",
    "\n",
    "            if n_occ1[ptcls_x]==1:\n",
    "                continue\n",
    "\n",
    "            val_orig = np.min( np.abs(inds_A- coords1[moved_elec]) )\n",
    "            val_dest = np.min( np.abs(inds_A-ptcls_x) )\n",
    "            if val_orig==0 and val_dest==0 :\n",
    "                inside_A=1\n",
    "                val_1 = np.min(np.abs(wf_inds-moved_elec))\n",
    "                ind_1 = np.argmin(np.abs(wf_inds-moved_elec))\n",
    "                assert val_1==0 , 'wrong inds_A for wf'\n",
    "            elif val_orig!=0 and val_dest!=0 :\n",
    "                inside_A=0\n",
    "            else:\n",
    "                continue\n",
    "            \n",
    "            pt_wf_1=np.transpose(V1[ptcls_x,:])\n",
    "\n",
    "            rel=np.dot(inv_wf[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",
    "                u_1=np.reshape(pt_wf_1,(N_pt,1))  - np.reshape(wf_1[:,moved_elec],(N_pt,1))\n",
    "                v=np.zeros((N_pt,1))\n",
    "                v[moved_elec]=1\n",
    "                inv_wf=inv_wf - np.dot(np.dot(np.dot(inv_wf,u_1),v.T),inv_wf) \\\n",
    "                                 /(1+np.dot(v.T,np.dot(inv_wf,u_1))) \n",
    "                wf_1[:,moved_elec]=np.reshape(pt_wf_1,(N_pt,))\n",
    "#                 wf_inv=np.linalg.inv(wf_1)\n",
    "    \n",
    "                if inside_A==1:\n",
    "#                     print(wf_1_swap_anc[:,wf_inds_anc[ind_1]].shape)\n",
    "#                     print(pt_wf_1.shape)\n",
    "                    wf_1_swap_anc[:,wf_inds_anc[ind_1]]= np.reshape(pt_wf_1,(N_pt,))\n",
    "                    inside_A=0\n",
    "                else:\n",
    "#                     print(wf_1_swap[:,moved_elec].shape)\n",
    "#                     print(pt_wf_1.shape)\n",
    "                    wf_1_swap[:,moved_elec]= np.reshape(pt_wf_1,(N_pt,))\n",
    "                \n",
    "                move_accepted=move_accepted+1\n",
    "                n_occ1[coords1[moved_elec]]= n_occ1[coords1[moved_elec]]-1\n",
    "                n_pos1[coords1[moved_elec]]= 0\n",
    "\n",
    "                coords1[moved_elec]=ptcls_x\n",
    "                n_occ1[ptcls_x]= n_occ1[ptcls_x]+1\n",
    "                n_pos1[ptcls_x]= moved_elec+1\n",
    "                \n",
    "        x=np.argwhere(n_pos1>0)\n",
    "        assert len(x)== N_pt, 'no of ptcle is %d' % (len(x))\n",
    "        assert np.sum(n_occ1[inds_A])== N_pt_A, 'n_occ_A anc ptcle is %d' % (np.sum(n_occ1[inds_A]))\n",
    "        assert np.sum(n_occ1)== N_pt, 'n_occ anc ptcle is %d' % (np.sum(n_occ1))\n",
    "\n",
    "####################### ancilla ############################\n",
    "\n",
    "        for moved_elec in range(N_pt):\n",
    "            move_attempted=move_attempted+1\n",
    "\n",
    "             # random walk of Ns steps left or right\n",
    "            stepx=random.sample(step_vals,1)[0]\n",
    "            ptcls_x= np.mod( coords1_anc[moved_elec]+stepx, N) # new configuration\n",
    "\n",
    "            if n_occ1_anc[ptcls_x]==1:\n",
    "                continue\n",
    "            \n",
    "            val_orig = np.min( np.abs(inds_A- coords1_anc[moved_elec]) )\n",
    "            val_dest = np.min( np.abs(inds_A-ptcls_x) )\n",
    "            if val_orig==0 and val_dest==0 :\n",
    "                inside_A = 1\n",
    "                val_1 = np.min(np.abs(wf_inds_anc-moved_elec))\n",
    "                ind_1 = np.argmin(np.abs(wf_inds_anc-moved_elec))\n",
    "#                 print(ind_1,val_1)\n",
    "                assert val_1==0 , 'wrong inds_A for anc wf'\n",
    "            elif val_orig!=0 and val_dest!=0 :\n",
    "                inside_A=0\n",
    "            else:\n",
    "                continue\n",
    "                \n",
    "            pt_wf_1=np.transpose(V1[ptcls_x,:])\n",
    "\n",
    "            rel=np.dot(inv_wf_anc[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",
    "                u_1=np.reshape(pt_wf_1,(N_pt,1))  - np.reshape(wf_1_anc[:,moved_elec],(N_pt,1))\n",
    "                v=np.zeros((N_pt,1))\n",
    "                v[moved_elec]=1\n",
    "                inv_wf_anc=inv_wf_anc - np.dot(np.dot(np.dot(inv_wf_anc,u_1),v.T),inv_wf_anc) \\\n",
    "                         /(1+np.dot(v.T,np.dot(inv_wf_anc,u_1))) \n",
    "                wf_1_anc[:,moved_elec]=np.reshape(pt_wf_1,(N_pt,))\n",
    "#                 wf_inv_anc=np.linalg.inv(wf_1_anc)\n",
    "\n",
    "                if inside_A==1:\n",
    "                    wf_1_swap[:,wf_inds[ind_1]]= np.reshape(pt_wf_1,(N_pt,))\n",
    "                    inside_A=0\n",
    "                else:\n",
    "                    wf_1_swap_anc[:,moved_elec]= np.reshape(pt_wf_1,(N_pt,))\n",
    "                \n",
    "                move_accepted=move_accepted+1\n",
    "                n_occ1_anc[coords1_anc[moved_elec]]= n_occ1_anc[coords1_anc[moved_elec]]-1\n",
    "                n_pos1_anc[coords1_anc[moved_elec]]= 0\n",
    "\n",
    "                coords1_anc[moved_elec]=ptcls_x\n",
    "                n_occ1_anc[ptcls_x]= n_occ1_anc[ptcls_x]+1\n",
    "                n_pos1_anc[ptcls_x]= moved_elec+1\n",
    "\n",
    "\n",
    "        x=np.argwhere(n_pos1_anc>0)\n",
    "        assert len(x)== N_pt, 'no of anc ptcle is %d' % (len(x))\n",
    "        assert np.sum(n_occ1_anc[inds_A])== N_pt_A, 'n_occ_A anc ptcle is %d' % (np.sum(n_occ1_anc[inds_A]))\n",
    "        assert np.sum(n_occ1_anc)== N_pt, 'n_occ anc ptcle is %d' % (np.sum(n_occ1_anc))\n",
    "#         print(n_pos1)\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",
    "                ent_ratio[count] = np.abs( np.linalg.det(wf_1_swap)*np.linalg.det(wf_1_swap_anc)/\\\n",
    "                                        (np.linalg.det(wf_1)*np.linalg.det(wf_1_anc)) )\n",
    "                count+=1\n",
    "        \n",
    "        if (step%500) ==0:\n",
    "            inv_wf=np.linalg.inv(wf_1)\n",
    "            inv_wf_anc=np.linalg.inv(wf_1_anc)\n",
    "\n",
    "\n",
    "    acc_ratio=move_accepted/move_attempted\n",
    "    \n",
    "    \n",
    "    print(\"Amplitude acceptance rate=\", acc_ratio)\n",
    "    return np.mean(ent_ratio)\n",
    "\n",
    "\n",
    "# VMC function \n",
    "def VMC_phase_ratio(numconfig,V1,wf_1,wf_1_anc,N,inds_A,N_pt,N_pt_A,\\\n",
    "             n_occ1,n_pos1,coords1,\\\n",
    "             n_occ1_anc,n_pos1_anc,coords1_anc):\n",
    "\n",
    "    move_attempted = 0\n",
    "    move_accepted = 0\n",
    "\n",
    "    inv_wf=np.linalg.inv(wf_1)\n",
    "    inv_wf_anc=np.linalg.inv(wf_1_anc)\n",
    "\n",
    "    inside_A = 0\n",
    "    \n",
    "    pt_num_inside = np.argwhere( n_occ1[inds_A]>0 )\n",
    "    pt_num_inside = np.reshape( pt_num_inside, (N_pt_A,)).tolist()\n",
    "    wf_inds=np.sort( n_pos1[ inds_A[pt_num_inside] ] )-1\n",
    "\n",
    "    pt_num_inside_anc = np.argwhere( n_occ1_anc[inds_A]>0 )\n",
    "    pt_num_inside_anc = np.reshape( pt_num_inside_anc, (N_pt_A,)).tolist()\n",
    "    wf_inds_anc = np.sort( n_pos1_anc[inds_A[pt_num_inside_anc] ] )-1\n",
    "\n",
    "    # swapping between original wf and ancilla wf\n",
    "    wf_1_swap = np.copy(wf_1)\n",
    "    wf_1_swap[:,wf_inds] = wf_1_anc[:,wf_inds_anc]\n",
    "    wf_1_swap_anc = np.copy(wf_1_anc)  \n",
    "    wf_1_swap_anc[:,wf_inds_anc] = wf_1[:,wf_inds] \n",
    "    \n",
    "    inv_wf_1_swap = np.linalg.inv(wf_1_swap)\n",
    "    inv_wf_1_swap_anc = np.linalg.inv(wf_1_swap_anc)\n",
    "\n",
    "    \n",
    "    min_step=500\n",
    "    \n",
    "    phi_0=np.angle( np.conj(np.linalg.det(wf_1)*np.linalg.det(wf_1_anc))*\\\n",
    "                   np.linalg.det(wf_1_swap)*np.linalg.det(wf_1_swap_anc) )    \n",
    "    phi=phi_0\n",
    "#     phi=np.exp(1j*phi_0)\n",
    "\n",
    "    ent_ratio=np.zeros(numconfig,dtype=np.complex64) \n",
    "    \n",
    "    # possible steps\n",
    "    step_abs=np.arange(1,Ns+1)\n",
    "    step_vals=np.sort(np.concatenate((-step_abs,step_abs),axis=0)).tolist()\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 Ns steps left or right\n",
    "            stepx=random.sample(step_vals,1)[0]\n",
    "            ptcls_x= np.mod( coords1[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",
    "            val_orig = np.min( np.abs(inds_A- coords1[moved_elec]) )\n",
    "            val_dest = np.min( np.abs(inds_A-ptcls_x) )\n",
    "            if val_orig==0 and val_dest==0 :\n",
    "                inside_A=1\n",
    "                val_1 = np.min(np.abs(wf_inds-moved_elec))\n",
    "                ind_1 = np.argmin(np.abs(wf_inds-moved_elec))\n",
    "                assert val_1==0 , 'wrong inds_A for wf'\n",
    "                rel = np.dot(inv_wf_1_swap_anc[wf_inds_anc[ind_1],:],pt_wf_1)\n",
    "            elif val_orig!=0 and val_dest!=0 :\n",
    "                inside_A=0\n",
    "                rel = np.dot(inv_wf_1_swap[moved_elec,:],pt_wf_1)\n",
    "            else:\n",
    "                continue\n",
    "            \n",
    "\n",
    "            rel=rel*np.conj(np.dot(inv_wf[moved_elec,:],pt_wf_1))\n",
    "\n",
    "            alpha=min(1, np.abs(rel))\n",
    "            \n",
    "            random_num=random.random()\n",
    "\n",
    "            if random_num <= alpha:\n",
    "                u_1=np.reshape(pt_wf_1,(N_pt,1))  - np.reshape(wf_1[:,moved_elec],(N_pt,1))\n",
    "                v=np.zeros((N_pt,1))\n",
    "                v[moved_elec]=1\n",
    "                inv_wf=inv_wf - np.dot(np.dot(np.dot(inv_wf,u_1),v.T),inv_wf) \\\n",
    "                                 /(1+np.dot(v.T,np.dot(inv_wf,u_1))) \n",
    "                wf_1[:,moved_elec]=np.reshape(pt_wf_1,(N_pt,))\n",
    "#                 wf_inv=np.linalg.inv(wf_1)\n",
    "\n",
    "                if inside_A==1:\n",
    "                    u_1=np.reshape(pt_wf_1,(N_pt,1)) - np.reshape(wf_1_swap_anc[:,wf_inds_anc[ind_1]],(N_pt,1))\n",
    "                    v=np.zeros((N_pt,1))\n",
    "                    v[wf_inds_anc[ind_1]]=1\n",
    "                    inv_wf_1_swap_anc=inv_wf_1_swap_anc - np.dot(np.dot(np.dot(inv_wf_1_swap_anc,u_1),v.T),inv_wf_1_swap_anc) \\\n",
    "                                 /(1+np.dot(v.T,np.dot(inv_wf_1_swap_anc,u_1))) \n",
    "                    wf_1_swap_anc[:,wf_inds_anc[ind_1]]= np.reshape(pt_wf_1,(N_pt,))\n",
    "                    inside_A=0\n",
    "                else:\n",
    "                    u_1=np.reshape(pt_wf_1,(N_pt,1)) - np.reshape(wf_1_swap[:,moved_elec],(N_pt,1))\n",
    "                    inv_wf_1_swap=inv_wf_1_swap - np.dot(np.dot(np.dot(inv_wf_1_swap,u_1),v.T),inv_wf_1_swap) \\\n",
    "                                 /(1+np.dot(v.T,np.dot(inv_wf_1_swap,u_1))) \n",
    "\n",
    "                    wf_1_swap[:,moved_elec]= np.reshape(pt_wf_1,(N_pt,))\n",
    "                \n",
    "#                 phi=phi*np.exp(1j*np.angle(rel))\n",
    "                phi=phi+np.angle(rel)\n",
    "                \n",
    "                move_accepted=move_accepted+1\n",
    "                n_occ1[coords1[moved_elec]]= n_occ1[coords1[moved_elec]]-1\n",
    "                n_pos1[coords1[moved_elec]]= 0\n",
    "\n",
    "                coords1[moved_elec]=ptcls_x\n",
    "                n_occ1[ptcls_x]= n_occ1[ptcls_x]+1\n",
    "                n_pos1[ptcls_x]= moved_elec+1\n",
    "                \n",
    "        x=np.argwhere(n_pos1>0)\n",
    "        assert len(x)== N_pt, 'no of ptcle is %d' % (len(x))\n",
    "        assert np.sum(n_occ1)== N_pt, 'n_occ anc ptcle is %d' % (np.sum(n_occ1))\n",
    "\n",
    "####################### ancilla ############################\n",
    "\n",
    "        for moved_elec in range(N_pt):\n",
    "            move_attempted=move_attempted+1\n",
    "\n",
    "             # random walk of Ns steps left or right\n",
    "            stepx=random.sample(step_vals,1)[0]\n",
    "            ptcls_x= np.mod( coords1_anc[moved_elec]+stepx, N) # new configuration\n",
    "\n",
    "            if n_occ1_anc[ptcls_x]==1:\n",
    "                continue\n",
    "            \n",
    "            pt_wf_1=np.transpose(V1[ptcls_x,:])\n",
    "\n",
    "            val_orig = np.min( np.abs(inds_A- coords1_anc[moved_elec]) )\n",
    "            val_dest = np.min( np.abs(inds_A-ptcls_x) )\n",
    "            if val_orig==0 and val_dest==0 :\n",
    "                inside_A = 1\n",
    "                val_1 = np.min(np.abs(wf_inds_anc-moved_elec))\n",
    "                ind_1 = np.argmin(np.abs(wf_inds_anc-moved_elec))\n",
    "#                 print(ind_1,val_1)\n",
    "                assert val_1==0 , 'wrong inds_A for anc wf'\n",
    "                rel = np.dot(inv_wf_1_swap[wf_inds[ind_1],:],pt_wf_1)\n",
    "            elif val_orig!=0 and val_dest!=0 :\n",
    "                inside_A=0\n",
    "                rel = np.dot(inv_wf_1_swap_anc[moved_elec,:],pt_wf_1)\n",
    "            else:\n",
    "                continue\n",
    "                \n",
    "\n",
    "            rel=rel*np.conj(np.dot(inv_wf_anc[moved_elec,:],pt_wf_1))\n",
    "\n",
    "            alpha=min(1, np.abs(rel))\n",
    "            \n",
    "            random_num=random.random()\n",
    "\n",
    "            if random_num <= alpha:\n",
    "                u_1=np.reshape(pt_wf_1,(N_pt,1))  - np.reshape(wf_1_anc[:,moved_elec],(N_pt,1))\n",
    "                v=np.zeros((N_pt,1))\n",
    "                v[moved_elec]=1\n",
    "                inv_wf_anc=inv_wf_anc - np.dot(np.dot(np.dot(inv_wf_anc,u_1),v.T),inv_wf_anc) \\\n",
    "                                 /(1+np.dot(v.T,np.dot(inv_wf_anc,u_1))) \n",
    "                wf_1_anc[:,moved_elec]=np.reshape(pt_wf_1,(N_pt,))\n",
    "#                 wf_inv_anc=np.linalg.inv(wf_1_anc)\n",
    "\n",
    "                if inside_A==1:\n",
    "                    u_1=np.reshape(pt_wf_1,(N_pt,1)) - np.reshape(wf_1_swap[:,wf_inds_anc[ind_1]],(N_pt,1))\n",
    "                    v=np.zeros((N_pt,1))\n",
    "                    v[wf_inds[ind_1]]=1\n",
    "                    inv_wf_1_swap=inv_wf_1_swap - np.dot(np.dot(np.dot(inv_wf_1_swap,u_1),v.T),inv_wf_1_swap) \\\n",
    "                                /(1+np.dot(v.T,np.dot(inv_wf_1_swap,u_1))) \n",
    "\n",
    "                    wf_1_swap[:,wf_inds[ind_1]]= np.reshape(pt_wf_1,(N_pt,))\n",
    "                    inside_A=0\n",
    "\n",
    "                else:\n",
    "                    u_1=np.reshape(pt_wf_1,(N_pt,1)) - np.reshape(wf_1_swap_anc[:,moved_elec],(N_pt,1))\n",
    "                    inv_wf_1_swap_anc=inv_wf_1_swap_anc - np.dot(np.dot(np.dot(inv_wf_1_swap_anc,u_1),v.T),inv_wf_1_swap_anc) \\\n",
    "                                 /(1+np.dot(v.T,np.dot(inv_wf_1_swap_anc,u_1))) \n",
    "\n",
    "                    wf_1_swap_anc[:,moved_elec]= np.reshape(pt_wf_1,(N_pt,))\n",
    "                \n",
    "#                 phi=phi*np.exp(1j*np.angle(rel))\n",
    "                phi=phi+np.angle(rel)\n",
    "\n",
    "\n",
    "                move_accepted=move_accepted+1\n",
    "                n_occ1_anc[coords1_anc[moved_elec]]= n_occ1_anc[coords1_anc[moved_elec]]-1\n",
    "                n_pos1_anc[coords1_anc[moved_elec]]= 0\n",
    "\n",
    "                coords1_anc[moved_elec]=ptcls_x\n",
    "                n_occ1_anc[ptcls_x]= n_occ1_anc[ptcls_x]+1\n",
    "                n_pos1_anc[ptcls_x]= moved_elec+1\n",
    "\n",
    "\n",
    "        x=np.argwhere(n_pos1_anc>0)\n",
    "        assert len(x)== N_pt, 'no of anc ptcle is %d' % (len(x))\n",
    "        assert np.sum(n_occ1_anc)== N_pt, 'n_occ anc ptcle is %d' % (np.sum(n_occ1_anc))\n",
    "\n",
    "\n",
    "# ##############################################################\n",
    "        assert n_occ1.sum()==N_pt, \"total particle number changes!!\"      \n",
    "            \n",
    "        if step> (min_step-1):\n",
    "            ent_ratio[step-min_step] = phi\n",
    "            \n",
    "        if (step%500) ==0:\n",
    "            inv_wf=np.linalg.inv(wf_1)\n",
    "            inv_wf_anc=np.linalg.inv(wf_1_anc)\n",
    "            inv_wf_1_swap=np.linalg.inv(wf_1_swap)\n",
    "            inv_wf_1_swap_anc=np.linalg.inv(wf_1_swap_anc)\n",
    "\n",
    "    acc_ratio=move_accepted/move_attempted\n",
    "        \n",
    "    print(\"Angle acceptance rate=\", acc_ratio)\n",
    "    return np.mean(np.exp(1j*ent_ratio))\n",
    "#     return np.mean(ent_ratio)\n",
    "\n",
    "\n",
    "def Renyi_vmc_runner(numconfig,V1,N, N_pt,Lsub):\n",
    "    \n",
    "    inds_A= np.arange(0,Lsub)\n",
    "\n",
    "    ent_amp=np.zeros(Lsub+1)\n",
    "    ent_phase=np.zeros(Lsub+1,dtype=np.complex64)\n",
    "\n",
    "    Na_list=np.arange(int(Lsub/2)-2,int(Lsub/2)+3)\n",
    "\n",
    "    for N_pt_A in Na_list:    \n",
    "        print(N_pt_A)\n",
    "        n_occ0, n_pos0, coords0=initialize_wf(N, N_pt, inds_A, N_pt_A)\n",
    "        wf_0=np.transpose(V1[coords0,:])\n",
    "        n_occ0_anc,n_pos0_anc, coords0_anc=initialize_wf(N, N_pt, inds_A, N_pt_A)\n",
    "        wf_0_anc=V1[coords0_anc,:]\n",
    "\n",
    "        epsilon= 1e-8\n",
    "        counter=0\n",
    "        while np.linalg.cond(wf_0) > 1/epsilon or np.linalg.cond(wf_0_anc) > 1/epsilon:\n",
    "            counter+=1\n",
    "            assert counter<=40, \"wave function cannot be constructed!\"\n",
    "            n_occ0, n_pos0, coords0=initialize_wf(N, N_pt, inds_A, N_pt_A)\n",
    "            wf_0=np.transpose(V1[coords0,:])\n",
    "            n_occ0_anc,n_pos0_anc, coords0_anc=initialize_wf(N, N_pt, inds_A, N_pt_A)\n",
    "            wf_0_anc=V1[coords0_anc,:]\n",
    "\n",
    "        ent_phase[N_pt_A]=VMC_phase_ratio(numconfig,V1,np.copy(wf_0),np.copy(wf_0_anc),\\\n",
    "                                          N,inds_A,N_pt,N_pt_A,\\\n",
    "                                          np.copy(n_occ0),np.copy(n_pos0),np.copy(coords0),\\\n",
    "                                          np.copy(n_occ0_anc),np.copy(n_pos0_anc),np.copy(coords0_anc))\n",
    "\n",
    "        ent_amp[N_pt_A]=VMC_amplitude_ratio(numconfig,V1,np.copy(wf_0),np.copy(wf_0_anc),\\\n",
    "                                          N,inds_A,N_pt,N_pt_A,\\\n",
    "                                          np.copy(n_occ0),np.copy(n_pos0),np.copy(coords0),\\\n",
    "                                          np.copy(n_occ0_anc),np.copy(n_pos0_anc),np.copy(coords0_anc))\n",
    "\n",
    "\n",
    "    n_occ0, coords0=initialize_wf_free(N, N_pt)\n",
    "    wf_0=np.transpose(V1[coords0,:])\n",
    "    n_occ0_anc, coords0_anc=initialize_wf_free(N, N_pt)\n",
    "    wf_0_anc=V1[coords0_anc,:]\n",
    "    ent_fr=VMC_equal_number(numconfig,V1,np.copy(wf_0),np.copy(wf_0_anc),\\\n",
    "                            N,inds_A,N_pt,\\\n",
    "                            np.copy(n_occ0),np.copy(coords0),\\\n",
    "                            np.copy(n_occ0_anc),np.copy(coords0_anc))\n",
    "\n",
    "\n",
    "\n",
    "#     ent_vmc= -np.log(np.sum(ent_fr*ent_amp*ent_phase))\n",
    "\n",
    "#     return ent_vmc\n",
    "    return ent_fr, ent_amp, ent_phase\n",
    "\n",
    "# initialize wavefunction \n",
    "def initialize_wf(N, N_pt, inds_A, N_pt_A):\n",
    "    \n",
    "    inds_A = inds_A.tolist()\n",
    "    inds_outsideA = np.arange(N)\n",
    "    inds_outsideA = np.delete(inds_outsideA, inds_A).tolist()\n",
    "    N_pt_remainder = N_pt-N_pt_A\n",
    "\n",
    "    assert N_pt_A <= len(inds_A), 'N_pt_A > subsystem A size'\n",
    "    assert N_pt_remainder <= N-len(inds_A), 'N_pt_B > subsystem B size'\n",
    "\n",
    "    coords = [0]*N_pt       \n",
    "    coords[:N_pt_A] = np.sort(random.sample(inds_A, N_pt_A))\n",
    "    coords[N_pt_A:] = np.sort(random.sample(inds_outsideA, N_pt_remainder))\n",
    "#     print(coords)\n",
    "    n_occ = np.zeros((N,))        \n",
    "    n_occ[coords] = np.ones((N_pt,))    \n",
    "    n_pos = np.zeros((N,),dtype=int)        \n",
    "    n_pos[coords] = np.arange(1,N_pt+1)\n",
    "    \n",
    "    return n_occ, n_pos, coords[:N_pt]\n",
    "\n",
    "# initialize wavefunction \n",
    "def initialize_wf_free(N,N_pt):\n",
    "    coords=np.sort(random.sample(range(N),N_pt))\n",
    "    n_occ=np.zeros((N,))        \n",
    "    n_occ[coords]=np.ones((N_pt,))        \n",
    "    n_pos = np.zeros((N,),dtype=int)        \n",
    "    n_pos[coords] = np.arange(1,N_pt+1)\n",
    "    \n",
    "    return n_occ, coords[:N_pt]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "subsystem size =  4\n",
      "0\n",
      "Angle acceptance rate= 0.19226666666666667\n",
      "Amplitude acceptance rate= 0.19103333333333333\n",
      "1\n",
      "Angle acceptance rate= 0.2414238095238095\n",
      "Amplitude acceptance rate= 0.23924761904761904\n",
      "2\n",
      "Angle acceptance rate= 0.26193809523809525\n",
      "Amplitude acceptance rate= 0.2575904761904762\n",
      "3\n",
      "Angle acceptance rate= 0.23927142857142858\n",
      "Amplitude acceptance rate= 0.23809523809523808\n",
      "4\n",
      "Angle acceptance rate= 0.1837952380952381\n",
      "Amplitude acceptance rate= 0.18453809523809525\n",
      "fraction acceptance rate= 0.30329047619047617\n",
      "VMC finished, elapsed time = 128.44628405570984 sec\n"
     ]
    }
   ],
   "source": [
    "scratch=\"output_files/\"\n",
    "\n",
    "# system size\n",
    "N=20\n",
    "# Lsub_list=np.arange(2,N-1)\n",
    "Lsub_list=[4]\n",
    "# Lsub_list=np.arange(4,10,2)\n",
    "N_pt = int(N/2)\n",
    "Ns=3 # number of steps of random walker\n",
    "\n",
    "# reference slater determinant\n",
    "t= 0.\n",
    "# hopping amplitudes\n",
    "t1= 1-t\n",
    "t2= 1+t\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",
    "V1=wf_gen(N,N_pt,BC,t1,t2) # eigenvectors in \n",
    "\n",
    "numconfig=10000\n",
    "t_timer=time.time()\n",
    "\n",
    "random.seed(time.time())\n",
    "# #################################################################\n",
    "# ##### Execution Part\n",
    "# #################################################################\n",
    "R2_ex=np.zeros(len(Lsub_list))\n",
    "R2_vmc=np.zeros(len(Lsub_list), dtype=np.complex64)\n",
    "Gmat=np.dot(V1,np.matrix(V1).H)\n",
    "\n",
    "for i_L in range(len(Lsub_list)):\n",
    "    Lsub=Lsub_list[i_L]\n",
    "    print('subsystem size = ', Lsub)\n",
    "    inds_A= np.arange(0,Lsub_list[i_L])\n",
    "\n",
    "    R2_ex[i_L]=exact_renyi_calc(2,Gmat[np.ix_(inds_A,inds_A)])\n",
    "#     R2_vmc[i_L]=Renyi_vmc_runner(numconfig,V1,N, N_pt,Lsub)\n",
    "    ent_fr, ent_amp, ent_phase=Renyi_vmc_runner(numconfig,V1,N, N_pt,Lsub)\n",
    "    R2_vmc[i_L]=-np.log(np.sum(ent_fr*ent_amp*ent_phase))\n",
    "    \n",
    "# print('VMC Renyi is', ent_vmc)\n",
    "# print('exact result is ', R2_ex)\n",
    "\n",
    "elapsed = time.time() - t_timer\n",
    "print(\"VMC finished, elapsed time =\", elapsed, \"sec\")\n",
    "\n",
    "# plt.figure()\n",
    "# plt.plot(Lsub_list,R2_ex,'o-')\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "vmc:  [0.86276156-4.732823e-08j]\n",
      "exact:  [0.87987479]\n",
      "error: [0.01944962]\n"
     ]
    }
   ],
   "source": [
    "print('vmc: ', R2_vmc)\n",
    "print('exact: ', R2_ex)\n",
    "print('error:', np.abs(R2_vmc-R2_ex)/R2_ex)\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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}