Merge lp:~nickpapior/siesta/4.1-scf-first into lp:siesta/4.1

[Nick Papior]

Fixing SCF-first for re-reading DM files

Previously one could not restart and perform mixing for the
first SCF step. This helps for big systems when needing restart
effects.

In some cases (for really large systems) it *could* be advantegeous
to also allow mixing in the first SCF step since any perturbation
cause huge dDmax/dHmax. In those cases we probably need some kind
of handler to allow mixing first step anyhow.

  SCF.Mix.First.MD <T|F> ?

This is a result of https://bugs.launchpad.net/siesta/+bug/1834949 where restart is not working properly or large systems.

Initially we had disabled scf.mix.first for all simulations which did a restart. This commit fixes this in the sense that iff the sparsity pattern is exactly the same, it will allow one to do mixing in the first iteration, otherwise not.

However, this *could* potentially still be a problem for very large systems where minor displacements are present. In that case it could still be viable to do small linear mixing for a few steps using the final DM from a different geometry.
Perhaps we need to introduce a new keyword (see above).

[Alberto Garcia (albertog)]

The issue is complex enough that in-program heuristics might not be able to treat it in all cases. I think it is better to allow the user to request explicitly the behavior desired, and keeping the current (including this patch) heuristics as default behaviour.

So I would introduce a new flag, something like "SCF.Force.Mix.First.Iter" and let it override all the internal heuristics (even for "dangerous" cases), with proper warnings for the user.

But:

-- Do the worries about non-idempotency or memory effects in the mixed DM apply when one is really mixing the hamiltonian?

-- I missed the orginal discussion of the above bug report, but what I can gather now is that the dH found by the user was unavoidable, and other features (reading H, re-reading more mixing history) would be more helpful.

[Nick Papior (nickpapior)]

I agree.
I will amend this with a SCF.Mix.First.Force (then it is equivalent to the SCF.Mix.First with another flag).

I will make SCF.Mix.First.Force default to true since this is typically what people *may* expect when going large scale.

-- non-idempotency: probably this has less effects when doing H-mixing. But ideally we should have a simple test that we can use to uncover some of these details. Effectively this is unexplored (from my side?) We should probably have a small test that ensures the DM retains idempotency.

-- Yes.

[Nick Papior (nickpapior)]

Sorry, will make SCF.Mix.First.Force default false. Just in case ;)

[lp:~nickpapior/siesta/4.1-scf-first]

1085. By Nick Papior on 2019-08-16

Added SCF.Mix.First.Force to overrule siesta heuristics

Now the user can decide to *always* override the Siesta heuristics
regarding mixing the first SCF step. This may be useful for very large
systems.

Added a print-out to users regarding this forced behaviour when
the heuristic would have chosen another method.

[Nick Papior (nickpapior)]

I have amended the requested changes.

I'll merge once you approve. :)

[lp:~nickpapior/siesta/4.1-scf-first]

1086. By Alberto Garcia on 2019-08-19

Longer explanation in manual. Change in warning

[Nick Papior (nickpapior)]

Is this now ready for merge?

[Alberto Garcia (albertog)]

Yes. This is ready to merge. I suggest the following commit message:

Added SCF.Mix.First.Force to overrule siesta heuristics

Please see the manual for the rationale and an explanation of the interactions
between this flag and the existing SCF.Mix.First.

review: Approve

=== modified file 'Docs/siesta.tex'
--- Docs/siesta.tex	2019-08-15 23:47:04 +0000
+++ Docs/siesta.tex	2019-08-19 08:18:39 +0000
@@ -4102,13 +4102,48 @@
 
 \begin{fdflogicalT}{SCF.Mix!First}
   \fdfindex*{DM.MixSCF1}[|see SCF.Mix.First]
+  \fdfdepend{SCF.Mix!First.Force}%
+  
+  This flag is used to decide whether mixing (of the DM or H) should
+  be done in the first SCF step. If mixing is not performed the output
+  DM or H generated in the first SCF step is used as input in the next
+  SCF step. When mixing the DM, this ``reset'' has the effect of
+  avoiding potentially undesirable memory effects: for example, a DM
+  read from file which corresponds to a different structure might not
+  satisfy the correct symmetry, and mixing will not fix it. On the
+  other hand, when reusing a DM for a restart of an interrupted
+  calculation, a full reset might not be advised.
+  
+  The value of this flag is one of the ingredients used by \siesta\ to
+  decide what to do. If \fdftrue\ (the default), mixing will be
+  performed in all cases, except when a DM has been read from file and
+  the sparsity pattern of the DM on file is different from the current
+  one. To ensure that a first-step mixing is done even in this case,
+  \fdf{SCF.Mix!First.Force} should be set to \fdftrue.
+  
+  If the flag is \fdffalse, no mixing in the first step will be
+  performed, except if overridden by \fdf{SCF.Mix!First.Force}.
+  
+  \note that the default value for this flag has changed from the old
+  (pre-version 4) setting in \siesta. The new setting is most
+  appropriate for the case of restarting calculations. On the other
+  hand, it means that mixing in the first SCF step will also be
+  performed for the standard case in which the initial DM is built as
+  a (diagonal) superposition of atomic orbital occupation values. In
+  some cases (e.g. spin-orbit calculations) better results might be
+  obtained by avoiding this mixing.
 
