Merge lp:~nickpapior/siesta/4.0-nrtab into lp:siesta/4.0

[Nick Papior]

Enabled old NRTAB behaviour via pre-processing flags.

By specifying SIESTA__LEGACY_NRTAB one recaptures the old
NRTAB behaviour which may be used to compare against older
siesta versions.

The compatibility.tex document has been updated to reflect this.

I simply added what we discussed at Slack (changing the pre-processor name).

[Alberto Garcia (albertog)]

Please forgive me for having second (third?) thoughts about this. An increasing number of users now depend on having SIESTA installed by a third person, be it a supercomputer center's staff or some grumpy sysadmin who cannot be counted on to provide a full range of SIESTA versions with all possible preprocessing options...

Thus, those users will be most likely locked into the new default value for NRTAB, and will be unable to reproduce older results with the new version.

This problem is not restricted to this issue, of course.

I have tentatively produced a patch that will make NRTAB configurable at two levels: a "legacy" switch to reproduce old results, and a full setter for people who want to try different values (not all of them would be allowed, as NQ must be compatible with NRTAB...).

(I will attach the patch somehow)

review: Needs Fixing

[lp:~nickpapior/siesta/4.0-nrtab]

561. By Nick Papior on 2018-04-05

Implemented NRTAB compatibility as an option

This enables users to change NRTAB with fdf files
rather than recompiling siesta.

[Nick Papior (nickpapior)]

I have amended our discussion.

I have also documented the fdf-option.

=== modified file 'Docs/compatibility.tex'
--- Docs/compatibility.tex	2018-02-19 08:15:36 +0000
+++ Docs/compatibility.tex	2018-04-05 07:11:38 +0000
@@ -3,8 +3,8 @@
 \usepackage[a4paper]{geometry}
 
 \newcommand\siesta{\textsc{SIESTA}}
