Merge lp:~nickpapior/siesta/4.1-init-DM into lp:siesta/4.1

[Nick Papior]

Added a new initialization method for DM

The current DM initialization *only* specifies the charge
on the truly diagonal DM entries. I.e. io == col(ind).
However, if one considers that atomic charges are actually
atomic *states* i.e. a state-vector with element in the
diagonal one will also get DM entries in off-diagonal
but io - 1 == col(ind) % no_u which is probably more correct.

I have added two new methods for DM initialization:

1. The DM.Init flag determines one of these methods:
   a) atomic
      the default on-site only atomic density
   b) expand to a state with on-site density
      This enables small auxiliary cell calculations
      to have the density spread amongst overlapping
      orbitals.
      In some systems this seem to behave better
      (graphite/graphene) whereas others it doesn't
      (si64).

2. Added DM.Init.RandomStates <integer> to remove
   N electrons from the initial density matrix
   and add N electrons using a randomized state.
   This state is ensured to be in the non-orthogonal
   basis but is not necessarily orthogonal to any
   other state generated (nor will it be orthogonal
   to the atomic states).
   Currently this method scales orbital coefficients with
   the atomic charges. This will match the atomic charges to
   some degree.
   This is done to limit randoized states with charge in
   high-lying shells which seems dubious in many cases.
   Generally this does not work too good, but may be good
   in some cases I don't know off.

Currently the defaults are as they have always been, but
perhaps the atomic-state would be good in the future?
This needs more heuristics.

I have tried atomic vs. atomic-state on the following tests:

The numbers shows the number of SCF steps needed to converge and the first step dDmax, dHmax

graphite_c6:
atomic: 9, 1.879391, 1.748376
atomic-state: 7, 1.568393, 0.207440

batio3:
atomic: 19, 1.961004, 22.268878
atomic-state: 21, 2.144922, 20.791337

sic-slab:
COMMENT: fun-fact, the SCF cycle indicates that atomic-states is better, but then suddenly atomic converges faster at iSCF 46! History has lots to say here since there is an electric field in this example.
atomic: 69, 1.746005, 13.382736
atomic-state: 144, 2.007246, 13.363516

nanotube-c-5-0:
atomic: 11, 0.741924, 4.242417
atomic-state: 11, 2.037193, 4.270371

[Alberto Garcia (albertog)]

I am trying to understand this a bit better.

The initial DM is used to compute the first H, but *only* the charge density computed from that DM is relevant (H0 does not use the DM, and dhscf only uses rho(r)). The initial charge density then
is the key for initialization. Currently the heuristic is that a superposition of atomic charges is an appropriate method. To improve on that we would need some extra heuristics about how bonding effects in that particular system might influence the charge density. But this cannot be guessed simply from manipulation of 'atomic' concepts.

Does your alternate method give as starting rho(r) the same atomic-charges superposition? I guess not, since the convergence results are different. So you are changing the charge heuristic, but not in a physical way, so I think that the results will be erratic and depend very much on whether your trick leads to a "closer" rho(r) for a particular system.

Besides this, the initial DM will indeed have an influence on the delta_DM found after the first scf step. For example, the 'diagonal' method disregards the nzeta-1 orbitals beyond the first in DM_in, and these are likely to be occupied in DM_out. But delta_DM(scf=1) should be discarded anyway if we are not mixing in the first scf step, as is the case when using atomic initialization.

There are other possibilites for a more targeted initialization: one could extract from a set of calculations in the same or similar system some kind of "deformation charge" estimate and add it to the initial charge density. Note that this is already implemented, but not really well explored.

Or a "physically motivated" initial H could be read. Maybe an extended-Hückel form.

[Nick Papior (nickpapior)]

1. Agreed on your first paragraph.

