TABLE OF CONTENTS


ABINIT/m_paw_uj [ Modules ]

[ Top ] [ 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