Run all tests, solved merge problems and updated convergence criteria

- Re-run all tests.

  This showed a couple of interesting things.

  1. There was a mistake in the compute_dm code after
  the PEXSI merge when spin-orbit coupling was introduced
  2. I have removed Hprev as it essentially is the same as
  Hold. Either way it shouldn't produce a huge difference
  in the tracking of the dEbs, etc. (after all they are
  not physically used other than for convergence criteria)
  3. The change to SCF.Mix Hamiltonian resulted in a huge
  number of changes in the output. This is because the first
  step prints out the energies at INIT. However, the Hamiltonian
  is different because it is initialized after the compute_dm step.

- Changed the logic in convergence criteria.
  Now the convergence criterias are additive and may be fully controlled.
  However, at least one convergence criteria must be used.

  Now the default convergence criteria is both the Hamiltonian and the density
  matrix.

  This is update...

Docs/compatibility.tex

@@ -39,6 +39,11 @@ This list describes the compatibility issues when using different versions of \s
     Nearly all flags regarding the mixing options have changed name. However, the old
     flags are still read (and used). In future versions they may be fully removed.
 
+    \begin{itemize}
+      \item All mixing routines are re-written. This means that one cannot compare with
+      prior versions in terms of the convergence path. Sometimes they will differ due to numerics.
+    \end{itemize}
+
     \item \fdf{NumberOfSpecies} and \fdf{NumberOfAtoms} are now defaulted to \# of lines
     in their respective blocks. If they are specified, however, they will be honoured.
 

Docs/siesta.tex

@@ -3518,7 +3518,7 @@ The current version introduces several changes:
   ``right'' on the basis of DM(out), exhibit somehow better convergence
   as a function of the scf step. In order to gain some more data and
   heuristics on this we have implemented a force-monitoring option,
-  activated by setting to \texttt{.true.} the variable \textbf{MonitorForcesInSCF}. The program will then print the maximum
+  activated by setting to \texttt{.true.} the variable \fdf{SCF.MonitorForces}. The program will then print the maximum
   absolute value of the change in forces from one step to the
   next. Other statistics could be implemented.
 
@@ -4906,104 +4906,127 @@ option only works with \tsiesta.
 \subsubsection{Convergence criteria}
 \index{SCF convergence criteria}
 
-\textbf{NOTE}: Some of the options below have a \texttt{DM} prefix. This is
-for historical reasons, as the density-matrix convergence was
-typically the only criterion for termination of the scf cycle.
+\textbf{NOTE}: The older options with a \texttt{DM} prefix is still
+working for backwards compatibility. However, these flags has
+precedence.
+
+Note that all convergence criteria are additive and may thus be used
+simultaneously for complete control.
+
+\begin{fdflogicalT}{SCF.Converge!DM}
+  \index{SCF!mixing!Density matrix convergence}
+
+  Logical variable to use the density matrix elements as monitor
+  of self-consistency.
+  
+\end{fdflogicalT}
+
+\begin{fdfentry}{SCF.Tolerance!DM}[real]<$10^{-4}$>
+  \fdfindex*{DM.Tolerance}
+
+  Tolerance of Density Matrix.
+  %
+  When the maximum difference between the output and the input on each
+  element of the DM in a SCF cycle is smaller than
+  \fdf{SCF.Tolerance!DM}, the self-consistency has been achieved.
 
-\begin{description}
-\item[\textbf{DM.Tolerance}] (\textit{real}):
-\index{DM.Tolerance@\textbf{DM.Tolerance}}
-Tolerance of Density Matrix.
-When the maximum difference between the output and the
-input on each element of the DM
-in a SCF cycle is smaller than DM.Tolerance,
-the selfconsistency has been achieved.
 
-\textit{Default value:} {$10^{-4}$}
+  \note \fdf{DM.Tolerance} is the actual default for this flag.
 
+\end{fdfentry}
+
+\begin{fdfentry}{DM.Normalization.Tolerance}[real]<$10^{-5}$>
+
+  Tolerance for unnormalized density matrices (typically the product
+  of solvers such as PEXSI which have a built-in electron-count
+  tolerance). If this tolerance is exceeded, the program stops. It is
+  understood as a fractional tolerance. For example, the default will
+  allow an excess or shorfall of 0.01 electrons in a 1000-electron
+  system.
+
+\end{fdfentry}
+
+
+
+\begin{fdflogicalT}{SCF.Converge!H}
+  \index{SCF!mixing!Hamiltonian convergence}
 
