scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Self interaction correction (SIC) 

tentatively added to FPMD,
work still in progress


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@758 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
cavazzon 2004-04-01 16:30:59 +00:00
parent e9575b3b67
commit a3e02b9d06
6 changed files with 179 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
Modules/control_flags.f90
View File

@ -52,7 +52,7 @@
tcp, tcap, tnodump, tdamp, tdampions, tconvthrs, &
convergence_criteria, tionstep, nstepe, tsteepdesc, &
ionic_conjugate_gradient, tconjgrad_ion, &
tatomicwfc, tscreen, gamma_only, ekin_conv_thr, ekin_maxiter
tatomicwfc, tscreen, gamma_only, ekin_conv_thr, ekin_maxiter, force_pairing
PUBLIC :: fix_dependencies, check_flags
@ -99,6 +99,10 @@
LOGICAL :: tscreen = .FALSE. ! Use screened coulomb potentials for cluster calculations
LOGICAL :: twfcollect = .FALSE. ! Collect wave function in the restart file at the end of run.
! ... Force pairing
LOGICAL :: force_pairing
TYPE (convergence_criteria) :: tconvthrs
! thresholds used to check GS convergence

23
Modules/energies.f90
View File

@ -45,16 +45,23 @@
REAL(dbl) :: EXC = 0.0_dbl
REAL(dbl) :: VXC = 0.0_dbl
REAL(dbl) :: EBAND = 0.0_dbl
REAL(dbl) :: self_ehte = 0.0d0
REAL(dbl) :: self_sxc = 0.0d0
REAL(dbl) :: self_vxc = 0.0d0
CONTAINS
! ---------------------------------------------------------------------------- !
SUBROUTINE total_energy(edft, omega, eexc, vvxc, eh, eps, nnr)
SUBROUTINE total_energy(edft, omega, eexc, vvxc, eh, eps, &
self_ehte_in, self_sxc_in, self_vxc_in, nnr)
TYPE (dft_energy_type) :: edft
REAL(dbl), INTENT(IN) :: OMEGA, EEXC, VVXC
REAL(dbl) :: VXC
REAL(dbl) :: self_sxc_in, self_vxc_in
COMPLEX(dbl) :: self_ehte_in
COMPLEX(dbl), INTENT(IN) :: EH, EPS
INTEGER, INTENT(IN) :: nnr
@ -65,6 +72,10 @@
esr = edft%esr
ekin = edft%ekin
self_ehte = REAL( self_ehte_in )
self_sxc = self_sxc_in
self_vxc = self_vxc_in
EXC = EEXC * omega / REAL(NNR)
VXC = VVXC * omega / REAL(NNR)
edft%exc = exc
@ -112,7 +123,8 @@
TYPE (dft_energy_type), OPTIONAL, INTENT(IN) :: edft
IF( PRESENT ( edft ) ) THEN
WRITE( stdout,1) edft%ETOT, edft%EKIN, edft%EHT, edft%ESELF, edft%ESR, &
edft%EPSEU, edft%ENL, edft%EXC, edft%EVDW, edft%emkin
edft%EPSEU, edft%ENL, edft%EXC, edft%EVDW, edft%emkin, &
self_ehte, self_sxc, self_vxc, edft%ehti
ELSE
WRITE( stdout,1) ETOT, EKIN, EHT, ESELF, ESR, EPSEU, ENL, EXC, EVDW
END IF
@ -125,7 +137,12 @@
,6X,' N-L PSEUDOPOTENTIAL ENERGY = ',F18.10,' A.U.'/ &
,6X,' EXCHANGE-CORRELATION ENERGY = ',F18.10,' A.U.'/ &
,6X,' VAN DER WAALS ENERGY = ',F18.10,' A.U.'/ &
,6X,' EMASS KINETIC ENERGY = ',F18.10,' A.U.'/,/)
,6X,' EMASS KINETIC ENERGY = ',F18.10,' A.U.'/ &
,6X,' SIC HARTREE ELECTRON ENERGY = ',F18.10,' A.U.'/ &
,6X,' SIC EXCHANGE-CORRELA ENERGY = ',F18.10,' A.U.'/ &
,6X,' SIC EXCHANGE-CORRELA POTENT = ',F18.10,' A.U.'/ &
,6X,' IONIC ELECTROSTATIC ENERGY = ',F18.10,' A.U.'/,/)
RETURN
END SUBROUTINE print_energies

41
Modules/input_parameters.f90
View File

