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