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

[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).

[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.

1086. By Alberto Garcia on 2019-08-19

Longer explanation in manual. Change in warning

Unmerged revisions

=== modified file 'Src/m_new_dm.F90'
--- Src/m_new_dm.F90	2018-11-23 10:32:56 +0000
+++ Src/m_new_dm.F90	2019-07-29 12:02:45 +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, regardless of user inquiry!
     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
@@ -337,7 +347,9 @@
     ! 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
+    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 +397,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 +469,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 +571,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 +604,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 +627,7 @@
 
        if ( TSDE_found ) then
           ! Signal we have read TSDE
-          init_method = 2
+          init_method = 3
 
           DM_found = .true.
 
@@ -670,6 +683,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 'version.info'
--- version.info	2019-06-27 10:51:17 +0000
+++ version.info	2019-07-29 12:02:45 +0000
@@ -1,1 +1,5 @@
+<<<<<<< TREE
 trunk-766
+=======
+siesta-4.1-1083--scf-first-1
+>>>>>>> MERGE-SOURCE