TABLE OF CONTENTS
ABINIT/m_outscfcv [ Modules ]
NAME
m_outscfcv
FUNCTION
COPYRIGHT
Copyright (C) 2005-2024 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 .
SOURCE
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 #ifdef HAVE_NETCDF 32 use netcdf 33 #endif 34 use m_nctk 35 use m_hdr 36 use m_plowannier 37 use m_splines 38 use m_ebands 39 use m_dtset 40 use m_dtfil 41 42 use defs_datatypes, only : pseudopotential_type, ebands_t 43 use defs_abitypes, only : MPI_type 44 use m_time, only : timab 45 use m_io_tools, only : open_file 46 use m_fstrings, only : strcat, endswith 47 use m_geometry, only : bonds_lgth_angles 48 use m_electronpositron, only : electronpositron_type,electronpositron_calctype 49 use m_oper, only : oper_type,init_oper,destroy_oper 50 use m_crystal, only : crystal_init, crystal_t, prt_cif 51 use m_results_gs, only : results_gs_type, results_gs_ncwrite 52 use m_ioarr, only : ioarr, fftdatar_write 53 use m_nucprop, only : calc_efg,calc_fc 54 use m_outwant, only : outwant 55 use m_pawang, only : pawang_type 56 use m_pawrad, only : pawrad_type, simp_gen, bound_deriv 57 use m_pawtab, only : pawtab_type 58 use m_paw_an, only : paw_an_type 59 use m_paw_ij, only : paw_ij_type 60 use m_paw_mkrho, only : denfgr 61 use m_pawfgrtab, only : pawfgrtab_type 62 use m_pawrhoij, only : pawrhoij_type, pawrhoij_nullify, pawrhoij_copy, pawrhoij_free 63 use m_pawcprj, only : pawcprj_type 64 use m_pawfgr, only : pawfgr_type 65 use m_paw_dmft, only : paw_dmft_type,init_dmft,destroy_dmft,print_dmft 66 use m_paw_optics, only : optics_paw,optics_paw_core 67 use m_paw_tools, only : pawprt 68 use m_numeric_tools, only : simpson_int 69 use m_epjdos, only : dos_calcnwrite, partial_dos_fractions, partial_dos_fractions_paw, & 70 epjdos_t, epjdos_new, prtfatbands, fatbands_ncwrite 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_datafordmft, only : datafordmft 79 use m_mkrho, only : read_atomden 80 use m_positron, only : poslifetime, posdoppler 81 use m_optics_vloc, only : optics_vloc 82 use m_green, only : green_type,compute_green,& 83 fourier_green,print_green,init_green,destroy_green,init_green_tau 84 use m_self, only : self_type,initialize_self,rw_self,destroy_self,destroy_self,selfreal2imag_self 85 use m_paw_correlations, only : loc_orbmom_cal 86 87 implicit none 88 89 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 intgres(nspden,natom)=integrated residuals from constrained DFT. They are also Lagrange parameters, or gradients with respect to constraints. 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.
SOURCE
195 subroutine outscfcv(atindx1,cg,compch_fft,compch_sph,cprj,dimcprj,dmatpawu,dtfil,dtset,& 196 & ecut,eigen,electronpositron,elfr,etotal,gmet,gprimd,grhor,hdr,intgres,kg,& 197 & lrhor,mband,mcg,mcprj,mgfftc,mkmem,mpi_enreg,mpsang,mpw,my_natom,natom,& 198 & nattyp,nfft,ngfft,nhat,nkpt,npwarr,nspden,nsppol,nsym,ntypat,n3xccc,occ,& 199 & paw_dmft,pawang,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,paw_an,paw_ij,& 200 & prtvol,psps,results_gs,rhor,rprimd,& 201 & taur,ucvol,usecprj,vhartr,vpsp,vtrial,vxc,wvl_den,xccc3d,xred) 202 203 !Arguments ------------------------------------ 204 !scalars 205 integer,intent(in) :: mband,mcg,mcprj,mgfftc,mkmem,mpsang,mpw,n3xccc,my_natom,natom,nfft 206 integer,intent(in) :: nkpt,nspden,nsppol,nsym,ntypat,prtvol,usecprj 207 real(dp),intent(in) :: compch_fft,compch_sph,ecut,ucvol 208 real(dp),intent(inout) :: etotal 209 type(electronpositron_type),pointer :: electronpositron 210 type(MPI_type),intent(inout) :: mpi_enreg 211 type(datafiles_type),intent(in) :: dtfil 212 type(dataset_type),intent(in) :: dtset 213 type(hdr_type),intent(inout) :: hdr 214 type(paw_dmft_type), intent(inout) :: paw_dmft 215 type(pawang_type),intent(in) :: pawang 216 type(pawfgr_type),intent(in) :: pawfgr 217 type(pseudopotential_type),intent(in) :: psps 218 type(results_gs_type),intent(in) :: results_gs 219 type(wvl_denspot_type), intent(in) :: wvl_den 220 !arrays 221 integer,intent(in) :: atindx1(natom),dimcprj(natom*usecprj) 222 integer,intent(in) :: kg(3,mpw*mkmem),nattyp(ntypat),ngfft(18),npwarr(nkpt) 223 real(dp),intent(in) :: dmatpawu(:,:,:,:),eigen(mband*nkpt*nsppol) 224 real(dp),intent(in) :: gmet(3,3),gprimd(3,3) 225 real(dp),intent(in) :: intgres(:,:) ! (nspden,natom) if constrainedDFT otherwise (nspden,0) 226 real(dp),intent(in) :: occ(mband*nkpt*nsppol) 227 real(dp),intent(in) :: rprimd(3,3),vhartr(nfft),xccc3d(n3xccc) 228 real(dp),intent(in) :: vpsp(nfft) 229 real(dp),intent(inout) :: cg(2,mcg) 230 real(dp),intent(inout) :: nhat(nfft,nspden*psps%usepaw) 231 real(dp),intent(inout),target :: rhor(nfft,nspden),vtrial(nfft,nspden) 232 real(dp),intent(inout) :: vxc(nfft,nspden),xred(3,natom) 233 real(dp),pointer :: elfr(:,:),grhor(:,:,:),lrhor(:,:),taur(:,:) 234 type(pawcprj_type),intent(inout) :: cprj(natom,mcprj*usecprj) 235 type(paw_an_type),intent(inout) :: paw_an(my_natom*psps%usepaw) 236 type(pawfgrtab_type),intent(in) :: pawfgrtab(my_natom*psps%usepaw) 237 type(paw_ij_type),intent(inout) :: paw_ij(my_natom*psps%usepaw) 238 type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw) 239 type(pawrhoij_type),target,intent(inout) :: pawrhoij(my_natom*psps%usepaw) 240 type(pawtab_type),intent(inout) :: pawtab(ntypat*psps%usepaw) 241 242 !Local variables------------------------------- 243 !scalars 244 integer,parameter :: master=0,cplex1=1,fform_den=52,rdwr2=2,rdwrpaw0=0 245 integer :: bantot,fform,collect,timrev 246 integer :: accessfil,coordn 247 integer :: ii,ierr,ifft,ikpt,ispden,isppol,itypat 248 integer :: me_fft,n1,n2,n3 249 integer :: ifgd, iatom, iatom_tot,nradint 250 integer :: me,my_natom_tmp 251 integer :: occopt 252 integer :: prtnabla 253 integer :: pawprtden 254 integer :: iband,nocc,spacecomm,comm_fft,tmp_unt,nfft_tot 255 integer :: my_comm_atom 256 integer :: opt_imagonly 257 integer :: indsym(4,dtset%nsym,dtset%natom) 258 #ifdef HAVE_NETCDF 259 integer :: ncid 260 #endif 261 real(dp) :: norm,occ_norm,unocc_norm 262 real(dp) :: rate_dum,rate_dum2 263 real(dp) :: yp1, ypn, dr 264 character(len=500) :: msg 265 character(len=fnlen) :: fname 266 !arrays 267 integer, allocatable :: isort(:) 268 integer, pointer :: my_atmtab(:) 269 real(dp) :: tsec(2),nt_ntone_norm(nspden),rhomag(2,nspden) 270 real(dp),allocatable :: eigen2(:) 271 real(dp),allocatable :: elfr_down(:,:),elfr_up(:,:),intgden(:,:) 272 real(dp),allocatable :: rhor_paw(:,:),rhor_paw_core(:,:),rhor_paw_val(:,:),vpaw(:,:),vwork(:,:) 273 real(dp),allocatable :: rhor_n_one(:,:),rhor_nt_one(:,:),ps_norms(:,:,:) 274 real(dp), allocatable :: doccde(:) 275 real(dp), allocatable :: vh1spl(:) 276 real(dp), allocatable :: vh1_interp(:) 277 real(dp), allocatable :: vh1_integ(:) 278 real(dp), allocatable :: vh1_corrector(:) 279 real(dp), allocatable :: radii(:) 280 real(dp), ABI_CONTIGUOUS pointer :: rho_ptr(:,:) 281 type(pawrhoij_type) :: pawrhoij_dum(1) 282 !type(pawrhoij_type) :: pawrhoij_dum(0) 283 type(pawrhoij_type),pointer :: pawrhoij_all(:) 284 logical :: remove_inv 285 logical :: paral_atom, paral_fft, my_atmtab_allocated 286 real(dp) :: e_zeeman 287 real(dp) :: dmatdum(0,0,0,0) 288 real(dp) :: e_fermie, e_fermih 289 type(oper_type) :: dft_occup 290 type(crystal_t) :: crystal 291 type(ebands_t) :: ebands 292 type(epjdos_t) :: dos 293 type(plowannier_type) :: wan 294 type(self_type) :: selfr 295 type(self_type) :: self 296 type(green_type) :: greenr 297 298 ! ************************************************************************* 299 300 DBG_ENTER("COLL") 301 302 call timab(1150,1,tsec) ! outscfcv 303 call timab(1151,1,tsec) ! outscfcv(preparation) 304 305 if ((usecprj==0.or.mcprj==0).and.psps%usepaw==1.and. & 306 (dtset%prtwant==2.or.dtset%prtwant==3.or.dtset%prtnabla>0.or.dtset%prtdos==3 & 307 .or.dtset%kssform==3.or.dtset%pawfatbnd>0.or.dtset%pawprtwf>0)) then 308 write (msg,'(5a)')& 309 & 'cprj datastructure must be allocated',ch10,& 310 & 'with options prtwant=2,3, prtnabla>0, prtdos>3, kssform==3, pawfatbnd>0, pawprtwf>0',ch10,& 311 & 'Action: change pawusecp input keyword.' 312 ABI_ERROR(msg) 313 end if 314 315 ! Parameters for MPI-FFT 316 n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3); nfft_tot = product(ngfft(1:3)) 317 comm_fft = mpi_enreg%comm_fft 318 me_fft = xmpi_comm_rank(comm_fft) 319 paral_fft = (mpi_enreg%paral_kgb==1) 320 321 spacecomm = mpi_enreg%comm_cell 322 me = xmpi_comm_rank(spacecomm) 323 324 paral_atom=(my_natom/=natom) 325 my_comm_atom = mpi_enreg%comm_atom 326 nullify(my_atmtab) 327 if (paral_atom) then 328 call get_my_atmtab(mpi_enreg%comm_atom, my_atmtab, my_atmtab_allocated, paral_atom,natom,my_natom_ref=my_natom) 329 else 330 ABI_MALLOC(my_atmtab, (natom)) 331 my_atmtab = (/ (iatom, iatom=1, natom) /) 332 my_atmtab_allocated = .true. 333 end if 334 335 ! Initialize two objects to facilitate the propagation of info. 336 ! These objects should used more frequently, actually they should 337 ! become basic objects used in abinit. 338 339 ! Crystalline structure. 340 remove_inv=.false. 341 if (dtset%nspden==4 .and. dtset%usedmft==1) remove_inv=.true. ! MG: why this? 342 343 timrev = 2; if (any(dtset%kptopt == [3, 4])) timrev= 1 344 call crystal_init(dtset%amu_orig(:,1),crystal,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 dtset%symrel,dtset%tnons,dtset%symafm) 348 349 ! Electron band energies. 350 bantot= dtset%mband*dtset%nkpt*dtset%nsppol 351 ABI_CALLOC(doccde, (bantot)) 352 call ebands_init(bantot,ebands,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 hdr%kptrlatt, hdr%nshiftk, hdr%shiftk) 357 358 ABI_FREE(doccde) 359 360 ebands%fermie = results_gs%energies%e_fermie 361 e_fermie = results_gs%energies%e_fermie 362 ebands%fermih = results_gs%energies%e_fermih 363 e_fermih = results_gs%energies%e_fermih 364 ebands%entropy = results_gs%energies%entropy 365 366 ! YAML output 367 if (me == master) then 368 call results_gs%yaml_write(ab_out, cryst=crystal, info="Summary of ground state results",& 369 & occopt=dtset%occopt, with_conv=(dtset%nstep > 0) ) 370 end if 371 372 call timab(1151,2,tsec) 373 374 !wannier interface 375 call timab(1152,1,tsec) 376 377 if (dtset%prtwant==2) then 378 379 call mlwfovlp(crystal, ebands, hdr, atindx1,cg,cprj,dtset,dtfil,eigen,gprimd,kg,& 380 & mband,mcg,mcprj,mgfftc,mkmem,mpi_enreg,mpw,natom,& 381 & nattyp,nfft,ngfft,nkpt,npwarr,nsppol,ntypat,occ,& 382 & pawang,pawrad,pawtab,prtvol,psps,rprimd,ucvol,xred) 383 384 else if (dtset%prtwant==3) then 385 386 ! Convert cg and eigen to GW quasiparticle wave functions and eigenvalues in mlwfovlp_qp 387 ABI_MALLOC(eigen2,(mband*nkpt*nsppol)) 388 eigen2=eigen 389 390 call mlwfovlp_qp(cg,cprj,dtset,dtfil,eigen2,mband,mcg,mcprj,mkmem,mpw,natom,& 391 & nkpt,npwarr,nspden,nsppol,ntypat,Hdr,pawtab,rprimd,MPI_enreg) 392 393 ! Call Wannier90 394 call mlwfovlp(crystal, ebands, hdr, atindx1,cg,cprj,dtset,dtfil,eigen2,gprimd,kg,& 395 & mband,mcg,mcprj,mgfftc,mkmem,mpi_enreg,mpw,natom,& 396 & nattyp,nfft,ngfft,nkpt,npwarr,nsppol,ntypat,occ,& 397 & pawang,pawrad,pawtab,prtvol,psps,rprimd,ucvol,xred) 398 399 ! this is the old implementation, risky due to unpredictable size effects 400 ! now eigen is not overwritten, one should use other ways to print the GW corrections 401 ! eigen=eigen2 402 ABI_FREE(eigen2) 403 end if !prtwant 404 405 call timab(1152,2,tsec) 406 call timab(1153,1,tsec) 407 408 occopt=dtset%occopt 409 410 prtnabla=dtset%prtnabla 411 pawprtden=dtset%prtden-1 412 413 spacecomm=mpi_enreg%comm_cell; me=xmpi_comm_rank(spacecomm) 414 comm_fft=mpi_enreg%comm_fft 415 paral_atom=(my_natom/=natom) 416 417 !Warnings : 418 !- core charge is excluded from the charge density; 419 !- the potential is the INPUT vtrial. 420 421 if (iwrite_fftdatar(mpi_enreg) .and. dtset%usewvl==0) then 422 423 ! output the density. 424 if (dtset%prtden/=0) then 425 if (dtset%positron/=1) rho_ptr => rhor 426 if (dtset%positron==1) rho_ptr => electronpositron%rhor_ep 427 call fftdatar_write("density",dtfil%fnameabo_app_den,dtset%iomode,hdr,& 428 crystal,ngfft,cplex1,nfft,nspden,rho_ptr,mpi_enreg,ebands=ebands) 429 430 if (dtset%positron/=0) then 431 if (dtset%positron/=1) rho_ptr => electronpositron%rhor_ep 432 if (dtset%positron==1) rho_ptr => rhor 433 fname = trim(dtfil%fnameabo_app_den)//'_POSITRON' 434 if (dtset%iomode == IO_MODE_ETSF) fname = strcat(fname, ".nc") 435 call fftdatar_write("positron_density",fname,dtset%iomode,hdr,& 436 crystal,ngfft,cplex1,nfft,nspden,rho_ptr,mpi_enreg,ebands=ebands) 437 end if 438 end if 439 440 else if (dtset%usewvl == 1 .and. dtset%prtden /= 0) then 441 !if iomode == 2 then set all outputs to netcdf format 442 !if iomode == 3 then set all outputs to ETSF format 443 accessfil = 0 444 if (dtset%iomode == IO_MODE_ETSF) accessfil = 3 445 if (dtset%iomode == IO_MODE_MPI) accessfil = 4 446 fform = fform_den 447 ! Write wavelet DEN. Note however that this should be delegate to separated Bigdft routines. 448 ! a lot of stuff written in outscf does not make sense if usewvl==0 449 call ioarr(accessfil,rhor,dtset,etotal,fform,dtfil%fnameabo_app_den, & 450 hdr,mpi_enreg,ngfft,cplex1,nfft,pawrhoij_dum,rdwr2,rdwrpaw0,wvl_den) 451 end if ! if master 452 453 !! MS - Printing of PAWDEN parallellised and several possible options included 454 !We output the total electron density in the PAW case 455 !this requires removing nhat from rhor and making PAW on-site corrections 456 if (pawprtden>0 .and. psps%usepaw==1) then 457 ! pawprtden 1 --> output PAW valence density 458 ! " 2 --> output PAW valence+core density 459 ! " 3 --> output core, valence and full atomic protodensity 460 ! " 4 --> options 1+3 461 ! " 5 --> options 2+3 462 ! " 6 --> output all individual PAW density contributions 463 if (pawprtden/=3) then ! calc PAW valence density 464 ABI_MALLOC(rhor_paw,(pawfgr%nfft,nspden)) 465 ABI_MALLOC(rhor_n_one,(pawfgr%nfft,nspden)) 466 ABI_MALLOC(rhor_nt_one,(pawfgr%nfft,nspden)) 467 ! If the communicator used for denfgr is kpt_comm, it is not compatible with paral_atom 468 if (mpi_enreg%paral_kgb==0.and.my_natom/=natom) then 469 my_natom_tmp=natom 470 ABI_MALLOC(pawrhoij_all,(natom)) 471 call pawrhoij_nullify(pawrhoij_all) 472 call pawrhoij_copy(pawrhoij,pawrhoij_all,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 473 & keep_cplex=.false.,keep_qphase=.false.,keep_itypat=.false.,keep_nspden=.false.) 474 else 475 my_natom_tmp=my_natom 476 pawrhoij_all => pawrhoij 477 end if 478 if (pawprtden/=6) then 479 call denfgr(atindx1,gmet,comm_fft,my_natom_tmp,natom,nattyp,ngfft,nhat,dtset%nspinor,nsppol,nspden,& 480 & ntypat,pawfgr,pawrad,pawrhoij_all,pawtab,prtvol,rhor,rhor_paw,rhor_n_one,& 481 & rhor_nt_one,rprimd,dtset%typat,ucvol,xred,& 482 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 483 else 484 call denfgr(atindx1,gmet,comm_fft,my_natom_tmp,natom,nattyp,ngfft,nhat,dtset%nspinor,nsppol,nspden,& 485 & ntypat,pawfgr,pawrad,pawrhoij_all,pawtab,prtvol,rhor,rhor_paw,rhor_n_one,& 486 & rhor_nt_one,rprimd,dtset%typat,ucvol,xred,& 487 & abs_n_tilde_nt_diff=nt_ntone_norm,znucl=dtset%znucl,& 488 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 489 end if 490 if (mpi_enreg%paral_kgb==0.and.my_natom/=natom) then 491 call pawrhoij_free(pawrhoij_all) 492 ABI_FREE(pawrhoij_all) 493 end if 494 495 if (prtvol>9) then ! Check normalisation 496 norm = SUM(rhor_paw(:,1))*ucvol/PRODUCT(pawfgr%ngfft(1:3)) 497 call xmpi_sum(norm,comm_fft,ierr) 498 write(msg,'(a,F8.4)') ' PAWDEN - NORM OF DENSITY: ',norm 499 call wrtout(std_out, msg) 500 end if 501 end if 502 503 if (pawprtden>1.AND.pawprtden<6) then ! We will need the core density 504 ABI_MALLOC(rhor_paw_core,(pawfgr%nfft,nspden)) 505 call read_atomden(mpi_enreg,natom,pawfgr%nfft,pawfgr%ngfft,nspden,ntypat,rhor_paw_core,& 506 & dtset%typat,rprimd,xred,prtvol,file_prefix='core ') 507 508 if (prtvol>9) then ! Check normalisation 509 norm = SUM(rhor_paw_core(:,1))*ucvol/PRODUCT(pawfgr%ngfft(1:3)) 510 call xmpi_sum(norm,comm_fft,ierr) 511 write(msg,'(a,F8.4)') ' ATMDEN - NORM OF CORE DENSITY: ', norm 512 call wrtout(std_out, msg) 513 end if 514 end if 515 516 if (pawprtden>2.AND.pawprtden<6) then ! We will need the valence protodensity 517 ABI_MALLOC(rhor_paw_val,(pawfgr%nfft,nspden)) 518 call read_atomden(mpi_enreg,natom,pawfgr%nfft,pawfgr%ngfft,nspden,ntypat,rhor_paw_val,& 519 & dtset%typat,rprimd,xred,prtvol,file_prefix='valence') 520 521 if (prtvol>9) then ! Check normalisation 522 norm = SUM(rhor_paw_val(:,1))*ucvol/PRODUCT(pawfgr%ngfft(1:3)) 523 call xmpi_sum(norm,comm_fft,ierr) 524 write(msg,'(a,F8.4)') ' ATMDEN - NORM OF VALENCE PROTODENSITY: ', norm 525 call wrtout(std_out, msg) 526 end if 527 end if 528 529 if (iwrite_fftdatar(mpi_enreg)) then 530 if (pawprtden/=3) then 531 if (pawprtden==2.or.pawprtden==5) rhor_paw = rhor_paw + rhor_paw_core 532 ! PAWDEN 533 call fftdatar_write("pawrhor",dtfil%fnameabo_app_pawden,dtset%iomode,hdr,& 534 crystal,ngfft,cplex1,nfft,nspden,rhor_paw,mpi_enreg,ebands=ebands) 535 end if 536 537 if (pawprtden>2.AND.pawprtden<6) then 538 ! ATMDEN_CORE 539 call fftdatar_write("pawrhor_core",dtfil%fnameabo_app_atmden_core,dtset%iomode,hdr,& 540 crystal,ngfft,cplex1,nfft,nspden,rhor_paw_core,mpi_enreg,ebands=ebands) 541 542 ! valence protodensity. ATMDEN_VAL 543 call fftdatar_write("pawrhor_val",dtfil%fnameabo_app_atmden_val,dtset%iomode,hdr,& 544 crystal,ngfft,cplex1,nfft,nspden,rhor_paw_val,mpi_enreg,ebands=ebands) 545 546 ! full protodensity. ATMDEN_FULL 547 rhor_paw_val = rhor_paw_val + rhor_paw_core 548 call fftdatar_write("pawrhor_full",dtfil%fnameabo_app_atmden_full,dtset%iomode,hdr,& 549 crystal,ngfft,cplex1,nfft,nspden,rhor_paw_val,mpi_enreg,ebands=ebands) 550 end if 551 552 if (pawprtden==6) then ! Print all individual contributions to the density 553 ! N_TILDE - N_HAT 554 ! Use rhor_paw_val as temporary array 555 if (.not.allocated(rhor_paw_val)) then 556 ABI_MALLOC(rhor_paw_val,(pawfgr%nfft,nspden)) 557 end if 558 rhor_paw_val = rhor - nhat 559 560 call fftdatar_write("pawrhor_ntilde_minus_nhat",dtfil%fnameabo_app_n_tilde,dtset%iomode,hdr,& 561 crystal,ngfft,cplex1,nfft,nspden,rhor_paw_val,mpi_enreg,ebands=ebands) 562 563 ! N_ONSITE 564 call fftdatar_write("pawrhor_n_one",dtfil%fnameabo_app_n_one,dtset%iomode,hdr,& 565 crystal,ngfft,cplex1,nfft,nspden,rhor_n_one,mpi_enreg,ebands=ebands) 566 567 ! N_TILDE_ONSITE 568 call fftdatar_write("pawrhor_nt_one",dtfil%fnameabo_app_nt_one,dtset%iomode,hdr,& 569 crystal,ngfft,cplex1,nfft,nspden,rhor_nt_one,mpi_enreg,ebands=ebands) 570 571 end if ! All indivdual density cont. 572 end if ! if master 573 574 ABI_SFREE(rhor_paw) 575 ABI_SFREE(rhor_paw_core) 576 ABI_SFREE(rhor_paw_val) 577 ABI_SFREE(rhor_n_one) 578 ABI_SFREE(rhor_nt_one) 579 580 end if ! if paw+pawprtden 581 582 call timab(1153,2,tsec) 583 call timab(1154,1,tsec) 584 585 ! Output of the GSR file (except when we are inside mover) 586 #ifdef HAVE_NETCDF 587 ! Temporarily disable for CRAY 588 #ifndef FC_CRAY 589 if (me == master .and. dtset%prtgsr == 1 .and. dtset%usewvl == 0) then 590 !.and. (dtset%ionmov /= 0 .or. dtset%optcell /= 0)) then 591 fname = strcat(dtfil%filnam_ds(4), "_GSR.nc") 592 593 ! Write crystal and band structure energies. 594 !call timab(1190,1,tsec) 595 NCF_CHECK(nctk_open_create(ncid, fname, xmpi_comm_self)) 596 !call timab(1190,2,tsec) 597 598 !call timab(1191,1,tsec) 599 NCF_CHECK(hdr%ncwrite(ncid, fform_den, spinat=dtset%spinat, nc_define=.True.)) 600 !call timab(1191,2,tsec) 601 602 !call timab(1192,1,tsec) 603 NCF_CHECK(crystal%ncwrite(ncid)) 604 !call timab(1192,2,tsec) 605 606 !call timab(1193,1,tsec) 607 NCF_CHECK(ebands_ncwrite(ebands, ncid)) 608 !call timab(1193,2,tsec) 609 610 ! Add energy, forces, stresses 611 !call timab(1194,1,tsec) 612 NCF_CHECK(results_gs_ncwrite(results_gs, ncid, dtset%ecut, dtset%pawecutdg)) 613 !call timab(1194,2,tsec) 614 615 !call timab(1195,1,tsec) 616 NCF_CHECK(nf90_close(ncid)) 617 !call timab(1195,2,tsec) 618 end if 619 #endif 620 #endif 621 622 call timab(1154,2,tsec) 623 call timab(1155,1,tsec) 624 625 ! Output of VCLMB file 626 ! The PAW correction has to be computed here (all processors contribute) 627 if (psps%usepaw > 0 .AND. dtset%prtvclmb>0) then 628 nradint = 1000 ! radial integration grid density 629 ABI_MALLOC(vpaw,(nfft,nspden)) 630 vpaw(:,:)=zero 631 if (me == master .and. my_natom > 0) then 632 if (paw_an(1)%cplex > 1) then 633 ABI_WARNING('cplex = 2 : complex hartree potential in PAW spheres. This is not coded yet. Imag part ignored') 634 end if 635 end if 636 637 do ispden=1,nspden 638 ! for points inside spheres, replace with full AE hartree potential. 639 ! In principle the correction could be more subtle (not spherical) 640 do iatom=1,my_natom 641 iatom_tot=iatom;if (paral_atom) iatom_tot=mpi_enreg%my_atmtab(iatom) 642 itypat=dtset%typat(iatom_tot) 643 644 ABI_MALLOC(vh1spl,(paw_an(iatom)%mesh_size)) 645 ABI_MALLOC(vh1_corrector,(paw_an(iatom)%mesh_size)) 646 ABI_MALLOC(vh1_interp,(pawfgrtab(iatom)%nfgd)) 647 ABI_MALLOC(radii,(pawfgrtab(iatom)%nfgd)) 648 ABI_MALLOC(isort,(pawfgrtab(iatom)%nfgd)) 649 ! vh1 vht1 contain the spherical first moments of the Hartree potentials, so re-divide by Y_00 = sqrt(four_pi) 650 vh1_corrector(:) = (paw_an(iatom)%vh1(:,1,ispden)-paw_an(iatom)%vht1(:,1,ispden)) / sqrt(four_pi) 651 652 ! get end point derivatives 653 call bound_deriv(vh1_corrector, pawrad(itypat), pawrad(itypat)%mesh_size, yp1, ypn) 654 ! spline the vh1 function 655 ! NB for second argument of vh1: only first moment lm_size appears to be used 656 ! NB2: vh1 can in principle be complex - not sure what to do with the imaginary part. Ignored for now. 657 call spline(pawrad(itypat)%rad, vh1_corrector, paw_an(iatom)%mesh_size, yp1, ypn, vh1spl) 658 659 do ifgd = 1, pawfgrtab(iatom)%nfgd 660 ! get radii for this point 661 isort(ifgd) = ifgd 662 radii(ifgd) = sqrt(sum(pawfgrtab(iatom)%rfgd(:,ifgd)**2)) 663 end do 664 665 if (pawfgrtab(iatom)%nfgd/=0) then 666 ! spline interpolate the vh1 value for current radii 667 call sort_dp(pawfgrtab(iatom)%nfgd, radii, isort, tol12) 668 call splint(pawrad(itypat)%mesh_size, pawrad(itypat)%rad, & 669 & vh1_corrector, vh1spl, pawfgrtab(iatom)%nfgd, radii, vh1_interp, ierr) 670 end if 671 672 norm=SUM(vh1_interp)*ucvol/PRODUCT(ngfft(1:3)) 673 call xmpi_sum(norm,comm_fft,ierr) 674 write(msg,'(a,i6,a,E20.10)') ' sum of Hartree correction term on fft grid of atom : ', iatom, ' = ', norm 675 call wrtout(std_out, msg) 676 677 if (pawfgrtab(iatom)%nfgd/=0) then 678 vpaw(pawfgrtab(iatom)%ifftsph(isort(1:pawfgrtab(iatom)%nfgd)),ispden) = & 679 & vpaw(pawfgrtab(iatom)%ifftsph(isort(1:pawfgrtab(iatom)%nfgd)),ispden) + & 680 & vh1_interp(1:pawfgrtab(iatom)%nfgd) 681 end if 682 683 ! get integral of correction term in whole sphere 684 ABI_FREE(radii) 685 ABI_FREE(vh1_interp) 686 687 ABI_MALLOC(radii,(nradint)) 688 ABI_MALLOC(vh1_interp,(nradint)) 689 690 ABI_MALLOC(vh1_integ,(nradint)) 691 dr = pawrad(itypat)%rad(paw_an(iatom)%mesh_size) / dble(nradint) 692 do ifgd = 1, nradint 693 radii(ifgd) = dble(ifgd-1)*dr 694 end do 695 696 ! spline interpolate the vh1 value for current radii 697 call splint(pawrad(itypat)%mesh_size, pawrad(itypat)%rad, & 698 & vh1_corrector, vh1spl, nradint, radii, vh1_interp, ierr) 699 700 do ifgd = 1, nradint 701 vh1_interp(ifgd) = vh1_interp(ifgd)*radii(ifgd)**2 702 end do 703 704 call simpson_int(nradint, dr, vh1_interp, vh1_integ) 705 write(msg,'(a,i6,a,E20.10)') ' integral of Hartree correction term in sphere of atom: ', iatom, & 706 & ' = ', vh1_integ(nradint)*four*pi 707 call wrtout(std_out, msg) 708 709 ABI_FREE(vh1spl) 710 ABI_FREE(vh1_corrector) 711 ABI_FREE(vh1_interp) 712 ABI_FREE(vh1_integ) 713 ABI_FREE(radii) 714 ABI_FREE(isort) 715 end do ! iatom 716 end do !ispden 717 call xmpi_sum_master(vpaw,master,mpi_enreg%comm_atom,ierr) 718 if (.not.iwrite_fftdatar(mpi_enreg)) then 719 ABI_FREE(vpaw) 720 end if 721 end if ! if paw - add all electron vhartree in spheres 722 723 call timab(1155,2,tsec) 724 725 if (iwrite_fftdatar(mpi_enreg)) then 726 727 call timab(1156,1,tsec) 728 729 ! output the electron localization function ELF 730 if (dtset%prtelf/=0) then 731 call fftdatar_write("elfr",dtfil%fnameabo_app_elf,dtset%iomode,hdr,& 732 crystal,ngfft,cplex1,nfft,nspden,elfr,mpi_enreg,ebands=ebands) 733 734 if (nspden==2)then 735 ABI_MALLOC(elfr_up,(nfft,nspden)) 736 elfr_up(:,:) = zero 737 do ifft=1,nfft 738 elfr_up(ifft,1) = elfr(ifft,2) 739 end do 740 ! ELF_UP 741 call fftdatar_write("elfr_up",dtfil%fnameabo_app_elf_up,dtset%iomode,hdr,& 742 crystal,ngfft,cplex1,nfft,nspden,elfr_up,mpi_enreg,ebands=ebands) 743 744 ABI_MALLOC(elfr_down,(nfft,nspden)) 745 elfr_down(:,:) = zero 746 do ifft=1,nfft 747 elfr_down(ifft,1) = elfr(ifft,3) 748 end do 749 ! ELF_DOWN' 750 call fftdatar_write("elfr_down",dtfil%fnameabo_app_elf_down,dtset%iomode,hdr,& 751 crystal,ngfft,cplex1,nfft,nspden,elfr_down,mpi_enreg,ebands=ebands) 752 753 ABI_FREE(elfr_up) 754 ABI_FREE(elfr_down) 755 end if 756 end if 757 758 call timab(1156,2,tsec) 759 call timab(1157,1,tsec) 760 761 ! We output the gradient of density 762 if (dtset%prtgden/=0) then 763 764 call fftdatar_write("grhor_1",dtfil%fnameabo_app_gden1,dtset%iomode,hdr,& 765 crystal,ngfft,cplex1,nfft,nspden,grhor(:,:,1),mpi_enreg,ebands=ebands) 766 767 call fftdatar_write("grhor_2",dtfil%fnameabo_app_gden2,dtset%iomode,hdr,& 768 crystal,ngfft,cplex1,nfft,nspden,grhor(:,:,2),mpi_enreg,ebands=ebands) 769 770 call fftdatar_write("grhor_3",dtfil%fnameabo_app_gden3,dtset%iomode,hdr,& 771 crystal,ngfft,cplex1,nfft,nspden,grhor(:,:,3),mpi_enreg,ebands=ebands) 772 end if 773 774 call timab(1157,2,tsec) 775 call timab(1158,1,tsec) 776 777 ! We output the total kinetic energy density KDEN 778 if (dtset%prtkden/=0) then 779 call fftdatar_write("kinedr",dtfil%fnameabo_app_kden,dtset%iomode,hdr,& 780 crystal,ngfft,cplex1,nfft,nspden,taur,mpi_enreg,ebands=ebands) 781 end if 782 783 call timab(1158,2,tsec) 784 call timab(1159,1,tsec) 785 786 787 ! We output the Laplacian of density 788 if (dtset%prtlden/=0) then 789 call fftdatar_write("laprhor",dtfil%fnameabo_app_lden,dtset%iomode,hdr,& 790 crystal,ngfft,cplex1,nfft,nspden,lrhor,mpi_enreg,ebands=ebands) 791 end if 792 793 call timab(1159,2,tsec) 794 call timab(1160,1,tsec) 795 796 ! POT 797 if (dtset%prtpot>0) then 798 call fftdatar_write("vtrial",dtfil%fnameabo_app_pot,dtset%iomode,hdr,& 799 crystal,ngfft,cplex1,nfft,nspden,vtrial,mpi_enreg,ebands=ebands) 800 end if 801 802 ! EIG 803 #if defined HAVE_NETCDF 804 if (dtset%prteig==2 .and. me == master) then 805 fname=trim(dtfil%fnameabo_app_eig)//'.nc' 806 call write_eig(eigen,e_fermie,fname,dtset%kptns,dtset%mband,dtset%nband,dtset%nkpt,dtset%nsppol,& 807 & results_gs%shiftfactor_extfpmd) ! Optional arguments 808 end if 809 #endif 810 811 call timab(1160,2,tsec) 812 call timab(1161,1,tsec) 813 814 if (dtset%prtgeo>0) then 815 coordn=dtset%prtgeo 816 call bonds_lgth_angles(coordn,dtfil%fnameabo_app_geo,natom,psps%ntypat,& 817 rprimd,dtset%typat,xred,dtset%znucl) 818 end if 819 820 if (dtset%prtcif > 0) then 821 call prt_cif(dtset%brvltt, dtfil%fnameabo_app_cif, natom, dtset%nsym, dtset%ntypat, rprimd, & 822 dtset%spgaxor, dtset%spgroup, dtset%spgorig, dtset%symrel, dtset%tnons, dtset%typat, xred, dtset%znucl) 823 end if 824 825 call timab(1161,2,tsec) 826 call timab(1162,1,tsec) 827 828 ! STM 829 if (dtset%prtstm/=0) then 830 call fftdatar_write("stm",dtfil%fnameabo_app_stm,dtset%iomode,hdr,& 831 crystal,ngfft,cplex1,nfft,nspden,rhor,mpi_enreg,ebands=ebands) 832 end if 833 834 call timab(1162,2,tsec) 835 call timab(1163,1,tsec) 836 837 if (dtset%prt1dm>0) then 838 call out1dm(dtfil%fnameabo_app_1dm,mpi_enreg,natom,nfft,ngfft,nspden,psps%ntypat,& 839 & rhor,rprimd,dtset%typat,ucvol,vtrial,xred,dtset%znucl) 840 end if 841 842 call timab(1163,2,tsec) 843 call timab(1164,1,tsec) 844 845 ! VHA 846 if (dtset%prtvha>0) then 847 ABI_MALLOC(vwork,(nfft,nspden)) 848 do ispden=1,nspden 849 vwork(:,ispden)=vhartr(:) 850 end do 851 852 call fftdatar_write("vhartree",dtfil%fnameabo_app_vha,dtset%iomode,hdr,& 853 crystal,ngfft,cplex1,nfft,nspden,vwork,mpi_enreg,ebands=ebands) 854 855 ABI_FREE(vwork) 856 end if 857 858 ! VPSP 859 if (dtset%prtvpsp>0) then 860 ABI_MALLOC(vwork,(nfft,nspden)) 861 do ispden=1,nspden 862 vwork(:,ispden)=vpsp(:) 863 end do 864 865 call fftdatar_write("vpsp",dtfil%fnameabo_app_vpsp,dtset%iomode,hdr,& 866 crystal,ngfft,cplex1,nfft,nspden,vwork,mpi_enreg,ebands=ebands) 867 868 ABI_FREE(vwork) 869 end if 870 871 ! VCouLoMB 872 if (dtset%prtvclmb>0) then 873 874 ABI_MALLOC(vwork,(nfft,nspden)) 875 do ispden=1,nspden 876 vwork(:,ispden)=vpsp(:)+vhartr(:) 877 end do 878 if (psps%usepaw==1) then 879 do ispden=1,nspden 880 vwork(:,ispden)=vwork(:,ispden)+vpaw(:,ispden) 881 end do 882 ABI_FREE(vpaw) 883 end if 884 885 call fftdatar_write("vhartree_vloc",dtfil%fnameabo_app_vclmb,dtset%iomode,hdr,& 886 & crystal,ngfft,cplex1,nfft,nspden,vwork,mpi_enreg,ebands=ebands) 887 888 !TODO: find out why this combination of calls with fftdatar_write then out1dm fails on buda with 4 mpi-fft procs (np_spkpt 1). 889 ! For the moment comment it out. Only DS2 of mpiio test 27 fails 890 ! call out1dm(dtfil%fnameabo_app_vclmb_1dm,mpi_enreg,natom,nfft,ngfft,nspden,psps%ntypat,& 891 !& rhor,rprimd,dtset%typat,ucvol,vwork,xred,dtset%znucl) 892 893 ! 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 894 ! where E is energy of electron, E0 rest mass, lambda the relativistic wavelength 895 ! values of CE at 200 300 and 1000 kV: 7.29e6 6.53e6 5.39e6 rad / V / m 896 ! vertical integral of vclmb * c / ngfft(3) / cross sectional area factor (= sin(gamma)) 897 ! * 0.5291772083e-10*27.2113834 to get to SI 898 ! * CE factor above 899 ! should be done for each plane perpendicular to the axes... 900 ABI_FREE(vwork) 901 end if ! prtvclmb 902 903 904 ! VHXC 905 if (dtset%prtvhxc>0) then 906 ABI_MALLOC(vwork,(nfft,nspden)) 907 do ispden=1,nspden 908 vwork(:,ispden)=vhartr(:)+vxc(:,ispden) 909 end do 910 911 call fftdatar_write("vhxc",dtfil%fnameabo_app_vhxc,dtset%iomode,hdr,& 912 & crystal,ngfft,cplex1,nfft,nspden,vwork,mpi_enreg,ebands=ebands) 913 914 ABI_FREE(vwork) 915 end if 916 917 ! VXC 918 if (dtset%prtvxc>0) then 919 call fftdatar_write("exchange_correlation_potential",dtfil%fnameabo_app_vxc,dtset%iomode,hdr,& 920 crystal,ngfft,cplex1,nfft,nspden,vxc,mpi_enreg,ebands=ebands) 921 end if 922 923 call timab(1164,2,tsec) 924 925 end if ! if iwrite_fftdatar 926 927 call timab(1165,1,tsec) 928 929 !Generate DOS using the tetrahedron method or using Gaussians 930 !FIXME: Should centralize all calculations of DOS here in outscfcv 931 if (dtset%prtdos>=2.or.dtset%pawfatbnd>0) then 932 dos = epjdos_new(dtset, psps, pawtab) 933 934 if (dos%partial_dos_flag>=1 .or. dos%fatbands_flag==1)then 935 ! Generate fractions for partial DOSs if needed partial_dos 1,2,3,4 give different decompositions 936 collect = 1 !; if (psps%usepaw==1 .and. dos%partial_dos_flag /= 2) collect = 0 937 if ((psps%usepaw==0.or.dtset%pawprtdos/=2) .and. dos%partial_dos_flag>=1) then 938 call partial_dos_fractions(dos,crystal,dtset,eigen,occ,npwarr,kg,cg,mcg,collect,mpi_enreg) 939 end if 940 941 if (psps%usepaw==1 .and. dos%partial_dos_flag /= 2) then 942 ! TODO: update partial_dos_fractions_paw for extra atoms - no PAW contribution normally, but check bounds and so on. 943 call partial_dos_fractions_paw(dos,cprj,dimcprj,dtset,mcprj,mkmem,mpi_enreg,pawrad,pawtab) 944 end if 945 946 else 947 dos%fractions(:,:,:,1)=one 948 end if 949 950 ! Here, print out fatbands for the k-points given in file appended _FATBANDS 951 if (me == master .and. dtset%pawfatbnd>0 .and. dos%fatbands_flag==1) then 952 call prtfatbands(dos,dtset,ebands,dtfil%fnameabo_app_fatbands,dtset%pawfatbnd,pawtab) 953 end if 954 955 ! Here, computation and output of DOS and partial DOS _DOS 956 if (dos%fatbands_flag == 0 .and. dos%prtdos /= 4) then 957 call dos_calcnwrite(dos,dtset,crystal,ebands,dtfil%fnameabo_app_dos,spacecomm) 958 end if 959 960 #ifdef HAVE_NETCDF 961 ! Write netcdf file with dos% results. 962 if (me == master) then 963 fname = trim(dtfil%filnam_ds(4))//'_FATBANDS.nc' 964 NCF_CHECK(nctk_open_create(ncid, fname, xmpi_comm_self)) 965 call fatbands_ncwrite(dos, crystal, ebands, hdr, dtset, psps, pawtab, ncid) 966 NCF_CHECK(nf90_close(ncid)) 967 end if 968 #endif 969 970 call dos%free() 971 end if ! prtdos > 1 972 973 call timab(1165,2,tsec) 974 call timab(1166,1,tsec) 975 976 !Output of integrated density inside atomic spheres 977 if ((dtset%prtdensph==1.and.dtset%usewvl==0) .or. sum(abs(dtset%zeemanfield)) > tol10) then 978 ABI_MALLOC(intgden,(nspden,natom)) 979 call calcdenmagsph(mpi_enreg,natom,nfft,ngfft,nspden,& 980 & ntypat,dtset%ratsm,dtset%ratsph,rhor,rprimd,dtset%typat,xred,1,cplex1,intgden=intgden,rhomag=rhomag) 981 ! for rhomag: 982 ! in collinear case component 1 is total density and 2 is _magnetization_ up-down 983 ! in non collinear case component 1 is total density, and 2:4 are the magnetization vector 984 985 if (dtset%prtdensph==1.and.dtset%usewvl==0) then 986 if(all(dtset%constraint_kind(:)==0))then 987 call prtdenmagsph(cplex1,intgden,natom,nspden,ntypat,ab_out,1,dtset%ratsm,dtset%ratsph,rhomag,dtset%typat) 988 call prtdenmagsph(cplex1,intgden,natom,nspden,ntypat,std_out,1,dtset%ratsm,dtset%ratsph,rhomag,dtset%typat) 989 else 990 call prtdenmagsph(cplex1,intgden,natom,nspden,ntypat,ab_out,1,dtset%ratsm,dtset%ratsph,rhomag,dtset%typat,dtset%ziontypat) 991 call prtdenmagsph(cplex1,intgden,natom,nspden,ntypat,std_out,1,dtset%ratsm,dtset%ratsph,rhomag,dtset%typat,dtset%ziontypat) 992 endif 993 if(any(dtset%constraint_kind(:)/=0))then 994 call prtdenmagsph(cplex1,intgres,natom,nspden,ntypat,ab_out,21,dtset%ratsm,dtset%ratsph,rhomag,dtset%typat) 995 call prtdenmagsph(cplex1,intgres,natom,nspden,ntypat,std_out,21,dtset%ratsm,dtset%ratsph,rhomag,dtset%typat) 996 endif 997 end if 998 999 !!!!!!!!!!!!!!!!!!!!!!!!if prt_lorbmag value is equal 1 and the calculations are noncollinear then the local orbital magnetic moments are calculated 1000 if (dtset%prt_lorbmag==1) then 1001 1002 if ((dtset%nspinor .ne. 2) .and. (dtset%nsppol .ne.4)) then 1003 write (msg,'(a)')" " 1004 call wrtout([std_out, ab_out], msg) 1005 write (msg,'(a)')"WARNING*" 1006 call wrtout([std_out, ab_out], msg) 1007 write (msg,'(a)')"prt_lorbmag=1, To calcualte orbital magnetisation, calculations need to be noncollinear" 1008 call wrtout([std_out, ab_out], msg) 1009 else 1010 if (dtset%usepawu .ne. 0)then 1011 call loc_orbmom_cal(1,0,dmatdum,0,0,indsym,my_natom,dtset%natom,dtset%natpawu,& 1012 & dtset%nspinor,dtset%nsppol,dtset%nsym,dtset%ntypat,paw_ij,pawang,pawrad,dtset%pawprtvol,& 1013 & pawrhoij,pawtab,dtset%spinat,dtset%symafm,dtset%typat,0,dtset%usepawu,& 1014 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 1015 else 1016 write (msg,'(a)')" " 1017 call wrtout([std_out, ab_out], msg) 1018 write (msg,'(a)')"WARNING*" 1019 call wrtout([std_out, ab_out], msg) 1020 write (msg,'(a)')"prt_lorbmag=1, To calcualte orbital magnetisation LDA+U calculations should be activated" 1021 call wrtout([std_out, ab_out], msg) 1022 end if 1023 endif 1024 end if 1025 1026 if (sum(abs(dtset%zeemanfield)) > tol10) then 1027 if(nspden==2)then 1028 e_zeeman = -half*rhomag(1,2)*dtset%zeemanfield(3) 1029 write (msg, "(a,E20.10,a)") " Collinear magnetization ", rhomag(1,2), & 1030 & " (in # of spins, without 1/2 for magnetic moment) " 1031 call wrtout([std_out, ab_out], msg) 1032 else if(nspden==4)then 1033 e_zeeman = -half * (dtset%zeemanfield(1)*rhomag(1,2)& ! x 1034 & +dtset%zeemanfield(2)*rhomag(1,3)& ! y 1035 & +dtset%zeemanfield(3)*rhomag(1,4)) ! z 1036 write (msg, "(a,3E20.10,a)") " Magnetization vector ", rhomag(1,2:4), & 1037 & " (in # of spins, without 1/2 for magnetic moment) " 1038 call wrtout([std_out, ab_out], msg) 1039 end if 1040 !TODO: this quantity should also be calculated in rhotov, and stored in 1041 ! results_gs%energies%e_zeeman, but for the moment it comes out 0 1042 write (msg, "(a,E20.10,a)") " Zeeman energy -m.B = ", e_zeeman, " Ha" 1043 call wrtout([std_out, ab_out], msg) 1044 end if 1045 ABI_SFREE(intgden) 1046 end if 1047 1048 call timab(1166,2,tsec) 1049 call timab(1167,1,tsec) 1050 1051 if (dtset%magconon /= 0) then 1052 ! calculate final value of terms for magnetic constraint: "energy" term, lagrange multiplier term, and atomic contributions 1053 call mag_penalty_e(dtset%magconon,dtset%magcon_lambda,mpi_enreg,& 1054 & natom,nfft,ngfft,nspden,ntypat,dtset%ratsm,dtset%ratsph,rhor,rprimd,dtset%spinat,dtset%typat,xred) 1055 end if 1056 1057 call timab(1167,2,tsec) 1058 call timab(1168,1,tsec) 1059 1060 !If PAW, provide additional outputs 1061 if (psps%usepaw==1) then 1062 ! Output of compensation charge 1063 if (dtset%nstep>0.or.dtfil%ireadwf/=0) then 1064 write(msg, '(4a)' )ch10,' PAW TEST:',ch10,& 1065 & ' ==== Compensation charge inside spheres ============' 1066 if (compch_sph>-1.d4.and.compch_fft>-1.d4) & 1067 & write(msg, '(3a)' ) trim(msg),ch10,' The following values must be close to each other ...' 1068 if (compch_sph>-1.d4) write(msg, '(3a,f22.15)' ) trim(msg),ch10,& 1069 & ' Compensation charge over spherical meshes = ',compch_sph 1070 if (compch_fft>-1.d4) then 1071 if (pawfgr%usefinegrid==1) then 1072 write(msg, '(3a,f22.15)' ) trim(msg),ch10,& 1073 & ' Compensation charge over fine fft grid = ',compch_fft 1074 else 1075 write(msg, '(3a,f22.15)' ) trim(msg),ch10,& 1076 & ' Compensation charge over fft grid = ',compch_fft 1077 end if 1078 end if 1079 call wrtout([std_out, ab_out], msg) 1080 end if 1081 ! Output of pseudopotential strength Dij and augmentation occupancies Rhoij 1082 call pawprt(dtset,my_natom,paw_ij,pawrhoij,pawtab,& 1083 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 1084 & electronpositron=electronpositron) 1085 end if 1086 1087 call timab(1168,2,tsec) 1088 call timab(1169,1,tsec) 1089 1090 1091 !PAW + output for optical conductivity _OPT and _OPT2 1092 if (psps%usepaw==1.and.prtnabla>0) then 1093 if (prtnabla==1.or.prtnabla==2) then 1094 call optics_paw(atindx1,cg,cprj,dimcprj,dtfil,dtset,eigen,gprimd,hdr,kg,& 1095 & mband,mcg,mcprj,mkmem,mpi_enreg,mpsang,mpw,natom,nkpt,npwarr,nsppol,pawang,& 1096 & pawrad,pawrhoij,pawtab,psps%znuclpsp) 1097 end if 1098 if (prtnabla==2.or.prtnabla==3) then 1099 call optics_paw_core(atindx1,cprj,dimcprj,dtfil,dtset,eigen,psps%filpsp,hdr,& 1100 & mband,mcprj,mkmem,mpi_enreg,mpsang,natom,nkpt,nsppol,pawang,pawrad,pawrhoij,pawtab,& 1101 & psps%znuclpsp) 1102 end if 1103 end if 1104 if (prtnabla<0) then 1105 ! TODO: This routine is not tested but it's used in production. 1106 call optics_vloc(cg,dtfil,dtset,eigen,gprimd,hdr,kg,& 1107 & mband,mcg,mkmem,mpi_enreg,mpw,nkpt,npwarr,nsppol) 1108 end if 1109 1110 call timab(1169,2,tsec) 1111 call timab(1170,1,tsec) 1112 1113 !Optionally provide output for AE wavefunctions (only for PAW) 1114 if (psps%usepaw==1 .and. dtset%pawprtwf==1) then 1115 ABI_MALLOC(ps_norms,(nsppol,nkpt,mband)) 1116 1117 call pawmkaewf(dtset,crystal,ebands,my_natom,mpw,mband,mcg,mcprj,nkpt,mkmem,nsppol,Dtset%nband,& 1118 & Dtset%istwfk,npwarr,Dtset%kptns,Dtset%ngfftdg,kg,dimcprj,pawfgrtab,& 1119 & Pawrad,Pawtab,Hdr,Dtfil,cg,Cprj,& 1120 & MPI_enreg,ierr,pseudo_norms=ps_norms,set_k=dtset%pawprt_k,set_band=dtset%pawprt_b,& 1121 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1122 1123 if (dtset%pawprt_b==0) then 1124 fname = strcat(dtfil%filnam_ds(4), '_PAWSTAT') 1125 if (open_file(fname, msg,newunit=tmp_unt,status='unknown',form='formatted') /= 0) then 1126 ABI_ERROR(msg) 1127 end if 1128 write(tmp_unt,'(5a)') '# This file contains the statistics on the cancellation of',ch10,& 1129 & '# the onsite pseudo component of the all-electron wavefunction',ch10,& 1130 & '# with the plane wave part' 1131 ii = 0 1132 do isppol=1,nsppol 1133 write(tmp_unt,'(a,i0)') '# isppol = ',isppol 1134 do ikpt=1,nkpt 1135 write(tmp_unt,'(a,i0)') '# ikpt = ',ikpt 1136 write(tmp_unt,'(a)') '# band norm' 1137 occ_norm = zero; unocc_norm = zero; nocc = 0 1138 do iband=1,dtset%nband(ikpt + (isppol-1)*nkpt) 1139 ii = ii + 1 1140 write(tmp_unt,'(i8,ES16.6)') iband,ps_norms(isppol,ikpt,iband) 1141 if (abs(occ(ii)) <= tol16) then 1142 unocc_norm = unocc_norm + ps_norms(isppol,ikpt,iband) 1143 else 1144 occ_norm = occ_norm + ps_norms(isppol,ikpt,iband) 1145 nocc = nocc + 1 1146 end if 1147 end do 1148 if(mband/=nocc)then 1149 write(tmp_unt,'(2(a,ES16.6))') '# occ average: ',occ_norm/real(nocc),& 1150 & ' unocc average: ',unocc_norm/real(mband-nocc) 1151 else 1152 write(tmp_unt,'(2(a,ES16.6))') '# occ average: ',occ_norm/real(nocc) 1153 end if 1154 end do 1155 end do 1156 close(tmp_unt) 1157 end if 1158 ABI_FREE(ps_norms) 1159 end if 1160 1161 call timab(1170,2,tsec) 1162 call timab(1171,1,tsec) 1163 1164 if(dtset%plowan_compute>0 .and. dtset%plowan_compute<10) then 1165 write(msg,'(2a,i3)') ch10,& 1166 & ' ====================================================================================== ' 1167 call wrtout([std_out, ab_out], msg) 1168 write(msg,'(2a,i3)') ch10,& 1169 & ' == Start computation of Projected Local Orbitals Wannier functions == ',dtset%nbandkss 1170 call wrtout([std_out, ab_out], msg) 1171 1172 ! == compute psichi 1173 1174 call init_plowannier(dtset%plowan_bandf,dtset%plowan_bandi,dtset%plowan_compute,& 1175 & dtset%plowan_iatom,dtset%plowan_it,dtset%plowan_lcalc,dtset%plowan_natom,& 1176 & dtset%plowan_nbl,dtset%plowan_nt,dtset%plowan_projcalc,dtset%acell_orig,& 1177 & dtset%kptns,dtset%nimage,dtset%nkpt,dtset%nspinor,dtset%nsppol,dtset%wtk,dtset%dmft_t2g,wan) 1178 call compute_coeff_plowannier(crystal,cprj,dimcprj,dtset,eigen,e_fermie,& 1179 & mpi_enreg,occ,wan,pawtab,psps,usecprj,dtfil%unpaw,pawrad,dtfil) 1180 if (me==master) then 1181 call print_plowannier(wan) 1182 endif 1183 call destroy_plowannier(wan) 1184 end if 1185 1186 call timab(1171,2,tsec) 1187 call timab(1172,1,tsec) 1188 1189 !Optionally provide output for the GW part of ABINIT 1190 if (dtset%nbandkss/=0) then 1191 ! Use DMFT to compute wannier function for cRPA calculation. 1192 if(dtset%usedmft==1) then 1193 write(msg,'(2a,i3)') ch10,& 1194 & ' Warning: Psichi are renormalized in datafordmft because nbandkss is used',dtset%nbandkss 1195 call wrtout(std_out, msg) 1196 call init_dmft(dmatpawu,dtset,e_fermie,dtfil%fnameabo_app,& 1197 & dtfil%filnam_ds(3),dtset%nspinor,paw_dmft,pawtab,psps,dtset%typat) 1198 call print_dmft(paw_dmft,dtset%pawprtvol) 1199 1200 ! == compute psichi 1201 call init_oper(paw_dmft,dft_occup) 1202 1203 call datafordmft(crystal,cprj,dimcprj,dtset,eigen,e_fermie & 1204 & ,dft_occup,dtset%mband,dtset%mband,dtset%mkmem,mpi_enreg,& 1205 & dtset%nkpt,dtset%nspinor,dtset%nsppol,occ,& 1206 & paw_dmft,paw_ij,pawang,pawtab,psps,usecprj,dtfil%unpaw,dtset%nbandkss) 1207 1208 opt_imagonly=0 1209 if(paw_dmft%dmft_solv>=5) opt_imagonly=1 1210 1211 1212 ! Compute k-resolved spectral function in DMFT. 1213 if(dtset%dmft_kspectralfunc==1) then 1214 ! Initialize self on real axis 1215 call initialize_self(selfr,paw_dmft,wtype='real') 1216 1217 ! Initialize self on imag axis 1218 call initialize_self(self,paw_dmft) 1219 1220 ! Initialize green on real axis 1221 call init_green(greenr,paw_dmft,opt_oper_ksloc=3,wtype='real') 1222 1223 ! Read self energy in imag. Matsubara freq (for double counting 1224 ! and limit at high frequency) 1225 call rw_self(self,paw_dmft,prtopt=5,opt_rw=1,opt_stop=1) 1226 1227 ! Read self energy on real axis obtained from Maxent 1228 call rw_self(selfr,paw_dmft,prtopt=5,opt_rw=1,opt_imagonly=opt_imagonly, & 1229 & opt_selflimit=self%oper(self%nw)%matlu,opt_hdc=self%hdc%matlu,pawang=pawang,cryst_struc=crystal) 1230 1231 ! Check: from self on real axis, recompute self on Imaginary axis. 1232 call selfreal2imag_self(selfr,self,paw_dmft%filapp) 1233 1234 ! paw_dmft%fermie=hdr%fermie ! for tests 1235 write(std_out,*) " Fermi level is",paw_dmft%fermie 1236 1237 ! selfr does not have any double couting in self%hdc 1238 ! hdc from self%hdc has been put in real part of self in rw_self. 1239 ! For the DFT BS: use opt_self=0 and fermie=fermie_dft 1240 1241 ! Compute green function on real axis 1242 call compute_green(crystal,greenr,paw_dmft,pawang,1,selfr,& 1243 & opt_self=1,opt_nonxsum=0) 1244 1245 !write(6,*) "compute green done" 1246 if(me==master) then 1247 if(dtset%kptopt<0) then 1248 ! k-resolved Spectral function 1249 call print_green("from_realaxisself",greenr,5,paw_dmft,& 1250 & pawprtvol=3,opt_wt=1) 1251 else 1252 ! DOS Calculation 1253 call print_green("from_realaxisself",greenr,4,paw_dmft,& 1254 & pawprtvol=3,opt_wt=1) 1255 endif 1256 !write(6,*) "print green done" 1257 endif 1258 1259 call destroy_green(greenr) 1260 call destroy_self(selfr) 1261 call destroy_self(self) 1262 endif 1263 call destroy_dmft(paw_dmft) 1264 call destroy_oper(dft_occup) 1265 end if 1266 1267 call outkss(crystal,dtfil,dtset,ecut,gmet,gprimd,hdr,& 1268 & dtset%kssform,mband,mcg,mcprj,mgfftc,mkmem,mpi_enreg,mpsang,mpw,natom,natom,& 1269 & nfft,nkpt,npwarr,nspden,nsppol,nsym,psps%ntypat,occ,pawtab,pawfgr,paw_ij,& 1270 & prtvol,psps,rprimd,vtrial,xred,cg,usecprj,cprj,eigen,ierr) 1271 if (ierr/=0) then 1272 ABI_WARNING("outkss returned a non zero status error, check log") 1273 end if 1274 end if 1275 1276 call timab(1172,2,tsec) ! outscfcv(gw) 1277 1278 if (electronpositron_calctype(electronpositron)/=0) then 1279 1280 ! Optionally provide output for positron life time calculation 1281 call timab(1173,1,tsec) 1282 call poslifetime(dtset,electronpositron,gprimd,my_natom,& 1283 & mpi_enreg,n3xccc,nfft,ngfft,nhat,1,pawang,& 1284 & pawrad,pawrhoij,pawtab,rate_dum,rate_dum2,& 1285 & rhor,ucvol,xccc3d) 1286 call timab(1173,2,tsec) 1287 1288 ! Optionally provide output for momentum distribution of annihilation radiation 1289 if (dtset%posdoppler>0) then 1290 call timab(1174,1,tsec) 1291 call posdoppler(cg,cprj,crystal,dimcprj,dtfil,dtset,electronpositron,psps%filpsp,& 1292 & kg,mcg,mcprj,mpi_enreg,my_natom,n3xccc,nfft,ngfft,nhat,npwarr,& 1293 & occ,pawang,pawrad,pawrhoij,pawtab,rhor,xccc3d) 1294 call timab(1174,2,tsec) 1295 end if 1296 end if 1297 1298 !Optionally provide output for WanT 1299 if (dtset%prtwant==1) then 1300 call timab(1175,1,tsec) 1301 ! WARNING: mpi_enreg not used --> MPI is not supported 1302 call outwant(dtset,eigen,cg,kg,npwarr,mband,mcg,nkpt,nsppol,mkmem,mpw,dtset%prtwant) 1303 call timab(1175,2,tsec) 1304 end if 1305 1306 !Optionally provide output for electric field gradient calculation 1307 if (dtset%nucefg > 0) then 1308 call timab(1176,1,tsec) 1309 call calc_efg(mpi_enreg,my_natom,natom,nfft,ngfft,nhat,nspden,dtset%nsym,dtset%nucefg,& 1310 & ntypat,paw_an,pawang,pawrad,pawrhoij,pawtab,& 1311 & dtset%ptcharge,dtset%quadmom,rhor,rprimd,dtset%symrel,& 1312 & dtset%tnons,dtset%typat,ucvol,psps%usepaw,xred,psps%zionpsp,& 1313 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1314 call timab(1176,2,tsec) 1315 end if 1316 1317 !Optionally provide output for Fermi-contact term at nuclear positions 1318 if (dtset%nucfc > 0) then 1319 call timab(1177,1,tsec) 1320 call calc_fc(my_natom,natom,nspden,ntypat,pawrad,pawrhoij,pawtab,dtset%typat,psps%usepaw,& 1321 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1322 call timab(1177,2,tsec) 1323 end if 1324 1325 ! Output electron bands. 1326 if (me == master .and. dtset%tfkinfunc==0) then 1327 call timab(1178,1,tsec) 1328 if (size(dtset%kptbounds, dim=2) > 0) then 1329 call ebands_write(ebands, dtset%prtebands, dtfil%filnam_ds(4), kptbounds=dtset%kptbounds) 1330 else 1331 call ebands_write(ebands, dtset%prtebands, dtfil%filnam_ds(4)) 1332 end if 1333 call timab(1178,2,tsec) 1334 end if 1335 1336 !Optionally provide Xcrysden output for the Fermi surface (Only master writes) 1337 if (me == master .and. dtset%prtfsurf == 1) then 1338 call timab(1179,1,tsec) 1339 if (ebands_write_bxsf(ebands,crystal,dtfil%fnameabo_app_bxsf) /= 0) then 1340 msg = "Cannot produce BXSF file with Fermi surface, see log file for more info" 1341 ABI_WARNING(msg) 1342 call wrtout(ab_out, msg) 1343 end if 1344 call timab(1179,2,tsec) 1345 end if ! prtfsurf==1 1346 1347 !output nesting factor for Fermi surface (requires ph_nqpath) 1348 if (me == master .and. dtset%prtnest>0 .and. dtset%ph_nqpath > 0) then 1349 call timab(1180,1,tsec) 1350 ierr = ebands_write_nesting(ebands,crystal,dtfil%fnameabo_app_nesting,dtset%prtnest,& 1351 dtset%tsmear,dtset%fermie_nest,dtset%ph_qpath(:,1:dtset%ph_nqpath),msg) 1352 if (ierr /= 0) then 1353 ABI_WARNING(msg) 1354 call wrtout(ab_out, msg) 1355 end if 1356 call timab(1180,2,tsec) 1357 end if ! prtnest=1 1358 1359 if (dtset%prtdipole == 1) then 1360 call timab(1181,1,tsec) 1361 call multipoles_out(rhor,mpi_enreg,natom,nfft,ngfft,dtset%nspden,dtset%ntypat,rprimd,& 1362 dtset%typat,ucvol,ab_out,xred,dtset%ziontypat) 1363 call timab(1181,2,tsec) 1364 end if 1365 1366 ! BoltzTraP output files in GENEric format 1367 if (dtset%prtbltztrp == 1 .and. me==master)then 1368 call timab(1182,1,tsec) 1369 call ebands_prtbltztrp(ebands, crystal, dtfil%filnam_ds(4)) 1370 call timab(1182,2,tsec) 1371 endif 1372 1373 ! Band structure interpolation from eigenvalues computed on the k-mesh. 1374 if (nint(dtset%einterp(1)) /= 0 .and. dtset%kptopt > 0) then 1375 call timab(1183,1,tsec) 1376 call ebands_interpolate_kpath(ebands, dtset, crystal, [0, 0], dtfil%filnam_ds(4), spacecomm) 1377 call timab(1183,2,tsec) 1378 end if 1379 1380 ABI_SFREE_PTR(elfr) 1381 ABI_SFREE_PTR(grhor) 1382 ABI_SFREE_PTR(lrhor) 1383 1384 call crystal%free() 1385 call ebands_free(ebands) 1386 1387 ! Destroy atom table used for parallelism 1388 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 1389 1390 call timab(1150,2,tsec) ! outscfcv 1391 1392 DBG_EXIT("COLL") 1393 1394 end subroutine outscfcv