2. Regarding the 2nd paragraph.
The *diagonal* Currently the initial DM is made from *only* diagonal elements in the DM. I.e. for supercell systems where orbital io connects with mirror io only io, io is initialized.
For my *state* option the DM is populated as though the identity matrix is the set of eigenvectors (rescaled for the non-orthogonal basis). And thus in this case one will find DM values in both io, io, AND io, mirror io.
My main motivation was that an eigenstate with {1,0,0,0,...} would expand to also mirror io couplings as for calculations without auxiliary supercells.
Do you think this is *very* unphysical? If so then non-auxiliary supercells should always be turned off for diagonal indices since this leads to non atomic initial DM charges. This is substantiated by the sic-slab calculation performed in 3 different settings (in all cases the SCF has only the Gamma-point):
diagonal-forced-auxiliary: 83, 1.820250, 13.346312
diagonal-no-auxiliary: 96, 1.767353, 13.326340
state-forced-auxiliary: 139, 1.767353, 13.326336

Here it can be seen that the initial DM change for the diagonal-no-auxiliary vs. state-forced-auxiliary are the same (as expected via the above).

3. This was my *first* and most basic attempt to try and investigate implications of initial DM. As I have told you previously this could be moved to sisl (I have tried this with success) to do more advanced things quicker, but for now I thought that these would be nice to have for certain cases.
Note that my implementation for *randomized* states was inspired by PW codes which does this (apparently with great success?) and hence I thought it could be used as well here. However, they are much worse, but sometimes using 1 randomized state helps...

4. Agreed that some better heuristic and more thorough investigations could be needed here.

[Alberto Garcia (albertog)]

I was trying to frame the issue in a more general setting to understand it better.
And my "unphysical" label refers just to not including system-specific 'bonding' effects.

From what you say I might have missed an important more concrete point: is the current
algorithm not giving atomic-charge superpositions in all cases?
I thought we had fixed this (in a zeroth order approximation) in rev 576 of 4.0, which deals with the (arguably pathological) case of severe folding (diagonal images) and gamma-point.

Are there any other cases in which rho(r)[DM_initial] is not a perfect superposition of atomic charges?

[Nick Papior (nickpapior)]

1. Ok :)

2. No, the current algorithm will *not* give atomic-only-charges IFF 1) the auxiliary is not used AND 2) if the cell is so small that orbitals have interactions with its own mirror orbitals (diagonal images).

If we always *only* want atomic-only charges (i.e. no supercell image overlap charges), then my state solution should not go in.

3. Let me however note that rarely would people do NO k-point sampling for cells needing an auxiliary supercell.

This is probably also a notice for us to consider why we get different results (in SCF cycles) for small cells with/without auxiliary supercells.

[lp:~nickpapior/siesta/4.1-init-DM]

1114. By Nick Papior on 2019-09-19

Removed the atomic-state method which should never be used

Retained the randomized initialization for experimental
use.

[Nick Papior (nickpapior)]

Merged on GitLab.

=== modified file 'Docs/siesta.tex'
--- Docs/siesta.tex	2019-08-22 07:58:14 +0000
+++ Docs/siesta.tex	2019-09-19 13:04:58 +0000
@@ -5149,6 +5149,25 @@
   
 \end{fdflogicalF}
 
+\begin{fdfentry}{DM!Init}<atomic>
+  \index{spin!initialization}%
+
+  Specify the initial density matrix composition. All these methods
+  also obey \fdf{DM!InitSpin!AF}.
+
+  \begin{fdfoptions}
+    
+    \option[atomic]%
+    \fdfindex*{DM!Init!atomic}
+
+    Only initialize the diagonal (on-site) density matrices according
+    to the initial charges on the atomic orbital.
+
+  \end{fdfoptions}
+
+\end{fdfentry}
+
+
 \begin{fdflogicalF}{DM!InitSpin!AF}
   \index{spin!initialization}%
   \index{ferromagnetic initial DM}%
@@ -5205,6 +5224,28 @@
 
 \end{fdfentry}
 
