Line data Source code
1 : !!****m* ABINIT/m_gwr_driver
2 : !! NAME
3 : !! m_gwr_driver
4 : !!
5 : !! FUNCTION
6 : !! Driver for GWR calculations
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2021-2026 ABINIT group (MG)
10 : !! This file is distributed under the terms of the
11 : !! APACHE license version 2.0, see ~abinit/COPYING
12 : !! or https://www.apache.org/licenses/LICENSE-2.0 .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_gwr_driver
23 :
24 : use, intrinsic :: iso_c_binding
25 : USE_MPI
26 : use defs_basis
27 : use defs_wvltypes
28 : use m_errors
29 : use m_abicore
30 : use m_xmpi
31 : use m_xomp
32 : use m_hdr
33 : use libxc_functionals
34 : use m_crystal
35 : use m_ebands
36 : use m_dtset
37 : use m_dtfil
38 : use m_wfk
39 : use netcdf
40 : use m_nctk
41 :
42 : use defs_datatypes, only : pseudopotential_type
43 : use defs_abitypes, only : MPI_type
44 : use m_time, only : timab
45 : use m_io_tools, only : file_exists, open_file, get_unit, iomode_from_fname
46 : use m_time, only : cwtime, cwtime_report, sec2str
47 : use m_fstrings, only : strcat, sjoin, ftoa, itoa, string_in, ltoa
48 : use m_fftcore, only : print_ngfft, get_kg
49 : use m_fft, only : fourdp
50 : use m_ioarr, only : read_rhor
51 : use m_energies, only : energies_type
52 : use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq
53 : use m_pawang, only : pawang_type
54 : use m_pawrad, only : pawrad_type
55 : use m_pawtab, only : pawtab_type, pawtab_print, pawtab_get_lsize
56 : use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify
57 : use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify
58 : use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free, pawfgrtab_print
59 : use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free, &
60 : pawrhoij_inquire_dim, pawrhoij_symrhoij, pawrhoij_unpack
61 : use m_pawdij, only : pawdij, symdij_all
62 : use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy
63 : use m_paw_pwaves_lmn, only : paw_pwaves_lmn_t, paw_pwaves_lmn_init, paw_pwaves_lmn_free
64 : use m_pawpwij, only : pawpwff_t, pawpwff_init, pawpwff_free, paw_rho_tw_g
65 : use m_kg, only : getph
66 : use m_wfd, only : wfd_t, test_charge
67 : use m_pspini, only : pspini
68 : use m_paw_correlations,only : pawpuxinit
69 : use m_paw_dmft, only : paw_dmft_type
70 : use m_paw_sphharm, only : setsym_ylm
71 : use m_paw_mkrho, only : denfgr
72 : use m_paw_nhat, only : nhatgrid, pawmknhat
73 : use m_paw_tools, only : chkpawovlp, pawprt
74 : use m_paw_denpot, only : pawdenpot
75 : use m_paw_init, only : pawinit, paw_gencond
76 : use m_pawcprj, only : pawcprj_type, pawcprj_free, pawcprj_alloc ! , paw_overlap
77 : use m_pawxc, only : pawxc_get_usekden
78 : use m_ksdiago, only : ugb_t, hyb_t
79 : use m_mkrho, only : prtrhomxmn
80 : use m_melemts, only : melflags_t
81 : use m_setvtr, only : setvtr
82 : use m_vhxc_me, only : calc_vhxc_me
83 : use m_gwr, only : gwr_t
84 : use m_vcoul, only : vcgen_t
85 : use m_ephtk, only : ephtk_update_ebands
86 : use m_pstat, only : pstat_proc
87 :
88 : implicit none
89 :
90 : private
91 : !!***
92 :
93 : public :: gwr_driver
94 :
95 : #ifdef HAVE_MPI1
96 : include 'mpif.h'
97 : #endif
98 : !!***
99 :
100 : contains
101 : !!***
102 :
103 : !!****f* m_gwr_driver/gwr_driver
104 : !! NAME
105 : !! gwr_driver
106 : !!
107 : !! FUNCTION
108 : !! Main routine for GWR calculations.
109 : !!
110 : !! INPUTS
111 : !! acell(3)=Length scales of primitive translations (bohr)
112 : !! codvsn=Code version
113 : !! dtfil<datafiles_type>=Variables related to files.
114 : !! dtset<dataset_type>=All input variables for this dataset.
115 : !! pawang<pawang_type)>=PAW angular mesh and related data.
116 : !! pawrad(ntypat*usepaw)<pawrad_type>=Paw radial mesh and related data.
117 : !! pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data.
118 : !! psps<pseudopotential_type>=Variables related to pseudopotentials.
119 : !! Before entering the first time in the routine, a significant part of Psps has been initialized :
120 : !! the integers dimekb,lmnmax,lnmax,mpssang,mpssoang,mpsso,mgrid,ntypat,n1xccc,usepaw,useylm,
121 : !! and the arrays dimensioned to npsp. All the remaining components of Psps are to be initialized in
122 : !! the call to pspini. The next time the code enters bethe_salpeter, Psps might be identical to the
123 : !! one of the previous Dtset, in which case, no reinitialisation is scheduled in pspini.F90.
124 : !! xred(3,natom)=Reduced atomic coordinates.
125 : !!
126 : !! NOTES
127 : !!
128 : !! ON THE USE OF FFT GRIDS:
129 : !! =================
130 : !! In case of PAW:
131 : !! ---------------
132 : !! Two FFT grids are used:
133 : !! - A "coarse" FFT grid (defined by ecut) for the application of the Hamiltonian on the plane waves basis.
134 : !! It is defined by nfft, ngfft, mgfft, ...
135 : !! Hamiltonian, wave-functions, density related to WFs (rhor here), ... are expressed on this grid.
136 : !! - A "fine" FFT grid (defined) by ecutdg) for the computation of the density inside PAW spheres.
137 : !! It is defined by nfftf, ngfftf, mgfftf, ... Total density, potentials, ... are expressed on this grid.
138 : !! In case of norm-conserving:
139 : !! ---------------------------
140 : !! - Only the usual FFT grid (defined by ecut) is used. It is defined by nfft, ngfft, mgfft, ...
141 : !! For compatibility reasons, (nfftf,ngfftf,mgfftf) are set equal to (nfft,ngfft,mgfft) in that case.
142 : !!
143 : !! SOURCE
144 :
145 0 : subroutine gwr_driver(codvsn, dtfil, dtset, pawang, pawrad, pawtab, psps, xred)
146 :
147 : !Arguments ------------------------------------
148 : !scalars
149 : character(len=8),intent(in) :: codvsn
150 : type(datafiles_type),intent(in) :: dtfil
151 : type(dataset_type),intent(inout) :: dtset
152 : type(pawang_type),intent(inout) :: pawang
153 : type(pseudopotential_type),intent(inout) :: psps
154 : !arrays
155 : real(dp),intent(in) :: xred(3,dtset%natom)
156 : type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw)
157 : type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
158 :
159 : !Local variables ------------------------------
160 : !scalars
161 : integer,parameter :: master = 0, cplex1 = 1, ipert0 = 0, idir0 = 0, optrhoij1 = 1
162 : integer :: ii, comm, nprocs, my_rank, mgfftf, nfftf, omp_ncpus, work_size, nks_per_proc
163 : integer :: ierr, spin, ik_ibz, nband_k, iomode__, color, io_comm, usevxctau_paw !, kg_varid
164 : real(dp) :: eff, mempercpu_mb, max_wfsmem_mb, nonscal_mem, el_temp
165 : real(dp) :: ecore, ecut_eff, ecutdg_eff, cpu, wall, gflops, diago_cpu, diago_wall, diago_gflops
166 : logical, parameter :: is_dfpt = .false.
167 : logical :: read_wfk, write_wfk, cc4s_task, rectangular, with_pools, rdm_update, call_pawinit, cc4s_from_wfk
168 : character(len=500) :: msg
169 : character(len=fnlen) :: wfk_path, den_path, kden_path, out_path
170 0 : type(hdr_type) :: wfk_hdr, den_hdr, kden_hdr, owfk_hdr
171 0 : type(crystal_t) :: cryst, den_cryst, wfk_cryst
172 0 : type(ebands_t) :: ks_ebands, owfk_ebands, tmp_ebands
173 : type(pawfgr_type) :: pawfgr
174 0 : type(wvl_data) :: wvl
175 0 : type(mpi_type) :: mpi_enreg_seq
176 0 : type(gwr_t) :: gwr
177 0 : type(wfk_t) :: owfk
178 : !arrays
179 : real(dp), parameter :: k0(3) = zero
180 : integer :: cplex, cplex_dij, cplex_rhoij
181 : integer :: gnt_option,has_dijU,has_dijso,ider,izero
182 : integer :: istep, moved_atm_inside, moved_rhor, n3xccc, sc_mode
183 : !integer :: ngrvdw,nhatgrdim,nkxc,nkxc1,nprocs,nscf,nspden_rhoij,nzlmopt,optene
184 : integer :: ndij !,ndim,nfftf,nfftf_tot,nkcalc,gwc_nfft,gwc_nfftot,gwx_nfft,gwx_nfftot
185 : integer :: ngrvdw, nhatgrdim, nkxc, nkxc1, nspden_rhoij, optene, nzlmopt
186 : integer :: optcut, optgr0, optgr1, optgr2, optrad, psp_gencond, option
187 : integer :: rhoxsp_method, usexcnhat !, use_umklp
188 : real(dp) :: compch_fft, compch_sph !,r_s,rhoav,alpha
189 : !real(dp) :: drude_plsmf !,my_plsmf,ecut_eff,ecutdg_eff,ehartree
190 : real(dp) :: gsqcutc_eff, gsqcutf_eff, gsqcut_shp, vxcavg, gw_gsq, gs_fermie
191 : type(energies_type) :: KS_energies
192 : type(melflags_t) :: KS_mflags
193 0 : type(paw_dmft_type) :: Paw_dmft
194 0 : type(ugb_t) :: ugb
195 0 : type(hyb_t) :: hyb
196 0 : type(xmpi_pool2d_t) :: diago_pool
197 : !arrays
198 : integer :: ngfftc(18),ngfftf(18),units(2) !, grid3(3)
199 0 : integer,allocatable :: nq_spl(:), l_size_atm(:)
200 0 : integer,allocatable :: tmp_kstab(:,:,:), npwarr_ik(:), gvec_(:,:), istwfk_ik(:), nband_iks(:,:)
201 0 : real(dp) :: strsxc(6), diago_info(3, dtset%nkpt, dtset%nsppol),tsec(2)
202 0 : real(dp),allocatable :: grchempottn(:,:),grewtn(:,:),grvdw(:,:),qmax(:)
203 0 : real(dp),allocatable :: ks_nhat(:,:),ks_nhatgr(:,:,:),ks_rhog(:,:)
204 0 : real(dp),allocatable :: ks_rhor(:,:),ks_vhartr(:), ks_vtrial(:,:), ks_vxc(:,:)
205 0 : real(dp),allocatable :: ks_taur(:,:) !, ks_vxctau(:,:), xcctau3d(:)
206 0 : real(dp),allocatable :: kxc(:,:), ph1d(:,:), ph1df(:,:) !qp_kxc(:,:),
207 0 : real(dp),allocatable :: vpsp(:), xccc3d(:), dijexc_core(:,:,:) !, dij_hf(:,:,:)
208 0 : real(dp),allocatable :: eig_k(:), occ_k(:)
209 0 : real(dp),contiguous,pointer :: cg_k_ptr(:,:)
210 0 : type(paw_an_type),allocatable :: KS_paw_an(:)
211 0 : type(paw_ij_type),allocatable :: KS_paw_ij(:)
212 0 : type(pawfgrtab_type),allocatable :: Pawfgrtab(:)
213 0 : type(pawrhoij_type),allocatable :: KS_Pawrhoij(:)
214 0 : type(pawpwff_t),allocatable :: Paw_pwff(:)
215 : !type(pawcprj_type),allocatable :: cprj_k(:,:)
216 : !************************************************************************
217 :
218 : ! This part performs the initialization of the basic objects used to perform e-ph calculations:
219 : !
220 : ! 1) Crystal structure `cryst`
221 : ! 2) Ground state band energies: `ks_ebands`
222 : ! 5) Pseudos and PAW basic objects.
223 : !
224 : ! Once we have these objects, we can call specialized routines for e-ph calculations.
225 : ! Notes:
226 : !
227 : ! * Any modification to the basic objects mentioned above should be done here (e.g. change of efermi)
228 : ! * This routines shall not allocate big chunks of memory. The CPU-demanding sections should be
229 : ! performed in the subdriver that will employ different MPI distribution schemes optimized for that particular task.
230 :
231 : ! abirules!
232 : if (.False.) write(std_out,*)xred
233 :
234 0 : comm = xmpi_world; nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
235 0 : units(:) = [std_out, ab_out]
236 :
237 0 : call cwtime(cpu, wall, gflops, "start")
238 0 : call pstat_proc%print(_PSTAT_ARGS_)
239 :
240 : ! write(msg,'(a)')&
241 : ! ' GWR: Calculation of the GW corrections with GWR code ',ch10,ch10,&
242 : ! call wrtout(units, msg)
243 : !
244 : #if defined HAVE_GW_DPC
245 0 : write(msg,'(a,i2,a)')'.Using double precision arithmetic; gwpc = ',gwp,ch10
246 : #else
247 : write(msg,'(a,i2,a)')'.Using single precision arithmetic; gwpc = ',gwp,ch10
248 : #endif
249 0 : call wrtout(units, msg)
250 :
251 : ! autoparal section
252 : ! TODO: This just to activate autoparal in AbiPy. Lot of things should be improved.
253 0 : if (dtset%max_ncpus /= 0) then
254 0 : write(ab_out,'(a)')"--- !Autoparal"
255 0 : write(ab_out,"(a)")"# Autoparal section for GWR runs"
256 0 : write(ab_out,"(a)") "info:"
257 0 : write(ab_out,"(a,i0)")" autoparal: ",dtset%autoparal
258 0 : write(ab_out,"(a,i0)")" max_ncpus: ",dtset%max_ncpus
259 0 : write(ab_out,"(a,i0)")" nkpt: ",dtset%nkpt
260 0 : write(ab_out,"(a,i0)")" nsppol: ",dtset%nsppol
261 0 : write(ab_out,"(a,i0)")" nspinor: ",dtset%nspinor
262 0 : write(ab_out,"(a,i0)")" mband: ",dtset%mband
263 0 : write(ab_out,"(3a)") " gwr_task: '",trim(dtset%gwr_task),"'"
264 :
265 0 : if (string_in(dtset%gwr_task, "HDIAGO, HDIAGO_FULL, CC4S, CC4S_FULL")) then
266 0 : work_size = dtset%nkpt * dtset%nsppol * dtset%mpw
267 0 : max_wfsmem_mb = (two * dp * dtset%mpw * dtset%mband * dtset%nkpt * dtset%nsppol * dtset%nspinor * b2Mb) * 1.1_dp
268 : else
269 0 : work_size = dtset%gwr_ntau * dtset%nkpt * dtset%nsppol * dtset%mpw
270 0 : max_wfsmem_mb = (two * dp * dtset%mpw * dtset%mband * dtset%nkpt * dtset%nsppol * dtset%nspinor * b2Mb) * 1.1_dp
271 : end if
272 : ! Non-scalable memory in Mb i.e. memory that is not distributed with MPI.
273 0 : nonscal_mem = zero
274 :
275 : ! List of configurations.
276 : ! Assuming an OpenMP implementation with perfect speedup!
277 0 : write(ab_out,"(a)")"configurations:"
278 :
279 0 : do ii=1,dtset%max_ncpus
280 0 : nks_per_proc = work_size / ii
281 0 : nks_per_proc = nks_per_proc + mod(work_size, ii)
282 0 : eff = (one * work_size) / (ii * nks_per_proc)
283 : ! Add the non-scalable part and increase by 10% to account for other datastructures.
284 0 : mempercpu_mb = (max_wfsmem_mb + nonscal_mem) * 1.1_dp
285 0 : do omp_ncpus=1,1 !xomp_get_max_threads()
286 0 : write(ab_out,"(a,i0)")" - tot_ncpus: ",ii * omp_ncpus
287 0 : write(ab_out,"(a,i0)")" mpi_ncpus: ",ii
288 0 : write(ab_out,"(a,i0)")" omp_ncpus: ",omp_ncpus
289 0 : write(ab_out,"(a,f12.9)")" efficiency: ",eff
290 0 : write(ab_out,"(a,f12.2)")" mem_per_cpu: ",mempercpu_mb
291 : end do
292 : end do
293 0 : write(ab_out,'(a)')"..."
294 0 : call wrtout(std_out, "Stopping now!")
295 0 : call abi_abort('PERS', exit_status=0, print_config=.false.)
296 : end if
297 :
298 0 : cryst = dtset%get_crystal(img=1)
299 :
300 : ! Some variables need to be initialized/nullify at start
301 0 : usexcnhat = 0
302 0 : call KS_energies%init()
303 :
304 : !Get electronic temperature from dtset
305 0 : el_temp = merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9)
306 :
307 0 : den_path = dtfil%fildensin; wfk_path = dtfil%fnamewffk; kden_path = dtfil%filkdensin
308 : !use_den = f (string_in(dtset%gwr_task, "CC4S_FROM_WFK")) then
309 :
310 0 : if (my_rank == master) then
311 : ! Initialize filenames. Accept files in Fortran or in netcdf format.
312 0 : if (nctk_try_fort_or_ncfile(den_path, msg) /= 0) then
313 0 : ABI_ERROR(sjoin("Cannot find DEN file:", den_path, ". Error:", msg))
314 : end if
315 0 : call wrtout(units, sjoin("- Reading GS density from: ", den_path))
316 :
317 0 : if (dtset%usekden == 1) then
318 0 : if (nctk_try_fort_or_ncfile(kden_path, msg) /= 0) then
319 0 : ABI_ERROR(sjoin("Cannot find KDEN file:", kden_path, ". Error:", msg))
320 : end if
321 0 : call wrtout(units, sjoin("- Reading KDEN kinetic energy density from: ", kden_path))
322 : end if
323 0 : call wrtout(ab_out, ch10//ch10)
324 : end if ! master
325 :
326 : ! Broadcast filenames (needed if we are using netcdf files)
327 0 : call xmpi_bcast(den_path, master, comm, ierr)
328 0 : call xmpi_bcast(kden_path, master, comm, ierr)
329 :
330 : ! TODO: FFT meshes for DEN/POT should be initialized from the DEN file instead of the dtset.
331 : ! Interpolating the DEN indeed breaks degeneracies in the vxc matrix elements.
332 : call pawfgr_init(pawfgr, dtset, mgfftf, nfftf, ecut_eff, ecutdg_eff, ngfftc, ngfftf, &
333 0 : gsqcutc_eff=gsqcutc_eff, gsqcutf_eff=gsqcutf_eff, gmet=cryst%gmet, k0=k0)
334 :
335 0 : call print_ngfft([std_out], ngfftc, header='Coarse FFT mesh for the wavefunctions')
336 0 : call print_ngfft([std_out], ngfftf, header='Dense FFT mesh for densities and potentials')
337 :
338 : ! Fake MPI_type for the sequential part.
339 0 : call initmpi_seq(mpi_enreg_seq)
340 0 : call mpi_enreg_seq%distribfft%init_seq('c', ngfftc(2), ngfftc(3), 'all')
341 0 : call mpi_enreg_seq%distribfft%init_seq('f', ngfftf(2), ngfftf(3), 'all')
342 :
343 : ! ===========================================
344 : ! === Open and read pseudopotential files ===
345 : ! ===========================================
346 0 : call pspini(dtset, dtfil, ecore, psp_gencond, gsqcutc_eff, gsqcutf_eff, pawrad, pawtab, psps, cryst%rprimd, comm_mpi=comm)
347 :
348 : ! ============================
349 : ! ==== PAW initialization ====
350 : ! ============================
351 0 : if (dtset%usepaw == 1) then
352 0 : call chkpawovlp(cryst%natom, cryst%ntypat, dtset%pawovlp, pawtab, cryst%rmet, cryst%typat, cryst%xred)
353 :
354 0 : cplex_dij = dtset%nspinor; cplex = 1; ndij = 1
355 :
356 0 : ABI_MALLOC(ks_pawrhoij, (cryst%natom))
357 : call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij, nspden_rhoij=nspden_rhoij, &
358 0 : nspden=dtset%nspden, spnorb=dtset%pawspnorb, cpxocc=dtset%pawcpxocc)
359 0 : call pawrhoij_alloc(ks_pawrhoij, cplex_rhoij, nspden_rhoij, dtset%nspinor, dtset%nsppol, cryst%typat, pawtab=pawtab)
360 :
361 : ! Test if we have to call pawinit
362 0 : gnt_option = 1; if (dtset%pawxcdev == 2 .or. (dtset%pawxcdev == 1 .and. dtset%positron /= 0)) gnt_option = 2
363 0 : call paw_gencond(dtset, gnt_option, "test", call_pawinit)
364 : !call_pawinit = .True.
365 :
366 0 : if (psp_gencond == 1 .or. call_pawinit) then
367 0 : call timab(553, 1, tsec)
368 0 : gsqcut_shp = two * abs(dtset%diecut) * dtset%dilatmx**2 / pi**2
369 : call pawinit(dtset%effmass_free, gnt_option, gsqcut_shp, zero, dtset%pawlcutd, dtset%pawlmix, &
370 : psps%mpsang, dtset%pawnphi, cryst%nsym, dtset%pawntheta, pawang, pawrad, &
371 0 : dtset%pawspnorb, pawtab, dtset%pawxcdev, dtset%ixc, dtset%usepotzero)
372 0 : call timab(553,2,tsec)
373 :
374 : ! Update internal values
375 0 : call paw_gencond(dtset, gnt_option, "save", call_pawinit)
376 : else
377 0 : if (pawtab(1)%has_kij ==1) pawtab(1:cryst%ntypat)%has_kij = 2
378 0 : if (pawtab(1)%has_nabla==1) pawtab(1:cryst%ntypat)%has_nabla = 2
379 : end if
380 :
381 0 : psps%n1xccc = maxval(pawtab(1:cryst%ntypat)%usetcore)
382 :
383 : ! Initialize optional flags in Pawtab to zero
384 : ! Cannot be done in Pawinit since the routine is called only if some parts. are changed
385 0 : pawtab(:)%has_nabla = 0
386 0 : pawtab(:)%lamb_shielding = zero
387 :
388 0 : call setsym_ylm(cryst%gprimd, pawang%l_max-1, cryst%nsym, dtset%pawprtvol, cryst%rprimd, cryst%symrec, pawang%zarot)
389 :
390 : ! Initialize and compute data for DFT+U
391 0 : Paw_dmft%use_dmft = Dtset%usedmft
392 : call pawpuxinit(Dtset%dmatpuopt,Dtset%exchmix,Dtset%f4of2_sla,Dtset%f6of2_sla, &
393 : is_dfpt,Dtset%jpawu,Dtset%lexexch,Dtset%lpawu,dtset%nspinor,Cryst%ntypat,dtset%optdcmagpawu,Pawang,Dtset%pawprtvol, &
394 0 : Pawrad,Pawtab,Dtset%upawu,Dtset%usedmft,Dtset%useexexch,Dtset%usepawu,dtset%ucrpa)
395 :
396 0 : if (my_rank == master) call pawtab_print(Pawtab)
397 :
398 : ! Get Pawrhoij from the header of the WFK file.
399 : !call pawrhoij_copy(wfk_hdr%pawrhoij, KS_Pawrhoij)
400 :
401 : ! Evaluate form factor of radial part of phi.phj-tphi.tphj.
402 0 : gw_gsq = max(Dtset%ecutsigx, Dtset%ecuteps) / (two*pi**2)
403 :
404 : ! Set up q-grid, make qmax 20% larger than largest expected.
405 0 : ABI_MALLOC(nq_spl, (Psps%ntypat))
406 0 : ABI_MALLOC(qmax, (Psps%ntypat))
407 0 : qmax = SQRT(gw_gsq)*1.2d0
408 0 : nq_spl = Psps%mqgrid_ff
409 : ! write(std_out,*)"using nq_spl",nq_spl,"qmax=",qmax
410 :
411 0 : rhoxsp_method = 1 ! Arnaud-Alouani (default in sigma)
412 : !rhoxsp_method = 2 ! Shiskin-Kresse
413 0 : if (Dtset%pawoptosc /= 0) rhoxsp_method = Dtset%pawoptosc
414 :
415 0 : ABI_MALLOC(paw_pwff, (psps%ntypat))
416 0 : call pawpwff_init(Paw_pwff, rhoxsp_method, nq_spl, qmax, cryst%gmet, pawrad, pawtab, psps)
417 :
418 0 : ABI_FREE(nq_spl)
419 0 : ABI_FREE(qmax)
420 :
421 : ! Variables/arrays related to the fine FFT grid
422 0 : ABI_CALLOC(ks_nhat, (nfftf, Dtset%nspden))
423 :
424 0 : ABI_MALLOC(pawfgrtab, (cryst%natom))
425 0 : call pawtab_get_lsize(pawtab, l_size_atm, cryst%natom, cryst%typat)
426 :
427 : cplex = 1
428 0 : call pawfgrtab_init(pawfgrtab, cplex, l_size_atm, dtset%nspden, dtset%typat)
429 0 : ABI_FREE(l_size_atm)
430 0 : compch_fft=greatest_real
431 0 : usexcnhat = maxval(Pawtab(:)%usexcnhat)
432 : ! * 0 if Vloc in atomic data is Vbare (Blochl's formulation)
433 : ! * 1 if Vloc in atomic data is VH(tnzc) (Kresse's formulation)
434 0 : call wrtout(std_out, sjoin(' using usexcnhat: ', itoa(usexcnhat)))
435 : !
436 : ! Identify parts of the rectangular grid where the density has to be calculated
437 0 : optcut = 0; optgr0 = Dtset%pawstgylm; optgr1 = 0; optgr2 = 0; optrad = 1 - Dtset%pawstgylm
438 0 : if (Dtset%pawcross==1) optrad=1
439 0 : if (Dtset%xclevel==2.and.usexcnhat>0) optgr1=Dtset%pawstgylm
440 :
441 : call nhatgrid(cryst%atindx1, cryst%gmet, cryst%natom, cryst%natom, cryst%nattyp, ngfftf, cryst%ntypat,&
442 0 : optcut,optgr0,optgr1,optgr2,optrad,Pawfgrtab,Pawtab,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred)
443 :
444 0 : call pawfgrtab_print(Pawfgrtab,Cryst%natom,unit=std_out,prtvol=Dtset%pawprtvol)
445 :
446 : else
447 0 : ABI_MALLOC(Paw_pwff, (0))
448 0 : ABI_MALLOC(Pawfgrtab, (0))
449 : end if ! End of PAW Initialization
450 :
451 : ! Allocate these arrays anyway, since they are passed to subroutines.
452 0 : ABI_MALLOC_IFNOT(ks_nhat, (nfftf, 0))
453 0 : ABI_MALLOC_IFNOT(dijexc_core, (1, 1, 0))
454 :
455 : !=============================================
456 : ! Read density and compare crystal structures
457 : ! ============================================
458 0 : ABI_MALLOC(ks_rhor, (nfftf, dtset%nspden))
459 :
460 : call read_rhor(den_path, cplex1, dtset%nspden, nfftf, ngfftf, dtset%usepaw, mpi_enreg_seq, ks_rhor, &
461 0 : den_hdr, ks_pawrhoij, comm, allow_interp=.False., want_varname="density")
462 :
463 0 : den_cryst = den_hdr%get_crystal()
464 0 : if (cryst%compare(den_cryst, header=" Comparing input crystal with DEN crystal") /= 0) then
465 0 : ABI_ERROR("Crystal structure from input and from DEN file do not agree! Check messages above!")
466 : end if
467 : ! Get fermie from the GS calculation.
468 : ! NB: It might understimate the real fermi level, especially if the den was computed on a shifted k-mesh
469 : ! at present it's only used to implement pseudobands
470 0 : gs_fermie = den_hdr%fermie
471 0 : call den_cryst%free(); call den_hdr%free()
472 :
473 0 : ABI_MALLOC(ks_taur, (nfftf, dtset%nspden * dtset%usekden))
474 0 : if (dtset%usekden == 1) then
475 : call read_rhor(kden_path, cplex1, dtset%nspden, nfftf, ngfftf, 0, mpi_enreg_seq, ks_taur, &
476 0 : kden_hdr, ks_pawrhoij, comm, allow_interp=.False., want_varname="kinedr")
477 0 : call kden_hdr%free()
478 0 : call prtrhomxmn(std_out, mpi_enreg_seq, nfftf, ngfftf, dtset%nspden, 1, ks_taur, optrhor=1, ucvol=cryst%ucvol)
479 : end if
480 :
481 : !========================================
482 : !==== Additional computation for PAW ====
483 : !========================================
484 0 : nhatgrdim = 0
485 0 : if (dtset%usepaw == 1) then
486 : ! Calculate the compensation charge nhat.
487 0 : if (Dtset%xclevel==2) nhatgrdim = usexcnhat * Dtset%pawnhatxc
488 0 : cplex = 1; ider = 2 * nhatgrdim; izero = 0
489 0 : if (nhatgrdim > 0) then
490 0 : ABI_MALLOC(ks_nhatgr,(nfftf,Dtset%nspden,3*nhatgrdim))
491 : end if
492 0 : if (nhatgrdim == 0) then
493 0 : ABI_MALLOC(ks_nhatgr,(0,0,0))
494 : end if
495 :
496 : call pawmknhat(compch_fft,cplex,ider,idir0,ipert0,izero,Cryst%gprimd,&
497 : Cryst%natom,Cryst%natom,nfftf,ngfftf,nhatgrdim,Dtset%nspden,Cryst%ntypat,Pawang,&
498 : Pawfgrtab,ks_nhatgr,ks_nhat,KS_Pawrhoij,KS_Pawrhoij,Pawtab,k0,Cryst%rprimd,&
499 0 : Cryst%ucvol,dtset%usewvl,Cryst%xred)
500 :
501 : ! === Evaluate onsite energies, potentials, densities ===
502 : ! Initialize variables/arrays related to the PAW spheres.
503 : ! Initialize also lmselect (index of non-zero LM-moments of densities).
504 0 : ABI_MALLOC(KS_paw_ij, (Cryst%natom))
505 0 : has_dijso = Dtset%pawspnorb; has_dijU = merge(0, 1, Dtset%usepawu == 0)
506 :
507 0 : call paw_ij_nullify(KS_paw_ij)
508 : call paw_ij_init(KS_paw_ij,cplex,Dtset%nspinor,Dtset%nsppol,&
509 : Dtset%nspden,Dtset%pawspnorb,Cryst%natom,Cryst%ntypat,Cryst%typat,Pawtab,&
510 : has_dij=1,has_dijhartree=1,has_dijhat=1,has_dijxc=1,has_dijxc_hat=1,has_dijxc_val=1,&
511 : has_dijso=has_dijso,has_dijU=has_dijU,has_exexch_pot=1,has_pawu_occ=1, &
512 0 : has_dijfock=dtset%usefock)
513 :
514 0 : nkxc1 = 0
515 0 : ABI_MALLOC(KS_paw_an, (Cryst%natom))
516 0 : call paw_an_nullify(KS_paw_an)
517 0 : usevxctau_paw=pawxc_get_usekden(dtset%ixc)
518 : call paw_an_init(KS_paw_an,Cryst%natom,Cryst%ntypat,nkxc1,0,Dtset%nspden,&
519 0 : cplex,Dtset%pawxcdev,Cryst%typat,Pawang,Pawtab,has_vxc=1,has_vxcval=1,has_vxctau=usevxctau_paw)
520 :
521 : ! Calculate onsite vxc with and without core charge.
522 0 : nzlmopt=-1; option=0; compch_sph=greatest_real
523 : call pawdenpot(compch_sph,el_temp,Cryst%gprimd,ipert0,Dtset%ixc,Cryst%natom,Cryst%natom,Dtset%nspden,&
524 : Cryst%ntypat,Dtset%nucdipmom,nzlmopt,option,KS_Paw_an,KS_Paw_an,KS_energies%paw,KS_paw_ij,&
525 : Pawang,Dtset%pawprtvol,Pawrad,KS_Pawrhoij,Dtset%pawspnorb,&
526 : Pawtab,Dtset%pawxcdev,Dtset%spnorbscl,Dtset%xclevel,Dtset%xc_denpos,Dtset%xc_taupos,&
527 0 : Cryst%xred,Cryst%ucvol,Psps%znuclpsp,Dtset%spinaxis)
528 :
529 : else
530 0 : ABI_MALLOC(ks_nhatgr, (0, 0, 0))
531 0 : ABI_MALLOC(ks_paw_ij, (0))
532 0 : ABI_MALLOC(ks_paw_an, (0))
533 : end if ! PAW
534 :
535 : !call test_charge(nfftf,ks_ebands%nelect,Dtset%nspden,ks_rhor,Cryst%ucvol,&
536 : ! Dtset%usepaw,usexcnhat,Pawfgr%usefinegrid,compch_sph,compch_fft,drude_plsmf)
537 :
538 : ! For PAW, add the compensation charge on the FFT mesh, then get rho(G).
539 : ! NB: ks_nhat is already included in the density stored on file so we don't need to add it.to ks_rhor
540 : !if (dtset%usepaw==1) ks_rhor = ks_rhor + ks_nhat
541 :
542 : ! TODO: Overloaded interface with units or just change the API to accept units
543 0 : call prtrhomxmn(std_out, mpi_enreg_seq, nfftf, ngfftf, dtset%nspden, 1, ks_rhor, ucvol=cryst%ucvol)
544 0 : call prtrhomxmn(ab_out, mpi_enreg_seq, nfftf, ngfftf, dtset%nspden, 1, ks_rhor, ucvol=cryst%ucvol)
545 :
546 0 : if (dtset%usekden==1) then
547 0 : call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,ks_taur,optrhor=1,ucvol=cryst%ucvol)
548 0 : call prtrhomxmn(ab_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,ks_taur,optrhor=1,ucvol=cryst%ucvol)
549 : end if
550 :
551 : ! FFT n(r) --> n(g)
552 0 : ABI_MALLOC(ks_rhog, (2, nfftf))
553 0 : call fourdp(cplex1, ks_rhog, ks_rhor(:, 1), -1, mpi_enreg_seq, nfftf, 1, ngfftf, 0)
554 :
555 : ! Compute structure factor phases and large sphere cutoff
556 0 : ABI_MALLOC(ph1d, (2, 3 * (2 * Dtset%mgfft + 1) * Cryst%natom))
557 0 : ABI_MALLOC(ph1df, (2, 3 * (2 * mgfftf + 1) * Cryst%natom))
558 :
559 0 : call getph(cryst%atindx, cryst%natom, ngfftc(1), ngfftc(2), ngfftc(3), ph1d, cryst%xred)
560 :
561 0 : if (psps%usepaw == 1 .and. pawfgr%usefinegrid == 1) then
562 0 : call getph(cryst%atindx, cryst%natom, ngfftf(1), ngfftf(2), ngfftf(3), ph1df, cryst%xred)
563 : else
564 0 : ph1df(:,:)=ph1d(:,:)
565 : end if
566 :
567 : ! The following steps have been gathered in the setvtr routine:
568 : ! - get Ewald energy and Ewald forces
569 : ! - compute local ionic pseudopotential vpsp
570 : ! - eventually compute 3D core electron density xccc3d
571 : ! - eventually compute vxc and vhartr
572 : ! - set up ks_vtrial
573 : !
574 : !*******************************************************************
575 : !**** NOTE THAT HERE Vxc CONTAINS THE CORE-DENSITY CONTRIBUTION ****
576 : !*******************************************************************
577 :
578 0 : ngrvdw = 0
579 0 : ABI_MALLOC(grvdw, (3, ngrvdw))
580 0 : ABI_MALLOC(grchempottn, (3, cryst%natom))
581 0 : ABI_MALLOC(grewtn, (3, cryst%natom))
582 0 : nkxc = 0
583 0 : if (dtset%nspden == 1) nkxc = 2
584 0 : if (dtset%nspden >= 2) nkxc = 3 ! check GGA and spinor, quite a messy part!!!
585 : ! In case of MGGA, fxc and kxc are not available and we dont need them (for now ...)
586 0 : if (dtset%ixc < 0 .and. libxc_functionals_ismgga()) nkxc = 0
587 0 : if (nkxc /= 0) then
588 0 : ABI_MALLOC(kxc, (nfftf, nkxc))
589 : end if
590 :
591 0 : n3xccc = 0; if (psps%n1xccc /= 0) n3xccc = nfftf
592 0 : ABI_MALLOC(xccc3d, (n3xccc))
593 0 : ABI_MALLOC(ks_vhartr, (nfftf))
594 0 : ABI_MALLOC(ks_vtrial, (nfftf, dtset%nspden))
595 0 : ABI_MALLOC(vpsp, (nfftf))
596 0 : ABI_MALLOC(ks_vxc, (nfftf, dtset%nspden))
597 :
598 : ! TODO: I don't think direct diago can be used with mega-GGA due to the functional derivative wrt KS states.
599 : ! TB-BK should be OK though.
600 :
601 : !usevxctau=merge(1,0,xc_need_kden(dtset%ixc))
602 : !ABI_MALLOC(ks_vxctau, (nfftf, dtset%nspden * usevxctau))
603 : !ABI_MALLOC(xcctau3d, (n3xccc * dtset%usekden))
604 : !ABI_FREE(ks_vxctau)
605 : !ABI_FREE(xcctau3d)
606 :
607 0 : optene = 4; moved_atm_inside = 0; moved_rhor = 0; istep = 1
608 :
609 : call setvtr(Cryst%atindx1,Dtset,KS_energies,cryst%gmet,cryst%gprimd,grchempottn,grewtn,grvdw,gsqcutf_eff,&
610 : istep,kxc,mgfftf,moved_atm_inside,moved_rhor,MPI_enreg_seq,&
611 : Cryst%nattyp,nfftf,ngfftf,ngrvdw,ks_nhat,ks_nhatgr,nhatgrdim,nkxc,Cryst%ntypat,Psps%n1xccc,n3xccc,&
612 : optene,Pawang,Pawrad,KS_pawrhoij,Pawtab,ph1df,Psps,ks_rhog,ks_rhor,cryst%rmet,cryst%rprimd,strsxc,&
613 : Cryst%ucvol,usexcnhat,ks_vhartr,vpsp,ks_vtrial,ks_vxc,vxcavg,Wvl,xccc3d,Cryst%xred, &
614 0 : taur=ks_taur) !xcctau3d=xcctau3d, vxctau=ks_vxctau)
615 :
616 0 : ABI_FREE(grvdw)
617 0 : ABI_FREE(grchempottn)
618 0 : ABI_FREE(grewtn)
619 :
620 : !============================
621 : !==== Compute KS PAW Dij ====
622 : !============================
623 0 : if (dtset%usepaw == 1) then
624 0 : call timab(561,1,tsec)
625 :
626 : ! Calculate the unsymmetrized Dij.
627 : call pawdij(cplex,Dtset%enunit,Cryst%gprimd,ipert0,&
628 : Cryst%natom,Cryst%natom,nfftf,ngfftf(1)*ngfftf(2)*ngfftf(3),&
629 : Dtset%nspden,Cryst%ntypat,KS_paw_an,KS_paw_ij,Pawang,Pawfgrtab,&
630 : Dtset%pawprtvol,Pawrad,KS_Pawrhoij,Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,&
631 : k0,Dtset%spnorbscl,Cryst%ucvol,dtset%cellcharge(1),ks_vtrial,&
632 : ks_vxc,Cryst%xred,Dtset%znucl,&
633 0 : nucdipmom=Dtset%nucdipmom,spinaxis=Dtset%spinaxis)
634 :
635 : ! Symmetrize KS Dij
636 : call symdij_all(Cryst%gprimd,Cryst%indsym,ipert0,&
637 : Cryst%natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,KS_paw_ij,Pawang,&
638 0 : Dtset%pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec)
639 :
640 : ! Output the pseudopotential strengths Dij and the augmentation occupancies Rhoij.
641 0 : call pawprt(Dtset,Cryst%natom,KS_paw_ij,KS_Pawrhoij,Pawtab)
642 0 : call timab(561,2,tsec)
643 : end if
644 :
645 0 : call pstat_proc%print(_PSTAT_ARGS_)
646 0 : call cwtime_report(" prepare gwr_driver_init", cpu, wall, gflops)
647 :
648 0 : if (string_in(dtset%gwr_task, "HDIAGO, HDIAGO_FULL, CC4S, CC4S_FULL")) then
649 : ! ==========================================
650 : ! Direct diagonalization of the Hamiltonian
651 : ! ==========================================
652 0 : ABI_MALLOC(nband_iks, (dtset%nkpt, dtset%nsppol))
653 0 : ABI_MALLOC(npwarr_ik, (dtset%nkpt))
654 0 : ABI_MALLOC(istwfk_ik, (dtset%nkpt))
655 0 : istwfk_ik = 1
656 :
657 : ! Compute npw_k from ecut so that we can update the header and redefine %mpw
658 0 : do ik_ibz=1,dtset%nkpt
659 : !if (dtset%istwfk(ik_ibz) == 2) istwfk_ik(ik_ibz) = 2 ! TODO: istwkf 2 is not yet supported.
660 0 : call get_kg(dtset%kptns(:,ik_ibz), istwfk_ik(ik_ibz), dtset%ecut, cryst%gmet, npwarr_ik(ik_ibz), gvec_)
661 0 : ABI_FREE(gvec_)
662 : end do
663 0 : dtset%mpw = maxval(npwarr_ik)
664 :
665 : ! CC4S does not need to output the WFK file.
666 0 : write_wfk = string_in(dtset%gwr_task, "HDIAGO, HDIAGO_FULL")
667 :
668 : ! Use input nband or min of npwarr_ik to set the number of bands.
669 0 : if (string_in(dtset%gwr_task, "HDIAGO, CC4S")) nband_iks(:,:) = maxval(dtset%nband)
670 0 : if (string_in(dtset%gwr_task, "HDIAGO_FULL, CC4S_FULL")) nband_iks(:,:) = minval(npwarr_ik)
671 0 : cc4s_task = string_in(dtset%gwr_task, "CC4S, CC4S_FULL")
672 0 : if (cc4s_task) then
673 0 : ABI_CHECK_IEQ(dtset%nkpt, 1, "CC4S interface does not support more than one k-point.")
674 : end if
675 :
676 : ! Build header with new npwarr and nband.
677 0 : call owfk_ebands%from_dtset(dtset, npwarr_ik, nband=nband_iks)
678 0 : owfk_ebands%eig = zero
679 0 : owfk_ebands%istwfk = istwfk_ik
680 : !print *, "owfk_ebands%npwarr:", owfk_ebands%npwarr; stop
681 0 : call owfk_hdr%init(owfk_ebands, codvsn, dtset, pawtab, 0, psps, wvl%descr)
682 :
683 : ! Change the value of istwfk taken from dtset and set the Fermie level from gs_fermie.
684 0 : ABI_REMALLOC(owfk_hdr%istwfk, (dtset%nkpt))
685 0 : owfk_hdr%istwfk(:) = istwfk_ik
686 0 : owfk_hdr%fermie = gs_fermie
687 :
688 : ! Build MPI pools to distribute (kpt, spin).
689 : ! Try to get rectangular grids in each pool to improve efficiency in slk diago.
690 0 : rectangular = .True.; if (dtset%nkpt == 1) rectangular = .False.
691 0 : with_pools = .True.
692 0 : if (any(dtset%gwr_np_kgts /= 0) .and. (nprocs == dtset%gwr_np_kgts(2))) then
693 0 : with_pools = .False.; call wrtout(std_out, " Deactivating pools in HDIAGO")
694 : end if
695 : !with_pools = .False.
696 0 : call diago_pool%from_dims(dtset%nkpt, dtset%nsppol, comm, with_pools, rectangular=rectangular)
697 0 : diago_info = zero
698 :
699 : ! TODO: Build hyb descriptor with hybrid orbitals from WFK file.
700 0 : cc4s_from_wfk = .False.
701 : if (dtset%usefock == 1 .and. cc4s_from_wfk) then
702 : call hyb%from_wfk_file(cryst, dtfil, dtset, psps, pawtab, ngfftc, diago_pool, comm)
703 : end if
704 :
705 0 : if (write_wfk) then
706 : ! Master writes the Abinit header and the Fortran record markers
707 0 : out_path = dtfil%fnameabo_wfk; if (dtset%iomode == IO_MODE_ETSF) out_path = nctk_ncify(out_path)
708 0 : iomode__ = iomode_from_fname(out_path)
709 0 : call wrtout(std_out, sjoin(" Writing wavefunctions to file:", out_path))
710 0 : if (my_rank == master) then
711 0 : call owfk%open_write(owfk_hdr, out_path, 0, iomode__, get_unit(), xmpi_comm_self, write_hdr=.True., write_frm=.True.)
712 0 : call owfk%close()
713 : end if
714 0 : call xmpi_barrier(comm)
715 : end if
716 :
717 : ! Build H_k(g,g') for each k-point and spin and diagonalize the Hamiltonian with Scalapack/ELPA.
718 0 : do spin=1,dtset%nsppol
719 0 : do ik_ibz=1,dtset%nkpt
720 0 : if (.not. diago_pool%treats(ik_ibz, spin)) cycle
721 0 : nband_k = nband_iks(ik_ibz, spin)
722 :
723 0 : if (dtset%usefock == 1 .and. .not. cc4s_from_wfk) then
724 0 : call wrtout(units, "Reading ugb datatype from WFK file")
725 0 : ABI_CHECK(.not. string_in(dtset%gwr_task, "CC4S_FULL"), "CC4S_FULL cannot be used with Fock, please specify nband")
726 : call ugb%from_wfk_file(ik_ibz, spin, istwfk_ik(ik_ibz), dtset%kptns(:,ik_ibz), nband_k, dtset, &
727 0 : dtfil, cryst, eig_k, diago_pool%comm%value)
728 :
729 : else
730 0 : call cwtime(diago_cpu, diago_wall, diago_gflops, "start")
731 : call ugb%from_diago(spin, istwfk_ik(ik_ibz), dtset%kptns(:,ik_ibz), dtset%ecut, gs_fermie, nband_k, ngfftc, nfftf, &
732 0 : dtset, pawtab, pawfgr, ks_paw_ij, cryst, psps, ks_vtrial, eig_k, hyb, diago_pool%comm%value)
733 0 : call cwtime(diago_cpu, diago_wall, diago_gflops, "stop")
734 :
735 0 : if (diago_pool%comm%me == 0) diago_info(1, ik_ibz, spin) = diago_wall
736 0 : call cwtime(diago_cpu, diago_wall, diago_gflops, "start")
737 : end if
738 :
739 0 : owfk_ebands%eig(1:nband_k, ik_ibz, spin) = eig_k(1:nband_k)
740 :
741 0 : if (write_wfk) then
742 : ! occupancies are set to zero. Client code is responsible for recomputing occ and fermie when reading this WFK.
743 0 : call wrtout(std_out, sjoin(" Writing kpoint ik_ibz:", itoa(ik_ibz), ", spin: ", itoa(spin), " to disk ..."))
744 0 : ABI_CALLOC(occ_k, (nband_k))
745 0 : color = merge(1, 0, ugb%my_nband > 0)
746 0 : call xmpi_comm_split(diago_pool%comm%value, color, diago_pool%comm%me, io_comm, ierr)
747 : call wrtout(std_out, sjoin(" Number of MPI procs involved in IO operation:", &
748 0 : itoa(xmpi_comm_size(io_comm)), "/", itoa(diago_pool%comm%nproc)))
749 :
750 0 : if (ugb%my_nband > 0) then
751 0 : ABI_CHECK(all(shape(ugb%cg_k) == [2, ugb%npwsp, ugb%my_nband]), "Wrong shape")
752 0 : ABI_CHECK_IEQ(ugb%npw_k, owfk_hdr%npwarr(ik_ibz), "Wrong npw_k")
753 0 : call c_f_pointer(c_loc(ugb%cg_k), cg_k_ptr, shape=[2, ugb%npwsp * ugb%my_nband])
754 :
755 : ! Reopen file inside io_comm.
756 0 : call owfk%open_write(owfk_hdr, out_path, 0, iomode__, get_unit(), io_comm, write_hdr=.False., write_frm=.False.)
757 : !call wrtout(std_out, "After owfk%open_write")
758 :
759 : ! Write my bands.
760 : !sc_mode = merge(xmpio_single, xmpio_collective, ugb%has_idle_procs)
761 : !sc_mode = xmpio_collective
762 0 : sc_mode = xmpio_single
763 : call owfk%write_band_block([ugb%my_bstart, ugb%my_bstop], ik_ibz, spin, sc_mode, &
764 0 : kg_k=ugb%kg_k, cg_k=cg_k_ptr, eig_k=owfk_ebands%eig(:, ik_ibz, spin), occ_k=occ_k)
765 0 : call owfk%close()
766 : end if
767 0 : call xmpi_comm_free(io_comm)
768 0 : ABI_FREE(occ_k)
769 : end if
770 :
771 0 : call cwtime(diago_cpu, diago_wall, diago_gflops, "stop")
772 0 : if (diago_pool%comm%me == 0) diago_info(2:3, ik_ibz, spin) = [diago_wall, dble(diago_pool%comm%nproc)]
773 :
774 : ! Compute and write matrix elements required by CC4S
775 0 : if (cc4s_task) call cc4s_gamma(spin, ik_ibz, dtset, dtfil, cryst, owfk_ebands, psps, pawtab, paw_pwff, ugb)
776 :
777 0 : ABI_FREE(eig_k)
778 0 : call ugb%free()
779 : end do ! ik_ibz
780 : end do ! spin
781 0 : call wrtout(std_out, " Direct diago completed by this MPI pool. Other pools might take more time if k != 0")
782 :
783 0 : call xmpi_sum_master(diago_info, master, comm, ierr)
784 0 : if (my_rank == master) then
785 0 : do spin=1,dtset%nsppol
786 0 : do ik_ibz=1,dtset%nkpt
787 0 : associate (info => diago_info(:, ik_ibz, spin))
788 0 : write(std_out, "(2(a,i0),5a,i0)") " ik_ibz: ", ik_ibz, ", spin: ", spin, &
789 0 : ", diago_wall: ", trim(sec2str(info(1))), ", io_wall: ", trim(sec2str(info(2))), ", nprocs: ", int(info(3))
790 : end associate
791 : end do
792 : end do
793 : end if
794 :
795 : ! Collect eigenvalues for the different k-points/spins.
796 0 : do spin=1,dtset%nsppol
797 0 : do ik_ibz=1,dtset%nkpt
798 0 : if (diago_pool%treats(ik_ibz, spin) .and. diago_pool%comm%me /= 0) owfk_ebands%eig(:, ik_ibz, spin) = zero
799 : end do
800 : end do
801 0 : call xmpi_sum(owfk_ebands%eig, comm, ierr)
802 :
803 0 : call owfk_ebands%update_occ(dtset%spinmagntarget, prtvol=dtset%prtvol, fermie_to_zero=.False.)
804 :
805 0 : if (my_rank == master) then
806 0 : if (write_wfk .and. iomode__ == IO_MODE_ETSF) then
807 0 : NCF_CHECK(owfk_ebands%ncwrite_path(cryst, out_path))
808 : end if
809 0 : call owfk_ebands%print_gaps(units, header="KS gaps after direct diagonalization")
810 0 : if (cc4s_task) call cc4s_write_eigens(owfk_ebands, dtfil)
811 : end if
812 :
813 0 : ABI_FREE(npwarr_ik)
814 0 : ABI_FREE(istwfk_ik)
815 0 : ABI_FREE(nband_iks)
816 0 : call owfk_hdr%free(); call owfk_ebands%free(); call hyb%free(); call diago_pool%free()
817 :
818 0 : else if (string_in(dtset%gwr_task, "CC4S_FROM_WFK")) then
819 : ! Read orbitals from an external WFK file and produce output files for CC4S.
820 :
821 : ! Construct crystal and ks_ebands from the GS WFK file.
822 0 : tmp_ebands = wfk_read_ebands(wfk_path, comm, out_hdr=wfk_hdr)
823 0 : ks_ebands = tmp_ebands%chop(1, maxval(dtset%nband))
824 :
825 : ! Make sure that ef is inside the gap if semiconductor.
826 : !call ks_ebands%update_occ(dtset%spinmagntarget, prtvol=dtset%prtvol, fermie_to_zero=.True.)
827 :
828 0 : call tmp_ebands%free()
829 0 : call wfk_hdr%vs_dtset(dtset)
830 :
831 0 : wfk_cryst = wfk_hdr%get_crystal()
832 0 : if (cryst%compare(wfk_cryst, header=" Comparing input crystal with WFK crystal") /= 0) then
833 0 : ABI_ERROR("Crystal structure from input and from WFK file do not agree! Check messages above!")
834 : end if
835 0 : call wfk_cryst%free()
836 :
837 0 : if (my_rank == master) call cc4s_write_eigens(ks_ebands, dtfil)
838 :
839 : ! Build MPI pools to distribute (kpt, spin).
840 0 : with_pools = .True.
841 0 : if (any(dtset%gwr_np_kgts /= 0) .and. (nprocs == dtset%gwr_np_kgts(2))) then
842 0 : with_pools = .False.; call wrtout(std_out, " Deactivating pools in HDIAGO")
843 : end if
844 0 : call diago_pool%from_dims(dtset%nkpt, dtset%nsppol, comm, with_pools, rectangular=.False.)
845 0 : ABI_CHECK_IEQ(dtset%nkpt, 1, "Only Gamma-point sampling is supported")
846 0 : ABI_CHECK_IEQ(dtset%nsppol, 1, "Only spin-unpolarized calculations are supported")
847 :
848 0 : do spin=1,dtset%nsppol
849 0 : do ik_ibz=1,dtset%nkpt
850 0 : if (.not. diago_pool%treats(ik_ibz, spin)) cycle
851 0 : nband_k = dtset%nband(ik_ibz + (spin-1)*dtset%nkpt)
852 : call ugb%from_wfk_file(ik_ibz, spin, dtset%istwfk(ik_ibz), dtset%kptns(:,ik_ibz), nband_k, dtset, &
853 0 : dtfil, cryst, eig_k, diago_pool%comm%value)
854 0 : ABI_FREE(eig_k)
855 0 : call cc4s_gamma(spin, ik_ibz, dtset, dtfil, cryst, ks_ebands, psps, pawtab, paw_pwff, ugb)
856 0 : call ugb%free()
857 : end do
858 : end do
859 0 : call wfk_hdr%free(); call ks_ebands%free(); call diago_pool%free()
860 :
861 : else
862 : ! ====================================================
863 : ! === This is the real GWR stuff once all is ready ===
864 : ! ====================================================
865 0 : ABI_CHECK(dtset%usepaw == 0, "PAW in GWR not yet implemented.")
866 0 : read_wfk = .True.
867 : if (read_wfk) then
868 0 : if (my_rank == master) then
869 0 : if (nctk_try_fort_or_ncfile(wfk_path, msg) /= 0) then
870 0 : ABI_ERROR(sjoin("Cannot find GS WFK file:", wfk_path, ". Error:", msg))
871 : end if
872 0 : call wrtout(units, sjoin("- Reading GS states from WFK file:", wfk_path))
873 : end if
874 :
875 : ! Broadcast filenames (needed because they might have been changed if we are using netcdf files)
876 0 : call xmpi_bcast(wfk_path, master, comm, ierr)
877 :
878 : ! Construct crystal and ks_ebands from the GS WFK file.
879 0 : tmp_ebands = wfk_read_ebands(wfk_path, comm, out_hdr=wfk_hdr)
880 0 : ks_ebands = tmp_ebands%chop(1, maxval(dtset%nband))
881 0 : call tmp_ebands%free()
882 0 : call wfk_hdr%vs_dtset(dtset)
883 :
884 0 : wfk_cryst = wfk_hdr%get_crystal()
885 0 : if (cryst%compare(wfk_cryst, header=" Comparing input crystal with WFK crystal") /= 0) then
886 0 : ABI_ERROR("Crystal structure from input and from WFK file do not agree! Check messages above!")
887 : end if
888 : !call wfk_cryst%print(header="crystal structure from WFK file")
889 0 : call wfk_cryst%free()
890 :
891 : ! Make sure that ef is inside the gap if semiconductor.
892 0 : call ks_ebands%update_occ(dtset%spinmagntarget, prtvol=dtset%prtvol, fermie_to_zero=.True.)
893 :
894 : ! Here we change the GS bands (Fermi level, scissors operator ...)
895 : ! All the modifications to ebands should be done here.
896 : !call ephtk_update_ebands(dtset, dtfil%filqpdatain, ks_ebands, "Ground state energies", comm)
897 : end if
898 0 : call pstat_proc%print(_PSTAT_ARGS_)
899 :
900 0 : call gwr%init(dtset, dtfil, cryst, psps, pawtab, ks_ebands, mpi_enreg_seq, comm)
901 0 : if (gwr%idle_proc) goto 100
902 :
903 : !=== Calculate Vxc(b1,b2,k,s)=<b1,k,s|v_{xc}|b2,k,s> for all the states included in GW ===
904 : ! * This part is parallelized within wfd%comm since each node has all GW wavefunctions.
905 : ! * Note that vH matrix elements are calculated using the true uncutted interaction.
906 :
907 0 : rdm_update = dtset%gwr_task == "GAMMA_GW"
908 :
909 0 : call KS_mflags%reset()
910 0 : if (rdm_update) then
911 0 : KS_mflags%has_hbare=1
912 0 : KS_mflags%has_kinetic=1
913 : end if
914 0 : KS_mflags%has_vhartree=1
915 0 : KS_mflags%has_vxc =1
916 0 : KS_mflags%has_vxcval =1
917 0 : if (Dtset%usepawu /= 0 ) KS_mflags%has_vu = 1
918 0 : if (Dtset%useexexch /= 0) KS_mflags%has_lexexch = 1
919 0 : if (Dtset%usepaw==1 .and. Dtset%gw_sigxcore == 1) KS_mflags%has_sxcore = 1
920 : ! off-diagonal elements only for SC on wavefunctions.
921 0 : KS_mflags%only_diago = 1
922 0 : if (rdm_update) KS_mflags%only_diago = 0
923 :
924 : ! Load wavefunctions for Sigma_nk in gwr%kcalc_wfd.
925 0 : call gwr%load_kcalc_wfd(wfk_path, tmp_kstab)
926 :
927 : ! Compute gwr%ks_me matrix elements.
928 0 : if (.not. string_in(dtset%gwr_task, "RPA_ENERGY")) then
929 : ! FIXME: This routine allocates (nband, nband) matrices and should be rewritten!
930 : call calc_vhxc_me(gwr%kcalc_wfd, ks_mflags, gwr%ks_me, cryst, dtset, nfftf, ngfftf, &
931 : ks_vtrial, ks_vhartr, ks_vxc, psps, pawtab, ks_paw_an, pawang, pawfgrtab, ks_paw_ij, dijexc_core, &
932 0 : ks_rhor, usexcnhat, ks_nhat, ks_nhatgr, nhatgrdim, tmp_kstab, taur=ks_taur)
933 0 : if (my_rank == master) call gwr%ks_me%print(header="KS matrix elements", unit=std_out)
934 : end if
935 :
936 0 : ABI_FREE(tmp_kstab)
937 :
938 : if (read_wfk) then
939 : ! Read wavefunctions from WFK file.
940 0 : call gwr%read_ugb_from_wfk(wfk_path)
941 : else
942 : ! Diagonalize H on the fly and
943 : !call gwr%get_ugb_from_vtrial(ngfftf, ks_vtrial)
944 : !gwr%wfk_hdr = ?
945 : end if
946 :
947 : ! Now call high-level routines depending on gwr_task.
948 0 : select case (dtset%gwr_task)
949 : case ("RPA_ENERGY")
950 0 : call gwr%rpa_energy()
951 : case ("GAMMA_GW")
952 0 : call gwr%gamma_gw(nfftf, ngfftf, vpsp)
953 : case ("CHI0")
954 0 : call gwr%run_chi0()
955 : case ("G0W0")
956 0 : call gwr%run_g0w0()
957 : case ("G0V")
958 0 : call gwr%build_sigxme(compute_qp=.True.)
959 : case ("EGEW", "EGW0", "G0EW")
960 0 : call gwr%run_energy_scf()
961 : case default
962 0 : ABI_ERROR(sjoin("Invalid value for gwr_task:", dtset%gwr_task))
963 : end select
964 : end if
965 :
966 : !=====================
967 : !==== Free memory ====
968 : !=====================
969 0 : 100 call xmpi_barrier(comm)
970 0 : ABI_FREE(ks_nhat)
971 0 : ABI_FREE(ks_nhatgr)
972 0 : ABI_FREE(dijexc_core)
973 0 : call pawfgr_destroy(pawfgr)
974 :
975 0 : if (dtset%usepaw == 1) then
976 : ! Deallocation for PAW.
977 0 : call pawrhoij_free(ks_pawrhoij)
978 0 : ABI_FREE(ks_pawrhoij)
979 0 : call pawfgrtab_free(pawfgrtab)
980 0 : call paw_ij_free(ks_paw_ij)
981 0 : ABI_FREE(ks_paw_ij)
982 0 : call paw_an_free(ks_paw_an)
983 0 : call pawpwff_free(Paw_pwff)
984 : end if
985 :
986 0 : ABI_FREE(ph1d)
987 0 : ABI_FREE(ph1df)
988 0 : ABI_FREE(ks_rhor)
989 0 : ABI_FREE(ks_rhog)
990 0 : ABI_FREE(ks_taur)
991 0 : ABI_FREE(kxc)
992 0 : ABI_FREE(xccc3d)
993 0 : ABI_FREE(ks_vhartr)
994 0 : ABI_FREE(ks_vtrial)
995 0 : ABI_FREE(vpsp)
996 0 : ABI_FREE(ks_vxc)
997 : ! PAW stuff
998 0 : ABI_SFREE(paw_pwff)
999 0 : ABI_SFREE(pawfgrtab)
1000 0 : ABI_SFREE(ks_paw_an)
1001 :
1002 0 : call cryst%free(); call wfk_hdr%free(); call ks_ebands%free(); call destroy_mpi_enreg(mpi_enreg_seq)
1003 : #if defined FC_NVHPC
1004 : call wrtout(units, "- Cannot deallocate gwr datatype if FC_NVHPC, DO NOT USE DATASETS!")
1005 : #else
1006 0 : call gwr%free()
1007 : #endif
1008 :
1009 0 : call pstat_proc%print(_PSTAT_ARGS_)
1010 :
1011 0 : end subroutine gwr_driver
1012 : !!***
1013 :
1014 : !!****f* m_gwr_driver/cc4s_write_eigens
1015 : !! NAME
1016 : !! cc4s_write_eigens
1017 : !!
1018 : !! FUNCTION
1019 : !! Write eigenvalues in CC4S format. Only master proc should call this routine.
1020 : !!
1021 : !! INPUTS
1022 :
1023 0 : subroutine cc4s_write_eigens(ebands, dtfil)
1024 :
1025 : !Arguments ------------------------------------
1026 : type(ebands_t),intent(in) :: ebands
1027 : type(datafiles_type),intent(in) :: dtfil
1028 :
1029 : !Local variables-------------------------------
1030 : integer :: unt, ik_ibz, spin, band, gap_err
1031 : real(dp) :: my_fermie
1032 : character(len=500) :: msg
1033 : character(len=fnlen) :: filepath
1034 0 : type(gaps_t) :: ks_gaps
1035 : ! *************************************************************************
1036 :
1037 0 : ks_gaps = ebands%get_gaps(gap_err)
1038 0 : call ks_gaps%print([std_out], header="Kohn-Sham gaps and band edges from IBZ mesh")
1039 0 : if (any(ks_gaps%ierr /= 0)) then
1040 0 : ABI_ERROR("Cannot compute gaps!")
1041 : end if
1042 :
1043 : ! Make sure fermi level is within the gap. I know, the case nsppol = 2 is treated in a dirty way!
1044 0 : my_fermie = zero
1045 0 : do spin=1,ebands%nsppol
1046 0 : my_fermie = ks_gaps%vb_max(spin) + (ks_gaps%cb_min(spin) - ks_gaps%vb_max(spin)) / two
1047 : end do
1048 0 : my_fermie = my_fermie / ebands%nsppol
1049 0 : call ks_gaps%free()
1050 :
1051 : ! See https://manuals.cc4s.org/user-manual/objects/EigenEnergies.html
1052 0 : filepath = trim(dtfil%filnam_ds(4))//'_EigenEnergies.yaml'
1053 0 : write(ab_out, "(3a)")ch10," Writing Eigenenergies metadata to file: ", trim(filepath)
1054 0 : if (open_file(filepath, msg, newunit=unt, access="stream", form="formatted", status="replace", action="write") /= 0) then
1055 0 : ABI_ERROR(msg)
1056 : end if
1057 0 : write(unt,'(A)') 'version: 100'
1058 0 : write(unt,'(A)') 'type: Tensor'
1059 0 : write(unt,'(A)') 'scalarType: Real64'
1060 0 : write(unt,'(A)') 'dimensions:'
1061 0 : write(unt,'(A,I0)') '- length: ',ebands%mband * ebands%nkpt * ebands%nsppol
1062 0 : write(unt,'(A)') ' type: State'
1063 0 : write(unt,'(A)') 'elements:'
1064 0 : write(unt,'(A)') ' type: TextFile'
1065 0 : write(unt,'(A)') 'unit: 1.0 # Hartree units'
1066 0 : write(unt,'(A)') 'metaData:'
1067 0 : write(unt,'(A,E22.15)') ' fermiEnergy: ',my_fermie
1068 0 : write(unt,'(A)') ' energies:'
1069 :
1070 0 : do spin=1,ebands%nsppol
1071 0 : do ik_ibz=1,ebands%nkpt
1072 0 : do band=1,ebands%nband(ik_ibz + (spin-1)*ebands%nkpt)
1073 0 : write(unt,"(a,e22.15)") ' - ',ebands%eig(band,ik_ibz,spin)
1074 : end do
1075 : end do
1076 : end do
1077 0 : close(unt)
1078 :
1079 0 : filepath = trim(dtfil%filnam_ds(4))//'_EigenEnergies.elements'
1080 0 : write(ab_out, "(3a)")ch10," Writing Eigenenergies to file: ", trim(filepath)
1081 0 : if (open_file(filepath, msg, newunit=unt, access="stream", form="formatted", status="replace", action="write") /= 0) then
1082 0 : ABI_ERROR(msg)
1083 : end if
1084 :
1085 0 : do spin=1,ebands%nsppol
1086 0 : do ik_ibz=1,ebands%nkpt
1087 0 : do band=1,ebands%nband(ik_ibz + (spin-1)*ebands%nkpt)
1088 0 : write(unt,"(e22.15)") ebands%eig(band,ik_ibz,spin)
1089 : end do
1090 : end do
1091 : end do
1092 :
1093 0 : close(unt)
1094 :
1095 0 : end subroutine cc4s_write_eigens
1096 : !!***
1097 :
1098 : !!****f* m_gwr_driver/cc4s_gamma
1099 : !! NAME
1100 : !! cc4s_gamma
1101 : !!
1102 : !! FUNCTION
1103 : !! Interface with CC4S code.
1104 : !! Compute <b1,k|e^{-iGr}|b2,k> matrix elements and store them to disk
1105 : !!
1106 : !! INPUTS
1107 :
1108 0 : subroutine cc4s_gamma(spin, ik_ibz, dtset, dtfil, cryst, ebands, psps, pawtab, paw_pwff, ugb)
1109 :
1110 : use m_numeric_tools, only : blocked_loop
1111 : use m_gwdefs, only : GW_Q0_DEFAULT
1112 : use m_fftcore, only : sphereboundary
1113 : use m_fft_mesh, only : setmesh
1114 : use m_fft, only : uplan_t
1115 : use m_vcoul, only : vcgen_t
1116 :
1117 : use m_sort, only : sort_gvecs
1118 : use m_pawpwij, only : pawpwij_t, pawpwij_init, pawpwij_free
1119 :
1120 : !Arguments ------------------------------------
1121 : integer,intent(in) :: spin, ik_ibz
1122 : type(dataset_type),intent(in) :: dtset
1123 : type(datafiles_type),intent(in) :: dtfil
1124 : type(crystal_t),intent(in) :: cryst
1125 : type(ebands_t),intent(in) :: ebands
1126 : type(Pseudopotential_type),intent(in) :: psps
1127 : type(pawpwff_t),intent(in) :: paw_pwff(dtset%ntypat*dtset%usepaw)
1128 : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
1129 : type(ugb_t),target,intent(in) :: ugb
1130 :
1131 : !Local variables-------------------------------
1132 : !scalars
1133 : integer,parameter :: mG0(3) = 0, master = 0
1134 : integer :: nproc, my_rank, my_ib2st, npw_k, nspinor, m_npw, npwvec, ig, mpierr, fh, comm, buf_size, ierr
1135 : integer :: band1, band1_start, batch1_size, n1dat, idat1, m_istwfk, iatom, dim_rtwg
1136 : integer :: band2, band2_start, batch2_size, n2dat, idat2, units(2), ii, unt, nqibz_, nqbz_, nkbz_, test_unt, M_
1137 : integer(XMPI_OFFSET_KIND) :: offset
1138 : real(dp) :: cpu, wall, gflops, qpt(3), qbz_(3,1), gcart(3), kpt(3), max_abs_err, abs_err, my_gw_qlwl(3), mem_mb
1139 : character(len=500) :: msg
1140 : character(len=fnlen) :: filepath, cvx_filepath
1141 : logical :: k_is_gamma
1142 : logical,parameter :: trust_no_one = .False.
1143 0 : type(uplan_t) :: uplan_1, uplan_2, uplan_m
1144 0 : type(vcgen_t) :: vcgen
1145 0 : integer :: u_ngfft(18), u_nfft, u_mgfft, enforce_sym, method, nlmn_atm(cryst%natom)
1146 0 : integer,pointer :: gvec_max(:,:)
1147 0 : integer,allocatable,target :: m_gvec(:,:), sorted_kg_k(:,:)
1148 0 : complex(dp),allocatable :: ug1_batch(:,:), ur1_batch(:,:), ur2_batch(:,:), ur12_batch(:,:), ug12_batch(:,:), cwork(:)
1149 0 : complex(gwp),allocatable :: sqrt_vc(:), paw_rhotwg(:)
1150 0 : type(pawpwij_t),allocatable :: pwij(:)
1151 0 : type(pawcprj_type),allocatable :: cprj1(:,:)
1152 : ! *************************************************************************
1153 :
1154 0 : call cwtime(cpu, wall, gflops, "start")
1155 :
1156 0 : comm = ugb%comm; nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
1157 0 : units = [std_out, ab_out]
1158 0 : npw_k = ugb%npw_k; nspinor = ugb%nspinor
1159 :
1160 0 : if (dtset%prtvol > 10) call ugb%print([std_out], dtset%prtvol, header="ugb for CC4S")
1161 :
1162 : ! m_gvec is the g-sphere for the oscillators M computed from ecuteps (half-sphere if wavefunctions have TR).
1163 : ! setmesh assumes g-vectors sorted by norm so use kin_sorted = True and sort ug%kg_k below.
1164 0 : kpt = dtset%kptns(:,ik_ibz); k_is_gamma = all(abs(kpt) < tol12)
1165 0 : m_istwfk = 1; if (ugb%istwf_k == 2) m_istwfk = 2
1166 0 : call get_kg(kpt, m_istwfk, dtset%ecuteps, cryst%gmet, m_npw, m_gvec, kin_sorted=.True.)
1167 :
1168 : ! Setup FFT mesh
1169 0 : u_ngfft = dtset%ngfft
1170 0 : method = 2
1171 0 : if (dtset%fftgw==00 .or. dtset%fftgw==01) method=0
1172 0 : if (dtset%fftgw==10 .or. dtset%fftgw==11) method=1
1173 0 : if (dtset%fftgw==20 .or. dtset%fftgw==21) method=2
1174 0 : if (dtset%fftgw==30 .or. dtset%fftgw==31) method=3
1175 0 : enforce_sym = mod(dtset%fftgw, 10)
1176 : ! Gamma only --> we don't need to rotate wavefunctions in the BZ
1177 0 : if (k_is_gamma) enforce_sym = 0
1178 :
1179 0 : call sort_gvecs(npw_k, kpt, cryst%gmet, ugb%kg_k, sorted_kg_k)
1180 :
1181 0 : npwvec = npw_k; gvec_max => sorted_kg_k
1182 0 : if (m_npw > npw_k) then
1183 0 : npwvec = m_npw; gvec_max => m_gvec
1184 : end if
1185 0 : call setmesh(cryst%gmet, gvec_max, u_ngfft, npwvec, m_npw, npw_k, u_nfft, method, mG0, cryst, enforce_sym, unit=std_out)
1186 0 : ABI_FREE(sorted_kg_k)
1187 :
1188 : u_mgfft = maxval(u_ngfft(1:3))
1189 0 : qpt = zero
1190 :
1191 0 : nqibz_ = 1; nqbz_ = 1; qbz_ = zero; nkbz_ = 1
1192 : ! TODO: MC technique does not seem to work as expected, even in the legacy code.
1193 : call vcgen%init(cryst, ebands%kptrlatt, nkbz_, nqibz_, nqbz_, qbz_, &
1194 0 : dtset%gw_rcut, dtset%gw_icutcoul, dtset%vcutgeo, dtset%ecuteps, comm)
1195 :
1196 : ! NB: npweps = m_npw
1197 0 : ABI_MALLOC(sqrt_vc, (m_npw))
1198 0 : my_gw_qlwl(:) = GW_Q0_DEFAULT; if (dtset%gw_nqlwl > 0) my_gw_qlwl = dtset%gw_qlwl(:,1)
1199 0 : call vcgen%get_vc_sqrt(qpt, m_npw, m_gvec, my_gw_qlwl, cryst, sqrt_vc, comm)
1200 : ! Override the G=0 component with the correct value set by gw_icutcoulomb method
1201 0 : sqrt_vc(1) = sqrt(vcgen%i_sz)
1202 : !sqrt_vc(1) = zero
1203 0 : call vcgen%free()
1204 :
1205 0 : if (my_rank == master) then
1206 0 : call wrtout(units, " Computing oscilator matrix elements for CC4S.")
1207 0 : call print_ngfft([std_out], u_ngfft, header='FFT mesh for wavefunctions')
1208 :
1209 : ! =====================
1210 : ! Write files for CC4S
1211 : ! =====================
1212 0 : if (ik_ibz == 1 .and. spin == 1) then
1213 : ! Write g-vector files, see https://manuals.cc4s.org/user-manual/objects/GridVectors.html
1214 0 : filepath = trim(dtfil%filnam_ds(4))//'_GridVectors.yaml'
1215 0 : write(ab_out, "(3a)")ch10," Writing Gridvectors metadata to file: ", trim(filepath)
1216 0 : if (open_file(filepath, msg, newunit=unt, access="stream", form="formatted", status="replace", action="write") /= 0) then
1217 0 : ABI_ERROR(msg)
1218 : end if
1219 0 : write(unt,'(a)') 'version: 100'
1220 0 : write(unt,'(a)') 'type: Tensor'
1221 0 : write(unt,'(a)') 'scalarType: Real64'
1222 0 : write(unt,'(a)') 'dimensions:'
1223 0 : write(unt,'(a)') ' - length: 3'
1224 0 : write(unt,'(a)') ' type: Vector'
1225 0 : write(unt,'(a,i0)') ' - length: ',m_npw
1226 0 : write(unt,'(a)') ' type: Momentum'
1227 0 : write(unt,'(a)') 'elements:'
1228 0 : write(unt,'(a)') ' type: TextFile'
1229 0 : write(unt,'(a)') 'unit: 1.0 # Bohr^-1'
1230 : ! The last three lines correspond to the reciprocal lattice vectors (including the factor 2pi)
1231 0 : write(unt,'(a)') 'metaData:'
1232 : write(unt,'(" Gi: [",E22.15,",",E22.15,",",E22.15,"]")') &
1233 0 : two_pi*cryst%gprimd(1,1), two_pi*cryst%gprimd(2,1), two_pi*cryst%gprimd(3,1)
1234 : write(unt,'(" Gj: [",E22.15,",",E22.15,",",E22.15,"]")') &
1235 0 : two_pi*cryst%gprimd(1,2), two_pi*cryst%gprimd(2,2), two_pi*cryst%gprimd(3,2)
1236 : write(unt,'(" Gk: [",E22.15,",",E22.15,",",E22.15,"]")') &
1237 0 : two_pi*cryst%gprimd(1,3), two_pi*cryst%gprimd(2,3), two_pi*cryst%gprimd(3,3)
1238 0 : close(unt)
1239 :
1240 : ! Write g-vectors (cart coords)
1241 0 : filepath = trim(dtfil%filnam_ds(4))//'_GridVectors.elements'
1242 0 : if (open_file(filepath, msg, newunit=unt, access="stream", form="formatted", status="replace", action="write") /= 0) then
1243 0 : ABI_ERROR(msg)
1244 : end if
1245 :
1246 : !TODO: k-points to be implemented.
1247 0 : do ig=1,m_npw
1248 0 : gcart = two_pi * matmul(cryst%gprimd, m_gvec(:,ig))
1249 0 : do ii=1,3
1250 0 : write(unt, *) gcart(ii)
1251 : end do
1252 : end do
1253 0 : close(unt)
1254 :
1255 : ! https://manuals.cc4s.org/user-manual/objects/CoulombVertex.html
1256 0 : filepath = trim(dtfil%filnam_ds(4))//'_CoulombVertex.yaml'
1257 0 : write(ab_out, "(3a)")ch10, ' Writing CoulombVertex metadata to file: ', trim(filepath)
1258 0 : if (open_file(filepath, msg, newunit=unt, access="stream", form="formatted", status="replace", action="write") /= 0) then
1259 0 : ABI_ERROR(msg)
1260 : end if
1261 0 : write(unt,'(a)') 'version: 100'
1262 0 : write(unt,'(a)') 'type: Tensor'
1263 0 : write(unt,'(a)') 'scalarType: Complex64'
1264 0 : write(unt,'(a)') 'dimensions:'
1265 0 : write(unt,'(a,i0)') '- length: ',m_npw
1266 0 : write(unt,'(a)') ' type: AuxiliaryField'
1267 0 : write(unt,'(a,i0)') '- length: ',ugb%nband_k
1268 0 : write(unt,'(a)') ' type: State'
1269 0 : write(unt,'(a,i0)') '- length: ',ugb%nband_k
1270 0 : write(unt,'(a)') ' type: State'
1271 0 : write(unt,'(a)') 'elements:'
1272 0 : write(unt,'(a)') ' type: IeeeBinaryFile'
1273 0 : write(unt,'(a)') 'unit: 1.0 # Atomic units'
1274 0 : write(unt,'(a)') 'metaData:'
1275 0 : if (m_istwfk == 2) then
1276 0 : write(unt,'(a)') ' halfGrid: 1'
1277 : else
1278 0 : write(unt,'(a)') ' halfGrid: 0'
1279 : end if
1280 0 : close(unt)
1281 :
1282 : ! https://manuals.cc4s.org/user-manual/objects/CoulombPotential.html
1283 0 : filepath = trim(dtfil%filnam_ds(4))//'_CoulombPotential.yaml'
1284 0 : write(ab_out, "(3a)")ch10, ' Writing CoulombPotential metadata to file: ', trim(filepath)
1285 0 : if (open_file(filepath, msg, newunit=unt, access="stream", form="formatted", status="replace", action="write") /= 0) then
1286 0 : ABI_ERROR(msg)
1287 : end if
1288 0 : write(unt,'(a)') 'version: 100'
1289 0 : write(unt,'(a)') 'type: Tensor'
1290 0 : write(unt,'(a)') 'scalarType: Real64'
1291 0 : write(unt,'(a)') 'dimensions:'
1292 0 : write(unt,'(a,i0)') ' - length: ',m_npw
1293 0 : write(unt,'(a)') ' type: Momentum'
1294 0 : write(unt,'(a)') 'elements:'
1295 0 : write(unt,'(a)') ' type: TextFile'
1296 0 : write(unt,'(a)') 'unit: 1.0 # Atomic units '
1297 0 : close(unt)
1298 :
1299 0 : filepath = trim(dtfil%filnam_ds(4))//'_CoulombPotential.elements'
1300 0 : write(ab_out, "(3a)")ch10, ' Writing CoulombPotential data to file: ', trim(filepath)
1301 0 : if (open_file(filepath, msg, newunit=unt, access="stream", form="formatted", status="replace", action="write") /= 0) then
1302 0 : ABI_ERROR(msg)
1303 : end if
1304 0 : do ig=1,m_npw
1305 0 : write(unt,*)real(sqrt_vc(ig) * conjg(sqrt_vc(ig)), kind=dp)
1306 : !write(unt,*)real(sqrt_vc(ig)**2)
1307 : end do
1308 0 : close(unt)
1309 : end if ! ik_ibz == 1 .and. spin == 1
1310 : end if ! my_rank == master
1311 :
1312 : ! Open binary file to store CoulombVertex
1313 0 : cvx_filepath = trim(dtfil%filnam_ds(4))//'_CoulombVertex.elements'
1314 0 : if (my_rank == master) write(ab_out, "(3a)")ch10, ' Writing CoulombVertex data to file: ', trim(cvx_filepath)
1315 :
1316 : #ifdef HAVE_MPI_IO
1317 0 : call MPI_FILE_OPEN(comm, cvx_filepath, MPI_MODE_CREATE + MPI_MODE_WRONLY, xmpio_info, fh, mpierr)
1318 0 : ABI_CHECK_MPI(mpierr, "MPI_FILE_OPEN")
1319 : #else
1320 : ABI_ERROR("CC4S interface requires MPI-IO!")
1321 : #endif
1322 :
1323 : if (trust_no_one .and. my_rank == 0) then
1324 : if (open_file("test_mg", msg, newunit=test_unt, form="formatted", status="replace", action="write") /= 0) then
1325 : ABI_ERROR(msg)
1326 : end if
1327 : end if
1328 :
1329 : ! Define batch sizes and allocate workspace arrays.
1330 : ! Increasing this value improves efficiency (less communication) at the price of more memory.
1331 0 : batch1_size = min(48, ugb%nband_k); batch2_size = min(48, ugb%nband_k) !; batch1_size = 1; batch2_size = 1
1332 0 : call wrtout(std_out, sjoin(" Using batch1_size:", itoa(batch1_size), ", batch2_size:", itoa(batch2_size)))
1333 :
1334 : mem_mb = (two * u_nfft * nspinor * batch1_size + &
1335 : two * u_nfft * nspinor * batch2_size * two + &
1336 0 : two * m_npw * nspinor * batch2_size) * dp * b2Mb
1337 0 : call wrtout(std_out, sjoin(" Memory for workspace arrays ", ftoa(mem_mb, fmt="f8.1"), "[Mb] <<< MEM"))
1338 :
1339 0 : ABI_MALLOC(ur1_batch, (u_nfft * nspinor, batch1_size))
1340 0 : ABI_MALLOC(ur2_batch, (u_nfft * nspinor, batch2_size))
1341 0 : ABI_MALLOC(ur12_batch, (u_nfft * nspinor, batch2_size))
1342 0 : ABI_MALLOC(ug12_batch, (m_npw * nspinor, batch2_size))
1343 :
1344 0 : if (psps%usepaw == 1) then
1345 : ! Evaluate oscillator matrix elements btw partial waves. Note q=Gamma
1346 0 : ABI_MALLOC(pwij, (psps%ntypat))
1347 0 : call pawpwij_init(pwij, m_npw, qpt, m_gvec, cryst%rprimd, psps, pawtab, paw_pwff)
1348 0 : do iatom=1,cryst%natom
1349 0 : nlmn_atm(iatom) = pawtab(cryst%typat(iatom))%lmn_size
1350 : end do
1351 0 : ABI_MALLOC(cprj1, (cryst%natom, nspinor*batch1_size))
1352 0 : call pawcprj_alloc(cprj1, 0, nlmn_atm)
1353 0 : dim_rtwg = nspinor
1354 0 : ABI_MALLOC(paw_rhotwg, (m_npw*dim_rtwg))
1355 : end if
1356 :
1357 : ! TODO:
1358 : ! 1) take advantage of M_{b1,b2}(g) = <b1|e^{-ig.r}|b2> => M_{b1,b2}(g) = M_{b2,b1}(-g)^*
1359 : ! once I have a better understanding of the fileformat expected by CC4S.
1360 : ! 2) Handle parallel IO if nsppol 2 (we are inside the spin loop that is already MPI distributed!)
1361 : ! 3) Clarify ordering of CoulombVertex (b1,b2 vs b2,b1) and eigenvalues (spin?)
1362 : ! 4) Treatment of q--> 0 in vc_coul
1363 : ! 4) See other TODOs below.
1364 :
1365 0 : call uplan_1%init(npw_k, nspinor, batch1_size, u_ngfft, ugb%istwf_k, ugb%kg_k, dp, dtset%gpu_option)
1366 0 : call uplan_2%init(npw_k, nspinor, batch2_size, u_ngfft, ugb%istwf_k, ugb%kg_k, dp, dtset%gpu_option)
1367 0 : call uplan_m%init(m_npw, nspinor, batch2_size, u_ngfft, m_istwfk, m_gvec, dp, dtset%gpu_option)
1368 :
1369 0 : M_ = m_npw
1370 :
1371 : ! Blocked loop over group of b1 indices. NB: Assuming bands distributed in contiguous blocks.
1372 0 : do band1_start=1, ugb%nband_k, batch1_size
1373 : ! Collect n1dat bands starting from band1_start on each proc.
1374 0 : n1dat = blocked_loop(band1_start, ugb%nband_k, batch1_size)
1375 :
1376 0 : call ugb%mat%collect_cplx(ugb%npwsp, n1dat, [1, band1_start], ug1_batch)
1377 0 : if (psps%usepaw == 1) call ugb%collect_cprj(nspinor, n1dat, band1_start, cprj1)
1378 :
1379 : ! FFT: ug1_batch --> ur1_batch
1380 0 : call uplan_1%execute_gr(n1dat, ug1_batch(:,1), ur1_batch(:,1))
1381 0 : if (ugb%istwf_k /= 2) ur1_batch = conjg(ur1_batch) ! Not needed if k == Gamma as ur1 is real.
1382 0 : ABI_FREE(ug1_batch)
1383 :
1384 : ! Blocked loop over MY group of b2 indices (contiguous blocks)
1385 0 : do band2_start=ugb%my_bstart, ugb%my_bstop, batch2_size
1386 0 : n2dat = blocked_loop(band2_start, ugb%my_bstop, batch2_size)
1387 0 : my_ib2st = band2_start - ugb%my_bstart + 1
1388 :
1389 : ! FFT: ugb%mat --> ur2_batch for n2dat states.
1390 0 : call uplan_2%execute_gr(n2dat, ugb%mat%buffer_cplx(:,my_ib2st), ur2_batch(:,1))
1391 :
1392 : ! For each row of the submatrix, build n2dat products (band1, idat2) in r-space, then r --> g.
1393 0 : do idat1=1,n1dat
1394 0 : band1 = band1_start + idat1 - 1
1395 :
1396 0 : do idat2=1,n2dat
1397 0 : ur12_batch(:,idat2) = ur1_batch(:,idat1) * ur2_batch(:,idat2)
1398 : end do
1399 0 : call uplan_m%execute_rg(n2dat, ur12_batch(:,1), ug12_batch(:,1))
1400 :
1401 0 : if (psps%usepaw == 1) then
1402 : ! Add PAW on-site contributions
1403 0 : do idat2=1,n2dat
1404 : associate (cprj1_kmq => cprj1(:, 1 + (idat1-1)*nspinor), &
1405 0 : cprj2_k => ugb%cprj_k(:, 1 + (my_ib2st+idat2-2)*nspinor)) ! NB: ugb%cprj_k(2, nspinor*my_nband)
1406 0 : paw_rhotwg = zero
1407 0 : call paw_rho_tw_g(cryst, pwij, m_npw, dim_rtwg, nspinor, m_gvec, cprj1_kmq, cprj2_k, paw_rhotwg)
1408 0 : ug12_batch(:,idat2) = ug12_batch(:,idat2) + paw_rhotwg
1409 : end associate
1410 : end do
1411 : end if
1412 :
1413 : ! This to zero the matrix elements for certain (band1, band2) entries.
1414 : ! in this case: band1 = 1, band2 = 1.
1415 : !if (band1 == 1) then
1416 : ! do idat2=1,n2dat
1417 : ! ! This is the condition for band2
1418 : ! if (band2_start + idat2 - 1 == 1) ug12_batch(:,idat2) = zero
1419 : ! end do
1420 : !end if
1421 :
1422 0 : if (nspinor == 2) then
1423 : ! Sum over spinors and repack data in the first n2dat positions to prepare IO operation.
1424 0 : do idat2=1,n2dat
1425 0 : ug12_batch(1:m_npw,idat2) = ug12_batch(1:m_npw,idat2) + ug12_batch(m_npw+1:,idat2)
1426 : end do
1427 0 : do idat2=2,n2dat,2
1428 0 : ug12_batch(m_npw+1:,idat2-1) = ug12_batch(1:m_npw,idat2)
1429 : end do
1430 : end if
1431 :
1432 0 : do idat2=1,n2dat
1433 : !if (band1 == band2_start + idat2 -1) then
1434 : ! write(std_out,*) " ug12_batch(g=0,band1,band2), band", ug12_batch(1,idat2), band1, band2_start + idat2 -1
1435 : !end if
1436 : ! Multiply by sqrt(vc(g))
1437 : ! This renormalization is needed to make CC4S converge (work done by AlejandroG, FabienB, MatteoG)
1438 0 : ug12_batch(:,idat2) = ug12_batch(:,idat2) * sqrt_vc(:) / sqrt(cryst%ucvol)
1439 : !FBru:
1440 : ! CC4S convention:
1441 : ! When using time-reversal, include a factor sqrt(2) for the non-zero G's
1442 : ! Then CC4S will perform the sum only over half the grid so all the G/=0 component should be counted twice.
1443 0 : if (m_istwfk == 2) ug12_batch(2:,idat2) = ug12_batch(2:,idat2) * sqrt(two)
1444 : end do
1445 : !write(std_out,*)" max(abs(ug12_batch)):", maxval(abs(ug12_batch(:,1:n2dat)))
1446 :
1447 : #ifdef HAVE_MPI_IO
1448 : ! Write ug12_batch using Stream-IO
1449 0 : buf_size = m_npw * n2dat
1450 :
1451 : ! FBru convert to longer integer before the multiplication
1452 : offset = INT(m_npw, KIND=XMPI_OFFSET_KIND) &
1453 : * INT( (band2_start-1) + (band1-1) * ugb%nband_k, KIND=XMPI_OFFSET_KIND) &
1454 0 : * INT(xmpi_bsize_dpc, KIND=XMPI_OFFSET_KIND)
1455 :
1456 0 : call MPI_FILE_WRITE_AT(fh, offset, ug12_batch, buf_size, MPI_DOUBLE_COMPLEX, MPI_STATUS_IGNORE, mpierr)
1457 0 : ABI_HANDLE_MPIERR(mpierr)
1458 :
1459 0 : if (my_rank == 0 .and. trust_no_one) then
1460 : do idat2=1,n2dat
1461 : band2 = band2_start + idat2 - 1
1462 : write(test_unt,*)band1, band2, ug12_batch(1:M_,idat2)
1463 : end do
1464 : end if
1465 : #endif
1466 : end do ! idat1
1467 :
1468 : end do ! band2_start
1469 : end do ! band1_start
1470 :
1471 : #ifdef HAVE_MPI_IO
1472 0 : call MPI_FILE_CLOSE(fh, mpierr)
1473 0 : ABI_CHECK_MPI(mpierr, "FILE_CLOSE!")
1474 0 : call xmpi_barrier(comm)
1475 :
1476 0 : if (my_rank == 0) then
1477 0 : buf_size = 4
1478 0 : call wrtout(units, sjoin(" Reading norm of Coulomb vertex for testing purposes with ng:", itoa(buf_size)), newlines=1, pre_newlines=1)
1479 0 : call wrtout(units, " List of g-vectors in the Coulomb vertex:")
1480 0 : do ig=1,buf_size
1481 0 : call wrtout(units, sjoin(itoa(ig), ":", ltoa(m_gvec(:,ig))))
1482 : end do
1483 0 : ABI_MALLOC(cwork, (buf_size))
1484 0 : call MPI_FILE_OPEN(xmpi_comm_self, cvx_filepath, MPI_MODE_RDONLY, xmpio_info, fh, mpierr)
1485 0 : ABI_CHECK_MPI(mpierr, "MPI_FILE_OPEN")
1486 0 : ierr = 0
1487 0 : band1_loop: do band1=1, ugb%nband_k
1488 0 : do band2=1, ugb%nband_k
1489 0 : ierr = ierr + 1; if (ierr == 6) exit band1_loop
1490 0 : offset = ((band2-1) * m_npw + (band1-1) * m_npw * ugb%nband_k) * xmpi_bsize_dpc
1491 0 : call MPI_FILE_READ_AT(fh, offset, cwork, buf_size, MPI_DOUBLE_COMPLEX, MPI_STATUS_IGNORE, mpierr)
1492 0 : ABI_HANDLE_MPIERR(mpierr)
1493 0 : call wrtout(units, sjoin(" For band1:", itoa(band1), ", band2:", itoa(band2)))
1494 0 : where (abs(cwork) < tol8)
1495 : cwork = zero
1496 : end where
1497 0 : write(msg, "(*(1x, es12.5))")abs(cwork(1:buf_size))
1498 0 : call wrtout(units, msg)
1499 : end do
1500 : end do band1_loop
1501 0 : ABI_FREE(cwork)
1502 : end if
1503 : #endif
1504 :
1505 : ! =============
1506 : ! DEBUG SECTION
1507 : ! =============
1508 : #ifdef HAVE_MPI_IO
1509 : if (my_rank == 0 .and. trust_no_one) then
1510 : close(test_unt)
1511 : if (open_file("test_mg", msg, newunit=test_unt, form="formatted", status="old", action="read") /= 0) then
1512 : ABI_ERROR(msg)
1513 : end if
1514 :
1515 : call MPI_FILE_OPEN(xmpi_comm_self, cvx_filepath, MPI_MODE_RDONLY, xmpio_info, fh, mpierr)
1516 : ABI_CHECK_MPI(mpierr, "MPI_FILE_OPEN")
1517 :
1518 : ABI_MALLOC(cwork, (m_npw))
1519 : max_abs_err = zero
1520 : do ig=1, ugb%nband_k**2
1521 : read(test_unt,*) band1, band2, ug12_batch(1:M_,1)
1522 : offset = ((band2-1) * m_npw + (band1-1) * m_npw * ugb%nband_k) * xmpi_bsize_dpc
1523 : call MPI_FILE_READ_AT(fh, offset, cwork, m_npw, MPI_DOUBLE_COMPLEX, MPI_STATUS_IGNORE, mpierr)
1524 : ABI_HANDLE_MPIERR(mpierr)
1525 :
1526 : abs_err = maxval(abs(ug12_batch(1:M_,1) - cwork(1:M_)))
1527 : max_abs_err = max(max_abs_err, abs_err)
1528 : if (abs_err > zero) write(std_out, *)" For ig:", ig, "/", ugb%nband_k**2, "abs_err", abs_err
1529 : !write(std_out, *)"1:", ug12_batch(1:M_,1); write(std_out, *)"2:", cwork(1:M_)
1530 : end do
1531 :
1532 : close(test_unt)
1533 : call MPI_FILE_CLOSE(fh, mpierr)
1534 : ABI_CHECK_MPI(mpierr, "FILE_CLOSE!")
1535 : ABI_FREE(cwork)
1536 :
1537 : write(std_out,*)" max_abs_err:", max_abs_err
1538 : ABI_CHECK(max_abs_err < tol16, sjoin("max_abs_err:", ftoa(max_abs_err)))
1539 : call wrtout(std_out, " Debugging section OK!!!")
1540 : end if
1541 : #endif
1542 :
1543 : ! Free memory
1544 0 : call uplan_1%free(); call uplan_2%free(); call uplan_m%free()
1545 :
1546 0 : ABI_FREE(m_gvec)
1547 0 : ABI_FREE(ur1_batch)
1548 0 : ABI_FREE(ur2_batch)
1549 0 : ABI_FREE(ur12_batch)
1550 0 : ABI_FREE(ug12_batch)
1551 0 : ABI_FREE(sqrt_vc)
1552 :
1553 0 : if (psps%usepaw == 1) then
1554 0 : call pawpwij_free(pwij)
1555 0 : ABI_FREE(pwij)
1556 0 : call pawcprj_free(cprj1)
1557 0 : ABI_FREE(cprj1)
1558 0 : ABI_FREE(paw_rhotwg)
1559 : end if
1560 :
1561 0 : call cwtime_report(" cc4s_gamma", cpu, wall, gflops)
1562 :
1563 0 : end subroutine cc4s_gamma
1564 : !!***
1565 :
1566 0 : end module m_gwr_driver
1567 : !!***
|
