-
Notifications
You must be signed in to change notification settings - Fork 1
/
VCA_DIAG.f90
550 lines (508 loc) · 19.8 KB
/
VCA_DIAG.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
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
module VCA_DIAG
USE VCA_INPUT_VARS
USE VCA_VARS_GLOBAL
USE VCA_EIGENSPACE
USE VCA_SETUP
USE VCA_HAMILTONIAN
!
USE SF_CONSTANTS
USE SF_STAT
USE SF_SP_LINALG
USE SF_LINALG, only: eigh
USE SF_TIMER, only: start_timer,stop_timer,eta
USE SF_IOTOOLS, only: reg,free_unit,file_length
USE SF_MISC, only: assert_shape
!
implicit none
private
!>diagonalize cluster Hamiltonian
public :: diagonalize_cluster
contains
!+-------------------------------------------------------------------+
!>Setup the Hilbert space, create the Hamiltonian, get the
! GS, build the Green's functions calling all the necessary routines
!+-------------------------------------------------------------------+
subroutine diagonalize_cluster
call vca_pre_diag
call vca_diag_c
!call vca_post_diag(.true.)
end subroutine diagonalize_cluster
!+-------------------------------------------------------------------+
!PURPOSE : diagonalize the Hamiltonian in each sector
!+------------------------------------------------------------------+
subroutine vca_diag_c
integer :: nup,ndw,isector,dim
integer :: DimUps(Ns_Ud),DimUp
integer :: DimDws(Ns_Ud),DimDw
integer :: Nups(Ns_Ud)
integer :: Ndws(Ns_Ud)
integer :: sz,nt
integer :: i,j,iter,unit,vecDim,PvecDim
integer :: Nitermax,Neigen,Nblock
integer :: iter_spectrum
real(8) :: oldzero,enemin,Ei,neigen_sector_error
real(8),allocatable :: eig_values(:)
complex(8),allocatable :: eig_basis(:,:),eig_basis_tmp(:,:)
logical :: lanc_solve,Tflag,lanc_verbose,bool,converged_spectrum
!
if(state_list%status)call es_delete_espace(state_list)
state_list=es_init_espace()
oldzero=1000.d0
if(MPIMASTER)then
write(LOGfile,"(A)")"Diagonalize cluster Hc+Hint:"
call start_timer()
endif
!
lanc_verbose=.false.
if(verbose>2)lanc_verbose=.true.
!
converged_spectrum=.false.
iter_spectrum=0
!
anneal: do while(.not.converged_spectrum.AND.iter_spectrum<lanc_niter_spectrum)
iter_spectrum=iter_spectrum+1
call es_free_espace(state_list)
!print*,"Spectrum annealing iteration:",iter_spectrum,"of ",lanc_niter_spectrum
iter=0
sector: do isector=1,Nsectors
if(.not.sectors_mask(isector))cycle sector
if(.not.twin_mask(isector))cycle sector !cycle loop if this sector should not be investigated
iter=iter+1
call get_Nup(isector,nups)
call get_Ndw(isector,ndws)
Tflag = twin_mask(isector).AND.vca_twin
bool=.true.
do i=1,Ns_ud
Bool=Bool.AND.(nups(i)/=ndws(i))
enddo
Tflag=Tflag.AND.Bool
!
Dim = getdim(isector)
!
select case (lanc_method)
case default !use P-ARPACK
Neigen = min(dim,neigen_sector(isector))
Nitermax = min(dim,lanc_niter)
Nblock = min(dim,lanc_ncv_factor*max(Neigen,lanc_nstates_sector) + lanc_ncv_add)
case ("lanczos")
Neigen = 1
Nitermax = min(dim,lanc_niter)
Nblock = 1
end select
!
lanc_solve = .true.
if(Neigen==dim)lanc_solve=.false.
if(dim<=max(lanc_dim_threshold,MPISIZE))lanc_solve=.false.
!
if(MPIMASTER)then
call get_DimUp(isector,DimUps) ; DimUp = product(DimUps)
call get_DimDw(isector,DimDws) ; DimDw = product(DimDws)
if(verbose>=3)then
if(lanc_solve)then
write(LOGfile,"(1X,I4,A,I4,A6,"&
//str(Ns_Ud)//"I3,A6,"&
//str(Ns_Ud)//"I3,A7,"&
//str(Ns_Ud)//"I6,"//str(Ns_Ud)//"I6,I15,A12,3I6)")&
iter,"-Solving sector:",isector,", nup:",nups,", ndw:",ndws,", dims=",&
DimUps,DimDws,getdim(isector),", Lanc Info:",Neigen,Nitermax,Nblock
else
write(LOGfile,"(1X,I4,A,I4,A6,"&
//str(Ns_Ud)//"I3,A6,"&
//str(Ns_Ud)//"I3,A7,"&
//str(Ns_Ud)//"I6,"//str(Ns_Ud)//"I6,I15)")&
iter,"-Solving sector:",isector,", nup:",nups,", ndw:",ndws,", dims=",&
DimUps,DimDws,getdim(isector)
endif
elseif(verbose==1.OR.verbose==2)then
call eta(iter,count(twin_mask),LOGfile)
endif
endif
!
if(allocated(eig_values))deallocate(eig_values)
if(allocated(eig_basis))deallocate(eig_basis)
if(lanc_solve)then
allocate(eig_values(Neigen))
eig_values=0d0
!
call build_Hv_sector(isector)
!
vecDim = vecDim_Hv_sector(isector)
allocate(eig_basis(vecDim,Neigen))
eig_basis=zero
!
select case (lanc_method)
case default !use P-ARPACK
#ifdef _MPI
if(MpiStatus)then
call sp_eigh(MpiComm,spHtimesV_p,eig_values,eig_basis,&
Nblock,&
Nitermax,&
tol=lanc_tolerance,&
iverbose=(verbose>3))
else
call sp_eigh(spHtimesV_p,eig_values,eig_basis,&
Nblock,&
Nitermax,&
tol=lanc_tolerance,&
iverbose=(verbose>3))
endif
#else
call sp_eigh(spHtimesV_p,eig_values,eig_basis,&
Nblock,&
Nitermax,&
tol=lanc_tolerance,&
iverbose=(verbose>=3))
#endif
case ("lanczos") !use Simple Lanczos
#ifdef _MPI
if(MpiStatus)then
call sp_lanc_eigh(MpiComm,spHtimesV_p,eig_values(1),eig_basis(:,1),Nitermax,&
iverbose=(verbose>3),threshold=lanc_tolerance)
else
call sp_lanc_eigh(spHtimesV_p,eig_values(1),eig_basis(:,1),Nitermax,&
iverbose=(verbose>3),threshold=lanc_tolerance)
endif
#else
call sp_lanc_eigh(spHtimesV_p,eig_values(1),eig_basis(:,1),Nitermax,&
iverbose=(verbose>3),threshold=lanc_tolerance)
#endif
end select
!
if(MpiMaster.AND.verbose>3)write(LOGfile,*)""
call delete_Hv_sector()
call Bcast_MPI(MpiComm,eig_values)
else
allocate(eig_values(Dim)) ; eig_values=0d0
allocate(eig_basis_tmp(Dim,Dim)) ; eig_basis_tmp=zero
call build_Hv_sector(isector,eig_basis_tmp)
if(MpiMaster)call eigh(eig_basis_tmp,eig_values,jobz='V',uplo='U')
if(dim==1)eig_basis_tmp(dim,dim)=one
!
call delete_Hv_sector()
#ifdef _MPI
if(MpiStatus)then
call Bcast_MPI(MpiComm,eig_values)
vecDim = vecDim_Hv_sector(isector)
allocate(eig_basis(vecDim,Neigen)) ; eig_basis=zero
call scatter_basis_MPI(MpiComm,eig_basis_tmp,eig_basis)
else
allocate(eig_basis(Dim,Neigen)) ; eig_basis=zero
eig_basis = eig_basis_tmp(:,1:Neigen)
endif
#else
allocate(eig_basis(Dim,Neigen)) ; eig_basis=zero
eig_basis = eig_basis_tmp(:,1:Neigen)
#endif
endif
!
if(verbose>=3)then
write(LOGfile,*)"EigValues: ",eig_values(:Neigen)
write(LOGfile,*)""
write(LOGfile,*)""
endif
!
if(finiteT)then
do i=1,Neigen
call es_add_state(state_list,eig_values(i),eig_basis(:,i),isector,twin=Tflag,size=lanc_nstates_total,verbose=(verbose>=3))
enddo
else
do i=1,Neigen
enemin = eig_values(i)
if (enemin < oldzero-10.d0*gs_threshold)then
oldzero=enemin
call es_free_espace(state_list)
call es_add_state(state_list,enemin,eig_basis(:,i),isector,twin=Tflag)
elseif(abs(enemin-oldzero) <= gs_threshold)then
oldzero=min(oldzero,enemin)
call es_add_state(state_list,enemin,eig_basis(:,i),isector,twin=Tflag)
endif
enddo
endif
!
if(MPIMASTER)then
unit=free_unit()
open(unit,file="eigenvalues_list"//reg(file_suffix)//".ed",position='append',action='write')
call print_eigenvalues_list(isector,eig_values(1:Neigen),unit,lanc_solve,mpiAllThreads)
close(unit)
endif
!
if(allocated(eig_values))deallocate(eig_values)
if(allocated(eig_basis_tmp))deallocate(eig_basis_tmp)
if(allocated(eig_basis))deallocate(eig_basis)
!
enddo sector
!
if(finiteT)then
converged_spectrum = (exp(-beta*(state_list%emax-state_list%emin)) < cutoff)
else
!> at T=0 you only need to solve one: there is no annealing of the spectrum
converged_spectrum=.true.
endif
!
call vca_post_diag(.true.)
!
!
enddo anneal
!
call stop_timer(unit=LOGfile)
!
close(unit)
!
end subroutine vca_diag_c
!###################################################################################################
!
! PRE-PROCESSING ROUTINES
!
!###################################################################################################
subroutine vca_pre_diag
integer :: Indices(2*Ns_Ud),Jndices(2*Ns_Ud)
integer :: Nups(Ns_ud),Ndws(Ns_ud)
integer :: Jups(Ns_ud),Jdws(Ns_ud)
integer :: i,iud,iorb
integer :: isector,jsector
integer :: unit,unit2,status,istate,ishift,isign
logical :: IOfile
integer :: list_len
integer,dimension(:),allocatable :: list_sector
!
sectors_mask=.true.
!
if(vca_sectors)then
inquire(file="sectors_list"//reg(file_suffix)//".restart",exist=IOfile)
if(IOfile)then
sectors_mask=.false.
write(LOGfile,"(A)")"Analysing sectors_list to reduce sectors scan:"
list_len=file_length("sectors_list"//reg(file_suffix)//".restart")
!
open(free_unit(unit),file="sectors_list"//reg(file_suffix)//".restart",status="old")
open(free_unit(unit2),file="list_of_sectors"//reg(file_suffix)//".ed")
do istate=1,list_len
read(unit,*,iostat=status)Indices
call get_Sector(Indices,Ns_Orb,isector)
sectors_mask(isector)=.true.
write(unit2,*)isector,sectors_mask(isector),Indices
!
do i=1,2*Ns_Ud
do ishift=1,vca_sectors_shift
do isign=-1,1,2
Jndices = Indices
Jndices(i) = Indices(i) + isign*ishift
call get_Sector(Jndices,Ns_Orb,jsector)
sectors_mask(jsector)=.true.
write(unit2,*)jsector,sectors_mask(jsector),Jndices
enddo
enddo
enddo
!
enddo
close(unit)
close(unit2)
!
endif
endif
!
end subroutine vca_pre_diag
!###################################################################################################
!
! POST-PROCESSING ROUTINES
!
!###################################################################################################
!+-------------------------------------------------------------------+
!PURPOSE : analyse the spectrum and print some information after
!lanczos diagonalization.
!+------------------------------------------------------------------+
subroutine vca_post_diag(iprint)
logical :: iprint
integer :: nup,ndw,sz,n,isector,dim
integer :: istate
integer :: i,unit
integer :: nups(Ns_Ud),ndws(Ns_Ud),Indices(2*Ns_Ud)
integer :: Nsize,NtoBremoved,nstates_below_cutoff
integer :: numgs
real(8) :: Egs,Ei,Ec,Etmp
type(histogram) :: hist
real(8) :: hist_a,hist_b,hist_w
integer :: hist_n
integer,allocatable :: list_sector(:),count_sector(:)
!
!POST PROCESSING:
if(MPIMASTER)then
open(free_unit(unit),file="state_list"//reg(file_suffix)//".ed")
call save_state_list(unit)
close(unit)
endif
if(verbose>=2)call print_state_list(LOGfile)
!
zeta_function=0d0
Egs = state_list%emin
if(finiteT)then
do i=1,state_list%size
ei = es_return_energy(state_list,i)
zeta_function = zeta_function + exp(-beta*(Ei-Egs))
enddo
else
zeta_function=real(state_list%size,8)
end if
!
!
numgs=es_return_gs_degeneracy(state_list,gs_threshold)
if(numgs>Nsectors)stop "diag: too many gs"
if(MPIMASTER.AND.verbose>=2)then
do istate=1,numgs
isector = es_return_sector(state_list,istate)
Egs = es_return_energy(state_list,istate)
call get_Nup(isector,Nups)
call get_Ndw(isector,Ndws)
write(LOGfile,"(A,F20.12,"//str(Ns_Ud)//"I4,"//str(Ns_Ud)//"I4)")'Egs =',Egs,nups,ndws
enddo
!omega_potential = -1d0/beta*log(zeta_function)
write(LOGfile,"(A,F20.12)")'Z =',zeta_function
!write(LOGfile,"(A,F20.12)")'Omega =',omega_potential
endif
!
!
!get histogram distribution of the sector contributing to the evaluated spectrum:
!go through states list and update the neigen_sector(isector) sector-by-sector
if(.not.finiteT)then
!generate a sector_list to be reused in case we want to reduce sectors scan
open(free_unit(unit),file="sectors_list"//reg(file_suffix)//".restart")
do istate=1,state_list%size
isector = es_return_sector(state_list,istate)
call get_Indices(isector,Ns_Orb,Indices)
write(unit,*)Indices
enddo
close(unit)
else
!get histogram distribution of the sector contributing to the evaluated spectrum:
!go through states list and update the neigen_sector(isector) sector-by-sector
if(MPIMASTER)then
unit=free_unit()
open(unit,file="histogram_states"//reg(file_suffix)//".ed",position='append')
hist_n = Nsectors
hist_a = 1d0
hist_b = dble(Nsectors)
hist_w = 1d0
hist = histogram_allocate(hist_n)
call histogram_set_range_uniform(hist,hist_a,hist_b)
do i=1,state_list%size
isector = es_return_sector(state_list,i)
call histogram_accumulate(hist,dble(isector),hist_w)
enddo
call histogram_print(hist,unit)
write(unit,*)""
close(unit)
endif
!
!
!
allocate(list_sector(state_list%size),count_sector(Nsectors))
!get the list of actual sectors contributing to the list
do i=1,state_list%size
list_sector(i) = es_return_sector(state_list,i)
enddo
!count how many times a sector appears in the list
do i=1,Nsectors
count_sector(i) = count(list_sector==i)
enddo
!adapt the number of required Neig for each sector based on how many
!appeared in the list.
do i=1,Nsectors
if(any(list_sector==i))then !if(count_sector(i)>1)then
neigen_sector(i)=neigen_sector(i)+1
else
neigen_sector(i)=neigen_sector(i)-1
endif
!prevent Neig(i) from growing unbounded but
!try to put another state in the list from sector i
if(neigen_sector(i) > count_sector(i))neigen_sector(i)=count_sector(i)+1
if(neigen_sector(i) <= 0)neigen_sector(i)=1
enddo
!check if the number of states is enough to reach the required accuracy:
!the condition to fullfill is:
! exp(-beta(Ec-Egs)) < \epsilon_c
! if this condition is violated then required number of states is increased
! if number of states is larger than those required to fullfill the cutoff:
! trim the list and number of states.
Egs = state_list%emin
Ec = state_list%emax
Nsize= state_list%size
if(exp(-beta*(Ec-Egs)) > cutoff)then
lanc_nstates_total=lanc_nstates_total + lanc_nstates_step
if(MPIMASTER)write(LOGfile,"(A,I4)")"Increasing lanc_nstates_total:",lanc_nstates_total
else
! !Find the energy level beyond which cutoff condition is verified & cut the list to that size
write(LOGfile,*)
isector = es_return_sector(state_list,state_list%size)
Ei = es_return_energy(state_list,state_list%size)
do while ( exp(-beta*(Ei-Egs)) <= cutoff )
if(verbose>=1.AND.MPIMASTER)write(LOGfile,"(A,I4,I5)")"Trimming state:",isector,state_list%size
call es_pop_state(state_list)
isector = es_return_sector(state_list,state_list%size)
Ei = es_return_energy(state_list,state_list%size)
enddo
if(verbose>=1.AND.MPIMASTER)then
write(LOGfile,*)"Trimmed state list:"
call print_state_list(LOGfile)
endif
!
lanc_nstates_total=max(state_list%size,lanc_nstates_step)+lanc_nstates_step
write(LOGfile,"(A,I4)")"Adjusting lanc_nstates_total to:",lanc_nstates_total
!
endif
endif
if(verbose>=1)write(LOGfile,*)""
if(verbose>=1)write(LOGfile,*)""
end subroutine vca_post_diag
subroutine print_state_list(unit)
integer :: indices(2*Ns_Ud),isector
integer :: istate
integer :: unit
real(8) :: Estate
if(MPIMASTER)then
write(unit,"(A1,A6,A18,2x,A19,1x,2A10,A12)")"#","i","E_i","exp(-(E-E0)/T)","Sect","Dim","Indices:"
do istate=1,state_list%size
Estate = es_return_energy(state_list,istate)
isector = es_return_sector(state_list,istate)
write(unit,"(i6,f18.12,2x,ES19.12,1x,2I10)",advance='no')&
istate,Estate,exp(-beta*(Estate-state_list%emin)),isector,getdim(isector)
call get_Indices(isector,Ns_Orb,Indices)
write(unit,"("//str(2*Ns_Ud)//"I4)")Indices
enddo
endif
end subroutine print_state_list
subroutine save_state_list(unit)
integer :: indices(2*Ns_Ud),isector
integer :: istate
integer :: unit
if(MPIMASTER)then
do istate=1,state_list%size
isector = es_return_sector(state_list,istate)
call get_Indices(isector,Ns_Orb,Indices)
write(unit,"(i8,i12,"//str(2*Ns_Ud)//"i8)")istate,isector,Indices
enddo
endif
end subroutine save_state_list
subroutine print_eigenvalues_list(isector,eig_values,unit,lanc,allt)
integer :: isector
real(8),dimension(:) :: eig_values
integer :: unit,i,indices(2*Ns_Ud)
logical :: lanc,allt
if(MPIMASTER)then
if(lanc)then
if(allt)then
write(unit,"(A9,A15)")" # Sector","Indices"
else
write(unit,"(A10,A15)")" #T Sector","Indices"
endif
else
write(unit,"(A10,A15)")" #X Sector","Indices"
endif
call get_Indices(isector,Ns_Orb,Indices)
write(unit,"(I9,"//str(2*Ns_Ud)//"I6)")isector,Indices
do i=1,size(eig_values)
write(unit,*)eig_values(i)
enddo
write(unit,*)""
endif
end subroutine print_eigenvalues_list
END MODULE VCA_DIAG