TABLE OF CONTENTS
ABINIT/setnoccmmp [ Functions ]
NAME
setnoccmmp
FUNCTION
PAW+U only:
Compute density matrix nocc_{m,m_prime}
or
Impose value of density matrix using dmatpawu input array, then symetrize it.
noccmmp^{\sigma}_{m,m'}=\sum_{ni,nj}[\rho^{\sigma}_{ni,nj}*phiphjint_{ni,nj}]
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (BA,MT) 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
compute_dmat= flag: if 1, nocc_{m,mp} is computed
dimdmat=first dimension of dmatpawu array
dmatpawu(dimdmat,dimdmat,nsppol*nspinor,natpawu)=input density matrix to be copied into noccmpp
dmatudiag= flag controlling the use of diagonalization:
0: no diagonalization of nocc_{m,mp}
1: diagonalized nocc_{m,mp} matrix is printed
2: dmatpawu matrix is expressed in the basis where nocc_(m,mp} is diagonal
impose_dmat= flag: if 1, nocc_{m,mp} is replaced by dmatpawu
indsym(4,nsym,natom)=indirect indexing array for atom labels
mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
comm_atom=--optional-- MPI communicator over atoms
my_natom=number of atoms treated by current processor
natom=number of atoms in cell
natpawu=number of atoms on which PAW+U is applied
nspinor=number of spinorial components of the wavefunctions
nsppol=number of independant spin components
nsym=number of symmetry elements in space group
ntypat=number of atom types
paw_ij(my_natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels
pawang <type(pawang_type)>=paw angular mesh and related data
pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
spinat(3,matom)=initial spin of each atom, in unit of hbar/2
symafm(nsym)=(anti)ferromagnetic part of symmetry operations
typat(natom)=type for each atom
useexexch=1 if local-exact-exchange is activated
usepawu=1 if PAW+U is activated
OUTPUT
paw_ij(natom)%noccmmp(cplex_dij,2*lpawu+1,2*lpawu+1,nsppol or ndij)=density matrix
NOTES
For non-collinear magnetism,
- nocc_{m,mp} is computed as:noccmmp(:,:,:,1)= n{m,mp}
noccmmp(:,:,:,2)= m_x{m,mp}
noccmmp(:,:,:,3)= m_y{m,mp}
noccmmp(:,:,:,4)= m_z{m,mp}
- but nocc_{m,mp} is stored as: noccmmp(:,:,:,1)= n^{up,up}_{m,mp}
noccmmp(:,:,:,2)= n^{dn,dn}_{m,mp}
noccmmp(:,:,:,3)= n^{up,dn}_{m,mp}
noccmmp(:,:,:,4)= n^{dn,up}_{m,mp}
We choose to have noccmmp complex when ndij=4 (ie nspinor=2)
If ndij=4 and pawspnorb=0, one could keep noccmmp real
with the n11, n22, Re(n12), Im(n21) representation, but it would
less clear to change the representation when pawspnorb is activated.
If ndij=4, nocc_{m,mp} is transformed to the Ylm basis
and then to the J, M_J basis (if cplex_dij==2)
Note that n_{m,mp}=<mp|hat(n)|m> because rhoij=<p_j|...|p_i>
PARENTS
afterscfloop,pawdenpot,pawprt,scfcv,vtorho
CHILDREN
dgemm,dsyev,free_my_atmtab,get_my_atmtab,mat_mlms2jmj,mat_slm2ylm wrtout,zgemm,zheev
SOURCE
81 #if defined HAVE_CONFIG_H 82 #include "config.h" 83 #endif 84 85 #include "abi_common.h" 86 87 subroutine setnoccmmp(compute_dmat,dimdmat,dmatpawu,dmatudiag,impose_dmat,indsym,my_natom,natom,& 88 & natpawu,nspinor,nsppol,nsym,ntypat,paw_ij,pawang,pawprtvol,pawrhoij,pawtab,& 89 & spinat,symafm,typat,useexexch,usepawu, & 90 & mpi_atmtab,comm_atom) ! optional arguments (parallelism) 91 92 use defs_basis 93 use m_profiling_abi 94 use m_errors 95 use m_xmpi, only : xmpi_comm_self 96 use m_linalg_interfaces 97 98 use m_pawang, only : pawang_type, mat_mlms2jmj, mat_slm2ylm 99 use m_pawtab, only : pawtab_type 100 use m_paw_ij, only : paw_ij_type 101 use m_pawrhoij, only : pawrhoij_type 102 use m_paw_sphharm, only : mat_mlms2jmj, mat_slm2ylm 103 use m_paral_atom, only : get_my_atmtab, free_my_atmtab 104 105 !This section has been created automatically by the script Abilint (TD). 106 !Do not modify the following lines by hand. 107 #undef ABI_FUNC 108 #define ABI_FUNC 'setnoccmmp' 109 use interfaces_14_hidewrite 110 !End of the abilint section 111 112 implicit none 113 114 !Arguments --------------------------------------------- 115 !scalars 116 integer,intent(in) :: compute_dmat,dimdmat,dmatudiag,impose_dmat,my_natom,natom,natpawu 117 integer,intent(in) :: nspinor,nsppol,nsym,ntypat,useexexch,usepawu 118 integer,optional,intent(in) :: comm_atom 119 type(pawang_type),intent(in) :: pawang 120 integer,intent(in) :: pawprtvol 121 !arrays 122 integer,intent(in) :: indsym(4,nsym,natom),symafm(nsym),typat(natom) 123 integer,optional,target,intent(in) :: mpi_atmtab(:) 124 real(dp),intent(in) :: dmatpawu(dimdmat,dimdmat,nspinor*nsppol,natpawu*impose_dmat) 125 !real(dp),intent(in) :: dmatpawu(:,:,:,:) 126 real(dp),intent(in) :: spinat(3,natom) 127 type(paw_ij_type),intent(inout) :: paw_ij(my_natom) 128 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom) 129 type(pawtab_type),intent(in) :: pawtab(ntypat) 130 131 !Local variables --------------------------------------- 132 !scalars 133 integer,parameter :: limp=0 ! could become an input variable 134 integer :: at_indx,cplex_dij,dmatudiag_loc,iafm,iatom,iatom_tot,iatpawu,icount 135 integer :: ilm,im1,im2,in1,in2,info,iplex,irot,ispden, irhoij,itypat,jlm,jrhoij 136 integer :: jspden,klmn,kspden,lcur,ldim,lmax,lmin,lpawu,lwork,my_comm_atom,ndij,nmat,nspden,nsploop 137 logical,parameter :: afm_noncoll=.true. ! TRUE if antiferro symmetries are used with non-collinear magnetism 138 logical :: antiferro,my_atmtab_allocated,noccsym_error,paral_atom,use_afm 139 real(dp),parameter :: invsqrt2=one/sqrt2 140 real(dp) :: factafm,mnorm,mx,my,mz,ntot,nup,ndn,snorm,sx,sy,szp,szm 141 character(len=4) :: wrt_mode 142 character(len=500) :: message 143 !arrays 144 integer :: nsym_used(2) 145 integer,pointer :: my_atmtab(:) 146 real(dp) :: sumocc(2) 147 real(dp),allocatable :: eig(:),hdp(:,:,:),hdp2(:,:),noccmmptemp(:,:,:,:),noccmmp_tmp(:,:,:,:) 148 real(dp),allocatable :: rwork(:),ro(:),noccmmp2(:,:,:,:),nocctot2(:) 149 complex(dpc),allocatable :: noccmmp_ylm(:,:,:),noccmmp_jmj(:,:),noccmmp_slm(:,:,:) 150 complex(dpc),allocatable :: zhdp(:,:),zhdp2(:,:),znoccmmp_tmp(:,:),zwork(:) 151 character(len=9),parameter :: dspin(6)= (/"up ","down ","up-up ","down-down","Re[up-dn]","Im[up-dn]"/) 152 character(len=9),parameter :: dspinc(6)= (/"up ","down ","up-up ","down-down","up-dn ","dn-up "/) 153 character(len=9),parameter :: dspinc2(6)=(/"up ","down ","dn-dn ","up-up ","dn-up ","up-dn "/) 154 character(len=9),parameter :: dspinm(6)= (/"dn ","up i ","n ","mx ","my ","mz "/) 155 type(coeff4_type),allocatable :: tmp_noccmmp(:) 156 157 !********************************************************************* 158 159 DBG_ENTER("COLL") 160 161 !Tests 162 if (my_natom>0) then 163 if (nsppol/=paw_ij(1)%nsppol) then 164 message=' inconsistent values for nsppol !' 165 MSG_BUG(message) 166 end if 167 if (compute_dmat>0) then 168 if (pawrhoij(1)%nspden/=paw_ij(1)%nspden.and.& 169 & pawrhoij(1)%nspden/=4.and.paw_ij(1)%nspden/=1) then 170 message=' inconsistent values for nspden !' 171 MSG_BUG(message) 172 end if 173 end if 174 end if 175 if (usepawu>0.and.useexexch>0) then 176 message=' usepawu>0 and useexexch>0 not allowed !' 177 MSG_BUG(message) 178 end if 179 if (impose_dmat/=0.and.dimdmat==0) then 180 message=' dmatpawu must be allocated when impose_dmat/=0 !' 181 MSG_BUG(message) 182 end if 183 if (usepawu>0.and.compute_dmat/=0.and.impose_dmat/=0.and.pawang%nsym==0) then 184 message=' pawang%zarot must be allocated !' 185 MSG_BUG(message) 186 end if 187 188 !Some inits 189 if (usepawu==0.and.useexexch==0) return 190 nspden=1;ndij=1;cplex_dij=1 191 if (my_natom>0) then 192 nspden=paw_ij(1)%nspden 193 ndij=paw_ij(1)%ndij 194 cplex_dij=paw_ij(1)%cplex_dij 195 end if 196 antiferro=(nspden==2.and.nsppol==1) 197 use_afm=((antiferro).or.((nspden==4).and.afm_noncoll)) 198 dmatudiag_loc=dmatudiag 199 if (dmatudiag==2.and.(dimdmat==0.or.impose_dmat==0)) dmatudiag_loc=1 200 201 !Set up parallelism over atoms 202 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 203 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 204 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 205 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) !vz_d 206 wrt_mode='COLL';if (paral_atom) wrt_mode='PERS' 207 208 !If needed, store dmatpu in suitable format in tmp_noccmmp 209 if (usepawu>0.and.impose_dmat/=0) then 210 iatpawu=0 211 ABI_DATATYPE_ALLOCATE(tmp_noccmmp,(natom)) 212 do iatom_tot=1,natom 213 itypat=typat(iatom_tot) 214 lpawu=pawtab(itypat)%lpawu 215 if (lpawu/=-1) then 216 iatpawu=iatpawu+1 217 if (ndij/=4) then 218 ABI_ALLOCATE(tmp_noccmmp(iatom_tot)%value,(cplex_dij,2*lpawu+1,2*lpawu+1,nsppol)) 219 tmp_noccmmp(iatom_tot)%value(1,1:2*lpawu+1,1:2*lpawu+1,1:nsppol)=& 220 & dmatpawu(1:2*lpawu+1,1:2*lpawu+1,1:nsppol,iatpawu) 221 else 222 ABI_ALLOCATE(tmp_noccmmp(iatom_tot)%value,(cplex_dij,2*lpawu+1,2*lpawu+1,ndij)) 223 tmp_noccmmp(iatom_tot)%value=zero 224 if(limp==0) then ! default reading 225 snorm=sqrt(spinat(1,iatom_tot)**2+spinat(1,iatom_tot)**2+spinat(3,iatom_tot)**2) 226 if (snorm>tol12) then 227 sx=half*spinat(1,iatom_tot)/snorm 228 sy=half*spinat(2,iatom_tot)/snorm 229 szp=half*(one+spinat(3,iatom_tot)/snorm) 230 szm=half*(one-spinat(3,iatom_tot)/snorm) 231 else 232 sx=zero;sy=zero 233 szp=one;szm=zero 234 end if 235 do im2=1,2*lpawu+1 236 do im1=1,2*lpawu+1 237 nup=dmatpawu(im1,im2,1,iatpawu);ndn=dmatpawu(im1,im2,2,iatpawu) 238 tmp_noccmmp(iatom_tot)%value(1,im1,im2,1)=nup*szp+ndn*szm 239 tmp_noccmmp(iatom_tot)%value(1,im1,im2,2)=nup*szm+ndn*szp 240 tmp_noccmmp(iatom_tot)%value(1,im1,im2,3)=(nup-ndn)*sx 241 tmp_noccmmp(iatom_tot)%value(1,im1,im2,4)=(ndn-nup)*sy 242 end do 243 end do 244 else if(limp>=1) then 245 ABI_ALLOCATE(noccmmp_ylm,(2*lpawu+1,2*lpawu+1,ndij)) 246 noccmmp_ylm=czero 247 ABI_ALLOCATE(noccmmp_slm,(2*lpawu+1,2*lpawu+1,ndij)) 248 noccmmp_slm=czero 249 ABI_ALLOCATE(noccmmp_jmj,(2*(2*lpawu+1),2*(2*lpawu+1))) 250 noccmmp_jmj=czero 251 if(limp==1) then ! read input matrix in J,M_J basis (l-1/2, then l+1/2) 252 noccmmp_jmj=czero 253 do im1=1,2*lpawu+1 254 noccmmp_jmj(im1,im1)=cmplx(dmatpawu(im1,im1,1,iatpawu),zero,kind=dp) 255 noccmmp_jmj(im1+lpawu,im1+lpawu)=cmplx(dmatpawu(im1+lpawu,im1+lpawu,2,iatpawu),zero,kind=dp) 256 end do 257 write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,& 258 & ' == Imposed occupation matrix (in the J M_J basis: L-1/2 and L+1/2 states)' 259 call wrtout(std_out,message,wrt_mode) 260 call mat_mlms2jmj(lpawu,noccmmp_ylm,noccmmp_jmj,ndij,& 261 & 2,2,pawprtvol,std_out,wrt_mode) ! optspin=1: up spin are first 262 end if 263 if(limp==2) then ! read input matrix in Ylm basis 264 noccmmp_ylm=czero 265 do im1=1,2*lpawu+1 266 noccmmp_ylm(im1,im1,1)=cmplx(dmatpawu(im1,im1,1,iatpawu),zero,kind=dp) 267 noccmmp_ylm(im1,im1,2)=cmplx(dmatpawu(im1,im1,2,iatpawu),zero,kind=dp) 268 end do 269 write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,& 270 & ' == Imposed occupation matrix (in the Ylm basis), for dn and up spin' 271 call wrtout(std_out,message,wrt_mode) 272 end if 273 call mat_slm2ylm(lpawu,noccmmp_ylm,noccmmp_slm,ndij,& 274 & 2,2,pawprtvol,std_out,wrt_mode) ! optspin=1 because up spin are first 275 ! interchange upup and dndn 276 if(limp>=1) then 277 tmp_noccmmp(iatom_tot)%value(1,:,:,1)=real(noccmmp_slm(:,:,2)) 278 tmp_noccmmp(iatom_tot)%value(2,:,:,1)=aimag(noccmmp_slm(:,:,2)) 279 tmp_noccmmp(iatom_tot)%value(1,:,:,2)=real(noccmmp_slm(:,:,1)) 280 tmp_noccmmp(iatom_tot)%value(2,:,:,2)=aimag(noccmmp_slm(:,:,1)) 281 tmp_noccmmp(iatom_tot)%value(1,:,:,3)=real(noccmmp_slm(:,:,4)) 282 tmp_noccmmp(iatom_tot)%value(2,:,:,3)=aimag(noccmmp_slm(:,:,4)) 283 tmp_noccmmp(iatom_tot)%value(1,:,:,4)=real(noccmmp_slm(:,:,3)) 284 tmp_noccmmp(iatom_tot)%value(2,:,:,4)=aimag(noccmmp_slm(:,:,3)) 285 end if 286 if(abs(pawprtvol)>2) then 287 write(message, '(2a)' ) ch10,& 288 & " Check Imposed density matrix in different basis" 289 call wrtout(std_out,message,wrt_mode) 290 call mat_slm2ylm(lpawu,noccmmp_slm,noccmmp_ylm,ndij,& 291 & 1,2,pawprtvol,std_out,wrt_mode) ! optspin=1 because up spin are first 292 call mat_mlms2jmj(lpawu,noccmmp_ylm,noccmmp_jmj,ndij,1,2,& 293 & pawprtvol,std_out,wrt_mode) ! optspin=1: up spin are first 294 end if 295 ABI_DEALLOCATE(noccmmp_ylm) 296 ABI_DEALLOCATE(noccmmp_jmj) 297 ABI_DEALLOCATE(noccmmp_slm) 298 end if 299 end if 300 end if 301 end do 302 end if ! impose_dmat/=0 303 304 !Print message 305 if (usepawu>0.and.impose_dmat/=0) then 306 if (dmatudiag_loc/=2) then 307 write(message,'(6a)') ch10,'Occupation matrix for correlated orbitals is kept constant',ch10,& 308 & 'and equal to dmatpawu from input file !',ch10,& 309 & '----------------------------------------------------------' 310 else 311 write(message,'(6a)') ch10,'Occupation matrix for correlated orbitals is imposed',ch10,& 312 & 'and equal to dmatpawu in the diagonal basis !',ch10,& 313 & '----------------------------------------------------------' 314 end if 315 call wrtout(std_out,message,'COLL') 316 end if 317 318 if (usepawu>0.and.dmatudiag_loc/=0.and.my_natom>0) then 319 write(message,'(4a)') ch10,'Diagonalized occupation matrix "noccmmp" is printed !',ch10,& 320 & '-------------------------------------------------------------' 321 call wrtout(std_out,message,wrt_mode) 322 end if 323 324 !Loops over atoms 325 do iatom=1,my_natom 326 iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom) 327 itypat=pawrhoij(iatom)%itypat 328 329 if (useexexch>0) then 330 lcur=pawtab(itypat)%lexexch 331 else if (usepawu>0) then 332 lcur=pawtab(itypat)%lpawu 333 end if 334 if (lcur/=-1) then 335 336 ! ######################################################################################## 337 ! # Compute nocc_mmp 338 ! ######################################################################################## 339 if ((usepawu>0.and.compute_dmat/=0).or.useexexch>0) then 340 341 342 paw_ij(iatom)%noccmmp(:,:,:,:)=zero 343 344 ! Loop over spin components 345 ABI_ALLOCATE(noccmmptemp,(cplex_dij,2*lcur+1,2*lcur+1,ndij)) 346 noccmmptemp(:,:,:,:)=zero 347 if(ndij==4) then 348 ABI_ALLOCATE(noccmmp2,(cplex_dij,2*lcur+1,2*lcur+1,ndij)) 349 end if 350 if(ndij==4) then 351 ABI_ALLOCATE(nocctot2,(ndij)) 352 end if 353 do ispden=1,ndij 354 jrhoij=1 355 do irhoij=1,pawrhoij(iatom)%nrhoijsel 356 klmn=pawrhoij(iatom)%rhoijselect(irhoij) 357 im1=pawtab(itypat)%klmntomn(1,klmn) 358 im2=pawtab(itypat)%klmntomn(2,klmn) 359 in1=pawtab(itypat)%klmntomn(3,klmn) 360 in2=pawtab(itypat)%klmntomn(4,klmn) 361 lmin=pawtab(itypat)%indklmn(3,klmn) 362 lmax=pawtab(itypat)%indklmn(4,klmn) 363 ABI_ALLOCATE(ro,(cplex_dij)) 364 if (ndij==1) then 365 ro(1)=half*pawrhoij(iatom)%rhoijp(jrhoij,1) 366 else if (ndij==2) then 367 ro(1)=pawrhoij(iatom)%rhoijp(jrhoij,ispden) 368 else ! ndij==4 369 ! Non-collinear magnetism: transfer rhoij to ro_c (keep n, m storage because 370 ! it is easier for the computation of noccmmp from rhoij) 371 ! cplex_dij has to be used here, because it is the dimension of rhoijp 372 ro(1:cplex_dij)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij-1+cplex_dij,ispden) 373 end if 374 if(lmin==0.and.lmax==2*lcur) then 375 icount=in1+(in2*(in2-1))/2 376 if(pawtab(itypat)%ij_proj<icount) then 377 message=' PAW+U: Problem in the loop for calculating noccmmp !' 378 MSG_BUG(message) 379 end if 380 if(in1/=in2) then 381 if(im2<=im1) then 382 noccmmptemp(:,im1,im2,ispden)=noccmmptemp(:,im1,im2,ispden)+ro(:)*pawtab(itypat)%phiphjint(icount) 383 end if 384 end if 385 if(im2>=im1) then 386 387 paw_ij(iatom)%noccmmp(:,im1,im2,ispden)=paw_ij(iatom)%noccmmp(:,im1,im2,ispden) & 388 & +ro(:)*pawtab(itypat)%phiphjint(icount) 389 end if 390 end if 391 jrhoij=jrhoij+pawrhoij(iatom)%cplex 392 ABI_DEALLOCATE(ro) 393 end do ! irhoij 394 do im2=1,2*lcur+1 395 do im1=1,im2 396 paw_ij(iatom)%noccmmp(1,im1,im2,ispden)=paw_ij(iatom)%noccmmp(1,im1,im2,ispden) & 397 & +noccmmptemp(1,im2,im1,ispden) 398 if(cplex_dij==2) paw_ij(iatom)%noccmmp(2,im1,im2,ispden)=paw_ij(iatom)%noccmmp(2,im1,im2,ispden) & 399 & -noccmmptemp(2,im2,im1,ispden) 400 end do 401 end do 402 do im1=1,2*lcur+1 403 do im2=1,im1 404 paw_ij(iatom)%noccmmp(1,im1,im2,ispden)=paw_ij(iatom)%noccmmp(1,im2,im1,ispden) 405 if(cplex_dij==2) paw_ij(iatom)%noccmmp(2,im1,im2,ispden)=-paw_ij(iatom)%noccmmp(2,im2,im1,ispden) 406 end do 407 end do 408 end do ! ispden 409 ABI_DEALLOCATE(noccmmptemp) 410 ! Compute noccmmp2, occupation matrix in the spin basis (upup, dndn, updn, dnup) 411 if(ndij==4) then 412 noccmmp2(:,:,:,:)=zero 413 do im1=1,2*lcur+1 414 do im2=1,2*lcur+1 415 noccmmp2(1,im1,im2,1)=half*(paw_ij(iatom)%noccmmp(1,im1,im2,1)+paw_ij(iatom)%noccmmp(1,im1,im2,4)) 416 noccmmp2(2,im1,im2,1)=half*(paw_ij(iatom)%noccmmp(2,im1,im2,1)+paw_ij(iatom)%noccmmp(2,im1,im2,4)) 417 noccmmp2(1,im1,im2,2)=half*(paw_ij(iatom)%noccmmp(1,im1,im2,1)-paw_ij(iatom)%noccmmp(1,im1,im2,4)) 418 noccmmp2(2,im1,im2,2)=half*(paw_ij(iatom)%noccmmp(2,im1,im2,1)-paw_ij(iatom)%noccmmp(2,im1,im2,4)) 419 noccmmp2(1,im1,im2,3)=half*(paw_ij(iatom)%noccmmp(1,im1,im2,2)+paw_ij(iatom)%noccmmp(2,im1,im2,3)) 420 noccmmp2(2,im1,im2,3)=half*(paw_ij(iatom)%noccmmp(2,im1,im2,2)-paw_ij(iatom)%noccmmp(1,im1,im2,3)) 421 noccmmp2(1,im1,im2,4)=half*(paw_ij(iatom)%noccmmp(1,im1,im2,2)-paw_ij(iatom)%noccmmp(2,im1,im2,3)) 422 noccmmp2(2,im1,im2,4)=half*(paw_ij(iatom)%noccmmp(2,im1,im2,2)+paw_ij(iatom)%noccmmp(1,im1,im2,3)) 423 end do 424 end do 425 if(abs(pawprtvol)>=1) then 426 write(message,'(2a)') ch10,"== Calculated occupation matrix for correlated orbitals in the n, m basis :" 427 call wrtout(std_out,message,wrt_mode) 428 do ispden=1,ndij 429 write(message,'(3a)') ch10,"Calculated occupation matrix for component ",trim(dspinm(ispden+2*(ndij/4))) 430 call wrtout(std_out,message,wrt_mode) 431 do im1=1,lcur*2+1 ! ( order of indices in noccmmp is exchanged in order to have the same convention as rhoij: transposition is done after ) 432 if(cplex_dij==1)& 433 & write(message,'(12(1x,9(1x,f10.5)))')& 434 & (paw_ij(iatom)%noccmmp(1,im2,im1,ispden),im2=1,lcur*2+1) 435 if(cplex_dij==2)& 436 ! & write(message,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))')& 437 & write(message,'(12(1x,9(1x,"(",f10.5,",",f10.5,")")))')& 438 & (paw_ij(iatom)%noccmmp(:,im2,im1,ispden),im2=1,lcur*2+1) 439 call wrtout(std_out,message,wrt_mode) 440 end do 441 end do 442 end if ! pawprtvol >=1 443 end if 444 445 ! Compute total number of electrons per spin 446 paw_ij(iatom)%nocctot(:)=zero ! contains nmmp in the n m representation 447 if(ndij==4) nocctot2(:)=zero ! contains nmmp in the upup dndn updn dnup representation 448 do ispden=1,ndij 449 do im1=1,2*lcur+1 450 if(ndij==4) then 451 paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)+paw_ij(iatom)%noccmmp(1,im1,im1,ispden) 452 nocctot2(ispden)=nocctot2(ispden)+noccmmp2(1,im1,im1,ispden) 453 else 454 paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)+paw_ij(iatom)%noccmmp(1,im1,im1,ispden) 455 end if 456 end do 457 end do 458 ! noccmmp will now be in the up up , dn dn... representation and now n_mmp=<m|n|mp> instead of <mp|n|m> ! 459 if(ndij==4) then 460 do ispden=1,ndij 461 do iplex=1,cplex_dij 462 do im1=1,2*lcur+1 463 do im2=1,2*lcur+1 464 paw_ij(iatom)%noccmmp(iplex,im1,im2,ispden)=noccmmp2(iplex,im2,im1,ispden) ! now, noccmmp is in the upup dndn updn dnup representation 465 end do 466 end do 467 end do 468 end do 469 ABI_DEALLOCATE(noccmmp2) 470 end if 471 ! Printing of new nocc_mmp 472 if ((usepawu>0.and.usepawu<10).or.(usepawu>=10.and.pawprtvol>=3)) & 473 write(message, '(2a)' ) ch10,'========== LDA+U DATA =================================================== ' 474 if (useexexch>0) write(message, '(2a)' )ch10,'======= Local ex-exchange (PBE0) DATA =================================== ' 475 if(((usepawu>0.and.usepawu<10).or.(usepawu>=10.and.pawprtvol>=3)).or.useexexch>0) call wrtout(std_out,message,wrt_mode) 476 if (usepawu>=10.and.pawprtvol>=3) then 477 write(message, '(6a)' ) ch10,' ( A DFT+DMFT calculation is carried out ',& 478 ch10,' Thus, the following LDA+U occupation matrices are not physical ',& 479 ch10,' and just informative )' 480 call wrtout(std_out,message,wrt_mode) 481 end if 482 if(usepawu<10.or.pawprtvol>=3) then ! Always write except if DMFT and pawprtvol low 483 write(message,'(2a,i5,a,i4,a)') ch10,"====== For Atom", iatom_tot,& 484 & ", occupations for correlated orbitals. l =",lcur,ch10 485 call wrtout(std_out,message,wrt_mode) 486 if(ndij==2) then 487 do ispden=1,2 488 write(message,'(a,i4,3a,f10.5)') "Atom", iatom_tot,". Occupations for spin ",& 489 & trim(dspin(ispden))," =",paw_ij(iatom)%nocctot(ispden) 490 call wrtout(std_out,message,wrt_mode) 491 end do 492 write(message,'(a,i4,a,2x,e16.8)') "=> On atom",iatom_tot,", local Mag. is ",& 493 & paw_ij(iatom)%nocctot(2)-paw_ij(iatom)%nocctot(1) 494 call wrtout(std_out,message,wrt_mode) 495 end if 496 if(ndij==4) then 497 ntot=paw_ij(iatom)%nocctot(1) 498 mx=paw_ij(iatom)%nocctot(2) 499 my=paw_ij(iatom)%nocctot(3) 500 mz=paw_ij(iatom)%nocctot(4) 501 mnorm=sqrt(mx*mx+my*my+mz*mz) 502 nup=nocctot2(1) 503 ndn=nocctot2(2) 504 write(message,'(a,i4,a,2x,e16.8)') "=> On atom",iatom_tot,", local Mag. x is ",mx 505 call wrtout(std_out,message,wrt_mode) 506 write(message,'(14x,a,2x,e16.8)') " local Mag. y is ",my 507 call wrtout(std_out,message,wrt_mode) 508 write(message,'(14x,a,2x,e16.8)') " local Mag. z is ",mz 509 call wrtout(std_out,message,wrt_mode) 510 write(message,'(14x,a,2x,e16.8)') " norm of Mag. is ",mnorm 511 call wrtout(std_out,message,wrt_mode) 512 write(message,'(14x,a,2x,f10.5)') " occ. of majority spin is ",half*(ntot+mnorm) ! to be checked versus direct calc from noccmmp 513 call wrtout(std_out,message,wrt_mode) 514 if(abs(pawprtvol)>=1) write(message,'(14x,a,2x,f10.5)') " occ. for spin up (along z) ",nup 515 if(abs(pawprtvol)>=1) then 516 call wrtout(std_out,message,wrt_mode) 517 end if 518 write(message,'(14x,a,2x,f10.5)') " occ. of minority spin is ",half*(ntot-mnorm) 519 call wrtout(std_out,message,wrt_mode) 520 if(abs(pawprtvol)>=1) write(message,'(14x,a,2x,f10.5)') " occ. for spin dn (along z) ",ndn 521 if(abs(pawprtvol)>=1) then 522 call wrtout(std_out,message,wrt_mode) 523 end if 524 if(ndij==4) then 525 ABI_DEALLOCATE(nocctot2) 526 end if 527 end if 528 write(message,'(2a)') ch10,"== Calculated occupation matrix for correlated orbitals:" 529 call wrtout(std_out,message,wrt_mode) 530 do ispden=1,ndij 531 write(message,'(3a)') ch10,"Calculated occupation matrix for component ",trim(dspinc(ispden+2*(ndij/4))) 532 call wrtout(std_out,message,wrt_mode) 533 do im1=1,lcur*2+1 534 if(cplex_dij==1)& 535 & write(message,'(12(1x,9(1x,f10.5)))')& 536 & (paw_ij(iatom)%noccmmp(1,im1,im2,ispden),im2=1,lcur*2+1) 537 if(cplex_dij==2)& 538 & write(message,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))')& 539 & (paw_ij(iatom)%noccmmp(:,im1,im2,ispden),im2=1,lcur*2+1) 540 call wrtout(std_out,message,wrt_mode) 541 end do 542 end do 543 end if 544 545 ! Transformation matrices: real->complex spherical harmonics (for test) 546 if(ndij==4.and.abs(pawprtvol)>=0) then 547 ABI_ALLOCATE(noccmmp_ylm,(2*lcur+1,2*lcur+1,ndij)) 548 noccmmp_ylm=czero 549 ABI_ALLOCATE(noccmmp_slm,(2*lcur+1,2*lcur+1,ndij)) 550 noccmmp_slm=czero 551 ABI_ALLOCATE(noccmmp_jmj,(2*(2*lcur+1),2*(2*lcur+1))) 552 noccmmp_jmj=czero 553 ! go from real notation for complex noccmmp to complex notation in noccmmp_slm 554 noccmmp_slm(:,:,:)=cmplx(paw_ij(iatom)%noccmmp(1,:,:,:)& 555 & ,paw_ij(iatom)%noccmmp(2,:,:,:),kind=dp) 556 call mat_slm2ylm(lcur,noccmmp_slm,noccmmp_ylm,ndij,1,1,pawprtvol,std_out,wrt_mode) ! optspin=1: up spin are first 557 558 do ispden=1,ndij 559 write(message,'(3a)') ch10,"Calculated Ylm occupation matrix for component ",trim(dspinc(ispden+2*(ndij/4))) 560 call wrtout(std_out,message,wrt_mode) 561 do im1=1,lcur*2+1 562 write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))') (noccmmp_ylm(im1,im2,ispden),im2=1,lcur*2+1) 563 call wrtout(std_out,message,wrt_mode) 564 end do 565 end do 566 call mat_mlms2jmj(lcur,noccmmp_ylm,noccmmp_jmj,ndij,1,1,pawprtvol,std_out,wrt_mode) ! optspin=1: up spin are first 567 ABI_DEALLOCATE(noccmmp_ylm) 568 ABI_DEALLOCATE(noccmmp_jmj) 569 ABI_DEALLOCATE(noccmmp_slm) 570 end if !ndij==4 571 572 end if ! impose_dmat==0 573 574 ! ######################################################################################## 575 ! # Diagonalize nocc_mmp 576 ! ######################################################################################## 577 if(usepawu>0.and.dmatudiag_loc>0) then 578 579 lpawu=lcur;ldim=2*lpawu+1 580 ABI_ALLOCATE(noccmmp_tmp,(1,ldim,ldim,ndij)) 581 if (ndij==4) then 582 ABI_ALLOCATE(znoccmmp_tmp,(2*ldim,2*ldim)) 583 end if 584 585 ! Select noccmmp for this atom 586 do ispden=1,ndij 587 noccmmp_tmp(1,:,:,ispden)=paw_ij(iatom)%noccmmp(1,:,:,ispden) 588 end do 589 if (ndij==4) then 590 do im2=1,ldim 591 do im1=1,ldim 592 znoccmmp_tmp(im1 , im2)=cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,1)& 593 & ,paw_ij(iatom)%noccmmp(2,im1,im2,1),kind=dp) 594 znoccmmp_tmp(ldim+im1,ldim+im2)=cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,2)& 595 & ,paw_ij(iatom)%noccmmp(2,im1,im2,2),kind=dp) 596 znoccmmp_tmp( im1,ldim+im2)=cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,3)& 597 & ,paw_ij(iatom)%noccmmp(2,im1,im2,3),kind=dp) 598 znoccmmp_tmp(ldim+im1, im2)=cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,4)& 599 & ,paw_ij(iatom)%noccmmp(2,im1,im2,4),kind=dp) 600 end do 601 end do 602 end if 603 604 ! Diagonalize nocc_mmp 605 if (ndij/=4) then 606 ABI_ALLOCATE(hdp,(ldim,ldim,ndij)) 607 hdp=zero 608 lwork=3*ldim-1 609 ABI_ALLOCATE(rwork,(lwork)) 610 ABI_ALLOCATE(eig,(ldim)) 611 do ispden=1,ndij 612 call dsyev('v','u',ldim,noccmmp_tmp(1,:,:,ispden),ldim,eig,rwork,lwork,info) 613 if(info/=0) then 614 message=' Error in diagonalization of noccmmp (DSYEV)!' 615 MSG_ERROR(message) 616 end if 617 do ilm=1,ldim 618 hdp(ilm,ilm,ispden)=eig(ilm) 619 end do 620 end do ! ispden 621 ABI_DEALLOCATE(rwork) 622 ABI_DEALLOCATE(eig) 623 else 624 ABI_ALLOCATE(hdp,(2*ldim,2*ldim,1)) 625 hdp=zero 626 lwork=4*ldim-1 627 ABI_ALLOCATE(rwork,(6*ldim-2)) 628 ABI_ALLOCATE(zwork,(lwork)) 629 ABI_ALLOCATE(eig,(2*ldim)) 630 call zheev('v','u',2*ldim,znoccmmp_tmp,2*ldim,eig,zwork,lwork,rwork,info) 631 if(info/=0) then 632 message=' Error in diagonalization of znoccmmp_tmp (zheev) !' 633 MSG_ERROR(message) 634 end if 635 do ilm=1,2*ldim 636 hdp(ilm,ilm,1)=eig(ilm) 637 end do 638 ABI_DEALLOCATE(rwork) 639 ABI_DEALLOCATE(zwork) 640 ABI_DEALLOCATE(eig) 641 end if 642 643 ! Print diagonalized matrix and eigenvectors 644 do ispden=1,size(hdp,3) 645 write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,' == Diagonalized Occupation matrix' 646 if (ndij==1) write(message,fmt='(2a)') trim(message)," for spin up ==" 647 if (ndij==2) write(message,fmt='(2a,i3,a)')trim(message)," for spin ",ispden," ==" 648 if (ndij==4) write(message,fmt='(2a,i3,a)')trim(message)," ==" 649 call wrtout(std_out,message,wrt_mode) 650 do ilm=1,size(hdp,1) 651 write(message,'(12(1x,9(1x,f10.5)))') (hdp(ilm,jlm,ispden),jlm=1,size(hdp,2)) 652 call wrtout(std_out,message,wrt_mode) 653 end do 654 end do ! ispden 655 if(abs(pawprtvol)>=1) then 656 if (ndij/=4) then 657 do ispden=1,ndij 658 write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,' == Eigenvectors' 659 if (ndij==1) write(message,fmt='(2a)') trim(message),' for spin up ==' 660 if (ndij==2) write(message,fmt='(2a,i3,a)')trim(message),' for spin ',ispden,' ==' 661 call wrtout(std_out,message,wrt_mode) 662 do ilm=1,ldim 663 write(message,'(12(1x,9(1x,f10.5)))') (noccmmp_tmp(1,ilm,jlm,ispden),jlm=1,ldim) 664 call wrtout(std_out,message,wrt_mode) 665 end do 666 end do 667 else 668 write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,' == Eigenvectors (spinors) in the real harmonics basis ==' 669 call wrtout(std_out,message,wrt_mode) 670 do ilm=1,2*ldim 671 write(message,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))') (znoccmmp_tmp(ilm,jlm),jlm=1,2*ldim) 672 call wrtout(std_out,message,wrt_mode) 673 end do 674 end if 675 end if 676 677 ! Back rotation of diagonalized matrix and printing 678 if(abs(pawprtvol)>=1) then 679 if (ndij/=4) then 680 ABI_ALLOCATE(hdp2,(ldim,ldim)) 681 do ispden=1,ndij 682 call dgemm('n','t',ldim,ldim,ldim,one,hdp(:,:,ispden),ldim,noccmmp_tmp(1,:,:,ispden),ldim,zero,hdp2,ldim) 683 call dgemm('n','n',ldim,ldim,ldim,one,noccmmp_tmp(1,:,:,ispden),ldim,hdp2,ldim,zero,hdp(:,:,ispden),ldim) 684 noccmmp_tmp(1,:,:,ispden)=hdp(:,:,ispden) 685 end do ! ispden 686 ABI_DEALLOCATE(hdp2) 687 else 688 ABI_ALLOCATE(zhdp,(2*ldim,2*ldim)) 689 ABI_ALLOCATE(zhdp2,(2*ldim,2*ldim)) 690 zhdp(:,:)=cmplx(hdp(:,:,1),zero,kind=dp) 691 zhdp2(:,:)=cmplx(zero,zero,kind=dp) 692 call zgemm('n','c',2*ldim,2*ldim,2*ldim,cone,zhdp,2*ldim,znoccmmp_tmp,2*ldim,czero,zhdp2,2*ldim) 693 zhdp(:,:)=cmplx(zero,zero,kind=dp) 694 call zgemm('n','n',2*ldim,2*ldim,2*ldim,cone,znoccmmp_tmp,2*ldim,zhdp2,2*ldim,czero,zhdp,2*ldim) 695 znoccmmp_tmp=zhdp 696 ABI_DEALLOCATE(zhdp) 697 ABI_DEALLOCATE(zhdp2) 698 end if 699 nmat=ndij ; if(ndij==4.and.cplex_dij==2) nmat=1 700 do ispden=1,nmat 701 write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,& 702 & ' == Rotated back diagonalized matrix' 703 if (ndij==1) write(message,fmt='(2a)') trim(message)," for spin up ==" 704 if (ndij==2) write(message,fmt='(2a,i3,a)')trim(message)," for spin ",ispden," ==" 705 if (ndij==4.and.cplex_dij==2) write(message,fmt='(4a)') trim(message)," for all component " 706 call wrtout(std_out,message,wrt_mode) 707 do ilm=1,ldim*cplex_dij 708 if(ndij==1.or.ndij==2)& 709 & write(message,'(12(1x,9(1x,f10.5)))')& 710 & (noccmmp_tmp(1,ilm,jlm,ispden),jlm=1,ldim) 711 if(ndij==4.and.cplex_dij==2)& 712 & write(message,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))')& 713 & (znoccmmp_tmp(ilm,jlm),jlm=1,ldim*cplex_dij) 714 call wrtout(std_out,message,wrt_mode) 715 end do 716 end do ! ispden 717 end if 718 ABI_DEALLOCATE(hdp) 719 720 end if ! dmatudiag_loc 721 722 ! ######################################################################################## 723 ! # Impose value of nocc_mmp from dmatpu; symetrize it 724 ! ######################################################################################## 725 if (usepawu>0.and.impose_dmat/=0) then 726 727 lpawu=lcur 728 nsploop=nsppol;if (ndij==4) nsploop=4 729 noccsym_error=.false. 730 731 ! Loop over spin components 732 do ispden=1,nsploop 733 if (ndij/=4) then 734 jspden=min(3-ispden,paw_ij(iatom)%nsppol) 735 else if (ispden<=2) then 736 jspden=3-ispden 737 else 738 jspden=ispden 739 end if 740 741 ! Loops over components of nocc_mmp 742 do jlm=1,2*lpawu+1 743 do ilm=1,2*lpawu+1 744 745 if(nsym>1.and.ndij<4) then 746 747 nsym_used(1:2)=0 748 sumocc(1:2)=zero 749 750 ! Accumulate values of nocc_mmp over symmetries 751 do irot=1,nsym 752 if ((symafm(irot)/=1).and.(.not.use_afm)) cycle 753 kspden=ispden;if (symafm(irot)==-1) kspden=jspden 754 factafm=one;if (ispden>3) factafm=dble(symafm(irot)) 755 iafm=1;if ((antiferro).and.(symafm(irot)==-1)) iafm=2 756 nsym_used(iafm)=nsym_used(iafm)+1 757 at_indx=indsym(4,irot,iatom_tot) 758 do im2=1,2*lpawu+1 759 do im1=1,2*lpawu+1 760 ! Be careful: use here R_rel^-1 in term of spherical harmonics 761 ! which is tR_rec in term of spherical harmonics 762 ! so, use transpose[zarot] 763 sumocc(iafm)=sumocc(iafm)+factafm*tmp_noccmmp(at_indx)%value(1,im1,im2,kspden) & 764 & *pawang%zarot(im1,ilm,lpawu+1,irot)& 765 & *pawang%zarot(im2,jlm,lpawu+1,irot) 766 ! sumocc(iafm)=sumocc(iafm)+factafm*tmp_noccmmp(at_indx)%value(im1,im2,kspden) & 767 ! & *pawang%zarot(ilm,im1,lpawu+1,irot)& 768 ! & *pawang%zarot(jlm,im2,lpawu+1,irot) 769 end do 770 end do 771 end do ! End loop over symmetries 772 773 ! Store new values of nocc_mmp 774 paw_ij(iatom)%noccmmp(1,ilm,jlm,ispden)=sumocc(1)/nsym_used(1) 775 if (.not.noccsym_error)& 776 & noccsym_error=(abs(paw_ij(iatom)%noccmmp(1,ilm,jlm,ispden) & 777 & -tmp_noccmmp(iatom_tot)%value(1,ilm,jlm,ispden))>tol5) 778 779 ! Antiferromagnetic case: has to fill up "down" component of nocc_mmp 780 if (antiferro.and.nsym_used(2)>0) paw_ij(iatom)%noccmmp(1,ilm,jlm,2)=sumocc(2)/nsym_used(2) 781 782 else ! nsym=1 783 784 ! Case without symetries 785 paw_ij(iatom)%noccmmp(:,ilm,jlm,ispden)= tmp_noccmmp(iatom_tot)%value(:,ilm,jlm,ispden) 786 end if 787 788 end do !ilm 789 end do !jlm 790 end do ! ispden 791 do ispden=1,nsploop 792 paw_ij(iatom)%nocctot(ispden)=zero ! contains nmmp in the n m representation 793 do im1=1,2*lcur+1 794 if(ndij==4.and.ispden==1) then 795 ! in this case, on computes total number or electron for double counting correction 796 paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)+& 797 & paw_ij(iatom)%noccmmp(1,im1,im1,1)+paw_ij(iatom)%noccmmp(1,im1,im1,2) 798 else if(ndij==4.and.ispden==2) then 799 paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)+& 800 & paw_ij(iatom)%noccmmp(1,im1,im1,3)+paw_ij(iatom)%noccmmp(1,im1,im1,4) 801 else if(ndij==4.and.ispden==3) then 802 paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)-& 803 & paw_ij(iatom)%noccmmp(2,im1,im1,3)+paw_ij(iatom)%noccmmp(2,im1,im1,4) 804 else if(ndij==4.and.ispden==4) then 805 paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)+& 806 & paw_ij(iatom)%noccmmp(2,im1,im1,1)-paw_ij(iatom)%noccmmp(2,im1,im1,2) 807 else 808 paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)+& 809 & paw_ij(iatom)%noccmmp(1,im1,im1,ispden) 810 end if 811 end do 812 end do ! ispden 813 814 ! Printing of new nocc_mmp 815 do ispden=1,ndij 816 if(dmatudiag_loc==2) then 817 write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,& 818 & ' == Imposed occupation matrix (in the basis of diagonalization!!)' 819 else 820 write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,& 821 & ' == Imposed occupation matrix' 822 end if 823 if (ndij==1) write(message,fmt='(2a)') trim(message)," for spin up ==" 824 if (ndij==2) write(message,fmt='(2a,i3,a)')trim(message)," for spin ",ispden," ==" 825 if (ndij==4) write(message,fmt='(4a)') trim(message)," for component ", & 826 & trim(dspinc(ispden+2*(ndij/4)))," ==" 827 call wrtout(std_out,message,wrt_mode) 828 do ilm=1,2*lpawu+1 829 if(cplex_dij==1)& 830 & write(message,'(12(1x,9(1x,f10.5)))')& 831 & (paw_ij(iatom)%noccmmp(1,ilm,jlm,ispden),jlm=1,2*lpawu+1) 832 if(cplex_dij==2)& 833 & write(message,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))')& 834 & (paw_ij(iatom)%noccmmp(:,ilm,jlm,ispden),jlm=1,2*lpawu+1) 835 call wrtout(std_out,message,wrt_mode) 836 end do 837 end do 838 839 ! WARNING if symmetrization changes the matrix 840 if (noccsym_error) then 841 write(message, '(a,i4,6a)' ) & 842 ' After symmetrization, imposed occupation matrix for atom ',iatom_tot,ch10,& 843 & ' is different from dmatpawu value set in input file !',ch10,& 844 & ' It is likely that dmatpawu does not match the symmetry operations of the system.',ch10,& 845 & ' Action: change dmatpawu in input file or increase precision until 0.00001' 846 MSG_WARNING(message) 847 end if 848 849 end if ! impose_dmat/=0 850 851 ! ######################################################################################## 852 ! # Rotate imposed occupation matrix in the non-diagonal basis 853 ! ######################################################################################## 854 if (usepawu>0.and.impose_dmat/=0.and.dmatudiag_loc==2) then 855 856 lpawu=lcur;ldim=2*lpawu+1 857 858 ! Rotation of imposed nocc_mmp 859 if (ndij/=4) then 860 ABI_ALLOCATE(hdp2,(ldim,ldim)) 861 do ispden=1,ndij 862 call dgemm('n','t',ldim,ldim,ldim,one,& 863 & paw_ij(iatom)%noccmmp(1,:,:,ispden),ldim,noccmmp_tmp(1,:,:,ispden),ldim,zero,hdp2,ldim) 864 call dgemm('n','n',ldim,ldim,ldim,one,& 865 & noccmmp_tmp(1,:,:,ispden),ldim,hdp2,ldim,zero,paw_ij(iatom)%noccmmp(1,:,:,ispden),ldim) 866 end do ! ispden 867 ABI_DEALLOCATE(hdp2) 868 else 869 ABI_ALLOCATE(zhdp,(2*ldim,2*ldim)) 870 ABI_ALLOCATE(zhdp2,(2*ldim,2*ldim)) 871 do im2=1,ldim 872 do im1=1,ldim 873 zhdp( im1, im2)=cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,1),zero,kind=dp) ! to be checked 874 zhdp(ldim+im1,ldim+im2)=cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,2),zero,kind=dp) ! to be checked 875 zhdp( im1,ldim+im2)=& 876 & cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,3),+paw_ij(iatom)%noccmmp(1,im1,im2,4),kind=dp) ! to be checked 877 zhdp(ldim+im1, im2)=& 878 & cmplx(paw_ij(iatom)%noccmmp(1,im2,im1,3),-paw_ij(iatom)%noccmmp(1,im2,im1,4),kind=dp) ! to be checked 879 end do 880 end do 881 call zgemm('n','c',2*ldim,2*ldim,2*ldim,cone,zhdp,2*ldim,znoccmmp_tmp,2*ldim,czero,zhdp2,2*ldim) 882 call zgemm('n','n',2*ldim,2*ldim,2*ldim,cone,znoccmmp_tmp,2*ldim,zhdp2,2*ldim,czero,zhdp,2*ldim) 883 do jlm=1,ldim 884 do ilm=1,ldim 885 paw_ij(iatom)%noccmmp(1,ilm,jlm,1)= real(znoccmmp_tmp( ilm, jlm)) ! to be checked 886 paw_ij(iatom)%noccmmp(1,ilm,jlm,2)= real(znoccmmp_tmp(ldim+ilm,ldim+jlm)) ! to be checked 887 paw_ij(iatom)%noccmmp(1,ilm,jlm,3)= real(znoccmmp_tmp( ilm,ldim+jlm)) ! to be checked 888 paw_ij(iatom)%noccmmp(1,ilm,jlm,4)=aimag(znoccmmp_tmp( ilm,ldim+jlm)) ! to be checked 889 end do 890 end do 891 ABI_DEALLOCATE(zhdp) 892 ABI_DEALLOCATE(zhdp2) 893 end if 894 895 ! Printing of rotated imposed matrix 896 do ispden=1,ndij 897 write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,& 898 & ' == Imposed density matrix in original basis' 899 if (ndij==1) write(message,fmt='(2a)') trim(message)," for spin up ==" 900 if (ndij==2) write(message,fmt='(2a,i3,a)')trim(message)," for spin ",ispden," ==" 901 if (ndij==4) write(message,fmt='(4a)') trim(message)," for component ", & 902 & trim(dspin(ispden+2*(ndij/4)))," ==" 903 call wrtout(std_out,message,wrt_mode) 904 do ilm=1,2*lpawu+1 905 write(message,'(12(1x,9(1x,f10.5)))') (paw_ij(iatom)%noccmmp(1,ilm,jlm,ispden),jlm=1,2*lpawu+1) ! to be checked 906 call wrtout(std_out,message,wrt_mode) 907 end do 908 end do ! ispden 909 910 end if ! dmatudiag_loc==2 911 912 if (usepawu>0.and.dmatudiag_loc>0) then 913 ABI_DEALLOCATE(noccmmp_tmp) 914 if (ndij==4) then 915 ABI_DEALLOCATE(znoccmmp_tmp) 916 end if 917 end if 918 919 paw_ij(iatom)%has_pawu_occ=2 920 921 end if ! lcur 922 end do ! iatom 923 924 !Memory deallocation 925 if (usepawu>0.and.impose_dmat/=0) then 926 do iatom_tot=1,natom 927 lpawu=pawtab(typat(iatom_tot))%lpawu 928 if (lpawu/=-1) then 929 ABI_DEALLOCATE(tmp_noccmmp(iatom_tot)%value) 930 end if 931 end do 932 ABI_DATATYPE_DEALLOCATE(tmp_noccmmp) 933 end if 934 935 !Destroy atom table used for parallelism 936 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 937 938 DBG_EXIT("COLL") 939 940 end subroutine setnoccmmp