-
Notifications
You must be signed in to change notification settings - Fork 0
/
tebd_callback_routines_imps.f90
356 lines (235 loc) · 13.7 KB
/
tebd_callback_routines_imps.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
MODULE tebd_callback_routines_imps
USE utility
USE iteration_helper
USE simulation_parameters
USE observables_imps
IMPLICIT NONE
interface create_initial_state_vec
module procedure create_psi_vec_from_input_vec
module procedure create_rho_vec_from_input_vec
module procedure create_rho_vec_from_rho_mat
end interface create_initial_state_vec
INTEGER :: DUMP_NUM, G2_SUMMARY
private create_psi_vec_from_input_vec, create_rho_vec_from_input_vec, create_rho_vec_from_rho_mat
CONTAINS
!!! Test if both 1P && 2P observables have converged !!!
SUBROUTINE test_tebd_convergence(is_converged, new_obs_2P, old_obs_2P, new_obs_1P, old_obs_1P, rho_onesite, rho_twosite, tebd, iter, time)
LOGICAL, INTENT(INOUT) :: is_converged !Whether observables have converged
COMPLEX(KIND=DP), ALLOCATABLE, INTENT(INOUT) :: new_obs_2P(:,:), old_obs_2P(:,:) !New && old values of observables (to be compared)
COMPLEX(KIND=DP), ALLOCATABLE, INTENT(INOUT) :: new_obs_1P(:), old_obs_1P(:) !New && old values of observables (to be compared)
COMPLEX(KIND=DP), INTENT(IN) :: rho_onesite(:), rho_twosite(:,:,:) !One-site && Two-site reduced density matrices
TYPE(tebd_params), INTENT(IN) :: tebd !Time evolution -- propagation params (dt, N_steps, test_interval)
INTEGER, INTENT(IN) :: iter !Time evolution -- iter we are at
REAL(KIND=DP), INTENT(IN) :: time !Time evolution -- time we are at (corresponds to iter)
LOGICAL :: is_converged_1P, is_converged_2P
!! Initialize
is_converged = .FALSE.
!! Check 1P observables
CALL test_tebd_conv_1P(is_converged_1P, new_obs_1P, old_obs_1P, &
& rho_onesite, tebd, iter, time)
!! Check 2P observables
CALL test_tebd_conv_2P(is_converged_2P, new_obs_2P, old_obs_2P, &
& rho_onesite, rho_twosite, tebd, iter, time)
!! Both 1P && 2P must be converged
is_converged = is_converged_1P .AND. is_converged_2P
END SUBROUTINE test_tebd_convergence
!!! Test if 1P observables have converged !!!
SUBROUTINE test_tebd_conv_1P(is_converged, new_obs, old_obs, rho_onesite, tebd, iter, time)
USE datafile_utility
LOGICAL, INTENT(INOUT) :: is_converged !whether observables have converged
COMPLEX(KIND=DP), ALLOCATABLE, INTENT(INOUT) :: new_obs(:), old_obs(:) !new && old values of observables (to be compared)
COMPLEX(KIND=DP), INTENT(IN) :: rho_onesite(:) !one-site reduced density matrix for computing 1P observables
TYPE(tebd_params), INTENT(IN) :: tebd !Time evolution -- propagation params (dt, N_steps, test_interval)
INTEGER, INTENT(IN) :: iter !Time evolution -- iter we are on
REAL(KIND=DP), INTENT(IN) :: time !Time evolution -- time we are at (corresponds to iter)
INTEGER :: Nop, i !operator indexing
LOGICAL, ALLOCATABLE :: obs_converged(:) !convergence tracking for individual obs
REAL(KIND=DP) :: etaR !rescaled eta
COMPLEX(KIND=DP) :: norm !Norm of TN
!number of operators
Nop = SIZE(TRACE_OP_1P, 1)
!Convergence tracking for diff obs
ALLOCATE(obs_converged(Nop)); obs_converged(:) = .FALSE.
!Allocate storage for transient observables
IF(.NOT. ALLOCATED(new_obs)) ALLOCATE(new_obs(Nop))
IF(.NOT. ALLOCATED(old_obs)) ALLOCATE(old_obs(Nop))
!Rescaled convergence cutoff
etaR = (tebd % eta) * (tebd % test_interval) * ABS(tebd % dt)
!Get norm of rho TN (i.e. trace)
norm = rho_onesite(1)
!Compute all observables on the list i={1,nop}
DO i=1,Nop
CALL compute_obs_1P(new_obs(i), rho_onesite, TRACE_OP_1P(i,:,:))
new_obs(i) = new_obs(i)/norm
CALL test_convergence(obs_converged(i), new_obs(i), old_obs(i), &
& etaR, TRACE_OUT_1P(i), iter, tebd % test_interval, time, norm)
end DO
!Check convergence
is_converged = isTrueVec(obs_converged, 'AND')
!Set old --> old = new for next iter
old_obs = new_obs
END SUBROUTINE test_tebd_conv_1P
!!! Test if 2P observables have converged !!!
SUBROUTINE test_tebd_conv_2P(is_converged, new_obs, old_obs, rho_onesite, rho_twosite, tebd, iter, time)
USE datafile_utility
LOGICAL, INTENT(INOUT) :: is_converged !Whether observables have converged
COMPLEX(KIND=DP), ALLOCATABLE, INTENT(INOUT) :: new_obs(:,:), old_obs(:,:) !New && old values of observables (to be compared)
COMPLEX(KIND=DP), INTENT(IN) :: rho_onesite(:) !One-site reduced density matrix
COMPLEX(KIND=DP), INTENT(IN) :: rho_twosite(:,:,:) !Two-site reduced density matrix
TYPE(tebd_params), INTENT(IN) :: tebd !Time evolution -- propagation params (dt, N_steps, test_interval)
INTEGER, INTENT(IN) :: iter !Time evolution -- iter we are at
REAL(KIND=DP), INTENT(IN) :: time !Time evolution -- time we are at (corresponds to iter)
INTEGER :: Nop, i !operator indexing
INTEGER :: min_sep, max_sep, Nsep, sep !separations
LOGICAL, ALLOCATABLE :: obs_converged(:,:) !convergence tracking for individual obs
REAL(KIND=DP) :: etaR !rescaled eta
COMPLEX(KIND=DP) :: norm !Norm of TN
!number of operators
Nop = SIZE(TRACE_OP_2P, 1)
!separations
min_sep = tebd % min_sep
max_sep = tebd % max_sep
Nsep = max_sep - min_sep + 1
!Convergence tracking for diff obs
ALLOCATE(obs_converged(Nop, Nsep)); obs_converged(:,:) = .FALSE.
!Allocate storage for transient observables
IF(.NOT. ALLOCATED(new_obs)) ALLOCATE(new_obs(Nop, Nsep))
IF(.NOT. ALLOCATED(old_obs)) ALLOCATE(old_obs(Nop, Nsep))
!Rescaled convergence cutoff
etaR = (tebd % eta) * (tebd % test_interval) * ABS(tebd % dt)
!Get norm of rho TN (i.e. trace)
norm = rho_onesite(1)
!Compute all 2P observables on the list i={1,nop}
DO i=1,Nop
!compute obs
CALL compute_obs_2P(new_obs(i,:), rho_twosite, rho_onesite, TRACE_OP_2P(i,:,:,:), min_sep, max_sep)
!normalize obs, and test convergence
DO sep = min_sep, max_sep
new_obs(i, sep+1-min_sep) = new_obs(i, sep+1-min_sep) / norm**(sep/2 + 1)
CALL test_convergence(obs_converged(i,sep), new_obs(i,sep), old_obs(i,sep), &
& etaR, TRACE_OUT_2P(i,sep), iter, tebd % test_interval, time, norm)
end DO
end DO
!Check convergence
is_converged = isTrueMat(obs_converged, 'AND')
!Set old --> old = new for next iter
old_obs = new_obs
END SUBROUTINE test_tebd_conv_2P
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Lambda Convergence !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Test if Lambdas have converged !!!
SUBROUTINE test_lambda_convergence(is_converged, newLambda, oldLambda, tebd, datafile, iter, time)
LOGICAL, INTENT(INOUT) :: is_converged !whether observables have converged
TYPE(block_lambda), INTENT(IN) :: newLambda(:) !new Lambdas
TYPE(block_lambda), ALLOCATABLE, INTENT(INOUT) :: oldLambda(:) !old Lambdas
TYPE(tebd_params), INTENT(IN) :: tebd !TEBD params
INTEGER, INTENT(INOUT) :: datafile !datafile for convergence data
INTEGER, INTENT(IN) :: iter !Time evolution -- iter we are on
REAL(KIND=DP), INTENT(IN) :: time !Time evolution -- time we are at (corresponds to iter)
INTEGER :: N_sites, site !! lambda sites
LOGICAL, ALLOCATABLE :: lambda_converged(:) !! convergence tracking for individual lambdas
REAL(KIND=DP) :: etaR !! Rescaled eta
!! Number of lambdas
N_sites = SIZE(newLambda)
!! Convergence tracking for diff lambdas
ALLOCATE(lambda_converged(N_sites)); lambda_converged(:) = .FALSE.
!! Rescaled convergence cutoff
etaR = (tebd % eta) * (tebd % test_interval) * ABS(tebd % dt)
!! Test Lambda convergence
DO site=1,N_sites
CALL test_convergence_vec(lambda_converged(site), newLambda(site) % m, oldLambda(site) % m, etaR, &
& datafile, iter, TEBD % test_interval, xval=TEBD % dt, PRINT_ON=.TRUE.)
end DO
!! Check for convergence
is_converged = isTrueVec(lambda_converged, 'AND')
!! Set old --> old = new for next iter
CALL copy_lambda_block(oldLambda, newLambda)
END SUBROUTINE test_lambda_convergence
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Decide which param to use in the param-loop, based on "xname" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE read_xname(x, xname, v1, v1_str, v2, v2_str, v3, v3_str, v4, v4_str)
REAL(KIND=DP), INTENT(IN) :: x
CHARACTER(LEN=*), INTENT(IN) :: xname
REAL(KIND=DP), INTENT(INOUT) :: v1
CHARACTER(LEN=*), INTENT(IN) :: v1_str
REAL(KIND=DP), INTENT(INOUT), OPTIONAL :: v2
CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: v2_str
REAL(KIND=DP), INTENT(INOUT), OPTIONAL :: v3
CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: v3_str
REAL(KIND=DP), INTENT(INOUT), OPTIONAL :: v4
CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: v4_str
!Local copies of optional string vars
CHARACTER(LEN=16) :: v2_strX, v3_strX, v4_strX
!Set default values (If var not present -- set its string to INVALID_STRING, to guarantee it's not selected)
v2_strX = OptArg(v2_str, "INVALID_STRING")
v3_strX = OptArg(v3_str, "INVALID_STRING")
v4_strX = OptArg(v4_str, "INVALID_STRING")
!Decide which param to use
IF(xname .EQ. v1_str) THEN
WRITE(*,*) "Setting ", v1_str, " = ", x
v1 = x
ELSEIF(xname .EQ. v2_strX) THEN
WRITE(*,*) "Setting ", v2_str, " = ", x
CALL setOptVar(v2, x, "set_xvar: argument must be present")
ELSEIF(xname .EQ. v3_strX) THEN
WRITE(*,*) "Setting ", v3_str, " = ", x
CALL setOptVar(v3, x, "set_xvar: argument must be present")
ELSEIF(xname .EQ. v4_strX) THEN
WRITE(*,*) "Setting ", v4_str, " = ", x
CALL setOptVar(v4, x, "set_xvar: argument must be present")
ELSE
CALL invalid_flag("set_xvar: invalid xname ", xname)
end IF
END SUBROUTINE read_xname
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Create initial state vector !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Create initial state vector directly by reading in the input vector
SUBROUTINE create_psi_vec_from_input_vec(psi_vec, s0, psi_flag)
USE simulation_parameters, ONLY: local_dim
COMPLEX(KIND=DP), ALLOCATABLE, INTENT(INOUT) :: psi_vec(:)
REAL(KIND=DP), INTENT(INOUT) :: s0(:)
CHARACTER(LEN=*), INTENT(IN) :: psi_flag
!Consistency checks
IF(local_dim .LT. 2) CALL invalid_value("create_psi_vec_from_input_vec: invalid local dim ", local_dim)
CALL check_sizes_equal(SIZE(s0), local_dim, "create_psi_vec_from_input_vec: s0 must have a size of local_dim ")
ALLOCATE(psi_vec(local_dim))
!Set psi_vec = sqrt of probability amplitude coefficients
s0 = s0/sqrt(sum(s0*s0))
psi_vec = sqrt(s0)
END SUBROUTINE create_psi_vec_from_input_vec
!Create initial rho vector directly by reading in the input vector
SUBROUTINE create_rho_vec_from_input_vec(rho_vec, s0)
USE simulation_parameters, ONLY: local_dim
COMPLEX(KIND=DP), ALLOCATABLE, INTENT(INOUT) :: rho_vec(:)
REAL(KIND=DP), INTENT(INOUT) :: s0(:)
!Consistency checks
IF(local_dim .LT. 2) CALL invalid_value("create_rho_vec_from_input_vec: invalid local dim ", local_dim)
CALL check_sizes_equal(SIZE(s0), local_dim - 1, "create_rho_vec_from_input_vec: s0 must have a size of local_dim - 1 ")
ALLOCATE(rho_vec(local_dim))
!Read rho_vec (NB. to ensure positivity, we must make sure s0 is not too outside the Bloch sphere)
IF(SUM(s0*s0) .GT. 1.0D0) s0 = s0/sqrt(sum(s0*s0))
rho_vec(1) = 1.0D0
rho_vec(2:local_dim) = s0
END SUBROUTINE create_rho_vec_from_input_vec
!Create initial rho vector from input density matrix
SUBROUTINE create_rho_vec_from_rho_mat(rho_vec, rho_mat, represent)
USE simulation_parameters, ONLY: local_dim
COMPLEX(KIND=DP), ALLOCATABLE, INTENT(INOUT) :: rho_vec(:)
COMPLEX(KIND=DP), INTENT(IN) :: rho_mat(:,:)
INTERFACE
FUNCTION represent(rho_mat)
USE utility
USE simulation_parameters, ONLY: local_dim
COMPLEX(KIND=DP), INTENT(IN) :: rho_mat(INT(sqrt(local_dim*1.0d0)), INT(sqrt(local_dim*1.0d0)))
COMPLEX(KIND=DP) :: represent(local_dim)
end FUNCTION represent
end INTERFACE
COMPLEX(KIND=DP) :: norm
!Consistency checks
IF(local_dim .LT. 2) CALL invalid_value("create_rho_vec_from_rho_mat: invalid local dim ", local_dim)
!Rho normalization (trace)
norm = TRACE(rho_mat)
!Represent rho matrix as a state vector using Gell-Mann basis matrices
ALLOCATE(rho_vec(local_dim))
rho_vec = represent(rho_mat / norm)
END SUBROUTINE create_rho_vec_from_rho_mat
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
END MODULE tebd_callback_routines_imps