Merge lp:~albertog/siesta/trunk-nlso into lp:siesta

[Alberto Garcia]

Implement option to avoid the splitting of the SR and SO contributions in SO 'offsite' calculations.

The recipe for splitting of SR and SO contributions when full lj projectors are used in a SOC calculation is quite fragile. It mostly works when the projectors are computed from semilocal potentials in a .psf file, but even then there might be problems when semicore states (and thus multiple projectors) are used.

An option _not_ to perform the splitting is implemented by this patch. If the line

    soc-split-sr-so F

is included in the fdf file, the program will avoid the split code in nlefsm_SO_off, will
assign to Enl the full lj contribution (done at the time of printing in write_subs), and
will set Eso to zero. New tags 'Enl(+so)' and 'Eso(nil)' are used for printing.

A new 'option' line is printing in the output file, and a new 'Split-SR-SO' parameter added to the CML file.

For generality, the same option is available for SO+onsite calculations.

[lp:~albertog/siesta/trunk-nlso]

701. By Alberto Garcia on 2018-06-04

Sync to trunk-699

[Alberto Garcia (albertog)]

Sorry. I had forgotten to sync to trunk...

I should add that this patch is really needed for the PSML work, as there the projectors might not come from a simple semilocal source.

[Nick Papior (nickpapior)]

The code looks good to me. However, having looked closer at the code there are some points I particularly don't like. E.g. nint on 0.5, the facpm as you mention.

In general these routines could do with a careful re-coding. I would e.g. be cautious on non-IEEE standard optimizations and this code? Does it actually do the right thing?

review: Approve

[lp:~albertog/siesta/trunk-nlso]

702. By Alberto Garcia on 2018-06-07

Safer arithmetic and more comments in nlefsm_SO_off

=== modified file 'Src/Makefile'
--- Src/Makefile	2018-05-30 09:15:30 +0000
+++ Src/Makefile	2018-06-07 08:38:56 +0000
@@ -1136,7 +1136,8 @@
 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
 nlefsm.o: alloc.o atm_types.o atmfuncs.o atomlist.o chemical.o mneighb.o
-nlefsm.o: new_matel.o parallel.o parallelsubs.o precision.o sparse_matrices.o
+nlefsm.o: new_matel.o parallel.o parallelsubs.o precision.o siesta_options.o
+nlefsm.o: sparse_matrices.o
 normalize_dm.o: atomlist.o m_mpi_utils.o m_spin.o parallel.o precision.o
 normalize_dm.o: siesta_options.o sparse_matrices.o sys.o
 obc.o: alloc.o precision.o

=== modified file 'Src/atom.F'
--- Src/atom.F	2018-04-26 13:23:03 +0000
+++ Src/atom.F	2018-06-07 08:38:56 +0000
@@ -3136,7 +3136,7 @@
              if (lj_projs .and. l>0 ) then
               ! Set up the appropriate potential and label
               if (jk == 1) then
-                 ! V(l-j/2) = Vdn - (l+1)/2 Vup = Vion - (l+1)/2 Vso
+                 ! V(l-1/2) = Vdn - (l+1)/2 Vup = Vion - (l+1)/2 Vso
                  vp(:) = vps(:,l) - 0.5_dp*(l+1)*vps_u(:,l)
                  jstr = 'j-'
               else

=== modified file 'Src/nlefsm.f'
--- Src/nlefsm.f	2018-04-24 14:13:24 +0000
+++ Src/nlefsm.f	2018-06-07 08:38:56 +0000
@@ -15,6 +15,7 @@
 
       private
 
