Line data Source code
1 : !!****m* ABINIT/m_tddft
2 : !! NAME
3 : !! m_tddft
4 : !!
5 : !! FUNCTION
6 : !! Routines for computing excitation energies within TDDFT
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1999-2026 ABINIT group (XG, JYR, MB, MBELAND, SHAMEL)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
14 : !!
15 : !! NOTES
16 : !!
17 : !! SOURCE
18 :
19 : #if defined HAVE_CONFIG_H
20 : #include "config.h"
21 : #endif
22 :
23 : #include "abi_common.h"
24 :
25 : module m_tddft
26 :
27 : use defs_basis
28 : use m_abicore
29 : use m_xmpi
30 : use m_errors
31 : use m_wffile
32 : use m_sort
33 : use m_dtset
34 : use m_dtfil
35 : use iso_c_binding, only : c_ptr,c_loc,c_f_pointer
36 :
37 : use defs_abitypes, only : MPI_type
38 : use m_io_tools, only : get_unit
39 : use m_matrix, only : matr3inv
40 : use m_time, only : timab
41 : use m_fftcore, only : sphereboundary
42 : use m_spacepar, only : hartre
43 : use m_mpinfo, only : proc_distrb_cycle
44 : use m_fft, only : fourwf, fourdp
45 :
46 : implicit none
47 :
48 : private
49 : !!***
50 :
51 : #if defined HAVE_MPI1
52 : include 'mpif.h'
53 : #endif
54 :
55 : public :: tddft
56 : !!***
57 :
58 : contains
59 :
60 : !!****f* m_tddft/tddft
61 : !! NAME
62 : !! tddft
63 : !!
64 : !! FUNCTION
65 : !! Compute the excitation energies within TDLDA
66 : !! from input wavefunctions, eigenenergies, and band occupations.
67 : !!
68 : !! INPUTS
69 : !! cg(2,mpw*nspinor*mband*mkmem*nsppol)=wf in G space
70 : !! dtfil <type(datafiles_type)>=variables related to files
71 : !! dtset <type(dataset_type)>=all input variables for this dataset
72 : !! eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
73 : !! etotal=total energy of the ground-state (Ha)
74 : !! gmet(3,3)=metrix tensor in G space in Bohr**-2.
75 : !! gprimd(3,3)=dimensional reciprocal space primitive translations
76 : !! gsqcut=cutoff on (k+G)^2 (bohr^-2)
77 : !! kg(3,mpw*mkmem)=reduced planewave coordinates.
78 : !! kxc(nfft,nkxc)=exchange-correlation kernel
79 : !! mband=maximum number of bands
80 : !! mgfftdiel=maximum size of 1D FFTs, for the computation of the dielectric matrix
81 : !! mkmem=number of k-points treated by this node.
82 : !! mpi_enreg=information about MPI parallelization
83 : !! mpw=maximum allowed value for npw
84 : !! nfft=(effective) number of FFT grid points (for this processor)
85 : !! WARNING about parallelization: see below
86 : !! ngfftdiel(18)=contain all needed information about 3D FFT, for dielectric matrix,
87 : !! see ~abinit/doc/variables/vargs.htm#ngfft
88 : !! nkpt=number of k points
89 : !! nkxc=second dimension of the array kxc (see rhotoxc for a description)
90 : !! npwarr(nkpt)=number of planewaves at each k point
91 : !! nspinor=number of spinorial components of the wavefunctions
92 : !! nsppol=1 for unpolarized, 2 for spin-polarized
93 : !! occ(mband*nkpt*nsppol)=
94 : !! occupation numbers for each band (usually 2.0) at each k point
95 : !! ucvol=unit cell volume (Bohr**3)
96 : !! wffnew=unit number for current wf disk file
97 : !!
98 : !! OUTPUT
99 : !! (only writing)
100 : !!
101 : !! WARNING:
102 : !! This routine should not be parallelized on space for the time being,
103 : !! because the already existing parallelisation is not the usual one, found
104 : !! in the majority of ABINIT routines.
105 : !!
106 : !! NOTES
107 : !! * Only accept nspinor=1, nsppol=1, nkpt=1 (Gamma point), and occopt<3
108 : !! (insulating occupation numbers).
109 : !! It is expected to make it work for nsppol=2 in the future.
110 : !!
111 : !! * For the oscillator strengths, see the paper
112 : !! ''Time-Dependent Density Functional Response Theory of Molecular
113 : !! systems: Theory, Computational Methods, and Functionals'', by M.E. Casida,
114 : !! in Recent Developments and Applications of Modern Density Functional
115 : !! Theory, edited by J.M. Seminario (Elsevier, Amsterdam, 1996).
116 : !!
117 : !! SOURCE
118 :
119 22 : subroutine tddft(cg,dtfil,dtset,eigen,etotal,gmet,gprimd,gsqcut,&
120 22 : & kg,kxc,mband,mgfftdiel,mkmem,mpi_enreg,mpw,nfft,ngfftdiel,nkpt,nkxc,&
121 22 : & npwarr,nspinor,nsppol,occ,ucvol,wffnew)
122 :
123 : !Arguments ------------------------------------
124 : integer, intent(in) :: mband,mgfftdiel,mkmem,mpw,nfft,nkpt,nkxc,nsppol
125 : integer, intent(in) :: nspinor
126 : real(dp), intent(in) :: etotal,gsqcut,ucvol
127 : type(datafiles_type), intent(in) :: dtfil
128 : type(dataset_type), intent(in) :: dtset
129 : type(MPI_type), intent(in) :: mpi_enreg
130 : type(wffile_type), intent(inout) :: wffnew
131 : integer, intent(in) :: kg(3,mpw*mkmem),ngfftdiel(18),npwarr(nkpt)
132 : real(dp), intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol),eigen(mband*nkpt*nsppol)
133 : real(dp), intent(in) :: gmet(3,3),gprimd(3,3),kxc(nfft,nkxc),occ(mband*nkpt*nsppol)
134 :
135 : !Local variables-------------------------------
136 : integer,parameter :: nexcitout=20
137 : integer :: cplex,i1,i2,i3,iband,idir,ier,ierr,iexcit,iexcit1,iexcit2,ifft
138 : integer :: old_iexcit,ii,jj,isppol,jsppol,isppol1,isppol2,isppol_l,isppol_n
139 : integer :: isppol_n1,isppol_n2,iocc_n1,iocc_n2,iunocc_n1,iunocc_n2,temp_unit2
140 : integer :: ikpt,index,iocc,iocc1,iocc2,iocc_l,iocc_n
141 : integer :: istwf_k,iunocc,iunocc1,iunocc2,iunocc_l,iunocc_n,istate,jexcit
142 : integer :: jexcit_cbase,master,mcg_disk,me_loc
143 44 : integer :: nband_k(nsppol), nband_occ(nsppol), nband_unocc(nsppol)
144 44 : integer :: nstate_k, nstate_occ, nstate_unocc, nexcit_pol(nsppol)
145 : integer :: nstate_win,ndiel,ndiel1,ndiel2,ndiel3,ndiel4
146 : integer :: ndiel5,ndiel6,nexcit,nexcit_max,nexcit_win,nfftdiel,nlargest,nnext
147 : integer :: nnext1,nnext2
148 : integer :: nproc_loc,npw_k,pole_approx,sing_trip,spaceComm,mtag,tim_fourwf
149 : integer :: tim_rwwf,save_iomode
150 : integer :: rec,recl,idummy,jdummy
151 : real(dp) :: buffer,buffer_inv,diffeig,eigunocc,emax_win
152 : real(dp) :: factor,ff,flargest,fnext,fr_invsquare,fr_power
153 : real(dp) :: fnext1,fnext2
154 : real(dp) :: normint,myproduct,saa,sab,sbb
155 : real(dp) :: sumx
156 44 : real(dp) :: sum_kernel(2/nsppol)
157 : real(dp) :: weight,xx
158 : logical :: am_master,file_exist
159 22 : logical, allocatable :: done_excit(:,:),done_sexc(:) !,done_sexc2(:)
160 : character(len=18) :: chain1,chain2
161 : character(len=500) :: message
162 22 : integer,allocatable :: flag_state_win(:),gbound(:,:),indarr(:),index_state(:)
163 22 : integer,allocatable :: kg_k(:,:)
164 22 : integer,allocatable :: excit_coords(:,:)
165 : integer :: count_to_do, count, displ, countmax, displmax
166 : integer :: ijexcit, ijexcit2, sendcount
167 66 : real(dp) :: f_sing_trip(2/nsppol),sendbuf(5-nsppol)
168 : real(dp) :: cauchy(7),poscart(3),rprimd(3,3),tsec(2),dummy(2,1)
169 : integer :: iomode,action,me,nmaster,sender,source,sread,sskip
170 : integer :: formeig,icg,ikg,nband_k_
171 : logical :: mydata, tmaster, swrite
172 22 : integer,allocatable :: kg_disk(:,:)
173 22 : integer,allocatable :: counts(:),displs(:),recvcounts(:),tmpbuf(:)
174 22 : real(dp),allocatable :: cg_disk(:,:),cg_tmp(:,:)
175 22 : real(dp),allocatable :: cwavef(:,:),eexcit(:)
176 22 : real(dp),allocatable :: eexcit2(:)
177 22 : real(dp),allocatable :: matr(:)
178 22 : real(dp),allocatable :: kxc_for_tddft(:,:,:,:,:,:),omega_tddft_casida(:,:,:,:,:,:,:)
179 22 : real(dp),allocatable :: osc_str(:,:),pos(:,:),rhoaug(:,:,:),rhog(:,:)
180 22 : real(dp),allocatable :: sexc(:,:),sqrtks(:),vec(:,:,:),vhartr(:),wfprod(:,:,:)
181 44 : real(dp) :: omega_tddft_casida_dummy(2/nsppol)
182 22 : real(dp),allocatable :: wfraug(:,:,:,:),wfrspa(:,:,:,:),work(:),zhpev1(:,:)
183 22 : real(dp),allocatable :: zhpev2(:)
184 :
185 : integer :: iproc
186 : integer :: ipwnbd
187 22 : real(dp), allocatable,target :: recvbuf(:,:)
188 22 : real(dp),pointer :: recvbuf_ptr(:)
189 : type(c_ptr) :: cptr
190 :
191 : ! *************************************************************************
192 :
193 : !Init mpi_comm
194 22 : spaceComm=mpi_enreg%comm_cell
195 :
196 22 : am_master=.true.
197 22 : master = 0
198 44 : nproc_loc = xmpi_comm_size(spaceComm) !Init ntot proc max
199 22 : me_loc = xmpi_comm_rank(spaceComm) !Define who i am
200 :
201 : #if defined HAVE_MPI
202 22 : if (me_loc/=0) then
203 4 : am_master=.FALSE.
204 : end if
205 22 : write(message, '(a,i3,a)' ) ' TDDFT ',nproc_loc,' CPU synchronized'
206 22 : call wrtout(std_out,message,'COLL')
207 22 : write(message, '(a,3D12.5,a,3D12.5,a,3D12.5)' ) ' gmet ',&
208 22 : & gmet(1,1),gmet(1,2),gmet(1,3),ch10,&
209 22 : & gmet(2,1),gmet(2,2),gmet(2,3),ch10,&
210 44 : & gmet(3,1),gmet(3,2),gmet(3,3)
211 22 : call wrtout(std_out,message,'COLL')
212 : #endif
213 :
214 :
215 : !COMMENT these values should become arguments
216 : !the two first define the energy window
217 :
218 22 : emax_win=greatest_real*tol6
219 22 : if(dtset%td_maxene>tol6)then
220 1 : emax_win = dtset%td_maxene
221 : end if
222 :
223 22 : call timab(95,1,tsec)
224 :
225 22 : istwf_k=dtset%istwfk(1)
226 :
227 22 : if(nkpt/=1 .or. &
228 : & abs(dtset%kptns(1,1))+abs(dtset%kptns(2,1))+abs(dtset%kptns(3,1))>1.0d-6 )then
229 : write(message, '(a,a,a,a,a,i4,a,3es14.6,a,a,a,a,a)' )&
230 0 : & 'The computation of excited states using TDDFT is only allowed',ch10,&
231 0 : & 'with nkpt=1, kpt=(0 0 0), but the following values are input:',ch10,&
232 0 : & 'nkpt=',nkpt,', kpt=',dtset%kptns(1:3,1),'.',ch10,&
233 0 : & 'Action: in the input file, set nkpt to 1 and kpt to 0 0 0 ,',ch10,&
234 0 : & 'or change iscf.'
235 0 : ABI_ERROR(message)
236 : end if
237 :
238 22 : if(nspinor/=1)then
239 : write(message, '(a,a,a,a,a,a,a)' )&
240 0 : & 'The computation of excited states using TDDFT is restricted',ch10,&
241 0 : & 'for the time being to nspinor=1, while input nspinor=2.',ch10,&
242 0 : & 'Action: if you want to compute excited states within TDDFT,',ch10,&
243 0 : & 'set nsppol to 1 in the input file. Otherwise, do not use iscf=-1.'
244 0 : ABI_ERROR(message)
245 : end if
246 :
247 :
248 22 : if(nsppol==2 .and. (dtset%ixc==22 .or. dtset%ixc==20))then
249 : write(message, '(a,a,a,a,a,a,a,a,a,a,a)' )&
250 0 : & 'The computation of excited states using TDDFT in the spin',ch10,&
251 0 : & 'polarized case for the time being cannot be used with ixc=20',ch10,&
252 0 : & 'or ixc=22',ch10,&
253 0 : & 'Action: if you want to compute excited states within TDDFT,',ch10,&
254 0 : & 'set ixc different from 20 or 22. Otherwise, do not use iscf=-1',ch10,&
255 0 : & 'with nsppol=2.'
256 0 : ABI_ERROR(message)
257 : end if
258 :
259 :
260 22 : if(dtset%occopt>2)then
261 : write(message, '(a,a,a,i2,a,a,a,a,a)' )&
262 0 : & 'The computation of excited states using TDDFT is only allowed',ch10,&
263 0 : & 'with occopt=0, 1, or 2, while input occopt=',dtset%occopt,'.',ch10,&
264 0 : & 'Action: if you want to compute excited states within TDDFT,',ch10,&
265 0 : & 'set occopt=0, 1, or 2 in the input file. Otherwise, do not use iscf=-1.'
266 0 : ABI_ERROR(message)
267 : end if
268 :
269 : !Examine the occupation numbers, and determine the number of
270 : !occupied and unoccupied states and band.
271 : !States are numerated as usual in Abinit, before all spin up band
272 : !and after all spin down bands.
273 : !Note that if nsppol==1 nstate=nband_k
274 54 : do isppol=1,nsppol
275 32 : nband_k(isppol)=dtset%nband(isppol)
276 32 : nband_occ(isppol)=0
277 470 : do iband=1,nband_k(isppol)
278 438 : if(abs(occ(iband+(isppol-1)*nband_k(1))-two/nsppol)<tol6) &
279 142 : & nband_occ(isppol)=nband_occ(isppol)+1
280 : end do
281 32 : nband_unocc(isppol)=nband_k(isppol)-nband_occ(isppol)
282 : ! next line make no sense if spin flip is taken into account
283 54 : nexcit_pol(isppol)=nband_occ(isppol)*nband_unocc(isppol)
284 : end do
285 22 : nstate_k=nband_k(1)+(nsppol-1)*nband_k(nsppol)
286 22 : nstate_occ=nband_occ(1)+(nsppol-1)*nband_occ(nsppol)
287 22 : nstate_unocc=nstate_k-nstate_occ
288 : !next line to be changed if spin fli is taken into account
289 22 : nexcit=nexcit_pol(1)+(nsppol-1)*nexcit_pol(nsppol)
290 :
291 : !number of plane wave (does it work even for nsppol=2 ??)
292 22 : npw_k=npwarr(1)
293 :
294 : !mux number of excitations that is taken into account
295 22 : if(dtset%td_mexcit==0)then
296 : nexcit_max=nexcit
297 : else
298 1 : nexcit_max =dtset%td_mexcit
299 : end if
300 :
301 : !DEBUG
302 : !write(std_out,*) nband_occ(1),nband_unocc(1)
303 : !write(std_out,*) nband_occ(nsppol),nband_unocc(nsppol)
304 : !END DEBUG
305 :
306 :
307 22 : if(nsppol==1)then
308 12 : write(message, '(a,a,a,a,i4,a,i4,a,a,i4,a,a,a,i6,a)' )ch10,&
309 12 : & ' *** TDDFT : computation of excited states *** ',ch10,&
310 12 : & ' Splitting of',dtset%nband(1),' states in',nband_occ(1),' occupied states,',&
311 12 : & ' and',nband_unocc(1),' unoccupied states,',ch10,&
312 24 : & ' giving',nexcit,' excitations.'
313 12 : call wrtout(std_out,message,'COLL')
314 12 : call wrtout(ab_out,message,'COLL')
315 : else
316 10 : write(message, '(a,a,a,a,i4,a,i4,a,a,i4,a,a,a,i6,a,a,a)' )ch10,&
317 10 : & ' *** TDDFT : computation of excited states *** ',ch10,&
318 10 : & ' Splitting of',nstate_k,' states in',nstate_occ,' occupied states,',&
319 10 : & ' and',nstate_unocc,' unoccupied states,',ch10,&
320 10 : & ' giving',nexcit,' excitations. Note that spin flip is not possible actually.',ch10,&
321 20 : & ' So the number of excitation is the half of the product of the number of state'
322 10 : call wrtout(std_out,message,'COLL')
323 10 : call wrtout(ab_out,message,'COLL')
324 : end if
325 :
326 : !Allocate the matrices to be diagonalized.
327 : !Use a simple storage mode, to be improved in the future.
328 22 : ii=max(nband_occ(1),nband_occ(nsppol))
329 22 : jj=max(nband_unocc(1),nband_unocc(nsppol))
330 198 : ABI_MALLOC(omega_tddft_casida,(ii,jj,nsppol,ii,jj,nsppol,2/nsppol))
331 66 : ABI_MALLOC(eexcit,(nexcit))
332 44 : ABI_MALLOC(sqrtks,(nexcit))
333 66 : ABI_MALLOC(flag_state_win,(nstate_k))
334 84196 : omega_tddft_casida(:,:,:,:,:,:,:)=zero
335 :
336 :
337 : !Fill the diagonal elements with square of differences of KS eigenvalues
338 : !(also not very efficient, but OK for the present first coding)
339 : !Also compute the square root of Kohn-Sham eigenvalue differences
340 54 : do isppol=1,nsppol
341 382 : do iunocc=1,nband_unocc(isppol)
342 328 : eigunocc=eigen(iunocc+nband_occ(isppol)+(isppol-1)*nband_k(1))
343 1256 : do iocc=1,nband_occ(isppol)
344 896 : iexcit=iocc+(isppol-1)*nexcit_pol(1)+nband_occ(isppol)*(iunocc-1)
345 896 : diffeig=eigunocc-eigen(iocc+(isppol-1)*nband_k(1))
346 2018 : do sing_trip=1,2/nsppol
347 2018 : omega_tddft_casida(iocc,iunocc,isppol,iocc,iunocc,isppol,sing_trip)=diffeig**2
348 : end do
349 896 : eexcit(iexcit)=diffeig
350 1224 : sqrtks(iexcit)=sqrt(diffeig)
351 : end do
352 : end do
353 : end do
354 :
355 :
356 :
357 : !Sort the excitation energies : note that the array eexcit is reordered
358 66 : ABI_MALLOC(indarr,(nexcit))
359 1814 : indarr(:)=(/ (ii,ii=1,nexcit) /)
360 22 : call sort_dp(nexcit,eexcit,indarr,tol14)
361 :
362 : !Determine an energy window for the excitations
363 : !to take into account. This is necessary for large systems
364 :
365 22 : nexcit_win = 0
366 918 : do iexcit = 1, nexcit
367 918 : if ((eexcit(iexcit) < emax_win ).and.(nexcit_win < nexcit_max)) then
368 890 : nexcit_win = nexcit_win + 1
369 :
370 : ! DEBUG
371 : ! write(message,'(a,F12.5,a,a,i2,a,a,i2)') 'excitation energy:', eexcit(indarr(iexcit)),ch10, &
372 : ! & 'excitation number:', indarr(iexcit),ch10, &
373 : ! & 'nexcit_win: ', nexcit_win
374 : ! call wrtout(std_out,message,'COLL')
375 : ! ENDDEBUG
376 :
377 : end if
378 : end do
379 :
380 : !identification of the bands contributing to the
381 : !nexcit_win excitations within the window
382 :
383 :
384 22 : nstate_win = 0
385 460 : flag_state_win(:) = 0
386 912 : do iexcit = 1, nexcit_win
387 890 : iexcit1 = indarr(iexcit)
388 890 : isppol1 = min((iexcit1-1)/nexcit_pol(1) +1,2)
389 890 : iunocc1 = (iexcit1-(isppol1-1)*nexcit_pol(1)-1)/nband_occ(isppol1)+1
390 890 : iocc1 = iexcit1-(isppol1-1)*nexcit_pol(1)-(iunocc1-1)*nband_occ(isppol1)
391 890 : if (flag_state_win(nband_occ(isppol1)+(isppol1-1)*nband_k(1)+iunocc1)==0) &
392 322 : & flag_state_win(nband_occ(isppol1)+(isppol1-1)*nband_k(1)+iunocc1) =1
393 890 : if (flag_state_win(iocc1+(isppol1-1)*nband_k(1))==0) &
394 132 : & flag_state_win(iocc1+(isppol1-1)*nband_k(1)) =1
395 : ! DEBUG
396 : ! write(message,'(a,i2,a,a,i2,i2,a,a,i2,a,a,i2)') 'isppol:', isppol1,ch10, &
397 : ! & 'iocc,iunocc:', iocc1,iunocc1,ch10, &
398 : ! & 'flag_state_win:', flag_state_win(iocc1+(isppol1-1)*nband_k(1)), &
399 : ! & ch10, 'flag_state_win:', flag_state_win(nband_occ(isppol1)+(isppol1-1)*nband_k(1)+iunocc1)
400 : ! call wrtout(std_out,message,'COLL')
401 : ! END DEBUG
402 :
403 : end do
404 :
405 54 : do isppol=1,nsppol
406 492 : do iband=1,nband_k(isppol)
407 470 : nstate_win=nstate_win+flag_state_win(iband+(isppol-1)*nband_k(1))
408 : end do
409 : end do
410 :
411 :
412 22 : write(message,'(a,a,i5)') ch10,'Nr of states to Fourier transform : ',nstate_win
413 22 : call wrtout(std_out,message,'COLL')
414 :
415 22 : ndiel1=ngfftdiel(1) ; ndiel2=ngfftdiel(2) ; ndiel3=ngfftdiel(3)
416 : !ndiel4,ndiel5,ndiel6 are FFT dimensions, modified to avoid cache trashing
417 22 : ndiel4=ngfftdiel(4) ; ndiel5=ngfftdiel(5) ; ndiel6=ngfftdiel(6)
418 :
419 : !The evaluation of integrals, later, needs the following factor
420 22 : normint=one/(ucvol*dble(ndiel1*ndiel2*ndiel3))
421 :
422 : !Setup the positions in real space for later integration
423 22 : call matr3inv(gprimd,rprimd)
424 :
425 66 : ABI_MALLOC(pos,(max(ndiel1,ndiel2,ndiel3),3))
426 :
427 : !Select the reduced position of the point with respect to the box center,
428 : !in the interval ]-0.5,0.5].
429 : buffer=0.05_dp ; buffer_inv=one/buffer
430 88 : do idir=1,3
431 66 : if(idir==1)ndiel=ndiel1
432 44 : if(idir==2)ndiel=ndiel2
433 44 : if(idir==3)ndiel=ndiel3
434 2734 : do ii=1,ndiel
435 : ! dtset%boxcenter(3)=reduced coordinates of the center of the box,
436 : ! in view of the computation of the oscillator strength
437 2646 : pos(ii,idir)=(ii-1)/(one*ndiel)-dtset%boxcenter(idir)
438 2646 : pos(ii,idir)=pos(ii,idir)-nint(pos(ii,idir)-tol12)
439 : ! The linear behaviour is cut-off when one becomes
440 : ! close to the boundaries : the buffer allows to match smoothly
441 : ! one side of the cell to the other. This is important
442 : ! to get rid of small breakings of symmetry, that are
443 : ! confusing in accurate tests
444 2712 : if(abs(pos(ii,idir))>half-buffer)then
445 : ! xx is always positive, and goes linearly from 1 to 0
446 : ! in the buffer region
447 242 : xx=(half-abs(pos(ii,idir)))*buffer_inv
448 : ! The cut-off is applied to pos(:,:)
449 242 : pos(ii,idir)=pos(ii,idir)*xx*(two-xx)
450 : ! DEBUG
451 : ! if (idir==1)then
452 : ! write(std_out,'(i2)') ndiel
453 : ! write(std_out,'(a,i2,a,F12.5,F12.5)')'idiel : ',ii,' x : ',pos(ii,idir),&
454 : ! & dtset%boxcenter(idir)
455 : ! endif
456 : ! ENDDEBUG
457 : end if
458 : end do ! ii
459 : end do ! idir
460 :
461 : !need to run in MPI I/O case
462 22 : if (wffnew%iomode == IO_MODE_MPI ) then
463 0 : save_iomode=wffnew%iomode
464 0 : wffnew%iomode = IO_MODE_FORTRAN
465 : else
466 : ! Do not store value but set to have save_iomode /= 1
467 : save_iomode = IO_MODE_FORTRAN
468 : end if
469 :
470 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
471 : !2009 Chunping Hu
472 : !need to collect wavefunctions from processors to master
473 22 : me=me_loc
474 :
475 22 : tim_rwwf =0
476 22 : source = master
477 22 : sread = master
478 22 : tmaster=(master==me)
479 22 : swrite=tmaster
480 22 : sender=-1
481 :
482 22 : iomode=wffnew%iomode
483 :
484 22 : if(am_master)then
485 : #if defined HAVE_MPI
486 54 : ABI_MALLOC(cg_tmp,(2,mpw*nspinor*mband*nsppol))
487 : #endif
488 : end if
489 :
490 66 : ABI_MALLOC(kg_disk,(3,mpw))
491 22 : mcg_disk=mpw*nspinor*mband
492 22 : formeig=0
493 :
494 : #if defined HAVE_MPI
495 22 : call xmpi_barrier(spaceComm)
496 66 : ABI_MALLOC(cg_disk,(2,mcg_disk))
497 : #endif
498 :
499 22 : icg=0
500 : if(mpi_enreg%paralbd==0) tim_rwwf=6
501 22 : if(mpi_enreg%paralbd==1)tim_rwwf=12
502 :
503 54 : do isppol=1,nsppol
504 32 : ikg=0
505 86 : do ikpt=1,nkpt
506 32 : nband_k_=dtset%nband(ikpt+(isppol-1)*nkpt)
507 32 : npw_k=npwarr(ikpt)
508 : #if defined HAVE_MPI
509 32 : if (dtset%usewvl == 0) then
510 32 : mtag=ikpt+(isppol-1)*nkpt
511 32 : call xmpi_barrier(spaceComm)
512 : ! Must transfer the wavefunctions to the master processor
513 : ! Separate sections for paralbd=1 or other values ; might be merged
514 32 : if(mpi_enreg%paralbd==0) then
515 32 : nmaster=0
516 470 : source=minval(mpi_enreg%proc_distrb(ikpt,1:nband_k_,isppol))
517 32 : mydata=.false.
518 32 : if(source==me)mydata=.true.
519 24 : action=0
520 : ! I am the master node, and I have the data in cg or cg_disk
521 24 : if((tmaster).and.(mydata))action=1
522 : ! I am not the master, and I have the data => send to master
523 32 : if((.not.tmaster).and.(mydata))action=2
524 : ! I am the master, and I receive the data
525 32 : if((tmaster).and.(.not.mydata))action=3
526 : ! I have the data in cg or cg_disk ( MPI_IO case)
527 32 : if (iomode==IO_MODE_MPI) then
528 0 : action = 0
529 0 : sender=-1
530 0 : swrite=.false.
531 0 : if (mydata)then
532 : action=1
533 22 : swrite=.true.
534 22 : sender=me
535 : end if
536 : end if
537 : ! I am the master node, and I have the data in cg or cg_disk
538 : ! I have the data in cg or cg_disk ( MPI_IO case)
539 32 : if(action==1)then
540 : ! Copy from kg to kg_disk
541 162074 : kg_disk(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
542 : ! Copy from cg to cg_disk
543 546050 : do ipwnbd=1,nband_k_*npw_k*nspinor
544 546028 : cg_disk(1,ipwnbd)=cg(1,ipwnbd+icg)
545 546050 : cg_disk(2,ipwnbd)=cg(2,ipwnbd+icg)
546 : end do
547 : end if
548 : ! I am not the master, and I have the data => send to master
549 : ! I am the master, and I receive the data
550 32 : if ( action==2.or.action==3) then
551 4 : call timab(48,1,tsec)
552 4 : if(action==2)then
553 2 : call xmpi_exch(kg(:,1+ikg:npw_k+ikg),3*npw_k,source,kg_disk,nmaster,spaceComm,2*(mtag-1)+1,ierr)
554 : call xmpi_exch(cg(:,icg+1:icg+nband_k_*npw_k*nspinor),2*nband_k_*npw_k*nspinor &
555 2 : & ,source,cg_disk,nmaster,spaceComm,2*(mtag-1)+2,ierr)
556 : else
557 2 : call xmpi_exch(kg_disk,3*npw_k,source,kg_disk,nmaster,spaceComm,2*(mtag-1)+1,ierr)
558 2 : call xmpi_exch(cg_disk,2*nband_k_*npw_k*nspinor,source,cg_disk,nmaster,spaceComm,2*(mtag-1)+2,ierr)
559 : end if
560 4 : call timab(48,2,tsec)
561 : end if
562 0 : else if(mpi_enreg%paralbd==1)then
563 0 : nmaster=0
564 : #if defined HAVE_MPI_IO
565 0 : sender=-1
566 0 : if( iomode ==IO_MODE_MPI ) then
567 0 : nmaster=mpi_enreg%proc_distrb(ikpt,1,isppol)
568 0 : sender=nmaster
569 : end if
570 : #endif
571 : ! Note the loop over bands
572 0 : do iband=1,nband_k_
573 : ! The message passing related to kg is counted as one band
574 0 : action=0
575 : ! I am the master node, and I have the data in cg or cg_disk
576 0 : if( mpi_enreg%proc_distrb(ikpt,iband,isppol)==nmaster .and. me==nmaster) then
577 : action=1
578 : ! I am not the master, and I have the data => send to master
579 0 : elseif( mpi_enreg%proc_distrb(ikpt,iband,isppol)==me .and. me/=nmaster ) then
580 : action = 2
581 : ! I am the master, and I receive the data
582 0 : elseif( mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me .and. me==nmaster ) then
583 : action=3
584 : end if
585 : if(action==1) then
586 : ! I am the master node, and I have the data in cg or cg_disk
587 : ! Copy from kg to kg_disk
588 0 : if(iband==1)kg_disk(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
589 : ! Copy from cg to cg_disk
590 0 : do ipwnbd=1,npw_k*nspinor
591 0 : cg_disk(1,(iband-1)*npw_k*nspinor+ipwnbd)= cg(1,(iband-1)*npw_k*nspinor+ipwnbd+icg)
592 0 : cg_disk(2,(iband-1)*npw_k*nspinor+ipwnbd)= cg(2,(iband-1)*npw_k*nspinor+ipwnbd+icg)
593 : end do
594 : end if ! action=1
595 0 : if ( action==2.or.action==3) then
596 : ! action=2 : I am not the master, and I have the data => send to master
597 : ! action=3 : I am the master, and I receive the data
598 0 : call timab(48,1,tsec)
599 0 : if ( iband == 1 ) then
600 0 : if (action==2) then
601 : call xmpi_exch(kg(:,1+ikg:npw_k+ikg),3*npw_k,mpi_enreg%proc_distrb(ikpt,iband,isppol) &
602 0 : & ,kg_disk,nmaster,spaceComm,iband*(mtag-1)+1,ierr)
603 : else
604 : call xmpi_exch(kg_disk,3*npw_k,mpi_enreg%proc_distrb(ikpt,iband,isppol) &
605 0 : & ,kg_disk,nmaster,spaceComm,iband*(mtag-1)+1,ierr)
606 : end if
607 : end if ! iband =1
608 0 : ipwnbd=(iband-1)*npw_k*nspinor
609 0 : if (action==2)then
610 : call xmpi_exch( cg(:,ipwnbd+icg+1:ipwnbd+icg+npw_k*nspinor),2*npw_k*nspinor &
611 : & ,mpi_enreg%proc_distrb(ikpt,iband,isppol) &
612 0 : & ,cg_disk(:,ipwnbd+1:ipwnbd+npw_k*nspinor),nmaster,spaceComm,iband*(mtag-1)+2,ierr)
613 : else
614 : call xmpi_exch( cg_disk(:,ipwnbd+1:ipwnbd+npw_k*nspinor),2*npw_k*nspinor &
615 : & ,mpi_enreg%proc_distrb(ikpt,iband,isppol) &
616 0 : & ,cg_disk(:,ipwnbd+1:ipwnbd+npw_k*nspinor),nmaster,spaceComm,iband*(mtag-1)+2,ierr)
617 : end if
618 0 : call timab(48,2,tsec)
619 : end if ! action=2 or action=3
620 0 : if(iomode ==IO_MODE_MPI) then
621 : ! I have the data in cg or cg_disk
622 : swrite=.false.
623 : if (nmaster == me) then
624 0 : swrite=.true.
625 : end if
626 : end if
627 : ! End of loop over bands
628 : end do
629 : ! End of paralbd=1
630 : end if
631 : end if
632 : #endif
633 :
634 : ! The wavefunctions for the present k point and spin are stored into cg_tmp
635 32 : if(am_master)then
636 : #if defined HAVE_MPI
637 1856340 : cg_tmp(:,icg+1:icg+nband_k_*npw_k*nspinor)=cg_disk(:,:)
638 : #endif
639 : end if
640 :
641 32 : sskip=1
642 : #if defined HAVE_MPI
643 32 : if (dtset%usewvl == 0) then
644 32 : sskip=0
645 32 : if(.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k_,isppol,me)))sskip=1
646 : end if
647 : #endif
648 32 : if(sskip==1)then
649 24 : icg=icg+npw_k*nspinor*nband_k_
650 24 : ikg=ikg+npw_k
651 : end if
652 :
653 : end do ! ikpt
654 : end do ! isppol
655 22 : ABI_FREE(kg_disk)
656 : #if defined HAVE_MPI
657 22 : ABI_FREE(cg_disk)
658 : #endif
659 : !!!!!!!end of collecting wavefunction to master!!!!!!
660 :
661 :
662 22 : if(am_master)then
663 : ! -----------------------------------------------------------
664 : ! The disk access is only done by master...
665 :
666 72 : ABI_MALLOC(gbound,(2*mgfftdiel+8,2))
667 54 : ABI_MALLOC(kg_k,(3,npw_k))
668 :
669 18 : ikpt=1
670 : ! Only one k point
671 : ! do isppol=1,nsppol
672 126038 : kg_k(:,1:npw_k)=kg(:,1:npw_k)
673 18 : call sphereboundary(gbound,istwf_k,kg_k,mgfftdiel,npw_k)
674 : ! enddo
675 :
676 : end if ! am_master
677 :
678 : !need to run in MPI I/O case
679 22 : if ( save_iomode == 1 ) wffnew%iomode = IO_MODE_MPI
680 : !call wrtout(std_out,'After reading the wavefunction','COLL')
681 :
682 : !Use a simple implementation for the computation of the kernel elements
683 22 : if (am_master) then
684 54 : ABI_MALLOC(cwavef,(2,mpw))
685 90 : ABI_MALLOC(rhoaug,(ndiel4,ndiel5,ndiel6))
686 90 : ABI_MALLOC(wfraug,(2,ndiel4,ndiel5,ndiel6))
687 : end if
688 44 : ABI_MALLOC(index_state,(nstate_k))
689 :
690 : ! all real-space states are kept in memory
691 132 : ABI_MALLOC(wfrspa,(ndiel4,ndiel5,ndiel6,nstate_win))
692 :
693 : !DEBUG
694 : !write(message,'(a)') 'After allocating wfrspa'
695 : !call wrtout(std_out,message,'COLL')
696 : !ENDDEBUG
697 :
698 22 : weight=zero
699 :
700 : !Generate states in real space, only for states contributing to excitations in window
701 22 : istate=0
702 :
703 54 : do isppol=1,nsppol
704 492 : do iband=1,nband_k(isppol)
705 :
706 470 : if(flag_state_win(iband+(isppol-1)*nband_k(1)) == 1) then
707 432 : istate=istate+1
708 432 : index_state(iband+(isppol-1)*nband_k(1))=istate
709 :
710 432 : if (am_master) then
711 : #if defined HAVE_MPI
712 : ! Obtain Fourier transform in fft box
713 : cwavef(:,1:npw_k)=cg_tmp(:,1+(iband-1)*npw_k+(isppol-1)* &
714 1840092 : & (npw_k*nband_k(1)) : iband*npw_k+(isppol-1)*(npw_k*nband_k(1)))
715 : #else
716 : cwavef(:,1:npw_k)=cg(:,1+(iband-1)*npw_k+(isppol-1)* (npw_k*nband_k(1)) : iband*npw_k+(isppol-1)*(npw_k*nband_k(1)))
717 : #endif
718 :
719 : ! write(std_out,*)' iband : ',iband, ' isppol', isppol, ' -> index ', &
720 : ! & istate,index_state(iband+(isppol-1)*nband_k(1))
721 :
722 336 : tim_fourwf=14
723 : ! This call should be made by master, and then the results be sent to the other procs
724 :
725 : call fourwf(1,rhoaug,cwavef,dummy,wfraug,gbound,gbound,&
726 : & istwf_k,kg_k,kg_k,mgfftdiel,mpi_enreg,1,ngfftdiel,npw_k,1,ndiel4,ndiel5,ndiel6,&
727 336 : & 0,tim_fourwf,weight,weight,gpu_option=dtset%gpu_option)
728 :
729 : ! write(std_out,'(a,i5)')' After Fourier proc ',me_loc
730 :
731 : ! Fix the phase, and checks that the wavefunction is real
732 : ! (should be merged with routine fxphas)
733 336 : saa=zero ; sab=zero ; sbb=zero
734 12964 : do i3=1,ndiel3
735 501000 : do i2=1,ndiel2
736 22377136 : do i1=1,ndiel1
737 21876472 : saa=saa+wfraug(1,i1,i2,i3)**2
738 21876472 : sbb=sbb+wfraug(2,i1,i2,i3)**2
739 22364508 : sab=sab+wfraug(1,i1,i2,i3)*wfraug(2,i1,i2,i3)
740 : end do
741 : end do
742 : end do
743 :
744 336 : if(sbb>5.0d-9)then
745 :
746 : write(message, '(a,a,a,es20.10,a,i4,a,i2,a)' )&
747 0 : & 'The imaginary part of wavefunctions should be practically zero.',ch10,&
748 0 : & 'This is not the case, since sbb=',sbb,' for iband=',iband,'with sppol=',1,'.'
749 0 : ABI_WARNING(message)
750 0 : if(sbb>1.0d-7)then
751 0 : ABI_ERROR("sbb>1.0d-7")
752 : end if
753 : end if
754 :
755 : ! Possibility of writing to disk
756 :
757 23103220 : wfrspa(:,:,:,istate)=wfraug(1,:,:,:)
758 : end if !am_master
759 :
760 : end if
761 :
762 : ! End loop on iband
763 : end do
764 : ! End loop on nsppol
765 : end do
766 :
767 22 : if (am_master) then
768 18 : ABI_FREE(gbound)
769 18 : ABI_FREE(kg_k)
770 18 : ABI_FREE(cwavef)
771 18 : ABI_FREE(rhoaug)
772 18 : ABI_FREE(wfraug)
773 : #if defined HAVE_MPI
774 18 : ABI_FREE(cg_tmp)
775 : #else
776 : #endif
777 : end if
778 22 : ABI_FREE(flag_state_win)
779 :
780 : ! send wfrspa from master to world
781 22 : call xmpi_bcast(wfrspa,master,spaceComm,ierr)
782 : !call MPI_BCAST(wfrspa,nbuf,MPI_DOUBLE_PRECISION,master,spaceComm,ierr)
783 :
784 :
785 : !DEBUG
786 : !#if defined HAVE_MPI
787 : !call xmpi_barrier(spaceComm)
788 : !write(message,'(a)')' after iband loop synchronization done...'
789 : !call wrtout(std_out,message,'COLL')
790 : !#endif
791 : !ENDDEBUG
792 :
793 : !Compute the xc kernel, in the form needed for the singlet or triplet
794 : !excitation energy.
795 : !In the case ixc=20, kxc vanishes, but no change is made here, for simplicity.
796 : !(ixc=20 implemented only in the not spin polarized case)
797 :
798 : !DEBUG
799 : !write(std_out,*)' tddft : xc kernel '
800 : !do ifft=1,nkxc,41
801 : !write(std_out,*)ifft,kxc(ifft,1),kxc(ifft,2),kxc(ifft,3)
802 : !enddo
803 : !stop
804 : !ENDDEBUG
805 :
806 176 : ABI_MALLOC(kxc_for_tddft,(ndiel1,ndiel2,ndiel3,nsppol,nsppol,2/nsppol))
807 22 : if(dtset%ixc/=22)then
808 52 : do isppol=1,nsppol
809 103 : do jsppol=1,nsppol
810 : index=1
811 2195 : do i3=1,ndiel3
812 91227 : do i2=1,ndiel2
813 4541914 : do i1=1,ndiel1
814 9479766 : do sing_trip=1,2/nsppol
815 : kxc_for_tddft(i1,i2,i3,isppol,jsppol,sing_trip)=two/nsppol* &
816 9479766 : & (kxc(index,isppol+jsppol-1)-(sing_trip-1)*kxc(index,2))
817 : end do
818 4539801 : index=index+1
819 : end do
820 : end do
821 : end do
822 : end do
823 : end do
824 : else
825 : ! This is for the Burke-Petersilka-Gross hybrid, with ixc=22
826 : ! However, the implementation in case of spin-polarized system should not be expected to be the correct one !
827 2 : do isppol=1,nsppol
828 3 : do jsppol=1,nsppol
829 : index=1
830 32 : do i3=1,ndiel3
831 931 : do i2=1,ndiel2
832 27930 : do i1=1,ndiel1
833 81000 : do sing_trip=1,2/nsppol
834 81000 : kxc_for_tddft(i1,i2,i3,isppol,jsppol,sing_trip)=((-1)**(sing_trip+1))*kxc(index,2)
835 : end do
836 27900 : index=index+1
837 : end do
838 : end do
839 : end do
840 : end do
841 : end do
842 : end if
843 :
844 22 : pole_approx=0
845 :
846 66 : ABI_MALLOC(excit_coords,(nexcit_win**2,2))
847 :
848 : if (xmpi_paral==1) then
849 66 : ABI_MALLOC(counts,(0:nproc_loc-1))
850 44 : ABI_MALLOC(displs,(0:nproc_loc-1))
851 44 : ABI_MALLOC(recvcounts,(0:nproc_loc-1))
852 94 : ABI_MALLOC(recvbuf,(5-nsppol,nproc_loc-1))
853 : end if
854 :
855 : !DEBUG
856 : !write(std_out,*)'before first loop'
857 : !ENDDEBUG
858 :
859 : !0000000000000000000000000000000000000000000000000000000
860 : !check if matrix file fname_tdexcit exists on disk
861 : !if the file is present,calculation is a continuation
862 22 : if (am_master) then
863 18 : inquire(file=trim(dtfil%fnametmp_tdexcit),exist=file_exist)
864 : ! for direct access to excitation file
865 18 : if(nsppol==1)then
866 12 : inquire(iolength=recl) omega_tddft_casida(1,1,1,1,1,1,1),omega_tddft_casida(1,1,1,1,1,1,1), iexcit,jexcit
867 : else
868 6 : inquire(iolength=recl) omega_tddft_casida(1,1,1,1,1,1,1), iexcit,jexcit
869 : end if
870 :
871 18 : temp_unit2 = get_unit()
872 18 : open(temp_unit2, file=trim(dtfil%fnametmp_tdexcit),form='unformatted', recl=recl, access='DIRECT')
873 :
874 72 : ABI_MALLOC(done_excit,(nexcit_win,nexcit_win))
875 :
876 18 : if(file_exist)then
877 0 : write(std_out,*)'TDDFT continues from a previous run'
878 0 : rec=0
879 0 : do iexcit=1,nexcit_win
880 0 : iexcit2 = indarr(iexcit)
881 0 : isppol2 = min((iexcit2-1)/nexcit_pol(1) +1,2)
882 0 : iunocc2 = (iexcit2-(isppol2-1)*nexcit_pol(1)-1)/nband_occ(isppol2)+1
883 0 : iocc2 = iexcit2-(isppol2-1)*nexcit_pol(1)-(iunocc2-1)*nband_occ(isppol2)
884 0 : do jexcit=1,nexcit_win
885 0 : iexcit1 = indarr(jexcit)
886 0 : isppol1 = min((iexcit1-1)/nexcit_pol(1) +1,2)
887 0 : iunocc1 = (iexcit1-(isppol1-1)*nexcit_pol(1)-1)/nband_occ(isppol1)+1
888 0 : iocc1 = iexcit1-(isppol1-1)*nexcit_pol(1)-(iunocc1-1)*nband_occ(isppol1)
889 :
890 0 : rec=rec+1 ! record of the entry in the excitation file
891 0 : if(nsppol==1)then
892 0 : read(temp_unit2,rec=rec) omega_tddft_casida_dummy(1), omega_tddft_casida_dummy(2), idummy, jdummy
893 : else
894 0 : read(temp_unit2,rec=rec) omega_tddft_casida_dummy(1), idummy, jdummy
895 : end if
896 0 : done_excit(jexcit,iexcit)= ( idummy /= -1 .and. jdummy /= -1 ) ! if true, eigsqr_singlet and eigsqr_triplet are ok
897 : ! and a true is marked in the logical array done_excit
898 0 : if (done_excit(jexcit,iexcit)) then
899 0 : do sing_trip=1,2/nsppol
900 : omega_tddft_casida(iocc1,iunocc1,isppol1,iocc2,iunocc2,isppol2,sing_trip)= &
901 0 : & omega_tddft_casida_dummy(sing_trip)
902 : end do
903 : end if
904 : end do
905 : end do
906 :
907 : else
908 18 : write(std_out,*)'no excitation matrix on disk'
909 18 : write(std_out,*)'TDDFT starts from scratch'
910 31712 : done_excit = .false. ! initialize the logical array to false
911 18 : rec=0
912 628 : do iexcit=1,nexcit_win
913 31712 : do jexcit=1,nexcit_win
914 31084 : rec=rec+1
915 31694 : if(nsppol==1)then
916 5434 : write(temp_unit2,rec=rec) zero,zero,-1,-1
917 : else
918 25650 : write(temp_unit2,rec=rec) zero,-1,-1
919 : end if
920 : end do
921 : end do
922 : end if
923 :
924 : ! Need to list the elements to compute, taking the symmetry into account: valid only for Gamma point but this is already the case
925 18 : count_to_do=0
926 628 : do iexcit=1,nexcit_win
927 16475 : do jexcit=1,iexcit
928 16457 : if (.not. done_excit(jexcit,iexcit)) then
929 15847 : count_to_do=count_to_do+1
930 15847 : excit_coords(count_to_do,1)=iexcit
931 15847 : excit_coords(count_to_do,2)=jexcit
932 : end if
933 : end do
934 : end do
935 :
936 18 : ABI_FREE(done_excit)
937 :
938 :
939 :
940 : if (xmpi_paral==1) then
941 : ! Compute limits for load balancing
942 40 : do iproc=0,nproc_loc-1
943 22 : displs(iproc)=(iproc*count_to_do)/nproc_loc
944 40 : counts(iproc)=min(((iproc+1)*count_to_do)/nproc_loc,count_to_do)-displs(iproc)
945 : end do
946 : end if
947 :
948 : end if ! am_master
949 :
950 22 : call xmpi_bcast(count_to_do,master,spaceComm,ierr)
951 :
952 22 : displ=(me_loc*count_to_do)/nproc_loc
953 22 : count=min(((me_loc+1)*count_to_do)/nproc_loc,count_to_do)-displ
954 22 : displmax=((nproc_loc-1)*count_to_do)/nproc_loc
955 22 : countmax=count_to_do-displmax
956 :
957 22 : write(message,'(A,I6)') 'Maximum number of matrix elements per processor = ',countmax
958 22 : call wrtout(std_out,message,'COLL')
959 :
960 : if (xmpi_paral==1) then
961 : ! Need to dispatch the elements to compute to the different processes
962 66 : ABI_MALLOC(tmpbuf,(nexcit_win**2))
963 50706 : tmpbuf=0
964 22 : call xmpi_scatterv(excit_coords(:,1),counts,displs,tmpbuf,count,0,spaceComm,ierr)
965 50706 : excit_coords(:,1)=tmpbuf(:)
966 50706 : tmpbuf=0
967 22 : call xmpi_scatterv(excit_coords(:,2),counts,displs,tmpbuf,count,0,spaceComm,ierr)
968 50706 : excit_coords(:,2)=tmpbuf(:)
969 22 : ABI_FREE(tmpbuf)
970 : end if
971 :
972 22 : nfftdiel=ndiel1*ndiel2*ndiel3
973 110 : ABI_MALLOC(wfprod,(ndiel1,ndiel2,ndiel3))
974 66 : ABI_MALLOC(work,(nfftdiel))
975 66 : ABI_MALLOC(sexc,(3,nexcit_win))
976 66 : ABI_MALLOC(done_sexc,(nexcit_win))
977 : !ABI_MALLOC(done_sexc2,(nexcit_win))
978 66 : ABI_MALLOC(rhog,(2,nfftdiel))
979 44 : ABI_MALLOC(vhartr,(nfftdiel))
980 :
981 3582 : sexc(:,:)=zero
982 912 : done_sexc(:)=.false.
983 :
984 : !----------------------------------------------------------
985 : !Main double loop
986 :
987 : old_iexcit=0
988 15873 : do ijexcit=1,countmax
989 : ! we really loop only through count, but we need to go through countmax
990 : ! to make sure that all processes execute MPI_Gatherv below
991 15851 : if (ijexcit <= count) then
992 15847 : iexcit=excit_coords(ijexcit,1)
993 15847 : jexcit=excit_coords(ijexcit,2)
994 :
995 15847 : iexcit2 = indarr(iexcit)
996 15847 : isppol2 = min((iexcit2-1)/nexcit_pol(1) +1,2)
997 15847 : iunocc2 = (iexcit2-(isppol2-1)*nexcit_pol(1)-1)/nband_occ(isppol2)+1
998 15847 : iocc2 = iexcit2-(isppol2-1)*nexcit_pol(1)-(iunocc2-1)*nband_occ(isppol2)
999 :
1000 15847 : iexcit1 = indarr(jexcit)
1001 15847 : isppol1 = min((iexcit1-1)/nexcit_pol(1) +1,2)
1002 15847 : iunocc1 = (iexcit1-(isppol1-1)*nexcit_pol(1)-1)/nband_occ(isppol1)+1
1003 15847 : iocc1 = iexcit1-(isppol1-1)*nexcit_pol(1)-(iunocc1-1)*nband_occ(isppol1)
1004 :
1005 15847 : if (old_iexcit /= iexcit) then
1006 : ! We start a new column of the matrix
1007 : ! DEBUG
1008 : ! write(message,'(a,i5,a,i3)')'treating iexcit = ',iexcit,&
1009 : ! &' with proc ',me_loc
1010 : ! call wrtout(std_out,message,'PERS')
1011 : ! ENDDEBUG
1012 :
1013 : ! DEBUG
1014 : ! write(message,'(a,i3)')'Multiplicating phi1 phi2, on proc ',me_loc
1015 : ! call wrtout(std_out,message,'PERS')
1016 : ! ENDDEBUG
1017 :
1018 : ifft=1
1019 25607 : do i3=1,ndiel3
1020 1064228 : do i2=1,ndiel2
1021 52165476 : do i1=1,ndiel1
1022 : wfprod(i1,i2,i3)=wfrspa(i1,i2,i3,index_state(iocc2+(isppol2-1)*nband_k(1))) &
1023 51101862 : & *wfrspa(i1,i2,i3,index_state(iunocc2+nband_occ(isppol2)+(isppol2-1)*nband_k(1)))
1024 51101862 : work(ifft)=wfprod(i1,i2,i3)
1025 52140483 : ifft=ifft+1
1026 : end do
1027 : end do
1028 : end do
1029 614 : if (jexcit == 1) then
1030 25423 : do i3=1,ndiel3
1031 1055334 : do i2=1,ndiel2
1032 51719796 : do i1=1,ndiel1
1033 203688369 : do idir=1,3
1034 : poscart(idir)=rprimd(idir,1)*pos(i1,1)+&
1035 : & rprimd(idir,2)*pos(i2,2)+&
1036 151993386 : & rprimd(idir,3)*pos(i3,3)
1037 202657848 : sexc(idir,iexcit)=sexc(idir,iexcit)+poscart(idir)*wfprod(i1,i2,i3)
1038 : end do
1039 : end do
1040 : end do
1041 25423 : done_sexc(iexcit)=.true.
1042 : end do
1043 : end if
1044 :
1045 : ! For the singlet correction, must compute the hartre potential created
1046 : ! by the product of wavefunctions
1047 614 : cplex=1
1048 :
1049 : ! DEBUG
1050 : ! write(message,'(a,i3)')'Before Fourdp, on proc ',me_loc
1051 : ! call wrtout(std_out,message,'PERS')
1052 : ! ENDDEBUG
1053 :
1054 614 : call fourdp(cplex,rhog,work,-1,mpi_enreg,nfftdiel,1,ngfftdiel,0)
1055 :
1056 : ! DEBUG
1057 : ! write(message,'(a,i3)')'Before Hartree, on proc ',me_loc
1058 : ! call wrtout(std_out,message,'PERS')
1059 : ! write(std_out,*)'CPU ',me_loc,ch10,&
1060 : ! & ' cplex : ',cplex,ch10,&
1061 : ! & ' gmet(3,3) : ',gmet(3,3),ch10,&
1062 : ! & ' gsqcut : ',gsqcut,ch10,&
1063 : ! & ' rhog(1,1) :,',rhog(1,1),ch10,&
1064 : ! & ' vhartr(1) :,',vhartr(1)
1065 : ! ENDDEBUG
1066 :
1067 : call hartre(cplex,gsqcut,dtset%icutcoul,0,mpi_enreg,nfftdiel,ngfftdiel,&
1068 614 : &dtset%nkpt,dtset%rcut,rhog,rprimd,dtset%vcutgeo,vhartr)
1069 :
1070 : ! DEBUG
1071 : ! write(message,'(a,i3)')'After Hartree, on proc ',me_loc
1072 : ! call wrtout(std_out,message,'PERS')
1073 : ! ENDDEBUG
1074 : end if
1075 34521 : old_iexcit=iexcit
1076 :
1077 : ! DEBUG
1078 : ! write(std_out,*)' treating iexcit = ',jexcit
1079 : ! write(std_out,*)' indarr(iexcit) =',iexcit1,iocc1,iunocc1
1080 : ! write(std_out,*)' index ',index_state(iocc1+(isppol1-1)*nband_k(1)), &
1081 : ! & index_state(iunocc1+nband_occ(isppol1)+(isppol1-1)*nband_k(1))
1082 : ! ENDDEBUG
1083 :
1084 : if(pole_approx==0 .or. (iunocc1==iunocc2 .and. iocc1==iocc2 .and. isppol1==isppol2))then
1085 34521 : sum_kernel(:)=zero
1086 15847 : f_sing_trip(1)=two/dble(nsppol)
1087 15847 : if(nsppol==1) f_sing_trip(2)=zero
1088 : ! For Fermi-Amaldi kxc, the xc contribution is -1/2 the Hartree contribution
1089 : ! to the triplet state. The following factors combines both contributions.
1090 15847 : if(dtset%ixc==20 .or. dtset%ixc==22)then
1091 570 : if(nsppol==1)then
1092 570 : f_sing_trip(1)= one
1093 570 : f_sing_trip(2)=-one
1094 : end if
1095 : end if
1096 : ifft=1
1097 687541 : do i3=1,ndiel3
1098 29500219 : do i2=1,ndiel2
1099 1494141288 : do i1=1,ndiel1
1100 : myproduct=wfrspa(i1,i2,i3,index_state(iocc1+(isppol1-1)*nband_k(1))) &
1101 1464656916 : & *wfrspa(i1,i2,i3,index_state(iunocc1+nband_occ(isppol1)+(isppol1-1)*nband_k(1)))
1102 3163331268 : do sing_trip=1,2/nsppol
1103 : sum_kernel(sing_trip)=sum_kernel(sing_trip)+&
1104 : & myproduct*(f_sing_trip(sing_trip)*vhartr(ifft)+kxc_for_tddft(i1,i2,i3,isppol1,isppol2,sing_trip) &
1105 3163331268 : *wfprod(i1,i2,i3))
1106 : end do
1107 1493469594 : ifft=ifft+1
1108 : end do ! i1
1109 : end do ! i2
1110 : end do ! i3
1111 :
1112 : ! The factor two is coherent with the formulas of Vasiliev et al
1113 15847 : factor=two*sqrtks(iexcit1)*sqrtks(iexcit2)*normint
1114 34521 : do sing_trip=1,2/nsppol
1115 : omega_tddft_casida(iocc1,iunocc1,isppol1,iocc2,iunocc2,isppol2,sing_trip)= &
1116 : & omega_tddft_casida(iocc1,iunocc1,isppol1,iocc2,iunocc2,isppol2,sing_trip)+ &
1117 34521 : & factor*sum_kernel(sing_trip)
1118 : end do
1119 :
1120 : ! End condition of being diagonal element if pole approximation
1121 : end if
1122 :
1123 : ! Continue writing excitation matrix
1124 15847 : if (am_master) then
1125 : ! the master writes its results to disk
1126 12740 : if(nsppol==1)then
1127 : write(temp_unit2, rec=(iexcit-1)*nexcit_win+jexcit ) &
1128 2827 : & omega_tddft_casida(iocc1,iunocc1,isppol1,iocc2,iunocc2,isppol2,1),&
1129 2827 : & omega_tddft_casida(iocc1,iunocc1,isppol1,iocc2,iunocc2,isppol2,2),&
1130 5654 : & iexcit, jexcit
1131 : else
1132 : write(temp_unit2, rec=(iexcit-1)*nexcit_win+jexcit ) &
1133 9913 : & omega_tddft_casida(iocc1,iunocc1,isppol1,iocc2,iunocc2,isppol2,1), &
1134 19826 : & iexcit, jexcit
1135 : end if
1136 :
1137 : ! DEBUG
1138 : ! if(nsppol==1)then
1139 : ! write(std_out,*)'singlet: ',omega_tddft_casida(iocc1,iunocc1,isppol1,iocc2,iunocc2,isppol2,1),&
1140 : ! & 'iexcit: ',iexcit,'jexcit :',jexcit
1141 : ! write(std_out,*)'triplet: ',omega_tddft_casida(iocc1,iunocc1,isppol1,iocc2,iunocc2,isppol2,2),&
1142 : ! & 'iexcit: ',iexcit,'jexcit :',jexcit
1143 : ! else
1144 : ! write(std_out,*)'excitation: ',omega_tddft_casida(iocc1,iunocc1,isppol1,iocc2,iunocc2,isppol2,1),&
1145 : ! & 'iexcit: ',iexcit,'jexcit :',jexcit
1146 : ! endif
1147 : ! ENDDEBUG
1148 :
1149 12740 : sendcount=0
1150 : else
1151 3107 : sendcount=5-nsppol
1152 :
1153 3107 : if(nsppol==1)then
1154 : sendbuf=(/ omega_tddft_casida(iocc1,iunocc1,isppol1,iocc2,iunocc2,isppol2,1), &
1155 : & omega_tddft_casida(iocc1,iunocc1,isppol1,iocc2,iunocc2,isppol2,2), &
1156 0 : & real(iexcit,dp), real(jexcit,dp) /)
1157 : else
1158 : sendbuf=(/ omega_tddft_casida(iocc1,iunocc1,isppol1,iocc2,iunocc2,isppol2,1), &
1159 12428 : & real(iexcit,dp), real(jexcit,dp) /)
1160 : end if
1161 : end if ! am_master
1162 : else
1163 : ! ijexcit > count
1164 :
1165 : ! done with local work, so send message of zero length
1166 4 : sendcount=0
1167 :
1168 : end if ! ijexcit <= count
1169 :
1170 22 : if (xmpi_paral==1) then
1171 15851 : if (am_master) then
1172 :
1173 : ! Compute displacements and counts for the gathering of the results
1174 12742 : displs(0)=0
1175 12742 : recvcounts(0)=0
1176 15851 : do iproc=1,nproc_loc-1
1177 3109 : recvcounts(iproc)=min(((iproc+1)*count_to_do)/nproc_loc,count_to_do)-(iproc*count_to_do)/nproc_loc
1178 3109 : if (recvcounts(iproc) < countmax .and. ijexcit==countmax) then
1179 2 : recvcounts(iproc)=0
1180 : else
1181 3107 : recvcounts(iproc)=5-nsppol
1182 : end if
1183 15851 : displs(iproc)=displs(iproc-1)+recvcounts(iproc-1)
1184 : end do
1185 : end if
1186 :
1187 15851 : if (nproc_loc>1) then
1188 19896 : cptr=c_loc(recvbuf) ; call c_f_pointer(cptr,recvbuf_ptr,[size(recvbuf)])
1189 4974 : call xmpi_gatherv(sendbuf,sendcount,recvbuf_ptr,recvcounts,displs,0,spaceComm,ierr)
1190 : end if
1191 :
1192 15851 : if (am_master) then
1193 :
1194 : ! Extract eigsqr_singlet, eigsqr_triplet, iexcit, jexcit from receive buffer and
1195 : ! write to file
1196 31700 : do ijexcit2=1,sum(recvcounts)/(5-nsppol)
1197 3107 : iexcit=int(recvbuf(4-nsppol,ijexcit2))
1198 3107 : jexcit=int(recvbuf(5-nsppol,ijexcit2))
1199 :
1200 3107 : iexcit2 = indarr(iexcit)
1201 3107 : isppol2 = min((iexcit2-1)/nexcit_pol(1) +1,2)
1202 3107 : iunocc2 = (iexcit2-(isppol2-1)*nexcit_pol(1)-1)/nband_occ(isppol2)+1
1203 3107 : iocc2 = iexcit2-(isppol2-1)*nexcit_pol(1)-(iunocc2-1)*nband_occ(isppol2)
1204 :
1205 3107 : iexcit1 = indarr(jexcit)
1206 3107 : isppol1 = min((iexcit1-1)/nexcit_pol(1) +1,2)
1207 3107 : iunocc1 = (iexcit1-(isppol1-1)*nexcit_pol(1)-1)/nband_occ(isppol1)+1
1208 3107 : iocc1 = iexcit1-(isppol1-1)*nexcit_pol(1)-(iunocc1-1)*nband_occ(isppol1)
1209 :
1210 6214 : do sing_trip=1,2/nsppol
1211 : omega_tddft_casida(iocc1,iunocc1,isppol1,iocc2,iunocc2,isppol2,sing_trip)= &
1212 6214 : & recvbuf(sing_trip,ijexcit2)
1213 : end do
1214 :
1215 15849 : if(nsppol==1)then
1216 : write(temp_unit2, rec=(iexcit-1)*nexcit_win+jexcit ) &
1217 0 : & omega_tddft_casida(iocc1,iunocc1,isppol1,iocc2,iunocc2,isppol2,1), &
1218 0 : & omega_tddft_casida(iocc1,iunocc1,isppol1,iocc2,iunocc2,isppol2,2), &
1219 0 : & iexcit, jexcit
1220 : else
1221 : write(temp_unit2, rec=(iexcit-1)*nexcit_win+jexcit ) &
1222 3107 : omega_tddft_casida(iocc1,iunocc1,isppol1,iocc2,iunocc2,isppol2,1), &
1223 6214 : & iexcit, jexcit
1224 : end if
1225 : ! DEBUG
1226 : ! if(nsppol==1)then
1227 : ! write(std_out,*)'singlet: ',omega_tddft_casida(iocc1,iunocc1,isppol1,iocc2,iunocc2,isppol2,1), &
1228 : ! & 'iexcit: ',iexcit,'jexcit :',jexcit
1229 : ! write(std_out,*)'triplet: ',omega_tddft_casida(iocc1,iunocc1,isppol1,iocc2,iunocc2,isppol2,2),
1230 : ! & 'iexcit: ',iexcit,'jexcit :',jexcit
1231 : ! else
1232 : ! write(std_out,*)'excitation: ',omega_tddft_casida(iocc1,iunocc1,isppol1,iocc2,iunocc2,isppol2,1), &
1233 : ! & 'iexcit: ',iexcit,'jexcit :',jexcit
1234 : ! endif
1235 : ! ENDDEBUG
1236 :
1237 : end do
1238 :
1239 : end if
1240 : end if ! MPI parallel
1241 :
1242 : ! End indices loops
1243 : end do ! ijexcit
1244 :
1245 : !End of the main double loop
1246 : !--------------------------------------------------------------------
1247 :
1248 :
1249 : if (xmpi_paral==1) then
1250 : ! sexc needs to be summed here since it used only by master
1251 22 : call xmpi_barrier(spaceComm)
1252 : ! call xmpi_sum_master(sexc,master,spaceComm,ierr) ! Does not work on some machines
1253 22 : call xmpi_sum(sexc,spaceComm,ierr)
1254 22 : call xmpi_lor(done_sexc,spaceComm)
1255 : !done_sexc2=done_sexc
1256 : !call MPI_Reduce(done_sexc2,done_sexc,nexcit_win,MPI_LOGICAL,MPI_LOR,master,spaceComm,ierr)
1257 : end if
1258 :
1259 :
1260 22 : if (am_master) then
1261 : ! We compute sexc again if it was not done. Will only be executed if
1262 : ! there was a restart from values read from logical unit temp_unit2.
1263 :
1264 628 : do iexcit=1,nexcit_win
1265 :
1266 : ! DEBUG
1267 : ! write(std_out,*)'do on excitation',iexcit
1268 : ! END DEBUG
1269 :
1270 628 : if (.not.done_sexc(iexcit)) then
1271 0 : iexcit2 = indarr(iexcit)
1272 0 : isppol2 = min((iexcit2-1)/nexcit_pol(1) +1,2)
1273 0 : iunocc2 = (iexcit1-(isppol2-1)*nexcit_pol(1)-1)/nband_occ(isppol2)+1
1274 0 : iocc2 = iexcit2-(isppol2-1)*nexcit_pol(1)-(iunocc2-1)*nband_occ(isppol2)
1275 0 : do i3=1,ndiel3
1276 0 : do i2=1,ndiel2
1277 0 : do i1=1,ndiel1
1278 : wfprod(i1,i2,i3)=wfrspa(i1,i2,i3,index_state(iocc2+(isppol2-1)*nband_k(1))) &
1279 0 : & *wfrspa(i1,i2,i3,index_state(iunocc2+nband_occ(isppol2)+ (isppol2-1)*nband_k(1)))
1280 0 : do idir=1,3
1281 : poscart(idir)=rprimd(idir,1)*pos(i1,1)+&
1282 : & rprimd(idir,2)*pos(i2,2)+&
1283 0 : & rprimd(idir,3)*pos(i3,3)
1284 0 : sexc(idir,iexcit)=sexc(idir,iexcit)+poscart(idir)*wfprod(i1,i2,i3)
1285 : end do
1286 : end do
1287 : end do
1288 : end do
1289 : end if
1290 : end do
1291 : end if
1292 :
1293 22 : ABI_FREE(work)
1294 22 : ABI_FREE(rhog)
1295 22 : ABI_FREE(pos)
1296 22 : ABI_FREE(vhartr)
1297 22 : ABI_FREE(kxc_for_tddft)
1298 22 : ABI_FREE(wfprod)
1299 22 : ABI_FREE(index_state)
1300 22 : ABI_FREE(excit_coords)
1301 22 : ABI_FREE(wfrspa)
1302 :
1303 : !Write the first excitation energies
1304 22 : write(message, '(a,a,es18.8,a,a,a,a,a,a,a,a,a)' )ch10,&
1305 22 : & ' Ground state total energy (Ha) :',etotal,ch10,ch10,&
1306 22 : & ' Kohn-Sham energy differences,',ch10,&
1307 22 : & ' corresponding total energies and oscillator strengths (X,Y,Z and average)-',ch10,&
1308 22 : & ' (oscillator strengths smaller than 1.e-6 are set to zero)',ch10,&
1309 44 : & ' Transition (Ha) and (eV) Tot. Ene. (Ha) Aver XX YY ZZ'
1310 22 : call wrtout(ab_out,message,'COLL')
1311 22 : call wrtout(std_out,message,'COLL')
1312 :
1313 : if (xmpi_paral==1) then
1314 22 : ABI_FREE(counts)
1315 22 : ABI_FREE(displs)
1316 22 : ABI_FREE(recvbuf)
1317 22 : ABI_FREE(recvcounts)
1318 : end if
1319 :
1320 22 : if (am_master) then
1321 :
1322 54 : ABI_MALLOC(osc_str,(7,nexcit))
1323 :
1324 628 : do iexcit=1,nexcit_win
1325 610 : iexcit2 = indarr(iexcit)
1326 610 : isppol = min((iexcit2-1)/nexcit_pol(1) +1,2)
1327 610 : iunocc = (iexcit2-(isppol-1)*nexcit_pol(1)-1)/nband_occ(isppol)+1
1328 610 : iocc = iexcit2-(isppol-1)*nexcit_pol(1)-(iunocc-1)*nband_occ(isppol)
1329 :
1330 610 : osc_str(1,iexcit)=zero
1331 2440 : do idir=1,3
1332 : ! One of the factor of two comes from the spin degeneracy,
1333 : ! the other comes from Eq.(40) of Casida
1334 : osc_str(idir+1,iexcit)=&
1335 1830 : & (sexc(idir,iexcit)*sqrtks(iexcit2)*normint*ucvol)**2*two*two/nsppol
1336 : osc_str(1,iexcit)=osc_str(1,iexcit)&
1337 2440 : & +osc_str(idir+1,iexcit)*third
1338 : end do
1339 3050 : do ii=1,4
1340 3050 : if(abs(osc_str(ii,iexcit))<tol6)osc_str(ii,iexcit)=zero
1341 : end do
1342 : ! Changed, the whole spectrum is written
1343 : ! The array eexcit has been reordered previously, the others also
1344 610 : if(nsppol==1)then
1345 : write(message, '(i4,a,i3,2es12.5,es13.5,es11.4,3es9.2)' ) &
1346 220 : & iocc,'->',iunocc+nband_occ(isppol), &
1347 220 : & eexcit(iexcit), eexcit(iexcit)*Ha_eV,eexcit(iexcit)+etotal, &
1348 : ! XG 020209 : Jean-Yves, I assume that the printout of sexc is for debugging ?!
1349 : ! & osc_str(1:4,iexcit),sexc(1:3,iexcit)
1350 440 : & osc_str(1:4,iexcit)
1351 : else
1352 : write(message, '(i4,a,i3,a,i1,2es12.5,es13.5,es11.4,3es9.2)' ) &
1353 390 : & iocc,'->',iunocc+nband_occ(isppol),' s:',isppol, &
1354 390 : & eexcit(iexcit), eexcit(iexcit)*Ha_eV,eexcit(iexcit)+etotal, &
1355 780 : & osc_str(1:4,iexcit)
1356 : end if
1357 610 : call wrtout(ab_out,message,'COLL')
1358 628 : call wrtout(std_out,message,'COLL')
1359 :
1360 : end do
1361 :
1362 : ! Check of the Sum rule for Casida eq.47,
1363 : ! only exact if complete basis of excitations, as well as local potentials only.
1364 18 : sumx=zero
1365 628 : do iexcit=1,nexcit_win
1366 628 : sumx=sumx+osc_str(1,iexcit)
1367 : end do
1368 18 : write(message, '(a,es16.6)' )' Sum of osc. strength : ',sumx
1369 18 : call wrtout(ab_out,message,'COLL')
1370 18 : call wrtout(std_out,message,'COLL')
1371 :
1372 : ! -Diagonalize the excitation matrices----------------------------
1373 :
1374 54 : ABI_MALLOC(eexcit2,(nexcit_win))
1375 72 : ABI_MALLOC(vec,(2,nexcit_win,nexcit_win))
1376 :
1377 48 : do sing_trip=1,2/nsppol
1378 :
1379 : if(pole_approx==0)then
1380 :
1381 90 : ABI_MALLOC(matr,(nexcit_win*(nexcit_win+1)))
1382 90 : ABI_MALLOC(zhpev1,(2,2*nexcit_win-1))
1383 90 : ABI_MALLOC(zhpev2,(3*nexcit_win-2))
1384 37378 : matr(:)=zero
1385 30 : ier=0
1386 : ! DEBUG
1387 : ! write(std_out,*)' after allocation matrices '
1388 : ! ENDDEBUG
1389 :
1390 :
1391 : ! Store the matrix in proper mode before calling zhpev
1392 :
1393 30 : index=1
1394 860 : do iexcit=1,nexcit_win
1395 830 : iexcit2 = indarr(iexcit)
1396 830 : isppol2 = min((iexcit2-1)/nexcit_pol(1) +1,2)
1397 830 : iunocc2 = (iexcit2-(isppol2-1)*nexcit_pol(1)-1)/nband_occ(isppol2)+1
1398 830 : iocc2 = iexcit2-(isppol2-1)*nexcit_pol(1)-(iunocc2-1)*nband_occ(isppol2)
1399 19534 : do jexcit=1,iexcit
1400 18674 : iexcit1 = indarr(jexcit)
1401 18674 : isppol1 = min((iexcit1-1)/nexcit_pol(1) +1,2)
1402 18674 : iunocc1 = (iexcit1-(isppol1-1)*nexcit_pol(1)-1)/nband_occ(isppol1)+1
1403 18674 : iocc1 = iexcit1-(isppol1-1)*nexcit_pol(1)-(iunocc1-1)*nband_occ(isppol1)
1404 :
1405 18674 : matr(index)=omega_tddft_casida(iocc1,iunocc1,isppol1,iocc2,iunocc2,isppol2,sing_trip)
1406 18674 : matr(index+1)=zero
1407 19504 : index=index+2
1408 : end do
1409 :
1410 : end do
1411 :
1412 : ! DEBUG
1413 : ! write(std_out,*)' after filling matrices '
1414 : ! ENDDEBUG
1415 :
1416 : call ZHPEV ('V','U',nexcit_win,matr,eexcit2,vec,nexcit_win,zhpev1,&
1417 30 : & zhpev2,ier)
1418 :
1419 30 : ABI_FREE(matr)
1420 30 : ABI_FREE(zhpev1)
1421 30 : ABI_FREE(zhpev2)
1422 : ! DEBUG
1423 : ! write(std_out,*)' after deallocating matrices '
1424 : ! ENDDEBUG
1425 :
1426 :
1427 : else
1428 :
1429 : vec(:,:,:)=zero
1430 : do isppol=1, nsppol
1431 : do iunocc=1,nband_k(isppol)
1432 : do iocc=1,nband_k(isppol)
1433 : index=iocc+nband_k(isppol)*(iunocc-1)+(isppol-1)*nexcit_pol(1)
1434 : eexcit2(index)=omega_tddft_casida(iocc,iunocc,isppol,iocc,iunocc,isppol,sing_trip)
1435 : vec(1,index,index)=one
1436 : end do
1437 : end do
1438 : end do
1439 :
1440 : end if
1441 :
1442 : ! Compute the excitation energies from the square root of eexcit2
1443 : ! eexcit(:)=sqrt(eexcit2(:)
1444 :
1445 30 : ABI_FREE(eexcit)
1446 90 : ABI_MALLOC(eexcit,(nexcit_win))
1447 :
1448 860 : eexcit(:)=sqrt(dabs(eexcit2(:))+tol10**2)
1449 : ! Write the first excitation energies
1450 30 : if(sing_trip==1)then
1451 18 : if(nsppol==1)then
1452 12 : write(message, '(a,a,a,a,a,a)' )ch10,&
1453 12 : & ' TDDFT singlet excitation energies (at most 20 of them are printed),',ch10,&
1454 12 : & ' and corresponding total energies. ',ch10,&
1455 24 : & ' Excit# (Ha) and (eV) total energy (Ha) major contributions '
1456 : else
1457 6 : write(message, '(a,a,a,a,a,a)' )ch10,&
1458 6 : & ' TDDFT mixed excitation energies (at most 40 of them are printed),',ch10,&
1459 6 : & ' and corresponding total energies. ',ch10,&
1460 12 : & ' Excit# (Ha) and (eV) total energy (Ha) major contributions '
1461 : end if
1462 : else
1463 12 : write(message, '(a,a,a,a,a,a)' )ch10,&
1464 12 : & ' TDDFT triplet excitation energies (at most 20 of them are printed),',ch10,&
1465 12 : & ' and corresponding total energies. ',ch10,&
1466 24 : & ' Excit# (Ha) and (eV) total energy (Ha) major contributions '
1467 : end if
1468 30 : call wrtout(ab_out,message,'COLL')
1469 30 : call wrtout(std_out,message,'COLL')
1470 30 : call wrtout(std_out,' tddft : before iexcit loop',"COLL")
1471 :
1472 500 : do iexcit=1,min(nexcit_win,nexcitout)
1473 470 : write(std_out,*)' tddft : iexcit=',iexcit
1474 : ! Select largest and next contributions
1475 470 : flargest=zero ; fnext=zero
1476 470 : nlargest=0 ; nnext=0
1477 470 : if(nsppol==2)then
1478 120 : fnext1=zero ; fnext2=zero
1479 120 : nnext1=0 ; nnext2=0
1480 : end if
1481 15988 : do jexcit=1,nexcit_win
1482 15518 : ff=vec(1,jexcit,iexcit)**2+vec(2,jexcit,iexcit)**2
1483 15988 : if(ff>flargest+tol12)then
1484 1475 : if(nsppol==2)then
1485 392 : nnext2=nnext1 ; fnext2=fnext1
1486 392 : nnext1=nnext ; fnext1=fnext
1487 : end if
1488 1475 : nnext=nlargest ; fnext=flargest
1489 1475 : nlargest=indarr(jexcit) ; flargest=ff
1490 14043 : else if(ff>fnext+tol12)then
1491 825 : if(nsppol==2)then
1492 256 : nnext2=nnext1 ; fnext2=fnext1
1493 256 : nnext1=nnext ; fnext1=fnext
1494 : end if
1495 825 : nnext=indarr(jexcit) ; fnext=ff
1496 13218 : else if(nsppol==2)then
1497 7152 : if(ff>fnext1+tol12)then
1498 228 : nnext2=nnext1 ; fnext2=fnext1
1499 228 : nnext1=indarr(jexcit) ; fnext1=ff
1500 6924 : else if(ff>fnext2+tol12)then
1501 224 : nnext2=indarr(jexcit) ; fnext2=ff
1502 : end if
1503 : end if
1504 :
1505 : end do
1506 :
1507 470 : isppol_l = min((nlargest-1)/nexcit_pol(1) +1,nsppol)
1508 470 : iunocc_l = (nlargest-(isppol_l-1)*nexcit_pol(1)-1)/nband_occ(isppol_l)+1
1509 470 : iocc_l = nlargest-(isppol_l-1)*nexcit_pol(1)-(iunocc_l-1)*nband_occ(isppol_l)
1510 470 : isppol_n = min((nnext-1)/nexcit_pol(1) +1,nsppol)
1511 470 : iunocc_n = (nnext-(isppol_n-1)*nexcit_pol(1)-1)/nband_occ(isppol_n)+1
1512 470 : iocc_n = nnext-(isppol_n-1)*nexcit_pol(1)-(iunocc_n-1)*nband_occ(isppol_n)
1513 470 : if(nsppol==2)then
1514 120 : isppol_n1 = min((nnext1-1)/nexcit_pol(1) +1,nsppol)
1515 120 : iunocc_n1 = (nnext1-(isppol_n1-1)*nexcit_pol(1)-1)/nband_occ(isppol_n1)+1
1516 120 : iocc_n1 = nnext1-(isppol_n1-1)*nexcit_pol(1)-(iunocc_n1-1)*nband_occ(isppol_n1)
1517 120 : isppol_n2 = min((nnext2-1)/nexcit_pol(1) +1,nsppol)
1518 120 : iunocc_n2 = (nnext2-(isppol_n2-1)*nexcit_pol(1)-1)/nband_occ(isppol_n2)+1
1519 120 : iocc_n2 = nnext2-(isppol_n2-1)*nexcit_pol(1)-(iunocc_n2-1)*nband_occ(isppol_n2)
1520 : end if
1521 :
1522 500 : if(nsppol==1)then
1523 : write(message,'(i4,es15.5,es14.5,es16.6,f8.2,a,i3,a,i3,a,f6.2,a,i3,a,i3,a)') &
1524 350 : & iexcit,eexcit(iexcit),&
1525 350 : & eexcit(iexcit)*Ha_eV,eexcit(iexcit)+etotal,&
1526 350 : & flargest,'(',iocc_l,'->',iunocc_l+nband_occ(1),')',&
1527 700 : & fnext, '(',iocc_n,'->',iunocc_n+nband_occ(1),')'
1528 350 : call wrtout(ab_out,message,'COLL')
1529 350 : call wrtout(std_out,message,'COLL')
1530 : else
1531 120 : write(chain1,'(f8.2,a,i3,a,i3,a)')flargest,'(',iocc_l,'->',iunocc_l+nband_occ(isppol_l),')'
1532 120 : write(chain2,'(f8.2,a,i3,a,i3,a)')fnext,'(',iocc_n,'->',iunocc_n+nband_occ(isppol_n),')'
1533 120 : if(trim(chain1)==trim(chain2))then
1534 : write(message,'(i4,es15.5,es14.5,es16.6,a,a,a,a)') &
1535 80 : & iexcit,eexcit(iexcit),&
1536 160 : & eexcit(iexcit)*Ha_eV,eexcit(iexcit)+etotal,trim(chain1),'(1)',trim(chain2),'(2)'
1537 : else
1538 : write(message,'(i4,es15.5,es14.5,es16.6,a,a,i1,a,a,a,i1,a)') &
1539 40 : & iexcit,eexcit(iexcit),&
1540 40 : & eexcit(iexcit)*Ha_eV,eexcit(iexcit)+etotal,trim(chain1),'(',isppol_l,')',&
1541 80 : & trim(chain2),'(',isppol_n,')'
1542 : end if
1543 120 : call wrtout(ab_out,message,'COLL')
1544 120 : call wrtout(std_out,message,'COLL')
1545 120 : write(chain1,'(f8.2,a,i3,a,i3,a)')fnext1,'(',iocc_n1,'->',iunocc_n1+nband_occ(isppol_n1),')'
1546 120 : write(chain2,'(f8.2,a,i3,a,i3,a)')fnext2,'(',iocc_n2,'->',iunocc_n2+nband_occ(isppol_n2),')'
1547 120 : if(trim(chain1)==trim(chain2))then
1548 : write(message,'(a,a,a,a,a)' ) &
1549 80 : & ' ',&
1550 160 : & chain1,'(1)',chain2,'(2)'
1551 : else
1552 : write(message,'(a,a,a,i1,a,a,a,i1,a)' ) &
1553 40 : & ' ',&
1554 40 : & chain1,'(',isppol_n1,')',&
1555 80 : & chain2,'(',isppol_n2,')'
1556 : end if
1557 120 : call wrtout(ab_out,message,'COLL')
1558 120 : call wrtout(std_out,message,'COLL')
1559 : end if
1560 : end do
1561 :
1562 : ! For each iexcit excitation, compute the oscillator strength (Casida, eq 47)
1563 30 : write(message, '(a,a,a,a)' )ch10,&
1564 30 : & ' Oscillator strengths : (elements smaller than 1.e-6 are set to zero)',ch10,&
1565 60 : & ' Excit# (Ha) Average XX YY ZZ XY XZ YZ'
1566 30 : call wrtout(ab_out,message,'COLL')
1567 30 : call wrtout(std_out,message,'COLL')
1568 :
1569 860 : do iexcit=1,nexcit_win
1570 :
1571 : ! One of the factor of two comes from the spin degeneracy,
1572 : ! the other comes from Eq.(40) of Casida
1573 830 : factor=(normint*ucvol)**2*two*two/nsppol
1574 :
1575 6640 : osc_str(:,iexcit)=zero
1576 :
1577 : ! factor=(normint*ucvol)**2*two*two/nsppol
1578 :
1579 37348 : do jexcit=1,nexcit_win
1580 36518 : jexcit_cbase=indarr(jexcit)
1581 146902 : do idir=1,3
1582 : osc_str(idir+1,iexcit)=osc_str(idir+1,iexcit)+ &
1583 : & sexc(idir,jexcit)*sqrtks(jexcit_cbase)*sqrt(factor)* &
1584 146072 : & vec(1,jexcit,iexcit)
1585 : end do ! idir
1586 : end do ! jexcit
1587 :
1588 : ! The "standard" definition of the oscillator strength is the square
1589 : ! of the matrix elements.
1590 : ! So, instead of the coding
1591 : ! do idir=1,3
1592 : ! osc_str(1,iexcit)=osc_str(1,iexcit)+osc_str(idir+1,iexcit)**2*third
1593 : ! enddo
1594 : ! I think that the following is more "standard"
1595 : ! Now, osc_str(2:4,iexcit) are the X, Y and Z matrix elements, not
1596 : ! yet the oscillator strengths
1597 830 : osc_str(5,iexcit)=osc_str(2,iexcit)*osc_str(3,iexcit) ! off diag XY
1598 830 : osc_str(6,iexcit)=osc_str(2,iexcit)*osc_str(4,iexcit) ! off diag XZ
1599 830 : osc_str(7,iexcit)=osc_str(3,iexcit)*osc_str(4,iexcit) ! off diag ZZ
1600 3320 : do idir=1,3
1601 : ! Here the X,Y, and Z matrix elements are combined to give diagonal osc. strengths
1602 2490 : osc_str(idir+1,iexcit)=osc_str(idir+1,iexcit)**2
1603 3320 : osc_str(1,iexcit)=osc_str(1,iexcit)+osc_str(idir+1,iexcit)*third ! compute the trace
1604 : end do
1605 : ! At this stage, osc_str(1,iexcit) is exactly the same as from your coding
1606 : ! ***End of section to be checked
1607 :
1608 6640 : do ii=1,7
1609 6640 : if(abs(osc_str(ii,iexcit))<tol6)osc_str(ii,iexcit)=zero
1610 : end do
1611 : ! XG 020209 : Jean-Yves, the off-diagonal oscillator strengths
1612 : ! can become negative. It is important for automatic
1613 : ! checking that the numbers are separated by a blank, even
1614 : ! if they are negative. So replace the following format, to have at least one blank.
1615 830 : write(message, '(i4,es12.5,es10.3,3es10.3,3es10.2)' )iexcit,eexcit(iexcit),osc_str(1:7,iexcit)
1616 830 : call wrtout(ab_out,message,'COLL')
1617 860 : call wrtout(std_out,message,'COLL')
1618 : end do
1619 :
1620 : ! Check of the Sum rule for Casida eq.47,
1621 : ! only exact if complete basis of excitations, as well as local potentials only.
1622 30 : sumx=zero
1623 860 : do iexcit=1,nexcit_win
1624 860 : sumx=sumx+osc_str(1,iexcit)
1625 : end do
1626 30 : write(message, '(a,es16.6)' )' Sum of osc. strength : ',sumx
1627 30 : call wrtout(ab_out,message,'COLL')
1628 30 : call wrtout(std_out,message,'COLL')
1629 :
1630 : ! If singlet, compute Cauchy coefficients
1631 48 : if(sing_trip==1.AND.nsppol==1)then
1632 12 : cauchy(:)=zero
1633 232 : do iexcit=1,nexcit_win
1634 220 : fr_invsquare=one/(eexcit(iexcit)**2)
1635 220 : fr_power=one
1636 1772 : do ii=1,7
1637 1540 : fr_power=fr_power*fr_invsquare
1638 1760 : cauchy(ii)=cauchy(ii)+osc_str(1,iexcit)*fr_power
1639 : end do
1640 : end do
1641 : write(message, '(a,es11.3,a,es11.3,a,es11.3,a,a,es11.3,a,es11.3,a,es11.3,a,es11.3)' ) &
1642 12 : & ' Cauchy coeffs (au) : ( -2)->',cauchy(1),&
1643 12 : & ', ( -4)->',cauchy(2),', ( -6)->',cauchy(3),ch10,&
1644 24 : & ' (-8)->',cauchy(4),', (-10)->',cauchy(5),', (-12)->',cauchy(6),', (-14)->',cauchy(7)
1645 12 : call wrtout(ab_out,message,'COLL')
1646 12 : call wrtout(std_out,message,'COLL')
1647 : end if
1648 :
1649 : ! End the loop on singlet or triplet
1650 : end do
1651 :
1652 18 : ABI_FREE(eexcit2)
1653 18 : ABI_FREE(vec)
1654 18 : ABI_FREE(osc_str)
1655 :
1656 : !! The temporary files should be deleted at the end of this routine
1657 18 : close(temp_unit2,status='delete')
1658 18 : call timab(95,2,tsec)
1659 : end if ! end of am_master
1660 :
1661 22 : ABI_FREE(omega_tddft_casida)
1662 22 : ABI_FREE(eexcit)
1663 22 : ABI_FREE(sqrtks)
1664 22 : ABI_FREE(sexc)
1665 22 : ABI_FREE(done_sexc)
1666 22 : ABI_FREE(indarr)
1667 : !ABI_FREE(done_sexc2)
1668 :
1669 110 : end subroutine tddft
1670 : !!***
1671 :
1672 : end module m_tddft
1673 : !!***
|
