Merge lp:~nickpapior/siesta/4.1-fix-spin into lp:~albertog/siesta/4.1-fix-spin

[Nick Papior]

Further "fixes" for the non-colinear calculations.
This is just to clarify (in the comments) and in the diag2g routine.

[Alberto Garcia (albertog)]

I need to discuss these things with Ramón in a larger context when we get back to ICMAB.

[Alberto Garcia (albertog)]

I am ready to accept the changes on simplicity arguments (i.e., having a single, compatible way to build the DM in the Noncol and SOC cases).

I will merge this changes in my 4.1-spin-fix branch (which also contains other clarifications), and we will mark this branch for merging into OSSO when it is ready.

review: Approve

=== modified file 'Src/diag2g.F'
--- Src/diag2g.F	2017-12-21 15:49:49 +0000
+++ Src/diag2g.F	2018-01-05 08:40:43 +0000
@@ -32,6 +32,7 @@
 C This version is for non-collinear spin at gamma point.
 C Writen by J.Soler, May and August 1998.
 C Reduced memory requirements, Nick, Aug. 2017
+C Clarified DM calculation, N.R. Papior, Jan. 2018      
 C **************************** INPUT **********************************
 C integer nuo                 : Number of basis orbitals on local node
 C integer no                  : Number of basis orbitals
@@ -107,13 +108,13 @@
       integer           BNode, BTest, ie, ierror, iie, iio
       integer           ind, io, j, jo, nd, iuo, juo
       real(dp)          ee, pipj, qe, t
-      complex(dp)       cicj, D11, D22, D12, D21
+      complex(dp) :: cicj, D11, D22, D12, D21
 #ifdef MPI
       integer MPIerror
 #endif
       external              cdiag
 #ifdef DEBUG
-      call write_debug( '    PRE diag2k' )
+      call write_debug( '    PRE diag2g' )
 #endif
 !***********************************************************************
 !     B E G I N
@@ -136,18 +137,16 @@
 ! 1. Hermiticity imposes H_{i,j}^{is,js}=H_{j,i}^{js,is}^*
 ! 2. Since wave functions are real, if there are no single P or L
 ! operators:
-!      (a) H_{i,j}^{is,js}=H_{j,i}^{is,js}
-!      (b) These imply spin-box hermiticity:
-!      H_{i,j}^{is,js}=H_{i,j}^{js,is}^*
-!
-!      (c) Hence
-!
-!                               | H(ind,1)                H(ind,3) - i
-!                               H(ind,4) |
-!           H_{j,i} = H_{i,j} = |
-!           |
-!                               | H(ind,3) + i H(ind,4)   H(ind,2)
-!                               |
+!     (a) H_{i,j}^{is,js}=H_{j,i}^{is,js}
+!     (b) These imply spin-box hermiticity:
+!
+!          H_{i,j}^{is,js}=H_{i,j}^{js,is}^*
+!
+!     (c) Hence
+!
+!                              | H(ind,1)                H(ind,3) - i H(ind,4) |
+!          H_{j,i} = H_{i,j} = |                                               |
+!                              | H(ind,3) + i H(ind,4)   H(ind,2)              |
 !
 !
 ! The spin-orbit interaction and the orbital part of the Zeeman
@@ -155,7 +154,6 @@
 ! break the  "spin-box hermiticity" of H. Hence the full format of H is
 ! needed
 ! in those cases. 
-!
       
       Saux = dcmplx(0.0_dp, 0.0_dp)
       Haux = dcmplx(0.0_dp, 0.0_dp)
@@ -219,22 +217,21 @@
 !***********************************************************************
 !
 !                 | ------- 1,1 -------     ------- 1,2 ------- |
-!                 | c_{i,up} c_{j,up}^*     c_{i,up} c_{j,down)^* |
-!     D_{i,j} =   |                                             |
+!                 | c_{j,up}^* c_{i,up}     c_{j,up}^* c_{i,dn) |
+!     D_{j,i) =   |                                             |
 !                 | ------- 2,1 -------     ------- 2,2 ------- |