+      
       CONTAINS
 
       subroutine nlefsm( scell, nua, na, isa, xa, indxua,
@@ -515,7 +516,10 @@
       use m_new_matel,     only : new_matel
       use atm_types,       only: species_info, species
       use sparse_matrices, only: Dscf, xijo
-     
+      use siesta_options,  only: split_sr_so
+      
+      use fdf
+      
       integer, intent(in) ::
      .   maxnh, na, maxnd, nspin, nua
 
@@ -568,14 +572,15 @@
       logical, dimension(:), pointer ::  listed, listedall
 
       complex(dp) :: E_offsiteSO(4)
-
+      
       type(species_info), pointer        :: spp
 
       integer :: nd, ndn
+      real(dp) :: Vit_saved
 
 C ------------------------------------------------------------
 
-C Start time counte
+C Start time counter
       call timer( 'nlefsm', 1 )
 
 C Find unit cell volume
@@ -659,11 +664,16 @@
 C     Initialize neighb subroutine
       call mneighb( scell, rmax, na, xa, 0, 0, nna )
 
-      nd= 0; ndn= 0
-      Enl = 0.0d0; E_offsiteSO(1:4)=dcmplx(0.0d0,0.0d0)
+      nd= 0
+      ndn= 0
+      Enl = 0.0d0
+      E_offsiteSO(1:4)=dcmplx(0.0d0,0.0d0)
       Enl_offsiteSO = 0.0d0
 C     Loop on atoms with KB projectors      
-      do ka = 1,na      ! Supercell atoms
+
+      ! Use info from hsparse to remove "far" atoms, as in nlefsm
+
+      do ka = 1,na              ! Supercell atoms
        kua = indxua(ka) ! Equivalent atom in the UC
        ks = isa(ka)     ! Specie index of atom ka
        nkb = lastkb(ka) - lastkb(ka-1) ! number of KB projs of atom ka
@@ -673,12 +683,16 @@
 
 C      Find neighbour orbitals
        Ski(:,:) = 0.0_dp
-       nno = 0; ; iano(:)=0; iono(:)=0
+       nno = 0
+       iano(:)=0
+       iono(:)=0
        do ina = 1,nna  ! Neighbour atoms
         ia = iana(ina) ! Atom index of ina (the neighbour to ka)
         is = isa(ia)   ! Specie index of atom ia
         rki = sqrt(r2ki(ina)) ! Square distance
- 
+
+        ! Use this to further filter
+        !!     if (rki - rkbmax(ks) - rorbmax(is) > 0.d0) CYCLE
         do io = lasto(ia-1)+1,lasto(ia) ! Orbitals of atom ia
 
 C        Only calculate if needed locally
@@ -756,7 +770,8 @@
             do ispin = 1,min(2,nspin) ! Only diagonal parts 
              Di(jo) = Di(jo) + Dscf(ind,ispin)
             enddo
-            Ds(1,1,jo) = dcmplx(Dscf(ind,1), Dscf(ind,5))  ! D(ju,iu)
+            ! Should we *add to* Ds, as above for Di?
+            Ds(1,1,jo) = dcmplx(Dscf(ind,1), Dscf(ind,5)) ! D(ju,iu)
             Ds(2,2,jo) = dcmplx(Dscf(ind,2), Dscf(ind,6))  ! D(jd,id)
             Ds(1,2,jo) = dcmplx(Dscf(ind,3), Dscf(ind,4))  ! D(ju,id)
             Ds(2,1,jo) = dcmplx(Dscf(ind,7),-Dscf(ind,8))  ! D(jd,iu)
@@ -789,10 +804,24 @@
 c----------- Compute Vion
              if ( l.eq.0 ) then
               epsk(1) = epskb(ks,koa)
-              Vit = epsk(1) * Ski(koa,ino) * Ski(koa,jno)
+              Vit_saved = epsk(1) * Ski(koa,ino) * Ski(koa,jno)
+              if (split_sr_so) then
+                 Vit = Vit_saved
+              else
+                 ! Move 'ionic' part to V_so diagonal
+                 V_so(1,1,jo)= V_so(1,1,jo) + Vit_saved
+                 V_so(2,2,jo)= V_so(2,2,jo) + Vit_saved
+                 Vit = 0.0_dp
+              endif
               Vi(jo) = Vi(jo) + Vit
               if (.not. matrix_elements_only) then
-               Enl = Enl +  Di(jo) * Vit
+               if (split_sr_so) then
+                  Enl = Enl +  Di(jo) * Vit
+               else
+                  ! Move energy contribution to "Eso"
+                  E_offsiteSO(1) = E_offsiteSO(1) + Vit_saved*Ds(1,1,jo) 
+                  E_offsiteSO(2) = E_offsiteSO(2) + Vit_saved*Ds(2,2,jo) 
+               endif
                CVj  = epsk(1) * Ski(koa,jno)
                Cijk = 2.0_dp * Di(jo) * CVj
                do ix = 1,3
@@ -808,15 +837,18 @@
               ko = ko + 1 
 
 c----------- Compute Vion from j+/-1/2 and V_so
-             else
+             else  ! l /= 0
+              ! First proj of l-1/2 block
               koa1 = -iphKB(ko+1)
+              ! Last proj of l+1/2 block
               koa2 = -iphKB(ko+2*(2*l+1))
               epsk(1) = epskb(ks,koa1)
               epsk(2) = epskb(ks,koa2)
 
               call calc_Vj_offsiteSO( l, epsk, Ski(koa1:koa2,ino), 
      &                       Ski(koa1:koa2,jno), grSki(:,koa1:koa2,ino),
-     &                       grSki(:,koa1:koa2,jno), Vit, V_sot, F_so )
+     &                       grSki(:,koa1:koa2,jno), Vit, V_sot, F_so,
+     &                       split_sr_so)
               Vi(jo) = Vi(jo) + Vit
               V_so(1:2,1:2,jo)= V_so(1:2,1:2,jo) + V_sot(1:2,1:2)
          
