Adding reshuffle in rows (#150)

* Adding reshuffle in rows

* fix

examples/twolevelmodel/input.in

@@ -1,13 +1,15 @@
 MHAM{
-  NOrbs= 1000
-  BMLType= Ellpack
-  EpsilonA= -10.0
-  EpsilonB= -0.0
-  DeltaAB= -1.0
+  NOrbs= 10
+  BMLType= Dense
+  EpsilonA= -1.0
+  EpsilonB= -2.0
+  DeltaAB= -0.0
   DeltaAiAj= -0.0
-  DeltaBiBj= -1.0
-  Decay= -1.0
+  DeltaBiBj= -0.0
+  Decay= -1000.01
   RCoeff= 0.0
+  Seed= 400
+  Reshuffle= T
 }
 
 

examples/twolevelmodel/input.in.metal

@@ -7,4 +7,6 @@ MHAM{
   DeltaAiAj= -1.0
   DeltaBiBj= -1.0
   Decay= -0.01
+  Seed= 100
+  Reshuffle= F
 }

examples/twolevelmodel/input.in.semicond

@@ -8,6 +8,8 @@ MHAM{
   DeltaBiBj= -1.0
   Decay= -1000.0
   RCoeff= 0.0
+  Seed= 100
+  Reshuffle= F
 }
 
 

examples/twolevelmodel/input.in.softmatt

@@ -8,6 +8,8 @@ MHAM{
   DeltaBiBj= -1.0
   Decay= -1.0
   RCoeff= 1.0
+  Seed= 100
+  Reshuffle= F
 }
 
 

examples/twolevelmodel/main.F90

@@ -45,7 +45,7 @@ program hmodel
 
   !Constructng the Hamiltonian
   call prg_twolevel_model(mham%ea, mham%eb, mham%dab, mham%daiaj, mham%dbibj, &
-  &mham%dec, mham%rcoeff, ham_bml, verbose)
+  &mham%dec, mham%rcoeff, mham%reshuffle, mham%seed, ham_bml, verbose)
 
   call bml_print_matrix("ham_bml",ham_bml,0,10,0,10)
 

src/prg_modelham_mod.F90

@@ -16,7 +16,7 @@ module prg_modelham_mod
 
   !> General ModelHam type
   type, public :: mham_type
-    integer :: norbs
+    integer :: norbs, seed
     character(100) :: jobname
     character(100) :: bml_type
     real(dp) :: ea
@@ -25,7 +25,7 @@ module prg_modelham_mod
     real(dp) :: daiaj
     real(dp) :: dbibj
     real(dp) :: dec, rcoeff    
-    logical :: log1
+    logical :: reshuffle
   end type mham_type  
 
   public :: prg_parse_mham, prg_twolevel_model
@@ -37,7 +37,7 @@ subroutine prg_parse_mham(mham,filename)
 
     implicit none
     type(mham_type), intent(inout) :: mham
-    integer, parameter :: nkey_char = 2, nkey_int = 1, nkey_re = 7, nkey_log = 1
+    integer, parameter :: nkey_char = 2, nkey_int = 2, nkey_re = 7, nkey_log = 1
     character(len=*) :: filename
 
     !Library of keywords with the respective defaults.
@@ -47,17 +47,17 @@ subroutine prg_parse_mham(mham,filename)
          'GetModelHam', 'Dense' ]
 
     character(len=50), parameter :: keyvector_int(nkey_int) = [character(len=50) :: &
-         'NOrbs=']
+         'NOrbs=', 'Seed=']
     integer :: valvector_int(nkey_int) = (/ &
-         10  /)
+         10, 100  /)
 
     character(len=50), parameter :: keyvector_re(nkey_re) = [character(len=50) :: &
          'EpsilonA=', 'EpsilonB=', 'DeltaAB=','DeltaAiAj=','DeltaBiBj=','Decay=','RCoeff=']
     real(dp) :: valvector_re(nkey_re) = (/&
          0.0, 0.0, -1.0, 0.0, -1.0, -100.0, 0.0 /)
 
     character(len=50), parameter :: keyvector_log(nkey_log) = [character(len=50) :: &
-         'Dummy=']
+         'Reshuffle=']
     logical :: valvector_log(nkey_log) = (/&
          .false./)
 
@@ -80,6 +80,7 @@ subroutine prg_parse_mham(mham,filename)
 
     !Integers
     mham%norbs = valvector_int(1)
+    mham%seed = valvector_int(2)
 
     !Reals
     mham%ea = valvector_re(1)
@@ -91,7 +92,7 @@ subroutine prg_parse_mham(mham,filename)
     mham%rcoeff = valvector_re(7)
 
     !Logicals
-    mham%log1 = valvector_log(1)
+    mham%reshuffle = valvector_log(1)
 
   end subroutine prg_parse_mham
 
@@ -104,15 +105,20 @@ end subroutine prg_parse_mham
   !! \param dbibj Intersite second level Hamiltonian elements
   !! \param dec Decay constant
   !! \param rcoeff Random coefficient
+  !! \param reshuffle If rows needs to be reshuffled
+  !! \param seed Random seed
   !! \param h_bml Output hamiltonian matrix 
   !! \param verbose Verbosity level