-!                 | c_{i,down} c_{j,up}^*     c_{i,dn} c_{j,dn)^* |
+!                 | c_{j,dn}^* c_{i,up}     c_{j,dn}^* c_{i,dn) |
 !
-!             =   | D_{i,j}(1)              D_{i,j}(3)-i D_{i,j}(4) |
-!                 | D_{i,j}(7)+i D_{i,j}(8) D_{i,j}(2)              |
+!             =   | D_{j,i}(1)              D_{j,i}(3)-i D_{j,i}(4) |
+!                 | D_{j,i}(7)+i D_{j,i}(8) D_{j,i}(2)              |
 
 ! The Energy is computed as E = Tr [ D_{i,j} H_{j,i} ]
 !
 ! The Density Matrix is not "spin box hermitian" even if H does not
-! contain P or
-! L operators. 
+! contain P or L operators. 
 !
 ! Spin-box symmetrization of D inside this present subroutine:
-! D_{i,j}(1,2) = 0.5 ( D_{i,j}(1,2) + D_{i,j}(2,1)^* )
+! D_{j,i}(1,2) = 0.5 ( D_{j,i}(1,2) + D_{j,i}(2,1)^* )
 ! does not affect the results in any case, since H is spin-box
 ! hermitian.
 
@@ -243,7 +240,6 @@
       
       BNode = 0
       iie   = 0
-      
       do ie = 1,2*nuotot
         qe = qo(ie)
         if (Node.eq.BNode) then
@@ -260,7 +256,7 @@
           endif
 #ifdef MPI
           call MPI_Bcast(caux(1,1),2*nuotot,MPI_double_complex,
-     .         BNode,MPI_Comm_World,MPIerror)
+     &         BNode,MPI_Comm_World,MPIerror)
 #endif
           ee = qo(ie) * eo(ie)
           do io = 1,nuo
@@ -270,30 +266,30 @@
               jo = listd(ind)
       
 !                 | ------- 1,1 -------     ------- 2,1 ------- |
-!                 | c_{j,up} c_{i,up}^*     c_{j,dn} c_{i,up)^* |
+!                 | c_{j,up}^* c_{i,up}     c_{j,dn}^* c_{i,up) |
 !     D_{j,i} =   |                                             |
 !                 | ------- 1,2 -------     ------- 2,2 ------- |
-!                 | c_{j,up} c_{i,dn}^*     c_{j,dn} c_{i,dn)^* |
+!                 | c_{j,up}^* c_{i,dn}     c_{j,dn}^* c_{i,dn) |
 !
 !------- 1,1 -----------------------------------------------------------
-              D11  = dconjg(caux(1,iio)) * caux(1,jo)
+              D11 = dconjg(caux(1,jo)) * caux(1,iio)
 !------- 2,2 -----------------------------------------------------------
-              D22  = dconjg(caux(2,iio)) * caux(2,jo)
+              D22 = dconjg(caux(2,jo)) * caux(2,iio)
 !------- 2,1 -----------------------------------------------------------
-              D12  = dconjg(caux(1,iio)) * caux(2,jo)
+              D12 = dconjg(caux(2,jo)) * caux(1,iio)
 !------- 1,2 -----------------------------------------------------------
-              D21  = dconjg(caux(2,iio)) * caux(1,jo)
+              D21 = dconjg(caux(1,jo)) * caux(2,iio)
       
 !------------ Density matrix has to be spin-box hermitian ----------------------
               D11 = 0.5_dp * (D11 + dconjg(D11))
               D22 = 0.5_dp * (D22 + dconjg(D22))
               D12 = 0.5_dp * (D12 + dconjg(D21))
-              D21 = dconjg(D12)
+              !D21 = dconjg(D12)
       
 !       Add contribution to density matrices of unit-cell orbitals
 !       ----------------------------------------------------------------
 !       | D11 = D_{j,i}(1)               D21 = D_{j,i}(3)-i D_{j,i}(4) |