@@ -843,7 +875,7 @@
                 enddo
                enddo
               endif
-              ko = ko+2*(2*l+1)
+              ko = ko+2*(2*l+1)    ! Point to next group of l-/+ 1/2 blocks
              endif
              if ( ko.ge.lastkb(ka) ) exit KB_loop
             enddo KB_loop
@@ -906,10 +938,21 @@
 c
 c-----------------------------------------------------------------------
       subroutine calc_Vj_offsiteSO( l, epskb, Ski, Skj, grSki, grSkj,
-     &                       V_ion, V_so, F_so )
+     &                       V_ion, V_so, F_so, split_sr_so )
 
       implicit none
 
+      ! Constructs the NL operator from the l+/- 1/2 information, which
+      ! is passed in two blocks. Hence the re-dimensioning of the dummy
+      ! arguments:
+      !      Ski(-l:l,2) means that there are '2' blocks of 2l+1 values
+      !      and so on.
+      ! The two epskb values correspond to the l-1/2 and l+1/2 blocks, 
+      ! respectively (all the projectors in a block have the same value)
+
+      ! On output, V_so and V_sr (V_ion) are separated, unless split_sr_so is .false.
+      ! The force contributions are not separated
+
       integer     , intent(in)  :: l
       real(dp)    , intent(in)  :: epskb(2)
       real(dp)    , intent(in)  :: Ski(-l:l,2), Skj(-l:l,2)
@@ -917,17 +960,17 @@
       real(dp)    , intent(out) :: V_ion
       complex(dp) , intent(out) :: F_so(3,2,2)
       complex(dp) , intent(out) :: V_so(2,2)
-
-
-      integer    :: J, ij, imj, m, is
+      logical     , intent(in)  :: split_sr_so
+
+      integer    :: J, ij, imj, m, is, facpm
       real(dp)   :: aj, amj, al, a2l1, fac, facm,
-     &              epskpm, V_iont, cp, cm, facpm
+     &              epskpm, V_iont, cp, cm
 
       real(dp)   :: cg(2*(2*l+1),2)
       complex(dp):: u(-l:l,-l:l)
       complex(dp):: SVi(2), SVj(2), grSVi(3,2)
 
