-
Notifications
You must be signed in to change notification settings - Fork 0
/
testing_dipole_response.jl
149 lines (125 loc) · 4 KB
/
testing_dipole_response.jl
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
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
println("\t\"He came riding fast, like a phoenix out of fire-flames.\" -- The Dancer, PJ Harvey" )
function randH(SiteEnergy, Edisorder, Jdisorder, modelJ, N)
# Random Trace / diagonal elements
S=SiteEnergy + Edisorder*randn(N)
# Random Off-diag elements
E=modelJ(0) + Jdisorder*randn(N-1)
return (S,E)
end
function prepare_model2()
S,E=randH(5.0,Edisorder, Jdisorder, modelJ, N)
H=diagm(E,-1)+diagm(S)+diagm(E,1) #build full matrix from diagonal elements; for comparison
psi=eigvecs(H)[:,1] # 1=gnd state
# Plot_H(H)
# TODO: UPDATE THIS FOR PLOTS.jl
## Construct initial dipoles structure
density = conj(psi).*psi
return S,E,H,psi,density
end
function FieldFromDensity(density)
Field=zeros(N)
for i in 1:N
Fsum=zero(eltype(density)) # Calculate total field at each
for j in 1:N
if (j==i)
continue # avoid infinite self energies
end
Fsum+=sign(i-j)*density[j]/(i-j)^2 # How much do the dipoles respond to the electron density?
end # check power of r for field from point charge.
Field[i]=Fsum/r^2
end
return Field
end
function FieldFromDipole(dipole)
FieldDI=zeros(N)
for i in 1:N
FDIsum=zero(eltype(dipole)) # Calculate total field at each
for j in 1:N
if (j==i)
continue # avoid infinite self energies
end
FDIsum+=-9*sign(j-i)*dipole[j]/((r*(j-i))^7) # How much do the dipoles respond to the electron density?
end
FieldDI[i]=FDIsum
end
return FieldDI
end
" Decompose Hamiltonian into Diagonal/S/PE and Off-diag/J/KE elements"
function Decompose_H(H)
S=eye(N).*H # elementwise to select for just diagonal terms
J=H-S
return S,J
end
function Propagate(H,psi,dt)
U=eye(H) # identiy matrix same size + type as H
S,J=Decompose_H(H) # split into diagonal and off-diag terms
#Trotter decomposition
U*=expm(-im*J*dt/2*hbar)*expm(-im*S*dt/hbar)*expm(-im*J*dt/2*hbar)
psi = U*psi
#psi=eigvecs(H)[:,1]
En = eigvals(H)
normalisation = norm(psi)
psi = psi/normalisation
density = conj(psi).*psi
return psi, density
end
function UpdateDipole(Field,density,dipole)
alphainv = -1*ones(N)#(221*N*density)#
M = diagm(alphainv)
n=0
n2=1
for i in 1:N
for j in 1:N
if (j==i)
#M[n+j] = alpha[i]
continue
end
M[n+j] = 1/(r*(j-i))^3
end
n+=N
end
new_dipole = \(M,Field)
dipole += 0.5*(new_dipole-dipole)
return dipole
end
function UpdateEnergy(dipole,density,FieldDI, Field,S,E)
Vpq = (dipole.*(Field)+ density.*(FieldDI))
norm_vpq = norm(Vpq)
#plot(real(Vpq), label = "Vpq")
Vpp = dipole.*(FieldDI)
norm_vpp = norm(Vpp)
#plot!(real(Vpp), label = "Vpp")
Vqq = N*density.*(Field)
norm_vqq = norm(Vqq)
#plot!(real(Vqq), label = "Vqq")
#KE = S - V
#E=?
S=Vpq+Vpp+Vqq
norm_S = norm(S)
#plot!(real(S), label = "Total")
println("norms: vpq", norm_vpq, "vpp", norm_vpp, "vqq", norm_vqq, "total", norm_S)
H=diagm(E,-1)+diagm(S)+diagm(E,1) #build full matrix from diagonal elements; for comparison
return S,E,H
end
function UpdateEnergy2(dipole,Field,S,E)
S=zeros(N)
for i in 1:N
for j in 1:N
if (j==i)
continue # avoid infinity in self energies
end
S[i]+=sign(j-i)*dipole[j]/(r*(j-i))^2 # Contribution to site energy (1 e- at site) from dipoles
end
# @printf("Site: i %d SiteEnergy: S[i] %f\n",i,S[i])
end
H=diagm(E,-1)+diagm(S)+diagm(E,1)
return S,E,H
end
function nondispersive_wavepacket(x0, λ)
ψ=[ exp(-(x-x0)^2)*(cos(2π*(x-x0)/λ) - im * sin(2π*(x-x0)/λ) ) for x=1:N ]
return ψ
end
"""
Think about which way the energy, dipoles and fields should be....
Think about magnitudes of dipoles and densities vs Fields (hence which have the most weight for energy)
"""