-!       | D12 = D_{i,j}(3)+i D_{i,j}(4)  D22 = D_{j,i}(2)              |
+!       | D12 = D_{j,i}(3)+i D_{j,i}(4)  D22 = D_{j,i}(2)              |
 !       ----------------------------------------------------------------
 !       ----------------------------------------------------------------
 !       | D11 = Dnew(1)                  D21 = Dnew(3)-i Dnew(4)       |
@@ -302,12 +298,12 @@
               Dnew(ind,1) = Dnew(ind,1) + dreal(D11) * qe
               Dnew(ind,2) = Dnew(ind,2) + dreal(D22) * qe
               Dnew(ind,3) = Dnew(ind,3) + dreal(D12) * qe
-              Dnew(ind,4) = Dnew(ind,4) + dimag(D12) * qe
+              Dnew(ind,4) = Dnew(ind,4) - dimag(D12) * qe
       
               Enew(ind,1) = Enew(ind,1) + dreal(D11) * ee
               Enew(ind,2) = Enew(ind,2) + dreal(D22) * ee
               Enew(ind,3) = Enew(ind,3) + dreal(D12) * ee
-              Enew(ind,4) = Enew(ind,4) + dimag(D12) * ee
+              Enew(ind,4) = Enew(ind,4) - dimag(D12) * ee
       
             enddo
           enddo
@@ -322,10 +318,8 @@
  1001 continue
       
 #ifdef DEBUG
-      call write_debug( '    POS diag2k' )
+      call write_debug( '    POS diag2g' )
 #endif
-!***********************************************************************
-      return
-      end
-!**********************************************************************
+
+      end subroutine diag2g
       

=== modified file 'Src/diag2k.F'
--- Src/diag2k.F	2017-08-06 11:17:47 +0000
+++ Src/diag2k.F	2018-01-05 08:40:43 +0000
@@ -37,6 +37,7 @@
 C by J.D. Gale November 2004.
 C Corrected by V.V. Maslyuk for parallel case, Mai 2007
 C Reduced memory requirements, Nick, Aug. 2017
+C Clarified DM calculation, N.R. Papior, Jan. 2018      
 C **************************** INPUT **********************************
 C integer nuo                 : Number of basis orbitals in unit cell
 C integer no                  : Number of basis orbitals in supercell
@@ -119,7 +120,7 @@
       complex(dp), dimension(2,nuotot) :: caux
       logical               getD
       
-!     TEMPOS, INTERNAL VARIABLES etc.
+!     Internal variables .............................................
       
       integer  BNode, BTest, ie, ierror, iie, ik, ind, io, iio
       integer  iuo, j, jo, juo, neigneeded
@@ -128,14 +129,12 @@
 !     Haux(js,juo,is,iuo) = <js,juo|H|is,iuo>
 !     Indices is and js are for spin components
 !     Indices iuo and juo are for orbital components
-      complex(dp)                                 :: cicj
-      complex(dp)                                 :: D11, D22, D12, D21
-      complex(dp)                                 :: kphs
-      real(dp)                                    :: kxij
+      complex(dp) :: cicj, D11, D22, D12, D21
+      complex(dp) :: kphs
+      real(dp) :: kxij
 #ifdef MPI
       integer            :: MPIerror
 #endif
-
       external              cdiag
 
 #ifdef DEBUG
@@ -165,24 +164,11 @@
 !  H_{jo,iuo} = |                                    |
 !               | H_{jo,iuo}^{d,u}  H_{jo,iuo}^{d,d} |
 !
-!            | H(ind,1) + i H(ind,5)   H(ind,3) - i H(ind,4) |
-!          = |                                               |
-!            | H(ind,7) + i H(ind,8)   H(ind,2) + i H(ind,6) |
+!               | H(ind,1) + i H(ind,5)   H(ind,3) - i H(ind,4) |
+!             = |                                               |
+!               | H(ind,7) + i H(ind,8)   H(ind,2) + i H(ind,6) |
 !
 ! 1. Hermiticity imposes H_{i,j}^{is,js}=H_{j,i}^{js,is}^*
