2003-01-20 05:58:50 +08:00
|
|
|
!
|
|
|
|
! Copyright (C) 2002 PWSCF group
|
|
|
|
! This file is distributed under the terms of the
|
|
|
|
! GNU General Public License. See the file `License'
|
|
|
|
! in the root directory of the present distribution,
|
|
|
|
! or http://www.gnu.org/copyleft/gpl.txt .
|
|
|
|
!
|
2004-06-26 01:25:37 +08:00
|
|
|
#include "f_defs.h"
|
2003-01-20 05:58:50 +08:00
|
|
|
!
|
|
|
|
!----------------------------------------------------------------------
|
2004-11-04 21:35:00 +08:00
|
|
|
SUBROUTINE stres_hub ( sigmah )
|
2003-01-20 05:58:50 +08:00
|
|
|
!----------------------------------------------------------------------
|
|
|
|
!
|
|
|
|
! This routines computes the Hubbard contribution to the internal stress
|
2003-02-08 00:04:36 +08:00
|
|
|
! tensor. It gives in output the array sigmah(i,j) which corresponds to
|
2003-01-20 05:58:50 +08:00
|
|
|
! the quantity -(1/\Omega)dE_{h}/d\epsilon_{i,j}
|
|
|
|
!
|
2004-06-12 21:44:18 +08:00
|
|
|
USE kinds, ONLY : DP
|
|
|
|
USE ions_base, ONLY : nat, ityp
|
|
|
|
USE cell_base, ONLY : omega, at, bg
|
|
|
|
USE ldaU, ONLY : hubbard_lmax, hubbard_l, hubbard_u, &
|
|
|
|
hubbard_alpha, ns, U_projection
|
|
|
|
USE lsda_mod, ONLY : nspin
|
|
|
|
USE symme, ONLY : s, nsym
|
|
|
|
USE io_files, ONLY : prefix, iunocc
|
|
|
|
USE wvfct, ONLY : gamma_only
|
2004-11-04 21:35:00 +08:00
|
|
|
USE io_global, ONLY : stdout, ionode
|
2004-06-12 21:44:18 +08:00
|
|
|
!
|
2004-11-04 21:35:00 +08:00
|
|
|
IMPLICIT NONE
|
2004-06-12 21:44:18 +08:00
|
|
|
!
|
2005-08-28 22:09:42 +08:00
|
|
|
REAL (DP) :: sigmah(3,3) ! output: the Hubbard stresses
|
2003-01-20 05:58:50 +08:00
|
|
|
|
2004-11-04 21:35:00 +08:00
|
|
|
INTEGER :: ipol, jpol, na, nt, is,isi, m1,m2,m3,m4
|
|
|
|
INTEGER :: ldim
|
2005-08-28 22:09:42 +08:00
|
|
|
REAL (DP) :: omin1, current_sum, inverse_sum, sum, temp, flag
|
2004-11-04 21:35:00 +08:00
|
|
|
LOGICAL :: exst
|
2005-08-28 22:09:42 +08:00
|
|
|
REAL (DP), ALLOCATABLE :: dns(:,:,:,:)
|
2003-07-31 20:54:09 +08:00
|
|
|
! dns(ldim,ldim,nspin,nat), ! the derivative of the atomic occupations
|
2003-02-10 16:58:33 +08:00
|
|
|
|
2004-11-04 21:35:00 +08:00
|
|
|
IF (U_projection .NE. "atomic") CALL errore("stres_hub", &
|
2004-04-03 00:05:17 +08:00
|
|
|
" stress for this U_projection_type not implemented",1)
|
|
|
|
|
2004-11-04 21:35:00 +08:00
|
|
|
IF (gamma_only) CALL errore('stres_hub',&
|
2004-02-14 04:51:21 +08:00
|
|
|
' LDA+U, stress AND gamma-only not implemented yet',1)
|
|
|
|
|
2003-01-20 05:58:50 +08:00
|
|
|
sigmah(:,:) = 0.d0
|
|
|
|
|
2003-02-10 16:58:33 +08:00
|
|
|
ldim = 2 * Hubbard_lmax + 1
|
2004-11-04 21:35:00 +08:00
|
|
|
ALLOCATE (dns(ldim,ldim,nspin,nat))
|
2003-01-20 05:58:50 +08:00
|
|
|
dns(:,:,:,:) = 0.d0
|
|
|
|
|
2004-11-04 21:35:00 +08:00
|
|
|
IF ( ionode ) THEN
|
|
|
|
|
2005-06-17 21:27:38 +08:00
|
|
|
CALL seqopn(iunocc,'occup','formatted',exst)
|
2004-11-04 21:35:00 +08:00
|
|
|
READ(iunocc,*) ns
|
|
|
|
CLOSE(unit=iunocc,status='keep')
|
|
|
|
|
|
|
|
END IF
|
|
|
|
|
2003-01-20 05:58:50 +08:00
|
|
|
|
|
|
|
#ifdef DEBUG
|
2004-11-04 21:35:00 +08:00
|
|
|
DO na=1,nat
|
|
|
|
DO is=1,nspin
|
2003-01-20 05:58:50 +08:00
|
|
|
nt = ityp(na)
|
2004-11-04 21:35:00 +08:00
|
|
|
IF (Hubbard_U(nt).NE.0.d0.OR.Hubbard_alpha(nt).NE.0.d0) THEN
|
2003-11-04 18:53:05 +08:00
|
|
|
WRITE( stdout,'(a,2i3)') 'NS(NA,IS) ', na,is
|
2004-11-04 21:35:00 +08:00
|
|
|
DO m1=1,ldim
|
2003-11-04 18:53:05 +08:00
|
|
|
WRITE( stdout,'(7f10.4)') (ns(m1,m2,is,na),m2=1,ldim)
|
2004-11-04 21:35:00 +08:00
|
|
|
END DO
|
|
|
|
END IF
|
|
|
|
END DO
|
|
|
|
END DO
|
2003-01-20 05:58:50 +08:00
|
|
|
#endif
|
|
|
|
omin1 = 1.d0/omega
|
2004-11-04 21:35:00 +08:00
|
|
|
DO ipol = 1,3
|
|
|
|
DO jpol = 1,ipol
|
|
|
|
CALL dndepsilon(dns,ldim,ipol,jpol)
|
|
|
|
DO na = 1,nat
|
2003-01-20 05:58:50 +08:00
|
|
|
nt = ityp(na)
|
2004-11-04 21:35:00 +08:00
|
|
|
IF (Hubbard_U(nt).NE.0.d0.OR.Hubbard_alpha(nt).NE.0.d0) THEN
|
|
|
|
DO is = 1,nspin
|
2003-01-20 05:58:50 +08:00
|
|
|
#ifdef DEBUG
|
2003-11-04 18:53:05 +08:00
|
|
|
WRITE( stdout,'(a,4i3)') 'DNS(IPOL,JPOL,NA,IS) ', ipol,jpol,na,is
|
|
|
|
WRITE( stdout,'(5f10.4)') ((dns(m1,m2,is,na),m2=1,5),m1=1,5)
|
2003-01-20 05:58:50 +08:00
|
|
|
#endif
|
2004-11-04 21:35:00 +08:00
|
|
|
DO m2 = 1, 2 * Hubbard_l(nt) + 1
|
2003-01-20 05:58:50 +08:00
|
|
|
sigmah(ipol,jpol) = sigmah(ipol,jpol) - omin1 * &
|
2003-07-31 20:54:09 +08:00
|
|
|
Hubbard_U(nt) * 0.5d0 * dns(m2,m2,is,na)
|
2004-11-04 21:35:00 +08:00
|
|
|
DO m1 = 1, 2 * Hubbard_l(nt) + 1
|
2003-01-20 05:58:50 +08:00
|
|
|
sigmah(ipol,jpol) = sigmah(ipol,jpol) + omin1 * &
|
2003-07-31 20:54:09 +08:00
|
|
|
Hubbard_U(nt) * ns(m2,m1,is,na) * dns(m1,m2,is,na)
|
2004-11-04 21:35:00 +08:00
|
|
|
END DO
|
|
|
|
END DO
|
|
|
|
END DO
|
|
|
|
END IF
|
|
|
|
END DO
|
|
|
|
END DO
|
|
|
|
END DO
|
|
|
|
IF (nspin.EQ.1) sigmah(:,:) = 2.d0 * sigmah(:,:)
|
2003-03-20 19:21:34 +08:00
|
|
|
|
2003-01-20 05:58:50 +08:00
|
|
|
!
|
|
|
|
! Symmetryze the stress tensor
|
|
|
|
!
|
2004-11-04 21:35:00 +08:00
|
|
|
DO ipol = 1,3
|
|
|
|
DO jpol = ipol,3
|
2003-01-20 05:58:50 +08:00
|
|
|
sigmah(ipol,jpol) = sigmah(jpol,ipol)
|
2004-11-04 21:35:00 +08:00
|
|
|
END DO
|
|
|
|
END DO
|
2003-03-20 19:21:34 +08:00
|
|
|
|
2004-11-04 21:35:00 +08:00
|
|
|
CALL trntns(sigmah,at,bg,-1)
|
|
|
|
CALL symtns(sigmah,nsym,s)
|
|
|
|
CALL trntns(sigmah,at,bg,1)
|
2003-01-20 05:58:50 +08:00
|
|
|
|
2004-11-04 21:35:00 +08:00
|
|
|
DEALLOCATE (dns)
|
2003-01-20 05:58:50 +08:00
|
|
|
|
2004-11-04 21:35:00 +08:00
|
|
|
RETURN
|
|
|
|
END SUBROUTINE stres_hub
|