+\begin{fdfentry}{DM!Init.RandomStates}[integer]<0>
+
+  Randomize $N$ number of states in the initial density-matrix. This
+  may be useful to give the system an initial kick. The randomized
+  states only has coefficients where the initial atomic charges are
+  non-zero. This prohibits initial charges in high-lying shells.
+
+  This only has an effect if the density matrix is started from an
+  atomic density and/or when using \fdf{DM!InitSpin}.
+
+  In case it is used together with \fdf{DM!InitSpin} it also
+  randomizes the spin-configuration which may be undesirable.
+
+  \note this is currently experimental since the randomized states are
+  not ensured orthogonal. This flag may be removed in later
+  revisions with some super seeding functions. If testing this, start
+  with a value of $1$ to see if it has an effect, any higher numbers
+  will probably be worse.
+
+\end{fdfentry}
+
+
 \begin{fdflogicalT}{DM!AllowReuse}
 
   Controls whether density matrix information from previous geometry

=== modified file 'Src/m_new_dm.F90'
--- Src/m_new_dm.F90	2019-08-22 07:58:14 +0000
+++ Src/m_new_dm.F90	2019-09-19 13:04:58 +0000
@@ -935,6 +935,8 @@
 
     ! *** Local variables
     type(block_fdf) :: bfdf
+    character(len=32) :: method
+    integer :: i_method ! 1 == atomic
     logical :: init_spin_block
     type(dData2D) :: arr_2D
     ! The pointers for the sparse pattern
@@ -942,14 +944,21 @@
     integer :: io, gio, i, ind, jo
     integer, pointer :: ncol(:), ptr(:), col(:)
     real(dp), pointer :: DM(:,:)
+
+    method = fdf_get('DM.Init', 'atomic')
+    if ( leqi(method, 'atomic') ) then
+      i_method = 1
+    else
+      call die('m_new_dm: erroneous option: DM.Init [atomic]')
+    end if
     
     ! First we need to figure out the options
     if ( spin%DM > 1 ) then
-       init_spin_block = fdf_block('DM.InitSpin', bfdf)
+      init_spin_block = fdf_block('DM.InitSpin', bfdf)
     else
-       ! Do not read the initspin block in case of
-       ! non-polarized calculation
-       init_spin_block = .false.
+      ! Do not read the initspin block in case of
+      ! non-polarized calculation
+      init_spin_block = .false.
     end if
 
     ! Retrieve pointers to data
@@ -963,10 +972,10 @@
     call delete(arr_2D)
     DM => val(DM_2D)
 
-    if ( .not. init_spin_block ) then
-       call init_atomic()
-    else
-       call init_user()
+    if ( init_spin_block ) then
+      call init_user()
+    else if ( i_method == 1 ) then
+      call init_atomic()
     end if
 
     ! We have initialized with atomic information. Correct in case we
@@ -975,18 +984,27 @@
     ! Gamma-point only. Otherwise S(diagonal) is always 1.0 and the
     ! simple atomic-orbital superposition works as intended.
 
+    if ( init_spin_block .or. i_method == 1 ) then
+      
       do io = 1, no_l
-         ! Retrieve global orbital index
-         gio = index_local_to_global(dit, io)
-         do i = 1 , ncol(io)
-            ind = ptr(io) + i
-            jo = col(ind)
-            ! Check for diagonal element
-            if ( gio == jo ) then
-               DM(ind,:) = DM(ind,:) / S(ind)
-            endif
-         end do
+        ! Retrieve global orbital index
+        gio = index_local_to_global(dit, io)
+        do i = 1 , ncol(io)
+          ind = ptr(io) + i
+          jo = col(ind)
+          ! Check for diagonal element
+          if ( gio == jo ) then
+            DM(ind,:) = DM(ind,:) / S(ind)
+          endif
+        end do
       end do
+     
+    end if
+      
+    ! Check if we should randomize some states
+    ! The default is *no* states are randomized
+    ! Also, this method does not obey anti-ferro!
+    call init_randomize()
 
   contains
 
@@ -1346,9 +1364,116 @@
 
     end subroutine init_user
 
