Branch data 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 : : !!***
|