-      external   :: die
+      external   :: die, message
 c-----------------------------------------------------------------------
 
 c---- set constants and factors
@@ -935,14 +978,29 @@
       a2l1 = dble( 2*l+1 )
 
 c---- load Clebsch-Gordan coefficients; cg(J,+-)
+      !
+      ! NOTE that this code will not work for l=0, since in this case
+      ! there is a single j subspace (j=0.5)
+      if ( l == 0 ) then
+         call die("Code in calc_Vj_offsiteSO does not work for l=0")
+      endif
+      
       J = 0
       cg(:,:) = 0.0_dp
+
+      ! ij ranges over the two j subspaces
+      ! The loop could have used: 'do ik = -1,1,2' for easier reading
+
       do ij = 1, 2
        aj = al + (2*ij-3)*0.5d0        ! j(ij=1)=l-1/2; j(ij=2)=l+1/2
-       facpm= (-1.0d0)**(aj-al-0.5d0)  ! +/- sign
-       do imj = 1, nint(2*aj)+1        ! Degeneracy for j
-        amj = -aj + dfloat(imj-1)      ! mj value
-        J = J+1                        ! (j,mj) index
+
+       ! This way of computing a sign is very fragile. Better below (and integer)
+       !facpm= (-1.0d0)**(aj-al-0.5d0) ! +/- sign: j=l-1/2: (-1)**(-1)=-1 ;  j = l+1/2: (-1)**0 = 1
+       facpm = 2*ij - 3
+
+       do imj = 1, nint(2*aj)+1 ! Degeneracy for j: 2j+1. Safe nint as aj is always half integer
+        amj = -aj + (imj-1)     ! mj value: -j, -j+1, ..., j-1, j
+        J = J+1                 ! Combined (j,mj) index  (great notation!)
 
         cp = sqrt( (al+0.5d0+amj)/a2l1 )
         cm = sqrt( (al+0.5d0-amj)/a2l1 )
@@ -969,26 +1027,28 @@
       enddo
 
 c---- Load V_so
-      V_so= cmplx(0.0d0,0.0d0); F_so= cmplx(0.0d0,0.0d0)
+      V_so= cmplx(0.0d0,0.0d0)
+      F_so= cmplx(0.0d0,0.0d0)
       J = 0
       do ij = 1, 2
        aj = al + (2*ij-3)*0.5d0        ! j value
        do imj = 1, nint(2*aj)+1        ! Degeneracy for j
-        amj = -aj + dfloat(imj-1)      ! mj value
-        J = J+1                        ! (j,mj) index
+        amj = -aj + imj-1              ! mj value
+        J = J+1                        ! Combined (j,mj) index
 
-        SVi(1:2)= cmplx(0.0d0,0.0d0); SVj(1:2)= cmplx(0.0d0,0.0d0)
+        SVi(1:2)= cmplx(0.0d0,0.0d0)
+        SVj(1:2)= cmplx(0.0d0,0.0d0)
         grSVi(1:3,1:2)= cmplx(0.0d0,0.0d0)
         do is = 1, 2  ! spin loop
 
-c        select correct m
+         ! select correct m. Safe nint as amj is always a half-integer
          if ( is.eq.1 ) then
           m = nint(amj-0.5d0)    ! up   => m=mj-1/2
          else
           m = nint(amj+0.5d0)    ! down => m=mj+1/2
          endif
      
-         if ( iabs(m).le.l ) then
+         if ( abs(m).le.l ) then
           SVi(is)= Ski(+M,ij)*u(+m,M)
           SVj(is)= Skj(+M,ij)*u(+m,M)
           grSVi(1:3,is)= grSki(1:3,+M,ij)*u(+m,M)
@@ -1005,6 +1065,8 @@
          endif
         enddo ! is
 
