1!-----------------------------------------------------------------------------------------------------------------------
2      subroutine get_sc_ensemble_energy(VEENSEMB,EDFT_orig,             &
3     &                                  cmo,nnashx,nnorbt,              &
4     &                                  ncroot,weights,lupri,           &
5     &                                  wrk,lwrk)
6
7      !> brief: calculate the ensemble-DFT energy in a self-consistent manner (see notes by E. Fromager)
8      !> author: Stefan Knecht and Manu Fromager
9!-----------------------------------------------------------------------------------------------------------------------
10      implicit none
11
12      integer, intent(in)    :: ncroot,lupri,lwrk,nnashx,nnorbt
13      real*8 , intent(in)    :: EDFT_orig
14      real*8 , intent(in)    :: weights(ncroot)
15      real*8 , intent(in)    :: cmo(*)
16      real*8 , intent(inout) :: wrk(lwrk)
17      real*8 , intent(inout) :: VEENSEMB
18!-----------------------------------------------------------------------------------------------------------------------
19      integer                :: kw2, KW2A, KUDVREF, KDVREF
20      integer                :: LW2A, ist
21      real*8                 :: uESRDV, uEDSR, uEJVSR, uEJCSR
22      real*8                 :: uEMYDFTAUX, uedft(3)
23      real*8                 :: UEJCVSR
24      real*8                 :: diff_ejvsr,diff_edsr
25
26      real*8 , allocatable   :: u_rho_ensemble(:)
27      real*8 , allocatable   :: ref_rho_ensemble(:)
28      real*8 , allocatable   :: diff_rho_ensemble(:)
29      real*8 , allocatable   :: ufsr(:)
30      real*8 , external      :: ddot
31
32      kw2      = 1
33
34      !> build new ensemble density matrix
35      allocate(u_rho_ensemble(nnashx),ref_rho_ensemble(nnashx),         &
36     &         diff_rho_ensemble(nnashx))
37      u_rho_ensemble    = 0
38      ref_rho_ensemble  = 0
39      diff_rho_ensemble = 0
40
41      do IST = 1,NCROOT
42        KUDVREF = KW2
43        KDVREF  = KUDVREF + NNASHX
44        KW2A    = KDVREF  + NNASHX
45        LW2A    = LWRK   - KW2A
46        CALL GETUDVREF('LUCITA   ',WRK(KUDVREF),IST,wrk,WRK(KW2A),LW2A)
47        CALL GETDVREF('LUCITA   ',WRK(KDVREF),IST,wrk,WRK(KW2A),LW2A)
48        call daxpy(nnashx,weights(ist),WRK(KUDVREF),1,u_rho_ensemble,1)
49        call daxpy(nnashx,weights(ist),WRK(KDVREF),1,ref_rho_ensemble,1)
50      end do
51      call dcopy(nnashx,u_rho_ensemble,1,diff_rho_ensemble,1)
52      call daxpy(nnashx,-1.0d0,ref_rho_ensemble,1,diff_rho_ensemble,1)
53
54      kw2 = 1
55
56      allocate(ufsr(nnorbt))
57      ufsr  = 0
58
59      CALL SRFMAT_vensemble(ufsr,cmo,u_rho_ensemble,ref_rho_ensemble,   &
60     &                      diff_rho_ensemble,                          &
61     &                      diff_ejvsr,diff_edsr,uEDFT,                 &
62     &                      wrk(kw2),lwrk,0)
63
64      VEENSEMB = VEENSEMB + uEDFT(1) - EDFT_orig                        &
65     &                    + diff_edsr - diff_ejvsr
66
67#ifdef DEBUG
68      write(lupri,*) 'contributions to the variational ensemble energy'
69      write(lupri,*) '------------------------------------------------'
70      write(lupri,*) 'uEDFT        --> ',uEDFT
71      write(lupri,*) ' EDFT_orig   --> ',EDFT_orig
72      write(lupri,*) ' diff_edsr   --> ',diff_edsr
73      write(lupri,*) ' diff_ejvsr  --> ',diff_ejvsr
74      write(lupri,*) '------------------------------------------------'
75#endif
76
77      write(lupri,806) '******************************************'
78      write(lupri,806) 'Variational ensemble CI-srDFT energy:    ',     &
79     &                  VEENSEMB
80      write(lupri,806) '*******************************************'
81!
82  805   FORMAT( 1X,A,F25.12)
83  806   FORMAT(/1X,A,F25.12)
84  807   FORMAT(/1X,A,I2,A,F20.12)
85
86      deallocate(u_rho_ensemble,ref_rho_ensemble,diff_rho_ensemble,ufsr)
87
88      end subroutine get_sc_ensemble_energy
89!-----------------------------------------------------------------------------------------------------------------------
90!-----------------------------------------------------------------------------------------------------------------------
91#ifdef backup
92      subroutine get_sc_ensemble_energy(emydft_orig,EJCSR_orig,         &
93     &                                  EJVSR_orig,EDSR_orig,           &
94     &                                  EDFT_orig,EMYDFTAUX_orig,       &
95     &                                  eci_orig,emy_ci_orig,           &
96     &                                  energy_root_orig,eensemb_orig,  &
97     &                                  cmo,srac_orig,nnashx,nnorbt,    &
98     &                                  ncroot,istaci,weights,lupri,    &
99     &                                  potnuc,                         &
100     &                                  wrk,lwrk)
101
102      !> brief: calculate the ensemble-DFT energy in a self-consistent manner (see notes by E. Fromager)
103      !> author: Stefan Knecht and Manu Fromager
104!-----------------------------------------------------------------------------------------------------------------------
105      implicit none
106
107      integer, intent(in)    :: ncroot,istaci,lupri,lwrk,nnashx,nnorbt
108      real*8 , intent(in)    :: emydft_orig, EJCSR_orig, EJVSR_orig
109      real*8 , intent(in)    :: EDFT_orig, EMYDFTAUX_orig, emy_ci_orig
110      real*8 , intent(in)    :: eensemb_orig, EDSR_orig
111      real*8 , intent(in)    :: eci_orig(ncroot)
112      real*8 , intent(in)    :: energy_root_orig(ncroot)
113      real*8 , intent(in)    :: potnuc
114      real*8 , intent(in)    :: weights(ncroot), srac_orig(nnashx)
115      real*8 , intent(in)    :: cmo(*)
116      real*8 , intent(inout) :: wrk(lwrk)
117!-----------------------------------------------------------------------------------------------------------------------
118      integer                :: kw2, KW2A, KUDVREF
119      integer                :: LW2A, ist
120      real*8                 :: uESRDV, uEDSR, uEJVSR, uEJCSR
121      real*8                 :: uEMYDFTAUX, uedft(3)
122      real*8                 :: UEJCVSR
123      real*8                 :: VEENSEMB
124
125      real*8 , allocatable   :: u_rho_ensemble(:)
126      real*8 , allocatable   :: ref_rho_ensemble(:)
127      real*8 , allocatable   :: diff_rho_ensemble(:)
128      real*8 , allocatable   :: ufsr(:)
129      real*8 , external      :: ddot
130
131      kw2 = 1
132
133      !> build new ensemble density matrix
134      allocate(u_rho_ensemble(nnashx),ref_rho_ensemble(nnashx),           &
135     ^         diff_rho_ensemble(nnashx))
136      u_rho_ensemble    = 0
137      ref_rho_ensemble  = 0
138      diff_rho_ensemble = 0
139
140      do IST = 1,NCROOT
141        KUDVREF = KW2
142        KDVREF  = KUDVREF + NNASHX
143        KW2A    = KDVREF  + NNASHX
144        LW2A    = LWRK   - KW2A
145        CALL GETUDVREF('LUCITA   ',WRK(KUDVREF),IST,wrk,WRK(KW2A),LW2A)
146        CALL GETDVREF('LUCITA   ',WRK(KDVREF),IST,wrk,WRK(KW2A),LW2A)
147        call daxpy(nnashx,weights(ist),WRK(KUDVREF),1,u_rho_ensemble,1)
148        call daxpy(nnashx,weights(ist),WRK(KDVREF),1,ref_rho_ensemble,1)
149        call dcopy(nnashx,u_rho_ensemble,1,diff_rho_ensemble,1)
150        call daxpy(nnashx,-1.0d0,ref_rho_ensemble,1,diff_rho_ensemble,1)
151      end do
152      kw2 = 1
153
154      allocate(ufsr(nnorbt))
155      ufsr  = 0
156
157      CALL SRFMAT(ufsr,cmo,u_rho_ensemble,                              &
158     &                      uEJCSR,uEJVSR,uEDSR,uEDFT,                  &
159     &                      uEMYDFTAUX,uEJCVSR,                         &
160     &                      wrk(kw2),lwrk,0)
161
162      uESRDV   = DDOT(NNASHX,u_rho_ensemble,1,srac_orig,1)
163
164      VEENSEMB = eensemb_orig + uEDFT(1)  - EDFT_orig - EDSR_orig       &
165     &         + uEJCVSR - uESRDV       - uEJVSR - EJVSR_orig
166
167      write(lupri,*) 'contributions to the variational ensemble energy'
168      write(lupri,*) '------------------------------------------------'
169      write(lupri,*) 'eensemb_orig --> ',eensemb_orig
170      write(lupri,*) 'uEDFT        --> ',uEDFT
171      write(lupri,*) ' EDFT_orig   --> ',EDFT_orig
172      write(lupri,*) ' EDSR_orig   --> ',EDSR_orig
173      write(lupri,*) 'uEJCVSR      --> ',uEJCVSR
174      write(lupri,*) 'uESRDV       --> ',uESRDV
175      write(lupri,*) 'uEJVSR       --> ',uEJVSR
176      write(lupri,*) ' EJVSR_orig  --> ',EJVSR_orig
177      write(lupri,*) '------------------------------------------------'
178
179      write(lupri,806) '*******************************'
180      write(lupri,806) 'Final variational ensemble CI-srDFT energy:',   &
181     &                  VEENSEMB
182      write(lupri,806) '*******************************'
183!
184  805   FORMAT( 1X,A,F25.12)
185  806   FORMAT(/1X,A,F25.12)
186  807   FORMAT(/1X,A,I2,A,F20.12)
187
188      deallocate(u_rho_ensemble,ref_rho_ensemble,diff_rho_ensemble,ufsr)
189
190      end subroutine get_sc_ensemble_energy
191!-----------------------------------------------------------------------------------------------------------------------
192#endif
193