PROGRESS  master
prg_modelham_mod.F90
Go to the documentation of this file.
1 
5 
6  use bml
9 
10 
11  implicit none
12 
13  private !Everything is private by default
14 
15  integer, parameter :: dp = kind(1.0d0)
16 
18  type, public :: mham_type
19  integer :: norbs, seed
20  character(100) :: jobname
21  character(100) :: bml_type
22  real(dp) :: ea
23  real(dp) :: eb
24  real(dp) :: dab
25  real(dp) :: daiaj
26  real(dp) :: dbibj
27  real(dp) :: dec, rcoeff
28  logical :: reshuffle
29  end type mham_type
30 
32 
33 contains
34 
36  subroutine prg_parse_mham(mham,filename)
37 
38  implicit none
39  type(mham_type), intent(inout) :: mham
40  integer, parameter :: nkey_char = 2, nkey_int = 2, nkey_re = 7, nkey_log = 1
41  character(len=*) :: filename
42 
43  !Library of keywords with the respective defaults.
44  character(len=50), parameter :: keyvector_char(nkey_char) = [character(len=50) :: &
45  'JobName=', 'BMLType=']
46  character(len=100) :: valvector_char(nkey_char) = [character(len=100) :: &
47  'GetModelHam', 'Dense' ]
48 
49  character(len=50), parameter :: keyvector_int(nkey_int) = [character(len=50) :: &
50  'NOrbs=', 'Seed=']
51  integer :: valvector_int(nkey_int) = (/ &
52  10, 100 /)
53 
54  character(len=50), parameter :: keyvector_re(nkey_re) = [character(len=50) :: &
55  'EpsilonA=', 'EpsilonB=', 'DeltaAB=','DeltaAiAj=','DeltaBiBj=','Decay=','RCoeff=']
56  real(dp) :: valvector_re(nkey_re) = (/&
57  0.0, 0.0, -1.0, 0.0, -1.0, -100.0, 0.0 /)
58 
59  character(len=50), parameter :: keyvector_log(nkey_log) = [character(len=50) :: &
60  'Reshuffle=']
61  logical :: valvector_log(nkey_log) = (/&
62  .false./)
63 
64  !Start and stop characters
65  character(len=50), parameter :: startstop(2) = [character(len=50) :: &
66  'MHAM{', '}']
67 
68  call prg_parsing_kernel(keyvector_char,valvector_char&
69  ,keyvector_int,valvector_int,keyvector_re,valvector_re,&
70  keyvector_log,valvector_log,trim(filename),startstop)
71 
72  !Characters
73  mham%JobName = valvector_char(1)
74 
75  if(valvector_char(2) == "Dense")then
76  mham%bml_type = bml_matrix_dense
77  elseif(valvector_char(2) == "Ellpack")then
78  mham%bml_type = bml_matrix_ellpack
79  elseif(valvector_char(2) == "CSR")then
80  mham%bml_type = bml_matrix_csr
81  elseif(valvector_char(2) == "Ellblock")then
82  mham%bml_type = bml_matrix_ellblock
83  endif
84 
85  !Integers
86  mham%norbs = valvector_int(1)
87  mham%seed = valvector_int(2)
88 
89  !Reals
90  mham%ea = valvector_re(1)
91  mham%eb = valvector_re(2)
92  mham%dab = valvector_re(3)
93  mham%daiaj = valvector_re(4)
94  mham%dbibj = valvector_re(5)
95  mham%dec = valvector_re(6)
96  mham%rcoeff = valvector_re(7)
97 
98  !Logicals
99  mham%reshuffle = valvector_log(1)
100 
101  end subroutine prg_parse_mham
102 
116  subroutine prg_twolevel_model(ea, eb, dab, daiaj, dbibj, dec, rcoeff, reshuffle, &
117  & seed, h_bml, verbose)
118  real(dp), intent(in) :: ea, eb, dab, daiaj, dbibj, rcoeff
119  integer, intent(in) :: verbose, seed
120  integer, allocatable :: seedin(:)
121  logical, intent(in) :: reshuffle
122  type(bml_matrix_t),intent(inout) :: h_bml
123  real(dp), allocatable :: diagonal(:), row(:), rowi(:), rowj(:)
124  type(bml_matrix_t) :: ht_bml
125  integer :: norbs, i, j, ssize
126  real(dp) :: dec, dist, ran
127 
128  norbs = bml_get_n(h_bml)
129  allocate(diagonal(norbs))
130  allocate(row(norbs))
131 
132  call random_seed()
133  call random_seed(size=ssize)
134  allocate(seedin(ssize))
135  seedin = seed
136  call random_seed(put=seedin)
137 
138  do i=1,norbs
139  if(mod(i,2) == 0)then
140  call random_number(ran)
141  diagonal(i) = ea + rcoeff*(2.0_dp*ran - 1.0_dp)
142  else
143  call random_number(ran)
144  diagonal(i) = eb + rcoeff*(2.0_dp*ran - 1.0_dp)
145  endif
146  enddo
147 
148  do i=1,norbs
149  do j=1,norbs
150  if(abs(real(i-j,dp)) <= norbs/2.0d0) then
151  dist = max(abs(real(i-j,dp))- 2.0_dp,0.0_dp)
152  else
153  dist = max((-abs(real(i-j,dp))+norbs) - 2.0_dp,0.0_dp)
154  endif
155  !A-A type
156  if((mod(i,2) .ne. 0) .and. (mod(j,2) .ne. 0))then
157  call random_number(ran)
158  row(j) = (daiaj + rcoeff*(2.0_dp*ran - 1.0_dp))*exp(dec*dist)
159  !A-B type
160  elseif((mod(i,2) == 0) .and. (mod(j,2) == 0))then
161  call random_number(ran)
162  row(j) = (dbibj + rcoeff*(2.0_dp*ran - 1.0_dp))*exp(dec*dist)
163  !B-B type
164  else
165  call random_number(ran)
166  row(j) = (dab + rcoeff*(2.0_dp*ran - 1.0_dp))*exp(dec*dist)
167  endif
168  ! write(*,*)i,j,row(j),mod(i,2),mod(j,2),abs(real(i-j,dp)+norbs),abs(real(i-j,dp)+norbs)
169 
170  enddo
171  call bml_set_row(h_bml,i,row)
172  enddo
173 
174  call bml_set_diagonal(h_bml,diagonal)
175 
176  !Symmetrization
177  call bml_copy_new(h_bml,ht_bml)
178  call bml_transpose_new(h_bml,ht_bml)
179  if(verbose.gt.0)then
180  call bml_print_matrix("h_bml",h_bml,0,10,0,10)
181  call bml_print_matrix("ht_bml",ht_bml,0,10,0,10)
182  endif
183  call bml_add(h_bml,ht_bml,0.5d0,0.5d0,0.0d0)
184 
185  call bml_deallocate(ht_bml)
186 
187  if(reshuffle)then
188  allocate(rowj(norbs))
189  allocate(rowi(norbs))
190  do i=1,norbs
191  call random_number(ran)
192  j = int(floor(ran*norbs+1))
193  call bml_get_row(h_bml,i,rowi)
194  call bml_get_row(h_bml,j,rowj)
195  call bml_set_row(h_bml,i,rowj)
196  call bml_set_row(h_bml,j,rowi)
197  enddo
198  deallocate(rowi)
199  deallocate(rowj)
200  endif
201 
202  end subroutine prg_twolevel_model
203 
217  subroutine prg_twolevel_model3d(ea, eb, dab, daiaj, dbibj, dec, rcoeff, reshuffle, &
218  & seed, h_bml, verbose)
219  real(dp), intent(in) :: ea, eb, dab, daiaj, dbibj, rcoeff
220  integer, intent(in) :: verbose, seed
221  integer, allocatable :: seedin(:)
222  logical, intent(in) :: reshuffle
223  type(bml_matrix_t),intent(inout) :: h_bml
224  real(dp), allocatable :: diagonal(:), row1(:), row2(:), rowi(:), rowj(:)
225  type(bml_matrix_t) :: ht_bml
226  integer :: norbs, i, j, ssize, i1, j1, k1, i2, j2, k2, n1d
227  real(dp) :: dec, dist, di, dj, dk, ran, tmp
228 
229  norbs = bml_get_n(h_bml)
230 
231  allocate(diagonal(norbs))
232  allocate(row1(norbs))
233  allocate(row2(norbs))
234 
235  tmp = norbs
236  tmp = tmp*0.5
237  tmp = tmp**(1./3.)
238  n1d = int(tmp)
239  print*,'n1d = ',n1d
240 
241  call random_seed()
242  call random_seed(size=ssize)
243  allocate(seedin(ssize))
244  seedin = seed
245  call random_seed(put=seedin)
246 
247  do i=1,norbs
248  if(mod(i,2) == 0)then
249  call random_number(ran)
250  diagonal(i) = ea + rcoeff*(2.0_dp*ran - 1.0_dp)
251  else
252  call random_number(ran)
253  diagonal(i) = eb + rcoeff*(2.0_dp*ran - 1.0_dp)
254  endif
255  enddo
256 
257  do i1=1,n1d
258  do j1=1,n1d
259  do k1=1,n1d
260  i = (i1-1) + n1d*(j1-1) + n1d*n1d*(k1-1)
261  do i2=1,n1d
262  do j2=1,n1d
263  do k2=1,n1d
264  j = (i2-1) + n1d*(j2-1) + n1d*n1d*(k2-1)
265  if(abs(real(i1-i2,dp)) <= n1d/2.0d0) then
266  di = max(abs(real(i1-i2,dp))- 2.0_dp,0.0_dp)
267  else
268  di = max((-abs(real(i1-i2,dp))+n1d)- 2.0_dp,0.0_dp)
269  endif
270  if(abs(real(j1-j2,dp)) <= n1d/2.0d0) then
271  dj = max(abs(real(j1-j2,dp))- 2.0_dp,0.0_dp)
272  else
273  dj = max((-abs(real(j1-j2,dp))+n1d)- 2.0_dp,0.0_dp)
274  endif
275  if(abs(real(k1-k2,dp)) <= n1d/2.0d0) then
276  dk = max(abs(real(k1-k2,dp))- 2.0_dp,0.0_dp)
277  else
278  dk = max((-abs(real(k1-k2,dp))+n1d)- 2.0_dp,0.0_dp)
279  endif
280  dist = dsqrt(di+di+dj*dj+dk*dk)
281 
282  !assign matrix elements for all interactions between two atoms
283  !A-A type
284  call random_number(ran)
285  row1(2*j+1) = (daiaj + rcoeff*(2.0_dp*ran - 1.0_dp))*exp(dec*dist)
286  !A-B type
287  call random_number(ran)
288  row1(2*j+2) = (dab + rcoeff*(2.0_dp*ran - 1.0_dp))*exp(dec*dist)
289  row2(2*j+1) = (dab + rcoeff*(2.0_dp*ran - 1.0_dp))*exp(dec*dist)
290  !B-B type
291  call random_number(ran)
292  row2(2*j+2) = (dbibj + rcoeff*(2.0_dp*ran - 1.0_dp))*exp(dec*dist)
293  ! write(*,*)i,j,row(j),mod(i,2),mod(j,2),abs(real(i-j,dp)+norbs),abs(real(i-j,dp)+norbs)
294  enddo
295  enddo
296  enddo
297  call bml_set_row(h_bml,2*i+1,row1)
298  call bml_set_row(h_bml,2*i+2,row2)
299  enddo
300  enddo
301  enddo
302 
303  call bml_set_diagonal(h_bml,diagonal)
304 
305  !Symmetrization (necessary when random factor used)
306  call bml_copy_new(h_bml,ht_bml)
307  call bml_transpose_new(h_bml,ht_bml)
308  if(verbose.gt.0)then
309  call bml_print_matrix("h_bml",h_bml,0,10,0,10)
310  call bml_print_matrix("ht_bml",ht_bml,0,10,0,10)
311  endif
312  call bml_add(h_bml,ht_bml,0.5d0,0.5d0,0.0d0)
313 
314  call bml_deallocate(ht_bml)
315 
316  if(reshuffle)then
317  allocate(rowj(norbs))
318  allocate(rowi(norbs))
319  do i=1,norbs
320  call random_number(ran)
321  j = int(floor(ran*norbs+1))
322  call bml_get_row(h_bml,i,rowi)
323  call bml_get_row(h_bml,j,rowj)
324  call bml_set_row(h_bml,i,rowj)
325  call bml_set_row(h_bml,j,rowi)
326  enddo
327  deallocate(rowi)
328  deallocate(rowj)
329  endif
330 
331  end subroutine prg_twolevel_model3d
332 
333 end module prg_modelham_mod
Some general parsing functions.
subroutine, public prg_parsing_kernel(keyvector_char, valvector_char, keyvector_int, valvector_int, keyvector_re, valvector_re, keyvector_log, valvector_log, filename, startstop)
The general parsing function. It is used to vectorize a set of "keywords" "value" pairs as included i...
The prg_hamiltonian module.
subroutine, public prg_parse_mham(mham, filename)
Model Ham parse.
subroutine, public prg_twolevel_model(ea, eb, dab, daiaj, dbibj, dec, rcoeff, reshuffle, seed, h_bml, verbose)
Construct a two-level model Hamiltonian.
subroutine, public prg_twolevel_model3d(ea, eb, dab, daiaj, dbibj, dec, rcoeff, reshuffle, seed, h_bml, verbose)
Construct a two-level model Hamiltonian.
A module to read and handle chemical systems.
General ModelHam type.