@ -419,13 +419,44 @@ MODULE input_parameters
INTEGER :: report = 1
CHARACTER(LEN=80) :: sic = 'none'
! sic = 'none' | 'sic_pz' | 'sic_mac' | 'only_sich' | 'only_sicxc_pz' | 'only_sicxc_mac' |
! Where:
! 'none' => self_interaction == 0 -> no USIC
! 'sic_pz' => self_interaction == 1 -> USIC = Exc[rhoup-rhodown,0] + Uhartree[rhoup-rhodown]
!! proposed by Perdew, Zunger, PRB 23 (1981) 5048
! 'sic_mac' => self_interaction == 2 -> USIC = Exc[rhoup,rhodown] - Exc[rhopaired, rhopaired] + Uhartree[rhoup-rhodown]
!! proposed by Lundin-Eriksson, IJQC 81 (2003) 247
!! implemented by Mauri-Calandra-d'Avezac (2003/04) for one electron/hole
! 'only_sich' => self_interaction == 3 -> USIC = Uhartree[rhoup-rhodown]
!!SIcorrection only on hartree part )
! 'only_sich_pz' => self_interaction == -1 -> USIC = Exc[rhoup-rhodown,0]
!!SIcorrection only on xc by PZ
! 'only_sich_mac' => self_interaction == -2 -> USIC = Exc[rhoup,rhodown] - Exc[rhopaired, rhopaired]
!!sIcorrection only on xc by MAC:
!!rhopaired==rhodown since that the charge more is defined in the highest level and up
!!rhounpaired==rhoup-rhodown
REAL(dbl) :: sic_epsilon = 1.0d0
LOGICAL :: force_pairing= .FALSE.
! FORCEPAIRING
! paires electrons forcedly in two spin channels calculation
! works for stuctures with at most one unpaired eletron
! just a comment: nelu is the number of el with spin up
! while neld== number of el with spin down
! define the unpaired el with spin up
NAMELIST / system / ibrav, celldm, a, b, c, cosab, cosac, cosbc, nat,&
ntyp, nbnd, nelec, ecutwfc, ecutrho, nr1, nr2, nr3, nr1s, nr2s, &
nr3s, nr1b, nr2b, nr3b, nosym, starting_magnetization, &
occupations, degauss, ngauss, nelup, neldw, nspin, ecfixed, &
qcutz, q2sigma, xc_type, lda_plus_U, Hubbard_U, Hubbard_alpha, &
edir, emaxpos, eopreg, eamp, smearing, starting_ns_eigenvalue, &
noncolin, mcons, lambda, i_cons, angle1, angle2, report
noncolin, mcons, lambda, i_cons, angle1, angle2, report, &
sic, sic_epsilon, force_pairing
!=----------------------------------------------------------------------------=!
@ -889,13 +920,15 @@ MODULE input_parameters
REAL(KIND=DP) :: w_1 = 0.5D-1
REAL(KIND=DP) :: w_2 = 0.5D0
REAL(KIND=DP) :: sic_rloc = 0.0d0
NAMELIST / ions / ion_dynamics, ion_radius, ion_damping, ion_positions, &
ion_velocities, ion_temperature, tempw, fnosep, tranp, amprp, greasp, &
tolp, ion_nstepe, ion_maxstep, upscale, potential_extrapolation, &
num_of_images, CI_scheme, VEC_scheme, minimization_scheme, &
optimization, damp, temp_req, ds, k_max, k_min, neb_thr, &
trust_radius_max, trust_radius_min, trust_radius_ini, trust_radius_end, &
w_1, w_2, lbfgs_ndim
w_1, w_2, lbfgs_ndim, sic_rloc
!
!=----------------------------------------------------------------------------=!
@ -1032,6 +1065,7 @@ MODULE input_parameters
REAL(dbl) :: rd_pos(3,natx) = 0.0d0 ! unsorted position from input
INTEGER :: sp_pos(natx) = 0
INTEGER :: if_pos(3,natx) = 1
INTEGER :: id_loc(natx) = 0
INTEGER :: na_inp(nsx) = 0 ! number of atom for each specie
LOGICAL :: tapos = .FALSE.
LOGICAL :: tscal = .TRUE.
@ -1045,6 +1079,8 @@ MODULE input_parameters
REAL (KIND=DP), ALLOCATABLE :: pos(:,:)
!
! ION_VELOCITIES
!
@ -1184,6 +1220,7 @@ MODULE input_parameters
LOGICAL :: tprnrho = .FALSE.
!
! CLIMBING_IMAGES
!

47
Modules/ions_base.f90
View File

