Merge lp:~nickpapior/siesta/4.1-naive-aux-cell into lp:siesta/4.1

[Nick Papior]

Reimplemented the NaiveAuxiliaryCell option

And further enhanced it by forcing the aux-cell size to
be the largest one for anything but single-shot calculations.

This does not fix *all* cases of error handling since when performing
variable cell calculations the auxiliary cell may again change.

One solution could be if naive aux-cell is used then the aux-cell is
forced to be the same for cell_can_change calculations?

This reimplements the NaiveAuxiliaryCell option which forces the auxiliary cell to be constant.

However, I would like to get your opinion (Alberto) on the cell_can_change comment above?

[Alberto Garcia (albertog)]

It seems that there are two types of "errors" that can appear:

1. A supercell is not big enough. This causes a crash in 'hsparse', and should be avoided at all costs.

2. The supercell for the new geometry is different (because of cell changes), and thus the DM cannot be reused safely. In this case we could just skip the reuse of the DM for this particular instance.
Hopefully the number of times that this happens is limited, as changes will tend to become smaller and smaller (?).

Now, I believe it is true that the "naive-aux-cell" algorithm gives a safe auxcell for a given unit cell, isn't it? Note that to be completely sure we can run *always* the tight graph algorithm, but I do not think it is necessary.

Then we can leave the code as it is, recomputing nsc at each geometry, but checking whether nscold /= nsc, and in that case forbid the reuse of the DM.

[Nick Papior (nickpapior)]

1) Yes, that error will (should) only happen in MD simulations when the lattice parameters reduces during the optimization (as far as I know).

2) By looking at m_new_dm I think this is already the case. The DM gets re-initialized if the supercell changes. If so the DM history is reset, and the DM gets atomicly initialized.

