Line data Source code
1 : !!****m* ABINIT/m_outscfcv
2 : !! NAME
3 : !! m_outscfcv
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 2005-2026 ABINIT group (XG)
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! SOURCE
14 :
15 : #if defined HAVE_CONFIG_H
16 : #include "config.h"
17 : #endif
18 :
19 : #include "abi_common.h"
20 :
21 : module m_outscfcv
22 :
23 : use defs_basis
24 : use defs_wvltypes
25 : use m_abicore
26 : use m_sort
27 : use m_efield
28 : use m_errors
29 : use m_xmpi
30 : use m_mpinfo
31 : use netcdf
32 : use m_nctk
33 : use m_hdr
34 : use m_plowannier
35 : use m_splines
36 : use m_ebands
37 : use m_dtset
38 : use m_dtfil
39 : use m_rcpaw
40 :
41 : use defs_datatypes, only : pseudopotential_type
42 : use defs_abitypes, only : MPI_type
43 : use m_time, only : timab
44 : use m_io_tools, only : open_file
45 : use m_fstrings, only : strcat, endswith
46 : use m_geometry, only : bonds_lgth_angles
47 : use m_electronpositron, only : electronpositron_type,electronpositron_calctype
48 : use m_oper, only : oper_type,init_oper,destroy_oper
49 : use m_crystal, only : crystal_t, prt_cif
50 : use m_results_gs, only : results_gs_type
51 : use m_ioarr, only : ioarr, fftdatar_write
52 : use m_matlu, only : copy_matlu,destroy_matlu,init_matlu,matlu_type
53 : use m_nucprop, only : calc_efg,calc_fc
54 : use m_outwant, only : outwant
55 : use m_rcpaw, only : rcpaw_type
56 : use m_pawang, only : pawang_type
57 : use m_pawrad, only : pawrad_type, simp_gen, bound_deriv
58 : use m_pawtab, only : pawtab_type
59 : use m_paw_an, only : paw_an_type
60 : use m_paw_ij, only : paw_ij_type
61 : use m_paw_mkrho, only : denfgr
62 : use m_pawfgrtab, only : pawfgrtab_type
63 : use m_pawrhoij, only : pawrhoij_type, pawrhoij_nullify, pawrhoij_copy, pawrhoij_free
64 : use m_pawcprj, only : pawcprj_type
65 : use m_pawfgr, only : pawfgr_type
66 : use m_paw_dmft, only : paw_dmft_type,init_dmft,destroy_dmft,print_dmft
67 : use m_paw_optics, only : optics_paw,optics_paw_core
68 : use m_paw_tools, only : pawprt
69 : use m_numeric_tools, only : simpson_int
70 : use m_epjdos, only : epjdos_t
71 : use m_paral_atom, only : get_my_atmtab, free_my_atmtab
72 : use m_io_kss, only : outkss
73 : use m_multipoles, only : multipoles_out, out1dm
74 : use m_mlwfovlp_qp, only : mlwfovlp_qp
75 : use m_paw_mkaewf, only : pawmkaewf
76 : use m_dens, only : mag_penalty_e, calcdenmagsph, prtdenmagsph
77 : !use m_mlwfovlp, only : mlwfovlp
78 : use m_wfd_wannier, only : wfd_run_wannier
79 : use m_datafordmft, only : datafordmft
80 : use m_mkrho, only : read_atomden
81 : use m_positron, only : poslifetime, posdoppler
82 : use m_optics_vloc, only : optics_vloc
83 : use m_green, only : green_type,compute_green,&
84 : fourier_green,print_green,init_green,destroy_green,init_green_tau
85 : use m_self, only : self_type,initialize_self,rw_self,destroy_self,destroy_self,selfreal2imag_self
86 : use m_paw_correlations, only : loc_orbmom_cal
87 : use m_paw_denpot, only : paw_relax_core
88 :
89 : implicit none
90 :
91 : private
92 : !!***
93 :
94 : public :: outscfcv
95 : !!***
96 :
97 : contains
98 : !!***
99 :
100 : !!****f* ABINIT/outscfcv
101 : !! NAME
102 : !! outscfcv
103 : !!
104 : !! FUNCTION
105 : !! Output routine for the scfcv.F90 routine
106 : !!
107 : !! INPUTS
108 : !! atindx(natom)=index table for atoms (see gstate.f)
109 : !! atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
110 : !! cg(2,mcg)=planewave coefficients of wavefunctions (see also side effects)
111 : !! compch_fft=compensation charge, from FFT grid
112 : !! compch_sph=compensation charge, from sphere
113 : !! cprj(natom,mcprj*usecprj)=<p_lmn|Cnk> coefficients for each WF |Cnk>
114 : !! and each |p_lmn> non-local projector. See also side effects
115 : !! dimcprj(natom*usecprj)=array of dimensions of array cprj (not ordered)
116 : !! dmatpawu= fixed occupation matrix of correlated orbitals (DFT+U or DMFT only)
117 : !! dtfil <type(datafiles_type)>=variables related to files
118 : !! dtset <type(dataset_type)>=all input variables for this dataset
119 : !! ecut=cut-off energy for plane wave basis sphere (Ha)
120 : !! eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
121 : !! elfr(nfft,nspden(+1))=electron localization function, real space.
122 : !! (+1) if spin-polarized in order to get total, spin up and spin down elf
123 : !! etotal=total energy
124 : !! gmet(3,3)=reciprocal space metric tensor in bohr**-2.
125 : !! gprimd(3,3)=dimensional reciprocal space primitive translations
126 : !! grhor(nfft,nspden,3)= gradient of electron density in electrons/bohr**4, real space
127 : !! hdr <type(hdr_type)>=the header of wf, den and pot files
128 : !! intgres(nspden,natom)=integrated residuals from constrained DFT. They are also Lagrange parameters, or gradients with respect to constraints.
129 : !! kg(3,mpw*mkmem)=reduced planewave coordinates.
130 : !! lrhor(nfft,nspden)= Laplacian of electron density in electrons/bohr**5, real space
131 : !! mband=maximum number of bands
132 : !! mcg=size of wave-functions array (cg) =mpw*my_nspinor*mband*mkmem*nsppol
133 : !! mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
134 : !! mgfftc=maximum size of 1D FFTs for the PAW coarse grid
135 : !! mkmem =number of k points treated by this node.
136 : !! mpi_enreg=information about MPI parallelization
137 : !! mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
138 : !! mpw=maximum dimensioned size of npw.
139 : !! my_natom=number of atoms treated by current processor
140 : !! natom=number of atoms in cell.
141 : !! nattyp(ntypat)= # atoms of each type.
142 : !! nfft=(effective) number of FFT grid points (for this processor) (see NOTES at beginning of scfcv)
143 : !! ngfft(18)=contain all needed information about 3D FFT (see NOTES at beginning of scfcv)
144 : !! nhat(nfft,nspden*usepaw)= compensation charge density (PAW)
145 : !! nkpt=number of k points.
146 : !! npwarr(nkpt)=number of planewaves in basis at this k point
147 : !! nspden=number of spin-density components
148 : !! nsppol=1 for unpolarized, 2 for spin-polarized
149 : !! nsym=number of symmetries in space group
150 : !! ntypat=number of types of atoms in unit cell.
151 : !! n3xccc=dimension of the xccc3d array (0 or nfft).
152 : !! occ(mband*nkpt*nsppol)=occupation number for each band (usually 2) for each k.
153 : !! paw_an(my_natom) <type(paw_an_type)>=paw arrays given on angular mesh
154 : !! pawang <type(pawang_type)>=paw angular mesh and related data
155 : !! pawfgr(natom) <type(pawfgr_type)>=fine grid parameters and related data
156 : !! pawfgrtab(my_natom) <type(pawfgrtab_type)> tables on PAW fine grid
157 : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
158 : !! pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
159 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
160 : !! paw_ij(my_natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels
161 : !! note:structure factors are given on the coarse grid for PAW
162 : !! prtvol=control print volume and debugging output
163 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
164 : !! results_gs <type(results_gs_type)>=results (energy and its components,
165 : !! forces and its components, the stress tensor) of a ground-state computation
166 : !! rhor(nfft,nspden)=total electron density in electrons/bohr**3, real space.
167 : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
168 : !! taur(nfft,nspden)=total kinetic energy density in bohr**(-5), real space.
169 : !! ucvol=unit cell volume (bohr**3)
170 : !! usecprj=1 if cprj datastructure has been allocated
171 : !! vhartr(nfft)=Hartree potential
172 : !! vxc(nfft,nspden)=xc potential
173 : !! vtrial(nfft,nspden)=the trial potential = vxc + vpsp + vhartr, roughly speaking
174 : !! xccc3d(n3xccc)=3D core electron density for XC core correction (bohr^-3)
175 : !! xred(3,natom)=reduced dimensionless atomic coordinates
176 : !!
177 : !! OUTPUT
178 : !! (only writing, printing)
179 : !!
180 : !! SIDE EFFECTS
181 : !! electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation
182 : !! If prtwant==3 the following quantitities are updated using the unitary transformation
183 : !! defining the QP amplitudes in terms of the KS basis set:
184 : !! cg(2,mcg)=planewave coefficients of wavefunctions.
185 : !! cprj(natom,mcprj*usecpyj)=<p_lmn|Cnk> coefficients for each WF |Cnk> and each |p_lmn> non-local projector
186 : !!
187 : !! NOTES
188 : !! The string passed to fftdatar_write (first argument) gives the name used to store the data in the netcdf file
189 : !! The function varname_from_fname defined in the module m_hdr.F90 gives the mapping between the Abinit
190 : !! file extension and the netcdf name e.g. foo_VHXC.nc --> vxc
191 : !! This function is used in cut3d so that we can immediately select the data to analyze without having
192 : !! to prompt the user. Remember to update varname_from_fname if you add a new file or if you change the
193 : !! name of the variable.
194 : !!
195 : !! SOURCE
196 :
197 20754 : subroutine outscfcv(atindx1,cg,compch_fft,compch_sph,cprj,dimcprj,dmatpawu,dtfil,dtset,&
198 13836 : & ecut,eigen,electronpositron,elfr,etotal,gmet,gprimd,grhor,hdr,intgres,kg,&
199 : & lrhor,mband,mcg,mcprj,mgfftc,mkmem,mpi_enreg,mpsang,mpw,my_natom,natom,&
200 6918 : & nattyp,nfft,ngfft,nhat,nkpt,npwarr,nspden,nsppol,nsym,ntypat,n3xccc,occ,&
201 6918 : & paw_dmft,pawang,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,paw_an,paw_ij,&
202 6918 : & prtvol,psps,results_gs,rhor,rprimd,&
203 13836 : & taur,ucvol,usecprj,vhartr,vpsp,vtrial,vxc,wvl_den,xccc3d,xred,rcpaw)
204 :
205 : !Arguments ------------------------------------
206 : !scalars
207 : integer,intent(in) :: mband,mcg,mcprj,mgfftc,mkmem,mpsang,mpw,n3xccc,my_natom,natom,nfft
208 : integer,intent(in) :: nkpt,nspden,nsppol,nsym,ntypat,prtvol,usecprj
209 : real(dp),intent(in) :: compch_fft,compch_sph,ecut,ucvol
210 : real(dp),intent(inout) :: etotal
211 : type(electronpositron_type),pointer :: electronpositron
212 : type(MPI_type),intent(inout) :: mpi_enreg
213 : type(datafiles_type),intent(in) :: dtfil
214 : type(dataset_type),intent(in) :: dtset
215 : type(hdr_type),intent(inout) :: hdr
216 : type(paw_dmft_type), intent(inout) :: paw_dmft
217 : type(pawang_type),intent(in) :: pawang
218 : type(pawfgr_type),intent(in) :: pawfgr
219 : type(pseudopotential_type),intent(inout) :: psps
220 : type(results_gs_type),intent(in) :: results_gs
221 : type(wvl_denspot_type), intent(in) :: wvl_den
222 : type(rcpaw_type),intent(in),pointer :: rcpaw
223 : !arrays
224 : integer,intent(in) :: atindx1(natom),dimcprj(natom*usecprj)
225 : integer,intent(in) :: kg(3,mpw*mkmem),nattyp(ntypat),ngfft(18),npwarr(nkpt)
226 : real(dp),intent(in) :: dmatpawu(:,:,:,:),eigen(mband*nkpt*nsppol)
227 : real(dp),intent(in) :: gmet(3,3),gprimd(3,3)
228 : real(dp),intent(in) :: intgres(:,:) ! (nspden,natom) if constrainedDFT otherwise (nspden,0)
229 : real(dp),intent(in) :: occ(mband*nkpt*nsppol)
230 : real(dp),intent(in) :: rprimd(3,3),vhartr(nfft),xccc3d(n3xccc)
231 : real(dp),intent(in) :: vpsp(nfft)
232 : real(dp),intent(inout) :: cg(2,mcg)
233 : real(dp),intent(inout) :: nhat(nfft,nspden*psps%usepaw)
234 : real(dp),intent(inout),target :: rhor(nfft,nspden),vtrial(nfft,nspden)
235 : real(dp),intent(inout) :: vxc(nfft,nspden),xred(3,natom)
236 : real(dp),pointer :: elfr(:,:),grhor(:,:,:),lrhor(:,:),taur(:,:)
237 : type(pawcprj_type),intent(inout) :: cprj(natom,mcprj*usecprj)
238 : type(paw_an_type),intent(inout) :: paw_an(my_natom*psps%usepaw)
239 : type(pawfgrtab_type),intent(in) :: pawfgrtab(my_natom*psps%usepaw)
240 : type(paw_ij_type),intent(inout) :: paw_ij(my_natom*psps%usepaw)
241 : type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw)
242 : type(pawrhoij_type),target,intent(inout) :: pawrhoij(my_natom*psps%usepaw)
243 : type(pawtab_type),intent(inout) :: pawtab(ntypat*psps%usepaw)
244 :
245 : !Local variables-------------------------------
246 : !scalars
247 : integer,parameter :: master=0,cplex1=1,fform_den=52,rdwr2=2,rdwrpaw0=0
248 : integer :: bantot,fform,collect,timrev, accessfil,coordn,option
249 : integer :: ii,ierr,ifft,ikpt,ispden,isppol,itypat, me_fft,n1,n2,n3
250 : integer :: ifgd, iatom, iatom_tot,nradint, me,my_natom_tmp
251 : integer :: n,nn,iln,iln2,l,ll
252 : integer :: occopt, opt_moments, prtnabla, pawprtden, ncid, ncerr,nphicor
253 : integer :: iband,nocc,comm,comm_fft,tmp_unt,nfft_tot, my_comm_atom, opt_imagonly
254 13836 : integer :: indsym(4,dtset%nsym,dtset%natom)
255 : real(dp) :: norm,occ_norm,unocc_norm, rate_dum,rate_dum2, yp1, ypn, dr
256 : character(len=500) :: msg
257 : character(len=fnlen) :: fname
258 : type(rcpaw_type),pointer :: rcpaw_dirac=>null()
259 : !arrays
260 : integer :: units(2)
261 6918 : integer, allocatable :: isort(:)
262 6918 : integer, pointer :: my_atmtab(:)
263 6918 : real(dp) :: tsec(2),nt_ntone_norm(nspden),rhomag(2,nspden)
264 6918 : real(dp),allocatable :: efg(:,:,:),eigen2(:)
265 6918 : real(dp),allocatable :: elfr_down(:,:),elfr_up(:,:),intgden(:,:)
266 6918 : real(dp),allocatable :: rhor_paw(:,:),rhor_paw_core(:,:),rhor_paw_val(:,:),vpaw(:,:),vwork(:,:)
267 6918 : real(dp),allocatable :: rhor_n_one(:,:),rhor_nt_one(:,:),ps_norms(:,:,:)
268 6918 : real(dp), allocatable :: doccde(:)
269 6918 : real(dp), allocatable :: vh1spl(:)
270 6918 : real(dp), allocatable :: vh1_interp(:)
271 6918 : real(dp), allocatable :: vh1_integ(:)
272 6918 : real(dp), allocatable :: vh1_corrector(:)
273 6918 : real(dp), allocatable :: radii(:)
274 6918 : real(dp), ABI_CONTIGUOUS pointer :: rho_ptr(:,:)
275 69180 : type(pawrhoij_type) :: pawrhoij_dum(1)
276 : !type(pawrhoij_type) :: pawrhoij_dum(0)
277 6918 : type(pawrhoij_type),pointer :: pawrhoij_all(:)
278 : logical :: remove_inv
279 : logical :: paral_atom, paral_fft, my_atmtab_allocated
280 : real(dp) :: e_hspinfield
281 : real(dp) :: dmatdum(0,0,0,0)
282 : real(dp) :: e_fermie, e_fermih
283 6918 : type(oper_type) :: dft_occup
284 359736 : type(crystal_t) :: crystal
285 6918 : type(ebands_t) :: ebands
286 6918 : type(epjdos_t) :: dos
287 6918 : type(plowannier_type) :: wan
288 6918 : type(self_type) :: selfr
289 6918 : type(self_type), target :: self
290 6918 : type(green_type) :: greenr
291 6918 : type(matlu_type), allocatable :: opt_selflimit(:)
292 6918 : integer,allocatable :: nphicor_arr(:),lcor(:,:)
293 6918 : real(dp),allocatable :: energy_cor(:,:),occ_cor(:,:)
294 :
295 : ! *************************************************************************
296 :
297 : DBG_ENTER("COLL")
298 :
299 20754 : units = [std_out, ab_out]
300 :
301 6918 : call timab(1150,1,tsec) ! outscfcv
302 6918 : call timab(1151,1,tsec) ! outscfcv(preparation)
303 :
304 6918 : if ((usecprj==0.or.mcprj==0).and.psps%usepaw==1.and. &
305 : (dtset%prtwant==2.or.dtset%prtwant==3.or.dtset%prtnabla>0.or.dtset%prtdos==3 &
306 : .or.dtset%kssform==3.or.dtset%pawfatbnd>0.or.dtset%pawprtwf>0)) then
307 : write (msg,'(5a)')&
308 0 : & 'cprj datastructure must be allocated',ch10,&
309 0 : & 'with options prtwant=2,3, prtnabla>0, prtdos>3, kssform==3, pawfatbnd>0, pawprtwf>0',ch10,&
310 0 : & 'Action: change pawusecp input keyword.'
311 0 : ABI_ERROR(msg)
312 : end if
313 :
314 : ! Parameters for MPI-FFT
315 6918 : n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3); nfft_tot = product(ngfft(1:3))
316 6918 : comm_fft = mpi_enreg%comm_fft
317 6918 : me_fft = xmpi_comm_rank(comm_fft)
318 6918 : paral_fft = (mpi_enreg%paral_kgb==1)
319 :
320 6918 : comm = mpi_enreg%comm_cell
321 6918 : me = xmpi_comm_rank(comm)
322 :
323 6918 : paral_atom=(my_natom/=natom)
324 6918 : my_comm_atom = mpi_enreg%comm_atom
325 6918 : nullify(my_atmtab)
326 6918 : if (paral_atom) then
327 558 : call get_my_atmtab(mpi_enreg%comm_atom, my_atmtab, my_atmtab_allocated, paral_atom,natom,my_natom_ref=my_natom)
328 : else
329 19080 : ABI_MALLOC(my_atmtab, (natom))
330 43992 : my_atmtab = (/ (iatom, iatom=1, natom) /)
331 6360 : my_atmtab_allocated = .true.
332 : end if
333 :
334 : ! Initialize two objects to facilitate the propagation of info.
335 : ! These objects should used more frequently, actually they should
336 : ! become basic objects used in abinit.
337 :
338 : ! Crystalline structure.
339 6918 : remove_inv=.false.
340 : ! What about Wannier90 DMFT?
341 6918 : if (dtset%nspden==4 .and. dtset%usedmft==1) remove_inv=.true. ! MG: why this?
342 :
343 6918 : timrev = 2; if (any(dtset%kptopt == [3, 4])) timrev= 1
344 : call crystal%init(dtset%amu_orig(:,1),dtset%spgroup,natom,dtset%npsp,ntypat, &
345 : dtset%nsym,rprimd,dtset%typat,xred,dtset%ziontypat,dtset%znucl,timrev,&
346 : dtset%nspden==2.and.dtset%nsppol==1,remove_inv,hdr%title,&
347 13806 : dtset%symrel,dtset%tnons,dtset%symafm)
348 :
349 : ! Electron band energies.
350 6918 : bantot= dtset%mband*dtset%nkpt*dtset%nsppol
351 511304 : ABI_CALLOC(doccde, (bantot))
352 : call ebands%init(bantot,dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%ivalence,&
353 : doccde,eigen,hdr%istwfk,hdr%kptns,hdr%nband,&
354 : hdr%nkpt,hdr%npwarr,hdr%nsppol,hdr%nspinor,hdr%tphysel,hdr%tsmear,hdr%occopt,hdr%occ,hdr%wtk,&
355 : hdr%cellcharge, hdr%kptopt, hdr%kptrlatt_orig, hdr%nshiftk_orig, hdr%shiftk_orig, &
356 6918 : hdr%kptrlatt, hdr%nshiftk, hdr%shiftk)
357 :
358 6918 : ABI_FREE(doccde)
359 :
360 6918 : ebands%fermie = results_gs%energies%e_fermie
361 6918 : e_fermie = results_gs%energies%e_fermie
362 6918 : ebands%fermih = results_gs%energies%e_fermih
363 6918 : e_fermih = results_gs%energies%e_fermih
364 6918 : ebands%entropy = results_gs%energies%entropy_ks
365 :
366 : ! YAML output
367 6918 : if (me == master) then
368 : call results_gs%yaml_write(ab_out, cryst=crystal, info="Summary of ground state results",&
369 5444 : occopt=dtset%occopt, with_conv=(dtset%nstep > 0))
370 : end if
371 :
372 6918 : call timab(1151,2,tsec)
373 :
374 : !wannier interface
375 6918 : call timab(1152,1,tsec)
376 :
377 6918 : if (dtset%prtwant==2 .and. ( dtset%usedmft /= 10 ) ) then
378 : call wfd_run_wannier(cryst=crystal, ebands=ebands, hdr=hdr, mpi_enreg=mpi_enreg, &
379 : & ngfftc=ngfft, ngfftf=ngfft, dtset=dtset, dtfil=dtfil, &
380 : & pawang=pawang, pawrad=pawrad, pawtab=pawtab, psps=psps, &
381 13 : & kg=kg, cg=cg, cprj=cprj)
382 : ! else
383 : !
384 : ! call mlwfovlp(crystal, ebands, hdr, atindx1,cg,cprj,dtset,dtfil,eigen,gprimd,kg,&
385 : !& mband,mcg,mcprj,mgfftc,mkmem,mpi_enreg,mpw,natom,&
386 : !& nattyp,nfft,ngfft,nkpt,npwarr,nsppol,ntypat,occ,&
387 : !& pawang,pawrad,pawtab,prtvol,psps,rprimd,ucvol,xred)
388 :
389 6905 : else if (dtset%prtwant==3) then
390 :
391 : ! Convert cg and eigen to GW quasiparticle wave functions and eigenvalues in mlwfovlp_qp
392 3 : ABI_MALLOC(eigen2,(mband*nkpt*nsppol))
393 642 : eigen2=eigen
394 :
395 : call mlwfovlp_qp(cg,cprj,dtset,dtfil,eigen2,mband,mcg,mcprj,mkmem,mpw,natom,&
396 1 : & nkpt,npwarr,nspden,nsppol,ntypat,Hdr,pawtab,rprimd,MPI_enreg)
397 :
398 : ! Call Wannier90
399 : ! call mlwfovlp(crystal, ebands, hdr, atindx1,cg,cprj,dtset,dtfil,eigen2,gprimd,kg,&
400 : !& mband,mcg,mcprj,mgfftc,mkmem,mpi_enreg,mpw,natom,&
401 : !& nattyp,nfft,ngfft,nkpt,npwarr,nsppol,ntypat,occ,&
402 : !& pawang,pawrad,pawtab,prtvol,psps,rprimd,ucvol,xred)
403 :
404 : call wfd_run_wannier(cryst=crystal, ebands=ebands, hdr=hdr, mpi_enreg=mpi_enreg, &
405 : & ngfftc=ngfft, ngfftf=ngfft, dtset=dtset, dtfil=dtfil, &
406 : & pawang=pawang, pawrad=pawrad, pawtab=pawtab, psps=psps, &
407 1 : & kg=kg, cg=cg, cprj=cprj)
408 :
409 :
410 : ! this is the old implementation, risky due to unpredictable size effects
411 : ! now eigen is not overwritten, one should use other ways to print the GW corrections
412 : ! eigen=eigen2
413 1 : ABI_FREE(eigen2)
414 : end if !prtwant
415 :
416 6918 : call timab(1152,2,tsec)
417 6918 : call timab(1153,1,tsec)
418 :
419 6918 : occopt=dtset%occopt
420 :
421 6918 : prtnabla=dtset%prtnabla
422 6918 : pawprtden=dtset%prtden-1
423 :
424 6918 : comm=mpi_enreg%comm_cell; me=xmpi_comm_rank(comm)
425 6918 : comm_fft=mpi_enreg%comm_fft
426 6918 : paral_atom=(my_natom/=natom)
427 :
428 : !Warnings :
429 : !- core charge is excluded from the charge density;
430 : !- the potential is the INPUT vtrial.
431 :
432 6918 : if (iwrite_fftdatar(mpi_enreg) .and. dtset%usewvl==0) then
433 :
434 : ! output the density.
435 5562 : if (dtset%prtden/=0) then
436 3574 : if (dtset%positron/=1) rho_ptr => rhor
437 3574 : if (dtset%positron==1) rho_ptr => electronpositron%rhor_ep
438 : call fftdatar_write("density",dtfil%fnameabo_app_den,dtset%iomode,hdr,&
439 3574 : crystal,ngfft,cplex1,nfft,nspden,rho_ptr,mpi_enreg,ebands=ebands)
440 :
441 3574 : if (dtset%positron/=0) then
442 11 : if (dtset%positron/=1) rho_ptr => electronpositron%rhor_ep
443 11 : if (dtset%positron==1) rho_ptr => rhor
444 11 : fname = trim(dtfil%fnameabo_app_den)//'_POSITRON'
445 11 : if (dtset%iomode == IO_MODE_ETSF) fname = strcat(fname, ".nc")
446 : call fftdatar_write("positron_density",fname,dtset%iomode,hdr,&
447 11 : crystal,ngfft,cplex1,nfft,nspden,rho_ptr,mpi_enreg,ebands=ebands)
448 : end if
449 : end if
450 :
451 1356 : else if (dtset%usewvl == 1 .and. dtset%prtden /= 0) then
452 : !if iomode == 2 then set all outputs to netcdf format
453 : !if iomode == 3 then set all outputs to ETSF format
454 0 : accessfil = 0
455 0 : if (dtset%iomode == IO_MODE_ETSF) accessfil = 3
456 0 : if (dtset%iomode == IO_MODE_MPI) accessfil = 4
457 0 : fform = fform_den
458 : ! Write wavelet DEN. Note however that this should be delegate to separated Bigdft routines.
459 : ! a lot of stuff written in outscf does not make sense if usewvl==0
460 : call ioarr(accessfil,rhor,dtset,etotal,fform,dtfil%fnameabo_app_den, &
461 0 : hdr,mpi_enreg,ngfft,cplex1,nfft,pawrhoij_dum,rdwr2,rdwrpaw0,wvl_den)
462 : end if ! if master
463 :
464 : !! MS - Printing of PAWDEN parallellised and several possible options included
465 : !We output the total electron density in the PAW case
466 : !this requires removing nhat from rhor and making PAW on-site corrections
467 6918 : if (pawprtden>0 .and. psps%usepaw==1) then
468 : ! pawprtden 1 --> output PAW valence density
469 : ! " 2 --> output PAW valence+core density
470 : ! " 3 --> output core, valence and full atomic protodensity
471 : ! " 4 --> options 1+3
472 : ! " 5 --> options 2+3
473 : ! " 6 --> output all individual PAW density contributions
474 1 : if (pawprtden/=3) then ! calc PAW valence density
475 4 : ABI_MALLOC(rhor_paw,(pawfgr%nfft,nspden))
476 3 : ABI_MALLOC(rhor_n_one,(pawfgr%nfft,nspden))
477 3 : ABI_MALLOC(rhor_nt_one,(pawfgr%nfft,nspden))
478 : ! If the communicator used for denfgr is kpt_comm, it is not compatible with paral_atom
479 1 : if (mpi_enreg%paral_kgb==0.and.my_natom/=natom) then
480 0 : my_natom_tmp=natom
481 0 : ABI_MALLOC(pawrhoij_all,(natom))
482 0 : call pawrhoij_nullify(pawrhoij_all)
483 : call pawrhoij_copy(pawrhoij,pawrhoij_all,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
484 0 : & keep_cplex=.false.,keep_qphase=.false.,keep_itypat=.false.,keep_nspden=.false.)
485 : else
486 1 : my_natom_tmp=my_natom
487 1 : pawrhoij_all => pawrhoij
488 : end if
489 1 : if (pawprtden/=6) then
490 : call denfgr(atindx1,gmet,comm_fft,my_natom_tmp,natom,nattyp,ngfft,nhat,dtset%nspinor,nsppol,nspden,&
491 : & ntypat,pawfgr,pawrad,pawrhoij_all,pawtab,prtvol,rhor,rhor_paw,rhor_n_one,&
492 : & rhor_nt_one,rprimd,dtset%typat,ucvol,xred,&
493 1 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
494 : else
495 : call denfgr(atindx1,gmet,comm_fft,my_natom_tmp,natom,nattyp,ngfft,nhat,dtset%nspinor,nsppol,nspden,&
496 : & ntypat,pawfgr,pawrad,pawrhoij_all,pawtab,prtvol,rhor,rhor_paw,rhor_n_one,&
497 : & rhor_nt_one,rprimd,dtset%typat,ucvol,xred,&
498 : & abs_n_tilde_nt_diff=nt_ntone_norm,znucl=dtset%znucl,&
499 0 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
500 : end if
501 1 : if (mpi_enreg%paral_kgb==0.and.my_natom/=natom) then
502 0 : call pawrhoij_free(pawrhoij_all)
503 0 : ABI_FREE(pawrhoij_all)
504 : end if
505 :
506 1 : if (prtvol>9) then ! Check normalisation
507 0 : norm = SUM(rhor_paw(:,1))*ucvol/PRODUCT(pawfgr%ngfft(1:3))
508 0 : call xmpi_sum(norm,comm_fft,ierr)
509 0 : write(msg,'(a,F8.4)') ' PAWDEN - NORM OF DENSITY: ',norm
510 0 : call wrtout(std_out, msg)
511 : end if
512 : end if
513 :
514 1 : if (pawprtden>1.AND.pawprtden<6) then ! We will need the core density
515 0 : ABI_MALLOC(rhor_paw_core,(pawfgr%nfft,nspden))
516 : call read_atomden(mpi_enreg,natom,pawfgr%nfft,pawfgr%ngfft,nspden,ntypat,rhor_paw_core,&
517 0 : & dtset%typat,rprimd,xred,prtvol,file_prefix='core ')
518 :
519 0 : if (prtvol>9) then ! Check normalisation
520 0 : norm = SUM(rhor_paw_core(:,1))*ucvol/PRODUCT(pawfgr%ngfft(1:3))
521 0 : call xmpi_sum(norm,comm_fft,ierr)
522 0 : write(msg,'(a,F8.4)') ' ATMDEN - NORM OF CORE DENSITY: ', norm
523 0 : call wrtout(std_out, msg)
524 : end if
525 : end if
526 :
527 0 : if (pawprtden>2.AND.pawprtden<6) then ! We will need the valence protodensity
528 0 : ABI_MALLOC(rhor_paw_val,(pawfgr%nfft,nspden))
529 : call read_atomden(mpi_enreg,natom,pawfgr%nfft,pawfgr%ngfft,nspden,ntypat,rhor_paw_val,&
530 0 : & dtset%typat,rprimd,xred,prtvol,file_prefix='valence')
531 :
532 0 : if (prtvol>9) then ! Check normalisation
533 0 : norm = SUM(rhor_paw_val(:,1))*ucvol/PRODUCT(pawfgr%ngfft(1:3))
534 0 : call xmpi_sum(norm,comm_fft,ierr)
535 0 : write(msg,'(a,F8.4)') ' ATMDEN - NORM OF VALENCE PROTODENSITY: ', norm
536 0 : call wrtout(std_out, msg)
537 : end if
538 : end if
539 :
540 1 : if (iwrite_fftdatar(mpi_enreg)) then
541 1 : if (pawprtden/=3) then
542 1 : if (pawprtden==2.or.pawprtden==5) rhor_paw = rhor_paw + rhor_paw_core
543 : ! PAWDEN
544 : call fftdatar_write("pawrhor",dtfil%fnameabo_app_pawden,dtset%iomode,hdr,&
545 1 : crystal,ngfft,cplex1,nfft,nspden,rhor_paw,mpi_enreg,ebands=ebands)
546 : end if
547 :
548 1 : if (pawprtden>2.AND.pawprtden<6) then
549 : ! ATMDEN_CORE
550 : call fftdatar_write("pawrhor_core",dtfil%fnameabo_app_atmden_core,dtset%iomode,hdr,&
551 0 : crystal,ngfft,cplex1,nfft,nspden,rhor_paw_core,mpi_enreg,ebands=ebands)
552 :
553 : ! valence protodensity. ATMDEN_VAL
554 : call fftdatar_write("pawrhor_val",dtfil%fnameabo_app_atmden_val,dtset%iomode,hdr,&
555 0 : crystal,ngfft,cplex1,nfft,nspden,rhor_paw_val,mpi_enreg,ebands=ebands)
556 :
557 : ! full protodensity. ATMDEN_FULL
558 0 : rhor_paw_val = rhor_paw_val + rhor_paw_core
559 : call fftdatar_write("pawrhor_full",dtfil%fnameabo_app_atmden_full,dtset%iomode,hdr,&
560 0 : crystal,ngfft,cplex1,nfft,nspden,rhor_paw_val,mpi_enreg,ebands=ebands)
561 : end if
562 :
563 1 : if (pawprtden==6) then ! Print all individual contributions to the density
564 : ! N_TILDE - N_HAT
565 : ! Use rhor_paw_val as temporary array
566 0 : if (.not.allocated(rhor_paw_val)) then
567 0 : ABI_MALLOC(rhor_paw_val,(pawfgr%nfft,nspden))
568 : end if
569 0 : rhor_paw_val = rhor - nhat
570 :
571 : call fftdatar_write("pawrhor_ntilde_minus_nhat",dtfil%fnameabo_app_n_tilde,dtset%iomode,hdr,&
572 0 : crystal,ngfft,cplex1,nfft,nspden,rhor_paw_val,mpi_enreg,ebands=ebands)
573 :
574 : ! N_ONSITE
575 : call fftdatar_write("pawrhor_n_one",dtfil%fnameabo_app_n_one,dtset%iomode,hdr,&
576 0 : crystal,ngfft,cplex1,nfft,nspden,rhor_n_one,mpi_enreg,ebands=ebands)
577 :
578 : ! N_TILDE_ONSITE
579 : call fftdatar_write("pawrhor_nt_one",dtfil%fnameabo_app_nt_one,dtset%iomode,hdr,&
580 0 : crystal,ngfft,cplex1,nfft,nspden,rhor_nt_one,mpi_enreg,ebands=ebands)
581 :
582 : end if ! All indivdual density cont.
583 : end if ! if master
584 :
585 1 : ABI_SFREE(rhor_paw)
586 1 : ABI_SFREE(rhor_paw_core)
587 1 : ABI_SFREE(rhor_paw_val)
588 1 : ABI_SFREE(rhor_n_one)
589 1 : ABI_SFREE(rhor_nt_one)
590 :
591 : end if ! if paw+pawprtden
592 :
593 6918 : call timab(1153,2,tsec)
594 :
595 6918 : call timab(1155,1,tsec)
596 :
597 : ! Output of VCLMB file
598 : ! The PAW correction has to be computed here (all processors contribute)
599 6918 : if (psps%usepaw > 0 .AND. dtset%prtvclmb>0) then
600 32 : nradint = 1000 ! radial integration grid density
601 128 : ABI_MALLOC(vpaw,(nfft,nspden))
602 938392 : vpaw(:,:)=zero
603 32 : if (me == master .and. my_natom > 0) then
604 8 : if (paw_an(1)%cplex > 1) then
605 0 : ABI_WARNING('cplex = 2 : complex hartree potential in PAW spheres. This is not coded yet. Imag part ignored')
606 : end if
607 : end if
608 :
609 88 : do ispden=1,nspden
610 : ! for points inside spheres, replace with full AE hartree potential.
611 : ! In principle the correction could be more subtle (not spherical)
612 164 : do iatom=1,my_natom
613 76 : iatom_tot=iatom;if (paral_atom) iatom_tot=mpi_enreg%my_atmtab(iatom)
614 76 : itypat=dtset%typat(iatom_tot)
615 :
616 228 : ABI_MALLOC(vh1spl,(paw_an(iatom)%mesh_size))
617 152 : ABI_MALLOC(vh1_corrector,(paw_an(iatom)%mesh_size))
618 228 : ABI_MALLOC(vh1_interp,(pawfgrtab(iatom)%nfgd))
619 152 : ABI_MALLOC(radii,(pawfgrtab(iatom)%nfgd))
620 228 : ABI_MALLOC(isort,(pawfgrtab(iatom)%nfgd))
621 : ! vh1 vht1 contain the spherical first moments of the Hartree potentials, so re-divide by Y_00 = sqrt(four_pi)
622 47040 : vh1_corrector(:) = (paw_an(iatom)%vh1(:,1,ispden)-paw_an(iatom)%vht1(:,1,ispden)) / sqrt(four_pi)
623 :
624 : ! get end point derivatives
625 76 : call bound_deriv(vh1_corrector, pawrad(itypat), pawrad(itypat)%mesh_size, yp1, ypn)
626 : ! spline the vh1 function
627 : ! NB for second argument of vh1: only first moment lm_size appears to be used
628 : ! NB2: vh1 can in principle be complex - not sure what to do with the imaginary part. Ignored for now.
629 76 : call spline(pawrad(itypat)%rad, vh1_corrector, paw_an(iatom)%mesh_size, yp1, ypn, vh1spl)
630 :
631 70184 : do ifgd = 1, pawfgrtab(iatom)%nfgd
632 : ! get radii for this point
633 70108 : isort(ifgd) = ifgd
634 280508 : radii(ifgd) = sqrt(sum(pawfgrtab(iatom)%rfgd(:,ifgd)**2))
635 : end do
636 :
637 76 : if (pawfgrtab(iatom)%nfgd/=0) then
638 : ! spline interpolate the vh1 value for current radii
639 74 : call sort_dp(pawfgrtab(iatom)%nfgd, radii, isort, tol12)
640 : call splint(pawrad(itypat)%mesh_size, pawrad(itypat)%rad, &
641 74 : & vh1_corrector, vh1spl, pawfgrtab(iatom)%nfgd, radii, vh1_interp, ierr)
642 : end if
643 :
644 70412 : norm=SUM(vh1_interp)*ucvol/PRODUCT(ngfft(1:3))
645 76 : call xmpi_sum(norm,comm_fft,ierr)
646 76 : write(msg,'(a,i6,a,E20.10)') ' sum of Hartree correction term on fft grid of atom : ', iatom, ' = ', norm
647 76 : call wrtout(std_out, msg)
648 :
649 76 : if (pawfgrtab(iatom)%nfgd/=0) then
650 : vpaw(pawfgrtab(iatom)%ifftsph(isort(1:pawfgrtab(iatom)%nfgd)),ispden) = &
651 : & vpaw(pawfgrtab(iatom)%ifftsph(isort(1:pawfgrtab(iatom)%nfgd)),ispden) + &
652 280506 : & vh1_interp(1:pawfgrtab(iatom)%nfgd)
653 : end if
654 :
655 : ! get integral of correction term in whole sphere
656 76 : ABI_FREE(radii)
657 76 : ABI_FREE(vh1_interp)
658 :
659 228 : ABI_MALLOC(radii,(nradint))
660 152 : ABI_MALLOC(vh1_interp,(nradint))
661 :
662 152 : ABI_MALLOC(vh1_integ,(nradint))
663 76 : dr = pawrad(itypat)%rad(paw_an(iatom)%mesh_size) / dble(nradint)
664 76076 : do ifgd = 1, nradint
665 76076 : radii(ifgd) = dble(ifgd-1)*dr
666 : end do
667 :
668 : ! spline interpolate the vh1 value for current radii
669 : call splint(pawrad(itypat)%mesh_size, pawrad(itypat)%rad, &
670 76 : & vh1_corrector, vh1spl, nradint, radii, vh1_interp, ierr)
671 :
672 76076 : do ifgd = 1, nradint
673 76076 : vh1_interp(ifgd) = vh1_interp(ifgd)*radii(ifgd)**2
674 : end do
675 :
676 76 : call simpson_int(nradint, dr, vh1_interp, vh1_integ)
677 76 : write(msg,'(a,i6,a,E20.10)') ' integral of Hartree correction term in sphere of atom: ', iatom, &
678 152 : & ' = ', vh1_integ(nradint)*four*pi
679 76 : call wrtout(std_out, msg)
680 :
681 76 : ABI_FREE(vh1spl)
682 76 : ABI_FREE(vh1_corrector)
683 76 : ABI_FREE(vh1_interp)
684 76 : ABI_FREE(vh1_integ)
685 76 : ABI_FREE(radii)
686 360 : ABI_FREE(isort)
687 : end do ! iatom
688 : end do !ispden
689 32 : call xmpi_sum_master(vpaw,master,mpi_enreg%comm_atom,ierr)
690 32 : if (.not.iwrite_fftdatar(mpi_enreg)) then
691 14 : ABI_FREE(vpaw)
692 : end if
693 : end if ! if paw - add all electron vhartree in spheres
694 :
695 6918 : call timab(1155,2,tsec)
696 :
697 6918 : if (iwrite_fftdatar(mpi_enreg)) then
698 :
699 5562 : call timab(1156,1,tsec)
700 :
701 : ! output the electron localization function ELF
702 5562 : if (dtset%prtelf/=0) then
703 : call fftdatar_write("elfr",dtfil%fnameabo_app_elf,dtset%iomode,hdr,&
704 4 : crystal,ngfft,cplex1,nfft,nspden,elfr,mpi_enreg,ebands=ebands)
705 :
706 4 : if (nspden==2)then
707 8 : ABI_MALLOC(elfr_up,(nfft,nspden))
708 108006 : elfr_up(:,:) = zero
709 54002 : do ifft=1,nfft
710 54002 : elfr_up(ifft,1) = elfr(ifft,2)
711 : end do
712 : ! ELF_UP
713 : call fftdatar_write("elfr_up",dtfil%fnameabo_app_elf_up,dtset%iomode,hdr,&
714 2 : crystal,ngfft,cplex1,nfft,nspden,elfr_up,mpi_enreg,ebands=ebands)
715 :
716 6 : ABI_MALLOC(elfr_down,(nfft,nspden))
717 108006 : elfr_down(:,:) = zero
718 54002 : do ifft=1,nfft
719 54002 : elfr_down(ifft,1) = elfr(ifft,3)
720 : end do
721 : ! ELF_DOWN'
722 : call fftdatar_write("elfr_down",dtfil%fnameabo_app_elf_down,dtset%iomode,hdr,&
723 2 : crystal,ngfft,cplex1,nfft,nspden,elfr_down,mpi_enreg,ebands=ebands)
724 :
725 2 : ABI_FREE(elfr_up)
726 2 : ABI_FREE(elfr_down)
727 : end if
728 : end if
729 :
730 5562 : call timab(1156,2,tsec)
731 5562 : call timab(1157,1,tsec)
732 :
733 : ! We output the gradient of density
734 5562 : if (dtset%prtgden/=0) then
735 :
736 : call fftdatar_write("grhor_1",dtfil%fnameabo_app_gden1,dtset%iomode,hdr,&
737 3 : crystal,ngfft,cplex1,nfft,nspden,grhor(:,:,1),mpi_enreg,ebands=ebands)
738 :
739 : call fftdatar_write("grhor_2",dtfil%fnameabo_app_gden2,dtset%iomode,hdr,&
740 3 : crystal,ngfft,cplex1,nfft,nspden,grhor(:,:,2),mpi_enreg,ebands=ebands)
741 :
742 : call fftdatar_write("grhor_3",dtfil%fnameabo_app_gden3,dtset%iomode,hdr,&
743 3 : crystal,ngfft,cplex1,nfft,nspden,grhor(:,:,3),mpi_enreg,ebands=ebands)
744 : end if
745 :
746 5562 : call timab(1157,2,tsec)
747 5562 : call timab(1158,1,tsec)
748 :
749 : ! We output the total kinetic energy density KDEN
750 5562 : if (dtset%prtkden/=0) then
751 : call fftdatar_write("kinedr",dtfil%fnameabo_app_kden,dtset%iomode,hdr,&
752 45 : crystal,ngfft,cplex1,nfft,nspden,taur,mpi_enreg,ebands=ebands)
753 : end if
754 :
755 5562 : call timab(1158,2,tsec)
756 5562 : call timab(1159,1,tsec)
757 :
758 :
759 : ! We output the Laplacian of density
760 5562 : if (dtset%prtlden/=0) then
761 : call fftdatar_write("laprhor",dtfil%fnameabo_app_lden,dtset%iomode,hdr,&
762 2 : crystal,ngfft,cplex1,nfft,nspden,lrhor,mpi_enreg,ebands=ebands)
763 : end if
764 :
765 5562 : call timab(1159,2,tsec)
766 5562 : call timab(1160,1,tsec)
767 :
768 : ! POT
769 5562 : if (dtset%prtpot>0) then
770 : call fftdatar_write("vtrial",dtfil%fnameabo_app_pot,dtset%iomode,hdr,&
771 58 : crystal,ngfft,cplex1,nfft,nspden,vtrial,mpi_enreg,ebands=ebands)
772 : end if
773 :
774 : ! EIG
775 5562 : if (dtset%prteig==2 .and. me == master) then
776 0 : fname=trim(dtfil%fnameabo_app_eig)//'.nc'
777 : call write_eig(eigen,e_fermie,fname,dtset%kptns,dtset%mband,dtset%nband,dtset%nkpt,dtset%nsppol,&
778 0 : results_gs%extfpmd_eshift) ! Optional arguments
779 : end if
780 :
781 5562 : call timab(1160,2,tsec)
782 5562 : call timab(1161,1,tsec)
783 :
784 5562 : if (dtset%prtgeo>0) then
785 45 : coordn=dtset%prtgeo
786 45 : call bonds_lgth_angles(coordn,dtfil%fnameabo_app_geo,natom,psps%ntypat, rprimd,dtset%typat,xred,dtset%znucl)
787 : end if
788 :
789 5562 : if (dtset%prtcif > 0) then
790 : call prt_cif(dtset%brvltt, dtfil%fnameabo_app_cif, natom, dtset%nsym, dtset%ntypat, rprimd, &
791 2 : dtset%spgaxor, dtset%spgroup, dtset%spgorig, dtset%symrel, dtset%tnons, dtset%typat, xred, dtset%znucl)
792 : end if
793 :
794 5562 : call timab(1161,2,tsec)
795 5562 : call timab(1162,1,tsec)
796 :
797 : ! STM
798 5562 : if (dtset%prtstm/=0) then
799 3 : call fftdatar_write("stm",dtfil%fnameabo_app_stm,dtset%iomode,hdr,crystal,ngfft,cplex1,nfft,nspden,rhor,mpi_enreg,ebands=ebands)
800 : end if
801 :
802 5562 : call timab(1162,2,tsec)
803 5562 : call timab(1163,1,tsec)
804 :
805 5562 : if (dtset%prt1dm>0) then
806 37 : call out1dm(dtfil%fnameabo_app_1dm,mpi_enreg,natom,nfft,ngfft,nspden,psps%ntypat,rhor,rprimd,dtset%typat,ucvol,vtrial,xred,dtset%znucl)
807 : end if
808 :
809 5562 : call timab(1163,2,tsec)
810 5562 : call timab(1164,1,tsec)
811 :
812 : ! VHA
813 5562 : if (dtset%prtvha>0) then
814 148 : ABI_MALLOC(vwork,(nfft,nspden))
815 86 : do ispden=1,nspden
816 688958 : vwork(:,ispden)=vhartr(:)
817 : end do
818 :
819 37 : call fftdatar_write("vhartree",dtfil%fnameabo_app_vha,dtset%iomode,hdr,crystal,ngfft,cplex1,nfft,nspden,vwork,mpi_enreg,ebands=ebands)
820 37 : ABI_FREE(vwork)
821 : end if
822 :
823 : ! VPSP
824 5562 : if (dtset%prtvpsp>0) then
825 132 : ABI_MALLOC(vwork,(nfft,nspden))
826 78 : do ispden=1,nspden
827 686030 : vwork(:,ispden)=vpsp(:)
828 : end do
829 :
830 33 : call fftdatar_write("vpsp",dtfil%fnameabo_app_vpsp,dtset%iomode,hdr,crystal,ngfft,cplex1,nfft,nspden,vwork,mpi_enreg,ebands=ebands)
831 :
832 33 : ABI_FREE(vwork)
833 : end if
834 :
835 : ! VCouLoMB
836 5562 : if (dtset%prtvclmb>0) then
837 :
838 140 : ABI_MALLOC(vwork,(nfft,nspden))
839 82 : do ispden=1,nspden
840 687314 : vwork(:,ispden)=vpsp(:)+vhartr(:)
841 : end do
842 35 : if (psps%usepaw==1) then
843 48 : do ispden=1,nspden
844 454512 : vwork(:,ispden)=vwork(:,ispden)+vpaw(:,ispden)
845 : end do
846 18 : ABI_FREE(vpaw)
847 : end if
848 :
849 : call fftdatar_write("vhartree_vloc",dtfil%fnameabo_app_vclmb,dtset%iomode,hdr,&
850 35 : crystal,ngfft,cplex1,nfft,nspden,vwork,mpi_enreg,ebands=ebands)
851 :
852 : !TODO: find out why this combination of calls with fftdatar_write then out1dm fails on buda with 4 mpi-fft procs (np_spkpt 1).
853 : ! For the moment comment it out. Only DS2 of mpiio test 27 fails
854 : ! call out1dm(dtfil%fnameabo_app_vclmb_1dm,mpi_enreg,natom,nfft,ngfft,nspden,psps%ntypat,&
855 : !& rhor,rprimd,dtset%typat,ucvol,vwork,xred,dtset%znucl)
856 :
857 : ! TODO: add TEM phase with CE = (2 pi / lambda) (E+E0)/(E(E+2E0)) from p.49 of RE Dunin Borkowski 2004 encyclopedia of nanoscience volume 3 pp 41-99
858 : ! where E is energy of electron, E0 rest mass, lambda the relativistic wavelength
859 : ! values of CE at 200 300 and 1000 kV: 7.29e6 6.53e6 5.39e6 rad / V / m
860 : ! vertical integral of vclmb * c / ngfft(3) / cross sectional area factor (= sin(gamma))
861 : ! * Bohr_Ang * 1.0e-10* Ha_eV to get to SI
862 : ! * CE factor above
863 : ! should be done for each plane perpendicular to the axes...
864 35 : ABI_FREE(vwork)
865 : end if ! prtvclmb
866 :
867 :
868 : ! VHXC
869 5562 : if (dtset%prtvhxc>0) then
870 144 : ABI_MALLOC(vwork,(nfft,nspden))
871 84 : do ispden=1,nspden
872 694876 : vwork(:,ispden)=vhartr(:)+vxc(:,ispden)
873 : end do
874 :
875 : call fftdatar_write("vhxc",dtfil%fnameabo_app_vhxc,dtset%iomode,hdr,&
876 36 : crystal,ngfft,cplex1,nfft,nspden,vwork,mpi_enreg,ebands=ebands)
877 36 : ABI_FREE(vwork)
878 : end if
879 :
880 : ! VXC
881 5562 : if (dtset%prtvxc>0) then
882 : call fftdatar_write("exchange_correlation_potential",dtfil%fnameabo_app_vxc,dtset%iomode,hdr,&
883 34 : crystal,ngfft,cplex1,nfft,nspden,vxc,mpi_enreg,ebands=ebands)
884 : end if
885 :
886 5562 : call timab(1164,2,tsec)
887 :
888 : end if ! if iwrite_fftdatar
889 :
890 6918 : call timab(1165,1,tsec)
891 :
892 : !Generate DOS using the tetrahedron method or using Gaussians
893 : !FIXME: Should centralize all calculations of DOS here in outscfcv
894 6918 : if (dtset%prtdos>=2.or.dtset%pawfatbnd>0) then
895 49 : call dos%init(dtset, psps, pawtab)
896 :
897 49 : if (dos%partial_dos_flag>=1 .or. dos%fatbands_flag==1)then
898 : ! Generate fractions for partial DOSs if needed partial_dos 1,2,3,4 give different decompositions
899 45 : collect = 1 !; if (psps%usepaw==1 .and. dos%partial_dos_flag /= 2) collect = 0
900 45 : if ((psps%usepaw==0.or.dtset%pawprtdos/=2) .and. dos%partial_dos_flag>=1) then
901 42 : call dos%partial_dos_fractions(crystal,dtset,eigen,occ,npwarr,kg,cg,mcg,collect,mpi_enreg)
902 : end if
903 :
904 45 : if (psps%usepaw==1 .and. dos%partial_dos_flag /= 2) then
905 : ! TODO: update partial_dos_fractions_paw for extra atoms - no PAW contribution normally, but check bounds and so on.
906 24 : call dos%partial_dos_fractions_paw(cprj,dimcprj,dtset,mcprj,mkmem,mpi_enreg,pawrad,pawtab)
907 : end if
908 :
909 : else
910 370 : dos%fractions(:,:,:,1)=one
911 : end if
912 :
913 : ! Here, print out fatbands for the k-points given in file appended _FATBANDS
914 49 : if (me == master .and. dtset%pawfatbnd>0 .and. dos%fatbands_flag==1) then
915 1 : call dos%prtfatbands(dtset,ebands,dtfil%fnameabo_app_fatbands,dtset%pawfatbnd,pawtab)
916 : end if
917 :
918 : ! Here, computation and output of DOS and partial DOS _DOS
919 49 : if (dos%fatbands_flag == 0 .and. dos%prtdos /= 4) then
920 46 : call dos%calcnwrite(dtset,crystal,ebands,dtfil%fnameabo_app_dos,comm)
921 : end if
922 :
923 : ! Write netcdf file with dos% results.
924 49 : if (me == master) then
925 28 : fname = trim(dtfil%filnam_ds(4))//'_FATBANDS.nc'
926 28 : NCF_CHECK(nctk_open_create(ncid, fname, xmpi_comm_self))
927 28 : call dos%ncwrite(crystal, ebands, hdr, dtset, psps, pawtab, ncid)
928 28 : NCF_CHECK(nf90_close(ncid))
929 : end if
930 :
931 : !TODO: do not free dos here, but use the fractions below in calcdenmagsph
932 49 : call dos%free()
933 : end if ! prtdos > 1
934 :
935 6918 : call timab(1165,2,tsec)
936 6918 : call timab(1166,1,tsec)
937 :
938 : !Output of integrated density inside atomic spheres
939 27672 : if ((dtset%prtdensph==1.and.dtset%usewvl==0) .or. sum(abs(dtset%hspinfield)) > tol10) then
940 27672 : ABI_MALLOC(intgden, (nspden, natom))
941 :
942 : call calcdenmagsph(mpi_enreg,natom,nfft,ngfft,nspden,&
943 6918 : ntypat,dtset%ratsm,dtset%ratsph,rhor,rprimd,dtset%typat,xred,1,cplex1,dtset%qgbt,dtset%use_gbt,intgden=intgden,rhomag=rhomag)
944 :
945 : ! for rhomag:
946 : ! in collinear case component 1 is total density and 2 is _magnetization_ up-down
947 : ! in non collinear case component 1 is total density, and 2:4 are the magnetization vector
948 :
949 6918 : if (dtset%prtdensph==1.and.dtset%usewvl==0) then
950 16395 : if(all(dtset%constraint_kind(:)==0))then
951 6864 : call prtdenmagsph(cplex1,intgden,natom,nspden,ntypat,units,1,dtset%qgbt,dtset%ratsm,dtset%ratsph,rhomag,dtset%typat,dtset%znucl,dtset%spinaxis)
952 : else
953 54 : call prtdenmagsph(cplex1,intgden,natom,nspden,ntypat,units,1,dtset%qgbt,dtset%ratsm,dtset%ratsph,rhomag,dtset%typat,dtset%znucl,dtset%spinaxis,dtset%ziontypat)
954 : endif
955 16395 : if(any(dtset%constraint_kind(:)/=0))then
956 54 : call prtdenmagsph(cplex1,intgres,natom,nspden,ntypat,units,21,dtset%qgbt,dtset%ratsm,dtset%ratsph,rhomag,dtset%typat,dtset%znucl,dtset%spinaxis)
957 : endif
958 : end if !end prtdensph==1 .and. usewvl==0
959 :
960 : !!!!!!!!!!!!!!!!!!!!!!!!if prt_lorbmag value is equal 1 and the calculations are noncollinear then the local orbital magnetic moments are calculated
961 6918 : if (dtset%prt_lorbmag==1) then
962 :
963 1 : if ((dtset%nspinor .ne. 2) .and. (dtset%nspden .ne.4)) then
964 0 : write (msg,'(a)')" "
965 0 : call wrtout(units, msg)
966 0 : write (msg,'(a)')"WARNING*"
967 0 : call wrtout(units, msg)
968 0 : write (msg,'(a)')"prt_lorbmag=1, To calculate orbital magnetisation, calculations need to be noncollinear"
969 0 : call wrtout(units, msg)
970 : else
971 1 : if (dtset%usepawu .ne. 0)then
972 : call loc_orbmom_cal(1,0,dmatdum,0,0,indsym,my_natom,dtset%natom,dtset%natpawu,&
973 : & dtset%nspinor,dtset%nsppol,dtset%nsym,dtset%ntypat,paw_ij,pawang,pawrad,dtset%pawprtvol,&
974 : & pawrhoij,pawtab,dtset%spinat,dtset%symafm,dtset%typat,0,dtset%usepawu,dtset%znucl,&
975 1 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
976 : else
977 0 : write (msg,'(a)')" "
978 0 : call wrtout(units, msg)
979 0 : write (msg,'(a)')"WARNING*"
980 0 : call wrtout(units, msg)
981 0 : write (msg,'(a)')"prt_lorbmag=1, To calculate orbital magnetisation LDA+U calculations should be activated"
982 0 : call wrtout(units, msg)
983 : end if
984 : endif
985 : end if
986 :
987 27672 : if (sum(abs(dtset%hspinfield)) > tol10) then
988 12 : if(nspden==2)then
989 2 : e_hspinfield = -half*rhomag(1,2)*dtset%hspinfield(3)
990 2 : write (msg, "(a,E20.10,a)") " Collinear magnetization ", rhomag(1,2), &
991 4 : " (in # of spins, without 1/2 for magnetic moment) "
992 2 : call wrtout(units, msg)
993 10 : else if(nspden==4)then
994 : e_hspinfield = -half * (dtset%hspinfield(1)*rhomag(1,2)& ! x
995 : & +dtset%hspinfield(2)*rhomag(1,3)& ! y
996 10 : & +dtset%hspinfield(3)*rhomag(1,4)) ! z
997 10 : write (msg, "(a,3E20.10,a)") " Magnetization vector ", rhomag(1,2:4), &
998 20 : & " (in # of spins, without 1/2 for magnetic moment) "
999 10 : call wrtout(units, msg)
1000 : end if
1001 : !TODO: this quantity should also be calculated in rhotov, and stored in
1002 : ! results_gs%energies%e_hspinfield, but for the moment it comes out 0
1003 12 : write (msg, "(a,E20.10,a)") " Spin magnetic energy -m.B = ", e_hspinfield, " Ha"
1004 12 : call wrtout(units, msg)
1005 : end if
1006 : end if ! end if prtdensph or magnetic field
1007 :
1008 6918 : call timab(1166,2,tsec)
1009 6918 : call timab(1167,1,tsec)
1010 :
1011 6918 : if (dtset%magconon /= 0) then
1012 : ! calculate final value of terms for magnetic constraint: "energy" term, lagrange multiplier term, and atomic contributions
1013 : call mag_penalty_e(dtset%magconon,dtset%magcon_lambda,mpi_enreg,&
1014 5 : & natom,nfft,ngfft,nspden,ntypat,dtset%ratsm,dtset%ratsph,rhor,rprimd,dtset%spinat,dtset%typat,xred,dtset%znucl,dtset%qgbt,dtset%use_gbt,dtset%spinaxis)
1015 : end if
1016 :
1017 6918 : call timab(1167,2,tsec)
1018 6918 : call timab(1168,1,tsec)
1019 :
1020 : !If PAW, provide additional outputs
1021 6918 : if (psps%usepaw==1) then
1022 : ! Output of compensation charge
1023 1500 : if (dtset%nstep>0.or.dtfil%ireadwf/=0) then
1024 1499 : write(msg, '(4a)' )ch10,' PAW TEST:',ch10,&
1025 2998 : & ' ==== Compensation charge inside spheres ============'
1026 1499 : if (compch_sph>-1.d4.and.compch_fft>-1.d4) &
1027 1479 : & write(msg, '(3a)' ) trim(msg),ch10,' The following values must be close to each other ...'
1028 2998 : if (compch_sph>-1.d4) write(msg, '(3a,f22.15)' ) trim(msg),ch10,&
1029 2998 : & ' Compensation charge over spherical meshes = ',compch_sph
1030 1499 : if (compch_fft>-1.d4) then
1031 1479 : if (pawfgr%usefinegrid==1) then
1032 1410 : write(msg, '(3a,f22.15)' ) trim(msg),ch10,&
1033 2820 : & ' Compensation charge over fine fft grid = ',compch_fft
1034 : else
1035 69 : write(msg, '(3a,f22.15)' ) trim(msg),ch10,&
1036 138 : & ' Compensation charge over fft grid = ',compch_fft
1037 : end if
1038 : end if
1039 1499 : call wrtout(units, msg)
1040 : end if
1041 : ! Output of pseudopotential strength Dij and augmentation occupancies Rhoij
1042 : call pawprt(dtset,my_natom,paw_ij,pawrhoij,pawtab,&
1043 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
1044 1500 : & electronpositron=electronpositron)
1045 : end if
1046 :
1047 6918 : call timab(1168,2,tsec)
1048 6918 : call timab(1169,1,tsec)
1049 :
1050 :
1051 : !PAW + output for optical conductivity _OPT and _OPT2
1052 6918 : if (psps%usepaw==1.and.prtnabla>0) then
1053 28 : if (prtnabla==1.or.prtnabla==2) then
1054 : call optics_paw(atindx1,cg,cprj,dimcprj,dtfil,dtset,eigen,gprimd,hdr,kg,&
1055 : & mband,mcg,mcprj,mkmem,mpi_enreg,mpsang,mpw,natom,nkpt,npwarr,nsppol,pawang,&
1056 26 : & pawrad,pawrhoij,pawtab,psps%znuclpsp)
1057 : end if
1058 28 : if (prtnabla==2.or.prtnabla==3) then
1059 3 : if(associated(rcpaw)) then
1060 1 : if(dtset%nspinor==2) then
1061 1 : ABI_MALLOC(rcpaw_dirac,)
1062 : call rcpaw_init(rcpaw_dirac,dtset,psps%filpsp,pawrad,pawtab,dtset%ntypat,1,.True.,&
1063 1 : & my_natom,mpi_enreg%comm_atom,mpi_enreg%my_atmtab)
1064 1 : option=0;if (dtset%iscf>0.and.dtset%iscf<10) option=1
1065 1 : if(dtset%extfpmd_pawsph>0) then
1066 0 : ABI_ERROR('Extfpmd pawsph > 0 with nspinor=2 is work in progress for optics')
1067 : endif
1068 2 : do itypat=1,rcpaw_dirac%ntypat
1069 6 : do iln=1,rcpaw_dirac%atm(itypat)%ln_size
1070 4 : l=rcpaw_dirac%atm(itypat)%indln(1,iln)
1071 4 : n=rcpaw_dirac%atm(itypat)%indln(2,iln)
1072 10 : do iln2=1,rcpaw%atm(itypat)%ln_size
1073 9 : ll=rcpaw%atm(itypat)%indln(1,iln2)
1074 9 : nn=rcpaw%atm(itypat)%indln(2,iln2)
1075 9 : if(nn==n.and.ll==l) then
1076 : rcpaw_dirac%atm(itypat)%occ(iln,1)=rcpaw_dirac%atm(itypat)%occ(iln,1)*&
1077 4 : rcpaw%atm(itypat)%occ(iln2,1)/rcpaw%atm(itypat)%max_occ(iln2,1)
1078 4 : exit
1079 : endif
1080 : enddo
1081 : enddo
1082 : enddo
1083 : call paw_relax_core(pawtab,pawrad,pawang,pawrhoij,dtset%ntypat,rcpaw_dirac,psps,dtset,&
1084 : & 1,dtset%pawnzlm,option,ucvol,paw_an,my_natom,&
1085 1 : & mpi_enreg%my_atmtab,mpi_enreg%comm_atom)
1086 2 : do itypat=1,rcpaw%ntypat
1087 7 : rcpaw_dirac%atm(itypat)%eig=rcpaw_dirac%atm(itypat)%eig+rcpaw%atm(itypat)%eigshift
1088 : enddo
1089 : call optics_paw_core(atindx1,cprj,dimcprj,dtfil,dtset,eigen,psps%filpsp,hdr,&
1090 : & mband,mcprj,mkmem,mpi_enreg,mpsang,natom,nkpt,nsppol,pawang,pawrad,pawrhoij,pawtab,&
1091 1 : & psps%znuclpsp,rcpaw=rcpaw_dirac)
1092 1 : call rcpaw_destroy(rcpaw_dirac)
1093 1 : ABI_FREE(rcpaw_dirac)
1094 : else
1095 : call optics_paw_core(atindx1,cprj,dimcprj,dtfil,dtset,eigen,psps%filpsp,hdr,&
1096 : & mband,mcprj,mkmem,mpi_enreg,mpsang,natom,nkpt,nsppol,pawang,pawrad,pawrhoij,pawtab,&
1097 0 : & psps%znuclpsp,rcpaw=rcpaw)
1098 : endif
1099 : else
1100 : call optics_paw_core(atindx1,cprj,dimcprj,dtfil,dtset,eigen,psps%filpsp,hdr,&
1101 : & mband,mcprj,mkmem,mpi_enreg,mpsang,natom,nkpt,nsppol,pawang,pawrad,pawrhoij,pawtab,&
1102 2 : & psps%znuclpsp)
1103 : endif
1104 : end if
1105 : end if
1106 6918 : if (prtnabla<0) then
1107 : ! TODO: This routine is not tested but it's used in production.
1108 : call optics_vloc(cg,dtfil,dtset,eigen,gprimd,hdr,kg,&
1109 0 : & mband,mcg,mkmem,mpi_enreg,mpw,nkpt,npwarr,nsppol)
1110 : end if
1111 :
1112 6918 : call timab(1169,2,tsec)
1113 6918 : call timab(1170,1,tsec)
1114 :
1115 : !Optionally provide output for AE wavefunctions (only for PAW)
1116 6918 : if (psps%usepaw==1 .and. dtset%pawprtwf>=1) then
1117 5 : ABI_MALLOC(ps_norms,(nsppol,nkpt,mband))
1118 :
1119 : call pawmkaewf(dtset,crystal,ebands,my_natom,mpw,mband,mcg,mcprj,nkpt,mkmem,nsppol,Dtset%nband,&
1120 : & Dtset%istwfk,npwarr,Dtset%kptns,Dtset%ngfftdg,kg,dimcprj,pawfgrtab,&
1121 : & Pawrad,Pawtab,Hdr,Dtfil,cg,Cprj,&
1122 : & MPI_enreg,ierr,pseudo_norms=ps_norms,set_k=dtset%pawprt_k,set_band=dtset%pawprt_b,&
1123 1 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1124 :
1125 1 : if (dtset%pawprt_b==0) then
1126 1 : fname = strcat(dtfil%filnam_ds(4), '_PAWSTAT')
1127 1 : if (open_file(fname, msg,newunit=tmp_unt,status='unknown',form='formatted') /= 0) then
1128 0 : ABI_ERROR(msg)
1129 : end if
1130 1 : write(tmp_unt,'(5a)') '# This file contains the statistics on the cancellation of',ch10,&
1131 1 : & '# the onsite pseudo component of the all-electron wavefunction',ch10,&
1132 2 : & '# with the plane wave part'
1133 1 : ii = 0
1134 2 : do isppol=1,nsppol
1135 1 : write(tmp_unt,'(a,i0)') '# isppol = ',isppol
1136 3 : do ikpt=1,nkpt
1137 1 : write(tmp_unt,'(a,i0)') '# ikpt = ',ikpt
1138 1 : write(tmp_unt,'(a)') '# band norm'
1139 1 : occ_norm = zero; unocc_norm = zero; nocc = 0
1140 7 : do iband=1,dtset%nband(ikpt + (isppol-1)*nkpt)
1141 6 : ii = ii + 1
1142 6 : write(tmp_unt,'(i8,ES16.6)') iband,ps_norms(isppol,ikpt,iband)
1143 7 : if (abs(occ(ii)) <= tol16) then
1144 0 : unocc_norm = unocc_norm + ps_norms(isppol,ikpt,iband)
1145 : else
1146 6 : occ_norm = occ_norm + ps_norms(isppol,ikpt,iband)
1147 6 : nocc = nocc + 1
1148 : end if
1149 : end do
1150 2 : if(mband/=nocc)then
1151 0 : write(tmp_unt,'(2(a,ES16.6))') '# occ average: ',occ_norm/real(nocc),&
1152 0 : & ' unocc average: ',unocc_norm/real(mband-nocc)
1153 : else
1154 1 : write(tmp_unt,'(2(a,ES16.6))') '# occ average: ',occ_norm/real(nocc)
1155 : end if
1156 : end do
1157 : end do
1158 1 : close(tmp_unt)
1159 : end if
1160 1 : ABI_FREE(ps_norms)
1161 : end if
1162 :
1163 6918 : call timab(1170,2,tsec)
1164 6918 : call timab(1171,1,tsec)
1165 :
1166 6918 : if(dtset%plowan_compute>0 .and. dtset%plowan_compute<10) then
1167 0 : write(msg,'(2a,i3)') ch10,&
1168 0 : & ' ====================================================================================== '
1169 0 : call wrtout(units, msg)
1170 0 : write(msg,'(2a,i3)') ch10,&
1171 0 : & ' == Start computation of Projected Local Orbitals Wannier functions == ',dtset%nbandkss
1172 0 : call wrtout(units, msg)
1173 :
1174 : ! == compute chipsi
1175 :
1176 : call init_plowannier(dtset%plowan_bandf,dtset%plowan_bandi,dtset%plowan_compute,&
1177 : & dtset%plowan_iatom,dtset%plowan_it,dtset%plowan_lcalc,dtset%plowan_natom,&
1178 : & dtset%plowan_nbl,dtset%plowan_nt,dtset%plowan_projcalc,dtset%acell_orig,&
1179 0 : & dtset%kptns,sum(dtset%plowan_nbl),dtset%nimage,dtset%nkpt,dtset%nspinor,dtset%nsppol,dtset%wtk,dtset%dmft_t2g,wan)
1180 : call compute_coeff_plowannier(crystal,cprj,dimcprj,dtset,eigen,e_fermie,&
1181 0 : & mpi_enreg,occ,wan,pawtab,psps,usecprj,dtfil%unpaw,pawrad,dtfil)
1182 0 : if (me==master) then
1183 0 : call print_plowannier(wan)
1184 : endif
1185 0 : call destroy_plowannier(wan)
1186 : end if
1187 :
1188 6918 : call timab(1171,2,tsec)
1189 6918 : call timab(1172,1,tsec)
1190 :
1191 : !Optionally provide output for the GW part of ABINIT
1192 6918 : if (dtset%nbandkss/=0) then
1193 : ! Use DMFT to compute wannier function for cRPA calculation.
1194 6 : if(dtset%usedmft==1) then
1195 4 : write(msg,'(2a,i3)') ch10,&
1196 8 : & ' Warning: Chipsi are orthonormalized in the DMFT code because nbandkss is used, with the value ',dtset%nbandkss
1197 4 : call wrtout(std_out, msg)
1198 : call init_dmft(crystal,dmatpawu(:,:,:,:),dtset,e_fermie,dtfil%filctqmcdatain,dtfil%filselfin, &
1199 4 : & dtfil%filnam_ds(3),dtfil%fnameabo_app,dtfil%ireadctqmcdata,dtfil%ireadself,paw_dmft,pawtab(:))
1200 4 : call print_dmft(paw_dmft,dtset%pawprtvol)
1201 :
1202 : ! == compute chipsi
1203 4 : call init_oper(paw_dmft,dft_occup,opt_ksloc=2)
1204 :
1205 : call datafordmft(cg(:,:),cprj(:,:),crystal,dft_occup,dimcprj,dtset,eigen(:),dtset%mband,mcg,&
1206 4 : & mpi_enreg,dtset%nspinor,occ(:),paw_dmft,paw_ij,pawtab(:),usecprj,dtset%nbandkss)
1207 :
1208 4 : opt_imagonly=0
1209 4 : if(paw_dmft%dmft_solv>=5) opt_imagonly=1
1210 :
1211 :
1212 : ! Compute k-resolved spectral function in DMFT.
1213 4 : if(dtset%dmft_kspectralfunc==1) then
1214 : ! Initialize self on real axis
1215 4 : call initialize_self(selfr,paw_dmft,wtype='real')
1216 :
1217 4 : opt_moments = 0
1218 4 : if (paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7) opt_moments = 1
1219 :
1220 : ! Initialize self on imag axis
1221 4 : call initialize_self(self,paw_dmft,opt_moments=opt_moments)
1222 :
1223 : ! Initialize green on real axis
1224 4 : call init_green(greenr,paw_dmft,opt_oper_ksloc=3,wtype='real')
1225 :
1226 : ! Read self energy in imag. Matsubara freq (for double counting and asymptotic value)
1227 4 : call rw_self(self,paw_dmft,prtopt=5,opt_rw=1,opt_stop=1)
1228 :
1229 16 : ABI_MALLOC(opt_selflimit,(paw_dmft%natom))
1230 4 : call init_matlu(paw_dmft%natom,paw_dmft%nspinor,paw_dmft%nsppol,paw_dmft%lpawu(:),opt_selflimit(:))
1231 :
1232 4 : if (opt_moments == 0) then
1233 4 : call copy_matlu(self%oper(self%nw)%matlu(:),opt_selflimit(:),paw_dmft%natom,opt_re=1)
1234 : else
1235 0 : call copy_matlu(self%moments(1)%matlu(:),opt_selflimit(:),paw_dmft%natom)
1236 : end if ! moments
1237 :
1238 : ! Read self energy on real axis obtained from Maxent
1239 : call rw_self(selfr,paw_dmft,prtopt=5,opt_rw=1,opt_imagonly=opt_imagonly, &
1240 4 : & opt_selflimit=opt_selflimit(:),opt_hdc=self%hdc%matlu(:),opt_maxent=1)
1241 :
1242 4 : call destroy_matlu(opt_selflimit(:),paw_dmft%natom)
1243 8 : ABI_FREE(opt_selflimit)
1244 :
1245 : ! Check: from self on real axis, recompute self on Imaginary axis.
1246 4 : call selfreal2imag_self(selfr,self,paw_dmft%filapp,paw_dmft)
1247 :
1248 : ! paw_dmft%fermie=hdr%fermie ! for tests
1249 4 : write(std_out,*) " Fermi level is",paw_dmft%fermie
1250 :
1251 : ! For the DFT BS: use opt_self=0 and fermie=fermie_dft
1252 :
1253 : ! Compute Green's function on real axis
1254 4 : call compute_green(greenr,paw_dmft,1,selfr,opt_self=1,opt_nonxsum=0)
1255 :
1256 : !write(6,*) "compute green done"
1257 4 : if(me==master) then
1258 1 : if(dtset%kptopt<0) then
1259 : ! k-resolved Spectral function
1260 1 : call print_green("from_realaxisself",greenr,5,paw_dmft,opt_wt=1)
1261 : else
1262 : ! DOS Calculation
1263 0 : call print_green("from_realaxisself",greenr,4,paw_dmft,opt_wt=1)
1264 : endif
1265 : !write(6,*) "print green done"
1266 : endif
1267 :
1268 4 : call destroy_green(greenr)
1269 4 : call destroy_self(selfr)
1270 4 : call destroy_self(self)
1271 : endif
1272 4 : call destroy_dmft(paw_dmft)
1273 4 : call destroy_oper(dft_occup)
1274 : end if
1275 :
1276 : call outkss(crystal,dtfil,dtset,ecut,gmet,gprimd,hdr,&
1277 : & dtset%kssform,mband,mcg,mcprj,mgfftc,mkmem,mpi_enreg,mpsang,mpw,natom,natom,&
1278 : & nfft,nkpt,npwarr,nspden,nsppol,nsym,psps%ntypat,occ,pawtab,pawfgr,paw_ij,&
1279 6 : & prtvol,psps,rprimd,vtrial,xred,cg,usecprj,cprj,eigen,ierr)
1280 6 : if (ierr/=0) then
1281 0 : ABI_WARNING("outkss returned a non zero status error, check log")
1282 : end if
1283 : end if
1284 :
1285 6918 : call timab(1172,2,tsec) ! outscfcv(gw)
1286 :
1287 6918 : if (electronpositron_calctype(electronpositron)/=0) then
1288 :
1289 : ! Optionally provide output for positron life time calculation
1290 17 : call timab(1173,1,tsec)
1291 : call poslifetime(dtset,electronpositron,gprimd,my_natom,&
1292 : & mpi_enreg,n3xccc,nfft,ngfft,nhat,1,pawang,&
1293 : & pawrad,pawrhoij,pawtab,rate_dum,rate_dum2,&
1294 17 : & rhor,ucvol,xccc3d)
1295 17 : call timab(1173,2,tsec)
1296 :
1297 : ! Optionally provide output for momentum distribution of annihilation radiation
1298 17 : if (dtset%posdoppler>0) then
1299 3 : call timab(1174,1,tsec)
1300 : call posdoppler(cg,cprj,crystal,dimcprj,dtfil,dtset,electronpositron,psps%filpsp,&
1301 : & kg,mcg,mcprj,mpi_enreg,my_natom,n3xccc,nfft,ngfft,nhat,npwarr,&
1302 3 : & occ,pawang,pawrad,pawrhoij,pawtab,rhor,xccc3d)
1303 3 : call timab(1174,2,tsec)
1304 : end if
1305 : end if
1306 :
1307 : !Optionally provide output for WanT
1308 6918 : if (dtset%prtwant==1) then
1309 1 : call timab(1175,1,tsec)
1310 : ! WARNING: mpi_enreg not used --> MPI is not supported
1311 1 : call outwant(dtset,eigen,cg,kg,npwarr,mband,mcg,nkpt,nsppol,mkmem,mpw,dtset%prtwant)
1312 1 : call timab(1175,2,tsec)
1313 : end if
1314 :
1315 : !Optionally provide output for electric field gradient calculation
1316 6918 : if (dtset%nucefg > 0) then
1317 :
1318 9 : if(allocated(efg)) then
1319 0 : ABI_FREE(efg)
1320 : end if
1321 27 : ABI_MALLOC(efg,(3,3,natom))
1322 542 : efg=zero
1323 9 : call timab(1176,1,tsec)
1324 : call calc_efg(efg,mpi_enreg,my_natom,natom,nfft,ngfft,nhat,nspden,dtset%nsym,dtset%nucefg,&
1325 : & ntypat,paw_an,pawang,pawrad,pawrhoij,pawtab,&
1326 : & dtset%ptcharge,dtset%quadmom,rhor,rprimd,dtset%symrel,&
1327 : & dtset%tnons,dtset%typat,ucvol,psps%usepaw,xred,psps%zionpsp,&
1328 9 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1329 9 : call timab(1176,2,tsec)
1330 : end if
1331 :
1332 : !Optionally provide output for Fermi-contact term at nuclear positions
1333 6918 : if (dtset%nucfc > 0) then
1334 5 : call timab(1177,1,tsec)
1335 : call calc_fc(my_natom,natom,nspden,ntypat,pawrad,pawrhoij,pawtab,dtset%typat,psps%usepaw,&
1336 5 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1337 5 : call timab(1177,2,tsec)
1338 : end if
1339 :
1340 : ! Output electron bands.
1341 6918 : if (me == master .and. dtset%tfkinfunc==0) then
1342 5407 : call timab(1178,1,tsec)
1343 5407 : if (size(dtset%kptbounds, dim=2) > 0) then
1344 132 : call ebands%write(dtset%prtebands, dtfil%filnam_ds(4), kptbounds=dtset%kptbounds)
1345 : else
1346 5275 : call ebands%write(dtset%prtebands, dtfil%filnam_ds(4))
1347 : end if
1348 5407 : call timab(1178,2,tsec)
1349 : end if
1350 :
1351 : !Optionally provide Xcrysden output for the Fermi surface (Only master writes)
1352 5444 : if (me == master .and. dtset%prtfsurf == 1) then
1353 1 : call timab(1179,1,tsec)
1354 1 : if (ebands%write_bxsf(crystal,dtfil%fnameabo_app_bxsf) /= 0) then
1355 0 : msg = "Cannot produce BXSF file with Fermi surface, see log file for more info"
1356 0 : ABI_WARNING(msg)
1357 0 : call wrtout(ab_out, msg)
1358 : end if
1359 1 : call timab(1179,2,tsec)
1360 : end if ! prtfsurf==1
1361 :
1362 : !output nesting factor for Fermi surface (requires ph_nqpath)
1363 5444 : if (me == master .and. dtset%prtnest>0 .and. dtset%ph_nqpath > 0) then
1364 0 : call timab(1180,1,tsec)
1365 : ierr = ebands%write_nesting(crystal,dtfil%fnameabo_app_nesting,dtset%prtnest,&
1366 0 : dtset%tsmear,dtset%fermie_nest,dtset%ph_qpath(:,1:dtset%ph_nqpath),msg)
1367 0 : if (ierr /= 0) then
1368 0 : ABI_WARNING(msg)
1369 0 : call wrtout(ab_out, msg)
1370 : end if
1371 0 : call timab(1180,2,tsec)
1372 : end if ! prtnest=1
1373 :
1374 6918 : if (dtset%prtdipole == 1) then
1375 1 : call timab(1181,1,tsec)
1376 : call multipoles_out(rhor,mpi_enreg,natom,nfft,ngfft,dtset%nspden,dtset%ntypat,rprimd,&
1377 1 : dtset%typat,ucvol,ab_out,xred,dtset%ziontypat)
1378 1 : call timab(1181,2,tsec)
1379 : end if
1380 :
1381 : ! BoltzTraP output files in GENEric format
1382 6918 : if (dtset%prtbltztrp == 1 .and. me==master)then
1383 1 : call timab(1182,1,tsec)
1384 1 : call ebands%prtbltztrp(crystal, dtfil%filnam_ds(4))
1385 1 : call timab(1182,2,tsec)
1386 : endif
1387 :
1388 : ! Band structure interpolation from eigenvalues computed on the k-mesh.
1389 6918 : if (nint(dtset%einterp(1)) /= 0 .and. dtset%kptopt > 0) then
1390 1 : call timab(1183,1,tsec)
1391 1 : call ebands%interpolate_kpath(dtset, crystal, [0, 0], dtfil%filnam_ds(4), comm)
1392 1 : call timab(1183,2,tsec)
1393 : end if
1394 :
1395 6918 : call timab(1154,1,tsec)
1396 :
1397 : ! Output of the GSR file (except when we are inside mover)
1398 6918 : if (me == master .and. dtset%prtgsr == 1 .and. dtset%usewvl == 0) then
1399 : !.and. (dtset%ionmov /= 0 .or. dtset%optcell /= 0)) then
1400 4147 : fname = strcat(dtfil%filnam_ds(4), "_GSR.nc")
1401 : ! Write crystal and band structure energies.
1402 4147 : NCF_CHECK(nctk_open_create(ncid, fname, xmpi_comm_self))
1403 4147 : NCF_CHECK(hdr%ncwrite(ncid, fform_den, spinat=dtset%spinat, nc_define=.True.))
1404 4147 : NCF_CHECK(crystal%ncwrite(ncid))
1405 4147 : NCF_CHECK(ebands%ncwrite(ncid))
1406 : ! Add energy, forces, stresses
1407 4147 : NCF_CHECK(results_gs%ncwrite(ncid, dtset%ecut, dtset%pawecutdg))
1408 :
1409 : ! Add info on GBT.
1410 : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: &
1411 : "use_gbt" &
1412 8294 : ], defmode=.True.)
1413 4147 : NCF_CHECK(ncerr)
1414 :
1415 : ncerr = nctk_def_arrays(ncid, [ &
1416 : nctkarr_t("qgbt", "dp", "three") &
1417 8294 : ])
1418 4147 : NCF_CHECK(ncerr)
1419 :
1420 4147 : NCF_CHECK(nctk_set_datamode(ncid))
1421 : ncerr = nctk_write_iscalars(ncid, [character(len=nctk_slen) :: &
1422 : "use_gbt"], &
1423 12441 : [dtset%use_gbt])
1424 4147 : NCF_CHECK(ncerr)
1425 4147 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "qgbt"), dtset%qgbt))
1426 :
1427 4147 : if (allocated(intgden)) then
1428 : ! Write integrated density inside atomic spheres and ratsph(ntypat)=radius of spheres around atoms
1429 : ncerr = nctk_def_arrays(ncid, [ &
1430 : nctkarr_t("intgden", "dp", "number_of_components, number_of_atoms"), &
1431 : nctkarr_t("ratsph", "dp", "number_of_atom_species"), &
1432 : nctkarr_t("rhomag", "dp", "two, number_of_components") &
1433 16588 : ], defmode=.True.)
1434 4147 : NCF_CHECK(ncerr)
1435 4147 : NCF_CHECK(nctk_set_datamode(ncid))
1436 4147 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "intgden"), intgden))
1437 4147 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ratsph"), dtset%ratsph))
1438 4147 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "rhomag"), rhomag))
1439 : end if
1440 :
1441 4147 : if(allocated(efg)) then
1442 : ! write EFG tensors to GSR if available
1443 : ncerr = nctk_def_dims(ncid, [ &
1444 : nctkdim_t("ndir",3),&
1445 : nctkdim_t("natom",dtset%natom),&
1446 36 : nctkdim_t("ntypat",dtset%ntypat)],defmode=.True.)
1447 9 : NCF_CHECK(ncerr)
1448 : ncerr = nctk_def_arrays(ncid, [&
1449 : nctkarr_t("quadmom", "dp", "ntypat"),&
1450 27 : nctkarr_t("efg", "dp", "ndir, ndir, natom")])
1451 9 : NCF_CHECK(ncerr)
1452 9 : NCF_CHECK(nctk_set_datamode(ncid))
1453 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "quadmom"), dtset%quadmom))
1454 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "efg"), efg))
1455 : end if
1456 :
1457 4147 : if(associated(rcpaw)) then
1458 3 : nphicor=0
1459 6 : do itypat=1,dtset%ntypat
1460 6 : nphicor=max(nphicor,rcpaw%atm(itypat)%ln_size)
1461 : enddo
1462 9 : ABI_MALLOC(nphicor_arr,(dtset%ntypat))
1463 12 : ABI_MALLOC(energy_cor,(nphicor,dtset%ntypat))
1464 9 : ABI_MALLOC(occ_cor,(nphicor,dtset%ntypat))
1465 12 : ABI_MALLOC(lcor,(nphicor,dtset%ntypat))
1466 15 : lcor=0
1467 15 : energy_cor=zero
1468 15 : occ_cor=one
1469 6 : do itypat=1,dtset%ntypat
1470 3 : nphicor_arr(itypat)=rcpaw%atm(itypat)%ln_size
1471 15 : do ii=1,rcpaw%atm(itypat)%ln_size
1472 9 : lcor(ii,itypat)=rcpaw%atm(itypat)%indln(1,ii)
1473 9 : energy_cor(ii,itypat)=rcpaw%atm(itypat)%eig(ii,1)
1474 12 : occ_cor(ii,itypat)=rcpaw%atm(itypat)%occ(ii,1)
1475 : enddo
1476 : enddo
1477 : ncerr = nctk_def_dims(ncid, [ &
1478 : nctkdim_t("max_number_of_core_states",nphicor),&
1479 9 : nctkdim_t("number_of_atom_types",dtset%ntypat)],defmode=.True.)
1480 3 : NCF_CHECK(ncerr)
1481 : ncerr = nctk_def_arrays(ncid, [&
1482 : nctkarr_t("eigenvalues_core", "dp", "max_number_of_core_states,number_of_atom_types"),&
1483 : nctkarr_t("occupation_core", "dp", "max_number_of_core_states,number_of_atom_types"),&
1484 : nctkarr_t("number_of_core_states", "int", "number_of_atom_types"),&
1485 15 : nctkarr_t("l_quantum_number_core", "int","max_number_of_core_states,number_of_atom_types")])
1486 3 : NCF_CHECK(ncerr)
1487 3 : NCF_CHECK(nctk_set_datamode(ncid))
1488 3 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "eigenvalues_core"),energy_cor))
1489 3 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "occupation_core"),occ_cor))
1490 3 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "number_of_core_states"),nphicor_arr))
1491 3 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid,"l_quantum_number_core"),lcor))
1492 3 : ABI_FREE(nphicor_arr)
1493 3 : ABI_FREE(energy_cor)
1494 3 : ABI_FREE(occ_cor)
1495 3 : ABI_FREE(lcor)
1496 : endif
1497 :
1498 4147 : NCF_CHECK(nf90_close(ncid))
1499 : end if
1500 :
1501 6918 : call timab(1154,2,tsec)
1502 :
1503 6918 : ABI_SFREE_PTR(elfr)
1504 6918 : ABI_SFREE_PTR(grhor)
1505 6918 : ABI_SFREE_PTR(lrhor)
1506 :
1507 6918 : ABI_SFREE(efg)
1508 6918 : ABI_SFREE(intgden)
1509 :
1510 6918 : call crystal%free()
1511 6918 : call ebands%free()
1512 :
1513 : ! Destroy atom table used for parallelism
1514 6918 : call free_my_atmtab(my_atmtab, my_atmtab_allocated)
1515 :
1516 6918 : call timab(1150,2,tsec) ! outscfcv
1517 :
1518 : DBG_EXIT("COLL")
1519 :
1520 27672 : end subroutine outscfcv
1521 : !!***
1522 :
1523 : end module m_outscfcv
1524 : !!***
|
