-
Notifications
You must be signed in to change notification settings - Fork 170
Expand file tree
/
Copy pathexthubbard.cc
More file actions
123 lines (106 loc) · 2.99 KB
/
exthubbard.cc
File metadata and controls
123 lines (106 loc) · 2.99 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
#include "itensor/all.h"
using namespace itensor;
int main(int argc, char* argv[])
{
//Parse the input file
if(argc != 2) { printfln("Usage: %s inputfile_exthubbard",argv[0]); return 0; }
auto input = InputGroup(argv[1],"input");
auto N = input.getInt("N");
auto Npart = input.getInt("Npart",N); //number of particles, default is N (half filling)
auto nsweeps = input.getInt("nsweeps");
auto t1 = input.getReal("t1",1);
auto t2 = input.getReal("t2",0);
auto U = input.getReal("U",0);
auto V1 = input.getReal("V1",0);
auto quiet = input.getYesNo("quiet",false);
auto table = InputGroup(input,"sweeps");
auto sweeps = Sweeps(nsweeps,table);
println(sweeps);
//
// Initialize the site degrees of freedom.
//
auto sites = Electron(N);
//
// Create the Hamiltonian using AutoMPO
//
auto ampo = AutoMPO(sites);
for(int i = 1; i <= N; ++i)
{
ampo += U,"Nupdn",i;
}
for(int b = 1; b < N; ++b)
{
ampo += -t1,"Cdagup",b,"Cup",b+1;
ampo += -t1,"Cdagup",b+1,"Cup",b;
ampo += -t1,"Cdagdn",b,"Cdn",b+1;
ampo += -t1,"Cdagdn",b+1,"Cdn",b;
ampo += V1,"Ntot",b,"Ntot",b+1;
}
for(int b = 1; b < N-1; ++b)
{
ampo += -t2,"Cdagup",b,"Cup",b+2;
ampo += -t2,"Cdagup",b+2,"Cup",b;
ampo += -t2,"Cdagdn",b,"Cdn",b+2;
ampo += -t2,"Cdagdn",b+2,"Cdn",b;
}
auto H = toMPO(ampo);
//
// Set the initial wavefunction matrix product state
// to be a Neel state.
//
auto state = InitState(sites);
int p = Npart;
for(int i = N; i >= 1; --i)
{
if(p > i)
{
println("Doubly occupying site ",i);
state.set(i,"UpDn");
p -= 2;
}
else
if(p > 0)
{
println("Singly occupying site ",i);
state.set(i,(i%2==1 ? "Up" : "Dn"));
p -= 1;
}
else
{
state.set(i,"Emp");
}
}
auto psi0 = MPS(state);
Print(totalQN(psi0));
//
// Begin the DMRG calculation
//
auto [energy,psi] = dmrg(H,psi0,sweeps,{"Quiet",quiet});
//
// Measure spin densities
//
Vector upd(N),dnd(N);
for(int j = 1; j <= N; ++j)
{
psi.position(j);
upd(j-1) = elt(dag(prime(psi(j),"Site"))*op(sites,"Nup",j)*psi(j));
dnd(j-1) = elt(dag(prime(psi(j),"Site"))*op(sites,"Ndn",j)*psi(j));
}
println("Up Density:");
for(int j = 0; j < N; ++j)
printfln("%d %.10f",1+j,upd(j));
println();
println("Dn Density:");
for(int j = 0; j < N; ++j)
printfln("%d %.10f",1+j,dnd(j));
println();
println("Total Density:");
for(int j = 0; j < N; ++j)
printfln("%d %.10f",1+j,(upd(j)+dnd(j)));
println();
//
// Print the final energy reported by DMRG
//
printfln("\nGround State Energy = %.10f",energy);
return 0;
}