Yes, the Naive-aux-cell gives the safest auxiliary cell with a minimum of 3x3x3 because it does not take into account actual atomic positions. This is rather unfortunate but will be handled once we have a DM(nsc) -> DM(nsc') routine.

=== modified file 'Docs/siesta.tex'
--- Docs/siesta.tex	2018-06-19 11:00:00 +0000
+++ Docs/siesta.tex	2018-06-25 12:35:49 +0000
@@ -5762,7 +5762,6 @@
  
 \begin{fdflogicalF}{ForceAuxCell}
 
-
 If \fdftrue, the program uses an auxiliary cell even for gamma-point-only
 calculations. This might be needed for COOP/COHP calculations, as
 noted above, \index{COOP/COHP curves!Folding in Gamma-point
@@ -5781,6 +5780,32 @@
 
 \end{fdflogicalF}
 
+\begin{fdfentry}{NaiveAuxiliaryCell}[logical]<\fdfvalue{MD.TypeOfRun}>
+
+  The auxiliary supercell is determined using the atomic positions,
+  the orbital extend and the lattice parameters. For a fixed geometry
+  configuration the auxiliary supercell is uniquely defined. However,
+  when performing MD simulations (or relaxations) the displacement of
+  certain atoms may result in another required auxiliary supercell. In
+  those cases \siesta\ will currently crash unless this option is
+  \fdftrue. When this is \fdftrue\ the used supercell will always be
+  maximized based on the orbital extend (fixed for \emph{any}
+  calculation) and the lattice vectors (depends on
+  \fdf{MD.VariableCell}).
+
+  By default this will be \fdffalse\ when performing
+  \fdf{MD.TypeOfRun:CG} with \fdf{MD.NumCGsteps:1} or below. Otherwise
+  this will default to \fdftrue.
+
+  \note one may try to set this to \fdffalse\ for better performance,
+  however at the risk of \siesta\ crashing after each MD step.
+
+  \note this option may made obsolete in later versions when a
+  supercell conversion routine is in place.
+
+\end{fdfentry}
+
+
 
 \subsection{Calculation of the electronic structure}
 

=== modified file 'Src/read_options.F90'
--- Src/read_options.F90	2018-05-15 09:48:18 +0000
+++ Src/read_options.F90	2018-06-25 12:35:49 +0000
@@ -97,6 +97,8 @@
   !                             6 = Force constants
   !                             7 = Deprecated (Forces for PHONON program)
   !                             8 = Force evaluation
+  !                             9 = Explicit set of coordinates
+  !                            10 = Lua controlled dynamics
   ! integer istart           : Initial time step for MD
   ! integer ifinal           : Final time step for MD
   ! integer nmove            : Number of CG steps in CG optimization
@@ -1579,6 +1581,14 @@
      end select
   end if
 
+
+  ! Default variable depending on the action required
+  if ( idyn == 0 .and. nmove <= 1 ) then
+    naive_aux_cell = fdf_get( 'NaiveAuxiliaryCell', .false. )
+  else
+    naive_aux_cell = fdf_get( 'NaiveAuxiliaryCell', .true. )
+  end if
+  
   writec                 = fdf_get( 'WriteCoorStep', outlng )
   savehs                 = fdf_get( 'SaveHS', .false. )
   initdmaux              = fdf_get( 'ReInitialiseDM', .TRUE. )

=== modified file 'Src/siesta_init.F'
--- Src/siesta_init.F	2018-06-07 20:08:23 +0000
+++ Src/siesta_init.F	2018-06-25 12:35:49 +0000
@@ -609,16 +609,17 @@
 !     Find required supercell
 !     2*rmaxh is used to guarantee that two given orbitals in the
 !     supercell can only overlap once
-      if (gamma) then
+      if ( gamma ) then
         nsc(1:3) = 1
-      else
-        ! Discussion with Jose-Soler we get this size ! NPA
+      else if ( naive_aux_cell ) then
+        ! The auxiliary supercell is symmetric in the connections,
+        ! hence we need a factor of 2.
         do i = 1 , 3
           veclen = sqrt(ucell(1,i)**2+ucell(2,i)**2+ucell(3,i)**2)
           nsc(i) = 1 + 2 * ceiling( rmaxh / veclen )
         end do
-        ! The above is kept for historical reasons,
-        ! but a tight supercell can be found from the atom-graph info:
+      else
+        ! Determine the actual required number of supercells
         call exact_sc_ag(negl,ucell,na_u,isa,xa,nsc)
       end if
 

=== modified file 'Src/siesta_options.F90'
--- Src/siesta_options.F90	2018-05-15 09:48:18 +0000
+++ Src/siesta_options.F90	2018-06-25 12:35:49 +0000
@@ -181,6 +181,9 @@
   integer :: call_diagon_default    ! Default number of SCF steps for which to use diagonalization before OMM
   integer :: call_diagon_first_step ! Number of SCF steps for which to use diagonalization before OMM (first MD step)
 
+  logical :: naive_aux_cell ! Control whether the auxiliary superc-cell will always be the maximal allowed cell size
+                            ! This will default to true for MD, FC and Lua controlled calculations since atoms may move
+  
   real(dp) :: beta          ! Inverse temperature for Chebishev expansion.
   real(dp) :: bulkm         ! Bulk modulus
   real(dp) :: charnet       ! Net electric charge

=== modified file 'Src/state_init.F'
--- Src/state_init.F	2018-05-15 09:48:18 +0000
+++ Src/state_init.F	2018-06-25 12:35:49 +0000
@@ -258,16 +258,16 @@
      &                 'state_init' )
 
 !       Find required supercell
-        if (gamma) then
+        if ( gamma ) then
            nsc(1:3) = 1
+        else if ( naive_aux_cell ) then
+          do i = 1 , 3
+            veclen = sqrt(ucell(1,i)**2+ucell(2,i)**2+ucell(3,i)**2)
+            nsc(i) = 1 + 2 * ceiling( rmaxh / veclen )
+          end do
         else
-           do i = 1 , 3
-              veclen = sqrt(ucell(1,i)**2+ucell(2,i)**2+ucell(3,i)**2)
-              nsc(i) = 1 + 2 * ceiling( rmaxh / veclen )
-           end do
-           ! The above is kept for historical reasons,
-           ! but a tight supercell can be found from the atom-graph info:
-           call exact_sc_ag(negl,ucell,na_u,isa,xa,nsc)
+          ! Determine the actual required number of supercells
+          call exact_sc_ag(negl,ucell,na_u,isa,xa,nsc)
         endif
 
         mscell = 0.0_dp

=== modified file 'version.info'
--- version.info	2018-06-19 11:00:00 +0000
+++ version.info	2018-06-25 12:35:49 +0000
@@ -1,1 +1,1 @@
-siesta-4.1--933
+siesta-4.1--933--naive-aux-cell-1