enh: minor changes to finalize for the release

This adds a link to the latest article describing siesta
developments and allows users to see this.

I also note that users should clearly state which version
of siesta they are using to clarify this for readers.

Added MODULOP which makes mod operations on negative indices
much easier.

Docs/siesta.bib

@@ -189,3 +189,19 @@ url = {https://link.aps.org/doi/10.1103/PhysRevB.100.195417},
 volume = {100},
 year = {2019}
 }
+
+
+@article{doi:10.1063/5.0005077,
+author = {García,Alberto  and Papior,Nick  and Akhtar,Arsalan  and Artacho,Emilio  and Blum,Volker  and Bosoni,Emanuele  and Brandimarte,Pedro  and Brandbyge,Mads  and Cerdá,J. I.  and Corsetti,Fabiano  and Cuadrado,Ramón  and Dikan,Vladimir  and Ferrer,Jaime  and Gale,Julian  and García-Fernández,Pablo  and García-Suárez,V. M.  and García,Sandra  and Huhs,Georg  and Illera,Sergio  and Korytár,Richard  and Koval,Peter  and Lebedeva,Irina  and Lin,Lin  and López-Tarifa,Pablo  and Mayo,Sara G.  and Mohr,Stephan  and Ordejón,Pablo  and Postnikov,Andrei  and Pouillon,Yann  and Pruneda,Miguel  and Robles,Roberto  and Sánchez-Portal,Daniel  and Soler,Jose M.  and Ullah,Rafi  and Yu,Victor Wen-zhe  and Junquera,Javier },
+title = {Siesta: Recent developments and applications},
+journal = {The Journal of Chemical Physics},
+volume = {152},
+number = {20},
+pages = {204108},
+year = {2020},
+doi = {10.1063/5.0005077},
+URL = {https://doi.org/10.1063/5.0005077},
+eprint = {https://doi.org/10.1063/5.0005077}
+}
+
+

Docs/siesta.tex

@@ -342,7 +342,7 @@ time and memory scale linearly with the number of atoms.
 Simulations with several hundred atoms are feasible with
 modest workstations.
 \item
-It is written in Fortran 95 and memory is allocated dynamically.
+It is written in Fortran 2003 and memory is allocated dynamically.
 \item
 It may be compiled for serial or parallel execution (under MPI).
 
@@ -468,6 +468,13 @@ Extensive review of applications as of summer 2003.
 
  Description of the original \tsiesta\ method (prior to 4.1).
 
+ \item
+ ``Siesta: Recent developments and applications'',
+ Alberto Garc\'\ia, \emph{et al.}, J. Chem. Phys. \textbf{152}, 204108
+ (2020).
+
+ Extensive review of applications and developments as of 2020.
+
 \end{itemize}
 
 For more information you can visit the web page

Src/dhscf.F

@@ -301,7 +301,7 @@ C     Vscf and Vaux are not used here
      &    .false., .false. )
         call reord( rhopcc, rhopcc, nml, nsm, TO_SEQUENTIAL )
         if ( debug_dhscf ) then
-           write(*,debug_fmt) Node,'rhopcc',sqrt(sum(rhopcc**2))
+           write(*,debug_fmt) Node,'rhopcc',sum(rhopcc)*dvol
         end if
       endif
 
@@ -381,7 +381,7 @@ C     Create Rhoatm using UNIFORM distr, in sequential form
      &     UNIFORM, rhoatm, TO_SEQUENTIAL )
 
       if ( debug_dhscf ) then
-         write(*,debug_fmt) Node,'rhoatm', sqrt(sum(rhoatm**2))
+         write(*,debug_fmt) Node,'rhoatm', sum(rhoatm)*dvol
       end if
 
 !
@@ -948,7 +948,7 @@ C       Set the UNIFORM distribution again and copy DRho to it
 
         if ( debug_dhscf ) then
            write(*,debug_fmt) Node,'Rho',
-     &          (sqrt(sum(DRho(:,ispin)**2)),ispin=1,nspin)
+     &          (sum(DRho(:,ispin))*dvol,ispin=1,nspin)
         end if
 
         if (save_initial_charge_density) then
@@ -1055,7 +1055,7 @@ C       Find diffuse ionic charge on mesh
         call reord( Chlocal, Chlocal, nml, nsm, TO_SEQUENTIAL )
 
         if ( debug_dhscf ) then
-           write(*,debug_fmt) Node,'Chlocal',sqrt(sum(Chlocal**2))
+           write(*,debug_fmt) Node,'Chlocal',sum(Chlocal)*dvol
         end if
 
 C       Save diffuse ionic charge 
@@ -1141,11 +1141,6 @@ C Find difference between selfconsistent and atomic densities
       enddo
 !$OMP end parallel
 
-      if ( debug_dhscf ) then
-        write(*,debug_fmt) Node,'DRho',
-     &      (sqrt(sum(DRho(:,ispin)**2)),ispin=1,nspin)
-      end if
-
 C ----------------------------------------------------------------------
 C Save electron density difference
 C ----------------------------------------------------------------------
@@ -1201,6 +1196,11 @@ C     Transform spin density into sum and difference
 !$OMP end parallel do
       endif
 