-! 2. Since wave functions are real, if there are no single P or L
-! operators:
-!      (a) H_{i,j}^{is,js}=H_{j,i}^{is,js}
-!      (b) These imply spin-box hermiticity:
-!      H_{i,j}^{is,js}=H_{i,j}^{js,is}^*
-!
-!      (c) Hence
-!
-!                               | H(ind,1)                H(ind,3) - i H(ind,4) |
-!           H_{j,i} = H_{i,j} = |                                               |
-!                               | H(ind,3) + i H(ind,4)   H(ind,2)              |
-!
-!
 
 !     Find eigenvalues at every k point
       do ik = 1,nk
@@ -334,37 +320,34 @@
 !***********************************************************************
 !
 !                 | ------- 1,1 -------     ------- 1,2 ------- |
-!                 | c_{i,up} c_{j,up}^*     c_{i,up} c_{j,down)^* |
-!     D_{i,j} =   |                                             |
+!                 | c_{j,up}^* c_{i,up}     c_{j,up}^* c_{i,dn) |
+!     D_{j,i) =   |                                             |
 !                 | ------- 2,1 -------     ------- 2,2 ------- |
-!                 | c_{i,down} c_{j,up}^*     c_{i,dn} c_{j,dn)^* |
+!                 | c_{j,dn}^* c_{i,up}     c_{j,dn}^* c_{i,dn) |
 !
-!             =   | D_{i,j}(1)              D_{i,j}(3)-i D_{i,j}(4) |
-!                 | D_{i,j}(7)+i D_{i,j}(8) D_{i,j}(2)              |
+!             =   | D_{j,i}(1)              D_{j,i}(3)-i D_{j,i}(4) |
+!                 | D_{j,i}(7)+i D_{j,i}(8) D_{j,i}(2)              |
 
 ! The Energy is computed as E = Tr [ D_{i,j} H_{j,i} ]
 !
 ! The Density Matrix is not "spin box hermitian" even if H does not
-! contain P or
-! L operators.
+! contain P or L operators.
 !
 ! Spin-box symmetrization of D inside this present subroutine:
-! D_{i,j}(1,2) = 0.5 ( D_{i,j}(1,2) + D_{i,j}(2,1)^* )
+! D_{j,i}(1,2) = 0.5 ( D_{j,i}(1,2) + D_{j,i}(2,1)^* )
 ! does not affect the results in any case, since H is spin-box
 ! hermitian.
 !
-! The wave function coefficients c_{i,sigma}(ie,k) only satisfy
+! The wave function coefficients c_{j,sigma}(ie,k) only satisfy
 ! "TR"-related 
 ! symmetry relations if the state is collinear (nspin = 2). In this case
-! c_{i,sigma}(ie,k) = c_{i,sigma}(ie, -k)^*
+! c_{j,sigma}(ie,k) = c_{j,sigma}(ie, -k)^*
 !
 ! For a non-collinear state, the above relation does not hold
-! generically. Hence the
-! full BZ must be used to compute the DM.
-
+! generically. Hence the full BZ must be used to compute the DM.
       
-        Dk  = dcmplx(0.0_dp,0.0_dp)
-        Ek  = dcmplx(0.0_dp,0.0_dp)
+        Dk = dcmplx(0.0_dp,0.0_dp)
+        Ek = dcmplx(0.0_dp,0.0_dp)
       
         BNode = 0
         iie = 0
@@ -384,31 +367,26 @@
             endif
 #ifdef MPI
             call MPI_Bcast(caux(1,1),2*nuotot,MPI_double_complex,BNode,
-     .           MPI_Comm_World,MPIerror)
+     &           MPI_Comm_World,MPIerror)
 #endif
-            ee = qo(ie,ik) * eo(ie,ik)
+            ee = qe * eo(ie,ik)
             do iuo = 1,nuo
               call LocalToGlobalOrb(iuo,Node,Nodes,iio)
               do juo = 1,nuotot