-  Whether the first SCF should be mixed or it uses the output as input
-  in the next SCF step. It is generally advised to set this to true,
-  at least when restarting calculations.
-  
 \end{fdflogicalT}
 
+\begin{fdflogicalF}{SCF.Mix!First.Force}
+
+  Force the mixing (of DM or H) in the first SCF step, regardless of
+  what \siesta\ may heuristically decide.
+
+  This overrules \fdf{SCF.Mix!First}.
+  
+\end{fdflogicalF}
+
 
 In the following the density matrix ($\DM$) will be used in the
 equations, while for Hamiltonian mixing, $\DM$, should be replaced by
 
=== modified file 'Src/m_new_dm.F90'
--- Src/m_new_dm.F90	2018-09-14 12:29:11 +0000
+++ Src/m_new_dm.F90	2019-08-19 08:18:39 +0000
@@ -79,6 +79,16 @@
     logical :: DM_init ! Initialize density matrix from file or from atomic density
     logical :: read_DM
     integer :: DM_in_history, n_depth
+
+    ! Variable to signal how the initialization was performed.
+    ! The following table lists the (current) possibilities:
+    !  0: DM is filled from atomic information (no information present)
+    !  1: .DM file is read, Siesta continuation
+    !  2: An extrapolation of previous geometries.
+    !     The DM's from the previous coordinates are kept in memory
+    !  3: .TSDE file is read, TranSiesta continuation
+    ! If the value is negative it means that one cannot mix in the
+    ! first step, unless SCF.Mix.First.Force is true!
     integer :: init_method
 
     real(dp), pointer :: DM(:,:), EDM(:,:)
@@ -192,8 +202,8 @@
           call print_type(DM_2D)
        end if
 
-       ! Defines the initialiazation to be re-using
-       init_method = -1
+       ! A negative value specifies that one cannot mix in the initial step
+       init_method = -3
 
     end if
 
@@ -290,7 +300,7 @@
 
           old_Ef = Ef
           Ef = fdf_get('TS.Fermi.Initial',Ef,'Ry')
-          if ( init_method == 2 ) then
+          if ( abs(init_method) == 3 ) then
              ! As the fermi-level has been read in from a previous
              ! calculation (TSDE), the EDM should only be shifted by the difference
              diff_Ef = Ef - old_Ef
@@ -301,10 +311,10 @@
 
           if ( IONode ) then
              write(*,*) ! new-line
-             if ( init_method < 2 ) then
+             if ( abs(init_method) < 3 ) then
                 write(*,'(a,f9.5,a)')'transiesta: Setting the Fermi-level to: ', &
                      Ef / eV,' eV'
-             else if ( init_method == 2 ) then
+             else if ( abs(init_method) == 3 ) then
                 write(*,'(a,2(f10.6,a))')'transiesta: Changing Fermi-level from -> to: ', &
                      old_Ef / eV,' -> ',Ef / eV, ' eV'
              end if
@@ -325,7 +335,7 @@
     ! Determine how the mixing of the first step should be performed
     !
     ! The idea is that one will allow mixing on the first SCF step
-    ! only for atomicly filled orbitals. This will in most cases
+    ! only for atomically filled orbitals. This will in most cases
     ! be a good initial choice, while in certain systems (spin-orbit)
     ! it may not be so good.
     ! Subsequently for any MD steps beyond the initial step, we will not
@@ -336,8 +346,19 @@
     ! iterations. For instance when performing FC runs (externally driven)
     ! this would lead to non-degenerate transverse eigenmodes for simple molecules
     if ( init_method == 0 ) then ! atomicly filled data