+  Logical variable to use the Hamiltonian matrix elements as monitor
+  of self-consistency: this is considered achieved when the maximum
+  absolute change (dHmax) in the H matrix elements is below
+  \fdf{SCF.Tolerance!H}. The actual meaning of dHmax depends on
+  whether DM or H mixing is in effect: if mixing the DM, dHmax refers
+  to the change in H(in) with respect to the previous step; if mixing
+  H, dHmax refers to H(out)-H(in) in the previous(?) step. 
+  
+\end{fdflogicalT}
 
-\item[\fdf{DM.Normalization.Tolerance}] (\textit{real}):
-\index{DM.Normalization.Tolerance@\fdf*{DM.Normalization.Tolerance}}
-Tolerance for unnormalized density matrices (typically the product of
-solvers such as PEXSI which have a built-in electron-count tolerance).
-If this tolerance is exceeded, the program stops. It is understood as a
-fractional tolerance. For example, the default will allow an excess or
-shorfall of 0.01 electrons in a 1000-electron system.
+\begin{fdfentry}{SCF.Tolerance!H}[energy]<$10^{-3}\,\mathrm{eV}$>
 
-\textit{Default value:} $10^{-5}$
+  If \fdf{SCF.Converge!H} is \fdftrue, then self-consistency is
+  achieved when the maximum absolute change in the Hamiltonian matrix
+  elements is below this value.
+  
+\end{fdfentry}
 
 
 
-\item[\fdf{DM.Require.Energy.Convergence}] (\textit{logical}):
+\begin{fdflogicalF}{SCF.Converge!FreeE}
   \index{SCF!mixing!energy convergence}
-  \index{DM.Require.Energy.Convergence@\textbf{DM.Require.Energy.Convergence}} Logical variable to request an
-  additional requirement for self-consistency: it is considered
-  achieved when the change in the total (free) energy between cycles
-  of the SCF procedure is below \textbf{DM.EnergyTolerance} and the
-  density matrix change criterion is also satisfied.
+  \fdfindex*{DM.RequireEnergyConvergence}
 
-\textit{Default value:} \texttt{.false.}
+  Logical variable to request an additional requirement for
+  self-consistency: it is considered achieved when the change in the
+  total (free) energy between cycles of the SCF procedure is below
+  \fdf{SCF.Tolerance!FreeE} and the density matrix change criterion is
+  also satisfied.
+
+\end{fdflogicalF}
 
-\item[\textbf{DM.Energy.Tolerance}] (\textit{real energy}):
-  \index{DM.EnergyTolerance@\textbf{DM.EnergyTolerance}} If \fdf{DM.Require.Energy.Convergence} is \texttt{.true.}, then
-  self-consistency is achieved when the change in the total (free)
-  energy between cycles of the SCF procedure is below this value and
-  the density matrix change criterion is also satisfied.
-
-\textit{Default value:} {$10^{-4}$ eV}
-
-\item[\textbf{DM.Require.Harris.Convergence}] (\textit{logical}):
-\index{SCF!mixing!harris energy convergence}
-\index{DM.Require.Harris.Convergence@\textbf{DM.Require.Harris.Convergence}}
-Logical variable to use the Harris energy as monitor of
-self-consistency: this is considered achieved when the change in the Harris energy between cycles
-of the SCF procedure is below \textbf{DM.Harris.Tolerance}. No density
-matrix change criterion is used.
-This is useful if only energies are needed, as the Harris energy tends
-to converge faster than the Kohn-Sham energy.
-The user is responsible for using the correct energies in further
-processing, e.g., the Harris energy if the Harris criterion is used.
-
-To help in basis-optimization tasks, a new file BASIS\_HARRIS\_ENTHALPY
-is provided, holding the same information as BASIS\_ENTHALPY but using
-the Harris energy instead of the Kohn-Sham energy.
+\begin{fdfentry}{SCF.Tolerance!FreeE}[energy]<$10^{-4}\,\mathrm{eV}$>
+  \fdfindex*{DM.EnergyTolerance}
 