-CC RC
-C    In the JF 3.1.1 version the following matrix elements are
-C    calculated by means of dconj(caux(1,juo)) instead of caux(1,juo)
-C    I will leave it as it is. maybe it needs further revision...
-      
 !------- 1,1 -----------------------------------------------------------
-               cicj = caux(1,iio) * dconjg(caux(1,juo))
+               cicj = dconjg(caux(1,juo)) * caux(1,iio)
                Dk(1,juo,1,iuo) = Dk(1,juo,1,iuo) + qe * cicj
                Ek(1,juo,1,iuo) = Ek(1,juo,1,iuo) + ee * cicj
 !------- 2,2 -----------------------------------------------------------
-               cicj = caux(2,iio) * dconjg(caux(2,juo))
+               cicj = dconjg(caux(2,juo)) * caux(2,iio)
                Dk(2,juo,2,iuo) = Dk(2,juo,2,iuo) + qe * cicj
                Ek(2,juo,2,iuo) = Ek(2,juo,2,iuo) + ee * cicj
 !------- 2,1 -----------------------------------------------------------
-               cicj = caux(1,iio) * dconjg(caux(2,juo))
+               cicj = dconjg(caux(2,juo)) * caux(1,iio)
                Dk(1,juo,2,iuo) = Dk(1,juo,2,iuo) + qe * cicj
                Ek(1,juo,2,iuo) = Ek(1,juo,2,iuo) + ee * cicj
 !------- 1,2 -----------------------------------------------------------
-               cicj = caux(2,iio) * dconjg(caux(1,juo))
+               cicj = dconjg(caux(1,juo)) * caux(2,iio)
                Dk(2,juo,1,iuo) = Dk(2,juo,1,iuo) + qe * cicj
                Ek(2,juo,1,iuo) = Ek(2,juo,1,iuo) + ee * cicj
               enddo
@@ -424,7 +402,7 @@
 !       Add contribution to density matrices of unit-cell orbitals
 !       ------------------------------------------------------
 !       | D_{j,i}(1)                 D_{j,i}(3)-i D_{j,i}(4) |
-!       | D_{i,j}(3)+i D_{i,j}(4)    D_{j,i}(2)              |
+!       | D_{j,i}(3)+i D_{j,i}(4)    D_{j,i}(2)              |
 !       ------------------------------------------------------
       
         do iuo = 1,nuo
@@ -433,74 +411,37 @@
           jo = listd(ind)
           juo = indxuo(jo)
           kxij = kpoint(1,ik) * xij(1,ind) +
-     .           kpoint(2,ik) * xij(2,ind) +
-     .           kpoint(3,ik) * xij(3,ind)
-! JFR: the sign in the phase of H_ji is reflected in the DM
+     &           kpoint(2,ik) * xij(2,ind) +
+     &           kpoint(3,ik) * xij(3,ind)
           kphs = cdexp(dcmplx(0.0_dp,-1.0_dp)*kxij)
      
-CC RC  
-C      In the JF's version D11, D22, etc are not time-reversal averaged
-C      I will leave as it is here but maybe it will need major revision
- 
-!         Average k and -k solutions because time-reversal symetry
-!          D11 = 0.5_dp * (Dkc(1,juo,1,iuo) * kphs
-!     .        +    dconjg(Dkc(1,juo,1,iuo) * kphs))
-!      
-!          D22 = 0.5_dp * (Dkc(2,juo,2,iuo) * kphs
-!     .        +    dconjg(Dkc(2,juo,2,iuo) * kphs))
-!          
-!          D12 = 0.5_dp * (Dkc(1,juo,2,iuo) * kphs
-!     .        +    dconjg(Dkc(2,juo,1,iuo) * kphs))
-!      
-!          D21 = 0.5_dp * (Dkc(2,juo,1,iuo) * kphs
-!     .        +    dconjg(Dkc(1,juo,2,iuo) * kphs))
-
           D11 = Dk(1,juo,1,iuo) * kphs
           D22 = Dk(2,juo,2,iuo) * kphs
           D12 = Dk(1,juo,2,iuo) * kphs
           D21 = Dk(2,juo,1,iuo) * kphs