+    subroutine init_randomize()
+
+      use parallel, only: Node
+      use geom_helper, only: ucorb
+#ifdef MPI
+      use mpi_siesta
+#endif
+
+      integer :: io, gio, i, ind, jo, i1, i2
+      integer :: n_random, ir
+      real(dp), allocatable :: rs(:)
+      real(dp) :: qtot, qinit, scale
+#ifdef MPI
+      integer :: ierr
+#endif
+      
+      n_random = fdf_get('DM.Init.RandomStates', 0)
+
+      ! If the user does not request any randomization, simply return
+      if ( n_random == 0 ) return
+
+      ! Calculate the number of charges in the current DM (this will be our target!)
+      qtot = 0._dp
+      do i2 = 1 , spin%DM
+        do i1 = 1 , nnz
+          qtot = qtot + DM(i1,i2) * S(i1)
+        end do
+      end do
+#ifdef MPI
+      call MPI_AllReduce(qtot, qinit, 1, mpi_double_precision, MPI_SUM, &
+          MPI_COMM_WORLD, ierr)
+      qtot = qinit
+#endif
+
+      ! Truncate number of random states (NOT nint)
+      if ( qtot < n_random ) n_random = int(qtot)
+
+      ! Calculate the target *initial DM*
+      qinit = qtot - n_random
+      scale = qinit / qtot
+
+      ! Rescale DM to remove n_random electrons
+      do i2 = 1 , spin%spinor
+        do i1 = 1 , nnz
+          DM(i1,i2) = DM(i1,i2) * scale
+        end do
+      end do
+      
+      ! Allocate for random states
+      ! TODO allocate all states and make them orthogonal
+      allocate(rs(no_u))
+
+      do ir = 1, n_random
+
+        if ( Node == 0 ) then
+          call random_number(rs)
+          ! Scale with atomic weights to not have *random* states
+          ! on, say, f-shells ;)
+          do io = 1 , no_u
+            rs(io) = (rs(io) - 0.5_dp) * DM_atom(io)
+          end do
+        end if
+#ifdef MPI
+        call MPI_Bcast(rs,no_u,MPI_Double_Precision,0,Mpi_Comm_World,ierr)
+#endif
+
+        ! Calculate inner product to rescale for 1 electron for this state
+        scale = 0._dp
+        do io = 1, no_l
+          gio = index_local_to_global(dit, io)
+          do ind = ptr(io) + 1 , ptr(io) + ncol(io)
+            jo = ucorb(col(ind), no_u)
+            scale = scale + rs(gio) * S(ind) * rs(jo)
+          end do
+        end do
+#ifdef MPI
+        call MPI_AllReduce(scale, qinit, 1, mpi_double_precision, MPI_SUM, &
+            MPI_COMM_WORLD, ierr)
+        scale = qinit
+#endif
+        ! Inner product value for getting a charge of 1
+        scale = sqrt(scale)
+        if ( spin%spinor > 1 ) then
+          ! Same state for spin-up/down
+          rs(:) = rs(:) / (scale * 2)
+        else
+          rs(:) = rs(:) / scale
+        end if
+
+        do io = 1, no_l
+          gio = index_local_to_global(dit, io)
+          do ind = ptr(io) + 1 , ptr(io) + ncol(io)
+            jo = ucorb(col(ind), no_u)
+            DM(ind,1:spin%spinor) = DM(ind,1:spin%spinor) + rs(gio) * rs(jo)
+          end do
+        end do
+      end do
+
+      ! Clean-up
+      deallocate(rs)
+
+      if ( IONode ) then
+        write(*,'(a,i0,a)') 'DM randomized ', n_random, ' states'
+      end if
+      
+    end subroutine init_randomize
+
   end subroutine init_DM_atomic
 
-  
+
   subroutine extrapolate_dm_with_coords(DM_history,na_u,xa,sparse_pattern,DMnew)
     use class_Sparsity
     use class_dData2D

=== modified file 'version.info'
--- version.info	2019-09-19 12:58:46 +0000
+++ version.info	2019-09-19 13:04:58 +0000
@@ -1,1 +1,5 @@
+<<<<<<< TREE
 siesta-4.1-1120
+=======
+siesta-4.1-1112--init-DM-2
+>>>>>>> MERGE-SOURCE