-\textit{Default value:} \texttt{.false.}
+  If \fdf{SCF.Converge!FreeE} is \fdftrue, then self-consistency is
+  achieved when the change in the total (free) energy between cycles
+  of the SCF procedure is below this value and the density matrix
+  change criterion is also satisfied.
 
+\end{fdfentry}
 
-\item[\textbf{DM.Harris.Tolerance}] (\textit{real energy}):
-\index{DM.Harris.Tolerance@\textbf{DM.Harris.Tolerance}}
-If \textbf{DM.Require.Harris.Convergence} is \texttt{.true.}, then
-self-consistency is achieved when the change in the Harris energy between cycles
-of the SCF procedure is below this value.
-This is useful if only energies are needed, as the Harris energy tends
-to converge faster than the Kohn-Sham energy.
+\begin{fdflogicalF}{SCF.Converge.Harris}
+  \index{SCF!mixing!harris energy convergence}
+  \fdfindex*{DM.Require.Harris.Convergence}
 
-\textit{Default value:} {$10^{-4}$ eV}
+  Logical variable to use the Harris energy as monitor of
+  self-consistency: this is considered achieved when the change in the
+  Harris energy between cycles of the SCF procedure is below
+  \fdf{SCF.Tolerance!Harris}. This is useful if only
+  energies are needed, as the Harris energy tends to converge faster
+  than the Kohn-Sham energy. The user is responsible for using the
+  correct energies in further processing, e.g., the Harris energy if
+  the Harris criterion is used.
 
-\item[\textbf{SCF.Require.Hamiltonian.Convergence}] (\textit{logical}):
-  \index{SCF!mixing!Hamiltonian convergence}
-  \index{DM.Require.Hamiltonian.Convergence@\textbf{DM.Require.Hamiltonian.Convergence}} Logical variable to use the
-  Hamiltonian matrix elements as monitor of self-consistency: this is
-  considered achieved when the maximum absolute change (dHmax) in the
-  H matrix elements is below \textbf{SCF.H.Tolerance}.  The actual
-  meaning of dHmax depends on whether DM or H mixing is in effect: if
-  mixing the DM, dHmax refers to the change in H(in) with respect to
-  the previous step; if mixing H, dHmax refers to H(out)-H(in) in the
-  previous(?) step. No density matrix change criterion is used.
+  To help in basis-optimization tasks, a new file
+  \file{BASIS\_HARRIS\_ENTHALPY} is provided, holding the same
+  information as \file{BASIS\_ENTHALPY} but using the Harris energy
+  instead of the Kohn-Sham energy.
 
-\textit{Default value:} \texttt{.false.}
+  \note setting this to \fdftrue\ makes \fdf{SCF.Converge!DM}
+  \fdf{SCF.Converge!H} default to \fdffalse.
 
+\end{fdflogicalF}
 
-\item[\textbf{SCF.H.Tolerance}] (\textit{real energy}):
-  \index{SCF.H.Tolerance@\textbf{SCF.H.Tolerance}} If \fdf{SCF.Require.Hamiltonian.Convergence} is \texttt{.true.}, then
-  self-consistency is achieved when the maximum absolute change in the
-  H matrix elements is below this value.
+\begin{fdfentry}{SCF.Tolerance!Harris}[energy]<$10^{-4}\,\mathrm{eV}$>
+
+  If \fdf{SCF.Converge!Harris} is \fdftrue, then self-consistency is
+  achieved when the change in the Harris energy between cycles of the
+  SCF procedure is below this value. This is useful if only energies
+  are needed, as the Harris energy tends to converge faster than the
+  Kohn-Sham energy.
+  
+\end{fdfentry}
 
-\textit{Default value:} \texttt{0.1 mRy}
 
-\end{description}
 
 \vspace{5pt}
 \subsection{The real-space grid and the eggbox-effect}
@@ -8574,6 +8597,8 @@ their values at the previous step for the analysis of the electronic
 structure (band structure calculations, DOS, LDOS, etc).
 \item Some output files reflect the values of the ``un-moved''
   coordinates.
+\item The default convergence criteria is now \emph{both} density
+and Hamiltonian convergence, see \fdf{SCF.Converge!DM} and \fdf{SCF.Converge!H}.
 \end{itemize}
 
 To recover the previous behavior, the user can turn on the

