-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathVCA_MAIN.f90
330 lines (304 loc) · 10.7 KB
/
VCA_MAIN.f90
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
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
module VCA_MAIN
USE VCA_VARS_GLOBAL
USE VCA_AUX_FUNX
USE VCA_SETUP
USE VCA_DIAG
USE VCA_OBSERVABLES
USE VCA_GREENS_FUNCTIONS
USE VCA_BATH_SETUP
USE VCA_EIGENSPACE
USE VCA_HAMILTONIAN
USE VCA_OMEGA
!
USE SF_LINALG, only: eigh,diag,kron,eye
USE SF_IOTOOLS, only: str,free_unit
USE SF_TIMER,only: start_timer,stop_timer
USE SF_MISC, only: assert_shape,sort_array
implicit none
private
!>INIT VCA SOLVER
interface vca_init_solver
module procedure :: vca_init_solver_serial
#ifdef _MPI
module procedure :: vca_init_solver_mpi
#endif
end interface vca_init_solver
!>
public :: vca_init_solver
!>VCA SOLVER
interface vca_solve
module procedure :: vca_solve_serial
#ifdef _MPI
module procedure :: vca_solve_mpi
#endif
end interface vca_solve
!>
public :: vca_solve
contains
!+-----------------------------------------------------------------------------+!
! PURPOSE: allocate and initialize one or multiple baths -+!
!+-----------------------------------------------------------------------------+!
subroutine vca_init_solver_serial(bath_h,bath_v)
complex(8),intent(inout),optional :: bath_h(:,:,:,:,:,:)
complex(8),intent(inout),optional :: bath_v(:,:,:,:,:,:)
logical,save :: isetup=.true.
!
write(LOGfile,"(A)")"INIT SOLVER FOR "//trim(file_suffix)
!
if(present(bath_h).and.present(bath_v))then
if(Nlat_bath<1 .or. Norb_bath<1)stop "VCA_INIT_SOLVER error: wrong bath dimensions"
call vca_allocate_bath(vca_bath)
call vca_init_bath(vca_bath)
else
write(LOGfile,"(A)") "Bath not present, setting Nbath to 0"
Nlat_bath=0
Norb_bath=0
endif
!Init Structure & memory
if(isetup)call init_cluster_structure()
!
if(isetup)call setup_global
!
if(vca_bath%status)call vca_deallocate_bath(vca_bath)
!
isetup=.false.
!
end subroutine vca_init_solver_serial
#ifdef _MPI
subroutine vca_init_solver_mpi(MpiComm,bath_h,bath_v)
integer :: MpiComm
complex(8),intent(inout),optional :: bath_h(:,:,:,:,:,:)
complex(8),intent(inout),optional :: bath_v(:,:,:,:,:,:)
logical :: check
logical,save :: isetup=.true.
integer :: i
!
call vca_set_mpicomm(MpiComm)
!
write(LOGfile,"(A)")"INIT SOLVER FOR "//trim(file_suffix)
!
if(present(bath_h).and.present(bath_v))then
if(Nlat_bath<1 .or. Norb_bath<1)stop "VCA_INIT_SOLVER error: wrong bath dimensions"
call vca_allocate_bath(vca_bath)
call vca_init_bath(vca_bath)
else
write(LOGfile,"(A)") "Bath not present, setting Nbath to 0"
Nlat_bath=0
Norb_bath=0
endif
!Init Structure & memory
if(isetup)call init_cluster_structure()
!
!Init bath:
!call set_hloc(Hloc)
!
!
if(isetup)call setup_global
!
if(vca_bath%status)call vca_deallocate_bath(vca_bath)
!
isetup=.false.
!
call vca_del_MpiComm()
!
end subroutine vca_init_solver_mpi
#endif
!+-----------------------------------------------------------------------------+!
!PURPOSE: Diag the cluster, reference system
!+-----------------------------------------------------------------------------+!
subroutine vca_solve_serial(Hloc,Hk,bath_h,bath_v)
complex(8),intent(in),dimension(:,:,:,:,:,:) :: Hloc ![Nlat,Nlat,Nspin,Nspin,Norb,Norb]
complex(8),intent(in),dimension(:,:,:,:,:,:,:) :: Hk ![Nlat,Nlat,Nspin,Nspin,Norb,Norb,Nktot]
complex(8),intent(inout),optional :: bath_h(:,:,:,:,:,:)
complex(8),intent(inout),optional :: bath_v(:,:,:,:,:,:)
!
integer :: Lk,Nsites
integer :: ilat,jlat
integer :: iorb,jorb
integer :: ispin
integer :: i,j
real(8) :: Tr,omega_integral,omegaprime
integer :: unit
!
write(LOGfile,"(A)")"SOLVING VCA"
!
if(rank(Hloc) .ne. 6) STOP "STOP: wrong cluster matrix dimensions"
if(rank(Hk) .ne. 7) STOP "STOP: wrong lattice matrix dimensions"
!
if(present(bath_h).and.present(bath_v))then
call assert_shape(bath_h,[Nlat_bath,Nlat_bath,Nspin,Nspin,Norb_bath,Norb_bath],"vca_solve","bath_h")
call assert_shape(bath_v,[Nlat ,Nlat_bath,Nspin,Nspin,Norb ,Norb_bath],"vca_solve","bath_v")
call vca_allocate_bath(vca_bath)
call vca_set_bath(bath_h,bath_v,vca_bath)
endif
!
!
select case(vca_sparse_H)
case (.true.)
spHtimesV_p => spMatVec_main
case (.false.)
spHtimesV_p => directMatVec_main
case default
stop "vca_solve_single ERROR: vca_sparse_H undefined"
end select
!
!
!GENERATE THE CLUSTER HAMILTONIAN AND THE HOPPING MATRIX FOR THE LATTICE
!
allocate(impHloc(Nlat,Nlat,Nspin,Nspin,Norb,Norb))
allocate(impHk(Nlat,Nlat,Nspin,Nspin,Norb,Norb,size(Hk,7)))
impHloc=zero
impHk=zero
call vca_set_Hcluster(Hloc)
call vca_set_Hk(Hk)
!
!GET CLUSTER GREEN'S FUNCTION AND GROUND STATE ENERGY
!
call diagonalize_cluster() !find target states by digonalization of Hamiltonian
call build_gf_cluster() !build the one-particle Green's functions and Self-Energies
if(print_observables)then
call observables_cluster() !obtain impurity observables as thermal averages.
call get_custom_observables() !obtain impurity observables as thermal averages.
endif
!
!CALCULATE THE VARIATIONAL GRAND POTENTIAL
!
omegaprime=0.d0
omega_integral=0.d0
!
if(finiteT)then
do i=1,state_list%size
omegaprime=omegaprime+exp(-beta*(es_return_energy(state_list,i)-state_list%emin))
enddo
omegaprime=state_list%emin-(1.d0/beta)*log(omegaprime)
omega_integral=frequency_integration_finite_t()
else
omegaprime=state_list%emin
omega_integral=frequency_integration()
endif
!
sft_potential = omegaprime-omega_integral
!
write(LOGfile,"(A,10f18.12,A)")"EGS PER SITE",state_list%emin/NLAT
write(LOGfile,"(A,10f18.12,A)")"OMEGA POTENTIAL PER SITE=",(state_list%emin-omega_integral)/NLAT
open(free_unit(unit),file="SFT_potential.vca",position='append')
write(unit,*)sft_potential,omegaprime,-omega_integral
close(unit)
!
!CLEAN UP
!
call es_delete_espace(state_list)
nullify(spHtimesV_p)
!
!
deallocate(impHloc)
deallocate(impHk)
if(vca_bath%status)call vca_deallocate_bath(vca_bath)
!
end subroutine vca_solve_serial
#ifdef _MPI
subroutine vca_solve_mpi(MpiComm,Hloc,Hk,bath_h,bath_v)
complex(8),intent(in),dimension(:,:,:,:,:,:) :: Hloc ![Nlat,Nlat,Nspin,Nspin,Norb,Norb]
complex(8),intent(in),dimension(:,:,:,:,:,:,:) :: Hk ![Nlat,Nlat,Nspin,Nspin,Norb,Norb,Nktot]
complex(8),intent(inout),optional :: bath_h(:,:,:,:,:,:),bath_v(:,:,:,:,:,:)
!
integer :: Lk,Nsites
integer :: ilat,jlat
integer :: iorb,jorb
integer :: ispin
integer :: i,j
real(8) :: Tr,omega_integral,omegaprime
integer :: unit
integer :: MpiComm
logical :: check
!
!SET THE LOCAL MPI COMMUNICATOR :
!
call vca_set_MpiComm(MpiComm)
!
if(rank(Hloc) .ne. 6) STOP "STOP: wrong cluster matrix dimensions"
if(rank(Hk) .ne. 7) STOP "STOP: wrong lattice matrix dimensions"
!
!
!
if(present(bath_h).and.present(bath_v))then
call assert_shape(bath_h,[Nlat_bath,Nlat_bath,Nspin,Nspin,Norb_bath,Norb_bath],"vca_solve","bath_h")
call assert_shape(bath_v,[Nlat ,Nlat_bath,Nspin,Nspin,Norb ,Norb_bath],"vca_solve","bath_v")
call vca_allocate_bath(vca_bath)
call vca_set_bath(bath_h,bath_v,vca_bath)
call vca_write_bath(vca_bath)
endif
!
select case(vca_sparse_H)
case (.true.)
spHtimesV_p => spMatVec_main
case (.false.)
spHtimesV_p => directMatVec_main
case default
stop "vca_solve_single ERROR: vca_sparse_H undefined"
end select
!
!GENERATE THE CLUSTER HAMILTONIAN AND THE HOPPING MATRIX FOR THE LATTICE
!
allocate(impHloc(Nlat,Nlat,Nspin,Nspin,Norb,Norb))
allocate(impHk(Nlat,Nlat,Nspin,Nspin,Norb,Norb,size(Hk,7)))
impHloc=zero
impHk=zero
call vca_set_Hcluster(Hloc)
call vca_set_Hk(Hk)
!
!GET CLUSTER GREEN'S FUNCTION AND GROUND STATE ENERGY
!
call diagonalize_cluster() !find target states by digonalization of Hamiltonian
call build_gf_cluster() !build the one-particle Green's functions and Self-Energies
if(print_observables)then
call observables_cluster() !obtain impurity observables as thermal averages.
call get_custom_observables() !obtain impurity observables as thermal averages.
endif
!call save_gfprime("gfprime",use_formatted=.true.)
!call read_gfprime("gfprime",use_formatted=.true.)
!call reconstruct_g()
!
!CALCULATE THE VARIATIONAL GRAND POTENTIAL
!
omegaprime=0.d0
omega_integral=0.d0
!
if(finiteT)then
do i=1,state_list%size
omegaprime=omegaprime+exp(-beta*(es_return_energy(state_list,i)-state_list%emin))
enddo
omegaprime=state_list%emin-(1.d0/beta)*log(omegaprime)
if(MPIMASTER) omega_integral=frequency_integration_finite_t()
else
omegaprime=state_list%emin
if(MPIMASTER) omega_integral=frequency_integration()
endif
!
if(MPIMASTER) sft_potential = omegaprime-omega_integral
call Bcast_MPI(MpiComm,sft_potential)
!
if(MPIMASTER)then
write(LOGfile,"(A,10f18.12,A)")"EGS PER SITE",omegaprime/NLAT
write(LOGfile,"(A,10f18.12,A)")"OMEGA POTENTIAL PER SITE=",(omegaprime-omega_integral)/NLAT
open(free_unit(unit),file="SFT_potential.vca",position='append')
write(unit,*)sft_potential,omegaprime,-omega_integral
close(unit)
endif
!
call barrier_mpi(MpiComm)
!
!CLEAN UP
!
call es_delete_espace(state_list)
nullify(spHtimesV_p)
!
deallocate(impHloc)
deallocate(impHk)
if(vca_bath%status)call vca_deallocate_bath(vca_bath)
call vca_del_MpiComm()
!
end subroutine vca_solve_mpi
#endif
end module VCA_MAIN
! !