+        ! Note that these involve just one epskb at a time.
+        
 c       up-up = <i,+|V,J><V,J|j,+>
         V_so(1,1)  = V_so(1,1)  + SVi(1) * epskb(ij) * conjg(SVj(1))
         F_so(:,1,1)= F_so(:,1,1)+ grSVi(:,1) * epskb(ij) * conjg(SVj(1))
@@ -1031,10 +1093,36 @@
        call die('calc_Vj_LS: ERROR')
       endif
 
-c---- substract out V_ion
-      epskpm = sqrt( epskb(1)*epskb(2) )
-      epskpm = sign(epskpm,epskb(1))
- 
+      if (split_sr_so) then
+
+      ! Attempt to get SR and SO parts of the total lj non-local contribution
+      ! (Note that there is an extra contribution to SR coming from l=0, computed
+      ! directly in the caller routine)
+
+      ! This is some kind of average, following Hemstreet, of the l+/- 1/2 
+      ! components. (In contrast to Hamann's more involved procedure)
+
+      !  v_sr = ( (l+1) v_j+ + l v_j- ) / (2l+1)
+      !  and the v_j+ and v_j- are "square roots" of the KB projectors      
+
+      !!  epskpm = sqrt( epskb(1)*epskb(2) ) ! This sqrt should be guarded with an abs()
+      !!  epskpm = sign(epskpm,epskb(1))   ! The value of epskpm with the sign of epskb(1) WHY?
+
+      ! This is the only acceptable solution from symmetry arguments
+      ! (and within the fragility of the approach)
+      epskpm  = sqrt(abs(epskb(1)*epskb(2)))
+      if (epskb(1)*epskb(2) > 0) then
+         ! same sign
+         epskpm  = epskpm * sign(1.0_dp,epskb(1)) 
+      else
+         ! different sign: this possibility was not taken into account in previous versions,
+         ! and probably implies that the Hemstreet ansatz is ill-defined
+         ! We should flag this instead of accepting this heuristic blindly.
+         epskpm  = - epskpm
+         call message("WARNING",
+     $        "The Enl-Eso split in energies might not be accurate")
+      endif
+      
       V_ion = 0.0d0
       do M = -l, l
        V_iont = ( l**2     * Ski(M,1)*epskb(1)*Skj(M,1)
@@ -1044,10 +1132,15 @@
        V_ion = V_ion + V_iont
       enddo
    
+      !---- substract out V_ion from V_so
       V_so(1,1) = V_so(1,1) - cmplx(1.0d0,0.0d0)*V_ion
       V_so(2,2) = V_so(2,2) - cmplx(1.0d0,0.0d0)*V_ion
 
-      return
+      else
+         !---- Keep all the NL contribution in V_so
+         V_ion = 0.0_dp
+      endif
+
       end subroutine calc_Vj_offsiteSO
 
       end module m_nlefsm

=== modified file 'Src/read_options.F90'
--- Src/read_options.F90	2018-05-15 13:01:22 +0000
+++ Src/read_options.F90	2018-06-07 08:38:56 +0000
@@ -27,7 +27,7 @@
   use units,     only : eV, Ang, Kelvin
   use siesta_cml
   use m_target_stress, only: set_target_stress
-  use m_spin, only: print_spin_options
+  use m_spin, only: print_spin_options, spin
 
   use m_charge_add, only : read_charge_add
   use m_hartree_add, only : read_hartree_add
@@ -843,6 +843,17 @@
           units='siestaUnits:eSpin' )
   endif
 
+  ! For SOC calculations: If .false., Enl will contain
+  ! the SO part of the energy.
+  if (spin%SO) then
+     split_sr_so = fdf_get('SOC.Split.SR.SO',.true.)
+     write(6,1) 'redata: Split SR and SO contributions', split_sr_so
+     if (cml_p) then 
+        call cmlAddParameter( xf=mainXML, name='Split-SR-SO', &
+             value=split_sr_so, dictref='siesta:split_sr_so' )
+     endif
+  endif
+  
   ! Order-N solution parameters ...
   !     Maximum number of CG minimization iterations
   ncgmax = fdf_get('ON.MaxNumIter',1000)

=== modified file 'Src/siesta_options.F90'
--- Src/siesta_options.F90	2018-05-15 13:01:22 +0000
+++ Src/siesta_options.F90	2018-06-07 08:38:56 +0000
@@ -126,6 +126,7 @@
   logical :: muldeb        ! Write Mulliken populations at every SCF step?
   logical :: spndeb        ! Write spin-polarization information at every SCF step?
   logical :: orbmoms       ! Write orbital moments?
+  logical :: split_sr_so   ! Cosmetic: split full lj NL energies into SR and SO parts
 
   ! Convergence options
   logical :: converge_FreeE   ! free Energy conv. to finish SCF iteration?

=== modified file 'Src/write_subs.F'
--- Src/write_subs.F	2018-04-19 10:08:47 +0000
+++ Src/write_subs.F	2018-06-07 08:38:56 +0000
@@ -304,7 +304,8 @@
       real(dp), intent(in) :: dHmax ! Max. change in H elements
 
       character(len=64) :: fmt
-      character(len=6) :: scf_name
+      character(len=6)  :: scf_name
+      character(len=17) :: enl_str, eso_str
 
       logical              :: first_scf_step
       integer              :: i
@@ -313,6 +314,16 @@
       if ( TSrun ) then
          scf_name = 'ts-scf'
       end if
+
+      if (spin%SO .and. (.not. split_sr_so)) then
+         Enl = Enl + Eso 
+         Eso = 0.0_dp
+         enl_str = 'siesta: Enl(+so)='
+         eso_str = 'siesta: Eso(nil)='
+      else
+         enl_str = 'siesta: Enl     ='
+         eso_str = 'siesta: Eso     ='
+      endif
       
       first_scf_step = (iscf == 1)
       ! Only print out full decomposition at very beginning and end.
@@ -323,8 +334,8 @@
      .     'siesta: Eions   =', Eions/eV,
      .     'siesta: Ena     =', Ena/eV,
      .     'siesta: Ekin    =', Ekin/eV,
-     .     'siesta: Enl     =', Enl/eV,
-     .     'siesta: Eso     =', Eso/eV,
+     .      enl_str, Enl/eV,
+     .      eso_str, Eso/eV,
      .     'siesta: Eldau   =', Eldau/eV,
      .     'siesta: DEna    =', DEna/eV,
      .     'siesta: DUscf   =', DUscf/eV,
@@ -496,13 +507,20 @@
         
       else !final
 !       Print out additional information in finalization.
+
+        if (spin%SO .and. (.not. split_sr_so)) then
+           eso_str = '     Eso(nil) ='
+        else
+           eso_str = '      Eso     ='
+        endif
+
         write(6,'(/,a)') 'siesta: Final energy (eV):'
         write(6,'(a,a15,f15.6)')
      .    'siesta: ', 'Band Struct. =', Ebs/eV,
      .    'siesta: ',      'Kinetic =', Ekin/eV,
      .    'siesta: ',      'Hartree =', Uscf/eV,
      .    'siesta: ',      'Eldau   =', Eldau/eV,
-     .    'siesta: ',      'Eso     =', Eso/eV,
+     .    'siesta: ',      eso_str    , Eso/eV,
      .    'siesta: ',   'Ext. field =', DUext/eV,
      .    'siesta: ',      'Enegf   =', DE_NEGF/eV,
      .    'siesta: ',  'Exch.-corr. =', Exc/eV,

=== modified file 'version.info'
--- version.info	2018-05-30 09:15:30 +0000
+++ version.info	2018-06-07 08:38:56 +0000
@@ -1,1 +1,1 @@
-trunk-699
+trunk-699--nlso-702