@ -52,6 +52,13 @@
INTEGER, ALLOCATABLE :: if_pos(:,:)
INTEGER :: fixatom !!! to be removed
INTEGER :: ind_localisation(natx) = 0 ! true if we want to know the localization arount the atom
INTEGER :: nat_localisation = 0
LOGICAL :: print_localisation = .FALSE. ! Calculates hartree energy around specified atoms
INTEGER :: self_interaction = 0
REAL(dbl) :: si_epsilon = 0.0d0
REAL(dbl) :: rad_localisation = 0.0d0
REAL(dbl), ALLOCATABLE :: pos_localisation(:,:)
LOGICAL :: tions_base_init = .FALSE.
@ -96,7 +103,7 @@
SUBROUTINE ions_base_init( nsp_ , nat_ , na_ , ityp_ , tau_ , amass_ , &
atm_ , if_pos_ )
atm_ , if_pos_ , id_loc_ , sic_ , sic_epsilon_, sic_rloc_ )
USE constants, ONLY: scmass
IMPLICIT NONE
INTEGER, INTENT(IN) :: nsp_ , nat_ , na_ (:) , ityp_ (:)
@ -104,6 +111,10 @@
REAL(dbl), INTENT(IN) :: amass_(:)
CHARACTER(LEN=*), INTENT(IN) :: atm_ (:)
INTEGER, INTENT(IN) :: if_pos_ (:,:)
INTEGER, OPTIONAL, INTENT(IN) :: id_loc_ (:)
CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sic_
REAL(dbl), OPTIONAL, INTENT(IN) :: sic_epsilon_
REAL(dbl), OPTIONAL, INTENT(IN) :: sic_rloc_
INTEGER :: i, ia
nsp = nsp_
@ -139,6 +150,39 @@
if_pos = 1
if_pos(:,:) = if_pos_ (:,1:nat)
IF( PRESENT( sic_ ) ) THEN
select case ( TRIM( sic_ ) )
case ( 'sic_pz' )
self_interaction = 1
case ( 'sic_mac' )
self_interaction = 2
case ( 'only_sich' )
self_interaction = 3
case ( 'only_sicxc_pz' )
self_interaction = -1
case ( 'only_sicxc_mac' )
self_interaction = -2
case default
self_interaction = 0
end select
END IF
IF( PRESENT( sic_epsilon_ ) ) THEN
si_epsilon = sic_epsilon_
END IF
IF( PRESENT( sic_rloc_ ) ) THEN
rad_localisation = sic_rloc_
END IF
IF( PRESENT( id_loc_ ) ) THEN
ind_localisation(1:nat) = id_loc_ ( 1:nat )
nat_localisation = COUNT( ind_localisation > 0 )
ALLOCATE( pos_localisation( 4, nat_localisation ) )
!counting the atoms around which i want to calculate the charge localization
ELSE
ind_localisation(1:nat) = 0
nat_localisation = 0
END IF
!
IF( nat_localisation > 0 ) print_localisation = .TRUE.
!
! ... TEMP: calculate fixatom (to be removed)
!
@ -166,6 +210,7 @@
IF( ALLOCATED( ityp ) ) DEALLOCATE( ityp )
IF( ALLOCATED( tau ) ) DEALLOCATE( tau )
IF( ALLOCATED( if_pos ) ) DEALLOCATE( if_pos )
IF( ALLOCATED( pos_localisation ) ) DEALLOCATE( pos_localisation )
RETURN
END SUBROUTINE

28
Modules/read_cards.f90
View File