-       ! allow mix_scf_first
-    else if ( mix_scf_first ) then
+      ! allow mix_scf_first
+    else if ( mix_scf_first_force ) then
+      ! user requested a mix in the first step
+      mix_scf_first = .true.
+      if ( IONode .and. init_method < 0 ) then
+        ! Warn the user about this
+        write(*,"(a)") "Mixing first scf iteration is *forced* although the sparsity pattern has changed..."
+        write(*,"(a)") "Mixing non-compatible matrices might result in problems."
+        write(*,"(a)") "Please use a linear mixer for the first few steps to remove history effects."
+      end if
+    else if ( mix_scf_first .and. init_method < 0 ) then
+       ! Do not allow mixing first SCF step since we are reusing a DM
+       ! from another geometry.
        mix_scf_first = .false.
     end if
     
@@ -385,7 +406,7 @@
     ! integer init_method           : returns method it has read the data with
     !                                   0 == atomic filling (possibly user-defined
     !                                   1 == .DM read
-    !                                   2 == .TSDE read
+    !                                   3 == .TSDE read
     ! *******************************************************************
 
     ! The spin-configuration that is used to determine the spin-order.
@@ -457,7 +478,7 @@
     else 
 
        ! Print-out whether transiesta is starting, or siesta is starting
-       call ts_method_init( init_method == 2 )
+       call ts_method_init( abs(init_method) == 3 )
 
     end if
 
@@ -559,7 +580,7 @@
     !                                 method used to read the DM/EDM
     !                                   0 == not read
     !                                   1 == .DM read
-    !                                   2 == .TSDE read
+    !                                   3 == .TSDE read
     ! *******************************************************************
 
     !> The spin-configuration that is used to determine the spin-order.
@@ -592,6 +613,7 @@
 
     ! The currently read stuff
     type(dSpData2D) :: DM_read
+    type(Sparsity), pointer :: sp_read
     type(dSpData2D) :: EDM_read
 
     ! Signal the file has not been found.
@@ -614,7 +636,7 @@
 
        if ( TSDE_found ) then
           ! Signal we have read TSDE
-          init_method = 2
+          init_method = 3
 
           DM_found = .true.
 
@@ -670,6 +692,16 @@
     ! Density matrix size checks
     if ( DM_found ) then
 
+      ! Specify the *correct* init_method
+      ! In cases where the sparsity pattern changes, we assume this has to do with
+      ! MD simulations and thus a geometry change.
+      ! By default we do not allow mixing the first step when re-using a DM from another
+      ! geometry. This is because of idempotency.
+      sp_read => spar(DM_read)
+      if ( .not. equivalent(sp, sp_read) ) then
+        init_method = -init_method
+      end if
+
       corrected_nsc = .false.
       if ( nsc_read(1) /= 0 .and. any(nsc /= nsc_read) ) then
 
 
=== modified file 'Src/read_options.F90'
--- Src/read_options.F90	2019-07-29 10:25:09 +0000
+++ Src/read_options.F90	2019-08-19 08:18:39 +0000
@@ -413,8 +413,14 @@
   mix_scf_first = fdf_get('DM.MixSCF1', &
        .not. compat_pre_v4_DM_H)
   mix_scf_first = fdf_get('SCF.Mix.First', mix_scf_first)
+  mix_scf_first_force = fdf_get('SCF.Mix.First.Force', .false.)
+  if ( mix_scf_first_force ) then
+    ! Also set this, to note the user of mixing first SCF regardless
+    ! of flag.
+    mix_scf_first = .true.
+  end if
   if (ionode) then
-     write(6,1) 'redata: Mix DM in first SCF step',mix_scf_first
+    write(6,1) 'redata: Mix DM in first SCF step',mix_scf_first
   endif
 
   if (cml_p) then
 
=== modified file 'Src/siesta_options.F90'
--- Src/siesta_options.F90	2018-08-24 13:23:21 +0000
+++ Src/siesta_options.F90	2019-08-19 08:18:39 +0000
@@ -29,6 +29,7 @@
 
   
   logical :: mix_scf_first ! Mix first SCF step?
+  logical :: mix_scf_first_force ! Mix first SCF step? and force it!
   logical :: mix_charge    ! New: mix fourier components of rho
   logical :: mixH          ! Mix H instead of DM
   logical :: h_setup_only  ! H Setup only
 
=== modified file 'version.info'
--- version.info	2019-08-19 07:06:27 +0000
+++ version.info	2019-08-19 08:18:39 +0000
@@ -1,1 +1,5 @@
+<<<<<<< TREE
 siesta-4.1-1093
+=======
+siesta-4.1-1083--scf-first-3
+>>>>>>> MERGE-SOURCE