TABLE OF CONTENTS
- ABINIT/m_paw_uj
- m_paw_uj/blow_pawuj
- m_paw_uj/chiscwrt
- m_paw_uj/lcalcu
- m_paw_uj/linvmat
- m_paw_uj/lprtmat
- m_paw_uj/macro_uj_type
- m_paw_uj/pawuj_det
- m_paw_uj/pawuj_free
- m_paw_uj/pawuj_ini
- m_paw_uj/pawuj_red
ABINIT/m_paw_uj [ Modules ]
NAME
m_paw_uj
FUNCTION
This module contains several routines relevant only for automatic determination of U in PAW+U context (linear response method according to Phys. Rev. B 71, 035105)
COPYRIGHT
Copyright (C) 2018-2022 ABINIT group (DJA) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 MODULE m_paw_uj 24 25 use defs_basis 26 use m_abicore 27 use m_errors 28 use m_linalg_interfaces 29 use m_xmpi 30 use m_dtset 31 use netcdf 32 use m_nctk 33 34 use m_fstrings, only : strcat 35 use m_pptools, only : prmat 36 use m_special_funcs, only : iradfnh 37 use m_geometry, only : shellstruct,ioniondist 38 use m_parser, only : prttagm 39 use m_supercell, only : mksupercell 40 use m_pawrad, only : pawrad_type 41 use m_pawtab, only : pawtab_type 42 use m_paw_ij, only : paw_ij_type 43 use m_paral_atom, only : get_my_atmtab, free_my_atmtab 44 use m_dtfil, only : datafiles_type 45 use m_crystal, only : crystal_t 46 47 implicit none 48 49 private
m_paw_uj/blow_pawuj [ Functions ]
[ Top ] [ m_paw_uj ] [ Functions ]
NAME
blow_pawuj
FUNCTION
This subroutine reads a real nxn matrice and appends lines n+1 and clumn n+1 containing the sum of the lines
INPUTS
mat(nj,nj) matrix to be completed
OUTPUT
matt(nj+1,nj+1) completed matrix
SOURCE
1424 subroutine blow_pawuj(mat,nj,matt) 1425 1426 !Arguments ------------------------------------ 1427 !scalars 1428 integer,intent(in) :: nj 1429 !arrays 1430 real(dp),intent(in) :: mat(nj,nj) 1431 real(dp),intent(out) :: matt(nj+1,nj+1) 1432 1433 !Local variables------------------------------- 1434 !scalars 1435 integer :: ii 1436 !arrays 1437 1438 ! ************************************************************************* 1439 1440 matt(1:nj,1:nj)=mat 1441 do ii = 1,nj 1442 matt(ii,nj+1)=-sum(matt(ii,1:nj)) 1443 end do 1444 1445 do ii = 1,nj+1 1446 matt(nj+1,ii)=-sum(matt(1:nj,ii)) 1447 end do 1448 1449 end subroutine blow_pawuj
m_paw_uj/chiscwrt [ Functions ]
[ Top ] [ m_paw_uj ] [ Functions ]
NAME
chiscwrt
FUNCTION
PRIVATE ROUTINE Distributes values of chi_org on chi_sc according to ion-ion distances and multiplicities in shells
INPUTS
chi_org = response matrix in original unit cell disv_org = distances (multiplied by magntization) in original unit cell nat_org = number of atoms in original unit cell sdisv_org = radii of shells in original unit cell smult_org = multiplicities of shells in original unit cell nsh_org = number of shells in original unit cell disv_sc = distances (multiplied by magntization) in super-cell nat_sc = number of atoms in super-cell sdisv_sc = radii of shells in super-cell (was unused, so suppressed) smult_sc = multiplicities of shells in super-cell nsh_sc = number of shells in super-cell opt =
OUTPUT
chi_sc = first line of reponse matrix in super-cell
SIDE EFFECTS
SOURCE
1055 subroutine chiscwrt(chi_org,disv_org,nat_org,sdisv_org,smult_org,nsh_org,chi_sc,& 1056 & disv_sc,nat_sc,smult_sc,nsh_sc,opt,prtvol) 1057 1058 !Arguments ------------------------------------ 1059 !scalars 1060 integer,intent(in) :: nat_org,nsh_org,nsh_sc 1061 integer,intent(in) :: nat_sc 1062 integer,intent(in),optional :: prtvol 1063 integer,intent(in),optional :: opt 1064 !arrays 1065 real(dp),intent(in) :: chi_org(nat_org),disv_org(nat_org),sdisv_org(nsh_org) 1066 integer,intent(in) :: smult_org(nsh_org), smult_sc(nsh_sc) 1067 real(dp),intent(in) :: disv_sc(nat_sc) 1068 real(dp),intent(out) :: chi_sc(nat_sc) 1069 1070 !Local variables------------------------------- 1071 !scalars 1072 integer :: iatom,jatom,jsh,optt 1073 character(len=500) :: message 1074 !arrays 1075 real(dp) :: chi_orgl(nat_org) 1076 1077 ! ************************************************************************* 1078 1079 if (present(opt)) then 1080 optt=opt 1081 else 1082 optt=1 1083 end if 1084 1085 1086 do iatom=1,nat_org 1087 do jsh=1,nsh_org 1088 if (disv_org(iatom)==sdisv_org(jsh)) then 1089 if (opt==2) then 1090 chi_orgl(iatom)=chi_org(iatom) 1091 else if (opt==1) then 1092 chi_orgl(iatom)=chi_org(iatom)*smult_org(jsh)/smult_sc(jsh) 1093 end if 1094 exit 1095 end if 1096 end do !iatom 1097 end do !jsh 1098 1099 if (prtvol>1) then 1100 write(message,fmt='(a,150f10.5)')' chiscwrt: chi at input ',chi_org 1101 call wrtout(std_out,message,'COLL') 1102 write(message,fmt='(a,150f10.5)')' chiscwrt: chi after division ',chi_orgl 1103 call wrtout(std_out,message,'COLL') 1104 end if 1105 1106 do iatom=1,nat_sc 1107 do jatom=1,nat_org 1108 if (disv_org(jatom)==disv_sc(iatom)) then 1109 chi_sc(iatom)=chi_orgl(jatom) 1110 exit 1111 else if (jatom==nat_org) then 1112 chi_sc(iatom)=0_dp 1113 end if 1114 end do 1115 end do 1116 1117 if (prtvol>1) then 1118 write(message,'(a)')' chiscwrt, chi_sc ' 1119 call wrtout(std_out,message,'COLL') 1120 call prmat(chi_sc,1,nat_sc,1,std_out) 1121 end if 1122 1123 end subroutine chiscwrt
m_paw_uj/lcalcu [ Functions ]
[ Top ] [ m_paw_uj ] [ Functions ]
NAME
lcalcu
FUNCTION
prints out real the real matrice mmat
INPUTS
magv=magnetic ordering of the ions (-1 of down/1 for up) natom=number of atoms rprimd(3,3)=lattic vectors of unit cell xred(3,natom)=positions of atoms chi(natom)=full response of atoms due to shift on atom pawujat chi0(natom)= response of atoms due to shift on atom pawujat pawujat=specifies on which atom the potential shift was done prtvol=controls output to files (see subroutine lprtmat) gam=gamma to be used for inversion of matrices (see subroutine livmat) opt=wether to use charge bath (1 or 3) or not (else)
OUTPUT
ures=resulting U (in eV) on atom pawujat
SOURCE
1323 subroutine lcalcu(magv,natom,rprimd,xred,chi,chi0,pawujat,ures,prtvol,gam,opt) 1324 1325 !Arguments ------------------------------------ 1326 !scalars 1327 integer,intent(in) :: natom 1328 integer,intent(in),optional :: opt,pawujat,prtvol 1329 real(dp),intent(in),optional:: gam 1330 real(dp),intent(out) :: ures 1331 !arrays 1332 integer,intent(in) :: magv(natom) 1333 real(dp),intent(in) :: rprimd(3,3),xred(3,natom),chi(natom),chi0(natom) 1334 1335 !Local variables------------------------------- 1336 !scalars 1337 character(len=500) :: message 1338 integer :: optt,prtvoll,nnatom 1339 real(dp) :: gamm 1340 !arrays 1341 real(dp),allocatable :: chi0matrix1(:,:),chimatrix2(:,:),auginvchi0matrix3(:,:) 1342 real(dp),allocatable :: auginvchimatrix4(:,:),uresmatrix(:,:) 1343 1344 ! ********************************************************************* 1345 1346 if (present(opt)) then 1347 optt=opt 1348 else 1349 optt=1 1350 end if 1351 1352 if (present(prtvol)) then 1353 prtvoll=prtvol 1354 else 1355 prtvoll=1 1356 end if 1357 1358 if (present(gam)) then 1359 gamm=gam 1360 else 1361 gamm=1_dp 1362 end if 1363 1364 if (optt==1.or.optt==3) then 1365 nnatom=natom+1 1366 else 1367 nnatom=natom 1368 end if 1369 1370 ABI_MALLOC(chi0matrix1,(natom,natom)) 1371 ABI_MALLOC(chimatrix2,(natom,natom)) 1372 ABI_MALLOC(auginvchi0matrix3,(nnatom,nnatom)) 1373 ABI_MALLOC(auginvchimatrix4,(nnatom,nnatom)) 1374 ABI_MALLOC(uresmatrix,(nnatom,nnatom)) 1375 1376 call ioniondist(natom,rprimd,xred,chi0matrix1,3,chi0,magv,pawujat,prtvoll) 1377 call ioniondist(natom,rprimd,xred,chimatrix2,3,chi,magv,pawujat,prtvoll) 1378 1379 write(message,fmt='(a)')'response chi_0' 1380 call linvmat(chi0matrix1,auginvchi0matrix3,natom,message,optt,gamm,prtvoll) 1381 call wrtout(std_out,message,'COLL') 1382 1383 write(message,fmt='(a)')'response chi' 1384 call linvmat(chimatrix2,auginvchimatrix4,natom,message,optt,gamm,prtvoll) 1385 call wrtout(std_out,message,'COLL') 1386 1387 uresmatrix=(auginvchi0matrix3-auginvchimatrix4)*Ha_eV 1388 1389 write(message,fmt='(a,i3,a)')' (chi_0)^(-1)-(chi)^(-1) (eV)' 1390 call lprtmat(message,2,prtvoll,uresmatrix,nnatom) 1391 call wrtout(std_out,message,'COLL') 1392 1393 ures=uresmatrix(1,pawujat) 1394 1395 ABI_FREE(chi0matrix1) 1396 ABI_FREE(chimatrix2) 1397 ABI_FREE(auginvchi0matrix3) 1398 ABI_FREE(auginvchimatrix4) 1399 ABI_FREE(uresmatrix) 1400 1401 1402 end subroutine lcalcu
m_paw_uj/linvmat [ Functions ]
[ Top ] [ m_paw_uj ] [ Functions ]
NAME
linvmat
FUNCTION
inverts real matrix inmat
INPUTS
inmat(1:nat,1:nat)=matrix to be inverted nat=dimension of inmat nam=comment specifiying the input matrix (to be printed in output) option=how to invert inmat option=1 or 3 add charge bath to matrix and add gam for inversion option=2 simply invert matrix gam=gamma added to inmat before inversion in case charge bath is used (allows inversion of otherwise singular matrix) prtvol=controls output to files (see subroutine lprtmat)
OUTPUT
oumat(nnat,nnat)=inverse of inmat, nnat=nat+1 for option=1 or option=3; nnat=nat for option=2
SOURCE
1151 subroutine linvmat(inmat,oumat,nat,nam,option,gam,prtvol) 1152 1153 !Arguments ------------------------------- 1154 1155 integer,intent(in) :: nat 1156 real(dp),intent(in) :: gam,inmat(nat,nat) 1157 real(dp),intent(inout) :: oumat(:,:) ! nat+1,nat+1 for option=1 or 2 1158 ! nat,nat for option=2 1159 character(len=500),intent(in) :: nam 1160 integer,intent(in),optional :: prtvol,option 1161 1162 !Local variables ------------------------- 1163 character(len=500) :: message 1164 character(len=500) :: bastrin,gastrin 1165 integer :: info,nnat,optionn !,ii,jj 1166 integer,allocatable :: ipvt(:) 1167 real(dp),allocatable :: hma(:,:),work(:) 1168 1169 ! ********************************************************************* 1170 1171 if (present(option)) then 1172 optionn=option 1173 else 1174 optionn=1 1175 end if 1176 1177 if (option==1.or.option==3) then 1178 write(bastrin,'(a)')'+ charge bath ' 1179 write(gastrin,'(a,d10.2,a)')'+ gamma (=',gam,') ' 1180 else 1181 write(bastrin,'(a)')'' 1182 write(gastrin,'(a)')'' 1183 end if 1184 1185 write(message,fmt='(a)')' matrix '//trim(nam) 1186 call lprtmat(message,1,prtvol,inmat,nat) 1187 1188 if (option==1.or.option==3) then 1189 call blow_pawuj(inmat,nat,oumat) 1190 write(message,fmt='(a,a)')' ',trim(nam)//trim(bastrin) 1191 call lprtmat(message,1,prtvol,oumat,nat+1) 1192 oumat=oumat+gam 1193 nnat=nat+1 1194 else 1195 nnat=nat 1196 oumat=inmat 1197 oumat(1,1)=inmat(1,1) 1198 end if 1199 1200 1201 ABI_MALLOC(hma,(nnat,nnat)) 1202 ABI_MALLOC(work,(nnat)) 1203 ABI_MALLOC(ipvt,(nnat)) 1204 work=0_dp 1205 hma(:,:)=oumat 1206 1207 call dgetrf(nnat,nnat,hma,nnat,ipvt,info) 1208 if (.not.info==0) then 1209 write(message, '(3a)' ) 'Matrix '//trim(nam)//' is singular',ch10,'Probably too many symmetries kept' 1210 call wrtout(ab_out,message,'COLL') 1211 return 1212 end if 1213 1214 call dgetri(nnat,hma,nnat,ipvt,work,nnat,info) 1215 oumat=hma(:,:) 1216 1217 write(message,fmt='(2a,a)')' ('//trim(nam)//trim(bastrin)//trim(gastrin)//')^(-1)' 1218 call lprtmat(message,1,prtvol,oumat,nnat) 1219 1220 ABI_FREE(hma) 1221 ABI_FREE(work) 1222 ABI_FREE(ipvt) 1223 1224 end subroutine linvmat
m_paw_uj/lprtmat [ Functions ]
[ Top ] [ m_paw_uj ] [ Functions ]
NAME
lprtmat
FUNCTION
prints out the real matrix mmat
INPUTS
mmat(nat,nat)=matrix to be printed nat=dimension of mmat prtvol specifies the volume of printing 3: print the whole matrix 2: print the first line 1: do not print anything chan specifies the output files 1: output only to std_out 2: output also to ab_out commnt=comment specifying matirix
OUTPUT
oumat(nat+1,nat+1)=inverse of inmat
SOURCE
1253 subroutine lprtmat(commnt,chan,prtvol,mmat,nat) 1254 1255 !Arguments ------------------------------- 1256 integer,intent(in) :: nat,chan,prtvol 1257 real(dp),intent(in) :: mmat(nat,nat) 1258 character(len=500),intent(in) :: commnt 1259 1260 !Local variables ------------------------- 1261 character(len=500) :: message 1262 ! ********************************************************************* 1263 1264 if (prtvol==3) then 1265 write(message,fmt='(a)') trim(commnt) 1266 call wrtout(std_out,message,'COLL') 1267 call prmat(mmat,nat,nat,nat,std_out) 1268 if (chan==2) then 1269 call wrtout(ab_out,message,'COLL') 1270 call prmat(mmat,nat,nat,nat,ab_out) 1271 end if 1272 write(message,*)ch10 1273 call wrtout(std_out,message,'COLL') 1274 if (chan==2) then 1275 call wrtout(ab_out,message,'COLL') 1276 end if 1277 end if 1278 1279 if (prtvol==2) then 1280 write(message,fmt='(a)') trim(commnt) 1281 call wrtout(std_out,message,'COLL') 1282 call prmat(mmat,1,nat,nat,std_out) 1283 if (chan==2) then 1284 call wrtout(ab_out,message,'COLL') 1285 call prmat(mmat,1,nat,nat,ab_out) 1286 end if 1287 write(message,*)ch10 1288 call wrtout(std_out,message,'COLL') 1289 if (chan==2) then 1290 call wrtout(ab_out,message,'COLL') 1291 end if 1292 end if 1293 1294 end subroutine lprtmat
m_paw_uj/macro_uj_type [ Types ]
[ Top ] [ m_paw_uj ] [ Types ]
NAME
dtmacro_uj
FUNCTION
This data type contains the potential shifts and the occupations for the determination of U and J for the DFT+U calculations. iuj=1,3: non-selfconsistent calculations. iuj=2,4 selfconsistent calculations.
SOURCE
65 type, public :: macro_uj_type 66 67 ! Integer 68 integer :: iuj ! dataset treated 69 integer :: nat ! number of atoms U (J) is determined on 70 integer :: ndtset ! total number of datasets 71 integer :: nspden ! number of densities treated 72 integer :: macro_uj ! which mode the determination runs in 73 integer :: pawujat ! which atom U (J) is determined on 74 integer :: pawprtvol ! controlling amount of output 75 integer :: option ! controls the determination of U (1 with compensating charge bath, 2 without) 76 integer :: dmatpuopt ! controls the renormalisation of the PAW projectors 77 78 ! Real 79 real(dp) :: diemix ! mixing parameter 80 real(dp) :: diemixmag ! magnetic mixing parameter 81 real(dp) :: mdist ! maximal distance of ions 82 real(dp) :: pawujga ! gamma for inversion of singular matrices 83 real(dp) :: ph0phiint ! integral of phi(:,1)*phi(:,1) 84 real(dp) :: pawujrad ! radius to which U should be extrapolated. 85 real(dp) :: pawrad ! radius of the paw atomic data 86 87 ! Integer arrays 88 integer , allocatable :: scdim(:) 89 ! size of supercell 90 91 ! Real arrays 92 real(dp) , allocatable :: occ(:,:) 93 ! occupancies after a potential shift: occ(ispden,nat) 94 95 real(dp) , allocatable :: rprimd(:,:) 96 ! unit cell for symmetrization 97 98 real(dp) , allocatable :: vsh(:,:) 99 ! potential shifts on atoms, dimensions: nspden,nat 100 101 real(dp) , allocatable :: xred(:,:) 102 ! atomic position for symmetrization 103 104 real(dp) , allocatable :: wfchr(:) 105 ! wfchr(1:3): zion, n and l of atom on which projection is done 106 ! wfchr(4:6): coefficients ai of a0+a1*r+a2*r^2, fit to the wfc for r< r_paw 107 108 real(dp), allocatable :: zioneff(:) 109 ! zioneff(ij_proj), "Effective charge"*n "seen" at r_paw, deduced from Phi at r_paw, n: 110 ! pricipal quantum number; good approximation to model wave function outside PAW-sphere through 111 112 end type macro_uj_type
m_paw_uj/pawuj_det [ Functions ]
[ Top ] [ m_paw_uj ] [ Functions ]
NAME
pawuj_det
FUNCTION
From the complete dtpawuj-dataset determines U (or J) parameter for PAW+U calculations Relevant only for automatic determination of Hubbard Parameters in PAW+U context. Hubbard U = diemix/chi0-1/chi Hund's J = 1/chi-diemixmag/chi0
INPUTS
dtpawuj=potential shifts (vsh) and atomic occupations (occ)
OUTPUT
only printing (among other things a section in the ab.out that can be used for input in ujdet)
SOURCE
269 subroutine pawuj_det(dtpawuj,ndtpawuj,dtset,dtfil,ures,comm) 270 271 !Arguments ------------------------------------ 272 !scalars 273 !arrays 274 integer :: ndtpawuj, comm 275 type(macro_uj_type),intent(in) :: dtpawuj(0:ndtpawuj) 276 type(dataset_type),intent(in) :: dtset 277 type(datafiles_type),intent(in) :: dtfil 278 real(dp),intent(out) :: ures 279 280 !Local variables------------------------------- 281 !scalars 282 integer,parameter :: natmax=2,nwfchr=6 283 integer :: ii,jj,nat_org,natom,jdtset,nspden,macro_uj,marr 284 integer :: im1,ndtuj,nsh_org, nsh_sc,nat_sc,maxnat 285 integer :: pawujat,pawprtvol,pawujoption 286 integer :: dmatpuopt,invopt,ipert 287 integer :: my_rank, ncid, ncerr 288 real(dp) :: pawujga,ph0phiint,intg,fcorr,eyp 289 real(dp) :: diem,signum,scalarHP !LMac quantities 290 291 character(len=500) :: message 292 !character(len=2) :: hstr 293 character(len=5) :: pertname 294 character(len=1) :: parname 295 character(len=14) :: occmag 296 type(crystal_t) :: cryst 297 !arrays 298 integer :: ext(3) 299 real(dp) :: rprimd_sc(3,3),vsh(ndtpawuj),a(5),b(5) 300 integer,allocatable :: jdtset_(:),smult_org(:),smult_sc(:),intmagv_org(:),intmagv_sc(:) 301 real(dp),allocatable :: luocc(:,:),xred_org(:,:) 302 real(dp),allocatable :: magv_org(:),magv_sc(:),chi(:),chi0(:),chi0_sc(:), chi_sc(:), xred_sc(:,:) 303 real(dp),allocatable :: sdistv_org(:),sdistv_sc(:),distv_org(:),distv_sc(:) 304 integer,parameter :: ncid0 = 0, master = 0 305 integer :: units(2) 306 ! ********************************************************************* 307 308 DBG_ENTER("COLL") 309 310 my_rank = xmpi_comm_rank(comm) 311 units = [std_out, ab_out] 312 313 !write(std_out,*) 'pawuj 01' 314 !########################################################### 315 !### 01. Allocations 316 317 !Initializations 318 ndtuj=count(dtpawuj(:)%iuj/=-1)-1 ! number of datasets initialized by pawuj_red 319 ABI_MALLOC(jdtset_,(0:ndtuj)) 320 jdtset_(0:ndtuj)=pack(dtpawuj(:)%iuj,dtpawuj(:)%iuj/=-1) 321 jdtset=maxval(dtpawuj(:)%iuj) 322 323 !DEBUG 324 write(message,'(10(a,i3))')'pawuj_det jdtset ',jdtset,& 325 & ' ndtuj ', ndtuj,' ndtpawuj ',ndtpawuj 326 call wrtout(std_out,message,'COLL') 327 !call flush_unit(6) 328 !END DEBUG 329 330 nspden=dtpawuj(jdtset)%nspden 331 nat_org=dtpawuj(jdtset)%nat 332 natom=dtset%natom 333 write(message,'(a,i3)')'natom,nat_org',nat_org 334 call wrtout(std_out,message,'COLL') 335 ! ABI_CHECK(nat_org == 1, "MG: I'm not sure we access data with the right index if nat_org > 1") 336 337 macro_uj=dtpawuj(jdtset)%macro_uj 338 pawujat=dtpawuj(jdtset)%pawujat 339 pawprtvol=dtpawuj(jdtset)%pawprtvol 340 pawujga=dtpawuj(jdtset)%pawujga 341 pawujoption=dtpawuj(jdtset)%option 342 ph0phiint=dtpawuj(jdtset)%ph0phiint 343 dmatpuopt=dtpawuj(jdtset)%dmatpuopt 344 marr=maxval((/ 9, nspden*nat_org ,nat_org*3 /)) 345 eyp=2.5_dp ! for dmatpuopt==2 and 3 346 if (dmatpuopt==1) eyp=eyp+3.0_dp 347 if (dmatpuopt>=3) eyp=(eyp+3.0_dp-dmatpuopt) 348 349 !DEBUG 350 write(message,'(a,i3)')'pawujdet dmatpuopt',dmatpuopt 351 call wrtout(std_out,message,'COLL') 352 !END DEBUG 353 354 ABI_MALLOC(luocc,(ndtpawuj,nat_org)) 355 ABI_MALLOC(magv_org,(nat_org)) 356 ABI_MALLOC(intmagv_org,(nat_org)) 357 ABI_MALLOC(xred_org,(3,nat_org)) 358 ABI_MALLOC(chi0,(nat_org)) 359 ABI_MALLOC(chi,(nat_org)) 360 ABI_MALLOC(distv_org,(nat_org)) 361 !DEBUG 362 !write(message,fmt='((a,i3,a))')'pawuj_det init sg' 363 !call wrtout(std_out,message,'COLL') 364 !END DEBUG 365 luocc=zero 366 367 !########################################################### 368 !### 02. Testing consistency of parameters and outputting info 369 370 !LMac 371 ! if (ndtuj/=4) return 372 373 write(message, '(3a)' ) ch10,' ---------- calculate U, (J) start ---------- ',ch10 374 call wrtout(ab_out,message,'COLL') 375 376 !Tests if perturbed atom is consistent with ujdet atom 377 if (all(dtpawuj(1:ndtpawuj)%pawujat==pawujat)) then 378 write (message,fmt='(a,i3)') ' All pawujat ok and equal to ',pawujat 379 call wrtout(ab_out,message,'COLL') 380 else 381 write (message,fmt='(a,4i3,2a)') ' Differing values of pawujat were found: ',dtpawuj(1:ndtuj)%pawujat,ch10,& 382 & 'No determination of U.' 383 call wrtout(ab_out,message,'COLL') 384 return 385 end if 386 387 !Tests consistency of macro_uj, then writes message about macro_uj procedure selected. LMac 388 if (nspden==1) then 389 write(message,fmt='(2a)') ' nspden==1, determination',& 390 & ' of U-parameter for unpolarized structure (non standard)' 391 else if (macro_uj==1.and.nspden==2) then 392 write(message,fmt='(2a)') ' macro_uj=1 and nspden=2:',& 393 & ' standard determination of Hubbard U-parameter' 394 else if (macro_uj==2.and.nspden==2) then 395 write(message,fmt='(2a)') ' macro_uj=2 and nspden=2:',& 396 & ' determination of parameter on single spin channel (experimental)' 397 else if (macro_uj==3.and.nspden==2) then 398 write(message,fmt='(2a)') ' macro_uj=3 and nspden=2,',& 399 & ' determination of (not Hunds) J-parameter on single spin channel (experimental)' 400 else if (macro_uj==4.and.nspden==2) then 401 write(message,fmt='(2a)') ' macro_uj=4 and nspden=2,',& 402 & ' Hunds J determination -- L. MacEnulty August 2021' 403 end if 404 call wrtout(ab_out,message,'COLL') 405 406 !Tests compatibility of nspden and macro_uj 407 if (macro_uj>1.and.nspden==1) then 408 write (message,'(4a,2a)') ' U on a single spin channel (or J) can only be determined for nspden=2 ,',ch10,& 409 & 'No determination of U.' 410 call wrtout(ab_out,message,'COLL') 411 return 412 end if 413 414 !Calculation of response matrix 415 diem=dtpawuj(1)%diemix !Unscreened response in Hubbard U impacted by diemix 416 pertname='alpha' !Hubbard U perturbation; applied equally to spins up and down 417 parname='U' 418 occmag=' Occupations' 419 signum=1.0d0 !Hubbard U is signum*(1/chi0-1/chi) 420 do jdtset=1,4 421 if (nspden==1) then 422 luocc(jdtset,1:nat_org)=dtpawuj(jdtset)%occ(1,:) 423 !Hubbard U: uses sum of spin up and spin down occupancies 424 else if (macro_uj==1.and.nspden==2) then 425 luocc(jdtset,1:nat_org)=dtpawuj(jdtset)%occ(1,:)+dtpawuj(jdtset)%occ(2,:) !Total occupation 426 else if (macro_uj==2.and.nspden==2) then 427 luocc(jdtset,1:nat_org)=dtpawuj(jdtset)%occ(1,:) 428 else if (macro_uj==3.and.nspden==2) then 429 luocc(jdtset,1:nat_org)=dtpawuj(jdtset)%occ(2,:) 430 parname='J' 431 !Hund's J: uses difference of spin up and spin down occupancies 432 else if (macro_uj==4.and.nspden==2) then 433 luocc(jdtset,1:nat_org)=dtpawuj(jdtset)%occ(1,:)-dtpawuj(jdtset)%occ(2,:) !Magnetization 434 diem=dtpawuj(1)%diemixmag !Unscreened response in Hund's J impacted by diemixmag 435 pertname='beta ' !Hund's J perturbation: +beta to spin up, -beta to down 436 parname='J' 437 occmag='Magnetizations' 438 signum=-1.0d0 !Hund's J is -1*(1/chi0-1/chi) 439 end if 440 ! MG: Cannot use pawujat index to extract vsh as we have used: 441 ! 442 ! dtpawuj(icyc)%vsh=reshape(pack(atvshift,atvshmusk),(/ nspden,nnat /)) 443 444 ! thus the pawujat atom in the atvshift array becomes the first one in %vsh 445 ! 446 !vsh(jdtset)=dtpawuj(jdtset)%vsh(1,pawujat) 447 vsh(jdtset)=dtpawuj(jdtset)%vsh(1,1) 448 if (pawprtvol==-3) then 449 write(message,fmt='(2a,i3,a,f15.12)') ch10,' Potential shift vsh(',jdtset,') =',vsh(jdtset) 450 call wrtout(std_out,message,'COLL') 451 write(message,fmt='( a,i3,a,120f15.9)') ' Occupations occ(',jdtset,') ',( luocc(jdtset,ii),ii=1,nat_org ) 452 call wrtout(std_out,message,'COLL') 453 end if 454 end do 455 456 write(message,fmt='( a)') 'Occupations assigned.' 457 call wrtout(std_out,message,'COLL') 458 459 !Two-point linear regression of response matrices. 460 chi0=(luocc(1,1:nat_org)-luocc(3,1:nat_org))/(vsh(1)-vsh(3))/diem 461 chi=(luocc(2,1:nat_org)-luocc(4,1:nat_org))/(vsh(2)-vsh(4)) 462 463 ! MG: pawujat replaced by 1 because arrays are dimensioned with nat_org (usually 1) and not natom! 464 if ((abs(chi0(1))<0.0000001).or.(abs(chi(1))<0.0000001)) then 465 write(message, '(2a,2f12.5,a)' ) ch10,'Chi0 or Chi is too small for inversion.',& 466 &chi0(1),chi(1),ch10 467 call wrtout(ab_out,message,'COLL') 468 return 469 end if 470 471 !LMac: Scalar Hubbard Parameter 472 scalarHP=signum*(1.0d0/chi0(1)-1.0d0/chi(1))*Ha_eV 473 474 write(message,fmt='(a)')': ' 475 if (nspden==2) then 476 magv_org=dtpawuj(1)%occ(1,:)-dtpawuj(1)%occ(2,:) 477 if (all(abs(magv_org)<0.001)) then 478 magv_org=(/(1,im1=1,nat_org)/) 479 else 480 magv_org=abs(magv_org)/magv_org 481 if (magv_org(1).lt.0) magv_org=magv_org*(-1_dp) 482 if (all(magv_org(2:nat_org).lt.0)) then 483 magv_org=abs(magv_org) 484 write(message,'(a)')', (reset to fm): ' 485 end if 486 end if 487 else 488 magv_org=(/(1,im1=1,nat_org)/) 489 end if 490 intmagv_org=(/(int(magv_org(im1)),im1=1,nat_org)/) 491 492 ! if (pawprtvol==-3) then 493 write(message,fmt='(3a, 150f4.1)') ch10,' Magnetization',trim(message),magv_org 494 call wrtout(std_out,message,'COLL') 495 ! end if 496 497 !Case of extrapolation to larger r_paw: calculate intg 498 if (all(dtpawuj(1)%wfchr(:)/=0).and.ph0phiint/=1) then 499 if (dtpawuj(1)%pawujrad<20.and.dtpawuj(1)%pawujrad>dtpawuj(1)%pawrad) then 500 fcorr=(1-ph0phiint)/(IRadFnH(dtpawuj(1)%pawrad,20.0_dp,nint(dtpawuj(1)%wfchr(2)),& 501 & nint(dtpawuj(1)%wfchr(3)),dtpawuj(1)%wfchr(1))) 502 intg=ph0phiint/(1-fcorr*IRadFnH(dtpawuj(1)%pawujrad,20.0_dp,nint(dtpawuj(1)%wfchr(2)),& 503 & nint(dtpawuj(1)%wfchr(3)),dtpawuj(1)%wfchr(1))) 504 write(message, fmt='(a,f12.5,a,f12.5)') ' pawuj_det: met2 extrapolation to ', dtpawuj(1)%pawujrad,' using factor ',intg 505 call wrtout(std_out,message,'COLL') 506 else if (dtpawuj(1)%pawujrad<dtpawuj(1)%pawrad) then 507 a=0 ; a(1:3)=dtpawuj(1)%wfchr(4:6) 508 a=(/a(1)**2,a(1)*a(2),a(2)**2/3.0_dp+2.0_dp/3.0_dp*a(1)*a(3),a(2)*a(3)/2.0_dp,a(3)**2/5.0_dp/) 509 b=(/(dtpawuj(1)%pawujrad**(im1)-dtpawuj(1)%pawrad**(im1),im1=1,5)/) 510 intg=dot_product(a,b) 511 intg=ph0phiint/(ph0phiint+intg) 512 write(message, fmt='(a,f12.5,a,f12.5)') ' pawuj_det: met1 extrapolation to ', dtpawuj(1)%pawujrad,' using factor ',intg 513 call wrtout(std_out,message,'COLL') 514 else if (dtpawuj(1)%pawujrad==dtpawuj(1)%pawrad) then 515 intg=one 516 write(message, fmt='(a,2i7,3f12.5)') ' pawuj_det: no extrapolation (pawujrad=pawrad)' 517 call wrtout(std_out,message,'COLL') 518 else 519 intg=ph0phiint 520 write(message, fmt='(a,3f12.5)') ' pawuj_det: extrapolation to r->\inf using factor ',intg 521 call wrtout(std_out,message,'COLL') 522 end if 523 else 524 write(message, fmt='(a,2i7,3f12.5)') ' pawuj_det: no extrapolation (ph0phiint=1 or wfchr=0)' 525 call wrtout(std_out,message,'COLL') 526 intg=one 527 end if 528 529 !Determine U in primitive cell 530 write(message,fmt='(a)')' pawuj_det: determine U in primitive cell' 531 call wrtout(std_out,message,'COLL') 532 533 call lcalcu(intmagv_org,nat_org,dtpawuj(1)%rprimd,dtpawuj(1)%xred,chi,chi0,pawujat,ures,pawprtvol,pawujga,pawujoption) 534 !Begin calculate U in supercell 535 536 !Analyze shell structure of primitive cell 537 !and atomic distances in primitive cell 538 ABI_MALLOC(smult_org,(nat_org)) 539 ABI_MALLOC(sdistv_org,(nat_org)) 540 call shellstruct(dtpawuj(1)%xred,dtpawuj(1)%rprimd,nat_org,& 541 & intmagv_org,distv_org,smult_org,sdistv_org,nsh_org,pawujat,pawprtvol) 542 543 544 !LMac: Printing relevant information about the Hubbard parameter just calculated. 545 write(message,'(3a)') ch10,ch10,'*********************************************************************' 546 call wrtout(units,message) 547 write(message,'(4a)') '************************ Linear Response ',parname,' ************************',ch10 548 call wrtout(units,message) 549 write(message, '(a,i4,a)' ) ' Info printed for perturbed atom: ',pawujat,ch10 550 call wrtout(units,message) 551 write(message, fmt='(10a)')' Perturbations ',occmag,ch10,& 552 ' --------------- -----------------------------',ch10,& 553 ' ',pertname,' [eV] Unscreened Screened',ch10,& 554 ' --------------- -----------------------------' 555 call wrtout(units,message) 556 ! MG: pawujat --> 1. 557 do ipert=1,2 558 write(message, fmt='(3f15.10)') vsh(ipert*2-1)*Ha_eV,luocc(ipert*2-1,1),luocc(ipert*2,1) 559 call wrtout(units,message) 560 end do 561 write(message,'(2a)') ch10,' Scalar response functions:' 562 call wrtout(units,message) 563 write(message,fmt='(a,f12.5)') ' Chi0 [eV^-1]: ',chi0(1)/Ha_eV 564 call wrtout(units,message) 565 write(message,fmt='(a,f12.5)') ' Chi [eV^-1]: ',chi(1)/Ha_eV 566 call wrtout(units,message) 567 write(message,'(4a,f9.5,a)') ch10,' The scalar ',parname,' from the two-point regression scheme is ',scalarHP,' eV.' 568 call wrtout(units,message) 569 write(message,'(3a)') '*********************************************************************',ch10,& 570 '*********************************************************************' 571 call wrtout(units,message) 572 write(message,'(7a)') 'Note: For more reliable linear regressions of the response',ch10,& 573 'matrices, it is advised that you have more than two points.',ch10,& 574 'See the LRUJ protocol for more information.',ch10,ch10 575 call wrtout(std_out,message,'COLL') 576 call wrtout(ab_out,message,'COLL') 577 578 !########################################################### 579 !### Create the file LRUJ.nc for the LRUJ ujdet utility 580 581 if (my_rank == master) then 582 ! First call: 583 ! - Create NC file, define dimensions, scalars and arrays. 584 ! - Add crystal structure and metadata required to post-process the data. 585 NCF_CHECK(nctk_open_create(ncid, strcat(dtfil%filnam_ds(4), "_LRUJ.nc"), xmpi_comm_self)) 586 587 cryst = dtset%get_crystal(1) 588 NCF_CHECK(cryst%ncwrite(ncid)) 589 call cryst%free() 590 591 !Define dimensions. 592 ncerr = nctk_def_dims(ncid, [ & 593 nctkdim_t("natom", natom), & 594 nctkdim_t("nnat", nat_org), & 595 nctkdim_t("ndtpawuj", ndtpawuj), & 596 nctkdim_t("nspden", dtset%nspden), & 597 nctkdim_t("nsppol", dtset%nsppol) ], & 598 defmode=.True.) 599 NCF_CHECK(ncerr) 600 601 ! Define integer scalars 602 ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: & 603 "usepaw", "macro_uj", "pawujat", "nspden", "dmatpuopt" & 604 ]) 605 606 ! Define double precision scalars 607 ! @lmacenul: Here I write pawujv so that one can order the results by alpha in the post-processing tool. 608 ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: & 609 "diemix", "diemixmag", "ph0phiint", "uj_pert" & 610 ]) 611 NCF_CHECK(ncerr) 612 613 ! Define arrays with results. 614 ! TODO: Here I need an extra dimension to store results with different iuj and/or different names. 615 ncerr = nctk_def_arrays(ncid, [ & 616 nctkarr_t("luocc", "dp", "ndtpawuj, nnat") & 617 !nctkarr_t("vsh", "dp", "nspden, nnat") & 618 ]) 619 NCF_CHECK(ncerr) 620 621 ! =========================================== 622 ! Write metadata that does not depend on icyc 623 ! =========================================== 624 NCF_CHECK(nctk_set_datamode(ncid)) 625 626 ncerr = nctk_write_iscalars(ncid, [character(len=nctk_slen) :: & 627 "usepaw", "macro_uj", "pawujat", "nspden", "dmatpuopt"], & 628 [dtset%usepaw, macro_uj, pawujat, nspden, dmatpuopt]) 629 NCF_CHECK(ncerr) 630 631 associate (dt => dtpawuj(1)) 632 ncerr = nctk_write_dpscalars(ncid, [character(len=nctk_slen) :: & 633 "diemix", "diemixmag", "ph0phiint", "uj_pert" ], & 634 [dt%diemix, dt%diemixmag, ph0phiint, vsh(3)]) 635 NCF_CHECK(ncerr) 636 end associate 637 638 ! Write arrays 639 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "luocc"), luocc)) 640 641 NCF_CHECK(nf90_close(ncid)) 642 end if 643 644 ii=1 645 write(message, fmt='(8a)') ' URES ',' ii',' nat',' r_max',' U(J)[eV]',' U_ASA[eV]',' U_inf[eV]',ch10 646 write(message, fmt='(a,2i7,4f12.5)') & 647 trim(message)//' URES ',ii,nat_org,maxval(abs(distv_org)),signum*ures,signum*ures*exp(log(intg)*eyp),& 648 signum*ures*exp(log(ph0phiint)*eyp) 649 call wrtout(units,message) 650 651 if (pawprtvol>1) then 652 write(message,fmt='(a, 150f10.5)')' pawuj_det: ionic distances in original cell (distv_org) ', distv_org 653 call wrtout(std_out,message,'COLL') 654 end if 655 656 !Construct supercell, calculate limit dimensions of supercell 657 ii=1 658 maxnat=product(dtpawuj(1)%scdim(:))*nat_org 659 if (maxnat==0) then 660 maxnat=dtpawuj(1)%scdim(1) 661 ext=(/ii, ii, ii/) 662 else 663 jj=1 664 do while (jj<=3) 665 ext(jj)=minval( (/ii, dtpawuj(1)%scdim(jj) /) ) 666 jj=jj+1 667 end do 668 end if 669 ext=ext+(/ 1,1,1 /) 670 ii=ii+1 671 672 673 nat_sc=product(ext)*nat_org 674 675 !DEBUG 676 !write(message,fmt='(a,3i4)')'pawuj_det: ext ',ext 677 !call wrtout(std_out,message,'COLL') 678 !END DEBUG 679 680 do while (nat_sc<=maxnat) 681 ABI_MALLOC(chi0_sc,(nat_sc)) 682 ABI_MALLOC(chi_sc,(nat_sc)) 683 ABI_MALLOC(distv_sc,(nat_sc)) 684 ABI_MALLOC(magv_sc,(nat_sc)) 685 ABI_MALLOC(intmagv_sc,(nat_sc)) 686 ABI_MALLOC(sdistv_sc,(nat_sc)) 687 ABI_MALLOC(smult_sc,(nat_sc)) 688 ABI_MALLOC(xred_sc,(3,nat_sc)) 689 690 ! Determine positions=xred_sc and supercell dimensions=rpimd_sc 691 call mksupercell(dtpawuj(1)%xred,intmagv_org,dtpawuj(1)%rprimd,nat_org,& 692 & nat_sc,xred_sc,magv_sc,rprimd_sc,ext,pawprtvol) 693 694 ! Determine shell structure of supercell: multiplicities (smult_sc), radii of shells (sdistv_sc) 695 ! number of shells (nsh_sc) and atomic distances in supercell (distv_sc) 696 697 write(message,fmt='(a,3i3,a)')' pawuj_det: determine shell structure of ',ext(1:3),' supercell' 698 call wrtout(std_out,message,'COLL') 699 700 intmagv_sc=(/(int(magv_sc(im1)),im1=1,nat_sc)/) 701 call shellstruct(xred_sc,rprimd_sc,nat_sc,intmagv_sc,distv_sc,smult_sc,sdistv_sc,nsh_sc,pawujat,pawprtvol) 702 703 if (pawprtvol>=2) then 704 write(message,fmt='(a)')' pawuj_det: ionic distances in supercell (distv_sc) ' 705 call wrtout(std_out,message,'COLL') 706 call prmat(distv_sc,1,nat_sc,1,std_out) 707 end if 708 709 ! Determine chi and chi0 in supercell (chi0_sc, chi_sc) 710 ! DEBUG 711 ! write(message,fmt='(a)')'pawuj_det: chi and chi0 in supercell' 712 ! call wrtout(std_out,message,'COLL') 713 ! END DEBUG 714 715 if (pawujoption>2) then 716 invopt=2 717 else 718 invopt=1 719 end if 720 721 call chiscwrt(chi0,distv_org,nat_org,sdistv_org,smult_org,nsh_org,& 722 & chi0_sc,distv_sc,nat_sc,smult_sc,nsh_sc,invopt,pawprtvol) 723 call chiscwrt(chi,distv_org,nat_org,sdistv_org,smult_org,nsh_org,& 724 & chi_sc,distv_sc,nat_sc,smult_sc,nsh_sc,invopt,pawprtvol) 725 726 ! Calculate U in supercell 727 ! DEBUG 728 ! write(message,fmt='(a)')'pawuj_det: U in supercell' 729 ! call wrtout(std_out,message,'COLL') 730 ! END DEBUG 731 call lcalcu(intmagv_sc,nat_sc,rprimd_sc,xred_sc,chi_sc,chi0_sc,pawujat,ures,pawprtvol,pawujga,pawujoption) 732 write(message, fmt='(a,2i7,4f12.5)') ' URES ',ii,nat_sc,maxval(abs(distv_sc)),signum*ures,signum*ures*exp(log(intg)*eyp),& 733 & signum*ures*exp(log(ph0phiint)*eyp) 734 call wrtout(units,message) 735 736 ABI_FREE(chi0_sc) 737 ABI_FREE(chi_sc) 738 ABI_FREE(distv_sc) 739 ABI_FREE(magv_sc) 740 ABI_FREE(intmagv_sc) 741 ABI_FREE(sdistv_sc) 742 ABI_FREE(smult_sc) 743 ABI_FREE(xred_sc) 744 745 if (product(dtpawuj(1)%scdim(:))==0) then 746 ext=(/ii, ii, ii/) 747 else 748 jj=1 749 do while (jj<=3) 750 ext(jj)=minval( (/ii, dtpawuj(1)%scdim(jj) /) ) 751 jj=jj+1 752 end do 753 end if 754 ext=ext+(/ 1,1,1 /) 755 ii=ii+1 756 757 nat_sc=product(ext)*nat_org 758 759 end do 760 761 write(message, '(3a)' )ch10,' ---------- calculate U, (J) end -------------- ',ch10 762 call wrtout(ab_out,message,'COLL') 763 764 ABI_FREE(jdtset_) 765 ABI_FREE(chi) 766 ABI_FREE(chi0) 767 ABI_FREE(smult_org) 768 ABI_FREE(sdistv_org) 769 ABI_FREE(luocc) 770 ABI_FREE(magv_org) 771 ABI_FREE(intmagv_org) 772 ABI_FREE(xred_org) 773 ABI_FREE(distv_org) 774 775 DBG_EXIT("COLL") 776 777 end subroutine pawuj_det
m_paw_uj/pawuj_free [ Functions ]
[ Top ] [ m_paw_uj ] [ Functions ]
NAME
pawuj_free
FUNCTION
deallocate pawuj stuff
INPUTS
OUTPUT
SOURCE
215 subroutine pawuj_free(dtpawuj) 216 217 !Arguments ------------------------------- 218 type(macro_uj_type),intent(inout) :: dtpawuj 219 220 ! ********************************************************************* 221 222 if (allocated(dtpawuj%scdim)) then 223 ABI_FREE(dtpawuj%scdim) 224 end if 225 if (allocated(dtpawuj%occ)) then 226 ABI_FREE(dtpawuj%occ) 227 end if 228 if (allocated(dtpawuj%rprimd)) then 229 ABI_FREE(dtpawuj%rprimd) 230 end if 231 if (allocated(dtpawuj%vsh)) then 232 ABI_FREE(dtpawuj%vsh) 233 end if 234 if (allocated(dtpawuj%xred)) then 235 ABI_FREE(dtpawuj%xred) 236 end if 237 if (allocated(dtpawuj%wfchr)) then 238 ABI_FREE(dtpawuj%wfchr) 239 end if 240 if (allocated(dtpawuj%zioneff)) then 241 ABI_FREE(dtpawuj%zioneff) 242 end if 243 244 end subroutine pawuj_free
m_paw_uj/pawuj_ini [ Functions ]
[ Top ] [ m_paw_uj ] [ Functions ]
NAME
pawuj_ini
FUNCTION
Initialize dtpawuj datastructure for one SCF run Relevant only for automatic determination of U in PAW+U context
INPUTS
OUTPUT
SIDE EFFECTS
dtpawuj(0:ndtpawuj) (initialization of fields vsh, occ, iuj,nnat)
SOURCE
150 subroutine pawuj_ini(dtpawuj,ndtset) 151 152 !Arguments ------------------------------------ 153 !scalars 154 integer :: ndtset 155 type(macro_uj_type),intent(inout) :: dtpawuj(0:ndtset) 156 157 !Local variables ------------------------- 158 !Variables for partial dos calculation 159 !scalars 160 integer, parameter :: nwfchr=6 161 integer :: iuj,im1 162 ! ********************************************************************* 163 164 DBG_ENTER("COLL") 165 do iuj=0,ndtset 166 !write(std_out,*)'pawuj_ini iuj ',iuj 167 dtpawuj(iuj)%diemix=0.45_dp 168 dtpawuj(iuj)%diemixmag=0.45_dp 169 dtpawuj(iuj)%iuj=0 170 dtpawuj(iuj)%nat=0 171 dtpawuj(iuj)%ndtset=1 172 dtpawuj(iuj)%nspden=1 173 dtpawuj(iuj)%macro_uj=0 174 dtpawuj(iuj)%option=1 175 dtpawuj(iuj)%pawujat=1 176 dtpawuj(iuj)%pawujga=one 177 dtpawuj(iuj)%pawprtvol=1 178 dtpawuj(iuj)%ph0phiint=one 179 dtpawuj(iuj)%dmatpuopt=2 180 dtpawuj(iuj)%pawujrad=3.0_dp 181 dtpawuj(iuj)%pawrad=20.0_dp 182 !Allocate arrays 183 !write(std_out,*)'pawuj_ini before arrays' 184 ABI_MALLOC(dtpawuj(iuj)%rprimd,(3,3)) 185 ABI_MALLOC(dtpawuj(iuj)%scdim,(3)) 186 ABI_MALLOC(dtpawuj(iuj)%wfchr,(nwfchr)) 187 dtpawuj(iuj)%rprimd=reshape((/ 1,0,0,0,1,0,0,0,1/),(/ 3,3 /)) 188 dtpawuj(iuj)%scdim=reshape((/ 250,0,0 /),(/3 /)) 189 dtpawuj(iuj)%wfchr=(/ (0,im1=1,nwfchr) /) 190 if (iuj>0) then 191 dtpawuj(iuj)%iuj=-1 192 end if 193 end do 194 195 DBG_EXIT("COLL") 196 197 end subroutine pawuj_ini
m_paw_uj/pawuj_red [ Functions ]
[ Top ] [ m_paw_uj ] [ Functions ]
NAME
pawuj_red
FUNCTION
Store atomic occupancies, potential shift, positions in dtpawuj datastructure.
INPUTS
istep: SCF iteration step pert_state: 0 if the routine is called with unperturbed occupancies 1 if the routine is called after having applied the perturbation. fatvshift=factor that multiplies atvshift 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 ntypat = number of atom types paw_ij(my_natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels pawprtvol= printing volume pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data comm=MPI commuicator
OUTPUT
dtpawuj(0:ndtpawuj) (initialization of fields vsh, occ, occ0, iuj,nnat)
SOURCE
808 subroutine pawuj_red(istep, pert_state, dtfil, & 809 dtset,dtpawuj,fatvshift,my_natom,natom,ntypat,paw_ij,pawrad,pawtab,ndtpawuj, comm, & 810 & mpi_atmtab,comm_atom) ! optional arguments (parallelism) 811 812 !Arguments ------------------------------------ 813 !scalars 814 integer,intent(in) :: istep, pert_state, my_natom,natom,ntypat,ndtpawuj, comm 815 integer,optional,intent(in) :: comm_atom 816 real(dp),intent(in) :: fatvshift 817 !arrays 818 integer,optional,target,intent(in) :: mpi_atmtab(:) 819 type(paw_ij_type),intent(in) :: paw_ij(my_natom) 820 type(pawtab_type),intent(in) :: pawtab(ntypat) 821 type(pawrad_type),intent(in) :: pawrad(ntypat) 822 type(dataset_type),intent(in) :: dtset 823 type(datafiles_type),intent(in) :: dtfil 824 type(macro_uj_type),intent(inout) :: dtpawuj(0:ndtpawuj) 825 826 !Local variables------------------------------- 827 !scalars 828 integer,parameter :: natmax=2,ncoeff=3, master = 0 829 integer :: iatom,iatom_tot,ierr,im1,im2,ispden,itypat,ll,nspden,nsppol,iuj,ncyc,icyc 830 integer :: my_comm_atom,nnat,natpawu,natvshift,pawujat,ndtset,typawujat 831 !integer :: my_rank, ncid, ncerr 832 logical :: usepawu !antiferro, 833 logical :: my_atmtab_allocated,paral_atom 834 character(len=1000) :: message,hstr 835 character(len=500) :: messg 836 !arrays 837 logical :: musk(3,natom) 838 integer,pointer :: my_atmtab(:) 839 real(dp) :: rrtab(ncoeff),wftab(ncoeff),a(ncoeff,ncoeff),b(ncoeff,ncoeff)! ,a(ncoeff,ncoeff) 840 real(dp),allocatable :: nnocctot(:,:) !,magv(:) 841 real(dp),allocatable :: atvshift(:,:,:) ! atvshift(natvshift,2,natom) 842 logical,allocatable :: dmusk(:,:),atvshmusk(:,:,:) !atvshmusk(natvshift,2,natom) 843 844 ! ********************************************************************* 845 846 ABI_UNUSED((/dtfil%ireadwf, istep, pert_state, comm/)) 847 848 !Initializations 849 nspden=1;nsppol=1 850 if (my_natom>0) then 851 nspden=paw_ij(1)%nspden ; nsppol=paw_ij(1)%nsppol 852 end if 853 854 !my_rank = xmpi_comm_rank(comm) 855 856 natvshift=dtset%natvshift 857 pawujat=dtset%pawujat 858 natpawu=dtset%natpawu ; ndtset=dtset%ndtset 859 ABI_MALLOC(atvshift,(natvshift,nspden,natom)) 860 ABI_MALLOC(atvshmusk,(natvshift,nspden,natom)) 861 ABI_MALLOC(dmusk,(nspden,natom)) 862 ABI_MALLOC(nnocctot,(nspden,natom)) 863 musk=.false.; dmusk=.false. 864 atvshift=fatvshift*dtset%atvshift 865 atvshmusk=.false. 866 nnocctot=zero 867 typawujat=dtset%typat(pawujat) 868 usepawu=(count(pawtab(:)%usepawu/=0)>0) 869 870 !Set up parallelism over atoms 871 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 872 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 873 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 874 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) 875 876 nnat=0 877 if (usepawu) then 878 write(message,'(3a)') ch10, '---------- pawuj_red ------ ',ch10 879 call wrtout(std_out, message,'COLL'); 880 do iatom=1,my_natom 881 iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom) 882 itypat=paw_ij(iatom)%itypat;ll=pawtab(itypat)%lpawu 883 if ((ll>=0).and.(pawtab(itypat)%usepawu/=0).and.itypat==typawujat) then 884 musk(:,iatom_tot)=(/.true., .true., .true. /) 885 atvshmusk(:,:,iatom_tot)=reshape((/ (( (im1==1), im1=1,natvshift) ,im2=1,nspden ) /),(/natvshift,nspden/)) 886 do ispden=1,nspden 887 nnocctot(ispden,iatom_tot)=paw_ij(iatom)%nocctot(ispden) 888 dmusk(ispden,iatom_tot)=.true. 889 end do 890 nnat=nnat+1 891 end if 892 end do 893 894 ! Reduction in case of parallelism 895 if (paral_atom) then 896 call xmpi_sum(nnocctot ,my_comm_atom,ierr) 897 call xmpi_lor(atvshmusk,my_comm_atom) ! dim=natpawu ??? 898 call xmpi_lor(dmusk ,my_comm_atom) 899 call xmpi_lor(musk ,my_comm_atom) 900 call xmpi_sum(nnat ,my_comm_atom,ierr) 901 end if 902 903 iuj=maxval(dtpawuj(:)%iuj) 904 write(std_out,*)' pawuj_red: iuj',iuj 905 ! write(std_out,*)' pawuj_red: dtpawuj(:)%iuj ',(dtpawuj(ii)%iuj,ii=1,ndtpawuj) 906 907 !If this is the unperturbed state, then unscreened and screened occupancies 908 !are the same. Also set perturbation vsh to zero. 909 !If this is the unperturbed case, then do everything twice; once for iuj=1 910 !and the same for iuj=2. If this is the perturbed case, do everything only 911 !once. LMac 912 if (iuj==1) then 913 ABI_MALLOC(dtpawuj(0)%vsh,(nspden,nnat)) 914 ABI_MALLOC(dtpawuj(0)%occ,(nspden,nnat)) 915 ABI_MALLOC(dtpawuj(0)%xred,(3,nnat)) 916 dtpawuj(0)%vsh=0 917 dtpawuj(0)%occ=0 918 dtpawuj(0)%xred=0 919 ncyc=2 920 else 921 ncyc=iuj 922 end if 923 924 do icyc=iuj,ncyc 925 926 if (icyc==1) then ! 1 and 3: non-scf steps 927 dtpawuj(icyc+1)%iuj=icyc+1 928 dtpawuj(icyc+2)%iuj=icyc+2 929 else if (icyc==3) then 930 dtpawuj(icyc+1)%iuj=icyc+1 931 end if 932 933 ! TODO: check that this is correct: this point is passed several times 934 ! for a given value of iuj - should the stuff be accumulated instead of replaced? 935 if(allocated(dtpawuj(icyc)%vsh)) then 936 ABI_FREE(dtpawuj(icyc)%vsh) 937 end if 938 ABI_MALLOC(dtpawuj(icyc)%vsh,(nspden,nnat)) 939 !Allocate screened occupancy array: LMac 940 if(allocated(dtpawuj(icyc)%occ)) then 941 ABI_FREE(dtpawuj(icyc)%occ) 942 end if 943 ABI_MALLOC(dtpawuj(icyc)%occ,(nspden,nnat)) 944 if(allocated(dtpawuj(icyc)%xred)) then 945 ABI_FREE(dtpawuj(icyc)%xred) 946 end if 947 ABI_MALLOC(dtpawuj(icyc)%xred,(3,nnat)) 948 949 950 rrtab=(/0.75_dp,0.815_dp,1.0_dp/)*pawtab(typawujat)%rpaw 951 wftab=pawtab(typawujat)%phi(pawtab(typawujat)%mesh_size,pawtab(typawujat)%lnproju(1)) 952 953 do im1=1,ncoeff 954 if (pawrad(typawujat)%mesh_type==1) then 955 im2=nint(rrtab(im1)/pawrad(typawujat)%rstep+1) 956 else if (pawrad(typawujat)%mesh_type==2) then 957 im2=nint(log(rrtab(im1)/pawrad(typawujat)%rstep+1)/pawrad(typawujat)%lstep+1) 958 else if (pawrad(typawujat)%mesh_type==3) then 959 im2=nint(log(rrtab(im1)/pawrad(typawujat)%rstep)/pawrad(typawujat)%lstep+1) 960 else if (pawrad(typawujat)%mesh_type==4) then 961 im2=nint(pawtab(typawujat)%mesh_size*(1-exp((-one)*rrtab(im1)/pawrad(typawujat)%rstep))+1) 962 end if 963 964 rrtab(im1)=pawrad(typawujat)%rad(im2) 965 wftab(im1)=pawtab(typawujat)%phi(im2,pawtab(typawujat)%lnproju(1)) 966 end do 967 write(message,fmt='(a,i3,a,10f10.5)')' pawuj_red: mesh_type',pawrad(typawujat)%mesh_type,' rrtab:', rrtab 968 call wrtout(std_out,message,'COLL') 969 write(message,fmt='(a,10f10.5)')' pawuj_red: wftab', wftab 970 call wrtout(std_out,message,'COLL') 971 a=reshape((/ (( rrtab(im2)**(im1-1), im1=1,3) ,im2=1,3 )/),(/ncoeff,ncoeff/)) 972 write(messg,fmt='(a)')'A' 973 call linvmat(a,b,ncoeff,messg,2,0.0_dp,3) ! linvmat(inmat,oumat,nat,nam,option,gam,prtvol) 974 write(std_out,*) 'pawuj_red: a,b ', a,b 975 wftab=matmul(wftab,b) 976 write(std_out,*) 'pawuj_red: wftab ', wftab 977 dtpawuj(icyc)%wfchr(4:6)=wftab 978 979 dtpawuj(icyc)%nat=nnat 980 write(std_out,*) 'pawuj_red: m1' 981 dtpawuj(icyc)%vsh=reshape(pack(atvshift,atvshmusk),(/ nspden,nnat /)) 982 !factor in next line to compensate nocctot contains just occ of 1 spin channel for nspden=1 983 write(std_out,*) 'pawuj_red: m2' 984 dtpawuj(icyc)%occ=reshape(pack(nnocctot,dmusk),(/nspden,nnat/))*(3-nspden) 985 write(std_out,*) 'pawuj_red: m3' 986 987 write(std_out,*) 'pawuj_red: occ ', dtpawuj(icyc)%occ 988 989 dtpawuj(icyc)%xred=reshape(pack(dtset%xred_orig(:,:,1),musk),(/3,nnat/)) 990 dtpawuj(icyc)%ph0phiint=pawtab(typawujat)%ph0phiint(1) 991 dtpawuj(icyc)%wfchr(1:3)=(/ pawtab(typawujat)%zioneff(1)*(dtset%lpawu(typawujat)+2),& 992 & one*(dtset%lpawu(typawujat)+1),one*(dtset%lpawu(typawujat))/) 993 dtpawuj(icyc)%pawrad=pawtab(typawujat)%rpaw 994 995 write(std_out,*) 'pawuj_red: wfchr ',dtpawuj(icyc)%wfchr 996 997 if (icyc.le.2) dtpawuj(icyc)%vsh=0.0d0 998 999 write (hstr,'(I0)') icyc 1000 write(message,'(a,a,I3,a)') ch10, '---------- MARK ------ ',icyc,ch10 1001 call wrtout(std_out,message,'COLL') 1002 write(message,fmt='(a)') 'vsh'//trim(hstr) 1003 call wrtout(std_out,message,'COLL') 1004 call prmat(dtpawuj(icyc)%vsh(:,:),1,nnat*nspden,1) 1005 write(message,fmt='(a)') 'occ'//trim(hstr) 1006 call wrtout(std_out,message,'COLL') 1007 call prmat(dtpawuj(icyc)%occ(:,:),1,nnat*nspden,1) 1008 write(message, '(3a)' )'---------- MARK ---------- ',ch10 1009 call wrtout(std_out,message,'COLL') 1010 end do 1011 end if !usepawu 1012 1013 ABI_FREE(nnocctot) 1014 ABI_FREE(dmusk) 1015 ABI_FREE(atvshift) 1016 ABI_FREE(atvshmusk) 1017 1018 !Destroy atom table used for parallelism 1019 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 1020 1021 end subroutine pawuj_red