-\newcommand{\opt}[1]{{\bf #1}}
-\newcommand{\code}[1]{{\tt #1}}
+\newcommand{\opt}[1]{\textbf{#1}}
+\newcommand{\code}[1]{\texttt{#1}}
 
 \title{Compatibility issues across versions}
 \author{\siesta}
@@ -23,9 +23,7 @@
 
   \item[\emph{any} --- 4.0.2] %
   Internal tabulated values. We have increased the tabulated grid quantities related to
-  orbital basis functions. This change \emph{increases} precision for Siesta. To compare
-  with prior versions of Siesta set \texttt{integer, parameter :: NRTAB = 128} in
-  \texttt{Src/new\_matel.f}.
+  orbital basis functions. Please see \opt{Compat.Matel.NRTAB} for details.
 
   \item[\emph{any} --- 4.0-b2] The following compatibility issues should be remarked when
   comparing with any later version of \siesta.

=== modified file 'Docs/siesta.tex'
--- Docs/siesta.tex	2018-02-15 15:45:23 +0000
+++ Docs/siesta.tex	2018-04-05 07:11:38 +0000
@@ -4153,6 +4153,20 @@
 
 See also the {\bf WriteDMHS.NetCDF} and {\bf WriteDMHS.History.NetCDF}
 options.
+
+
+\item[{\bf Compat.Matel.NRTAB}] ({\it logical}):
+\index{Compat.Matel.NRTAB@{\bf Compat.Matel.NRTAB@}}
+Internally the matrix element calculations are tabulated with a preset number of
+elements. In versions 4.0.1 and prior this was default to 128. Since 4.0.2 the default
+value has been increased to 1024 which translates to higher precision matrix element
+calculations.
+
+This option is a compatibility option that enables comparison with 4.0.1 and prior
+versions. Should \emph{only} be set to true to compare. Otherwise use the default value.
+
+{\it Default value:} {\tt .false.}
+
 \end{description}
 
 \subsubsection{The auxiliary supercell}

=== modified file 'Src/Makefile'
--- Src/Makefile	2018-02-01 10:42:36 +0000
+++ Src/Makefile	2018-04-05 07:11:38 +0000
@@ -695,7 +695,7 @@
 new_dm.o: precision.o siesta_geom.o siesta_options.o sparse_matrices.o sys.o
 new_dm.o: units.o
 new_matel.o: alloc.o errorf.o interpolation.o matel_registry.o parallel.o
-new_matel.o: precision.o radfft.o spher_harm.o sys.o
+new_matel.o: precision.o radfft.o siesta_options.o spher_harm.o sys.o
 nlefsm.o: alloc.o atmfuncs.o mneighb.o new_matel.o parallel.o parallelsubs.o
 nlefsm.o: precision.o
 normalize_dm.o: atomlist.o m_mpi_utils.o m_spin.o parallel.o precision.o

=== renamed file 'Src/new_matel.f' => 'Src/new_matel.F'
--- Src/new_matel.f	2018-02-19 08:15:36 +0000
+++ Src/new_matel.F	2018-04-05 07:11:38 +0000
@@ -78,6 +78,7 @@
       use spher_harm, only: reset_spher_harm
       use sys,     only: die
       use m_radfft
+      use siesta_options, only: matel_NRTAB
 C -------------------------------------------------------------------
 
 C Argument types and dimensions -------------------------------------
@@ -98,14 +99,14 @@
 C  FFTOL is the tolerance for considering equal the radial part of
 C    two orbitals.
 C  TINY is a small number to avoid a zero denominator
-      INTEGER,          PARAMETER :: NRTAB  =  1024
+      INTEGER, save :: NRTAB = 1024
       INTEGER,          PARAMETER :: NQ     =  1024
       real(dp),         PARAMETER :: EXPAND =  1.20_dp
       real(dp),         PARAMETER :: Q2CUT  =  2.50e3_dp
       real(dp),         PARAMETER :: GWBYDR =  1.5_dp
       real(dp),         PARAMETER :: FFTOL  =  1.e-8_dp
       real(dp),         PARAMETER :: TINY   =  1.e-12_dp
-      CHARACTER(LEN=*), PARAMETER :: MYNAME =  'MATEL '
+      CHARACTER(LEN=*), PARAMETER :: MYNAME =  'MATEL'
 C -------------------------------------------------------------------
 
 C Internal variable types and dimensions ----------------------------
@@ -160,6 +161,12 @@
         CALL RE_ALLOC( INDFFY, 0,MFF, 'INDFFY', MYNAME )
         INDFFY(0) = 0
         NULLIFIED = .TRUE.
+
+!       Set only at specific times. If matel_nrtab is changed the new
+!       value will not take effect unless there is a re-initialization,
+!       to avoid corruption of the tables.
+        NRTAB = matel_NRTAB
+
       ENDIF
 
 C Check if tables must be re-initialized 
@@ -214,7 +221,7 @@
      &     MAX(IOPER,3).GT.UBOUND(INDF,2) ) THEN
         CALL RE_ALLOC( INDF, 1, MAX(UBOUND(INDF,1),IG1,IG2),
      &                       1,MAX(UBOUND(INDF,2),IOPER,3),
-     &                'INDF', MYNAME )
+     &                'INDF', MYNAME)
         CALL RE_ALLOC( NLM,  1,MAX(UBOUND(INDF,1),IG1,IG2),
      &                       1,MAX(UBOUND(INDF,2),IOPER,3),
      &                'NLM', MYNAME )

=== modified file 'Src/read_options.F90'
--- Src/read_options.F90	2017-06-30 08:31:34 +0000
+++ Src/read_options.F90	2018-04-05 07:11:38 +0000
@@ -292,12 +292,26 @@
   partial_charges_at_every_scf_step =  &
        fdf_get('PartialChargesAtEveryScfStep',.false.)
 
+  
+  if ( fdf_get('Compat.Matel.NRTAB', .false.) ) then
+    matel_NRTAB = 128
+  else
+    matel_NRTAB = 1024
+  end if
+  if ( IONode ) then
+    write(6,4) 'redata: Matel table size (NRTAB)         = ', matel_NRTAB
+  end if
+  if (cml_p) then
+    call cmlAddParameter( xf=mainXML, name='MatelNRTAB',value=matel_NRTAB, &
+        dictRef='siesta:matel_nrtab', units="cmlUnits:countable")
+  end if
 
   ! Planewave cutoff of the real space mesh ...
   g2cut = fdf_get('MeshCutoff',g2cut_default,'Ry')
   if (ionode) then
-     write(6,6) 'redata: Mesh Cutoff                      = ', g2cut,'  Ry'
-  endif
+     write(6,6)'redata: Mesh Cutoff                      = ', g2cut,' Ry'
+   endif
+   
 
   if (cml_p) then
      call cmlAddParameter( xf=mainXML, name='MeshCutOff', value=g2cut,     &

=== modified file 'Src/siesta_options.F90'
--- Src/siesta_options.F90	2017-03-05 19:42:40 +0000
+++ Src/siesta_options.F90	2018-04-05 07:11:38 +0000
@@ -9,6 +9,8 @@
   
   implicit none
   integer, parameter, private :: dp = selected_real_kind(10,100)
+
+  save
   
   PUBLIC
 
@@ -181,7 +183,10 @@
   real(dp) :: wmix          ! Mixing weight for DM in SCF iteration
   real(dp) :: wmixkick       ! Mixing weight for DM in special 'kick' SCF steps
 
-  character(len=150) :: sname   ! System name, used to initialise read
+  ! Matrix element compatibility variable
+  integer :: matel_NRTAB = 1024
+  
+  character(len=164) :: sname   ! System name, used to initialise read
 
   integer,  parameter :: SOLVE_DIAGON = 0
   integer,  parameter :: SOLVE_ORDERN = 1

=== modified file 'version.info'
--- version.info	2018-02-22 10:50:23 +0000
+++ version.info	2018-04-05 07:11:38 +0000
@@ -1,1 +1,1 @@
-siesta-4.0--559
+siesta-4.0--559--nrtab-2