-
       
 ! Make DM Spin-box hermitian 
           D12 = 0.5_dp * (D12 + dconjg(D21))
-          D21 = dconjg(D12)
+          !D21 = dconjg(D12)
 
           Dnew(ind,1) = Dnew(ind,1) + dreal(D11)
           Dnew(ind,2) = Dnew(ind,2) + dreal(D22)
           Dnew(ind,3) = Dnew(ind,3) + dreal(D12)
-          Dnew(ind,4) = Dnew(ind,4) + dimag(D21)
-
-CC RC  
-C      In the JF's version D11, D22, etc are not time-reversal averaged
-C      I will leave as it is here but maybe it will need major revision
-
-!         Average k and -k solutions because time-reversal symetry
-!          D11 = 0.5_dp * (Ekc(1,juo,1,iuo) * kphs
-!     .        +    dconjg(Ekc(1,juo,1,iuo) * kphs))
-!      
-!          D22 = 0.5_dp * (Ekc(2,juo,2,iuo) * kphs
-!     .        +    dconjg(Ekc(2,juo,2,iuo) * kphs))
-!          
-!          D12 = 0.5_dp * (Ekc(1,juo,2,iuo) * kphs
-!     .        +    dconjg(Ekc(2,juo,1,iuo) * kphs))
-!      
-!          D21 = 0.5_dp * (Ekc(2,juo,1,iuo) * kphs
-!     .        +    dconjg(Ekc(1,juo,2,iuo) * kphs))
+          Dnew(ind,4) = Dnew(ind,4) - dimag(D12)
 
           D11 = Ek(1,juo,1,iuo) * kphs
           D22 = Ek(2,juo,2,iuo) * kphs
           D12 = Ek(1,juo,2,iuo) * kphs
           D21 = Ek(2,juo,1,iuo) * kphs
-
-     
-! Make DM Spin-box hermitian  
+      
+! Make EDM Spin-box hermitian 
           D12 = 0.5_dp * (D12 + dconjg(D21))
-          D21 = dconjg(D12)
+          !D21 = dconjg(D12)
 
           Enew(ind,1) = Enew(ind,1) + dreal(D11)
           Enew(ind,2) = Enew(ind,2) + dreal(D22)
-          Enew(ind,3) = Enew(ind,3) + dreal(D12) 
-          Enew(ind,4) = Enew(ind,4) + dimag(D21) 
+          Enew(ind,3) = Enew(ind,3) + dreal(D12)
+          Enew(ind,4) = Enew(ind,4) - dimag(D12)
 
         enddo
         enddo
@@ -508,27 +449,9 @@
 
  1001 continue
       
-!        write(160+node,*)'ITERATION'
-!        write(170+node,*)'ITERATION'
-!        do iuo = 1,nuo
-!          do j = 1,numd(iuo)
-!            ind = listdptr(iuo) + j
-!            jo = listd(ind)
-!            juo = indxuo(jo)
-!             write(160+node,'(3i5,4f15.9)')iuo,juo,jo,Dnew(ind,1:4)
-!             write(170+node,'(3i5,4f15.9)')iuo,juo,jo,
-!     .    (Dnew(ind,1)+Dnew(ind,2))*0.5,
-!     .    (Dnew(ind,1)-Dnew(ind,2))*0.5,Dnew(ind,3:4)
-!          enddo
-!        enddo
-!      
-!      call die('maslyuk')
-!***********************************************************************
 #ifdef DEBUG
       call write_debug( '    POS diag2k' )
 #endif
 
-      return
       end subroutine diag2k
-!***********************************************************************      
       

=== modified file 'version.info'
--- version.info	2017-12-21 15:51:22 +0000
+++ version.info	2018-01-05 08:40:43 +0000
@@ -1,3 +1,3 @@
-siesta-4.1--829--spin-fix-2
+siesta-4.1--829--spin-fix-2-nicpa-1