TABLE OF CONTENTS
ABINIT/m_outscfcv [ Modules ]
NAME
m_outscfcv
FUNCTION
COPYRIGHT
Copyright (C) 2005-2018 ABINIT group (XG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
PARENTS
CHILDREN
SOURCE
20 #if defined HAVE_CONFIG_H 21 #include "config.h" 22 #endif 23 24 #include "abi_common.h" 25 26 module m_outscfcv 27 28 implicit none 29 30 private
ABINIT/outscfcv [ Functions ]
NAME
outscfcv
FUNCTION
Output routine for the scfcv.F90 routine
INPUTS
atindx(natom)=index table for atoms (see gstate.f) atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f) cg(2,mcg)=planewave coefficients of wavefunctions (see also side effects) compch_fft=compensation charge, from FFT grid compch_sph=compensation charge, from sphere cprj(natom,mcprj*usecprj)=<p_lmn|Cnk> coefficients for each WF |Cnk> and each |p_lmn> non-local projector. See also side effects dimcprj(natom*usecprj)=array of dimensions of array cprj (not ordered) dmatpawu= fixed occupation matrix of correlated orbitals (DFT+U or DMFT only) dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables for this dataset ecut=cut-off energy for plane wave basis sphere (Ha) eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree) elfr(nfft,nspden(+1))=electron localization function, real space. (+1) if spin-polarized in order to get total, spin up and spin down elf etotal=total energy gmet(3,3)=reciprocal space metric tensor in bohr**-2. gprimd(3,3)=dimensional reciprocal space primitive translations grhor(nfft,nspden,3)= gradient of electron density in electrons/bohr**4, real space hdr <type(hdr_type)>=the header of wf, den and pot files kg(3,mpw*mkmem)=reduced planewave coordinates. lrhor(nfft,nspden)= Laplacian of electron density in electrons/bohr**5, real space mband=maximum number of bands mcg=size of wave-functions array (cg) =mpw*my_nspinor*mband*mkmem*nsppol mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol mgfftc=maximum size of 1D FFTs for the PAW coarse grid mkmem =number of k points treated by this node. mpi_enreg=information about MPI parallelization mpsang= 1+maximum angular momentum for nonlocal pseudopotentials mpw=maximum dimensioned size of npw. my_natom=number of atoms treated by current processor natom=number of atoms in cell. nattyp(ntypat)= # atoms of each type. nfft=(effective) number of FFT grid points (for this processor) (see NOTES at beginning of scfcv) ngfft(18)=contain all needed information about 3D FFT (see NOTES at beginning of scfcv) nhat(nfft,nspden*usepaw)= compensation charge density (PAW) nkpt=number of k points. npwarr(nkpt)=number of planewaves in basis at this k point nspden=number of spin-density components nsppol=1 for unpolarized, 2 for spin-polarized nsym=number of symmetries in space group ntypat=number of types of atoms in unit cell. n3xccc=dimension of the xccc3d array (0 or nfft). occ(mband*nkpt*nsppol)=occupation number for each band (usually 2) for each k. paw_an(my_natom) <type(paw_an_type)>=paw arrays given on angular mesh pawang <type(pawang_type)>=paw angular mesh and related data pawfgr(natom) <type(pawfgr_type)>=fine grid parameters and related data pawfgrtab(my_natom) <type(pawfgrtab_type)> tables on PAW fine grid pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data paw_ij(my_natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels note:structure factors are given on the coarse grid for PAW prtvol=control print volume and debugging output psps <type(pseudopotential_type)>=variables related to pseudopotentials results_gs <type(results_gs_type)>=results (energy and its components, forces and its components, the stress tensor) of a ground-state computation rhor(nfft,nspden)=total electron density in electrons/bohr**3, real space. rprimd(3,3)=dimensional primitive translations for real space (bohr) taur(nfft,nspden)=total kinetic energy density in bohr**(-5), real space. ucvol=unit cell volume (bohr**3) usecprj=1 if cprj datastructure has been allocated vhartr(nfft)=Hartree potential vxc(nfft,nspden)=xc potential vtrial(nfft,nspden)=the trial potential = vxc + vpsp + vhartr, roughly speaking xccc3d(n3xccc)=3D core electron density for XC core correction (bohr^-3) xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
(only writing, printing)
SIDE EFFECTS
electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation If prtwant==3 the following quantitities are updated using the unitary transformation defining the QP amplitudes in terms of the KS basis set: cg(2,mcg)=planewave coefficients of wavefunctions. cprj(natom,mcprj*usecpyj)=<p_lmn|Cnk> coefficients for each WF |Cnk> and each |p_lmn> non-local projector
NOTES
The string passed to fftdatar_write (first argument) gives the name used to store the data in the netcdf file The function varname_from_fname defined in the module m_hdr.F90 gives the mapping between the Abinit file extension and the netcdf name e.g. foo_VHXC.nc --> vxc This function is used in cut3d so that we can immediately select the data to analyze without having to prompt the user. Remember to update varname_from_fname if you add a new file or if you change the name of the variable.
PARENTS
scfcv
CHILDREN
bonds_lgth_angles,bound_deriv,calc_efg,calc_fc,calcdensph compute_coeff_plowannier,crystal_free,crystal_init,datafordmft,denfgr destroy_dmft,destroy_oper,destroy_plowannier,dos_calcnwrite,ebands_free ebands_init,ebands_interpolate_kpath,ebands_prtbltztrp,ebands_write epjdos_free,fatbands_ncwrite,fftdatar_write,free_my_atmtab get_my_atmtab,init_dmft,init_oper,init_plowannier,ioarr,mag_constr_e mlwfovlp,mlwfovlp_qp,multipoles_out,optics_paw,optics_paw_core optics_vloc,out1dm,outkss,outwant,partial_dos_fractions partial_dos_fractions_paw,pawmkaewf,pawprt,pawrhoij_copy pawrhoij_nullify,posdoppler,poslifetime,print_dmft,prt_cif,prtfatbands read_atomden,simpson_int,sort_dp,spline,splint,timab,wrtout,xmpi_sum xmpi_sum_master
SOURCE
152 subroutine outscfcv(atindx1,cg,compch_fft,compch_sph,cprj,dimcprj,dmatpawu,dtfil,dtset,& 153 & ecut,eigen,electronpositron,elfr,etotal,gmet,gprimd,grhor,hdr,kg,& 154 & lrhor,mband,mcg,mcprj,mgfftc,mkmem,mpi_enreg,mpsang,mpw,my_natom,natom,& 155 & nattyp,nfft,ngfft,nhat,nkpt,npwarr,nspden,nsppol,nsym,ntypat,n3xccc,occ,& 156 & paw_dmft,pawang,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,paw_an,paw_ij,& 157 & prtvol,psps,results_gs,rhor,rprimd,& 158 & taur,ucvol,usecprj,vhartr,vpsp,vtrial,vxc,wvl_den,xccc3d,xred) 159 160 use defs_basis 161 use defs_datatypes 162 use defs_wvltypes 163 use defs_abitypes 164 use m_abicore 165 use m_sort 166 use m_efield 167 use m_errors 168 use m_xmpi 169 use m_mpinfo 170 #ifdef HAVE_NETCDF 171 use netcdf 172 #endif 173 use m_abi_etsf 174 use m_nctk 175 use m_hdr 176 use m_plowannier 177 use m_splines 178 use m_ebands 179 180 use m_time, only : timab 181 use m_io_tools, only : open_file 182 use m_fstrings, only : strcat, endswith 183 use m_geometry, only : bonds_lgth_angles 184 use m_electronpositron, only : electronpositron_type,electronpositron_calctype 185 use m_oper, only : oper_type,init_oper,destroy_oper 186 use m_crystal, only : crystal_init, crystal_free, crystal_t, prt_cif 187 use m_crystal_io, only : crystal_ncwrite 188 use m_results_gs, only : results_gs_type, results_gs_ncwrite 189 use m_ioarr, only : ioarr, fftdatar_write 190 use m_nucprop, only : calc_efg,calc_fc 191 use m_outwant, only : outwant 192 use m_pawang, only : pawang_type 193 use m_pawrad, only : pawrad_type, simp_gen, bound_deriv 194 use m_pawtab, only : pawtab_type 195 use m_paw_an, only : paw_an_type 196 use m_paw_ij, only : paw_ij_type 197 use m_paw_mkrho, only : denfgr 198 use m_pawfgrtab, only : pawfgrtab_type 199 use m_pawrhoij, only : pawrhoij_type, pawrhoij_nullify, pawrhoij_copy, pawrhoij_free 200 use m_pawcprj, only : pawcprj_type 201 use m_pawfgr, only : pawfgr_type 202 use m_paw_dmft, only : paw_dmft_type,init_dmft,destroy_dmft,print_dmft 203 use m_paw_optics, only : optics_paw,optics_paw_core 204 use m_paw_tools, only : pawprt 205 use m_numeric_tools, only : simpson_int 206 use m_epjdos, only : dos_calcnwrite, partial_dos_fractions, partial_dos_fractions_paw, & 207 epjdos_t, epjdos_new, epjdos_free, prtfatbands, fatbands_ncwrite 208 use m_paral_atom, only : get_my_atmtab, free_my_atmtab 209 use m_io_kss, only : outkss 210 use m_multipoles, only : multipoles_out, out1dm 211 use m_mlwfovlp_qp, only : mlwfovlp_qp 212 use m_paw_mkaewf, only : pawmkaewf 213 use m_dens, only : mag_constr_e, calcdensph 214 use m_mlwfovlp, only : mlwfovlp 215 use m_datafordmft, only : datafordmft 216 use m_mkrho, only : read_atomden 217 use m_positron, only : poslifetime, posdoppler 218 use m_optics_vloc, only : optics_vloc 219 220 !This section has been created automatically by the script Abilint (TD). 221 !Do not modify the following lines by hand. 222 #undef ABI_FUNC 223 #define ABI_FUNC 'outscfcv' 224 !End of the abilint section 225 226 implicit none 227 228 !Arguments ------------------------------------ 229 !scalars 230 integer,intent(in) :: mband,mcg,mcprj,mgfftc,mkmem,mpsang,mpw,n3xccc,my_natom,natom,nfft 231 integer,intent(in) :: nkpt,nspden,nsppol,nsym,ntypat,prtvol,usecprj 232 real(dp),intent(in) :: compch_fft,compch_sph,ecut,ucvol 233 real(dp),intent(inout) :: etotal 234 type(electronpositron_type),pointer :: electronpositron 235 type(MPI_type),intent(inout) :: mpi_enreg 236 type(datafiles_type),intent(in) :: dtfil 237 type(dataset_type),intent(in) :: dtset 238 type(hdr_type),intent(inout) :: hdr 239 type(paw_dmft_type), intent(inout) :: paw_dmft 240 type(pawang_type),intent(in) :: pawang 241 type(pawfgr_type),intent(in) :: pawfgr 242 type(pseudopotential_type),intent(in) :: psps 243 type(results_gs_type),intent(in) :: results_gs 244 type(wvl_denspot_type), intent(in) :: wvl_den 245 !arrays 246 integer,intent(in) :: atindx1(natom),dimcprj(natom*usecprj) 247 integer,intent(in) :: kg(3,mpw*mkmem),nattyp(ntypat),ngfft(18),npwarr(nkpt) 248 real(dp),intent(in) :: dmatpawu(:,:,:,:),eigen(mband*nkpt*nsppol) 249 real(dp),intent(in) :: gmet(3,3),gprimd(3,3) 250 real(dp),intent(in) :: occ(mband*nkpt*nsppol) 251 real(dp),intent(in) :: rprimd(3,3),vhartr(nfft),xccc3d(n3xccc) 252 real(dp),intent(in) :: vpsp(nfft) 253 real(dp),intent(inout) :: cg(2,mcg) 254 real(dp),intent(inout) :: nhat(nfft,nspden*psps%usepaw) 255 real(dp),intent(inout),target :: rhor(nfft,nspden),vtrial(nfft,nspden) 256 real(dp),intent(inout) :: vxc(nfft,nspden),xred(3,natom) 257 real(dp),pointer :: elfr(:,:),grhor(:,:,:),lrhor(:,:),taur(:,:) 258 type(pawcprj_type),intent(inout) :: cprj(natom,mcprj*usecprj) 259 type(paw_an_type),intent(inout) :: paw_an(my_natom*psps%usepaw) 260 type(pawfgrtab_type),intent(in) :: pawfgrtab(my_natom*psps%usepaw) 261 type(paw_ij_type),intent(inout) :: paw_ij(my_natom*psps%usepaw) 262 type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw) 263 type(pawrhoij_type),target,intent(inout) :: pawrhoij(my_natom*psps%usepaw) 264 type(pawtab_type),intent(inout) :: pawtab(ntypat*psps%usepaw) 265 266 !Local variables------------------------------- 267 !scalars 268 integer,parameter :: master=0,cplex1=1,fform_den=52,rdwr2=2,rdwrpaw0=0 269 integer :: bantot,fform,collect,timrev 270 integer :: accessfil,coordn 271 integer :: ii,ierr,ifft,ikpt,ispden,isppol,itypat 272 integer :: me_fft,n1,n2,n3 273 integer :: ifgd, iatom, iatom_tot,nradint 274 integer :: me,my_natom_tmp 275 integer :: occopt 276 integer :: prtnabla 277 integer :: pawprtden 278 integer :: iband,nocc,spacecomm,comm_fft,tmp_unt,nfft_tot 279 integer :: my_comm_atom 280 #ifdef HAVE_NETCDF 281 integer :: ncid 282 #endif 283 real(dp) :: norm,occ_norm,unocc_norm 284 real(dp) :: rate_dum,rate_dum2 285 real(dp) :: yp1, ypn, dr 286 character(len=500) :: message 287 character(len=fnlen) :: fname 288 !arrays 289 integer, allocatable :: isort(:) 290 integer, pointer :: my_atmtab(:) 291 real(dp) :: tsec(2),nt_ntone_norm(nspden) 292 real(dp),allocatable :: eigen2(:) 293 real(dp),allocatable :: elfr_down(:,:),elfr_up(:,:) 294 real(dp),allocatable :: rhor_paw(:,:),rhor_paw_core(:,:),rhor_paw_val(:,:),vpaw(:,:),vwork(:,:) 295 real(dp),allocatable :: rhor_n_one(:,:),rhor_nt_one(:,:),ps_norms(:,:,:) 296 real(dp), allocatable :: doccde(:) 297 real(dp), allocatable :: vh1spl(:) 298 real(dp), allocatable :: vh1_interp(:) 299 real(dp), allocatable :: vh1_integ(:) 300 real(dp), allocatable :: vh1_corrector(:) 301 real(dp), allocatable :: radii(:) 302 real(dp), ABI_CONTIGUOUS pointer :: rho_ptr(:,:) 303 type(pawrhoij_type) :: pawrhoij_dum(0) 304 type(pawrhoij_type),pointer :: pawrhoij_all(:) 305 logical :: remove_inv 306 logical :: paral_atom, paral_fft, my_atmtab_allocated 307 real(dp) :: e_fermie 308 type(oper_type) :: lda_occup 309 type(crystal_t) :: crystal 310 type(ebands_t) :: ebands 311 type(epjdos_t) :: dos 312 type(plowannier_type) :: wan 313 314 ! ************************************************************************* 315 316 DBG_ENTER("COLL") 317 call timab(950,1,tsec) ! outscfcv 318 319 if ((usecprj==0.or.mcprj==0).and.psps%usepaw==1.and. & 320 & (dtset%prtwant==2.or.dtset%prtwant==3.or.dtset%prtnabla>0.or.dtset%prtdos==3 & 321 & .or.dtset%kssform==3.or.dtset%pawfatbnd>0.or.dtset%pawprtwf>0)) then 322 write (message,'(5a)')& 323 & 'cprj datastructure must be allocated',ch10,& 324 & 'with options prtwant=2,3, prtnabla>0, prtdos>3, kssform==3, pawfatbnd>0, pawprtwf>0',ch10,& 325 & 'Action: change pawusecp input keyword.' 326 MSG_ERROR(message) 327 end if 328 329 !Initialize two objects to facilitate the propagation of info. 330 !These objects should used more frequently, actually they should 331 !become basic objects used in abinit. 332 333 !Crystalline structure. 334 remove_inv=.false. 335 if(dtset%nspden==4 .and. dtset%usedmft==1) remove_inv=.true. ! MG: why this? 336 337 timrev = 2; if (any(dtset%kptopt == [3, 4])) timrev= 1 338 call crystal_init(dtset%amu_orig(:,1),crystal,dtset%spgroup,natom,dtset%npsp,ntypat, & 339 & dtset%nsym,rprimd,dtset%typat,xred,dtset%ziontypat,dtset%znucl,timrev,& 340 & dtset%nspden==2.and.dtset%nsppol==1,remove_inv,hdr%title,& 341 & dtset%symrel,dtset%tnons,dtset%symafm) 342 343 !Electron band energies. 344 bantot= dtset%mband*dtset%nkpt*dtset%nsppol 345 ABI_MALLOC(doccde,(bantot)) 346 doccde=zero 347 348 call ebands_init(bantot,ebands,dtset%nelect,doccde,eigen,hdr%istwfk,hdr%kptns,hdr%nband,& 349 & hdr%nkpt,hdr%npwarr,hdr%nsppol,hdr%nspinor,hdr%tphysel,hdr%tsmear,hdr%occopt,hdr%occ,hdr%wtk,& 350 & hdr%charge, hdr%kptopt, hdr%kptrlatt_orig, hdr%nshiftk_orig, hdr%shiftk_orig, & 351 & hdr%kptrlatt, hdr%nshiftk, hdr%shiftk) 352 353 ABI_FREE(doccde) 354 355 ebands%fermie = results_gs%energies%e_fermie 356 e_fermie = results_gs%energies%e_fermie 357 ebands%entropy = results_gs%energies%entropy 358 !write(std_out,*)"ebands%efermi in outscfcv",ebands%fermie 359 !write(std_out,*)"results_gs%energies%e_fermie in outscfcv",e_fermie 360 !write(std_out,*)"results_gs%fermie in outscfcv",results_gs%fermie 361 !write(std_out,*)"hdr%efermi in outscfcv",hdr%fermie 362 363 ! Parameters for MPI-FFT 364 n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3); nfft_tot = product(ngfft(1:3)) 365 comm_fft = mpi_enreg%comm_fft 366 me_fft = xmpi_comm_rank(comm_fft) 367 paral_fft = (mpi_enreg%paral_kgb==1) 368 369 spacecomm = mpi_enreg%comm_cell 370 me = xmpi_comm_rank(spacecomm) 371 372 paral_atom=(my_natom/=natom) 373 my_comm_atom = mpi_enreg%comm_atom 374 nullify(my_atmtab) 375 if (paral_atom) then 376 call get_my_atmtab(mpi_enreg%comm_atom, my_atmtab, my_atmtab_allocated, paral_atom,natom,my_natom_ref=my_natom) 377 else 378 ABI_ALLOCATE(my_atmtab, (natom)) 379 my_atmtab = (/ (iatom, iatom=1, natom) /) 380 my_atmtab_allocated = .true. 381 end if 382 383 !wannier interface 384 call timab(951,1,tsec) 385 if (dtset%prtwant==2) then 386 387 call mlwfovlp(crystal, ebands, hdr, atindx1,cg,cprj,dtset,dtfil,eigen,gprimd,kg,& 388 & mband,mcg,mcprj,mgfftc,mkmem,mpi_enreg,mpw,natom,& 389 & nattyp,nfft,ngfft,nkpt,npwarr,nsppol,ntypat,occ,& 390 & pawang,pawrad,pawtab,prtvol,psps,rprimd,ucvol,xred) 391 392 else if (dtset%prtwant==3) then 393 394 ! Convert cg and eigen to GW quasiparticle wave functions and eigenvalues in mlwfovlp_qp 395 ABI_ALLOCATE(eigen2,(mband*nkpt*nsppol)) 396 eigen2=eigen 397 398 call mlwfovlp_qp(cg,cprj,dtset,dtfil,eigen2,mband,mcg,mcprj,mkmem,mpw,natom,& 399 & nkpt,npwarr,nspden,nsppol,ntypat,Hdr,pawtab,rprimd,MPI_enreg) 400 401 ! Call Wannier90 402 call mlwfovlp(crystal, ebands, hdr, atindx1,cg,cprj,dtset,dtfil,eigen2,gprimd,kg,& 403 & mband,mcg,mcprj,mgfftc,mkmem,mpi_enreg,mpw,natom,& 404 & nattyp,nfft,ngfft,nkpt,npwarr,nsppol,ntypat,occ,& 405 & pawang,pawrad,pawtab,prtvol,psps,rprimd,ucvol,xred) 406 407 ! this is the old implementation, risky due to unpredictable size effects 408 ! now eigen is not overwritten, one should use other ways to print the GW corrections 409 ! eigen=eigen2 410 ABI_DEALLOCATE(eigen2) 411 end if !prtwant 412 call timab(951,2,tsec) 413 414 occopt=dtset%occopt 415 416 prtnabla=dtset%prtnabla 417 pawprtden=dtset%prtden-1 418 419 call timab(952,1,tsec) 420 421 spacecomm=mpi_enreg%comm_cell; me=xmpi_comm_rank(spacecomm) 422 comm_fft=mpi_enreg%comm_fft 423 paral_atom=(my_natom/=natom) 424 425 !Warnings : 426 !- core charge is excluded from the charge density; 427 !- the potential is the INPUT vtrial. 428 429 if (iwrite_fftdatar(mpi_enreg) .and. dtset%usewvl==0) then 430 431 ! output the density. 432 if (dtset%prtden/=0) then 433 if (dtset%positron/=1) rho_ptr => rhor 434 if (dtset%positron==1) rho_ptr => electronpositron%rhor_ep 435 call fftdatar_write("density",dtfil%fnameabo_app_den,dtset%iomode,hdr,& 436 crystal,ngfft,cplex1,nfft,nspden,rho_ptr,mpi_enreg,ebands=ebands) 437 438 if (dtset%positron/=0) then 439 if (dtset%positron/=1) rho_ptr => electronpositron%rhor_ep 440 if (dtset%positron==1) rho_ptr => rhor 441 fname = trim(dtfil%fnameabo_app_den)//'_POSITRON' 442 if (dtset%iomode == IO_MODE_ETSF) fname = strcat(fname, ".nc") 443 call fftdatar_write("positron_density",fname,dtset%iomode,hdr,& 444 crystal,ngfft,cplex1,nfft,nspden,rho_ptr,mpi_enreg,ebands=ebands) 445 end if 446 end if 447 448 else if (dtset%usewvl == 1 .and. dtset%prtden /= 0) then 449 !if iomode == 2 then set all outputs to netcdf format 450 !if iomode == 3 then set all outputs to ETSF format 451 accessfil = 0 452 if (dtset%iomode == IO_MODE_ETSF) accessfil = 3 453 if (dtset%iomode == IO_MODE_MPI) accessfil = 4 454 fform = fform_den 455 ! Write wavelet DEN. Note however that this should be delegate to separated Bigdft routines. 456 ! a lot of stuff written in outscf does not make sense if usewvl==0 457 call ioarr(accessfil,rhor,dtset,etotal,fform,dtfil%fnameabo_app_den, & 458 hdr,mpi_enreg,ngfft,cplex1,nfft,pawrhoij_dum,rdwr2,rdwrpaw0,wvl_den) 459 end if ! if master 460 461 !! MS - Printing of PAWDEN parallellised and several possible options included 462 !We output the total electron density in the PAW case 463 !this requires removing nhat from rhor and making PAW on-site corrections 464 if (pawprtden>0 .and. psps%usepaw==1) then 465 ! pawprtden 1 --> output PAW valence density 466 ! " 2 --> output PAW valence+core density 467 ! " 3 --> output core, valence and full atomic protodensity 468 ! " 4 --> options 1+3 469 ! " 5 --> options 2+3 470 ! " 6 --> output all individual PAW density contributions 471 if (pawprtden/=3) then ! calc PAW valence density 472 ABI_ALLOCATE(rhor_paw,(pawfgr%nfft,nspden)) 473 ABI_ALLOCATE(rhor_n_one,(pawfgr%nfft,nspden)) 474 ABI_ALLOCATE(rhor_nt_one,(pawfgr%nfft,nspden)) 475 ! If the communicator used for denfgr is kpt_comm, it is not compatible with paral_atom 476 if (mpi_enreg%paral_kgb==0.and.my_natom/=natom) then 477 my_natom_tmp=natom 478 ABI_DATATYPE_ALLOCATE(pawrhoij_all,(natom)) 479 call pawrhoij_nullify(pawrhoij_all) 480 call pawrhoij_copy(pawrhoij,pawrhoij_all,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 481 & keep_cplex=.false.,keep_itypat=.false.,keep_nspden=.false.) 482 else 483 my_natom_tmp=my_natom 484 pawrhoij_all => pawrhoij 485 end if 486 if (pawprtden/=6) then 487 call denfgr(atindx1,gmet,comm_fft,my_natom_tmp,natom,nattyp,ngfft,nhat,dtset%nspinor,nsppol,nspden,& 488 & ntypat,pawfgr,pawrad,pawrhoij_all,pawtab,prtvol,rhor,rhor_paw,rhor_n_one,& 489 & rhor_nt_one,rprimd,dtset%typat,ucvol,xred,& 490 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 491 else 492 call denfgr(atindx1,gmet,comm_fft,my_natom_tmp,natom,nattyp,ngfft,nhat,dtset%nspinor,nsppol,nspden,& 493 & ntypat,pawfgr,pawrad,pawrhoij_all,pawtab,prtvol,rhor,rhor_paw,rhor_n_one,& 494 & rhor_nt_one,rprimd,dtset%typat,ucvol,xred,& 495 & abs_n_tilde_nt_diff=nt_ntone_norm,znucl=dtset%znucl,& 496 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 497 end if 498 if (mpi_enreg%paral_kgb==0.and.my_natom/=natom) then 499 call pawrhoij_free(pawrhoij_all) 500 ABI_DATATYPE_DEALLOCATE(pawrhoij_all) 501 end if 502 503 if (prtvol>9) then ! Check normalisation 504 norm = SUM(rhor_paw(:,1))*ucvol/PRODUCT(pawfgr%ngfft(1:3)) 505 call xmpi_sum(norm,comm_fft,ierr) 506 write(message,'(a,F8.4)') ' PAWDEN - NORM OF DENSITY: ',norm 507 call wrtout(std_out,message,'COLL') 508 end if 509 end if 510 511 if (pawprtden>1.AND.pawprtden<6) then ! We will need the core density 512 ABI_ALLOCATE(rhor_paw_core,(pawfgr%nfft,nspden)) 513 call read_atomden(mpi_enreg,natom,pawfgr%nfft,pawfgr%ngfft,nspden,ntypat,rhor_paw_core,& 514 & dtset%typat,rprimd,xred,prtvol,file_prefix='core ') 515 516 if (prtvol>9) then ! Check normalisation 517 norm = SUM(rhor_paw_core(:,1))*ucvol/PRODUCT(pawfgr%ngfft(1:3)) 518 call xmpi_sum(norm,comm_fft,ierr) 519 write(message,'(a,F8.4)') ' ATMDEN - NORM OF CORE DENSITY: ', norm 520 call wrtout(std_out,message,'COLL') 521 end if 522 end if 523 524 if (pawprtden>2.AND.pawprtden<6) then ! We will need the valence protodensity 525 ABI_ALLOCATE(rhor_paw_val,(pawfgr%nfft,nspden)) 526 call read_atomden(mpi_enreg,natom,pawfgr%nfft,pawfgr%ngfft,nspden,ntypat,rhor_paw_val,& 527 & dtset%typat,rprimd,xred,prtvol,file_prefix='valence') 528 529 if (prtvol>9) then ! Check normalisation 530 norm = SUM(rhor_paw_val(:,1))*ucvol/PRODUCT(pawfgr%ngfft(1:3)) 531 call xmpi_sum(norm,comm_fft,ierr) 532 write(message,'(a,F8.4)') ' ATMDEN - NORM OF VALENCE PROTODENSITY: ', norm 533 call wrtout(std_out,message,'COLL') 534 end if 535 end if 536 537 if (iwrite_fftdatar(mpi_enreg)) then 538 if (pawprtden/=3) then 539 if (pawprtden==2.or.pawprtden==5) rhor_paw = rhor_paw + rhor_paw_core 540 ! PAWDEN 541 call fftdatar_write("pawrhor",dtfil%fnameabo_app_pawden,dtset%iomode,hdr,& 542 crystal,ngfft,cplex1,nfft,nspden,rhor_paw,mpi_enreg,ebands=ebands) 543 end if 544 545 if (pawprtden>2.AND.pawprtden<6) then 546 ! ATMDEN_CORE 547 call fftdatar_write("pawrhor_core",dtfil%fnameabo_app_atmden_core,dtset%iomode,hdr,& 548 crystal,ngfft,cplex1,nfft,nspden,rhor_paw_core,mpi_enreg,ebands=ebands) 549 550 ! valence protodensity. ATMDEN_VAL 551 call fftdatar_write("pawrhor_val",dtfil%fnameabo_app_atmden_val,dtset%iomode,hdr,& 552 crystal,ngfft,cplex1,nfft,nspden,rhor_paw_val,mpi_enreg,ebands=ebands) 553 554 ! full protodensity. ATMDEN_FULL 555 rhor_paw_val = rhor_paw_val + rhor_paw_core 556 call fftdatar_write("pawrhor_full",dtfil%fnameabo_app_atmden_full,dtset%iomode,hdr,& 557 crystal,ngfft,cplex1,nfft,nspden,rhor_paw_val,mpi_enreg,ebands=ebands) 558 end if 559 560 if (pawprtden==6) then ! Print all individual contributions to the density 561 ! N_TILDE - N_HAT 562 ! Use rhor_paw_val as temporary array 563 if (.not.allocated(rhor_paw_val)) then 564 ABI_ALLOCATE(rhor_paw_val,(pawfgr%nfft,nspden)) 565 end if 566 rhor_paw_val = rhor - nhat 567 568 call fftdatar_write("pawrhor_ntilde_minus_nhat",dtfil%fnameabo_app_n_tilde,dtset%iomode,hdr,& 569 crystal,ngfft,cplex1,nfft,nspden,rhor_paw_val,mpi_enreg,ebands=ebands) 570 571 ! N_ONSITE 572 call fftdatar_write("pawrhor_n_one",dtfil%fnameabo_app_n_one,dtset%iomode,hdr,& 573 crystal,ngfft,cplex1,nfft,nspden,rhor_n_one,mpi_enreg,ebands=ebands) 574 575 ! N_TILDE_ONSITE 576 call fftdatar_write("pawrhor_nt_one",dtfil%fnameabo_app_nt_one,dtset%iomode,hdr,& 577 crystal,ngfft,cplex1,nfft,nspden,rhor_nt_one,mpi_enreg,ebands=ebands) 578 579 end if ! All indivdual density cont. 580 end if ! if master 581 582 if (allocated(rhor_paw)) then 583 ABI_DEALLOCATE(rhor_paw) 584 end if 585 if (allocated(rhor_paw_core)) then 586 ABI_DEALLOCATE(rhor_paw_core) 587 end if 588 if (allocated(rhor_paw_val)) then 589 ABI_DEALLOCATE(rhor_paw_val) 590 end if 591 if (allocated(rhor_n_one)) then 592 ABI_DEALLOCATE(rhor_n_one) 593 end if 594 if (allocated(rhor_nt_one)) then 595 ABI_DEALLOCATE(rhor_nt_one) 596 end if 597 598 end if ! if paw+pawprtden 599 600 call timab(952,2,tsec) 601 call timab(953,1,tsec) 602 603 ! Output of the GSR file (except when we are inside mover) 604 #ifdef HAVE_NETCDF 605 if (me == master .and. dtset%prtgsr==1 .and. dtset%usewvl == 0) then 606 !.and. (dtset%ionmov /= 0 .or. dtset%optcell /= 0)) then 607 fname = strcat(dtfil%filnam_ds(4), "_GSR.nc") 608 609 ! Write density, crystal, band structure energies. 610 NCF_CHECK(nctk_open_create(ncid, fname, xmpi_comm_self)) 611 NCF_CHECK(hdr_ncwrite(hdr, ncid, fform_den, nc_define=.True.)) 612 NCF_CHECK(crystal_ncwrite(crystal, ncid)) 613 NCF_CHECK(ebands_ncwrite(ebands, ncid)) 614 ! Add energy, forces, stresses 615 NCF_CHECK(results_gs_ncwrite(results_gs, ncid, dtset%ecut, dtset%pawecutdg)) 616 NCF_CHECK(nf90_close(ncid)) 617 end if 618 #endif 619 620 ! Output of VCLMB file 621 ! The PAW correction has to be computed here (all processors contribute) 622 if (psps%usepaw > 0 .AND. dtset%prtvclmb>0) then 623 nradint = 1000 ! radial integration grid density 624 ABI_ALLOCATE(vpaw,(nfft,nspden)) 625 vpaw(:,:)=zero 626 if (me == master .and. my_natom > 0) then 627 if (paw_an(1)%cplex > 1) then 628 MSG_WARNING('cplex = 2 : complex hartree potential in PAW spheres. This is not coded yet. Imag part ignored') 629 end if 630 end if 631 632 do ispden=1,nspden 633 ! for points inside spheres, replace with full AE hartree potential. 634 ! In principle the correction could be more subtle (not spherical) 635 do iatom=1,my_natom 636 iatom_tot=iatom;if (paral_atom) iatom_tot=mpi_enreg%my_atmtab(iatom) 637 itypat=dtset%typat(iatom_tot) 638 639 ABI_ALLOCATE(vh1spl,(paw_an(iatom)%mesh_size)) 640 ABI_ALLOCATE(vh1_corrector,(paw_an(iatom)%mesh_size)) 641 ABI_ALLOCATE(vh1_interp,(pawfgrtab(iatom)%nfgd)) 642 ABI_ALLOCATE(radii,(pawfgrtab(iatom)%nfgd)) 643 ABI_ALLOCATE(isort,(pawfgrtab(iatom)%nfgd)) 644 ! vh1 vht1 contain the spherical first moments of the Hartree potentials, so re-divide by Y_00 = sqrt(four_pi) 645 vh1_corrector(:) = (paw_an(iatom)%vh1(:,1,ispden)-paw_an(iatom)%vht1(:,1,ispden)) / sqrt(four_pi) 646 647 ! get end point derivatives 648 call bound_deriv(vh1_corrector, pawrad(itypat), pawrad(itypat)%mesh_size, yp1, ypn) 649 ! spline the vh1 function 650 ! NB for second argument of vh1: only first moment lm_size appears to be used 651 ! NB2: vh1 can in principle be complex - not sure what to do with the imaginary part. Ignored for now. 652 call spline(pawrad(itypat)%rad, vh1_corrector, paw_an(iatom)%mesh_size, yp1, ypn, vh1spl) 653 654 do ifgd = 1, pawfgrtab(iatom)%nfgd 655 ! get radii for this point 656 isort(ifgd) = ifgd 657 radii(ifgd) = sqrt(sum(pawfgrtab(iatom)%rfgd(:,ifgd)**2)) 658 end do 659 660 if (pawfgrtab(iatom)%nfgd/=0) then 661 ! spline interpolate the vh1 value for current radii 662 call sort_dp(pawfgrtab(iatom)%nfgd, radii, isort, tol12) 663 call splint(pawrad(itypat)%mesh_size, pawrad(itypat)%rad, & 664 & vh1_corrector, vh1spl, pawfgrtab(iatom)%nfgd, radii, vh1_interp, ierr) 665 end if 666 667 norm=SUM(vh1_interp)*ucvol/PRODUCT(ngfft(1:3)) 668 call xmpi_sum(norm,comm_fft,ierr) 669 write(message,'(a,i6,a,E20.10)') ' sum of Hartree correction term on fft grid of atom : ', iatom, & 670 & ' = ', norm 671 call wrtout(std_out,message,'COLL') 672 673 if (pawfgrtab(iatom)%nfgd/=0) then 674 vpaw(pawfgrtab(iatom)%ifftsph(isort(1:pawfgrtab(iatom)%nfgd)),ispden) = & 675 & vpaw(pawfgrtab(iatom)%ifftsph(isort(1:pawfgrtab(iatom)%nfgd)),ispden) + & 676 & vh1_interp(1:pawfgrtab(iatom)%nfgd) 677 end if 678 679 ! get integral of correction term in whole sphere 680 ABI_DEALLOCATE(radii) 681 ABI_DEALLOCATE(vh1_interp) 682 683 ABI_ALLOCATE(radii,(nradint)) 684 ABI_ALLOCATE(vh1_interp,(nradint)) 685 686 ABI_ALLOCATE(vh1_integ,(nradint)) 687 dr = pawrad(itypat)%rad(paw_an(iatom)%mesh_size) / dble(nradint) 688 do ifgd = 1, nradint 689 radii(ifgd) = dble(ifgd-1)*dr 690 end do 691 692 ! spline interpolate the vh1 value for current radii 693 call splint(pawrad(itypat)%mesh_size, pawrad(itypat)%rad, & 694 & vh1_corrector, vh1spl, nradint, radii, vh1_interp, ierr) 695 696 do ifgd = 1, nradint 697 vh1_interp(ifgd) = vh1_interp(ifgd)*radii(ifgd)**2 698 end do 699 700 call simpson_int(nradint, dr, vh1_interp, vh1_integ) 701 write(message,'(a,i6,a,E20.10)') ' integral of Hartree correction term in sphere of atom: ', iatom, & 702 & ' = ', vh1_integ(nradint)*four*pi 703 call wrtout(std_out,message,'COLL') 704 705 ABI_DEALLOCATE(vh1spl) 706 ABI_DEALLOCATE(vh1_corrector) 707 ABI_DEALLOCATE(vh1_interp) 708 ABI_DEALLOCATE(vh1_integ) 709 ABI_DEALLOCATE(radii) 710 ABI_DEALLOCATE(isort) 711 end do ! iatom 712 end do !ispden 713 call xmpi_sum_master(vpaw,master,mpi_enreg%comm_atom,ierr) 714 if (.not.iwrite_fftdatar(mpi_enreg)) then 715 ABI_DEALLOCATE(vpaw) 716 end if 717 end if ! if paw - add all electron vhartree in spheres 718 719 if (iwrite_fftdatar(mpi_enreg)) then 720 721 ! output the electron localization function ELF 722 if (dtset%prtelf/=0) then 723 call fftdatar_write("elfr",dtfil%fnameabo_app_elf,dtset%iomode,hdr,& 724 crystal,ngfft,cplex1,nfft,nspden,elfr,mpi_enreg,ebands=ebands) 725 726 if (nspden==2)then 727 ABI_ALLOCATE(elfr_up,(nfft,nspden)) 728 elfr_up(:,:) = zero 729 do ifft=1,nfft 730 elfr_up(ifft,1) = elfr(ifft,2) 731 end do 732 ! ELF_UP 733 call fftdatar_write("elfr_up",dtfil%fnameabo_app_elf_up,dtset%iomode,hdr,& 734 crystal,ngfft,cplex1,nfft,nspden,elfr_up,mpi_enreg,ebands=ebands) 735 736 ABI_ALLOCATE(elfr_down,(nfft,nspden)) 737 elfr_down(:,:) = zero 738 do ifft=1,nfft 739 elfr_down(ifft,1) = elfr(ifft,3) 740 end do 741 ! ELF_DOWN' 742 call fftdatar_write("elfr_down",dtfil%fnameabo_app_elf_down,dtset%iomode,hdr,& 743 crystal,ngfft,cplex1,nfft,nspden,elfr_down,mpi_enreg,ebands=ebands) 744 745 ABI_DEALLOCATE(elfr_up) 746 ABI_DEALLOCATE(elfr_down) 747 end if 748 end if 749 750 call timab(953,2,tsec) 751 call timab(954,1,tsec) 752 753 ! We output the gradient of density 754 if (dtset%prtgden/=0) then 755 756 call fftdatar_write("grhor_1",dtfil%fnameabo_app_gden1,dtset%iomode,hdr,& 757 crystal,ngfft,cplex1,nfft,nspden,grhor(:,:,1),mpi_enreg,ebands=ebands) 758 759 call fftdatar_write("grhor_2",dtfil%fnameabo_app_gden2,dtset%iomode,hdr,& 760 crystal,ngfft,cplex1,nfft,nspden,grhor(:,:,2),mpi_enreg,ebands=ebands) 761 762 call fftdatar_write("grhor_3",dtfil%fnameabo_app_gden3,dtset%iomode,hdr,& 763 crystal,ngfft,cplex1,nfft,nspden,grhor(:,:,3),mpi_enreg,ebands=ebands) 764 end if 765 766 ! We output the total kinetic energy density KDEN 767 if (dtset%prtkden/=0) then 768 call fftdatar_write("kinedr",dtfil%fnameabo_app_kden,dtset%iomode,hdr,& 769 crystal,ngfft,cplex1,nfft,nspden,taur,mpi_enreg,ebands=ebands) 770 end if 771 772 ! We output the Laplacian of density 773 if (dtset%prtlden/=0) then 774 call fftdatar_write("laprhor",dtfil%fnameabo_app_lden,dtset%iomode,hdr,& 775 crystal,ngfft,cplex1,nfft,nspden,lrhor,mpi_enreg,ebands=ebands) 776 end if 777 778 call timab(954,2,tsec) 779 call timab(955,1,tsec) 780 781 call timab(955,2,tsec) 782 call timab(956,1,tsec) 783 784 ! POT 785 if (dtset%prtpot>0) then 786 call fftdatar_write("vtrial",dtfil%fnameabo_app_pot,dtset%iomode,hdr,& 787 crystal,ngfft,cplex1,nfft,nspden,vtrial,mpi_enreg,ebands=ebands) 788 end if 789 790 call timab(956,2,tsec) 791 call timab(957,1,tsec) 792 793 if (dtset%prtgeo>0) then 794 coordn=dtset%prtgeo 795 call bonds_lgth_angles(coordn,dtfil%fnameabo_app_geo,natom,psps%ntypat,& 796 & rprimd,dtset%typat,xred,dtset%znucl) 797 end if 798 799 if (dtset%prtcif > 0) then 800 call prt_cif(dtset%brvltt, dtfil%fnameabo_app_cif, natom, dtset%nsym, dtset%ntypat, rprimd, & 801 & dtset%spgaxor, dtset%spgroup, dtset%spgorig, dtset%symrel, dtset%tnons, dtset%typat, xred, dtset%znucl) 802 end if 803 804 call timab(957,2,tsec) 805 call timab(958,1,tsec) 806 807 ! STM 808 if (dtset%prtstm>0) then 809 call fftdatar_write("stm",dtfil%fnameabo_app_stm,dtset%iomode,hdr,& 810 crystal,ngfft,cplex1,nfft,nspden,rhor,mpi_enreg,ebands=ebands) 811 end if 812 813 if (dtset%prt1dm>0) then 814 call out1dm(dtfil%fnameabo_app_1dm,mpi_enreg,natom,nfft,ngfft,nspden,psps%ntypat,& 815 & rhor,rprimd,dtset%typat,ucvol,vtrial,xred,dtset%znucl) 816 end if 817 818 ! VHA 819 if (dtset%prtvha>0) then 820 ABI_ALLOCATE(vwork,(nfft,nspden)) 821 do ispden=1,nspden 822 vwork(:,ispden)=vhartr(:) 823 end do 824 825 call fftdatar_write("vhartree",dtfil%fnameabo_app_vha,dtset%iomode,hdr,& 826 crystal,ngfft,cplex1,nfft,nspden,vwork,mpi_enreg,ebands=ebands) 827 828 ABI_DEALLOCATE(vwork) 829 end if 830 831 ! VPSP 832 if (dtset%prtvpsp>0) then 833 ABI_ALLOCATE(vwork,(nfft,nspden)) 834 do ispden=1,nspden 835 vwork(:,ispden)=vpsp(:) 836 end do 837 838 call fftdatar_write("vpsp",dtfil%fnameabo_app_vpsp,dtset%iomode,hdr,& 839 crystal,ngfft,cplex1,nfft,nspden,vwork,mpi_enreg,ebands=ebands) 840 841 ABI_DEALLOCATE(vwork) 842 end if 843 844 ! VCouLoMB 845 if (dtset%prtvclmb>0) then 846 847 ABI_ALLOCATE(vwork,(nfft,nspden)) 848 do ispden=1,nspden 849 vwork(:,ispden)=vpsp(:)+vhartr(:) 850 end do 851 if (psps%usepaw==1) then 852 do ispden=1,nspden 853 vwork(:,ispden)=vwork(:,ispden)+vpaw(:,ispden) 854 end do 855 ABI_DEALLOCATE(vpaw) 856 end if 857 858 call fftdatar_write("vhartree_vloc",dtfil%fnameabo_app_vclmb,dtset%iomode,hdr,& 859 & crystal,ngfft,cplex1,nfft,nspden,vwork,mpi_enreg,ebands=ebands) 860 861 !TODO: find out why this combination of calls with fftdatar_write then out1dm fails on buda with 4 mpi-fft procs (npkpt 1). 862 ! For the moment comment it out. Only DS2 of mpiio test 27 fails 863 ! call out1dm(dtfil%fnameabo_app_vclmb_1dm,mpi_enreg,natom,nfft,ngfft,nspden,psps%ntypat,& 864 !& rhor,rprimd,dtset%typat,ucvol,vwork,xred,dtset%znucl) 865 866 ! 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 867 ! where E is energy of electron, E0 rest mass, lambda the relativistic wavelength 868 ! values of CE at 200 300 and 1000 kV: 7.29e6 6.53e6 5.39e6 rad / V / m 869 ! vertical integral of vclmb * c / ngfft(3) / cross sectional area factor (= sin(gamma)) 870 ! * 0.5291772083e-10*27.2113834 to get to SI 871 ! * CE factor above 872 ! should be done for each plane perpendicular to the axes... 873 ABI_DEALLOCATE(vwork) 874 end if ! prtvclmb 875 876 877 ! VHXC 878 if (dtset%prtvhxc>0) then 879 ABI_ALLOCATE(vwork,(nfft,nspden)) 880 do ispden=1,nspden 881 vwork(:,ispden)=vhartr(:)+vxc(:,ispden) 882 end do 883 884 call fftdatar_write("vhxc",dtfil%fnameabo_app_vhxc,dtset%iomode,hdr,& 885 & crystal,ngfft,cplex1,nfft,nspden,vwork,mpi_enreg,ebands=ebands) 886 887 ABI_DEALLOCATE(vwork) 888 end if 889 890 ! VXC 891 if (dtset%prtvxc>0) then 892 call fftdatar_write("exchange_correlation_potential",dtfil%fnameabo_app_vxc,dtset%iomode,hdr,& 893 crystal,ngfft,cplex1,nfft,nspden,vxc,mpi_enreg,ebands=ebands) 894 end if 895 896 call timab(958,2,tsec) 897 end if ! if iwrite_fftdatar 898 899 call timab(959,1,tsec) 900 901 !Generate DOS using the tetrahedron method or using Gaussians 902 !FIXME: Should centralize all calculations of DOS here in outscfcv 903 if (dtset%prtdos>=2.or.dtset%pawfatbnd>0) then 904 dos = epjdos_new(dtset, psps, pawtab) 905 906 if (dos%partial_dos_flag>=1 .or. dos%fatbands_flag==1)then 907 ! Generate fractions for partial DOSs if needed partial_dos 1,2,3,4 give different decompositions 908 collect = 1 !; if (psps%usepaw==1 .and. dos%partial_dos_flag /= 2) collect = 0 909 if ((psps%usepaw==0.or.dtset%pawprtdos/=2) .and. dos%partial_dos_flag>=1) then 910 call partial_dos_fractions(dos,crystal,dtset,eigen,occ,npwarr,kg,cg,mcg,collect,mpi_enreg) 911 end if 912 913 if (psps%usepaw==1 .and. dos%partial_dos_flag /= 2) then 914 ! TODO: update partial_dos_fractions_paw for extra atoms - no PAW contribution normally, but check bounds and so on. 915 call partial_dos_fractions_paw(dos,cprj,dimcprj,dtset,mcprj,mkmem,mpi_enreg,pawrad,pawtab) 916 end if 917 918 else 919 dos%fractions(:,:,:,1)=one 920 end if 921 922 ! Here, print out fatbands for the k-points given in file appended _FATBANDS 923 if (me == master .and. dtset%pawfatbnd>0 .and. dos%fatbands_flag==1) then 924 call prtfatbands(dos,dtset,ebands,dtfil%fnameabo_app_fatbands,dtset%pawfatbnd,pawtab) 925 end if 926 927 ! Here, computation and output of DOS and partial DOS _DOS 928 if (dos%fatbands_flag == 0) then 929 if (dos%prtdos /= 4) then 930 call dos_calcnwrite(dos,dtset,crystal,ebands,dtfil%fnameabo_app_dos,spacecomm) 931 end if 932 end if 933 934 #ifdef HAVE_NETCDF 935 ! Write netcdf file with dos% results. 936 if (me == master) then 937 fname = trim(dtfil%filnam_ds(4))//'_FATBANDS.nc' 938 NCF_CHECK(nctk_open_create(ncid, fname, xmpi_comm_self)) 939 call fatbands_ncwrite(dos, crystal, ebands, hdr, dtset, psps, pawtab, ncid) 940 NCF_CHECK(nf90_close(ncid)) 941 end if 942 #endif 943 944 call epjdos_free(dos) 945 end if ! prtdos > 1 946 947 call timab(959,2,tsec) 948 call timab(960,1,tsec) 949 950 !Output of integrated density inside atomic spheres 951 if (dtset%prtdensph==1.and.dtset%usewvl==0)then 952 call calcdensph(gmet,mpi_enreg,natom,nfft,ngfft,nspden,& 953 & ntypat,ab_out,dtset%ratsph,rhor,rprimd,dtset%typat,ucvol,xred,1,cplex1) 954 end if 955 956 call timab(960,2,tsec) 957 958 if (dtset%magconon /= 0) then 959 ! calculate final value of terms for magnetic constraint: "energy" term, lagrange multiplier term, and atomic contributions 960 call mag_constr_e(dtset%magconon,dtset%magcon_lambda,mpi_enreg,& 961 & natom,nfft,ngfft,nspden,ntypat,dtset%ratsph,rhor,rprimd,dtset%spinat,dtset%typat,xred) 962 end if 963 964 call timab(961,1,tsec) 965 966 !If PAW, provide additional outputs 967 if (psps%usepaw==1) then 968 ! Output of compensation charge 969 if (dtset%nstep>0.or.dtfil%ireadwf/=0) then 970 write(message, '(4a)' )ch10,' PAW TEST:',ch10,& 971 & ' ==== Compensation charge inside spheres ============' 972 if (compch_sph>-1.d4.and.compch_fft>-1.d4) & 973 & write(message, '(3a)' ) trim(message),ch10,' The following values must be close to each other ...' 974 if (compch_sph>-1.d4) write(message, '(3a,f22.15)' ) trim(message),ch10,& 975 & ' Compensation charge over spherical meshes = ',compch_sph 976 if (compch_fft>-1.d4) then 977 if (pawfgr%usefinegrid==1) then 978 write(message, '(3a,f22.15)' ) trim(message),ch10,& 979 & ' Compensation charge over fine fft grid = ',compch_fft 980 else 981 write(message, '(3a,f22.15)' ) trim(message),ch10,& 982 & ' Compensation charge over fft grid = ',compch_fft 983 end if 984 end if 985 call wrtout(ab_out,message,'COLL') 986 call wrtout(std_out,message,'COLL') 987 end if 988 ! Output of pseudopotential strength Dij and augmentation occupancies Rhoij 989 call pawprt(dtset,my_natom,paw_ij,pawrhoij,pawtab,& 990 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 991 & electronpositron=electronpositron) 992 end if 993 994 call timab(961,2,tsec) 995 call timab(962,1,tsec) 996 997 998 !PAW + output for optical conductivity _OPT and _OPT2 999 if (psps%usepaw==1.and.prtnabla>0) then 1000 if (prtnabla==1.or.prtnabla==2) then 1001 call optics_paw(atindx1,cg,cprj,dimcprj,dtfil,dtset,eigen,gprimd,hdr,kg,& 1002 & mband,mcg,mcprj,mkmem,mpi_enreg,mpsang,mpw,natom,nkpt,npwarr,nsppol,pawrad,pawtab) 1003 end if 1004 if (prtnabla==2.or.prtnabla==3) then 1005 call optics_paw_core(atindx1,cprj,dimcprj,dtfil,dtset,eigen,psps%filpsp,hdr,& 1006 & mband,mcprj,mkmem,mpi_enreg,mpsang,natom,nkpt,nsppol,pawrad,pawtab) 1007 end if 1008 end if 1009 if (prtnabla<0) then 1010 ! TODO: This routine is not tested but it's used in production. 1011 call optics_vloc(cg,dtfil,dtset,eigen,gprimd,hdr,kg,& 1012 & mband,mcg,mkmem,mpi_enreg,mpw,nkpt,npwarr,nsppol) 1013 end if 1014 1015 call timab(962,2,tsec) 1016 call timab(963,1,tsec) 1017 1018 !Optionally provide output for AE wavefunctions (only for PAW) 1019 if (psps%usepaw==1 .and. dtset%pawprtwf==1) then 1020 ABI_ALLOCATE(ps_norms,(nsppol,nkpt,mband)) 1021 1022 call pawmkaewf(dtset,crystal,ebands,my_natom,mpw,mband,mcg,mcprj,nkpt,mkmem,nsppol,Dtset%nband,& 1023 & Dtset%istwfk,npwarr,Dtset%kptns,Dtset%ngfftdg,kg,dimcprj,pawfgrtab,& 1024 & Pawrad,Pawtab,Hdr,Dtfil,cg,Cprj,& 1025 & MPI_enreg,ierr,pseudo_norms=ps_norms,set_k=dtset%pawprt_k,set_band=dtset%pawprt_b,& 1026 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1027 1028 if (dtset%pawprt_b==0) then 1029 fname = strcat(dtfil%filnam_ds(4), '_PAWSTAT') 1030 if (open_file(fname, message,newunit=tmp_unt,status='unknown',form='formatted') /= 0) then 1031 MSG_ERROR(message) 1032 end if 1033 write(tmp_unt,'(5a)') '# This file contains the statistics on the cancellation of',ch10,& 1034 & '# the onsite pseudo component of the all-electron wavefunction',ch10,& 1035 & '# with the plane wave part' 1036 ii = 0 1037 do isppol=1,nsppol 1038 write(tmp_unt,'(a,i0)') '# isppol = ',isppol 1039 do ikpt=1,nkpt 1040 write(tmp_unt,'(a,i0)') '# ikpt = ',ikpt 1041 write(tmp_unt,'(a)') '# band norm' 1042 occ_norm = zero; unocc_norm = zero; nocc = 0 1043 do iband=1,dtset%nband(ikpt + (isppol-1)*nkpt) 1044 ii = ii + 1 1045 write(tmp_unt,'(i8,ES16.6)') iband,ps_norms(isppol,ikpt,iband) 1046 if (abs(occ(ii)) <= tol16) then 1047 unocc_norm = unocc_norm + ps_norms(isppol,ikpt,iband) 1048 else 1049 occ_norm = occ_norm + ps_norms(isppol,ikpt,iband) 1050 nocc = nocc + 1 1051 end if 1052 end do 1053 if(mband/=nocc)then 1054 write(tmp_unt,'(2(a,ES16.6))') '# occ average: ',occ_norm/real(nocc),& 1055 & ' unocc average: ',unocc_norm/real(mband-nocc) 1056 else 1057 write(tmp_unt,'(2(a,ES16.6))') '# occ average: ',occ_norm/real(nocc) 1058 end if 1059 end do 1060 end do 1061 close(tmp_unt) 1062 end if 1063 ABI_DEALLOCATE(ps_norms) 1064 end if 1065 1066 call timab(963,2,tsec) 1067 if(dtset%plowan_compute>0) then 1068 write(message,'(2a,i3)') ch10,& 1069 & ' ====================================================================================== ' 1070 call wrtout(std_out,message,'COLL') 1071 call wrtout(ab_out,message,'COLL') 1072 write(message,'(2a,i3)') ch10,& 1073 & ' == Start computation of Projected Local Orbitals Wannier functions == ',dtset%nbandkss 1074 call wrtout(std_out,message,'COLL') 1075 call wrtout(ab_out,message,'COLL') 1076 1077 ! == compute psichi 1078 1079 call init_plowannier(dtset,wan) 1080 call compute_coeff_plowannier(crystal,cprj,dimcprj,dtset,eigen,e_fermie,& 1081 & mpi_enreg,occ,wan,pawtab,psps,usecprj,dtfil%unpaw,pawrad,dtfil) 1082 call destroy_plowannier(wan) 1083 end if 1084 1085 !Optionally provide output for the GW part of ABINIT 1086 if (dtset%nbandkss/=0) then 1087 ! Use DMFT to compute wannier function for cRPA calculation. 1088 if(dtset%usedmft==1) then 1089 write(message,'(2a,i3)') ch10,& 1090 & ' Warning: Psichi are renormalized in datafordmft because nbandkss is used',dtset%nbandkss 1091 call wrtout(std_out,message,'COLL') 1092 call init_dmft(dmatpawu,dtset,e_fermie,dtfil%fnameabo_app,dtset%nspinor,paw_dmft,pawtab,psps,dtset%typat) 1093 call print_dmft(paw_dmft,dtset%pawprtvol) 1094 1095 ! == compute psichi 1096 call init_oper(paw_dmft,lda_occup) 1097 1098 call datafordmft(crystal,cprj,dimcprj,dtset,eigen,e_fermie,& 1099 & lda_occup,dtset%mband,dtset%mband,dtset%mkmem,mpi_enreg,& 1100 & dtset%nkpt,dtset%nspinor,dtset%nsppol,occ,& 1101 & paw_dmft,paw_ij,pawang,pawtab,psps,usecprj,dtfil%unpaw,dtset%nbandkss) 1102 1103 call destroy_dmft(paw_dmft) 1104 call destroy_oper(lda_occup) 1105 end if 1106 1107 call timab(964,1,tsec) ! outscfcv(outkss) 1108 1109 call outkss(crystal,dtfil,dtset,ecut,gmet,gprimd,hdr,& 1110 & dtset%kssform,mband,mcg,mcprj,mgfftc,mkmem,mpi_enreg,mpsang,mpw,natom,natom,& 1111 & nfft,nkpt,npwarr,nspden,nsppol,nsym,psps%ntypat,occ,pawtab,pawfgr,paw_ij,& 1112 & prtvol,psps,rprimd,vtrial,xred,cg,usecprj,cprj,eigen,ierr) 1113 call timab(964,2,tsec) ! outscfcv(outkss) 1114 if (ierr/=0) then 1115 MSG_WARNING("outkss returned a non zero status error, check log") 1116 end if 1117 end if 1118 1119 if (electronpositron_calctype(electronpositron)/=0) then 1120 1121 ! Optionally provide output for positron life time calculation 1122 call timab(965,1,tsec) 1123 call poslifetime(dtset,electronpositron,gprimd,my_natom,& 1124 & mpi_enreg,n3xccc,nfft,ngfft,nhat,1,pawang,& 1125 & pawrad,pawrhoij,pawtab,rate_dum,rate_dum2,& 1126 & rhor,ucvol,xccc3d) 1127 call timab(965,2,tsec) 1128 1129 ! Optionally provide output for momentum distribution of annihilation radiation 1130 if (dtset%posdoppler>0) then 1131 call posdoppler(cg,cprj,crystal,dimcprj,dtfil,dtset,electronpositron,psps%filpsp,& 1132 & kg,mcg,mcprj,mpi_enreg,my_natom,n3xccc,nfft,ngfft,nhat,npwarr,& 1133 & occ,pawang,pawrad,pawrhoij,pawtab,rhor,xccc3d) 1134 end if 1135 end if 1136 1137 !Optionally provide output for WanT 1138 if (dtset%prtwant==1) then 1139 call timab(966,1,tsec) 1140 ! WARNING: mpi_enreg not used --> MPI is not supported 1141 call outwant(dtset,eigen,cg,kg,npwarr,mband,mcg,nkpt,nsppol,mkmem,mpw,dtset%prtwant) 1142 call timab(966,2,tsec) 1143 end if 1144 1145 !Optionally provide output for electric field gradient calculation 1146 if (dtset%prtefg > 0) then 1147 call timab(967,1,tsec) 1148 call calc_efg(mpi_enreg,my_natom,natom,nfft,ngfft,nspden,dtset%nsym,ntypat,& 1149 & dtset%paral_kgb,paw_an,pawang,pawrad,pawrhoij,pawtab,& 1150 & dtset%ptcharge,dtset%prtefg,dtset%quadmom,rhor,rprimd,dtset%symrel,& 1151 & dtset%tnons,dtset%typat,ucvol,psps%usepaw,xred,psps%zionpsp,& 1152 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1153 call timab(967,2,tsec) 1154 end if 1155 1156 !Optionally provide output for Fermi-contact term at nuclear positions 1157 if (dtset%prtfc > 0) then 1158 call timab(967,1,tsec) 1159 call calc_fc(my_natom,natom,nspden,ntypat,pawrad,pawrhoij,pawtab,dtset%typat,psps%usepaw,& 1160 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1161 call timab(967,2,tsec) 1162 end if 1163 1164 ! Output electron bands. 1165 if (me == master .and. dtset%tfkinfunc==0) then 1166 if (size(dtset%kptbounds, dim=2) > 0) then 1167 call ebands_write(ebands, dtset%prtebands, dtfil%filnam_ds(4), kptbounds=dtset%kptbounds) 1168 else 1169 call ebands_write(ebands, dtset%prtebands, dtfil%filnam_ds(4)) 1170 end if 1171 end if 1172 1173 !Optionally provide Xcrysden output for the Fermi surface (Only master writes) 1174 if (me == master .and. dtset%prtfsurf == 1) then 1175 if (ebands_write_bxsf(ebands,crystal,dtfil%fnameabo_app_bxsf) /= 0) then 1176 message = "Cannot produce BXSF file with Fermi surface, see log file for more info" 1177 MSG_WARNING(message) 1178 call wrtout(ab_out, message) 1179 end if 1180 end if ! prtfsurf==1 1181 1182 !output nesting factor for Fermi surface (requires ph_nqpath) 1183 if (me == master .and. dtset%prtnest>0 .and. dtset%ph_nqpath > 0) then 1184 call timab(968,1,tsec) 1185 ierr = ebands_write_nesting(ebands,crystal,dtfil%fnameabo_app_nesting,dtset%prtnest,& 1186 & dtset%tsmear,dtset%fermie_nest,dtset%ph_qpath(:,1:dtset%ph_nqpath),message) 1187 if (ierr /=0) then 1188 MSG_WARNING(message) 1189 call wrtout(ab_out,message,'COLL') 1190 end if 1191 call timab(968,2,tsec) 1192 end if ! prtnest=1 1193 1194 call timab(969,1,tsec) 1195 1196 if (dtset%prtdipole == 1) then 1197 call multipoles_out(rhor,mpi_enreg,natom,nfft,ngfft,dtset%nspden,dtset%ntypat,rprimd,& 1198 & dtset%typat,ucvol,ab_out,xred,dtset%ziontypat) 1199 end if 1200 1201 ! BoltzTraP output files in GENEric format 1202 if (dtset%prtbltztrp == 1 .and. me==master) then 1203 call ebands_prtbltztrp(ebands, crystal, dtfil%filnam_ds(4)) 1204 end if 1205 1206 ! Band structure interpolation from eigenvalues computed on the k-mesh. 1207 if (nint(dtset%einterp(1)) /= 0) then 1208 call ebands_interpolate_kpath(ebands, dtset, crystal, [0,0], dtfil%filnam_ds(4), spacecomm) 1209 end if 1210 1211 call crystal_free(crystal) 1212 call ebands_free(ebands) 1213 1214 !Destroy atom table used for parallelism 1215 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 1216 1217 call timab(969,2,tsec) 1218 call timab(950,2,tsec) ! outscfcv 1219 1220 DBG_EXIT("COLL") 1221 1222 end subroutine outscfcv