1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
|
!
! Copyright (C) 2002 FPMD 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 .
!
!------------------------------------------------------------------------------!
MODULE sic_module
!------------------------------------------------------------------------------!
!
! The versions after 3.0 contain also the self-interaction-correction method
! has proposed by Mauri et al. (PRB 2005), taking also into account the 'comment'
! proposed by Sprik et al. (ICR 2005).
! Thus, we introduce the parameters sic_alpha and sic_epsilon to correct the
! exchange-correlation and the electronic hartree potentials, respectively.
! They are two empirical parameters, thus to remain in a ab-initio
! set them equal to 1.0_DP.
! Sprik et al. showed that, in same cases, i.e. OH radical, it should be better
! to under estimate the correction to ex-ch, since in same way the exch already
! corrects the electronic hartree part.
! HOW AND WHEN USE THE SIC::
! Fran's personal considerations:
! the SIC is a way to correct the self-interaction WHEN
! ONE and only ONE e- lives in an unpaired electronic level
! we have choosen for it the spin up
! Remember to select nspin == 2 and nelup = neldw + 1
! the other e- are fictitious calculate in a LSD approach:
! infact, even if the paired e- feel a different potential (for spin up and spin dw)
! we constrain them to have the same force, and the same eigenvalues, the same eigenstates
! When you applied this SIC scheme to a molecule or to an atom, which are neutral,
! remember that you have to consider another correction to the energy level as proposed
! by Landau: infact if you start from a neutral system and subtract the self-intereaction
! the unpaired e- feels a charge system. Thus remeber a correction term ~2.317(Madelung)/2L_box
USE kinds, ONLY: DP
!
IMPLICIT NONE
SAVE
INTEGER, ALLOCATABLE :: ind_localisation(:)
INTEGER :: nat_localisation = 0
LOGICAL :: print_localisation = .FALSE. ! Calculates hartree energy around specified atoms
INTEGER :: self_interaction = 0
REAL(DP) :: sic_epsilon = 0.0_DP
REAL(DP) :: sic_alpha = 0.0_DP
REAL(DP) :: sic_rloc = 0.0_DP
REAL(DP), ALLOCATABLE :: pos_localisation(:,:)
!------------------------------------------------------------------------------!
CONTAINS
!------------------------------------------------------------------------------!
SUBROUTINE sic_initval( nat_ , id_loc_ , sic_ , sic_epsilon_ , sic_alpha_, sic_rloc_ )
IMPLICIT NONE
INTEGER, INTENT(IN) :: nat_
INTEGER, INTENT(IN) :: id_loc_ (:)
CHARACTER(LEN=*), INTENT(IN) :: sic_
REAL(DP), INTENT(IN) :: sic_epsilon_
REAL(DP), INTENT(IN) :: sic_alpha_
REAL(DP), INTENT(IN) :: sic_rloc_
select case ( TRIM( sic_ ) )
case ( 'sic_mac' )
self_interaction = 2
case default
self_interaction = 0
end select
sic_epsilon = sic_epsilon_
sic_alpha = sic_alpha_
sic_rloc = sic_rloc_
! counting the atoms around which i want to calculate the charge localization
IF( ALLOCATED( ind_localisation ) ) DEALLOCATE( ind_localisation )
ALLOCATE( ind_localisation( nat_ ) )
ind_localisation( 1 : nat_ ) = id_loc_ ( 1 : nat_ )
nat_localisation = COUNT( ind_localisation > 0 )
IF( ALLOCATED( pos_localisation ) ) DEALLOCATE( pos_localisation )
ALLOCATE( pos_localisation( 4, MAX( nat_localisation, 1 ) ) )
!
IF( nat_localisation > 0 ) print_localisation = .TRUE.
!
RETURN
END SUBROUTINE sic_initval
!------------------------------------------------------------------------------!
SUBROUTINE deallocate_sic()
IMPLICIT NONE
IF( ALLOCATED( pos_localisation ) ) DEALLOCATE( pos_localisation )
IF( ALLOCATED( ind_localisation ) ) DEALLOCATE( ind_localisation )
RETURN
END SUBROUTINE deallocate_sic
!------------------------------------------------------------------------------!
SUBROUTINE sic_info( )
USE io_global, ONLY: stdout
IMPLICIT NONE
!
! prints the type of USIC we will do :
!
IF( self_interaction == 0 ) THEN
RETURN
END IF
WRITE(stdout, 591)
WRITE(stdout, 592) self_interaction
WRITE(stdout, 593)
!!select case (self_interaction)
IF ( self_interaction /= 0 ) THEN
write(stdout,*) &
' Unpaired-electron self-interaction correction by Mauri', self_interaction
write(stdout,*) &
' E_USIC_EHTE = U_hartree[rho_up + rho_dw]- sic_espilon * U_hartree[rho_up-rhp_down]'
write(stdout,*) &
' E_USIC_XC = E_xc[rho_up,rho_dw] - sic_alpha( E_xc[rho_up,rho_dw] + E_xc[rho_dw, rho_dw]) '
END IF !!select
591 FORMAT( 3X,' ')
592 FORMAT( 3X,'Introducing a Mauri Avezac Calandra Self_Interaction Correction: ', I3)
593 FORMAT( 3X,'----------------------------------------')
RETURN
END SUBROUTINE sic_info
!------------------------------------------------------------------------------!
END MODULE sic_module
!------------------------------------------------------------------------------!
|