+      if ( debug_dhscf ) then
+        write(*,debug_fmt) Node,'DRho-clean',
+     &      (sum(DRho(:,ispin))*dvol,ispin=1,nspin)
+      end if
+
 C Add a background charge to neutralize the net charge, to
 C model doped systems. It only adds the charge at points
 C where there are atoms (i.e., not in vacuum).          
@@ -1212,6 +1212,11 @@ C First, call with 'task=0' to add background charge
       ! Note that this routine will return immediately if no dopant is present
       call charge_add('+',cell, ntpl, DRho(:,1) )
 
+      if ( debug_dhscf ) then
+        write(*,debug_fmt) Node,'DRho-doped',
+     &      (sum(DRho(:,ispin))*dvol,ispin=1,nspin)
+      end if
+
 C ----------------------------------------------------------------------
 C Calculate the dipole moment
 C ----------------------------------------------------------------------
@@ -1316,6 +1321,7 @@ C ----------------------------------------------------------------------
       Enaatm = Enaatm * dVol
       Enascf = Enaatm + Enascf * dVol
 
+
       if ( SLABDIPOLECORR == SLABDIPOLECORR_VACUUM ) then
         ! In case the user requested a vacuum dipole, calculate it here
         ! REMARK

Src/dipole.F90

@@ -556,7 +556,7 @@ contains
     use mpi_siesta
 #endif
     use units, only: Pi, eV, Ang
-    use intrinsic_missing, only: MODP, vnorm
+    use intrinsic_missing, only: MODULOP, vnorm
 
     real(dp), intent(in) :: cell(3,3)
     integer, intent(in) :: na_u, isa(na_u)
@@ -626,13 +626,14 @@ contains
     if ( dip_vacuum%use_xyz ) then
       ! Find vacuum position
       rc = dot_product(dip_vacuum%xyz, rcell(:,idir))
+      rc = modulo(rc, 1._dp)
       if ( xmin <= rc .and. rc <= xmax ) then
         call die("dipole_potential: Vacuum point lies in atomic range.")
       end if
     else
       rc = (xmin + xmax) / 2 - 0.5_dp
     end if
-    i0 = MODP(nint(rc*ntm(idir)) + 10 * ntm(idir), ntm(idir))
+    i0 = MODULOP(nint(rc*ntm(idir)), ntm(idir))
 
     ! Initialize calculated field
     ! We will finally convert to dipole
@@ -702,7 +703,7 @@ contains
     deallocate(E)
 #endif
 
-    ! Now perform average calcualtion of field along idir
+    ! Now perform average calculation of field along idir
     ! This is a 3-point stencil
     if ( Node == 0 ) then
 

Src/intrinsic_missing.F90

@@ -42,6 +42,8 @@
 !             * NEAREST= 0: if not found, returns -1 (NORMAL behaviour)
 !  - MODP : Returns "FORTRAN" indexed MOD, i.e. if MOD would return 0,
 !           MODP returns the divisor.
+!  - MODULOP : Returns "FORTRAN" indexed MODULO, i.e. if MODULO would return 0,
+!           MODULOP returns the divisor.
 !  - EYE  : a routine to generate identity matrices
 !           this is a subroutine as I think a function would produce
 !           a large memory foot-print before it actually saves to the array.
@@ -95,7 +97,7 @@ module intrinsic_missing
 
 
 ! Elemental functions (can be called on arrays)
-  public :: MODP
+  public :: MODP, MODULOP
 
 ! Missing matrix stuff
   public :: EYE
@@ -234,6 +236,19 @@ contains
     MODP = MOD(a-1,p) + 1
   end function MODP
 
+! A MODULO function which behaves differently on the edge.
+! It will NEVER return 0, but instead return the divisor.
+! This makes it useful in FORTRAN do-loops which are not zero-based
+! but 1-based.
+! Thus we have the following scheme:
+!  MODULOP(x',x) == x         for x' == x
+!  MODULOP(x',x) == MODULOP(x',x) for x' /= x
+  elemental function MODULOP(a,p)
+    integer, intent(in) :: a,p
+    integer :: MODULOP
+    MODULOP = MODULO(a-1,p) + 1
+  end function MODULOP
+
 ! Function to return the unique COUNT of an integer array.
 ! Thus will return how many DIFFERENT entries there exists.
   pure function UNIQC(array)

Src/m_cite.F90

@@ -406,14 +406,17 @@ contains
   end subroutine write_citation
 
   subroutine announce_citations()
+    use version_info, only: version_str
     type(citation) :: cit
     integer :: ID
 
     ! Notify the user about which citations needs to be taken care of
     write(*,'(/,3a)')'cite: Please see "', trim(cite_file), &
         '" for an exhaustive BiBTeX file.'
-    write(*,'(a)') 'cite: This calculation has made use of the following articles.'
-    write(*,'(a)') 'cite: Articles are encouraged to be cited in a published work.'
+    write(*,'(2a)') 'cite: Please clearly indicate Siesta version in published work: ', &
+        trim(version_str)
+    write(*,'(a)') 'cite: This calculation has made use of the following articles'
+    write(*,'(a)') 'cite: which are encouraged to be cited in a published work.'
 
     do ID = 1, N_citations
       if ( used(ID) == 1 ) then ! has been used