-  subroutine prg_twolevel_model(ea, eb, dab, daiaj, dbibj, dec, rcoeff, h_bml, verbose)
+  subroutine prg_twolevel_model(ea, eb, dab, daiaj, dbibj, dec, rcoeff, reshuffle, &
+   & seed, h_bml, verbose)
     real(dp), intent(in) :: ea, eb, dab, daiaj, dbibj, rcoeff
     integer, intent(in) :: verbose
+    integer :: seed
+    logical, intent(in) :: reshuffle
     type(bml_matrix_t),intent(inout) :: h_bml
+    real(dp), allocatable :: diagonal(:), row(:), rowi(:), rowj(:)
     type(bml_matrix_t) :: ht_bml
-    integer :: norbs, i, j, seed
-    real(dp), allocatable :: diagonal(:), row(:)
+    integer :: norbs, i, j
     real(dp) :: dec, dist, ran
 
     norbs = bml_get_N(h_bml)
@@ -158,14 +164,29 @@ subroutine prg_twolevel_model(ea, eb, dab, daiaj, dbibj, dec, rcoeff, h_bml, ver
     enddo
 
     call bml_set_diagonal(h_bml,diagonal)
-    
+
     !Symmetrization
     call bml_copy_new(h_bml,ht_bml)
     call bml_transpose(h_bml,ht_bml)
     call bml_print_matrix("h_bml",h_bml,0,10,0,10)
     call bml_print_matrix("ht_bml",ht_bml,0,10,0,10)
     call bml_add(h_bml,ht_bml,0.5d0,0.5d0,0.0d0)
 
+    if(reshuffle)then
+      allocate(rowj(norbs))
+      allocate(rowi(norbs))
+      do i=1,norbs 
+        call random_number(ran)        
+        j = int(floor(ran*norbs+1))
+        call bml_get_row(h_bml,i,rowi)
+        call bml_get_row(h_bml,j,rowj)
+        call bml_set_row(h_bml,i,rowj)
+        call bml_set_row(h_bml,j,rowi)
+      enddo
+      deallocate(rowi)
+      deallocate(rowj)
+    endif  
+
   end subroutine prg_twolevel_model
 
 end module prg_modelham_mod

src/prg_sp2_mod.F90

@@ -796,20 +796,31 @@ subroutine prg_sp2_alg1(h_bml, rho_bml, threshold, bndfil, minsp2iter, maxsp2ite
        limdiff = abs(trx - trx2 - occ) - abs(trx + trx2 - occ)
 
        if (limdiff .ge. idemtol) then
-
+          write(*,*)"H1"
           ! X <- X + (X - X * X) <- 2 * X - X * X
           call bml_add_deprecated(1.0_dp, rho_bml, 1.0_dp, x2_bml, threshold)
 
           trx = trx + trx2
 
        elseif(limdiff .lt. -idemtol) then
+          write(*,*)"H2"
 
           ! X <- X - (X - X * X) <- X * X
           call bml_add_deprecated(1.0_dp, rho_bml, -1.0_dp, x2_bml, threshold)
 
           trx = trx - trx2
+
+       elseif((limdiff .eq. 0.0) .and. (iter .eq. 1)) then 
+          write(*,*)"H3"
+
+         ! X <- X - (X - X * X) <- X * X
+          call bml_add_deprecated(1.0_dp, rho_bml, -1.0_dp, x2_bml, threshold)
 
-       else
+          trx = trx - trx2
+
+       else  
+          write(*,*)"H4"
+          write(*,*)"limdiff,idemtol",limdiff, idemtol
 
           iter = iter - 1
           breakloop = 1
@@ -819,11 +830,14 @@ subroutine prg_sp2_alg1(h_bml, rho_bml, threshold, bndfil, minsp2iter, maxsp2ite
        idemperr2 = idemperr1
        idemperr1 = idemperr
        idemperr = abs(trx2)
-
-       if (sp2conv .eq. "Rel" .and. iter .ge. minsp2iter .and. &
+      write(*,*)sp2conv,iter,minsp2iter,idemtol
+       if (iter .ge. minsp2iter) then 
+          write(*,*)"sssss",iter,minsp2iter
+          if (sp2conv .eq. "Rel" .and. &
             (idemperr2 .le. idemperr .or. idemperr .lt. idemtol)) then
-          breakloop = 1
-       end if
+            breakloop = 1
+          end if
+       end if   
 
        if (iter .eq. maxsp2iter) then
           write(*,*) "SP2 purification is not converging: STOP!"

tests/src/main.F90

@@ -1057,7 +1057,7 @@ program main
      call bml_zero_matrix(mham%bml_type,bml_element_real,dp,mham%norbs,mham%norbs,aux_bml)
      call bml_zero_matrix(mham%bml_type,bml_element_real,dp,mham%norbs,mham%norbs,ham_bml)   
      call prg_twolevel_model(mham%ea, mham%eb, mham%dab, mham%daiaj, mham%dbibj, &
-            &mham%dec, mham%rcoeff, ham_bml, verbose)
+            &mham%dec, mham%rcoeff, mham%reshuffle, mham%seed, ham_bml, verbose)
      call bml_read_matrix(aux_bml,'hamiltonian-twolevel-ref.mtx')
      call bml_add_deprecated(-1.0_dp,aux_bml,1.0_dp,ham_bml,0.0_dp)
      error_calc = bml_fnorm(aux_bml)