@ -253,9 +253,7 @@ MODULE read_cards_module
ELSE IF ( TRIM(card) == 'CLIMBING_IMAGES' ) THEN
!
CALL card_climbing_images( input_line )
IF ( ( prog == 'FPMD' .OR. prog == 'CP' ) .AND. ionode ) &
WRITE( stdout,'(A)') 'Warning: card '//trim(input_line)//' ignored'
!
ELSE
!
IF ( ionode ) &
@ -557,6 +555,9 @@ MODULE read_cards_module
CALL read_line( input_line )
CALL field_count( nfield, input_line )
!
IF( sic /= 'none' .AND. nfield /= 8 ) &
CALL errore(' read_cards ', ' ATOMIC_POSITIONS with sic, 8 columns required ', 1 )
!
IF ( nfield == 4 ) THEN
!
READ(input_line,*) lb_pos, ( rd_pos(k,ia), k = 1, 3 )
@ -570,6 +571,16 @@ MODULE read_cards_module
if_pos(2,ia), &
if_pos(3,ia)
!
ELSE IF ( nfield == 8 ) THEN
!
READ(input_line,*) lb_pos, rd_pos(1,ia), &
rd_pos(2,ia), &
rd_pos(3,ia), &
if_pos(1,ia), &
if_pos(2,ia), &
if_pos(3,ia), &
id_loc(ia)
!
ELSE
!
CALL errore( ' read_cards ', ' wrong number of columns ' // &
@ -1720,11 +1731,9 @@ MODULE read_cards_module
!
! index1, ..., indexN are indeces of the images that have to climb
!
!
!
!------------------------------------------------------------------------
! BEGIN manual
!----------------------------------------------------------------------
! END manual
!------------------------------------------------------------------------
!
SUBROUTINE card_climbing_images( input_line )
!
@ -1779,6 +1788,11 @@ MODULE read_cards_module
END SUBROUTINE card_climbing_images
!
!
!------------------------------------------------------------------------
! BEGIN manual
!----------------------------------------------------------------------
!
!
!
! TEMPLATE
!

52
Modules/read_namelists.f90
View File

@ -265,6 +265,10 @@ MODULE read_namelists_module
diago_diis_ndim = 3
diago_diis_ethr = 1.0D-3
diago_diis_keep = .FALSE.
sic = 'none'
sic_epsilon = 1.d0
force_pairing = .false.
!
RETURN
!
@ -328,6 +332,8 @@ MODULE read_namelists_module
k_min = 0.1D0
neb_thr = 0.05D0
!
sic_rloc = 0.0d0
!
RETURN
!
END SUBROUTINE
@ -589,6 +595,9 @@ MODULE read_namelists_module
CALL mp_bcast( diago_diis_ndim, ionode_id )
CALL mp_bcast( diago_diis_ethr, ionode_id )
CALL mp_bcast( diago_diis_keep, ionode_id )
CALL mp_bcast( sic, ionode_id )
CALL mp_bcast( sic_epsilon , ionode_id )
CALL mp_bcast( force_pairing , ionode_id )
!
RETURN
!
@ -640,6 +649,8 @@ MODULE read_namelists_module
CALL mp_bcast( k_max, ionode_id )
CALL mp_bcast( k_min, ionode_id )
CALL mp_bcast( neb_thr, ionode_id )
CALL mp_bcast( sic_rloc, ionode_id )
!
RETURN
!
@ -885,9 +896,32 @@ MODULE read_namelists_module
IF ( i_cons < 0 .OR. i_cons > 2 ) &
CALL errore( sub_name ,' wrong i_cons ', -1 )
END IF
!
!!control on SIC variables
IF (sic /='none' ) then
IF (sic_epsilon > 1.d0 ) &
CALL errore( sub_name, &
& ' invalid sic_epsilon, greater than 1.',-1)
IF (sic_epsilon < 0.d0 ) &
CALL errore( sub_name, &
& ' invalid sic_epsilon, less than 0 ',-1)
IF ( force_pairing == .false.) &
CALL errore( sub_name, &
& ' invalid force_pairing with sic activated', -1)
IF ( nspin /= 2) &
CALL errore( sub_name, &
& ' invalid nspin with sic activated', -1)
IF ( (nelup /= 0) .or. (neldw /= 0) ) &
CALL errore( sub_name, &
& ' invalid nelup and neldwn spin with sic activated', -1)
IF ( nelup /= (neldw + 1) ) &
CALL errore( sub_name, &
& ' invalid nelup /= (neldwn +1) spin with sic activated', -1)
ENDIF !! on sic_control
RETURN
!
END SUBROUTINE
!
!=----------------------------------------------------------------------=!
@ -947,9 +981,9 @@ MODULE read_namelists_module
IF( empty_states_ethr < 0.0d0 ) &
CALL errore( sub_name, &
& ' invalid empty_states_ethr, less than 0 ',-1)
!
!
RETURN
!
END SUBROUTINE
!
!=----------------------------------------------------------------------=!
@ -1022,6 +1056,8 @@ MODULE read_namelists_module
IF ( .NOT. allowed ) &
CALL errore( sub_name, ' minimization_scheme '''// &
& TRIM( minimization_scheme )//''' not allowed ', -1 )
IF (sic /= 'none' .and. sic_rloc == 0.d0) &
CALL errore( sub_name, ' invalid sic_rloc with sic activated ', -1)
!
RETURN
!
@ -1185,6 +1221,7 @@ MODULE read_namelists_module
CALL errore( sub_name,' calculation '// &
& calculation//' not implemented ',1)
END SELECT
!
IF ( prog == 'PW' ) THEN
!
@ -1208,6 +1245,13 @@ MODULE read_namelists_module
END IF
!
END IF
IF ( TRIM( sic ) /= 'none' ) THEN
IF ( nspin == 2 .and. nelec > 1 .and. ( nelup == neldw .or. nelup == (neldw+1) ) ) THEN
force_pairing = .true.
END IF
END IF
!
RETURN
!