Src/Makefile

@@ -658,8 +658,7 @@ compute_energies.o: atomlist.o class_SpData1D.o class_SpData2D.o dhscf.o
 compute_energies.o: files.o m_dipol.o m_energies.o m_mpi_utils.o m_ntm.o
 compute_energies.o: m_rhog.o m_spin.o precision.o siesta_geom.o
 compute_energies.o: siesta_options.o sparse_matrices.o
-compute_max_diff.o: atomlist.o m_mpi_utils.o m_spin.o precision.o
-compute_max_diff.o: sparse_matrices.o
+compute_max_diff.o: m_mpi_utils.o precision.o
 compute_norm.o: m_mpi_utils.o m_spin.o precision.o sparse_matrices.o
 compute_pw_matrix.o: alloc.o m_planewavematrix.o m_planewavematrixvar.o
 compute_pw_matrix.o: parallel.o precision.o siesta2wannier90.o
@@ -736,14 +735,15 @@ fermid.o: errorf.o parallel.o precision.o sys.o
 fft.o: alloc.o fft1d.o m_timer.o mesh.o parallel.o parallelsubs.o precision.o
 fft.o: sys.o
 fft1d.o: parallel.o precision.o sys.o
-final_H_f_stress.o: alloc.o atomlist.o class_SpData2D.o dnaefs.o files.o
-final_H_f_stress.o: grdsam.o kinefsm.o ldau.o ldau_specs.o m_dipol.o
-final_H_f_stress.o: m_energies.o m_forces.o m_gamma.o m_hsx.o m_mpi_utils.o
-final_H_f_stress.o: m_ncdf_siesta.o m_ntm.o m_spin.o m_steps.o m_stress.o
-final_H_f_stress.o: m_ts_global_vars.o m_ts_io.o m_ts_kpoints.o m_ts_options.o
-final_H_f_stress.o: metaforce.o molecularmechanics.o naefs.o nlefsm.o overfsm.o
-final_H_f_stress.o: parallel.o siesta_geom.o siesta_options.o sparse_matrices.o
-final_H_f_stress.o: spinorbit.o sys.o
+final_H_f_stress.o: alloc.o atomlist.o class_SpData2D.o compute_max_diff.o
+final_H_f_stress.o: dnaefs.o files.o grdsam.o kinefsm.o ldau.o ldau_specs.o
+final_H_f_stress.o: m_dipol.o m_energies.o m_forces.o m_gamma.o m_hsx.o
+final_H_f_stress.o: m_mpi_utils.o m_ncdf_siesta.o m_ntm.o m_spin.o m_steps.o
+final_H_f_stress.o: m_stress.o m_ts_global_vars.o m_ts_io.o m_ts_kpoints.o
+final_H_f_stress.o: m_ts_options.o metaforce.o molecularmechanics.o naefs.o
+final_H_f_stress.o: nlefsm.o overfsm.o parallel.o siesta_geom.o
+final_H_f_stress.o: siesta_options.o sparse_matrices.o spinorbit.o sys.o
+final_H_f_stress.o: units.o
 find_kgrid.o: alloc.o minvec.o parallel.o precision.o units.o
 fire_optim.o: alloc.o m_fire.o m_mpi_utils.o parallel.o precision.o
 fire_optim.o: siesta_options.o units.o

Src/compute_dm.F

@@ -27,8 +27,7 @@
       use m_energies,   only: Ebs, Ecorrec, Ef, Entropy
       use m_rmaxh
       use m_eo
-      use m_spin,   only: NonCol, SpOrb, qs, efs
-      use m_spin,   only: nspin, h_spin_dim, spinor_dim
+      use m_spin,   only: spin, qs, efs
       use m_diagon, only: diagon
       use m_gamma
       use parallel, only: IONode
@@ -53,7 +52,9 @@
       use m_ts_global_vars, only : TSmode, TSinit, TSrun
       use m_transiesta,     only : transiesta
 #endif /* TRANSIESTA */
+      
       implicit none
+      
 !     Input variables
       integer, intent(in) :: iscf
 
@@ -93,7 +94,7 @@
 
       if (SIESTA_worker) then
          if (iscf > 1) then
