1!----------------------------------------------------------------------- 2 SUBROUTINE spinsq (c,bec,rhor) 3!----------------------------------------------------------------------- 4! 5! estimate of <S^2>=s(s+1) in two different ways. 6! 1) using as many-body wavefunction a single Slater determinant 7! constructed with Kohn-Sham orbitals: 8! 9! <S^2> = (Nup-Ndw)/2 * (Nup-Ndw)/2+1) + Ndw - 10! \sum_up\sum_dw < psi_up | psi_dw > 11! 12! where Nup, Ndw = number of up and down states, the sum is over 13! occupied states. Not suitable for fractionary occupancy. 14! In the ultrasoft scheme (c is the smooth part of \psi): 15! 16! < psi_up | psi_dw > = \sum_G c*_up(G) c_dw(G) + 17! \int Q_ij <c_up|beta_i><beta_j|c_dw> 18! 19! This is the usual formula, unsuitable for fractionary occupancy. 20! 2) using the "LSD model" of Wang, Becke, Smith, JCP 102, 3477 (1995): 21! 22! <S^2> = (Nup-Ndw)/2 * (Nup-Ndw)/2+1) + Ndw - 23! \int max(rhoup(r),rhodw(r)) dr 24! 25! Requires on input: c=psi, bec=<c|beta>, rhoup(r), rhodw(r) 26! Assumes real psi, with only half G vectors. 27! 28 USE kinds, only: dp 29 USE electrons_base, ONLY: nx => nbspx, n => nbsp, iupdwn, nupdwn, f, nel, nspin 30 USE io_global, ONLY: stdout 31 USE mp_global, ONLY: intra_bgrp_comm 32 USE mp, ONLY: mp_sum 33 USE fft_base, ONLY: dfftp 34 USE gvecw, ONLY: ngw 35 USE gvect, ONLY: gstart 36 USE cell_base, ONLY: omega 37 USE uspp, ONLY: nkb, nkbus, qq_nt, indv_ijkb0 38 USE uspp_param, ONLY: nh, upf 39 USE ions_base, ONLY: na, nat, ityp 40! 41 IMPLICIT NONE 42! input 43 REAL(dp) :: bec(nkb,n), rhor(dfftp%nnr,nspin) 44 COMPLEX(dp):: c(ngw,nx) 45! local variables 46 INTEGER :: nup, ndw, ir, i, j, jj, ig, ia, is, iv, jv, inl, jnl 47 REAL(dp) :: spin0, spin1, spin2, fup, fdw 48 REAL(dp), ALLOCATABLE:: overlap(:,:), temp(:) 49 LOGICAL :: frac 50! 51! 52 IF (nspin.EQ.1) RETURN 53! 54! find spin-up and spin-down states 55! 56 fup = 0.0d0 57 DO i=iupdwn(1),nupdwn(1) 58 fup = fup + f(i) 59 END DO 60 nup = NINT(fup) 61 ndw = nel(1)+nel(2) - nup 62! 63! paranoid checks 64! 65 frac= ABS(fup-nup).GT.1.0d-6 66 fup = 0.0d0 67 DO i=1,nup 68 fup = fup + f(i) 69 END DO 70 frac=frac.OR.ABS(fup-nup).GT.1.0d-6 71 fdw = 0.0d0 72 DO j=iupdwn(2),iupdwn(2)-1+ndw 73 fdw = fdw + f(j) 74 END DO 75 frac=frac.OR.ABS(fdw-ndw).GT.1.0d-6 76! 77 spin0 = ABS(fup-fdw)/2.d0 * ( ABS(fup-fdw)/2.d0 + 1.d0 ) + fdw 78! 79! Becke's formula for spin polarization 80! 81 spin1 = 0.0d0 82 DO ir=1,dfftp%nnr 83 spin1 = spin1 - MIN(rhor(ir,1),rhor(ir,2)) 84 END DO 85 CALL mp_sum( spin1, intra_bgrp_comm ) 86 spin1 = spin0 + omega/(dfftp%nr1*dfftp%nr2*dfftp%nr3)*spin1 87 IF (frac) THEN 88 WRITE( stdout,'(/" Spin contamination: s(s+1)=",f5.2," (Becke) ",& 89 & f5.2," (expected)")') & 90 & spin1, ABS(fup-fdw)/2.d0*(ABS(fup-fdw)/2.d0+1.d0) 91 RETURN 92 END IF 93! 94! Slater formula, smooth contribution to < psi_up | psi_dw > 95! 96 ALLOCATE (overlap(nup,ndw)) 97 ALLOCATE (temp(ngw)) 98 DO j=1,ndw 99 jj=j+iupdwn(2)-1 100 DO i=1,nup 101 overlap(i,j)=0.d0 102 DO ig=1,ngw 103 temp(ig)=DBLE(CONJG(c(ig,i))*c(ig,jj)) 104 END DO 105 overlap(i,j) = 2.d0*SUM(temp) 106 IF (gstart == 2) overlap(i,j) = overlap(i,j) - temp(1) 107 END DO 108 END DO 109 DEALLOCATE (temp) 110 CALL mp_sum( overlap, intra_bgrp_comm ) 111 DO j=1,ndw 112 jj=j+iupdwn(2)-1 113 DO i=1,nup 114! 115! vanderbilt contribution to < psi_up | psi_dw > 116! 117 DO ia=1,nat 118 is=ityp(ia) 119 IF( upf(is)%tvanp ) THEN 120 DO iv=1,nh(is) 121 DO jv=1,nh(is) 122 IF(ABS(qq_nt(iv,jv,is)).GT.1.e-5) THEN 123 inl = indv_ijkb0(ia) + iv 124 jnl = indv_ijkb0(ia) + jv 125 overlap(i,j) = overlap(i,j) + qq_nt(iv,jv,is)*bec(inl,i)*bec(jnl,jj) 126 ENDIF 127 END DO 128 END DO 129 END IF 130 END DO 131 END DO 132 END DO 133! 134 spin2 = spin0 135 DO j=1,ndw 136 DO i=1,nup 137 spin2 = spin2 - overlap(i,j)**2 138 END DO 139 END DO 140! 141 DEALLOCATE (overlap) 142! 143 WRITE( stdout,'(/" Spin contamination: s(s+1)=",f5.2," (Slater) ", & 144 & f5.2," (Becke) ",f5.2," (expected)")') & 145 & spin2,spin1, ABS(fup-fdw)/2.d0*(ABS(fup-fdw)/2.d0+1.d0) 146! 147 RETURN 148 END SUBROUTINE spinsq 149 150 151