TABLE OF CONTENTS
ABINIT/ks_ddiago [ Functions ]
NAME
ks_ddiago
FUNCTION
This routine performs the direct diagonalization of the Kohn-Sham Hamiltonian
for a given k-point and spin. The routine drives the following operations:
1) Re-computing <G|H|G_prim> matrix elements for all (G, G_prim).
starting from the knowledge of the local potential on the real-space FFT mesh.
2) Diagonalizing H in the plane-wave basis.
It is called in outkss.F90 during the generation of the KSS file
needed for a GW post-treatment. Since many-body calculations usually
require a large number of eigenstates eigen-functions, a direct
diagonalization of the Hamiltonian might reveal more stable than iterative
techniques that might be problematic when several high energy states are required.
The main drawback of the direct diagonalization is the bad scaling with the size
of the basis set (npw**3) and the large memory requirements.
At present, only norm-conserving pseudopotentials are implemented.
COPYRIGHT
Copyright (C) 2000-2018 ABINIT group (MG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
INPUTS
kpoint(3) prtvol=Integer Flags defining verbosity level ecut=cut-off energy for plane wave basis sphere (Ha) mgfftc=maximum size of 1D FFTs (coarse mesh). natom=number of atoms in cell. nfftf=(effective) number of FFT grid points in the dense FFT mesh (for this processor) (nfftf=nfft for norm-conserving potential runs) nspinor=number of spinorial components of the wavefunctions nsppol=1 for unpolarized, 2 for spin-polarized nspden=number of density components Pawtab(Psps%ntypat*Psps%usepaw) <type(pawtab_type)>=paw tabulated starting data Pawfgr<pawfgr_type>=fine grid parameters and related data Paw_ij(natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels Psps <type(pseudopotential_type)>=variables related to pseudopotentials rprimd(3,3)=dimensional primitive translations for real space (bohr) vtrial(nfftf,nspden)=the trial potential xred(3,natom)=reduced dimensionless atomic coordinates comm=MPI communicator. [Electronpositron] <electronpositron_type>=quantities for the electron-positron annihilation. nfftc=Number of points in the coarse FFT mesh. ngfftc(18)=Info about 3D FFT for the coarse mesh, see ~abinit/doc/variables/vargs.htm#ngfft Diago_ctl<ddiago_ctl_type>=Datatype storing variables and options controlling the direct diagonalization.
OUTPUT
ierr=Status error. onband_diago
SIDE EFFECTS
eig_ene(:)=Pointer used for allocating and storing the eigenvalues (hartree)
input: pointer to NULL
output: eig_ene(onband_diago)=The calculatated eigenvalues in ascending order.
eig_vec(:,:,:)=Pointer used for allocating and holding the wave functions at this k-point and spin.
input: pointer to NULL
output: eig_vec(2,npw_k*nspinor,onband_diago)=The calculated eigenvectors.
Cprj_k(natom,nspinor*onband_diago) PAW only===
input: pointer to NULL
output: Projected eigenstates <Proj_i|Cnk> from output eigenstates.
NOTES
* The routine can be time consuming (in particular when computing <G1|H|G2> elements for all (G1,G2)). So, it is recommended to call it once per run. * The routine RE-compute all Hamiltonian terms. So it is equivalent to an additional electronic SC cycle. (This has no effect is convergence was reach. If not, eigenvalues/vectors may differs from the conjugate gradient ones)
NOTES
Please, do NOT pass Dtset% to this routine. Either use a local variable properly initialized or add the additional variable to ddiago_ctl_type and change the creation method accordingly. ks_ddiago is designed such that it is possible to diagonalize the Hamiltonian at an arbitrary k-point or spin (not efficient but easy to code). Therefore ks_ddiago is useful non only for the KSS generation but also for testing more advanced iterative algorithms as well as interpolation techniques.
PARENTS
m_shirley,outkss
CHILDREN
destroy_hamiltonian,destroy_mpi_enreg,fftpac,getcprj,getghc init_distribfft_seq,init_hamiltonian,initmpi_seq,initylmg,kpgsph load_k_hamiltonian,load_spin_hamiltonian,metric,mkffnl,mkkin,mkkpg pawcprj_alloc,pawcprj_free,pawcprj_reorder,transgrid,wrtout,xheev xheevx,xhegv,xhegvx,xmpi_barrier
SOURCE
98 #if defined HAVE_CONFIG_H 99 #include "config.h" 100 #endif 101 102 #include "abi_common.h" 103 104 subroutine ks_ddiago(Diago_ctl,nband_k,nfftc,mgfftc,ngfftc,natom,& 105 & typat,nfftf,nspinor,nspden,nsppol,Pawtab,Pawfgr,Paw_ij,& 106 & Psps,rprimd,vtrial,xred,onband_diago,eig_ene,eig_vec,Cprj_k,comm,ierr,& 107 & Electronpositron) ! Optional arguments 108 109 use defs_basis 110 use defs_datatypes 111 use defs_abitypes 112 use m_profiling_abi 113 use m_xmpi 114 use m_errors 115 use m_hamiltonian 116 117 use m_geometry, only : metric 118 use m_abilasi, only : xheev, xhegv, xheevx, xhegvx 119 use m_electronpositron, only : electronpositron_type 120 use m_fftcore, only : kpgsph 121 use m_mpinfo, only : destroy_mpi_enreg 122 use m_pawtab, only : pawtab_type 123 use m_paw_ij, only : paw_ij_type 124 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, & 125 & pawcprj_reorder 126 use m_pawfgr, only : pawfgr_type 127 use m_kg, only : mkkin 128 129 !This section has been created automatically by the script Abilint (TD). 130 !Do not modify the following lines by hand. 131 #undef ABI_FUNC 132 #define ABI_FUNC 'ks_ddiago' 133 use interfaces_14_hidewrite 134 use interfaces_51_manage_mpi 135 use interfaces_53_ffts 136 use interfaces_56_recipspace 137 use interfaces_65_paw 138 use interfaces_66_nonlocal 139 use interfaces_66_wfs 140 !End of the abilint section 141 142 implicit none 143 144 !Arguments ------------------------------------ 145 !scalars 146 integer,intent(in) :: mgfftc,natom,comm,nband_k 147 integer,intent(in) :: nfftf,nsppol,nspden,nspinor,nfftc 148 integer,intent(out) :: ierr,onband_diago 149 type(pseudopotential_type),intent(in) :: Psps 150 type(pawfgr_type),intent(in) :: Pawfgr 151 type(ddiago_ctl_type),intent(in) :: Diago_ctl 152 !arrays 153 integer,intent(in) :: typat(natom) 154 integer,intent(in) :: ngfftc(18) 155 real(dp),intent(in) :: rprimd(3,3) 156 real(dp),intent(inout) :: vtrial(nfftf,nspden) 157 real(dp),intent(in) :: xred(3,natom) 158 real(dp),pointer :: eig_ene(:),eig_vec(:,:,:) 159 type(pawcprj_type),pointer :: Cprj_k(:,:) 160 type(pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Psps%usepaw) 161 type(paw_ij_type),intent(in) :: Paw_ij(natom*Psps%usepaw) 162 type(electronpositron_type),optional,pointer :: Electronpositron 163 164 !Local variables------------------------------- 165 !scalars 166 integer,parameter :: mkmem_=1,tim_getghc=4,paral_kgb=0 167 integer :: cprj_choice,cpopt,dimffnl,ib,ider,idir,isppol,npw_k 168 integer :: ikg,master,istwf_k,exchn2n3d,prtvol 169 integer :: jj,n1,n2,n3,n4,n5,n6,negv,nkpg,nprocs,npw_k_test 170 integer :: my_rank,optder 171 integer :: ispden,ndat,type_calc,sij_opt,igsp2,cplex_ghg 172 integer :: iband,iorder_cprj,ibs1,ibs2 173 real(dp) :: cfact,ucvol,ecutsm,effmass_free,lambda,size_mat,ecut 174 logical :: do_full_diago 175 character(len=50) :: jobz,range 176 character(len=80) :: frmt1,frmt2 177 character(len=10) :: stag(2) 178 character(len=500) :: msg 179 type(MPI_type) :: MPI_enreg_seq 180 type(gs_hamiltonian_type) :: gs_hamk 181 !arrays 182 integer :: nloalg(3) 183 integer,allocatable :: kg_k(:,:) 184 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),kptns_(3,1) 185 real(dp) :: kpoint(3),ylmgr_dum(1,1,1) 186 real(dp) :: rhodum(1) 187 real(dp),allocatable :: ph3d(:,:,:),pwave(:,:) 188 real(dp),allocatable :: ffnl(:,:,:,:),kinpw(:),kpg_k(:,:) 189 real(dp),allocatable :: vlocal(:,:,:,:),ylm_k(:,:),dum_ylm_gr_k(:,:,:),vlocal_tmp(:,:,:) 190 real(dp),allocatable :: ghc(:,:),gvnlc(:,:),gsc(:,:) 191 real(dp),allocatable :: ghg_mat(:,:,:),gtg_mat(:,:,:) 192 real(dp),allocatable :: cgrvtrial(:,:) 193 real(dp),pointer :: cwavef(:,:) 194 type(pawcprj_type),allocatable :: Cwaveprj(:,:) 195 196 ! ********************************************************************* 197 198 DBG_ENTER("COLL") 199 200 nprocs = xmpi_comm_size(comm) 201 my_rank = xmpi_comm_rank(comm) 202 master=0 203 204 if (nprocs>1) then 205 MSG_WARNING(" ks_ddiago not supported in parallel. Running in sequential.") 206 end if 207 208 call initmpi_seq(MPI_enreg_seq) ! Fake MPI_type for sequential part. 209 call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfftc(2),ngfftc(3),'all') 210 if (Pawfgr%usefinegrid/=0) then 211 call init_distribfft_seq(MPI_enreg_seq%distribfft,'f',Pawfgr%ngfft(2),pawfgr%ngfft(3),'all') 212 end if 213 214 isppol = Diago_ctl%isppol 215 kpoint = Diago_ctl%kpoint 216 istwf_k = Diago_ctl%istwf_k 217 !% nband_k = Diago_ctl%nband_k 218 npw_k = Diago_ctl%npw_k 219 nloalg = Diago_ctl%nloalg 220 ecut = Diago_ctl%ecut 221 ecutsm = Diago_ctl%ecutsm 222 effmass_free = Diago_ctl%effmass_free 223 prtvol = Diago_ctl%prtvol 224 225 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 226 227 if (nsppol==1) stag=(/' ',' '/) 228 if (nsppol==2) stag=(/'SPIN UP: ','SPIN DOWN:'/) 229 ! 230 !The coarse FFT mesh. 231 n1=ngfftc(1); n2=ngfftc(2); n3=ngfftc(3) 232 n4=ngfftc(4); n5=ngfftc(5); n6=ngfftc(6) 233 ! 234 !==================== 235 !=== Check input ==== 236 !==================== 237 ierr=0 238 ! 239 !* istwfk must be 1 for each k-point 240 if (istwf_k/=1) then 241 write(msg,'(7a)')& 242 & ' istwfk/=1 not allowed:',ch10,& 243 & ' States output not programmed for time-reversal symmetry.',ch10,& 244 & ' Action : change istwfk in input file (put it to 1 for all kpt).',ch10,& 245 & ' Program does not stop but _KSS file will not be created...' 246 MSG_WARNING(msg) 247 ierr=ierr+1 248 end if 249 250 if (paral_kgb/=0) then 251 write(msg,'(3a)')& 252 & ' paral_kgb/=0 not allowed:',ch10,& 253 & ' Program does not stop but _KSS file will not be created...' 254 MSG_WARNING(msg) 255 ierr=ierr+1 256 end if 257 258 if (ierr/=0) RETURN ! Houston we have a problem! 259 ! 260 !Initialize the Hamiltonian datatype on the coarse FFT mesh. 261 if (PRESENT(Electronpositron)) then 262 call init_hamiltonian(gs_hamk,Psps,pawtab,nspinor,nsppol,nspden,natom,typat,xred,nfftc,& 263 & mgfftc,ngfftc,rprimd,nloalg,paw_ij=Paw_ij,usecprj=0,Electronpositron=Electronpositron) 264 else 265 call init_hamiltonian(gs_hamk,Psps,pawtab,nspinor,nsppol,nspden,natom,typat,xred,nfftc,& 266 & mgfftc,ngfftc,rprimd,nloalg,paw_ij=Paw_ij,usecprj=0) 267 end if 268 269 !Check on the number of stored bands. 270 if (nband_k==-1.or.nband_k>=npw_k*nspinor) then 271 onband_diago=npw_k*nspinor 272 write(msg,'(4a)')ch10,& 273 & ' Since the number of bands to be computed was (-1) or',ch10,& 274 & ' too large, it has been set to the max. value npw_k*nspinor. ' 275 call wrtout(std_out,msg,'COLL') 276 else 277 onband_diago=nband_k 278 end if 279 280 !do_full_diago = (onband_diago==npw_k*nspinor) 281 do_full_diago = Diago_ctl%do_full_diago 282 283 if (do_full_diago) then 284 write(msg,'(6a)')ch10,& 285 & ' Since the number of bands to be computed',ch10,& 286 & ' is equal to the nb of G-vectors found for this k-pt,',ch10,& 287 & ' the program will perform complete diagonalizations.' 288 else 289 write(msg,'(6a)')ch10,& 290 & ' Since the number of bands to be computed',ch10,& 291 & ' is less than the number of G-vectors found,',ch10,& 292 & ' the program will perform partial diagonalizations.' 293 end if 294 if (prtvol>0) then 295 call wrtout(std_out,msg,'COLL') 296 end if 297 ! 298 299 !* Set up local potential vlocal with proper dimensioning, from vtrial. 300 !* Select spin component of interest if nspden<=2 as nvloc==1, for nspden==4, nvloc==4 301 !* option=2: vtrial(n1*n2*n3,ispden) --> vlocal(nd1,nd2,nd3) real case 302 303 !nvloc=1; if (nspden==4) nvloc=4 304 ABI_MALLOC(vlocal,(n4,n5,n6,gs_hamk%nvloc)) 305 306 if (nspden/=4)then 307 if (Psps%usepaw==0.or.Pawfgr%usefinegrid==0) then 308 call fftpac(isppol,MPI_enreg_seq,nspden,n1,n2,n3,n4,n5,n6,ngfftc,vtrial,vlocal,2) 309 else ! Move from fine to coarse FFT mesh (PAW) 310 ABI_MALLOC(cgrvtrial,(nfftc,nspden)) 311 call transgrid(1,MPI_enreg_seq,nspden,-1,0,0,paral_kgb,Pawfgr,rhodum,rhodum,cgrvtrial,vtrial) 312 call fftpac(isppol,MPI_enreg_seq,nspden,n1,n2,n3,n4,n5,n6,ngfftc,cgrvtrial,vlocal,2) 313 ABI_FREE(cgrvtrial) 314 end if 315 else 316 ABI_MALLOC(vlocal_tmp,(n4,n5,n6)) 317 if (Psps%usepaw==0.or.Pawfgr%usefinegrid==0) then 318 do ispden=1,nspden 319 call fftpac(ispden,MPI_enreg_seq,nspden,n1,n2,n3,n4,n5,n6,ngfftc,vtrial,vlocal_tmp,2) 320 vlocal(:,:,:,ispden)=vlocal_tmp(:,:,:) 321 end do 322 else ! Move from fine to coarse FFT mesh (PAW) 323 ABI_MALLOC(cgrvtrial,(nfftc,nspden)) 324 call transgrid(1,MPI_enreg_seq,nspden,-1,0,0,paral_kgb,Pawfgr,rhodum,rhodum,cgrvtrial,vtrial) 325 do ispden=1,nspden 326 call fftpac(ispden,MPI_enreg_seq,nspden,n1,n2,n3,n4,n5,n6,ngfftc,cgrvtrial,vlocal_tmp,2) 327 vlocal(:,:,:,ispden)=vlocal_tmp(:,:,:) 328 end do 329 ABI_FREE(cgrvtrial) 330 end if 331 ABI_FREE(vlocal_tmp) 332 end if 333 334 !Continue to initialize the Hamiltonian (spin-dependent part) 335 call load_spin_hamiltonian(gs_hamk,isppol,vlocal=vlocal,with_nonlocal=.true.) 336 ! 337 !* Calculate G-vectors, for this k-point. 338 !* Count the number of planewaves as a check. 339 exchn2n3d=0; ikg=0 340 ABI_MALLOC(kg_k,(3,npw_k)) 341 342 call kpgsph(ecut,exchn2n3d,gmet,ikg,0,istwf_k,kg_k,kpoint,0,MPI_enreg_seq,0,npw_k_test) 343 ABI_CHECK(npw_k_test==npw_k,"npw_k_test/=npw_k") 344 345 call kpgsph(ecut,exchn2n3d,gmet,ikg,0,istwf_k,kg_k,kpoint,mkmem_,MPI_enreg_seq,npw_k,npw_k_test) 346 347 !======================== 348 !==== Kinetic energy ==== 349 !======================== 350 ABI_MALLOC(kinpw,(npw_k)) 351 ! call mkkin(ecut,ecutsm,effmass_free,gmet,kg_k,kinpw,kpoint,npw_k) 352 call mkkin(ecut,ecutsm,effmass_free,gmet,kg_k,kinpw,kpoint,npw_k,0,0) 353 ! 354 !================================ 355 !==== Non-local form factors ==== 356 !================================ 357 ABI_MALLOC(ylm_k,(npw_k,Psps%mpsang**2*Psps%useylm)) 358 359 if (Psps%useylm==1) then 360 optder=0 361 ABI_MALLOC(dum_ylm_gr_k,(npw_k,3+6*(optder/2),Psps%mpsang**2)) 362 kptns_(:,1) = kpoint 363 364 ! Here mband is not used if paral_compil_kpt=0 365 call initylmg(gprimd,kg_k,kptns_,mkmem_,MPI_enreg_seq,Psps%mpsang,npw_k,(/nband_k/),1,& 366 & (/npw_k/),1,optder,rprimd,ylm_k,dum_ylm_gr_k) 367 368 ABI_FREE(dum_ylm_gr_k) 369 end if 370 371 !Compute (k+G) vectors (only if useylm=1) 372 nkpg=3*nloalg(3) 373 ABI_MALLOC(kpg_k,(npw_k,nkpg)) 374 if (nkpg>0) then 375 call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k) 376 end if 377 378 !Compute nonlocal form factors ffnl at all (k+G): 379 idir=0; ider=0; dimffnl=1+ider ! Now the derivative is not needed anymore. 380 ABI_MALLOC(ffnl,(npw_k,dimffnl,Psps%lmnmax,Psps%ntypat)) 381 382 call mkffnl(Psps%dimekb,dimffnl,Psps%ekb,ffnl,Psps%ffspl,gmet,gprimd,ider,idir,Psps%indlmn,& 383 & kg_k,kpg_k,kpoint,Psps%lmnmax,Psps%lnmax,Psps%mpsang,Psps%mqgrid_ff,nkpg,npw_k,& 384 & Psps%ntypat,Psps%pspso,Psps%qgrid_ff,rmet,Psps%usepaw,Psps%useylm,ylm_k,ylmgr_dum) 385 386 ABI_FREE(ylm_k) 387 388 !Load k-dependent part in the Hamiltonian datastructure 389 ABI_MALLOC(ph3d,(2,npw_k,gs_hamk%matblk)) 390 call load_k_hamiltonian(gs_hamk,kpt_k=kpoint,istwf_k=istwf_k,npw_k=npw_k,kinpw_k=kinpw,& 391 & kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl,ph3d_k=ph3d,& 392 & compute_ph3d=.true.,compute_gbound=.true.) 393 394 !Prepare the call to getghc. 395 ndat=1; lambda=zero; type_calc=0 ! For applying the whole Hamiltonian 396 sij_opt=0; if (Psps%usepaw==1) sij_opt=1 ! For PAW, <k+G|1+S|k+G"> is also needed. 397 398 cpopt=-1 ! If cpopt=-1, <p_lmn|in> (and derivatives) are computed here (and not saved) 399 if (Psps%usepaw==1.and..FALSE.) then ! TODO Calculate <p_lmn|k+G>. 400 cpopt = 0 ! <p_lmn|in> are computed here and saved 401 end if 402 403 ABI_MALLOC(ghc ,(2,npw_k*nspinor*ndat)) 404 ABI_MALLOC(gvnlc,(2,npw_k*nspinor*ndat)) 405 ABI_MALLOC(gsc ,(2,npw_k*nspinor*ndat*(sij_opt+1)/2)) 406 407 cplex_ghg=2 408 size_mat = cplex_ghg*(npw_k*nspinor)**2*dp*b2Mb 409 !; if (Psps%usepaw==1) size_mat=two*size_mat 410 !write(msg,'(a,f0.3,a)')" Memory required by the Hamiltonian matrix: ",size_mat," [Mb]." 411 !call wrtout(std_out,msg,"COLL") 412 413 write(msg,'(a,f0.3,a)')" Out-of-memory in ghg_mat. Memory required by the Hamiltonian matrix: ",size_mat," [Mb]." 414 ABI_STAT_MALLOC(ghg_mat,(cplex_ghg,npw_k*nspinor,npw_k*nspinor), ierr) 415 ABI_CHECK(ierr==0, msg) 416 417 write(msg,'(a,f0.3,a)')" Out-of-memory in gtg_mat. Memory required by the PAW overlap operator: ",size_mat," [Mb]." 418 ABI_STAT_MALLOC(gtg_mat,(cplex_ghg,npw_k*nspinor,npw_k*nspinor*Psps%usepaw), ierr) 419 ABI_CHECK(ierr==0, msg) 420 421 ABI_DT_MALLOC(Cwaveprj,(natom,nspinor*(1+cpopt)*gs_hamk%usepaw)) 422 if (cpopt==0) then ! Cwaveprj is ordered, see nonlop_ylm. 423 call pawcprj_alloc(Cwaveprj,0,gs_hamk%dimcprj) 424 end if 425 426 ABI_MALLOC(pwave,(2,npw_k*nspinor)) 427 pwave=zero ! Initialize plane-wave array: 428 429 if (prtvol>0) then 430 call wrtout(std_out,' Calculating <G|H|G''> elements','PERS') 431 end if 432 433 do igsp2=1,npw_k*nspinor ! Loop over the |beta,G''> component. 434 435 pwave(1,igsp2)=one ! Get <:|H|beta,G''> and <:|T_{PAW}|beta,G''> 436 437 call getghc(cpopt,pwave,Cwaveprj,ghc,gsc,gs_hamk,gvnlc,lambda,MPI_enreg_seq,ndat,& 438 & prtvol,sij_opt,tim_getghc,type_calc) 439 440 ! Fill the upper triangle. 441 ghg_mat(:,1:igsp2,igsp2) = ghc(:,1:igsp2) 442 if (Psps%usepaw==1) gtg_mat(:,1:igsp2,igsp2) = gsc(:,1:igsp2) 443 444 pwave(1,igsp2)=zero ! Reset the |G,beta> component that has been treated. 445 end do 446 447 !Free workspace memory allocated so far. 448 ABI_FREE(pwave) 449 ABI_FREE(kinpw) 450 ABI_FREE(vlocal) 451 ABI_FREE(ghc) 452 ABI_FREE(gvnlc) 453 ABI_FREE(gsc) 454 455 if (Psps%usepaw==1.and.cpopt==0) then 456 call pawcprj_free(Cwaveprj) 457 end if 458 ABI_DT_FREE(Cwaveprj) 459 ! 460 !=========================================== 461 !=== Diagonalization of <G|H|G''> matrix === 462 !=========================================== 463 ! 464 !*** Allocate the pointers *** 465 ABI_MALLOC(eig_ene,(onband_diago)) 466 ABI_MALLOC(eig_vec,(cplex_ghg,npw_k*nspinor,onband_diago)) 467 468 jobz =Diago_ctl%jobz !jobz="Vectors" 469 470 if (do_full_diago) then ! * Complete diagonalization 471 472 write(msg,'(2a,3es16.8,3x,3a,i5)')ch10,& 473 & ' Begin complete diagonalization for kpt= ',kpoint(:),stag(isppol),ch10,& 474 & ' - Size of mat.=',npw_k*nspinor 475 if (prtvol>0) then 476 call wrtout(std_out,msg,'PERS') 477 end if 478 479 if (Psps%usepaw==0) then 480 call xheev( jobz,"Upper",cplex_ghg,npw_k*nspinor,ghg_mat,eig_ene) 481 else 482 call xhegv(1,jobz,"Upper",cplex_ghg,npw_k*nspinor,ghg_mat,gtg_mat,eig_ene) 483 end if 484 485 eig_vec(:,:,:)= ghg_mat 486 487 else ! * Partial diagonalization 488 489 range=Diago_ctl%range !range="Irange" 490 491 write(msg,'(2a,3es16.8,3a,i5,a,i5)')ch10,& 492 & ' Begin partial diagonalization for kpt= ',kpoint,stag(isppol),ch10,& 493 & ' - Size of mat.=',npw_k*nspinor,' - # bnds=',onband_diago 494 if (prtvol>0) then 495 call wrtout(std_out,msg,'PERS') 496 end if 497 498 if (Psps%usepaw==0) then 499 call xheevx(jobz,range,"Upper",cplex_ghg,npw_k*nspinor,ghg_mat,zero,zero,& 500 & 1,onband_diago,-tol8,negv,eig_ene,eig_vec,npw_k*nspinor) 501 else 502 call xhegvx(1,jobz,range,"Upper",cplex_ghg,npw_k*nspinor,ghg_mat,gtg_mat,zero,zero,& 503 & 1,onband_diago,-tol8,negv,eig_ene,eig_vec,npw_k*nspinor) 504 end if 505 506 end if 507 508 ABI_FREE(ghg_mat) 509 ABI_FREE(gtg_mat) 510 ! 511 !======================================================== 512 !==== Calculate <Proj_i|Cnk> from output eigenstates ==== 513 !======================================================== 514 if (Psps%usepaw==1) then 515 516 iorder_cprj=1 ! Ordered (order does change wrt input file); will be changed later 517 ABI_DT_MALLOC(Cprj_k,(natom,nspinor*onband_diago)) 518 call pawcprj_alloc(Cprj_k,0,gs_hamk%dimcprj) 519 520 idir=0; cprj_choice=1 ! Only projected wave functions. 521 522 do iband=1,onband_diago 523 ibs1=nspinor*(iband-1)+1 524 ibs2=ibs1; if (nspinor==2) ibs2=ibs2+1 525 cwavef => eig_vec(1:2,1:npw_k,iband) 526 527 call getcprj(cprj_choice,0,cwavef,Cprj_k(:,ibs1:ibs2),& 528 & gs_hamk%ffnl_k,idir,gs_hamk%indlmn,gs_hamk%istwf_k,gs_hamk%kg_k,& 529 & gs_hamk%kpg_k,gs_hamk%kpt_k,gs_hamk%lmnmax,gs_hamk%mgfft,MPI_enreg_seq,& 530 & gs_hamk%natom,gs_hamk%nattyp,gs_hamk%ngfft,gs_hamk%nloalg,gs_hamk%npw_k,gs_hamk%nspinor,& 531 & gs_hamk%ntypat,gs_hamk%phkxred,gs_hamk%ph1d,gs_hamk%ph3d_k,gs_hamk%ucvol,gs_hamk%useylm) 532 end do 533 534 ! Reorder the cprj (order is now the same as in input file) 535 call pawcprj_reorder(Cprj_k,gs_hamk%atindx1) 536 537 ! deallocate(cwavef) 538 end if !usepaw==1 539 540 if (prtvol>0) then ! Write eigenvalues. 541 cfact=Ha_eV ; frmt1='(8x,9(1x,f7.2))' ; frmt2='(8x,9(1x,f7.2))' 542 write(msg,'(a,3es16.8,3x,a)')' Eigenvalues in eV for kpt= ',kpoint,stag(isppol) 543 call wrtout(std_out,msg,'COLL') 544 545 write(msg,frmt1)(eig_ene(ib)*cfact,ib=1,MIN(9,onband_diago)) 546 call wrtout(std_out,msg,'COLL') 547 if (onband_diago>9) then 548 do jj=10,onband_diago,9 549 write(msg,frmt2) (eig_ene(ib)*cfact,ib=jj,MIN(jj+8,onband_diago)) 550 call wrtout(std_out,msg,'COLL') 551 end do 552 end if 553 end if 554 555 ABI_FREE(kpg_k) 556 ABI_FREE(kg_k) 557 ABI_FREE(ph3d) 558 ABI_FREE(ffnl) 559 560 call destroy_mpi_enreg(MPI_enreg_seq) 561 call destroy_hamiltonian(gs_hamk) 562 563 call xmpi_barrier(comm) 564 565 DBG_EXIT("COLL") 566 567 end subroutine ks_ddiago