Line data Source code
1 : !!****m* ABINIT/m_rttddft_tdks
2 : !! NAME
3 : !! m_rttddft_tdks
4 : !!
5 : !! FUNCTION
6 : !! Contains the main object (tdks) to propagate
7 : !! the time-dependent Kohn-Sham equations in RT-TDDFT
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2021-2026 ABINIT group (FB)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : module m_rttddft_tdks
24 :
25 : use defs_basis
26 : use defs_abitypes, only: MPI_type
27 : use defs_datatypes, only: pseudopotential_type
28 : use defs_wvltypes, only: wvl_data, nullify_wvl_data
29 : use libxc_functionals, only: libxc_functionals_get_hybridparams
30 : use m_bandfft_kpt, only: bandfft_kpt, bandfft_kpt_init1, bandfft_kpt_destroy_array
31 : use m_cgprj, only: ctocprj
32 : use m_common, only: setup1
33 : use m_dtfil, only: datafiles_type
34 : use m_dtset, only: dataset_type
35 : use m_ebands, only: ebands_t, unpack_eneocc
36 : use m_energies, only: energies_type
37 : use m_errors, only: msg_hndl, assert
38 : use m_extfpmd, only: extfpmd_type
39 : use m_gemm_nonlop_projectors, only: init_gemm_nonlop, destroy_gemm_nonlop
40 : use m_geometry, only: fixsym
41 : use m_hdr, only: hdr_type
42 : use m_initylmg, only: initylmg
43 : use m_invovl, only: init_invovl, destroy_invovl
44 : use m_io_tools, only: open_file
45 : use m_inwffil, only: inwffil
46 : use m_kg, only: kpgio, getph, getcut
47 : use m_mpinfo, only: proc_distrb_cycle
48 : use m_occ, only: newocc
49 : use m_paw_an, only: paw_an_type, paw_an_init, paw_an_free, &
50 : paw_an_nullify
51 : use m_pawang, only: pawang_type
52 : use m_pawcprj, only: pawcprj_type,pawcprj_free,pawcprj_alloc, &
53 : pawcprj_getdim
54 : use m_paw_dmft, only: init_sc_dmft,destroy_sc_dmft,paw_dmft_type
55 : use m_pawfgr, only: pawfgr_type, pawfgr_init, pawfgr_destroy
56 : use m_pawfgrtab, only: pawfgrtab_type, pawfgrtab_init, pawfgrtab_free
57 : use m_paw_init, only: pawinit,paw_gencond
58 : use m_paw_ij, only: paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify
59 : use m_paw_nhat, only: nhatgrid
60 : use m_paw_occupancies, only: initrhoij
61 : use m_pawrad, only: pawrad_type
62 : use m_pawrhoij, only: pawrhoij_type, pawrhoij_copy, pawrhoij_free
63 : use m_paw_sphharm, only: setsym_ylm
64 : use m_pawtab, only: pawtab_type, pawtab_get_lsize
65 : use m_paw_tools, only: chkpawovlp
66 : use m_pawxc, only: pawxc_get_usekden
67 : use m_pspini, only: pspini
68 : use m_profiling_abi, only: abimem_record
69 : use m_rttddft_tdef, only: tdef_type
70 : use m_spacepar, only: setsym
71 : use m_specialmsg, only: wrtout
72 : use m_symtk, only: symmetrize_xred
73 : use m_wffile, only: wffile_type, WffClose
74 : use m_xmpi, only: xmpi_bcast, xmpi_sum
75 : use m_drivexc, only: xc_need_kden
76 :
77 : implicit none
78 :
79 : private
80 : !!***
81 :
82 : !! NAME
83 : !! tdks_type: Time Dependent Kohn-Sham type
84 : !! Object containing the TD KS orbitals and all other
85 : !! important variables required to run RT-TDDFT
86 : !!
87 : !! SOURCE
88 : type,public :: tdks_type
89 :
90 : !scalars
91 : integer :: bantot !total number of bands
92 : integer :: first_step !start propagation from first_step (for restart)
93 : integer :: mband_cprj !nb of band per proc (for cprj)
94 : integer :: mcg !nb of WFs (cg) coeffs
95 : integer :: mcprj !nb of cprj (projectors applied to WF)
96 : integer :: nfftf !nb of FFT grid pts (fine grid)
97 : integer :: nfft !nb of FFT grid pts (coarse grid)
98 : integer :: nhatgrdim !dimension of nhatgr array
99 : integer :: ngrvdw !dimension of grvdw array
100 : integer :: ntime !max nb of time steps
101 : integer :: current_unit!unit nb of the current density file
102 : integer :: tdener_unit !unit nb of the energy file
103 : integer :: tdef_unit !unit nb of the efield file
104 : integer :: tdrestart_unit !unit nb of the restart file
105 : integer :: unpaw !paw data tmp file unit
106 : integer :: usexcnhat !use nhat in the computation of the XC term
107 : real(dp) :: dt !propagation time step
108 : real(dp) :: ecore !core energy
109 : real(dp) :: etot !total energy
110 : real(dp) :: gsqcut !cut-off on G^2
111 : real(dp) :: ucvol !primitive cell volume
112 : real(dp) :: zion !total ionic charge
113 : logical :: gemm_nonlop_use_gemm !use efficient BLAS call
114 : !for computing non local potential
115 : type(energies_type) :: energies !contains various energy values
116 : type(hdr_type) :: hdr !header: contains various info
117 : type(paw_dmft_type) :: paw_dmft !paw_dmft object (unused but
118 : !required by various routines)
119 : type(pawfgr_type) :: pawfgr !FFT fine grid in PAW sphere
120 : type(pawang_type),pointer :: pawang => NULL() !angular grid in PAW sphere
121 : type(tdef_type) :: tdef !Object containing variables related to TD electric field
122 : type(wvl_data) :: wvl !wavelets ojects (unused but
123 : !required by various routines)
124 : character(len=fnlen) :: fname_current!Name of the TDCURRENT file
125 : character(len=fnlen) :: fname_tdener!Name of the TDENER file
126 : character(len=fnlen) :: fname_tdef !Name of the TDEFIELD file
127 : character(len=fnlen) :: fname_wfk0 !Name of the input WFK file containing
128 : !the intial (t=0) wfs
129 : !arrays
130 : integer,allocatable :: atindx(:) !index table of atom ordered by type
131 : integer,allocatable :: atindx1(:) !nb of the atom for each index in atindx
132 : integer,allocatable :: dimcprj(:) !Contains dimension for cprj array
133 : integer,allocatable :: dimcprj_srt(:) !Contains dimension for cprj array ordered by atom type
134 : integer,allocatable :: indsym(:,:,:) !atom indexing for symmetries
135 : integer,allocatable :: irrzon(:,:,:) !irreducible Brillouin zone
136 : integer,allocatable :: kg(:,:) !red. coord. of G vecs
137 : integer,allocatable :: nattyp(:) !nb of atoms of different types
138 : integer,allocatable :: npwarr(:) !number of PW at each k-point
139 : integer,allocatable :: symrec(:,:,:) !sym. operations in recip space
140 : real(dp) :: gprimd(3,3) !primitive cell vectors in recip space
141 : real(dp) :: gmet(3,3) !metric tensor in recip space
142 : real(dp) :: rprimd(3,3) !prim cell vectors in direct space
143 : real(dp) :: rmet(3,3) !metric tensor in direct space
144 : real(dp),allocatable :: cg(:,:) !WF coefficients in PW basis <k+G|psi_nk>
145 : real(dp),allocatable :: cg0(:,:) !Initial WF coefficients in PW basis <k+G|psi_nk>
146 : real(dp),allocatable :: current(:,:)!Current density
147 : real(dp),allocatable :: eigen(:) !eigen-energies
148 : real(dp),allocatable :: eigen0(:) !Initial eigen-energies (at t=0)
149 : real(dp),allocatable :: grvdw(:,:) !Gradient of the total energy coming
150 : !from VDW dispersion correction !FB: Needed?
151 : real(dp),allocatable :: nhat(:,:) !compensation charge density
152 : real(dp),allocatable :: nhatgr(:,:,:) !gradient of nhat
153 : real(dp),allocatable :: occ(:) !occupation numbers
154 : real(dp),allocatable :: occ0(:) !Initial occupation numbers
155 : real(dp),allocatable :: phnons(:,:,:) !For symmetries (nonsymmorphic translation phases)
156 : real(dp),allocatable :: ph1d(:,:) !Structure factor phase: exp(2Pi i G.xred)
157 : !on coarse grid
158 : real(dp),allocatable :: ph1df(:,:) !Structure factor phase: exp(2Pi i G.xred) for G
159 : !on fine grid
160 : real(dp),allocatable :: rhog(:,:) !charge density in recip space
161 : real(dp),allocatable :: rhor(:,:) !charge density in direct space
162 : real(dp),allocatable :: taug(:,:) !kin ener density in recip space !FB: Needed?
163 : real(dp),allocatable :: taur(:,:) !kin ener density in direct space !FB: Needed?
164 : real(dp),allocatable :: vhartr(:) !Hartree part of the potential
165 : real(dp),allocatable :: vpsp(:) !PSP part of the potential
166 : real(dp),allocatable :: vtrial(:,:) !"Trial" potential
167 : real(dp),allocatable :: vxc(:,:) !XC part of the potential
168 : real(dp),allocatable :: vxc_hybcomp(:,:) !Hybrid part of the xc potential !FB:Needed?
169 : real(dp),allocatable :: vxctau(:,:,:) !dV_{XC}/dtau (tau = kin. ener density)
170 : !for mGGAs !FB: Needed?
171 : real(dp),allocatable :: xred(:,:,:) !red. coord. of atoms
172 : real(dp),allocatable :: xccc3d(:) !3D core electron density
173 : !for XC core correction
174 : real(dp),allocatable :: xcctau3d(:) !3D core electron kin ener density
175 : !for XC core correction
176 : real(dp),allocatable :: ylm(:,:) !real spherical harmonics for each k+G
177 : real(dp),allocatable :: ylmgr(:,:,:)!real spherical harmonics gradients !FB: Needed?
178 : type(pawcprj_type), allocatable :: cprj(:,:) !projectors applied on WF <p_lmn|C_nk>
179 : type(pawcprj_type), allocatable :: cprj0(:,:) !projectors applied on WF <p_lmn|C_nk>
180 : type(paw_an_type), allocatable :: paw_an(:) !various arrays on angular mesh
181 : type(pawfgrtab_type),allocatable :: pawfgrtab(:) !PAW atomic data on fine grid
182 : type(paw_ij_type), allocatable :: paw_ij(:) !various arrays on partial waves (i,j channels)
183 : type(pawrad_type), pointer :: pawrad(:) => NULL() !radial grid in PAW sphere
184 : type(pawrhoij_type), pointer :: pawrhoij(:) => NULL() !operator rho_ij= <psi|p_i><p_j|psi>
185 : type(pawtab_type), pointer :: pawtab(:) => NULL() !tabulated PAW atomic data
186 :
187 : contains
188 :
189 : procedure :: init => tdks_init
190 : procedure :: free => tdks_free
191 :
192 : end type tdks_type
193 : !!***
194 :
195 : contains
196 : !!***
197 :
198 : !!****f* m_rttddft_tdks/tdks_init
199 : !!
200 : !! NAME
201 : !! tdks_init
202 : !!
203 : !! FUNCTION
204 : !! Initialize the tdks object
205 : !!
206 : !! INPUTS
207 : !! codvsn = code version
208 : !! dtfil <type datafiles_type> = infos about file names, file unit numbers
209 : !! dtset <type(dataset_type)> = all input variables for this dataset
210 : !! mpi_enreg <MPI_type> = MPI-parallelisation information
211 : !! pawang <type(pawang_type)> = paw angular mesh and related data
212 : !! pawrad(ntypat*usepaw) <type(pawrad_type)> = paw radial mesh and related data
213 : !! pawtab(ntypat*usepaw) <type(pawtab_type)> = paw tabulated starting data
214 : !! psps <type(pseudopotential_type)> = variables related to pseudopotentials
215 : !!
216 : !! OUTPUT
217 : !! tdks <class(tdks_type)> = the tdks object to initialize
218 : !!
219 : !! SOURCE
220 50 : subroutine tdks_init(tdks ,codvsn, dtfil, dtset, mpi_enreg, pawang, pawrad, pawtab, psps)
221 :
222 : !Arguments ------------------------------------
223 : !scalars
224 : class(tdks_type), intent(inout) :: tdks
225 : character(len=8), intent(in) :: codvsn
226 : type(datafiles_type), intent(in) :: dtfil
227 : type(dataset_type), intent(inout) :: dtset
228 : type(MPI_type), intent(inout) :: mpi_enreg
229 : type(pawang_type), intent(inout),target :: pawang
230 : type(pseudopotential_type), intent(inout) :: psps
231 : !arrays
232 : type(pawrad_type), intent(inout),target :: pawrad(psps%ntypat*psps%usepaw)
233 : type(pawtab_type), intent(inout),target :: pawtab(psps%ntypat*psps%usepaw)
234 :
235 : !Local variables-------------------------------
236 : !scalars
237 : integer :: ierr
238 : integer :: my_natom
239 : integer :: ncpgr
240 : integer :: psp_gencond
241 : real(dp) :: ecut_eff
242 : character(len=500) :: msg
243 : character(len=fnlen) :: fname_wfk
244 : type(extfpmd_type),pointer :: extfpmd => null()
245 : !arrays
246 50 : real(dp),allocatable :: doccde(:)
247 :
248 : ! ***********************************************************************
249 :
250 50 : my_natom=mpi_enreg%my_natom
251 :
252 : !1) Various initializations & checks (MPI, PW, FFT, PSP, Symmetry ...)
253 50 : call first_setup(codvsn,dtfil,dtset,ecut_eff,mpi_enreg,pawrad,pawtab,psps,psp_gencond,tdks)
254 :
255 : !2) Deals with restart and setup some basic variables and filenames
256 50 : tdks%first_step = 1
257 50 : tdks%fname_tdener = dtfil%fnameabo_td_ener
258 50 : tdks%fname_wfk0 = dtfil%fnamewffk
259 50 : fname_wfk = dtfil%fnamewffk
260 50 : tdks%fname_tdef = dtfil%fnameabo_td_ef
261 50 : tdks%fname_current = dtfil%fnameabo_td_current
262 50 : if (dtset%td_restart > 0) then
263 20 : if (mpi_enreg%me == 0) then
264 6 : if (open_file('TD_RESTART', msg, newunit=tdks%tdrestart_unit, status='old', form='formatted') /= 0) then
265 0 : write(msg,'(a,a,a)') 'Error while trying to open file TD_RESTART needed to restart the calculation.'
266 0 : ABI_ERROR(msg)
267 : end if
268 6 : read(tdks%tdrestart_unit,*) tdks%first_step
269 6 : tdks%first_step = tdks%first_step + 1
270 6 : read(tdks%tdrestart_unit,*) tdks%fname_wfk0
271 6 : read(tdks%tdrestart_unit,*) fname_wfk
272 6 : read(tdks%tdrestart_unit,*) tdks%fname_tdener
273 6 : if (dtset%td_ef_type /= 0) then
274 3 : read(tdks%tdrestart_unit,*) tdks%fname_tdef
275 : end if
276 6 : if (dtset%prtcurrent /= 0) then
277 3 : read(tdks%tdrestart_unit,*) tdks%fname_current
278 : end if
279 : end if
280 : !Send to all procs
281 20 : call xmpi_bcast(tdks%first_step,0,mpi_enreg%comm_world,ierr)
282 20 : call xmpi_bcast(tdks%fname_tdener,0,mpi_enreg%comm_world,ierr)
283 20 : call xmpi_bcast(tdks%fname_wfk0,0,mpi_enreg%comm_world,ierr)
284 20 : call xmpi_bcast(fname_wfk,0,mpi_enreg%comm_world,ierr)
285 20 : call xmpi_bcast(tdks%fname_tdef,0,mpi_enreg%comm_world,ierr)
286 20 : call xmpi_bcast(tdks%fname_current,0,mpi_enreg%comm_world,ierr)
287 : else
288 30 : if (mpi_enreg%me == 0) then
289 16 : if (open_file('TD_RESTART', msg, newunit=tdks%tdrestart_unit, status='replace', form='formatted') /= 0) then
290 0 : write(msg,'(a,a,a)') 'Error while trying to open file TD_RESTART.'
291 0 : ABI_ERROR(msg)
292 : end if
293 : end if
294 : end if
295 :
296 : !3) Reads initial KS orbitals from file (calls inwffil)
297 50 : call read_wfk(dtfil,dtset,ecut_eff,fname_wfk,mpi_enreg,tdks)
298 :
299 : !4) Init occupation numbers
300 150 : ABI_MALLOC(tdks%occ0,(dtset%mband*dtset%nkpt*dtset%nsppol))
301 9138 : tdks%occ0(:)=dtset%occ_orig(:,1)
302 : !calc occupation number with metallic occupation using the previously read WF
303 50 : if (dtset%occopt>=3.and.dtset%occopt<=9) then ! allowing for occopt 9
304 48 : ABI_MALLOC(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol))
305 : call newocc(doccde,tdks%eigen0,tdks%energies%entropy_ks,tdks%energies%e_fermie, &
306 : & tdks%energies%e_fermih,dtset%ivalence,dtset%spinmagntarget, &
307 : & dtset%mband,dtset%nband,dtset%nelect,dtset%ne_qFD,dtset%nh_qFD, &
308 : & dtset%nkpt,dtset%nspinor,dtset%nsppol,tdks%occ0,dtset%occopt, &
309 24 : & dtset%prtvol,dtset%tphysel,dtset%tsmear,dtset%wtk,extfpmd=extfpmd)
310 24 : ABI_FREE(doccde)
311 : end if
312 :
313 : !5) Some further initialization (Mainly for PAW and allocation of arrays for Hamiltonian and densities)
314 50 : call second_setup(dtset,mpi_enreg,pawang,pawrad,pawtab,psps,psp_gencond,tdks)
315 :
316 : !6) TD external elec. field perturbation
317 50 : if (dtset%td_ef_type/=0 .and. psps%useylm/=1) ABI_ERROR("TD Electric field only works with spherical harmonics (useylm=1)")
318 50 : if (dtset%td_ef_type/=0 .and. psps%usepaw/=1) ABI_ERROR("TD Electric field only works with PAW")
319 : !Init vector potential and associated constants
320 50 : if (dtset%td_ef_type/=0 .or. dtset%prtcurrent/=0) then
321 63 : ABI_MALLOC(tdks%current,(3,dtset%nsppol))
322 134 : tdks%current = zero
323 : end if
324 : call tdks%tdef%init(dtset%td_ef_type,dtset%td_ef_pol,dtset%td_ef_ezero,dtset%td_ef_tzero, &
325 50 : & dtset%td_ef_lambda,dtset%td_ef_tau,dtset%td_ef_induced_vecpot,dtset%nkpt,dtset%kptns)
326 50 : if (dtset%td_restart /=0) then
327 20 : call tdks%tdef%restart(mpi_enreg, tdks%tdrestart_unit)
328 : call tdks%tdef%update(dtset,mpi_enreg,(tdks%first_step-1)*dtset%dtele,tdks%rprimd,tdks%gprimd,tdks%kg, &
329 20 : & psps%mpsang,tdks%npwarr,tdks%ylm,tdks%ylmgr,tdks%current,update_vecpot_ind=.false.)
330 : else
331 : call tdks%tdef%update(dtset,mpi_enreg,(tdks%first_step-1)*dtset%dtele,tdks%rprimd,tdks%gprimd,tdks%kg, &
332 30 : & psps%mpsang,tdks%npwarr,tdks%ylm,tdks%ylmgr,tdks%current)
333 : end if
334 :
335 : !7) Keep initial cg and cproj in memory for occupations
336 : !Keep initial wavefunction in memory
337 50 : if (dtset%td_restart == 0) then
338 90 : ABI_MALLOC(tdks%cg0,(2,tdks%mcg))
339 3468894 : tdks%cg0(:,:) = tdks%cg(:,:)
340 : end if
341 : !and associated cprojs to compute occupations
342 50 : if (psps%usepaw ==1) then
343 30 : ncpgr=0
344 12216 : ABI_MALLOC(tdks%cprj0,(dtset%natom,tdks%mcprj))
345 30 : call pawcprj_alloc(tdks%cprj0,ncpgr,tdks%dimcprj)
346 : call ctocprj(tdks%atindx,tdks%cg0,1,tdks%cprj0,tdks%gmet,tdks%gprimd,0,0,0, &
347 : & dtset%istwfk,tdks%kg,tdks%tdef%kpa,tdks%mcg,tdks%mcprj,dtset%mgfft, &
348 : & dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,dtset%natom,tdks%nattyp, &
349 : & dtset%nband,dtset%natom,dtset%ngfft,dtset%nkpt,dtset%nloalg, &
350 : & tdks%npwarr,dtset%nspinor,dtset%nsppol,dtset%nsppol,psps%ntypat, &
351 : & dtset%paral_kgb,tdks%ph1d,psps,tdks%rmet,dtset%typat,tdks%ucvol, &
352 30 : & tdks%unpaw,tdks%xred,tdks%ylm,tdks%ylmgr)
353 : end if
354 150 : ABI_MALLOC(tdks%occ,(dtset%mband*dtset%nkpt*dtset%nsppol))
355 :
356 : !Keep some additional stuff in memory within the tdks object
357 50 : tdks%unpaw = dtfil%unpaw
358 50 : tdks%dt = dtset%dtele
359 50 : tdks%ntime = dtset%ntime
360 :
361 50 : tdks%pawang => pawang
362 50 : tdks%pawrad => pawrad
363 50 : tdks%pawtab => pawtab
364 :
365 50 : end subroutine tdks_init
366 : !!***
367 :
368 : !!****f* m_rttddft_tdks/tdks_free
369 : !!
370 : !! NAME
371 : !! tdks_free
372 : !!
373 : !! FUNCTION
374 : !! Free all the memory associated with the tdks object
375 : !!
376 : !! INPUTS
377 : !! tdks <class(tdks_type)> = the tdks object to free
378 : !! dtset <type(dataset_type)> = all input variables for this dataset
379 : !! mpi_enreg <MPI_type> = MPI-parallelisation information
380 : !! psps <type(pseudopotential_type)> = variables related to pseudopotentials
381 : !!
382 : !! OUTPUT
383 : !!
384 : !! SOURCE
385 49 : subroutine tdks_free(tdks,dtset,mpi_enreg,psps)
386 :
387 : !Arguments ------------------------------------
388 : !scalars
389 : class(tdks_type), intent(inout) :: tdks
390 : type(dataset_type), intent(inout) :: dtset
391 : type(MPI_type), intent(inout) :: mpi_enreg
392 : type(pseudopotential_type), intent(inout) :: psps
393 :
394 : ! ***********************************************************************
395 :
396 : !Destroy hidden save variables
397 49 : call bandfft_kpt_destroy_array(bandfft_kpt,mpi_enreg)
398 49 : if (psps%usepaw ==1) then
399 29 : call destroy_invovl(dtset%nkpt,dtset%gpu_option)
400 : end if
401 49 : if(tdks%gemm_nonlop_use_gemm .and. dtset%gpu_option==ABI_GPU_DISABLED) then
402 0 : call destroy_gemm_nonlop(dtset%gpu_option)
403 : end if
404 :
405 : !Call type destructors
406 49 : call destroy_sc_dmft(tdks%paw_dmft)
407 49 : call pawfgr_destroy(tdks%pawfgr)
408 49 : call tdks%hdr%free()
409 :
410 : !Nullify pointers
411 49 : if(associated(tdks%pawang)) tdks%pawang => null()
412 49 : if(associated(tdks%pawrad)) tdks%pawrad => null()
413 49 : if(associated(tdks%pawtab)) tdks%pawtab => null()
414 :
415 : !Deallocate allocatables
416 49 : ABI_SFREE(tdks%atindx)
417 49 : ABI_SFREE(tdks%atindx1)
418 49 : ABI_SFREE(tdks%cg)
419 49 : ABI_SFREE(tdks%cg0)
420 49 : ABI_SFREE(tdks%current)
421 49 : ABI_SFREE(tdks%dimcprj)
422 49 : ABI_SFREE(tdks%dimcprj_srt)
423 49 : ABI_SFREE(tdks%eigen)
424 49 : ABI_SFREE(tdks%eigen0)
425 49 : ABI_SFREE(tdks%grvdw)
426 49 : ABI_SFREE(tdks%indsym)
427 49 : ABI_SFREE(tdks%irrzon)
428 49 : ABI_SFREE(tdks%kg)
429 49 : ABI_SFREE(tdks%nattyp)
430 49 : ABI_SFREE(tdks%nhat)
431 49 : ABI_SFREE(tdks%nhatgr)
432 49 : ABI_SFREE(tdks%npwarr)
433 49 : ABI_SFREE(tdks%occ)
434 49 : ABI_SFREE(tdks%occ0)
435 49 : ABI_SFREE(tdks%ph1d)
436 49 : ABI_SFREE(tdks%ph1df)
437 49 : ABI_SFREE(tdks%phnons)
438 49 : ABI_SFREE(tdks%rhog)
439 49 : ABI_SFREE(tdks%rhor)
440 49 : ABI_SFREE(tdks%symrec)
441 49 : ABI_SFREE(tdks%taug)
442 49 : ABI_SFREE(tdks%taur)
443 49 : ABI_SFREE(tdks%vhartr)
444 49 : ABI_SFREE(tdks%vpsp)
445 49 : ABI_SFREE(tdks%vtrial)
446 49 : ABI_SFREE(tdks%vxc)
447 49 : ABI_SFREE(tdks%vxctau)
448 49 : ABI_SFREE(tdks%vxc_hybcomp)
449 49 : ABI_SFREE(tdks%xred)
450 49 : ABI_SFREE(tdks%xccc3d)
451 49 : ABI_SFREE(tdks%xcctau3d)
452 49 : ABI_SFREE(tdks%ylm)
453 49 : ABI_SFREE(tdks%ylmgr)
454 49 : ABI_SFREE(tdks%tdef%kpa)
455 :
456 49 : if(allocated(tdks%cprj)) then
457 29 : call pawcprj_free(tdks%cprj)
458 8093 : ABI_FREE(tdks%cprj)
459 : end if
460 49 : if(allocated(tdks%cprj0)) then
461 29 : call pawcprj_free(tdks%cprj0)
462 8093 : ABI_FREE(tdks%cprj0)
463 : end if
464 49 : if(allocated(tdks%paw_an)) then
465 29 : call paw_an_free(tdks%paw_an)
466 95 : ABI_FREE(tdks%paw_an)
467 : end if
468 49 : if(allocated(tdks%pawfgrtab)) then
469 29 : call pawfgrtab_free(tdks%pawfgrtab)
470 95 : ABI_FREE(tdks%pawfgrtab)
471 : end if
472 49 : if(allocated(tdks%paw_ij)) then
473 29 : call paw_ij_free(tdks%paw_ij)
474 95 : ABI_FREE(tdks%paw_ij)
475 : end if
476 49 : if(associated(tdks%pawrhoij)) then
477 49 : call pawrhoij_free(tdks%pawrhoij)
478 49 : ABI_FREE(tdks%pawrhoij)
479 : end if
480 :
481 49 : end subroutine tdks_free
482 : !!***
483 :
484 : !!****f* m_rttddft_tdks/first_setup
485 : !!
486 : !! NAME
487 : !! first_setup
488 : !!
489 : !! FUNCTION
490 : !! Intialize many important quantities before running RT-TDDFT
491 : !! (PW, FFT, PSP, Symmetry etc.)
492 : !!
493 : !! INPUTS
494 : !! codvsn = code version
495 : !! dtfil <type datafiles_type> = infos about file names, file unit numbers
496 : !! dtset <type(dataset_type)> = all input variables for this dataset
497 : !! mpi_enreg <MPI_type> = MPI-parallelisation information
498 : !! pawrad(ntypat*usepaw) <type(pawrad_type)> = paw radial mesh and related data
499 : !! pawtab(ntypat*usepaw) <type(pawtab_type)> = paw tabulated starting data
500 : !! psps <type(pseudopotential_type)> = variables related to pseudopotentials
501 : !! tdks <type(tdks_type)> = the tdks object to initialize
502 : !!
503 : !! OUTPUT
504 : !! psp_gencond <integer> = store conditions for generating psp
505 : !! ecut_eff <real(dp)> = effective PW cutoff energy
506 : !!
507 : !! NOTES
508 : !! USE OF FFT GRIDS:
509 : !! =================
510 : !! In case of PAW:
511 : !! ---------------
512 : !! Two FFT grids are used:
513 : !! - A "coarse" FFT grid (defined by ecut)
514 : !! for the application of the Hamiltonian on the plane waves basis.
515 : !! It is defined by nfft, ngfft, mgfft, ...
516 : !! Hamiltonian, wave-functions, density related to WFs (rhor here), ...
517 : !! are expressed on this grid.
518 : !! - A "fine" FFT grid (defined) by ecutdg)
519 : !! for the computation of the density inside PAW spheres.
520 : !! It is defined by nfftf, ngfftf, mgfftf, ...
521 : !! Total density, potentials, ...
522 : !! are expressed on this grid.
523 : !! In case of norm-conserving:
524 : !! ---------------------------
525 : !! - Only the usual FFT grid (defined by ecut) is used.
526 : !! It is defined by nfft, ngfft, mgfft, ...
527 : !! For compatibility reasons, (nfftf,ngfftf,mgfftf)
528 : !! are set equal to (nfft,ngfft,mgfft) in that case.
529 : !! In case of wavelets:
530 : !! --------------------
531 : !! - Only the usual FFT grid (defined by wvl_crmult) is used.
532 : !! It is defined by nfft, ngfft, mgfft, ... This is strictly not
533 : !! an FFT grid since its dimensions are not suited for FFTs. They are
534 : !! defined by wvl_setngfft().
535 : !! For compatibility reasons, (nfftf,ngfftf,mgfftf)
536 : !! are set equal to (nfft,ngfft,mgfft) in that case.
537 : !!
538 : !! SOURCE
539 50 : subroutine first_setup(codvsn,dtfil,dtset,ecut_eff,mpi_enreg,pawrad,pawtab,psps,psp_gencond,tdks)
540 :
541 : !Arguments ------------------------------------
542 : !scalars
543 : character(len=8), intent(in) :: codvsn
544 : integer, intent(out) :: psp_gencond
545 : real(dp), intent(out) :: ecut_eff
546 : type(datafiles_type), intent(in) :: dtfil
547 : type(dataset_type), intent(inout) :: dtset
548 : type(pseudopotential_type), intent(inout) :: psps
549 : type(MPI_type), intent(inout) :: mpi_enreg
550 : type(tdks_type), intent(inout) :: tdks
551 : !arrays
552 : type(pawrad_type), intent(inout) :: pawrad(psps%ntypat*psps%usepaw)
553 : type(pawtab_type), intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
554 :
555 : !Local variables-------------------------------
556 : !scalars
557 : integer,parameter :: response=0, cplex=1
558 : integer :: comm_psp
559 : integer :: gscase
560 : integer :: iatom, ierr, itypat, indx
561 : integer :: mgfftf, my_natom
562 : integer :: npwmin, nfftot
563 : real(dp) :: gsqcut_eff, gsqcutc_eff
564 : real(dp) :: ecutdg_eff
565 50 : type(ebands_t) :: bstruct
566 : !arrays
567 : character(len=500) :: msg
568 50 : integer, allocatable :: npwarr_(:)
569 : integer :: ngfft(18)
570 : integer :: ngfftf(18)
571 100 : integer :: npwtot(dtset%nkpt)
572 :
573 : ! ***********************************************************************
574 :
575 50 : my_natom=mpi_enreg%my_natom
576 :
577 : !** Init FFT grid(s) sizes (be careful !)
578 : !See NOTES in the comments at the beginning of this subroutine.
579 50 : tdks%nfft = dtset%nfft
580 : call pawfgr_init(tdks%pawfgr,dtset,mgfftf,tdks%nfftf,ecut_eff,ecutdg_eff, &
581 50 : & ngfft,ngfftf)
582 :
583 : !** Init to zero different energies
584 50 : call tdks%energies%init()
585 50 : tdks%ecore = zero
586 50 : tdks%etot = zero
587 :
588 : !** various additional setup mostly related to fft grids and the box (rprimd, metric..)
589 : call setup1(dtset%acell_orig,tdks%bantot,dtset,ecutdg_eff,ecut_eff,tdks%gmet, &
590 : & tdks%gprimd,gsqcut_eff,gsqcutc_eff,ngfftf,ngfft,dtset%nkpt, &
591 : & dtset%nsppol,response,tdks%rmet,dtset%rprim_orig,tdks%rprimd, &
592 50 : & tdks%ucvol,psps%usepaw)
593 :
594 : !** Set up the basis sphere of planewaves
595 150 : ABI_MALLOC(tdks%npwarr,(dtset%nkpt))
596 150 : ABI_MALLOC(tdks%kg,(3,dtset%mpw*dtset%mkmem))
597 : call kpgio(ecut_eff,dtset%exchn2n3d,tdks%gmet,dtset%istwfk,tdks%kg,dtset%kptns, &
598 : & dtset%mkmem,dtset%nband,dtset%nkpt,'PERS',mpi_enreg,dtset%mpw, &
599 50 : & tdks%npwarr,npwtot,dtset%nsppol)
600 : call bandfft_kpt_init1(bandfft_kpt,dtset%istwfk,tdks%kg,dtset%mgfft,dtset%mkmem, &
601 : & mpi_enreg,dtset%mpw,dtset%nband,dtset%nkpt,tdks%npwarr, &
602 50 : & dtset%nsppol)
603 :
604 : !** Use efficient BLAS calls for computing the non local potential (No GPU yet)
605 50 : if(dtset%use_gemm_nonlop == 1 .and. dtset%gpu_option==ABI_GPU_DISABLED) then
606 : ! set global variable
607 0 : tdks%gemm_nonlop_use_gemm = .true.
608 0 : call init_gemm_nonlop(dtset%gpu_option)
609 : else
610 50 : tdks%gemm_nonlop_use_gemm = .false.
611 : end if
612 :
613 : !** TODO: uncomment when gemm_nonlop can be used on GPU
614 : ! if(dtset%use_gemm_nonlop == 1 .and. dtset%gpu_option/=ABI_GPU_DISABLED) then
615 : ! ! set global variable
616 : ! tdks%gemm_nonlop_use_gemm_gpu = .true.
617 : ! !call init_gemm_nonlop_gpu(dtset%nkpt)
618 : ! else
619 : ! tdks%gemm_nonlop_use_gemm_gpu = .false.
620 : ! end if
621 :
622 : !** Initialize band structure datatype
623 150 : ABI_MALLOC(npwarr_,(dtset%nkpt))
624 958 : npwarr_(:)=tdks%npwarr(:)
625 50 : if (dtset%paral_kgb/=0) then
626 32 : call xmpi_sum(npwarr_,mpi_enreg%comm_bandfft,ierr)
627 : end if
628 50 : call bstruct%from_dtset(dtset, npwarr_)
629 50 : ABI_FREE(npwarr_)
630 50 : call unpack_eneocc(dtset%nkpt,dtset%nsppol,bstruct%mband,bstruct%nband,dtset%occ_orig(:,1),bstruct%occ,val=zero)
631 :
632 : !** Open and read pseudopotential files
633 50 : comm_psp=mpi_enreg%comm_cell
634 : call pspini(dtset,dtfil,tdks%ecore,psp_gencond,gsqcutc_eff,gsqcut_eff,pawrad, &
635 50 : & pawtab,psps,tdks%rprimd,comm_mpi=comm_psp)
636 :
637 : !In case of isolated computations, ecore must be set to zero
638 : !because its contribution is counted in the ewald energy as the ion-ion interaction.
639 50 : if (dtset%icoulomb == 1) tdks%ecore = zero
640 :
641 : !Include core energy?
642 100 : select case(dtset%usepotzero)
643 : case(0,1)
644 50 : tdks%energies%e_corepsp = tdks%ecore / tdks%ucvol
645 50 : tdks%energies%e_corepspdc = zero
646 : case(2)
647 : ! No need to include the PspCore energy since it is already included in the
648 : ! local pseudopotential (vpsp)
649 0 : tdks%energies%e_corepsp = zero
650 50 : tdks%energies%e_corepspdc = zero
651 : end select
652 :
653 : !** Initialize PAW atomic occupancies
654 218 : ABI_MALLOC(tdks%pawrhoij,(my_natom*psps%usepaw))
655 50 : if (psps%usepaw == 1) then
656 : call initrhoij(dtset%pawcpxocc,dtset%lexexch,dtset%lpawu,my_natom,dtset%natom, &
657 : & dtset%nspden,dtset%nspinor,dtset%nsppol,dtset%ntypat, &
658 : & tdks%pawrhoij,dtset%pawspnorb,pawtab,cplex,dtset%spinat, &
659 : & dtset%typat,comm_atom=mpi_enreg%comm_atom, &
660 30 : & mpi_atmtab=mpi_enreg%my_atmtab)
661 : end if
662 :
663 : !Nullify wvl_data. It is important to do so irregardless of the value of usewvl
664 : !Only needed here because hdr%init requires a wvl object for the wvl%descr input
665 50 : call nullify_wvl_data(tdks%wvl)
666 :
667 : !** Initialize header
668 50 : gscase=0
669 : call tdks%hdr%init(bstruct,codvsn,dtset,pawtab,gscase,psps,tdks%wvl%descr,&
670 50 : comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
671 :
672 : !Clean band structure datatype
673 50 : call bstruct%free()
674 :
675 : !** PW basis set: test if the problem is ill-defined.
676 958 : npwmin=minval(tdks%hdr%npwarr(:))
677 50 : if (dtset%mband > npwmin) then
678 : ! No way we can solve the problem. Abort now!
679 0 : write(msg,"(2(a,i0),4a)") "Number of bands nband= ",dtset%mband, &
680 0 : & " > number of planewaves npw= ",npwmin,ch10, &
681 0 : & "The number of eigenvectors cannot be greater that the size of the Hamiltonian!",&
682 0 : & ch10, "Action: decrease nband or, alternatively, increase ecut"
683 0 : if (dtset%ionmov/=23) then
684 0 : ABI_ERROR(msg)
685 : else
686 0 : ABI_WARNING(msg)
687 : end if
688 :
689 50 : else if (dtset%mband >= 0.9 * npwmin) then
690 : ! Warn the user
691 0 : write(msg,"(a,i0,a,f6.1,4a)") "Number of bands nband= ",dtset%mband, &
692 0 : & " >= 0.9 * maximum number of planewaves= ",0.9*npwmin,ch10,&
693 0 : & "This could lead to some instabilities, you might want to decrease nband or increase ecut!", &
694 0 : & ch10,"Assume experienced user. Execution will continue."
695 0 : ABI_WARNING(msg)
696 : end if
697 :
698 : !** Initialize symmetry
699 50 : nfftot=ngfft(1)*ngfft(2)*ngfft(3)
700 200 : ABI_MALLOC(tdks%irrzon,(nfftot**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
701 200 : ABI_MALLOC(tdks%phnons,(2,nfftot**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
702 200 : ABI_MALLOC(tdks%indsym,(4,dtset%nsym,dtset%natom))
703 150 : ABI_MALLOC(tdks%symrec,(3,3,dtset%nsym))
704 1360680 : tdks%irrzon(:,:,:)=0
705 2040820 : tdks%phnons(:,:,:)=zero
706 104838 : tdks%indsym(:,:,:)=0
707 76178 : tdks%symrec(:,:,:)=0
708 :
709 : !Do symmetry stuff if nsym>1
710 50 : if (dtset%nsym>1) then
711 : call setsym(tdks%indsym,tdks%irrzon,dtset%iscf,dtset%natom, &
712 : & nfftot,ngfft,dtset%nspden,dtset%nsppol,dtset%nsym, &
713 : & tdks%phnons,dtset%symafm,tdks%symrec,dtset%symrel, &
714 50 : & dtset%tnons,dtset%typat,dtset%xred_orig)
715 :
716 : !Make sure dtset%iatfix does not break symmetry
717 50 : call fixsym(dtset%iatfix,tdks%indsym,dtset%natom,dtset%nsym)
718 : else
719 : !The symrec array is used by initberry even in case nsym = 1
720 0 : tdks%symrec(:,:,1) = 0
721 0 : tdks%symrec(1,1,1) = 1 ; tdks%symrec(2,2,1) = 1 ; tdks%symrec(3,3,1) = 1
722 : end if
723 :
724 : !** Initialize and eventually symmetrize reduced atomic coordinates
725 200 : ABI_MALLOC(tdks%xred,(3,dtset%natom,dtset%nimage))
726 742 : tdks%xred = dtset%xred_orig
727 : !Eventually symmetrize atomic coordinates over space group elements
728 : call symmetrize_xred(dtset%natom,dtset%nsym,dtset%symrel,dtset%tnons,tdks%xred, &
729 50 : & indsym=tdks%indsym)
730 :
731 : !** Create the atindx array
732 : !** index table of atoms, in order for them to be used type after type.
733 150 : ABI_MALLOC(tdks%atindx,(dtset%natom))
734 100 : ABI_MALLOC(tdks%atindx1,(dtset%natom))
735 150 : ABI_MALLOC(tdks%nattyp,(psps%ntypat))
736 50 : indx=1
737 100 : do itypat=1,psps%ntypat
738 50 : tdks%nattyp(itypat)=0
739 248 : do iatom=1,dtset%natom
740 198 : if(dtset%typat(iatom)==itypat)then
741 148 : tdks%atindx(iatom)=indx
742 148 : tdks%atindx1(indx)=iatom
743 148 : indx=indx+1
744 148 : tdks%nattyp(itypat)=tdks%nattyp(itypat)+1
745 : end if
746 : end do
747 : end do
748 :
749 : !** Calculate zion: the total positive charge acting on the valence electrons
750 50 : tdks%zion=zero
751 198 : do iatom=1,dtset%natom
752 198 : tdks%zion=tdks%zion+psps%ziontypat(dtset%typat(iatom))
753 : end do
754 :
755 100 : end subroutine first_setup
756 : !!***
757 :
758 : !!****f* m_rttddft_tdks/second_setup
759 : !!
760 : !! NAME
761 : !! second_setup
762 : !!
763 : !! FUNCTION
764 : !! Further important initialization required after reading WFK and computing
765 : !! occupation numbers in paticular related to PAW
766 : !!
767 : !! INPUTS
768 : !! dtset <type(dataset_type)> = all input variables for this dataset
769 : !! mpi_enreg <MPI_type> = MPI-parallelisation information
770 : !! pawang <type(pawang_type)> = paw angular mesh and related data
771 : !! pawrad(ntypat*usepaw) <type(pawrad_type)> = paw radial mesh and related data
772 : !! pawtab(ntypat*usepaw) <type(pawtab_type)> = paw tabulated starting data
773 : !! psps <type(pseudopotential_type)> = variables related to pseudopotentials
774 : !! psp_gencond <integer> = store conditions for generating psp
775 : !! tdks <type(tdks_type)> = the tdks object to initialize
776 : !!
777 : !! OUTPUT
778 : !!
779 : !! SOURCE
780 50 : subroutine second_setup(dtset, mpi_enreg, pawang, pawrad, pawtab, psps, psp_gencond, tdks)
781 :
782 : !Arguments ------------------------------------
783 : !scalars
784 : integer, intent(in) :: psp_gencond
785 : type(pawang_type), intent(inout) :: pawang
786 : type(dataset_type), intent(inout) :: dtset
787 : type(pseudopotential_type), intent(inout) :: psps
788 : type(MPI_type), intent(inout) :: mpi_enreg
789 : !arrays
790 : type(pawrad_type), intent(inout) :: pawrad(psps%ntypat*psps%usepaw)
791 : type(pawtab_type), intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
792 : type(tdks_type), intent(inout) :: tdks
793 :
794 : !Local variables-------------------------------
795 : !scalars
796 : logical :: call_pawinit
797 : integer, parameter :: cplex = 1
798 : integer :: forces_needed
799 : integer :: gnt_option
800 : integer :: has_dijhat, has_vhartree, has_dijfock
801 : integer :: has_dijnd, has_dijU, has_vxctau
802 : integer :: my_natom, my_nspinor
803 : integer :: ncpgr
804 : integer :: optcut, optgr0, optgr1, optgr2, optrad
805 : integer :: stress_needed
806 : integer :: ylm_option
807 : integer :: use_hybcomp, usevxctau
808 : real(dp) :: boxcut
809 : real(dp) :: gsqcut_shp
810 : real(dp) :: hyb_range_fock
811 : real(dp),parameter :: k0(3)=(/zero,zero,zero/)
812 : !arrays
813 50 : integer,allocatable :: l_size_atm(:)
814 :
815 : ! ***********************************************************************
816 :
817 50 : my_natom=mpi_enreg%my_natom
818 50 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
819 :
820 : !FB: @MT needed?
821 50 : if (psps%usepaw==1) then
822 : call pawrhoij_copy(tdks%hdr%pawrhoij,tdks%pawrhoij,comm_atom=mpi_enreg%comm_atom, &
823 30 : & mpi_atmtab=mpi_enreg%my_atmtab)
824 : end if
825 :
826 : !FB: Needed because paw_dmft is required in mkrho
827 : !PAW related operations
828 : !Initialize paw_dmft, even if neither dmft not paw are used
829 50 : call init_sc_dmft(dtset,psps%mpsang,tdks%paw_dmft,mpi_enreg=mpi_enreg,use_sc_dmft=0)
830 :
831 :
832 : !*** Main PAW initialization ***
833 :
834 : !** Setup the Ylm for each k point
835 50 : if (psps%useylm==1) then
836 120 : ABI_MALLOC(tdks%ylm,(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm))
837 120 : ABI_MALLOC(tdks%ylmgr,(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm))
838 30 : ylm_option=0
839 : call initylmg(tdks%gprimd,tdks%kg,dtset%kptns,dtset%mkmem,mpi_enreg,&
840 : & psps%mpsang,dtset%mpw,dtset%nband,dtset%nkpt,&
841 30 : & tdks%npwarr,dtset%nsppol,ylm_option,tdks%rprimd,tdks%ylm,tdks%ylmgr)
842 : else
843 20 : ABI_MALLOC(tdks%ylm,(0,0))
844 20 : ABI_MALLOC(tdks%ylmgr,(0,0,0))
845 : end if
846 :
847 50 : tdks%mcprj=0;tdks%mband_cprj=0
848 50 : if(psps%usepaw==1) then
849 30 : gnt_option=1
850 30 : if (dtset%pawxcdev==2.or.(dtset%pawxcdev==1.and.dtset%positron/=0)) gnt_option=2
851 :
852 : !** Test if we have to call pawinit
853 : ! Some gen-cond have to be added...
854 30 : call paw_gencond(dtset,gnt_option,"test",call_pawinit)
855 :
856 30 : if (psp_gencond==1.or.call_pawinit) then
857 2 : gsqcut_shp=two*abs(dtset%diecut)*dtset%dilatmx**2/pi**2
858 2 : hyb_range_fock=zero
859 2 : if (dtset%ixc<0) then
860 2 : call libxc_functionals_get_hybridparams(hyb_range=hyb_range_fock)
861 : end if
862 : call pawinit(dtset%effmass_free,gnt_option,gsqcut_shp,hyb_range_fock, &
863 : & dtset%pawlcutd,dtset%pawlmix,psps%mpsang,dtset%pawnphi, &
864 : & dtset%nsym,dtset%pawntheta,pawang,pawrad,dtset%pawspnorb, &
865 2 : & pawtab,dtset%pawxcdev,dtset%ixc,dtset%usepotzero)
866 :
867 : ! Update internal values
868 2 : call paw_gencond(dtset,gnt_option,"save",call_pawinit)
869 : end if
870 60 : psps%n1xccc=maxval(pawtab(1:psps%ntypat)%usetcore)
871 : call setsym_ylm(tdks%gprimd,pawang%l_max-1,dtset%nsym,dtset%pawprtvol, &
872 30 : & tdks%rprimd,tdks%symrec,pawang%zarot)
873 :
874 : !** Initialisation of cprj
875 30 : tdks%mband_cprj=dtset%mband
876 30 : if (dtset%paral_kgb/=0) tdks%mband_cprj=tdks%mband_cprj/mpi_enreg%nproc_band
877 30 : tdks%mcprj=my_nspinor*tdks%mband_cprj*dtset%mkmem*dtset%nsppol
878 12216 : ABI_MALLOC(tdks%cprj,(dtset%natom,tdks%mcprj))
879 30 : ncpgr=0
880 : !FB: @MT dimcprj_srt needed?
881 90 : ABI_MALLOC(tdks%dimcprj,(dtset%natom))
882 60 : ABI_MALLOC(tdks%dimcprj_srt,(dtset%natom))
883 : call pawcprj_getdim(tdks%dimcprj,dtset%natom,tdks%nattyp,dtset%ntypat, &
884 30 : & dtset%typat,pawtab,'R')
885 : call pawcprj_getdim(tdks%dimcprj_srt,dtset%natom,tdks%nattyp,dtset%ntypat, &
886 30 : & dtset%typat,pawtab,'O')
887 : !call pawcprj_alloc(tdks%cprj,ncpgr,dimcprj_srt)
888 30 : call pawcprj_alloc(tdks%cprj,ncpgr,tdks%dimcprj)
889 : !ABI_FREE(dimcprj_srt)
890 :
891 : !** Variables/arrays related to the fine FFT grid
892 158 : ABI_MALLOC(tdks%pawfgrtab,(my_natom))
893 30 : if (my_natom>0) then
894 : call pawtab_get_lsize(pawtab,l_size_atm,my_natom,dtset%typat, &
895 30 : & mpi_atmtab=mpi_enreg%my_atmtab)
896 : call pawfgrtab_init(tdks%pawfgrtab,cplex,l_size_atm,dtset%nspden,dtset%typat, &
897 30 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
898 30 : ABI_FREE(l_size_atm)
899 : end if
900 60 : tdks%usexcnhat=maxval(pawtab(:)%usexcnhat)
901 :
902 : !** Variables/arrays related to the PAW spheres
903 158 : ABI_MALLOC(tdks%paw_ij,(my_natom))
904 158 : ABI_MALLOC(tdks%paw_an,(my_natom))
905 30 : call paw_an_nullify(tdks%paw_an)
906 30 : call paw_ij_nullify(tdks%paw_ij)
907 30 : has_dijhat=0; if (dtset%iscf==22) has_dijhat=1
908 30 : has_vhartree=0; if (dtset%prtvha > 0 .or. dtset%prtvclmb > 0) has_vhartree=1
909 302 : has_dijnd=0;if(any(abs(dtset%nucdipmom)>tol8)) has_dijnd=1
910 30 : has_dijfock=0
911 30 : has_dijU=merge(0,1,dtset%usepawu>0) !Be careful on this!
912 30 : has_vxctau=pawxc_get_usekden(dtset%ixc)
913 : call paw_an_init(tdks%paw_an,dtset%natom,dtset%ntypat,0,0,dtset%nspden, &
914 : & cplex,dtset%pawxcdev,dtset%typat,pawang,pawtab,has_vxc=1, &
915 : & has_vxctau=has_vxctau,has_vxc_ex=1,has_vhartree=has_vhartree, &
916 30 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
917 : call paw_ij_init(tdks%paw_ij,cplex,dtset%nspinor,dtset%nsppol,dtset%nspden, &
918 : & dtset%pawspnorb,dtset%natom,dtset%ntypat,dtset%typat,pawtab, &
919 : & has_dij=1,has_dijfock=has_dijfock,has_dijhartree=1, &
920 : & has_dijnd=has_dijnd,has_dijso=1,has_dijhat=has_dijhat, &
921 : & has_dijU=has_dijU,has_pawu_occ=1,has_exexch_pot=1, &
922 : & nucdipmom=dtset%nucdipmom,comm_atom=mpi_enreg%comm_atom, &
923 30 : & mpi_atmtab=mpi_enreg%my_atmtab)
924 :
925 : !** Check for non-overlapping spheres
926 : call chkpawovlp(dtset%natom,psps%ntypat,dtset%pawovlp,pawtab,tdks%rmet, &
927 30 : & dtset%typat,tdks%xred)
928 :
929 : !** Identify parts of the rectangular grid where the density has to be calculated
930 30 : optcut=0;optgr0=dtset%pawstgylm;optgr1=0;optgr2=0;optrad=1-dtset%pawstgylm
931 30 : forces_needed=0 !FB TODO Maybe needs to be changed if Ehrenfest?
932 30 : stress_needed=0
933 : if ((forces_needed==1) .or. &
934 30 : & (dtset%xclevel==2 .and. dtset%pawnhatxc>0 .and. tdks%usexcnhat>0) .or. &
935 : & (dtset%positron/=0.and.forces_needed==2)) then
936 0 : optgr1=dtset%pawstgylm
937 0 : if (stress_needed==1) optrad=1; if (dtset%pawprtwf==1) optrad=1
938 : end if
939 : call nhatgrid(tdks%atindx1,tdks%gmet,my_natom,dtset%natom, &
940 : & tdks%nattyp,tdks%pawfgr%ngfft,psps%ntypat,optcut,optgr0,optgr1, &
941 : & optgr2,optrad,tdks%pawfgrtab,pawtab,tdks%rprimd,dtset%typat, &
942 : & tdks%ucvol,tdks%xred,comm_atom=mpi_enreg%comm_atom, &
943 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_fft=mpi_enreg%comm_fft, &
944 30 : & distribfft=mpi_enreg%distribfft)
945 :
946 30 : tdks%nhatgrdim=0;if (dtset%xclevel==2) tdks%nhatgrdim=tdks%usexcnhat*dtset%pawnhatxc
947 30 : if (tdks%nhatgrdim>0) then
948 0 : ABI_MALLOC(tdks%nhatgr,(cplex*tdks%nfftf,dtset%nspden,3*tdks%nhatgrdim))
949 : else
950 30 : ABI_MALLOC(tdks%nhatgr,(0,0,0))
951 : end if
952 :
953 120 : ABI_MALLOC(tdks%nhat,(tdks%nfftf,dtset%nspden*psps%usepaw))
954 :
955 : !Required in the PAW case to compute the inverse of the overlap (invovl) operator
956 60 : call init_invovl(dtset%nkpt)
957 : else
958 20 : ABI_MALLOC(tdks%nhat,(0,0))
959 20 : ABI_MALLOC(tdks%nhatgr,(0,0,0))
960 20 : tdks%nhatgrdim=0
961 : end if
962 :
963 : !Allocate various required arrays for calculation of the Hamiltonian
964 : !Potentials
965 150 : ABI_MALLOC(tdks%vhartr,(tdks%nfftf))
966 989924 : tdks%vhartr=zero
967 100 : ABI_MALLOC(tdks%vpsp,(tdks%nfftf))
968 989924 : tdks%vpsp=zero
969 200 : ABI_MALLOC(tdks%vtrial,(tdks%nfftf,dtset%nspden))
970 989974 : tdks%vtrial=zero
971 150 : ABI_MALLOC(tdks%vxc,(tdks%nfftf,dtset%nspden))
972 989974 : tdks%vxc=zero
973 50 : if (psps%n1xccc/=0) then
974 100 : ABI_MALLOC(tdks%xccc3d,(tdks%nfftf))
975 : else
976 0 : ABI_MALLOC(tdks%xccc3d,(0))
977 : end if
978 989924 : tdks%xccc3d=zero
979 50 : if (psps%usepaw==1) then
980 90 : ABI_MALLOC(tdks%xcctau3d,(tdks%nfftf*dtset%usekden))
981 30 : tdks%xcctau3d=zero
982 : endif
983 : !For mGGA
984 50 : usevxctau=merge(1,0,xc_need_kden(dtset%ixc))
985 250 : ABI_MALLOC(tdks%vxctau,(tdks%nfftf,dtset%nspden,4*usevxctau))
986 50 : tdks%vxctau=zero
987 : !For hybrid functionals
988 50 : use_hybcomp=0
989 50 : if(mod(dtset%fockoptmix,100)==11) use_hybcomp=1
990 200 : ABI_MALLOC(tdks%vxc_hybcomp,(tdks%pawfgr%nfft,dtset%nspden*use_hybcomp))
991 50 : tdks%vxc_hybcomp=zero
992 : !For VDW corrected functionals
993 50 : tdks%ngrvdw=0
994 50 : if ((dtset%vdw_xc>=5.and.dtset%vdw_xc<=7)) then
995 0 : tdks%ngrvdw=dtset%natom
996 : end if
997 150 : ABI_MALLOC(tdks%grvdw,(3,tdks%ngrvdw))
998 50 : tdks%grvdw=zero
999 :
1000 : !Compute large sphere G^2 cut-off (gsqcut) and box / sphere ratio
1001 50 : if (psps%usepaw==1) then
1002 : call getcut(boxcut,dtset%pawecutdg,tdks%gmet,tdks%gsqcut,dtset%iboxcut, &
1003 30 : & std_out,k0,tdks%pawfgr%ngfft)
1004 : else
1005 : call getcut(boxcut,dtset%ecut,tdks%gmet,tdks%gsqcut,dtset%iboxcut, &
1006 20 : & std_out,k0,tdks%pawfgr%ngfft)
1007 : end if
1008 :
1009 : !Compute structure factor phases (exp(2Pi i G.xred)) on coarse and fine grid
1010 150 : ABI_MALLOC(tdks%ph1d,(2,3*(2*tdks%pawfgr%mgfftc+1)*dtset%natom))
1011 150 : ABI_MALLOC(tdks%ph1df,(2,3*(2*tdks%pawfgr%mgfft+1)*dtset%natom))
1012 : call getph(tdks%atindx,dtset%natom,tdks%pawfgr%ngfftc(1),tdks%pawfgr%ngfftc(2), &
1013 50 : & tdks%pawfgr%ngfftc(3),tdks%ph1d,tdks%xred)
1014 50 : if (psps%usepaw==1.and.tdks%pawfgr%usefinegrid==1) then
1015 : call getph(tdks%atindx,dtset%natom,tdks%pawfgr%ngfft(1),tdks%pawfgr%ngfft(2), &
1016 30 : & tdks%pawfgr%ngfft(3),tdks%ph1df,tdks%xred)
1017 : else
1018 43940 : tdks%ph1df(:,:)=tdks%ph1d(:,:)
1019 : end if
1020 :
1021 : !Allocate memory for density
1022 200 : ABI_MALLOC(tdks%rhor,(tdks%nfftf,dtset%nspden))
1023 200 : ABI_MALLOC(tdks%taur,(tdks%nfftf,dtset%nspden*dtset%usekden))
1024 150 : ABI_MALLOC(tdks%rhog,(2,tdks%nfftf))
1025 150 : ABI_MALLOC(tdks%taug,(2,tdks%nfftf*dtset%usekden))
1026 :
1027 50 : end subroutine second_setup
1028 : !!***
1029 :
1030 : !!****f* m_rttddft_tdks/read_wfk
1031 : !!
1032 : !! NAME
1033 : !! read_wfk
1034 : !!
1035 : !! FUNCTION
1036 : !! Reads initial wavefunctions (KS orbitals) in WFK file (call inwffil)
1037 : !!
1038 : !! INPUTS
1039 : !! dtfil <type datafiles_type> = infos about file names, file unit numbers
1040 : !! dtset <type(dataset_type)> = all input variables for this dataset
1041 : !! ecut_eff <real(dp)> = effective PW cutoff energy
1042 : !! mpi_enreg <MPI_type> = MPI-parallelisation information
1043 : !! tdks <type(tdks_type)> = the tdks object to initialize
1044 : !!
1045 : !! OUTPUT
1046 : !!
1047 : !! SOURCE
1048 50 : subroutine read_wfk(dtfil, dtset, ecut_eff, fname_wfk, mpi_enreg, tdks)
1049 :
1050 : !Arguments ------------------------------------
1051 : !scalars
1052 : character(len=fnlen), intent(in) :: fname_wfk
1053 : real(dp), intent(in) :: ecut_eff
1054 : type(datafiles_type), intent(in) :: dtfil
1055 : type(dataset_type), intent(inout) :: dtset
1056 : type(MPI_type), intent(inout) :: mpi_enreg
1057 : type(tdks_type), intent(inout) :: tdks
1058 :
1059 : !Local variables-------------------------------
1060 : !scalars
1061 : integer,parameter :: formeig=0
1062 : integer :: ask_accurate
1063 : integer :: band
1064 : integer :: cnt
1065 : integer :: ierr, ikpt
1066 : integer :: my_nspinor
1067 : integer :: optorth
1068 : integer :: spin
1069 : type(wffile_type) :: wff1, wffnow
1070 : !arrays
1071 : character(len=500) :: msg
1072 :
1073 : ! ***********************************************************************
1074 :
1075 : !If paral_kgb == 0, it may happen that some processors are idle (no entry in proc_distrb)
1076 : !but mkmem == nkpt and this can cause integer overflow in mcg or allocation error.
1077 : !Here we count the number of states treated by the proc. if cnt == 0, mcg is then set to 0.
1078 50 : cnt = 0
1079 100 : do spin=1,dtset%nsppol
1080 1008 : do ikpt=1,dtset%nkpt
1081 10046 : do band=1,dtset%nband(ikpt + (spin-1) * dtset%nkpt)
1082 9996 : if (.not. proc_distrb_cycle(mpi_enreg%proc_distrb, ikpt, band, band, spin, mpi_enreg%me_kpt)) cnt = cnt + 1
1083 : end do
1084 : end do
1085 : end do
1086 :
1087 50 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
1088 50 : tdks%mcg=dtset%mpw*my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol
1089 50 : if (cnt == 0) then
1090 0 : tdks%mcg = 0
1091 0 : write(msg,"(2(a,i0))")"rank: ",mpi_enreg%me, "does not have wavefunctions to treat. Setting mcg to: ",tdks%mcg
1092 0 : ABI_WARNING(msg)
1093 : end if
1094 :
1095 50 : if (dtset%usewvl == 0 .and. dtset%mpw > 0 .and. cnt /= 0)then
1096 50 : if (my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol > floor(real(HUGE(0))/real(dtset%mpw) )) then
1097 : ierr = 0
1098 : write (msg,'(9a)')&
1099 0 : & "Default integer is not wide enough to store the size of the wavefunction array (mcg).",ch10,&
1100 0 : & "This usually happens when paral_kgb == 0 and there are not enough procs to distribute kpts and spins",ch10,&
1101 0 : & "Action: if paral_kgb == 0, use nprocs = nkpt * nsppol to reduce the memory per node.",ch10,&
1102 0 : & "If tdks does not solve the problem, use paral_kgb 1 with nprocs > nkpt * nsppol and use npfft/npband/npspinor",ch10,&
1103 0 : & "to decrease the memory requirements. Consider also OpenMP threads."
1104 0 : ABI_ERROR_NOSTOP(msg,ierr)
1105 : write (msg,'(5(a,i0), 2a)')&
1106 0 : & "my_nspinor: ",my_nspinor, ", mpw: ",dtset%mpw, ", mband: ",dtset%mband,&
1107 0 : & ", mkmem: ",dtset%mkmem, ", nsppol: ",dtset%nsppol,ch10,&
1108 0 : & 'Note: Compiling with large int (int64) requires a full software stack (MPI/FFTW/BLAS...) compiled in int64 mode'
1109 0 : ABI_ERROR(msg)
1110 : end if
1111 : end if
1112 :
1113 : ! Alloc size for wfk and bands
1114 150 : ABI_MALLOC_OR_DIE(tdks%cg,(2,tdks%mcg),ierr)
1115 150 : ABI_MALLOC(tdks%eigen,(dtset%mband*dtset%nkpt*dtset%nsppol))
1116 100 : ABI_MALLOC(tdks%eigen0,(dtset%mband*dtset%nkpt*dtset%nsppol))
1117 :
1118 9138 : tdks%eigen(:) = zero
1119 50 : ask_accurate=1
1120 :
1121 : !Actually read the intial KS orbitals here
1122 50 : if (dtset%td_restart /= 1) then
1123 30 : write(msg,'(3a)') ch10,'------------------- Reading initial wavefunctions -------------------',ch10
1124 : else
1125 20 : write(msg,'(3a)') ch10,'------------------- Reading wavefunctions for restart ------------------',ch10
1126 : end if
1127 50 : call wrtout(ab_out,msg)
1128 50 : if (do_write_log) call wrtout(std_out,msg)
1129 50 : wff1%unwff=dtfil%unwff1
1130 50 : optorth=0 !No need to orthogonalize the wfk
1131 650 : tdks%hdr%rprimd=tdks%rprimd
1132 4724594 : tdks%cg=zero
1133 : call inwffil(ask_accurate,tdks%cg,dtset,dtset%ecut,ecut_eff,tdks%eigen, &
1134 : & dtset%exchn2n3d,formeig,tdks%hdr,1,dtset%istwfk,tdks%kg, &
1135 : & dtset%kptns,dtset%localrdwf,dtset%mband,tdks%mcg,dtset%mkmem, &
1136 : & mpi_enreg,dtset%mpw,dtset%nband,tdks%pawfgr%ngfft,dtset%nkpt, &
1137 : & tdks%npwarr,dtset%nsppol,dtset%nsym,dtset%occ_orig,optorth, &
1138 : & dtset%symafm,dtset%symrel,dtset%tnons,dtfil%unkg,wff1,wffnow, &
1139 50 : & dtfil%unwff1,fname_wfk,tdks%wvl)
1140 :
1141 : !Close file
1142 50 : call WffClose(wff1,ierr)
1143 :
1144 : !Keep initial eigenvalues in memory
1145 9138 : tdks%eigen0(:) = tdks%eigen(:)
1146 :
1147 : !In case of restart also read wfk file containing wave functions at t=0
1148 50 : if (dtset%td_restart == 1 .and. tdks%fname_wfk0 /= fname_wfk) then
1149 20 : write(msg,'(3a)') ch10,'------------------- Reading initial wavefunctions -------------------',ch10
1150 60 : ABI_MALLOC_OR_DIE(tdks%cg0,(2,tdks%mcg),ierr)
1151 20 : call wrtout(ab_out,msg)
1152 20 : if (do_write_log) call wrtout(std_out,msg)
1153 1255700 : tdks%cg0=zero
1154 : call inwffil(ask_accurate,tdks%cg0,dtset,dtset%ecut,ecut_eff,tdks%eigen0, &
1155 : & dtset%exchn2n3d,formeig,tdks%hdr,1,dtset%istwfk,tdks%kg, &
1156 : & dtset%kptns,dtset%localrdwf,dtset%mband,tdks%mcg,dtset%mkmem, &
1157 : & mpi_enreg,dtset%mpw,dtset%nband,tdks%pawfgr%ngfft,dtset%nkpt, &
1158 : & tdks%npwarr,dtset%nsppol,dtset%nsym,dtset%occ_orig,optorth, &
1159 : & dtset%symafm,dtset%symrel,dtset%tnons,dtfil%unkg,wff1,wffnow, &
1160 20 : & dtfil%unwff1,tdks%fname_wfk0,tdks%wvl)
1161 :
1162 : !Close file
1163 20 : call WffClose(wff1,ierr)
1164 :
1165 : end if
1166 :
1167 50 : end subroutine read_wfk
1168 : !!***
1169 :
1170 0 : end module m_rttddft_tdks
1171 : !!***
|
