Branch data 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 : 20853 : subroutine outscfcv(atindx1,cg,compch_fft,compch_sph,cprj,dimcprj,dmatpawu,dtfil,dtset,&
198 : 13902 : & 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 : 6951 : & nattyp,nfft,ngfft,nhat,nkpt,npwarr,nspden,nsppol,nsym,ntypat,n3xccc,occ,&
201 [ - + ]: 6951 : & paw_dmft,pawang,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,paw_an,paw_ij,&
202 : 6951 : & prtvol,psps,results_gs,rhor,rprimd,&
203 [ - + ]: 13902 : & 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 : 13902 : 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 : 6951 : integer, allocatable :: isort(:)
262 : 6951 : integer, pointer :: my_atmtab(:)
263 : 6951 : real(dp) :: tsec(2),nt_ntone_norm(nspden),rhomag(2,nspden)
264 : 6951 : real(dp),allocatable :: efg(:,:,:),eigen2(:)
265 : 6951 : real(dp),allocatable :: elfr_down(:,:),elfr_up(:,:),intgden(:,:)
266 : 6951 : real(dp),allocatable :: rhor_paw(:,:),rhor_paw_core(:,:),rhor_paw_val(:,:),vpaw(:,:),vwork(:,:)
267 : 6951 : real(dp),allocatable :: rhor_n_one(:,:),rhor_nt_one(:,:),ps_norms(:,:,:)
268 : 6951 : real(dp), allocatable :: doccde(:)
269 : 6951 : real(dp), allocatable :: vh1spl(:)
270 : 6951 : real(dp), allocatable :: vh1_interp(:)
271 : 6951 : real(dp), allocatable :: vh1_integ(:)
272 : 6951 : real(dp), allocatable :: vh1_corrector(:)
273 : 6951 : real(dp), allocatable :: radii(:)
274 : 6951 : real(dp), contiguous, pointer :: rho_ptr(:,:)
275 [ + + - + : 69510 : type(pawrhoij_type) :: pawrhoij_dum(1)
- + - + -
+ - + - +
+ + + - +
- + - + -
+ - + - ]
276 : : !type(pawrhoij_type) :: pawrhoij_dum(0)
277 : 6951 : 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 : 6951 : type(oper_type) :: dft_occup
284 [ + + + + : 361452 : type(crystal_t) :: crystal
+ + + + +
+ + + + +
+ + + + ]
285 : 6951 : type(ebands_t) :: ebands
286 : 6951 : type(epjdos_t) :: dos
287 : 6951 : type(plowannier_type) :: wan
288 : 6951 : type(self_type) :: selfr
289 : 6951 : type(self_type), target :: self
290 : 6951 : type(green_type) :: greenr
291 : 6951 : type(matlu_type), allocatable :: opt_selflimit(:)
292 : 6951 : integer,allocatable :: nphicor_arr(:),lcor(:,:)
293 : 6951 : real(dp),allocatable :: energy_cor(:,:),occ_cor(:,:)
294 : :
295 : : ! *************************************************************************
296 : :
297 : : DBG_ENTER("COLL")
298 : :
299 [ + + ]: 20853 : units = [std_out, ab_out]
300 : :
301 : 6951 : call timab(1150,1,tsec) ! outscfcv
302 : 6951 : call timab(1151,1,tsec) ! outscfcv(preparation)
303 : :
304 [ + + - + : 6951 : 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 : 6951 : n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3); nfft_tot = product(ngfft(1:3))
316 : 6951 : comm_fft = mpi_enreg%comm_fft
317 : 6951 : me_fft = xmpi_comm_rank(comm_fft)
318 : 6951 : paral_fft = (mpi_enreg%paral_kgb==1)
319 : :
320 : 6951 : comm = mpi_enreg%comm_cell
321 : 6951 : me = xmpi_comm_rank(comm)
322 : :
323 : 6951 : paral_atom=(my_natom/=natom)
324 : 6951 : my_comm_atom = mpi_enreg%comm_atom
325 : 6951 : nullify(my_atmtab)
326 [ + + ]: 6951 : 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 [ - + + - : 19179 : ABI_MALLOC(my_atmtab, (natom))
+ - + - -
+ - + ]
330 [ + - + + : 44141 : my_atmtab = (/ (iatom, iatom=1, natom) /)
+ + ]
331 : 6393 : 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 : 6951 : remove_inv=.false.
340 : : ! What about Wannier90 DMFT?
341 [ + + + + ]: 6951 : if (dtset%nspden==4 .and. dtset%usedmft==1) remove_inv=.true. ! MG: why this?
342 : :
343 [ + + ]: 6951 : 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 [ + + + + ]: 13872 : dtset%symrel,dtset%tnons,dtset%symafm)
348 : :
349 : : ! Electron band energies.
350 : 6951 : bantot= dtset%mband*dtset%nkpt*dtset%nsppol
351 [ - + + - : 514662 : 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 : 6951 : hdr%kptrlatt, hdr%nshiftk, hdr%shiftk)
357 : :
358 [ - + ]: 6951 : ABI_FREE(doccde)
359 : :
360 : 6951 : ebands%fermie = results_gs%energies%e_fermie
361 : 6951 : e_fermie = results_gs%energies%e_fermie
362 : 6951 : ebands%fermih = results_gs%energies%e_fermih
363 : 6951 : e_fermih = results_gs%energies%e_fermih
364 : 6951 : ebands%entropy = results_gs%energies%entropy_ks
365 : :
366 : : ! YAML output
367 [ + + ]: 6951 : if (me == master) then
368 : : call results_gs%yaml_write(ab_out, cryst=crystal, info="Summary of ground state results",&
369 : 5474 : occopt=dtset%occopt, with_conv=(dtset%nstep > 0))
370 : : end if
371 : :
372 : 6951 : call timab(1151,2,tsec)
373 : :
374 : : !wannier interface
375 : 6951 : call timab(1152,1,tsec)
376 : :
377 [ + + + - ]: 6951 : 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 [ + + ]: 6938 : 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 : 6951 : call timab(1152,2,tsec)
417 : 6951 : call timab(1153,1,tsec)
418 : :
419 : 6951 : occopt=dtset%occopt
420 : :
421 : 6951 : prtnabla=dtset%prtnabla
422 : 6951 : pawprtden=dtset%prtden-1
423 : :
424 : 6951 : comm=mpi_enreg%comm_cell; me=xmpi_comm_rank(comm)
425 : 6951 : comm_fft=mpi_enreg%comm_fft
426 : 6951 : 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 [ + + + - ]: 6951 : if (iwrite_fftdatar(mpi_enreg) .and. dtset%usewvl==0) then
433 : :
434 : : ! output the density.
435 [ + + ]: 5593 : if (dtset%prtden/=0) then
436 [ + + ]: 3588 : if (dtset%positron/=1) rho_ptr => rhor
437 [ + + ]: 3588 : if (dtset%positron==1) rho_ptr => electronpositron%rhor_ep
438 : : call fftdatar_write("density",dtfil%fnameabo_app_den,dtset%iomode,hdr,&
439 : 3588 : crystal,ngfft,cplex1,nfft,nspden,rho_ptr,mpi_enreg,ebands=ebands)
440 : :
441 [ + + ]: 3588 : 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 [ - + - - ]: 1358 : 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 [ + + + - ]: 6951 : 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 : 6951 : call timab(1153,2,tsec)
594 : :
595 : 6951 : call timab(1155,1,tsec)
596 : :
597 : : ! Output of VCLMB file
598 : : ! The PAW correction has to be computed here (all processors contribute)
599 [ + + + + ]: 6951 : 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 : 6951 : call timab(1155,2,tsec)
696 : :
697 [ + + ]: 6951 : if (iwrite_fftdatar(mpi_enreg)) then
698 : :
699 : 5593 : call timab(1156,1,tsec)
700 : :
701 : : ! output the electron localization function ELF
702 [ + + ]: 5593 : 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 : 5593 : call timab(1156,2,tsec)
731 : 5593 : call timab(1157,1,tsec)
732 : :
733 : : ! We output the gradient of density
734 [ + + ]: 5593 : 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 : 5593 : call timab(1157,2,tsec)
747 : 5593 : call timab(1158,1,tsec)
748 : :
749 : : ! We output the total kinetic energy density KDEN
750 [ + + ]: 5593 : 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 : 5593 : call timab(1158,2,tsec)
756 : 5593 : call timab(1159,1,tsec)
757 : :
758 : :
759 : : ! We output the Laplacian of density
760 [ + + ]: 5593 : 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 : 5593 : call timab(1159,2,tsec)
766 : 5593 : call timab(1160,1,tsec)
767 : :
768 : : ! POT
769 [ + + ]: 5593 : 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 [ - + - - ]: 5593 : 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 : 5593 : call timab(1160,2,tsec)
782 : 5593 : call timab(1161,1,tsec)
783 : :
784 [ + + ]: 5593 : 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 [ + + ]: 5593 : 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 : 5593 : call timab(1161,2,tsec)
795 : 5593 : call timab(1162,1,tsec)
796 : :
797 : : ! STM
798 [ + + ]: 5593 : 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 : 5593 : call timab(1162,2,tsec)
803 : 5593 : call timab(1163,1,tsec)
804 : :
805 [ + + ]: 5593 : 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 : 5593 : call timab(1163,2,tsec)
810 : 5593 : call timab(1164,1,tsec)
811 : :
812 : : ! VHA
813 [ + + ]: 5593 : 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 [ + + ]: 5593 : 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 [ + + ]: 5593 : 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 [ + + ]: 5593 : 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 [ + + ]: 5593 : 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 : 5593 : call timab(1164,2,tsec)
887 : :
888 : : end if ! if iwrite_fftdatar
889 : :
890 : 6951 : 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 [ + + + + ]: 6951 : 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 : 6951 : call timab(1165,2,tsec)
936 : 6951 : call timab(1166,1,tsec)
937 : :
938 : : !Output of integrated density inside atomic spheres
939 [ + + + - : 27804 : if ((dtset%prtdensph==1.and.dtset%usewvl==0) .or. sum(abs(dtset%hspinfield)) > tol10) then
- + - - ]
940 [ - + + - : 27804 : ABI_MALLOC(intgden, (nspden, natom))
- + + - +
- + - - +
- + - + ]
941 : :
942 : : call calcdenmagsph(mpi_enreg,natom,nfft,ngfft,nspden,&
943 : 6951 : 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 [ + - + - ]: 6951 : if (dtset%prtdensph==1.and.dtset%usewvl==0) then
950 [ + + + + : 16482 : if(all(dtset%constraint_kind(:)==0))then
+ + ]
951 : 6897 : 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 [ + + + + : 16482 : 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 [ + + ]: 6951 : 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 [ + + + + ]: 27804 : 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 : 6951 : call timab(1166,2,tsec)
1009 : 6951 : call timab(1167,1,tsec)
1010 : :
1011 [ + + ]: 6951 : 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 : 6951 : call timab(1167,2,tsec)
1018 : 6951 : call timab(1168,1,tsec)
1019 : :
1020 : : !If PAW, provide additional outputs
1021 [ + + ]: 6951 : if (psps%usepaw==1) then
1022 : : ! Output of compensation charge
1023 [ + + + + ]: 1523 : if (dtset%nstep>0.or.dtfil%ireadwf/=0) then
1024 : 1522 : write(msg, '(4a)' )ch10,' PAW TEST:',ch10,&
1025 : 3044 : & ' ==== Compensation charge inside spheres ============'
1026 [ + - + + ]: 1522 : if (compch_sph>-1.d4.and.compch_fft>-1.d4) &
1027 : 1502 : & write(msg, '(3a)' ) trim(msg),ch10,' The following values must be close to each other ...'
1028 [ + - ]: 3044 : if (compch_sph>-1.d4) write(msg, '(3a,f22.15)' ) trim(msg),ch10,&
1029 : 3044 : & ' Compensation charge over spherical meshes = ',compch_sph
1030 [ + + ]: 1522 : if (compch_fft>-1.d4) then
1031 [ + + ]: 1502 : if (pawfgr%usefinegrid==1) then
1032 : 1433 : write(msg, '(3a,f22.15)' ) trim(msg),ch10,&
1033 : 2866 : & ' 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 : 1522 : 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 : 1523 : & electronpositron=electronpositron)
1045 : : end if
1046 : :
1047 : 6951 : call timab(1168,2,tsec)
1048 : 6951 : call timab(1169,1,tsec)
1049 : :
1050 : :
1051 : : !PAW + output for optical conductivity _OPT and _OPT2
1052 [ + + + + ]: 6951 : 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 [ - + ]: 6951 : 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 : 6951 : call timab(1169,2,tsec)
1113 : 6951 : call timab(1170,1,tsec)
1114 : :
1115 : : !Optionally provide output for AE wavefunctions (only for PAW)
1116 [ + + + + ]: 6951 : 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 : 6951 : call timab(1170,2,tsec)
1164 : 6951 : call timab(1171,1,tsec)
1165 : :
1166 [ - + ]: 6951 : 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 : 6951 : call timab(1171,2,tsec)
1189 : 6951 : call timab(1172,1,tsec)
1190 : :
1191 : : !Optionally provide output for the GW part of ABINIT
1192 [ + + ]: 6951 : 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 : 6951 : call timab(1172,2,tsec) ! outscfcv(gw)
1286 : :
1287 [ + + ]: 6951 : 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 [ + + ]: 6951 : 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 [ + + ]: 6951 : 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 [ + + ]: 6951 : 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 [ + + + + ]: 6951 : if (me == master .and. dtset%tfkinfunc==0) then
1342 : 5437 : call timab(1178,1,tsec)
1343 [ + + ]: 5437 : 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 : 5305 : call ebands%write(dtset%prtebands, dtfil%filnam_ds(4))
1347 : : end if
1348 : 5437 : call timab(1178,2,tsec)
1349 : : end if
1350 : :
1351 : : !Optionally provide Xcrysden output for the Fermi surface (Only master writes)
1352 [ + + ]: 5474 : 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 [ + + - + ]: 5474 : 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 [ + + ]: 6951 : 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 [ + + + - ]: 6951 : 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 [ + + + - ]: 6951 : 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 : 6951 : call timab(1154,1,tsec)
1396 : :
1397 : : ! Output of the GSR file (except when we are inside mover)
1398 [ + + + + : 6951 : if (me == master .and. dtset%prtgsr == 1 .and. dtset%usewvl == 0) then
+ - ]
1399 : : !.and. (dtset%ionmov /= 0 .or. dtset%optcell /= 0)) then
1400 [ + - ]: 4177 : fname = strcat(dtfil%filnam_ds(4), "_GSR.nc")
1401 : : ! Write crystal and band structure energies.
1402 [ - + ]: 4177 : NCF_CHECK(nctk_open_create(ncid, fname, xmpi_comm_self))
1403 [ - + ]: 4177 : NCF_CHECK(hdr%ncwrite(ncid, fform_den, spinat=dtset%spinat, nc_define=.True.))
1404 [ - + ]: 4177 : NCF_CHECK(crystal%ncwrite(ncid))
1405 [ - + ]: 4177 : NCF_CHECK(ebands%ncwrite(ncid))
1406 : : ! Add energy, forces, stresses
1407 [ - + ]: 4177 : 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 [ + + ]: 8354 : ], defmode=.True.)
1413 [ - + ]: 4177 : NCF_CHECK(ncerr)
1414 : :
1415 : : ncerr = nctk_def_arrays(ncid, [ &
1416 : : nctkarr_t("qgbt", "dp", "three") &
1417 [ + + ]: 8354 : ])
1418 [ - + ]: 4177 : NCF_CHECK(ncerr)
1419 : :
1420 [ - + ]: 4177 : NCF_CHECK(nctk_set_datamode(ncid))
1421 : : ncerr = nctk_write_iscalars(ncid, [character(len=nctk_slen) :: &
1422 : : "use_gbt"], &
1423 [ + + + + ]: 12531 : [dtset%use_gbt])
1424 [ - + ]: 4177 : NCF_CHECK(ncerr)
1425 [ - + ]: 4177 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "qgbt"), dtset%qgbt))
1426 : :
1427 [ + - ]: 4177 : 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 [ + + ]: 16708 : ], defmode=.True.)
1434 [ - + ]: 4177 : NCF_CHECK(ncerr)
1435 [ - + ]: 4177 : NCF_CHECK(nctk_set_datamode(ncid))
1436 [ - + ]: 4177 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "intgden"), intgden))
1437 [ - + ]: 4177 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ratsph"), dtset%ratsph))
1438 [ - + ]: 4177 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "rhomag"), rhomag))
1439 : : end if
1440 : :
1441 [ + + ]: 4177 : 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 [ + + ]: 4177 : 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 [ - + ]: 4177 : NCF_CHECK(nf90_close(ncid))
1499 : : end if
1500 : :
1501 : 6951 : call timab(1154,2,tsec)
1502 : :
1503 [ + + ]: 6951 : ABI_SFREE_PTR(elfr)
1504 [ + + ]: 6951 : ABI_SFREE_PTR(grhor)
1505 [ + + ]: 6951 : ABI_SFREE_PTR(lrhor)
1506 : :
1507 [ + + ]: 6951 : ABI_SFREE(efg)
1508 [ + - ]: 6951 : ABI_SFREE(intgden)
1509 : :
1510 : 6951 : call crystal%free()
1511 : 6951 : call ebands%free()
1512 : :
1513 : : ! Destroy atom table used for parallelism
1514 : 6951 : call free_my_atmtab(my_atmtab, my_atmtab_allocated)
1515 : :
1516 : 6951 : call timab(1150,2,tsec) ! outscfcv
1517 : :
1518 : : DBG_EXIT("COLL")
1519 : :
1520 [ - + - + : 27804 : end subroutine outscfcv
- + - - -
- - + - +
- - - - -
+ - - - -
- + - - -
- - - - -
- - - + -
+ - + - +
- + - + -
- - - - +
- + - - -
- - - - -
- - - + -
- - - - -
- - - - -
+ - - - -
- + - + -
- - - - -
- - - - -
+ - - - -
- - - - -
- - + - +
- + - + -
+ - - - -
- + - + -
- - - - +
- + - + -
+ - + - +
- + - + -
+ - - - -
- - - - -
- - + - -
- - - - -
- - - - +
- - - - -
- - - - -
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - -
- - - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + -
+ ]
1521 : : !!***
1522 : :
1523 : : end module m_outscfcv
1524 : : !!***
|