-            call compute_Ebs_shift(Dscf,H,Hprev,delta_Ebs)
+            call compute_Ebs_shift(Dscf,H,Hold,delta_Ebs)
             delta_Ef = delta_Ebs / qtot
             if (ionode.and.isolve.eq.SOLVE_PEXSI) then
                write(6,"(a,f16.5)")
@@ -113,11 +114,11 @@
         ! This test done in node 0 since NonCol and SpOrb
         ! are not set for PEXSI-solver-only processes
          if (ionode) then
-            if (NonCol .or. SpOrb) call die(
+            if (spin%NCol .or. spin%SO) call die(
      $               "The PEXSI solver does not implement "//
      $               "non-coll spins or Spin-orbit yet")
          endif
-         call pexsi_solver(iscf, no_u, no_l, spinor_dim,
+         call pexsi_solver(iscf, no_u, no_l, spin%spinor,
      $              maxnh, numh, listhptr, listh,
      $              H, S, qtot, Dscf, Escf,
      $              ef, Entropy, temp, delta_Ef)
@@ -149,7 +150,7 @@
 !      &                                '.matrix'
 
 !        Note: only one-shot for now
-         call write_hs_formatted(no_u, nspin,
+         call write_hs_formatted(no_u, spin%H,
      $        maxnh, numh, listhptr, listh, H, S)
          call bye("End of run after writing H.matrix and S.matrix")
 
@@ -162,7 +163,7 @@ c$$$     &           no_u, no_s, maxnh, numh, listhptr, listh,
 c$$$     &           S, 'S.matrix')
 
       elseif ((isolve .eq. SOLVE_DIAGON) .or. (CallDiagon)) then
-        call diagon(no_s, spinor_dim, 
+        call diagon(no_s, spin%spinor, 
      &              no_l, maxnh, maxnh, no_u,
      &              numh, listhptr, listh, numh, listhptr, listh, 
      &              H, S, qtot, fixspin, qs, temp, e1, e2,
@@ -176,26 +177,26 @@ c$$$     &           S, 'S.matrix')
 #endif
       elseif (isolve .eq. SOLVE_ORDERN) then
         if (.not. gamma) call die("Cannot do O(N) with k-points.")
-        if (NonCol .or. SpOrb)
+        if ( spin%NCol .or. spin%SO )
      .      call die("Cannot do O(N) with non-coll spins or Spin-orbit")
         call ordern(usesavelwf, ioptlwf, na_u, no_u, no_l, lasto,
      &               isa, qa, rcoor, rmaxh, ucell, xa, iscf,
      &               istp, ncgmax, etol, eta, qtot, maxnh, numh,
      &               listhptr, listh, H, S, chebef, noeta, rcoorcp,
-     &               beta, pmax, Dscf, Escf, Ecorrec, nspin, qs )
+     &               beta, pmax, Dscf, Escf, Ecorrec, spin%H, qs )
         Entropy = 0.0_dp
       elseif (isolve .eq. SOLVE_MINIM) then
-        if (NonCol .or. SpOrb) 
+        if ( spin%NCol .or. spin%SO ) 
      &      call die('ERROR: Non-collinear spin calculations
      &                       not yet implemented with OMM!')
         H_kin => val(H_kin_1D)
         if (gamma) then
           call dminim(.false., PreviousCallDiagon, iscf, istp, no_l,
-     &                 nspin, no_u, maxnh, numh, listhptr, listh, Dscf,
+     &                 spin%H, no_u, maxnh, numh, listhptr, listh, Dscf,
      &                 eta, qs, H, S, H_kin)
         else
           call zminim(.false., PreviousCallDiagon, iscf, istp, no_l,
-     &                 nspin, no_u, maxnh, numh, listhptr, listh, Dscf,
+     &                 spin%H, no_u, maxnh, numh, listhptr, listh, Dscf,
      &                 eta, qs, no_s, xijo, indxuo, nkpnt, kpoint,
      &                 kweight, H, S, H_kin)
         end if
@@ -204,7 +205,7 @@ c$$$     &           S, 'S.matrix')
         PreviousCallDiagon=.false.
 #ifdef TRANSIESTA
       elseif (TSmode .and. TSinit) then
-        call diagon(no_s, spinor_dim, 
+        call diagon(no_s, spin%spinor, 
      &              no_l, maxnh, maxnh, no_u,
      &              numh, listhptr, listh, numh, listhptr, listh,
      &              H, S, qtot, fixspin, qs, temp, e1, e2,
@@ -216,7 +217,7 @@ c$$$     &           S, 'S.matrix')
 
       else if (TSrun) then
 
-         call transiesta(iscf,nspin, block_dist, sparse_pattern,
+         call transiesta(iscf,spin%H, block_dist, sparse_pattern,
      &        Gamma, ucell, nsc, isc_off, no_u, na_u, lasto, xa, maxnh,
      &        H, S, Dscf, Escf, Ef, Qtot, .false. )
 
@@ -229,9 +230,9 @@ c$$$     &           S, 'S.matrix')
 
 #ifdef CDF
       if ( writedmhs_cdf_history) then
-        call write_dmh_netcdf( no_l, maxnh, h_spin_dim, Dold, H, Dscf )
+        call write_dmh_netcdf( no_l, maxnh, spin%H, Dold, H, Dscf )
       else if (writedmhs_cdf) then
-        call write_dmh_netcdf( no_l, maxnh, h_spin_dim, Dold, H, Dscf,
+        call write_dmh_netcdf( no_l, maxnh, spin%H, Dold, H, Dscf,
      &                         overwrite=.true. )
       endif
 #endif
@@ -240,12 +241,12 @@ c$$$     &           S, 'S.matrix')
       if (muldeb) then
         if (ionode) write (6,"(/a)")
      &    'siesta: Mulliken populations (DM_out)'
-            if (SpOrb) then
+            if ( spin%SO ) then
               call moments( mullipop, na_u, no_u, maxnh, numh, listhptr,
      .                      listh, S, Dscf, isa, lasto, iaorb, iphorb,
      .                      indxuo )
             else
-             call mulliken( mullipop, nspin, na_u, no_u, maxnh,
+             call mulliken( mullipop, spin%DM, na_u, no_u, maxnh,
      &                 numh, listhptr, listh, S, Dscf, isa,
      &                 lasto, iaorb, iphorb )
             endif
@@ -273,8 +274,6 @@ c$$$     &           S, 'S.matrix')
       call normalize_dm(first=.false.)
 #endif
       
-! Save Hamiltonian used to generate the DM in this iteration
-      Hprev = H
 
       call timer( 'compute_dm', 2 )
 #ifdef SIESTA__PEXSI

Src/compute_energies.F90

@@ -63,7 +63,7 @@ CONTAINS
       real(dp) :: buffer1
 #endif
 
-      real(dp) :: const, Escf_out
+      real(dp) :: const
       real(dp) :: dummy_stress(3,3), dummy_fa(1,1)
       real(dp) :: dummy_E, g2max, dummy_H(1,1)
       logical  :: mixDM
@@ -84,9 +84,8 @@ CONTAINS
 
       ! E0 = Ena + Ekin + Enl + Eso - Eions
 
-      DEna = Enascf - Enaatm
-      Etot = E0 + DEna + DUscf + DUext + Exc + &
-           Ecorrec + Emad + Emm + Emeta + Eldau
+      call update_DEna()
+      call update_Etot()
 
 
 ! Harris energy
@@ -120,8 +119,8 @@ CONTAINS
       else if (mix_charge) then
          call compute_correct_EKS()
       endif
-         
-      FreeE  = Etot - Temp * Entropy
+
+      call update_FreeE(Temp)
 
  CONTAINS
 
@@ -242,12 +241,7 @@ CONTAINS
                   Exc, Dxc, dipol, dummy_stress, dummy_fa, dummy_stress)
       ! (when mix_charge is true, rhog will contain the output rho(G))
 
-      DEna     = Enascf - Enaatm
-
-      ! Clarify: 
-      !          DUext (external electric field) -- should it be in or out?
-
-      Escf_out = DEna + DUscf + DUext + Exc 
+      call update_DEna()
 
 !     Compute Tr[H_0*DM_out] = Ekin + Enl + Eso with DM_out
 
@@ -289,9 +283,7 @@ CONTAINS
       ! E0 = Ena + Ekin + Enl + Eso - Eions
 
       ! Clarify: Ecorrec (from O(N))
-      Etot = Ena + Ekin + Enl + Eso - Eions + &
-           Escf_out + Ecorrec + Emad + Emm + &
-           Emeta + Eldau
+      call update_Etot()
       
     end subroutine compute_correct_EKS
 

Src/compute_max_diff.F

-! ---
-! Copyright (C) 1996-2016	The SIESTA group
-!  This file is distributed under the terms of the
-!  GNU General Public License: see COPYING in the top directory
-!  or http://www.gnu.org/copyleft/gpl.txt .
-! See Docs/Contributors.txt for a list of contributors.
-! ---
-      module m_compute_max_diff
-      use precision,         only: dp
-      real(dp) :: dDmax_current
-      public   :: compute_max_diff
-      public   :: dDmax_current
-      CONTAINS
-      subroutine compute_max_diff(Xin,Xout,MaxDiff)
-
-      use sparse_matrices,   only: numh, listhptr
-      use m_spin,            only: nspin
-      use atomlist,          only: no_l
-#ifdef MPI
-      use m_mpi_utils,           only: globalize_max
-#endif
-
-      real(dp), intent(in)    :: Xin(:,:), Xout(:,:)
-      real(dp), intent(out)   :: MaxDiff
-
-      integer   :: is, i, ind, in
-#ifdef MPI
-      real(dp)              :: buffer1
-#endif
-
-      MaxDiff = 0.0_dp
-      do is = 1,nspin
-         do i = 1,no_l
-            do in = 1,numh(i)
-               ind = listhptr(i) + in
-               MaxDiff = max(MaxDiff, abs(Xout(ind,is) - Xin(ind,is)))
-            enddo
-         enddo
-      enddo
-
-#ifdef MPI
-!     Ensure that MaxDiff is the same on all nodes 
-      call globalize_max(MaxDiff,buffer1)
-      MaxDiff = buffer1
-#endif
-
-      dDmax_current = MaxDiff
-
-      end subroutine compute_max_diff
-      end module m_compute_max_diff

Src/compute_max_diff.F90

+! ---
+! Copyright (C) 1996-2016	The SIESTA group
+!  This file is distributed under the terms of the
+!  GNU General Public License: see COPYING in the top directory
+!  or http://www.gnu.org/copyleft/gpl.txt .
+! See Docs/Contributors.txt for a list of contributors.
+! ---
+module m_compute_max_diff
+  use precision,         only: dp
+
+  !> Temporary for storing the old maximum change
+  real(dp), public, save :: dDmax_current
+
+  interface compute_max_diff
+     module procedure compute_max_diff_1d
+     module procedure compute_max_diff_2d
+  end interface compute_max_diff
+  public   :: compute_max_diff
+  
+contains
+  
+  subroutine compute_max_diff_2d(X1,X2, max_diff)
+    
+#ifdef MPI
+    use m_mpi_utils, only: globalize_max
+#endif
+    
+    real(dp), intent(in)  :: X1(:,:), X2(:,:)
+    real(dp), intent(out) :: max_diff
+    
+    integer :: n1, n2
+    integer :: i1, i2
+#ifdef MPI
+    real(dp) :: buffer1
+#endif
+    
+    n1 = size(X1, 1)
+    n2 = size(X1, 2)
+    if ( size(X2, 1) /= n1 ) then
+       call die('compute_max_diff: Sizes of the arrays are not &
+            &conforming (1-D)')
+    end if
+    if ( size(X2, 2) /= n2 ) then
+       call die('compute_max_diff: Sizes of the arrays are not &
+            &conforming (2-D)')
+    end if
+    
+    max_diff = 0.0_dp
+!$OMP parallel do default(shared), private(i2,i1), &
+!$OMP& reduction(max:max_diff), collapse(2)
+    do i2 = 1 , n2
+       do i1 = 1 , n1
+          max_diff = max(max_diff, abs(X1(i1,i2) - X2(i1,i2)) )
+       end do
+    end do
+!$OMP end parallel do
+    
+#ifdef MPI
+    ! Ensure that max_diff is the same on all nodes 
+    call globalize_max(max_diff, buffer1)
+    max_diff = buffer1
+#endif
+
+    dDmax_current = max_diff
+
+  end subroutine compute_max_diff_2d