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-2018 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

18 #if defined HAVE_CONFIG_H
19 #include "config.h"
20 #endif
21 
22 #include "abi_common.h"
23 
24 MODULE m_paw_uj
25 
26  use defs_basis
27  use defs_abitypes
28  use m_abicore
29  use m_errors
30  use m_linalg_interfaces
31  use m_xmpi
32 
33  use m_pptools,       only : prmat
34  use m_special_funcs, only : iradfnh
35  use m_geometry,      only : shellstruct,ioniondist
36  use m_parser,        only : prttagm
37  use m_supercell,     only : mksupercell
38  use m_pawrad,        only : pawrad_type
39  use m_pawtab,        only : pawtab_type
40  use m_paw_ij,        only : paw_ij_type
41  use m_paral_atom,    only : get_my_atmtab, free_my_atmtab
42 
43  implicit none
44 
45  private
46 
47 !public procedures.
48  public :: pawuj_ini  ! Initialize dtpawuj datastructure
49  public :: pawuj_free ! Deallocate dtpawuj datastructure
50  public :: pawuj_det  ! Determine U (or J) parameter
51  public :: pawuj_red  ! Store atomic occupancies, potential shift, positions
52 
53 !private procedures.
54 !  chiscwrt:   distributes values of chi_org on chi_sc according to ion-ion distances
55 !  linvmat:    inverts real matrix inmat
56 !  lprtmat:    prints out the real matrix mmat
57 !  lcalcu:     prints out real the real matrice mmat
58 !  blow_pawuj: reads a real nxn matrice and appends lines n+1 and clumn n+1
59 
60 CONTAINS  !========================================================================================

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

PARENTS

      pawuj_utils

CHILDREN

SOURCE

1441 subroutine blow_pawuj(mat,nj,matt)
1442 
1443 
1444 !This section has been created automatically by the script Abilint (TD).
1445 !Do not modify the following lines by hand.
1446 #undef ABI_FUNC
1447 #define ABI_FUNC 'blow_pawuj'
1448 !End of the abilint section
1449 
1450  implicit none
1451 
1452 !Arguments ------------------------------------
1453 !scalars
1454  integer,intent(in)        :: nj
1455 !arrays
1456  real(dp),intent(in)       :: mat(nj,nj)
1457  real(dp),intent(out)      :: matt(nj+1,nj+1)
1458 
1459 !Local variables-------------------------------
1460 !scalars
1461  integer                   :: ii
1462 !arrays
1463 
1464 ! *************************************************************************
1465 
1466  matt(1:nj,1:nj)=mat
1467  do  ii = 1,nj
1468    matt(ii,nj+1)=-sum(matt(ii,1:nj))
1469  end do
1470 
1471  do  ii = 1,nj+1
1472    matt(nj+1,ii)=-sum(matt(1:nj,ii))
1473  end do
1474 
1475 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

PARENTS

      pawuj_det

CHILDREN

      prmat,wrtout

SOURCE

1026 subroutine chiscwrt(chi_org,disv_org,nat_org,sdisv_org,smult_org,nsh_org,chi_sc,&
1027 & disv_sc,nat_sc,smult_sc,nsh_sc,opt,prtvol)
1028 
1029 
1030 !This section has been created automatically by the script Abilint (TD).
1031 !Do not modify the following lines by hand.
1032 #undef ABI_FUNC
1033 #define ABI_FUNC 'chiscwrt'
1034 !End of the abilint section
1035 
1036  implicit none
1037 
1038 !Arguments ------------------------------------
1039 !scalars
1040  integer,intent(in)              :: nat_org,nsh_org,nsh_sc
1041  integer,intent(in)              :: nat_sc
1042  integer,intent(in),optional     :: prtvol
1043  integer,intent(in),optional     :: opt
1044 !arrays
1045  real(dp),intent(in)             :: chi_org(nat_org),disv_org(nat_org),sdisv_org(nsh_org)
1046  integer,intent(in)              :: smult_org(nsh_org), smult_sc(nsh_sc)
1047  real(dp),intent(in)             :: disv_sc(nat_sc)
1048  real(dp),intent(out)            :: chi_sc(nat_sc)
1049 
1050 !Local variables-------------------------------
1051 !scalars
1052  integer                      :: iatom,jatom,jsh,optt
1053  character(len=500)           :: message
1054 !arrays
1055  real(dp)                     :: chi_orgl(nat_org)
1056 
1057 ! *************************************************************************
1058 
1059  if (present(opt)) then
1060    optt=opt
1061  else
1062    optt=1
1063  end if
1064 
1065 
1066  do iatom=1,nat_org
1067    do jsh=1,nsh_org
1068      if (disv_org(iatom)==sdisv_org(jsh)) then
1069        if (opt==2) then
1070          chi_orgl(iatom)=chi_org(iatom)
1071        else if (opt==1) then
1072          chi_orgl(iatom)=chi_org(iatom)*smult_org(jsh)/smult_sc(jsh)
1073        end if
1074        exit
1075      end if
1076    end do !iatom
1077  end do  !jsh
1078 
1079  if (prtvol>1) then
1080    write(message,fmt='(a,150f10.5)')' chiscwrt: chi at input ',chi_org
1081    call wrtout(std_out,message,'COLL')
1082    write(message,fmt='(a,150f10.5)')' chiscwrt: chi after division ',chi_orgl
1083    call wrtout(std_out,message,'COLL')
1084  end if
1085 
1086  do iatom=1,nat_sc
1087    do jatom=1,nat_org
1088      if (disv_org(jatom)==disv_sc(iatom)) then
1089        chi_sc(iatom)=chi_orgl(jatom)
1090        exit
1091      else if (jatom==nat_org) then
1092        chi_sc(iatom)=0_dp
1093      end if
1094    end do
1095  end do
1096 
1097  if (prtvol>1) then
1098    write(message,'(a)')' chiscwrt, chi_sc '
1099    call wrtout(std_out,message,'COLL')
1100    call prmat(chi_sc,1,nat_sc,1,std_out)
1101  end if
1102 
1103 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

PARENTS

      pawuj_det

CHILDREN

SOURCE

1338 subroutine lcalcu(magv,natom,rprimd,xred,chi,chi0,pawujat,ures,prtvol,gam,opt)
1339 
1340 
1341 !This section has been created automatically by the script Abilint (TD).
1342 !Do not modify the following lines by hand.
1343 #undef ABI_FUNC
1344 #define ABI_FUNC 'lcalcu'
1345 !End of the abilint section
1346 
1347  implicit none
1348 
1349 !Arguments ------------------------------------
1350 !scalars
1351  integer,intent(in)          :: natom
1352  integer,intent(in),optional :: opt,pawujat,prtvol
1353  real(dp),intent(in),optional:: gam
1354  real(dp),intent(out)        :: ures
1355 !arrays
1356  integer,intent(in)          :: magv(natom)
1357  real(dp),intent(in)         :: rprimd(3,3),xred(3,natom),chi(natom),chi0(natom)
1358 
1359 !Local variables-------------------------------
1360 !scalars
1361  character(len=500)          :: message
1362  integer                     :: optt,prtvoll,nnatom
1363  real(dp)                    :: gamm
1364 !arrays
1365  real(dp),allocatable        :: tab(:,:,:)
1366 
1367 ! *********************************************************************
1368 
1369  if (present(opt)) then
1370    optt=opt
1371  else
1372    optt=1
1373  end if
1374 
1375  if (present(prtvol)) then
1376    prtvoll=prtvol
1377  else
1378    prtvoll=1
1379  end if
1380 
1381  if (present(gam)) then
1382    gamm=gam
1383  else
1384    gamm=1_dp
1385  end if
1386 
1387  if (optt==1.or.optt==3) then
1388    nnatom=natom+1
1389  else
1390    nnatom=natom
1391  end if
1392 
1393  ABI_ALLOCATE(tab,(4,nnatom,nnatom))
1394 
1395  call ioniondist(natom,rprimd,xred,tab(1,1:natom,1:natom),3,chi0,magv,pawujat,prtvoll)
1396  call ioniondist(natom,rprimd,xred,tab(2,1:natom,1:natom),3,chi,magv,pawujat,prtvoll)
1397 
1398 
1399  write(message,fmt='(a)')'response chi_0'
1400  call linvmat(tab(1,1:natom,1:natom),tab(3,1:nnatom,1:nnatom),natom,message,optt,gamm,prtvoll)
1401 
1402  write(message,fmt='(a)')'response chi'
1403  call linvmat(tab(2,1:natom,1:natom),tab(4,1:nnatom,1:nnatom),natom,message,optt,gamm,prtvoll)
1404 
1405  tab(1,1:nnatom,1:nnatom)=(tab(3,1:nnatom,1:nnatom)-tab(4,1:nnatom,1:nnatom))*Ha_eV
1406 
1407  write(message,fmt='(a,i3,a)')' (chi_0)^(-1)-(chi)^(-1) (eV)'
1408  call lprtmat(message,2,prtvoll,tab(1,1:nnatom,1:nnatom),nnatom)
1409 
1410  ures=tab(1,1,pawujat)
1411 
1412  ABI_DEALLOCATE(tab)
1413 
1414 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

PARENTS

      pawuj_red,pawuj_utils

CHILDREN

SOURCE

1136 subroutine linvmat(inmat,oumat,nat,nam,option,gam,prtvol)
1137 
1138 
1139 !This section has been created automatically by the script Abilint (TD).
1140 !Do not modify the following lines by hand.
1141 #undef ABI_FUNC
1142 #define ABI_FUNC 'linvmat'
1143 !End of the abilint section
1144 
1145  implicit none
1146 
1147 !Arguments -------------------------------
1148 
1149  integer,intent(in)              :: nat
1150  real(dp),intent(in)             :: gam,inmat(nat,nat)
1151  real(dp),intent(inout)          :: oumat(:,:)         ! nat+1,nat+1 for option=1 or 2
1152                                                        ! nat,nat for option=2
1153  character(len=500),intent(in)   :: nam
1154  integer,intent(in),optional     :: prtvol,option
1155 
1156 !Local variables -------------------------
1157  character(len=500)             :: message
1158  character(len=500)             :: bastrin,gastrin
1159  integer                        :: info,nnat,optionn
1160  integer,allocatable            :: ipvt(:)
1161  real(dp),allocatable           :: hma(:,:),work(:)
1162 
1163 ! *********************************************************************
1164 
1165  if (present(option)) then
1166    optionn=option
1167  else
1168    optionn=1
1169  end if
1170 
1171  if (option==1.or.option==3) then
1172    write(bastrin,'(a)')'+ charge bath '
1173    write(gastrin,'(a,d10.2,a)')'+ gamma  (=',gam,') '
1174  else
1175    write(bastrin,'(a)')''
1176    write(gastrin,'(a)')''
1177  end if
1178 
1179  write(message,fmt='(a)')' matrix '//trim(nam)
1180  call lprtmat(message,1,prtvol,inmat,nat)
1181 
1182  if (option==1.or.option==3) then
1183    call blow_pawuj(inmat,nat,oumat)
1184    write(message,fmt='(a,a)')' ',trim(nam)//trim(bastrin)
1185    call lprtmat(message,1,prtvol,oumat,nat+1)
1186    oumat=oumat+gam
1187    write(message,fmt='(a,a)')' ',trim(nam)//trim(bastrin) ! //trim(gastrin)
1188    call lprtmat(message,1,prtvol,oumat-gam,nat+1)
1189    nnat=nat+1
1190  else
1191    nnat=nat
1192    oumat=inmat
1193    oumat(1,1)=inmat(1,1)
1194  end if
1195 
1196 
1197  ABI_ALLOCATE(hma,(nnat,nnat))
1198  ABI_ALLOCATE(work,(nnat))
1199  ABI_ALLOCATE(ipvt,(nnat))
1200  work=0_dp
1201  hma(:,:)=oumat
1202 
1203  call dgetrf(nnat,nnat,hma,nnat,ipvt,info)
1204  if (.not.info==0) then
1205    write(message, '(3a)' ) 'Matrix '//trim(nam)//' is singular',ch10,'Probably too many symmetries kept'
1206    call wrtout(ab_out,message,'COLL')
1207    return
1208  end if
1209 
1210  call dgetri(nnat,hma,nnat,ipvt,work,nnat,info)
1211  oumat=hma(:,:)
1212 
1213  write(message,fmt='(2a,a)')' ('//trim(nam)//trim(bastrin)//trim(gastrin)//')^(-1)'
1214  call lprtmat(message,1,prtvol,oumat,nnat)
1215 
1216  ABI_DEALLOCATE(hma)
1217  ABI_DEALLOCATE(work)
1218  ABI_DEALLOCATE(ipvt)
1219 
1220 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

PARENTS

      pawuj_utils

CHILDREN

SOURCE

1254 subroutine lprtmat(commnt,chan,prtvol,mmat,nat)
1255 
1256 
1257 !This section has been created automatically by the script Abilint (TD).
1258 !Do not modify the following lines by hand.
1259 #undef ABI_FUNC
1260 #define ABI_FUNC 'lprtmat'
1261 !End of the abilint section
1262 
1263  implicit none
1264 
1265 !Arguments -------------------------------
1266  integer,intent(in)              :: nat,chan,prtvol
1267  real(dp),intent(in)             :: mmat(nat,nat)
1268  character(len=500),intent(in)  :: commnt
1269 
1270 !Local variables -------------------------
1271  character(len=500)             :: message
1272 ! *********************************************************************
1273 
1274  if (prtvol==3) then
1275    write(message,fmt='(a)') trim(commnt)
1276    call wrtout(std_out,message,'COLL')
1277    call prmat(mmat,nat,nat,nat,std_out)
1278    if (chan==2) then
1279      call wrtout(ab_out,message,'COLL')
1280      call prmat(mmat,nat,nat,nat,ab_out)
1281    end if
1282    write(message,*)ch10
1283    call wrtout(std_out,message,'COLL')
1284    if (chan==2) then
1285      call wrtout(ab_out,message,'COLL')
1286    end if
1287  end if
1288 
1289  if (prtvol==2) then
1290    write(message,fmt='(a)') trim(commnt)
1291    call wrtout(std_out,message,'COLL')
1292    call prmat(mmat,1,nat,nat,std_out)
1293    if (chan==2) then
1294      call wrtout(ab_out,message,'COLL')
1295      call prmat(mmat,1,nat,nat,ab_out)
1296    end if
1297    write(message,*)ch10
1298    call wrtout(std_out,message,'COLL')
1299    if (chan==2) then
1300      call wrtout(ab_out,message,'COLL')
1301    end if
1302  end if
1303 
1304 end subroutine lprtmat

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 U in PAW+U context

INPUTS

  dtpawuj=potential shifts (vsh) and atomic occupations (occ)
  ujdet_filename= Filename for write (Using NetCDF format)

OUTPUT

  only printing
  (among other things a section in the ab.out that can be used for input in ujdet)

PARENTS

      pawuj_drive,ujdet

CHILDREN

      chiscwrt,lcalcu,mksupercell,prmat,prttagm,shellstruct,wrtout

SOURCE

233 subroutine pawuj_det(dtpawuj,ndtpawuj,ujdet_filename,ures)
234 
235 
236 !This section has been created automatically by the script Abilint (TD).
237 !Do not modify the following lines by hand.
238 #undef ABI_FUNC
239 #define ABI_FUNC 'pawuj_det'
240 !End of the abilint section
241 
242  implicit none
243 
244 !Arguments ------------------------------------
245 !scalars
246 !arrays
247  integer                        :: ndtpawuj
248  type(macro_uj_type),intent(in) :: dtpawuj(0:ndtpawuj)
249  real(dp),intent(out)           :: ures
250  character(len=*),intent(in)    :: ujdet_filename
251 
252 !Local variables-------------------------------
253 !scalars
254  integer,parameter           :: natmax=2,nwfchr=6
255  integer                     :: ii,jj,nat_org,jdtset,nspden,macro_uj,kdtset,marr
256  integer                     :: im1,ndtuj,idtset, nsh_org, nsh_sc,nat_sc,maxnat
257  integer                     :: pawujat,pawprtvol,pawujoption
258  integer                     :: dmatpuopt,invopt
259  real(dp)                    :: pawujga,ph0phiint,intg,fcorr,eyp
260 
261  character(len=500)          :: message
262  character(len=2)            :: hstr
263 !arrays
264  integer                     :: ext(3)
265  real(dp)                    :: rprimd_sc(3,3),vsh(ndtpawuj),a(5),b(5)
266  integer,allocatable         :: narrm(:)
267  integer,allocatable         :: idum2(:,:),jdtset_(:),smult_org(:),smult_sc(:)
268  real(dp),allocatable        :: chih(:,:),dqarr(:,:),dqarrr(:,:),dparr(:,:),dparrr(:,:),xred_org(:,:),drarr(:,:)
269  real(dp),allocatable        :: magv_org(:),magv_sc(:),chi_org(:),chi0_org(:),chi0_sc(:), chi_sc(:), xred_sc(:,:)
270  real(dp),allocatable        :: sdistv_org(:),sdistv_sc(:),distv_org(:),distv_sc(:)
271  integer                     :: ncid=0
272 ! *********************************************************************
273 
274  DBG_ENTER("COLL")
275 
276 !write(std_out,*) 'pawuj 01'
277 !###########################################################
278 !### 01. Allocations
279 
280 !Initializations
281  ndtuj=count(dtpawuj(:)%iuj/=-1)-1 ! number of datasets initialized by pawuj_red
282  ABI_ALLOCATE(jdtset_,(0:ndtuj))
283  jdtset_(0:ndtuj)=pack(dtpawuj(:)%iuj,dtpawuj(:)%iuj/=-1)
284  jdtset=maxval(dtpawuj(:)%iuj)
285 
286 !DEBUG
287 !write(message,'(10(a,i3))')'pawuj_det jdtset ',jdtset,&
288 !& ' ndtuj ', ndtuj,' ndtpawuj ',ndtpawuj
289 !call wrtout(std_out,message,'COLL')
290 !call flush_unit(6)
291 !END DEBUG
292 
293  nspden=dtpawuj(jdtset)%nspden
294  nat_org=dtpawuj(jdtset)%nat
295  macro_uj=dtpawuj(jdtset)%macro_uj
296  pawujat=dtpawuj(jdtset)%pawujat
297  pawprtvol=dtpawuj(jdtset)%pawprtvol
298  pawujga=dtpawuj(jdtset)%pawujga
299  pawujoption=dtpawuj(jdtset)%option
300  ph0phiint=dtpawuj(jdtset)%ph0phiint
301  dmatpuopt=dtpawuj(jdtset)%dmatpuopt
302  marr=maxval((/ 9, nspden*nat_org ,nat_org*3 /))
303  eyp=2.5_dp ! for dmatpuopt==2 and 3
304  if (dmatpuopt==1) eyp=eyp+3.0_dp
305  if (dmatpuopt>=3) eyp=(eyp+3.0_dp-dmatpuopt)
306 
307  ABI_ALLOCATE(chih,(ndtpawuj,nat_org))
308  ABI_ALLOCATE(idum2,(marr,0:ndtuj))
309  ABI_ALLOCATE(drarr,(marr,0:ndtuj))
310  ABI_ALLOCATE(magv_org,(nat_org))
311  ABI_ALLOCATE(xred_org,(3,nat_org))
312  ABI_ALLOCATE(chi0_org,(nat_org))
313  ABI_ALLOCATE(chi_org,(nat_org))
314  ABI_ALLOCATE(dparr,(marr,0:ndtuj))
315  ABI_ALLOCATE(dparrr,(marr,0:ndtuj))
316  ABI_ALLOCATE(dqarr,(marr,0:ndtuj))
317  ABI_ALLOCATE(dqarrr,(marr,0:ndtuj))
318  ABI_ALLOCATE(distv_org,(nat_org))
319  ABI_ALLOCATE(narrm,(0:ndtuj))
320  dparr=-one ;  dparrr=-one ;  dqarr=-one ;  dqarrr=-one
321 !DEBUG
322 !write(message,fmt='((a,i3,a))')'pawuj_det init sg'
323 !call wrtout(std_out,message,'COLL')
324 !END DEBUG
325  idum2=1 ; drarr=one
326  chih=zero
327 
328 !write(std_out,*) 'pawuj 02'
329 !###########################################################
330 !### 02. Create the file UJDET.nc
331 
332  if(.false.)write(std_out,*)ujdet_filename ! This is for the abirules
333 
334 !write(std_out,*) 'pawuj 03'
335 !###########################################################
336 !### 03. Write out the Input for UJDET
337 
338  write(message, '(3a)' ) ch10,&
339 & ' # input for ujdet, cut it using ''sed -n "/MARK/,/MARK/p" abi.out > ujdet.in ''------- ',ch10
340  call wrtout(ab_out,message,'COLL')
341 
342  if (ndtuj/=4) then
343    write (hstr,'(I0)') jdtset              ! convert integer to string
344  else
345    hstr=''
346  end if
347  if (ndtuj>=2.or.jdtset==1) then
348    idum2(1,1:ndtuj)=4 !dtpawuj(:)%ndtuj
349    call prttagm(dparr,idum2,ab_out,jdtset_,1,marr,1,narrm,ncid,ndtuj,'ndtset','INT',0)
350  end if
351  idum2(1,0:ndtuj)=pack(dtpawuj(:)%nat,dtpawuj(:)%iuj/=-1)
352  call prttagm(dparr,idum2,ab_out,jdtset_,1,marr,1,narrm,ncid,ndtuj,'nat'//trim(hstr),'INT',0)
353 
354  idum2(1,0:ndtuj)=pack(dtpawuj(:)%nspden,dtpawuj(:)%iuj/=-1)
355  call prttagm(dparr,idum2,ab_out,jdtset_,1,marr,1,narrm,ncid,ndtuj,'nspden'//trim(hstr),'INT',0)
356 
357  idum2(1,0:ndtuj)=pack(dtpawuj(:)%macro_uj,dtpawuj(:)%iuj/=-1)
358  call prttagm(dparr,idum2,ab_out,jdtset_,1,marr,1,narrm,ncid,ndtuj,'macro_uj'//trim(hstr),'INT',0)
359 
360  idum2(1,0:ndtuj)=pack(dtpawuj(:)%pawujat,dtpawuj(:)%iuj/=-1)
361  call prttagm(dparr,idum2,ab_out,jdtset_,1,marr,1,narrm,ncid,ndtuj,'pawujat'//trim(hstr),'INT',0)
362 
363  dparr(1,0:ndtuj)=pack(dtpawuj(:)%pawujga,dtpawuj(:)%iuj/=-1)
364  call prttagm(dparr,idum2,ab_out,jdtset_,1,marr,1,narrm,ncid,ndtuj,'pawujga'//trim(hstr),'DPR',0)
365 
366  dparr(1,0:ndtuj)=pack(dtpawuj(:)%pawujrad,dtpawuj(:)%iuj/=-1)
367  call prttagm(dparr,idum2,ab_out,jdtset_,1,marr,1,narrm,ncid,ndtuj,'pawujrad'//trim(hstr),'DPR',0)
368 
369  dparr(1,0:ndtuj)=pack(dtpawuj(:)%pawrad,dtpawuj(:)%iuj/=-1)
370  call prttagm(dparr,idum2,ab_out,jdtset_,1,marr,1,narrm,ncid,ndtuj,'pawrad'//trim(hstr),'DPR',0)
371 
372  dparr(1,0:ndtuj)=pack(dtpawuj(:)%ph0phiint,dtpawuj(:)%iuj/=-1)
373  call prttagm(dparr,idum2,ab_out,jdtset_,1,marr,1,narrm,ncid,ndtuj,'ph0phiint'//trim(hstr),'DPR',0)
374 
375  idum2(1,0:ndtuj)=pack(dtpawuj(:)%pawprtvol,dtpawuj(:)%iuj/=-1)
376  call prttagm(dparr,idum2,ab_out,jdtset_,1,marr,1,narrm,ncid,ndtuj,'pawprtvol'//trim(hstr),'INT',0)
377 
378  idum2(1,0:ndtuj)=pack(dtpawuj(:)%option,dtpawuj(:)%iuj/=-1)
379  call prttagm(dparr,idum2,ab_out,jdtset_,1,marr,1,narrm,ncid,ndtuj,'pawujopt'//trim(hstr),'INT',0)
380 
381  idum2(1,0:ndtuj)=pack(dtpawuj(:)%dmatpuopt,dtpawuj(:)%iuj/=-1)
382  call prttagm(dparr,idum2,ab_out,jdtset_,1,marr,1,narrm,ncid,ndtuj,'dmatpuopt'//trim(hstr),'INT',0)
383 
384  kdtset=0
385 
386  do idtset=0,ndtpawuj
387 !  DEBUG
388 !  write(message,fmt='((a,i3,a))')'pawuj_det m2, idtset ',idtset,ch10
389 !  call wrtout(std_out,message,'COLL')
390 !  call flush_unit(6)
391 !  END DEBUG
392    if (dtpawuj(idtset)%iuj/=-1) then
393      dparr(1:nspden*nat_org,kdtset)=reshape(dtpawuj(idtset)%vsh,(/nspden*nat_org/))
394      dparrr(1:nspden*nat_org,kdtset)=reshape(dtpawuj(idtset)%occ,(/nspden*nat_org/))
395 !    DEBUG
396 !    write(message,fmt='((a,i3,a))')'pawuj_det m3, idtset ',idtset,ch10
397 !    call wrtout(std_out,message,'COLL')
398 !    write(std_out,*)' marr,narrm,ncid,ndtuj,nat_org,kdtset,idtset=',marr,narrm,ncid,ndtuj,nat_org,kdtset,idtset
399 !    write(std_out,*)' dtpawuj(idtset)%xred=',dtpawuj(idtset)%xred
400 !    call flush_unit(6)
401 !    END DEBUG
402      dqarr(1:nat_org*3,kdtset)=reshape(dtpawuj(idtset)%xred,(/nat_org*3/))
403      dqarrr(1:3*3,kdtset)=reshape(dtpawuj(idtset)%rprimd,(/3*3/))
404      idum2(1:3,kdtset)=reshape(dtpawuj(idtset)%scdim,(/3/))
405      drarr(1:nwfchr,kdtset)=reshape(dtpawuj(idtset)%wfchr,(/nwfchr/))
406 !    DEBUG
407 !    write(message,fmt='((a,i3,a))')'pawuj_det m4, idtset ',idtset,ch10
408 !    call wrtout(std_out,message,'COLL')
409 !    call flush_unit(6)
410 !    END DEBUG
411      kdtset=kdtset+1
412    end if
413  end do
414 
415 !DEBUG
416 !write(message,fmt='((a,i3,a))')'pawuj_det m5'
417 !call wrtout(std_out,message,'COLL')
418 !END DEBUG
419  call prttagm(dparr,idum2,ab_out,jdtset_,2,marr,nspden*nat_org,narrm,ncid,ndtuj,'vsh'//trim(hstr),'DPR',0)
420  call prttagm(dparrr,idum2,ab_out,jdtset_,2,marr,nspden*nat_org,narrm,ncid,ndtuj,'occ'//trim(hstr),'DPR',0)
421 !DEBUG
422 !write(message,fmt='((a,i3,a))')'pawuj_det m6'
423 !call wrtout(std_out,message,'COLL')
424 !END DEBUG
425  call prttagm(dqarr,idum2,ab_out,jdtset_,2,marr,nat_org*3,narrm,ncid,ndtuj,'xred'//trim(hstr),'DPR',0)
426  call prttagm(dqarrr,idum2,ab_out,jdtset_,2,marr,3*3,narrm,ncid,ndtuj,'rprimd'//trim(hstr),'DPR',0)
427  call prttagm(dqarrr,idum2,ab_out,jdtset_,2,marr,3,narrm,ncid,ndtuj,'scdim'//trim(hstr),'INT',0)
428  call prttagm(drarr,idum2,ab_out,jdtset_,2,marr,nwfchr,narrm,ncid,ndtuj,'wfchr'//trim(hstr),'DPR',0)
429  ABI_DEALLOCATE(narrm)
430 
431  write(message, '( 15a )'  ) ch10,' # further possible options: ',ch10,&
432 & ' #    scdim    2 2 2 ',ch10,&
433 & ' #    mdist    10.0  ',ch10,&
434 & ' #  pawujga    2 '    ,ch10,&
435 & ' # pawujopt    2 '    ,ch10,&
436 & ' # pawujrad    3.0'   ,ch10,&
437 & ' # ------- end input for ujdet: end-MARK  -------- ',ch10
438  call wrtout(ab_out,message,'COLL')
439 
440 !write(std_out,*) 'pawuj 04'
441 !###########################################################
442 !### 04. Test
443 
444  if (ndtuj/=4)  return
445 
446  write(message, '(3a)' ) ch10,' ---------- calculate U, (J) start ---------- ',ch10
447  call wrtout(ab_out,message,'COLL')
448 
449  if (all(dtpawuj(1:ndtpawuj)%pawujat==pawujat)) then
450    write (message,fmt='(a,i3)') ' All pawujat  ok and equal to ',pawujat
451    call wrtout(ab_out,message,'COLL')
452  else
453    write (message,fmt='(a,4i3,2a)') ' Differing values of pawujat were found: ',dtpawuj(1:ndtuj)%pawujat,ch10,&
454 &   'No determination of U.'
455    call wrtout(ab_out,message,'COLL')
456    return
457  end if
458 
459  if (all(dtpawuj(1:ndtpawuj)%macro_uj==macro_uj)) then
460    if (nspden==1) then
461      write(message,fmt='(2a)') ' pawuj_det found nspden==1, determination',&
462 &     ' of U-parameter for unpol. struct. (non standard)'
463    else if (macro_uj==1.and.nspden==2) then
464      write(message,fmt='(2a)') ' pawuj_det: found macro_uj=1 and nspden=2:',&
465 &     ' standard determination of U-parameter'
466    else if (macro_uj==2.and.nspden==2) then
467      write(message,fmt='(2a)') ' pawuj_det: found macro_uj=2 and nspden=2:',&
468 &     ' determination of U on single spin channel (experimental)'
469 
470    else if (macro_uj==3.and.nspden==2) then
471      write(message,fmt='(2a)') ' pawuj_det: found macro_uj=3 and nspden=2,',&
472 &     ' determination of J-parameter on single spin channel (experimental)'
473    end if
474    write (message,fmt='(a,i3,a,a)') ' All macro_uj ok and equal to ',macro_uj,ch10,trim(message)
475    write (message,'(a,i3)') ' All macro_uj ok and equal to ',macro_uj
476    call wrtout(ab_out,message,'COLL')
477  else
478    write (message,fmt='(a,10i3)') ' Differing values of macro_uj were found: ',dtpawuj(:)%macro_uj
479    write (message,fmt='(3a)')trim(message),ch10,' No determination of U.'
480    call wrtout(ab_out,message,'COLL')
481    return
482  end if
483 
484  if (macro_uj>1.and.nspden==1) then
485    write (message,'(4a,2a)') ' U on a single spin channel (or J) can only be determined for nspden=2 ,',ch10,&
486 &   'No determination of U.'
487    call wrtout(ab_out,message,'COLL')
488    return
489  end if
490 
491 !Calculation of response matrix
492 
493  do jdtset=1,4
494    if (nspden==1) then
495      chih(jdtset,1:nat_org)=dtpawuj(jdtset)%occ(1,:)
496    else if (macro_uj==1.and.nspden==2) then
497      chih(jdtset,1:nat_org)=dtpawuj(jdtset)%occ(1,:)+dtpawuj(jdtset)%occ(2,:)
498    else if (macro_uj==2.and.nspden==2) then
499      chih(jdtset,1:nat_org)=dtpawuj(jdtset)%occ(1,:)
500    else if (macro_uj==3.and.nspden==2) then
501      chih(jdtset,1:nat_org)=dtpawuj(jdtset)%occ(2,:)
502    end if
503    vsh(jdtset)=dtpawuj(jdtset)%vsh(1,pawujat)
504    if (pawprtvol==3) then
505      write(message,fmt='(2a,i3,a,f15.12)') ch10,' Potential shift vsh(',jdtset,') =',vsh(jdtset)
506      call wrtout(std_out,message,'COLL')
507      write(message,fmt='( a,i3,a,120f15.9)') ' Occupations occ(',jdtset,') ',chih(jdtset,1:nat_org)
508      call wrtout(std_out,message,'COLL')
509    end if
510  end do
511 
512  if (any(abs((/(vsh(ii)-vsh(ii+2), ii=1,2) /))<0.00000001)) then
513    write(message, '(2a,18f10.7,a)' )  ch10,' vshift is too small: ',abs((/(vsh(ii)-vsh(ii+2), ii=1,2) /))
514    call wrtout(ab_out,message,'COLL')
515    return
516  end if
517 
518 !DEBUG
519 !write(message,fmt='(a)')'pawuj_det: after test vsh'
520 !call wrtout(std_out,message,'COLL')
521 !END DEBUG
522 
523  chi0_org=(chih(1,1:nat_org)-chih(3,1:nat_org))/(vsh(1)-vsh(3))/dtpawuj(1)%diemix
524  chi_org=(chih(2,1:nat_org)-chih(4,1:nat_org))/(vsh(2)-vsh(4))
525 
526  if (pawprtvol==3) then
527    write(message,fmt='(2a, 150f15.10)') ch10,' Chi_0n ',chi0_org
528    call wrtout(std_out,message,'COLL')
529    write(message,fmt='(a, 150f15.10)') ' Chi_n ',chi_org
530    call wrtout(std_out,message,'COLL')
531  end if
532 
533  write(message,fmt='(a)')': '
534  if (nspden==2) then
535    magv_org=dtpawuj(1)%occ(1,:)-dtpawuj(1)%occ(2,:)
536    if (all(abs(magv_org)<0.001)) then
537      magv_org=(/(1,im1=1,nat_org)/)
538    else
539      magv_org=abs(magv_org)/magv_org
540      if (magv_org(1).lt.0) magv_org=magv_org*(-1_dp)
541      if (all(magv_org(2:nat_org).lt.0)) then
542        magv_org=abs(magv_org)
543        write(message,'(a)')', (reset to fm): '
544      end if
545    end if
546  else
547    magv_org=(/(1,im1=1,nat_org)/)
548  end if
549 
550  if (pawprtvol==3) then
551    write(message,fmt='(3a, 150f4.1)') ch10,' Magnetization',trim(message),magv_org
552    call wrtout(std_out,message,'COLL')
553  end if
554 
555 !Case of extrapolation to larger r_paw: calculate intg
556 
557  if (all(dtpawuj(1)%wfchr(:)/=0).and.ph0phiint/=1) then
558    if (dtpawuj(1)%pawujrad<20.and.dtpawuj(1)%pawujrad>dtpawuj(1)%pawrad) then
559      fcorr=(1-ph0phiint)/(IRadFnH(dtpawuj(1)%pawrad,20.0_dp,nint(dtpawuj(1)%wfchr(2)),&
560 &     nint(dtpawuj(1)%wfchr(3)),dtpawuj(1)%wfchr(1)))
561      intg=ph0phiint/(1-fcorr*IRadFnH(dtpawuj(1)%pawujrad,20.0_dp,nint(dtpawuj(1)%wfchr(2)),&
562 &     nint(dtpawuj(1)%wfchr(3)),dtpawuj(1)%wfchr(1)))
563      write(message, fmt='(a,f12.5,a,f12.5)') ' pawuj_det: met2 extrapolation to ', dtpawuj(1)%pawujrad,' using factor ',intg
564      call wrtout(std_out,message,'COLL')
565    else if (dtpawuj(1)%pawujrad<dtpawuj(1)%pawrad) then
566      a=0 ; a(1:3)=dtpawuj(1)%wfchr(4:6)
567      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/)
568      b=(/(dtpawuj(1)%pawujrad**(im1)-dtpawuj(1)%pawrad**(im1),im1=1,5)/)
569      intg=dot_product(a,b)
570      intg=ph0phiint/(ph0phiint+intg)
571      write(message, fmt='(a,f12.5,a,f12.5)') ' pawuj_det: met1 extrapolation to ', dtpawuj(1)%pawujrad,' using factor ',intg
572      call wrtout(std_out,message,'COLL')
573    else if (dtpawuj(1)%pawujrad==dtpawuj(1)%pawrad) then
574      intg=one
575      write(message, fmt='(a,2i7,3f12.5)') ' pawuj_det: no extrapolation (pawujrad=pawrad)'
576      call wrtout(std_out,message,'COLL')
577    else
578      intg=ph0phiint
579      write(message, fmt='(a,3f12.5)') ' pawuj_det: extrapolation to r->\inf using factor ',intg
580      call wrtout(std_out,message,'COLL')
581    end if
582  else
583    write(message, fmt='(a,2i7,3f12.5)') ' pawuj_det: no extrapolation (ph0phiint=1 or wfchr=0)'
584    call wrtout(std_out,message,'COLL')
585    intg=one
586  end if
587 
588 
589 !Determine U in primitive cell
590 
591  write(message,fmt='(a)')' pawuj_det: determine U in primitive cell'
592  call wrtout(std_out,message,'COLL')
593 
594  call lcalcu(int(magv_org),nat_org,dtpawuj(1)%rprimd,dtpawuj(1)%xred,chi_org,chi0_org,pawujat,ures,pawprtvol,pawujga,pawujoption)
595 
596 !Begin calculate U in supercell
597 
598 !Analize shell structure of primitive cell
599 !and atomic distances in primitive cell
600  ABI_ALLOCATE(smult_org,(nat_org))
601  ABI_ALLOCATE(sdistv_org,(nat_org))
602  call shellstruct(dtpawuj(1)%xred,dtpawuj(1)%rprimd,nat_org,&
603 & int(magv_org),distv_org,smult_org,sdistv_org,nsh_org,pawujat,pawprtvol)
604 
605  ii=1
606  write(message, fmt='(8a)') ' URES ','     ii','    nat','       r_max','    U(J)[eV]','   U_ASA[eV]','   U_inf[eV]',ch10
607  write(message, fmt='(a,2i7,4f12.5)') trim(message)//' URES ',ii,nat_org,maxval(abs(distv_org)),ures,ures*exp(log(intg)*eyp),&
608 & ures*exp(log(ph0phiint)*eyp)
609  call wrtout(std_out,message,'COLL')
610  call wrtout(ab_out,message,'COLL')
611 
612  if (pawprtvol>1) then
613    write(message,fmt='(a, 150f10.5)')' pawuj_det: ionic distances in original cell (distv_org) ', distv_org
614    call wrtout(std_out,message,'COLL')
615  end if
616 
617 !Construct supercell, calculate limit dimensions of supercell
618  ii=1
619  maxnat=product(dtpawuj(1)%scdim(:))*nat_org
620  if (maxnat==0) then
621    maxnat=dtpawuj(1)%scdim(1)
622    ext=(/ii, ii, ii/)
623  else
624    jj=1
625    do while (jj<=3)
626      ext(jj)=minval( (/ii, dtpawuj(1)%scdim(jj) /) )
627      jj=jj+1
628    end do
629  end if
630  ext=ext+(/ 1,1,1 /)
631  ii=ii+1
632 
633 
634  nat_sc=product(ext)*nat_org
635 
636 !DEBUG
637 !write(message,fmt='(a,3i4)')'pawuj_det: ext ',ext
638 !call wrtout(std_out,message,'COLL')
639 !END DEBUG
640 
641  do while (nat_sc<=maxnat)
642    ABI_ALLOCATE(chi0_sc,(nat_sc))
643    ABI_ALLOCATE(chi_sc,(nat_sc))
644    ABI_ALLOCATE(distv_sc,(nat_sc))
645    ABI_ALLOCATE(magv_sc,(nat_sc))
646    ABI_ALLOCATE(sdistv_sc,(nat_sc))
647    ABI_ALLOCATE(smult_sc,(nat_sc))
648    ABI_ALLOCATE(xred_sc,(3,nat_sc))
649 
650 !  Determine positions=xred_sc and supercell dimensions=rpimd_sc
651    call mksupercell(dtpawuj(1)%xred,int(magv_org),dtpawuj(1)%rprimd,nat_org,&
652 &   nat_sc,xred_sc,magv_sc,rprimd_sc,ext,pawprtvol)
653 
654 !  Determine shell structure of supercell: multiplicities (smult_sc), radii of shells (sdistv_sc)
655 !  number of shells (nsh_sc) and atomic distances in supercell (distv_sc)
656 
657    write(message,fmt='(a,3i3,a)')' pawuj_det: determine shell structure of ',ext(1:3),' supercell'
658    call wrtout(std_out,message,'COLL')
659 
660    call shellstruct(xred_sc,rprimd_sc,nat_sc,int(magv_sc),distv_sc,smult_sc,sdistv_sc,nsh_sc,pawujat,pawprtvol)
661 
662    if (pawprtvol>=2) then
663      write(message,fmt='(a)')' pawuj_det: ionic distances in supercell (distv_sc) '
664      call wrtout(std_out,message,'COLL')
665      call prmat(distv_sc,1,nat_sc,1,std_out)
666    end if
667 
668 !  Determine chi and chi0 in supercell (chi0_sc, chi_sc)
669 !  DEBUG
670 !  write(message,fmt='(a)')'pawuj_det:  chi and chi0 in supercell'
671 !  call wrtout(std_out,message,'COLL')
672 !  END DEBUG
673 
674    if (pawujoption>2) then
675      invopt=2
676    else
677      invopt=1
678    end if
679 
680    call chiscwrt(chi0_org,distv_org,nat_org,sdistv_org,smult_org,nsh_org,&
681 &   chi0_sc,distv_sc,nat_sc,smult_sc,nsh_sc,invopt,pawprtvol)
682    call chiscwrt(chi_org,distv_org,nat_org,sdistv_org,smult_org,nsh_org,&
683 &   chi_sc,distv_sc,nat_sc,smult_sc,nsh_sc,invopt,pawprtvol)
684 
685 !  Calculate U in supercell
686 !  DEBUG
687 !  write(message,fmt='(a)')'pawuj_det:   U in supercell'
688 !  call wrtout(std_out,message,'COLL')
689 !  END DEBUG
690    call lcalcu(int(magv_sc),nat_sc,rprimd_sc,xred_sc,chi_sc,chi0_sc,pawujat,ures,pawprtvol,pawujga,pawujoption)
691 
692    write(message, fmt='(a,2i7,4f12.5)') ' URES ',ii,nat_sc,maxval(abs(distv_sc)),ures,ures*exp(log(intg)*eyp),&
693 &   ures*exp(log(ph0phiint)*eyp)
694    call wrtout(std_out,message,'COLL')
695    call wrtout(ab_out,message,'COLL')
696 
697    ABI_DEALLOCATE(chi0_sc)
698    ABI_DEALLOCATE(chi_sc)
699    ABI_DEALLOCATE(distv_sc)
700    ABI_DEALLOCATE(magv_sc)
701    ABI_DEALLOCATE(sdistv_sc)
702    ABI_DEALLOCATE(smult_sc)
703    ABI_DEALLOCATE(xred_sc)
704 
705    if (product(dtpawuj(1)%scdim(:))==0) then
706      ext=(/ii, ii, ii/)
707    else
708      jj=1
709      do while (jj<=3)
710        ext(jj)=minval( (/ii, dtpawuj(1)%scdim(jj) /) )
711        jj=jj+1
712      end do
713    end if
714    ext=ext+(/ 1,1,1 /)
715    ii=ii+1
716 
717    nat_sc=product(ext)*nat_org
718 
719  end do
720 
721  write(message, '(3a)' )ch10,' ---------- calculate U, (J) end -------------- ',ch10
722  call wrtout(ab_out,message,'COLL')
723 
724  ABI_DEALLOCATE(dparrr)
725  ABI_DEALLOCATE(dqarr)
726  ABI_DEALLOCATE(dqarrr)
727  ABI_DEALLOCATE(jdtset_)
728  ABI_DEALLOCATE(chi_org)
729  ABI_DEALLOCATE(chi0_org)
730  ABI_DEALLOCATE(smult_org)
731  ABI_DEALLOCATE(sdistv_org)
732  ABI_DEALLOCATE(chih)
733  ABI_DEALLOCATE(idum2)
734  ABI_DEALLOCATE(drarr)
735  ABI_DEALLOCATE(dparr)
736  ABI_DEALLOCATE(magv_org)
737  ABI_DEALLOCATE(xred_org)
738  ABI_DEALLOCATE(distv_org)
739 
740  DBG_EXIT("COLL")
741 
742 end subroutine pawuj_det

m_paw_uj/pawuj_free [ Functions ]

[ Top ] [ m_paw_uj ] [ Functions ]

NAME

  pawuj_free

FUNCTION

   deallocate pawuj stuff

INPUTS

OUTPUT

PARENTS

      pawuj_drive,ujdet

CHILDREN

SOURCE

166 subroutine pawuj_free(dtpawuj)
167 
168 
169 !This section has been created automatically by the script Abilint (TD).
170 !Do not modify the following lines by hand.
171 #undef ABI_FUNC
172 #define ABI_FUNC 'pawuj_free'
173 !End of the abilint section
174 
175  implicit none
176 
177 !Arguments -------------------------------
178  type(macro_uj_type),intent(inout) :: dtpawuj
179 
180 ! *********************************************************************
181 
182  if (allocated(dtpawuj%scdim))    then
183    ABI_DEALLOCATE(dtpawuj%scdim)
184  end if
185  if (allocated(dtpawuj%occ))      then
186    ABI_DEALLOCATE(dtpawuj%occ)
187  end if
188  if (allocated(dtpawuj%rprimd))   then
189    ABI_DEALLOCATE(dtpawuj%rprimd)
190  end if
191  if (allocated(dtpawuj%vsh))      then
192    ABI_DEALLOCATE(dtpawuj%vsh)
193  end if
194  if (allocated(dtpawuj%xred))     then
195    ABI_DEALLOCATE(dtpawuj%xred)
196  end if
197  if (allocated(dtpawuj%wfchr))    then
198    ABI_DEALLOCATE(dtpawuj%wfchr)
199  end if
200  if (allocated(dtpawuj%zioneff))  then
201    ABI_DEALLOCATE(dtpawuj%zioneff)
202  end if
203 
204 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)

PARENTS

      pawuj_drive,ujdet

CHILDREN

SOURCE

 87 subroutine pawuj_ini(dtpawuj,ndtset)
 88 
 89 
 90 !This section has been created automatically by the script Abilint (TD).
 91 !Do not modify the following lines by hand.
 92 #undef ABI_FUNC
 93 #define ABI_FUNC 'pawuj_ini'
 94 !End of the abilint section
 95 
 96  implicit none
 97 
 98 !Arguments ------------------------------------
 99 !scalars
100  integer                           :: ndtset
101  type(macro_uj_type),intent(inout) :: dtpawuj(0:ndtset)
102 
103 !Local variables -------------------------
104 !Variables for partial dos calculation
105 !scalars
106  integer, parameter :: nwfchr=6
107  integer            :: iuj,im1
108 ! *********************************************************************
109 
110  DBG_ENTER("COLL")
111 
112  do iuj=0,ndtset
113    !write(std_out,*)'pawuj_ini iuj ',iuj
114    dtpawuj(iuj)%diemix=0.45_dp
115    dtpawuj(iuj)%iuj=0
116    dtpawuj(iuj)%nat=0
117    dtpawuj(iuj)%ndtset=1
118    dtpawuj(iuj)%nspden=1
119    dtpawuj(iuj)%macro_uj=0
120    dtpawuj(iuj)%option=1
121    dtpawuj(iuj)%pawujat=1
122    dtpawuj(iuj)%pawujga=one
123    dtpawuj(iuj)%pawprtvol=1
124    dtpawuj(iuj)%ph0phiint=one
125    dtpawuj(iuj)%dmatpuopt=3
126    dtpawuj(iuj)%pawujrad=3.0_dp
127    dtpawuj(iuj)%pawrad=20.0_dp
128    !Allocate arrays
129    !write(std_out,*)'pawuj_ini before arrays'
130    ABI_ALLOCATE(dtpawuj(iuj)%rprimd,(3,3))
131    ABI_ALLOCATE(dtpawuj(iuj)%scdim,(3))
132    ABI_ALLOCATE(dtpawuj(iuj)%wfchr,(nwfchr))
133    dtpawuj(iuj)%rprimd=reshape((/ 1,0,0,0,1,0,0,0,1/),(/ 3,3 /))
134    dtpawuj(iuj)%scdim=reshape((/ 250,0,0 /),(/3 /))
135    dtpawuj(iuj)%wfchr=(/ (0,im1=1,nwfchr) /)
136    if (iuj>0) then
137      dtpawuj(iuj)%iuj=-1
138    end if
139  end do
140 
141  DBG_EXIT("COLL")
142 
143 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

  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

OUTPUT

  dtpawuj(0:ndtpawuj) (initialization of fields vsh, occ, iuj,nnat)

PARENTS

      scfcv

CHILDREN

      free_my_atmtab,get_my_atmtab,linvmat,prmat,wrtout,xmpi_lor,xmpi_sum

SOURCE

775 subroutine pawuj_red(dtset,dtpawuj,fatvshift,my_natom,natom,ntypat,paw_ij,pawrad,pawtab,ndtpawuj,&
776 &                    mpi_atmtab,comm_atom) ! optional arguments (parallelism)
777 
778 
779 !This section has been created automatically by the script Abilint (TD).
780 !Do not modify the following lines by hand.
781 #undef ABI_FUNC
782 #define ABI_FUNC 'pawuj_red'
783 !End of the abilint section
784 
785  implicit none
786 
787 !Arguments ------------------------------------
788 !scalars
789  integer,intent(in)                 :: my_natom,natom,ntypat,ndtpawuj
790  integer,optional,intent(in)        :: comm_atom
791  real(dp),intent(in)                :: fatvshift
792 !arrays
793  integer,optional,target,intent(in) :: mpi_atmtab(:)
794  type(paw_ij_type),intent(in)       :: paw_ij(my_natom)
795  type(pawtab_type),intent(in)       :: pawtab(ntypat)
796  type(pawrad_type),intent(in)       :: pawrad(ntypat)
797  type(dataset_type),intent(in)      :: dtset
798  type(macro_uj_type),intent(inout)  :: dtpawuj(0:ndtpawuj)
799 
800 !Local variables-------------------------------
801 !scalars
802  integer,parameter           :: natmax=2,ncoeff=3
803  integer                     :: iatom,iatom_tot,ierr,im1,im2,ispden,itypat,ll,nspden,nsppol,iuj
804  integer                     :: my_comm_atom,nnat,natpawu,natvshift,pawujat,ndtset,typawujat
805  logical                     :: usepawu !antiferro,
806  logical                     :: my_atmtab_allocated,paral_atom
807  character(len=1000)         :: message,hstr
808  character(len=500)          :: messg
809 !arrays
810  logical                     :: musk(3,natom)
811  integer,pointer             :: my_atmtab(:)
812  real(dp)                    :: rrtab(ncoeff),wftab(ncoeff),a(ncoeff,ncoeff),b(ncoeff,ncoeff)! ,a(ncoeff,ncoeff)
813  real(dp),allocatable        :: nnocctot(:,:) !,magv(:)
814  real(dp),allocatable        :: atvshift(:,:,:) ! atvshift(natvshift,2,natom)
815  logical,allocatable         :: dmusk(:,:),atvshmusk(:,:,:) !atvshmusk(natvshift,2,natom)
816 
817 ! *********************************************************************
818 
819 !Initializations
820  nspden=1;nsppol=1
821  if (my_natom>0) then
822    nspden=paw_ij(1)%nspden ; nsppol=paw_ij(1)%nsppol
823  end if
824 
825  natvshift=dtset%natvshift
826  pawujat=dtset%pawujat
827  natpawu=dtset%natpawu   ; ndtset=dtset%ndtset
828  ABI_ALLOCATE(atvshift,(natvshift,nspden,natom))
829  ABI_ALLOCATE(atvshmusk,(natvshift,nspden,natom))
830  ABI_ALLOCATE(dmusk,(nspden,natom))
831  ABI_ALLOCATE(nnocctot,(nspden,natom))
832  musk=.false.; dmusk=.false.
833  atvshift=fatvshift*dtset%atvshift
834  atvshmusk=.false.
835  nnocctot=zero
836  typawujat=dtset%typat(pawujat)
837  usepawu=(count(pawtab(:)%usepawu>0)>0)
838 
839 !Set up parallelism over atoms
840  paral_atom=(present(comm_atom).and.(my_natom/=natom))
841  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
842  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
843  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
844 
845  nnat=0
846  if (usepawu) then
847    write(message,'(3a)') ch10, '---------- pawuj_red ------ ',ch10
848    call wrtout(std_out,  message,'COLL');
849    do iatom=1,my_natom
850      iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
851      itypat=paw_ij(iatom)%itypat;ll=pawtab(itypat)%lpawu
852      if ((ll>=0).and.(pawtab(itypat)%usepawu>0).and.itypat==typawujat) then
853        musk(:,iatom_tot)=(/.true., .true., .true. /)
854        atvshmusk(:,:,iatom_tot)=reshape((/ (( (im1==1), im1=1,natvshift)  ,im2=1,nspden ) /),(/natvshift,nspden/))
855        do ispden=1,nspden
856          nnocctot(ispden,iatom_tot)=paw_ij(iatom)%nocctot(ispden)
857          dmusk(ispden,iatom_tot)=.true.
858        end do
859        nnat=nnat+1
860      end if
861    end do
862 
863 !  Reduction in case of parallelism
864    if (paral_atom) then
865      call xmpi_sum(nnocctot ,my_comm_atom,ierr)
866      call xmpi_lor(atvshmusk,my_comm_atom) ! dim=natpawu ???
867      call xmpi_lor(dmusk    ,my_comm_atom)
868      call xmpi_lor(musk     ,my_comm_atom)
869      call xmpi_sum(nnat     ,my_comm_atom,ierr)
870    end if
871 
872    iuj=maxval(dtpawuj(:)%iuj)
873    !write(std_out,*)' pawuj_red: iuj',iuj
874    !write(std_out,*)' pawuj_red: dtpawuj(:)%iuj ',dtpawuj(:)%iuj
875 
876    if (iuj==1.or.iuj==3) then  ! 1 and 3: non-scf steps
877      dtpawuj(iuj+1)%iuj=iuj+1
878    end if
879 
880    ! TODO: check that this is correct: this point is passed several times
881    ! for a given value of iuj - should the stuff be accumulated instead of replaced?
882    if(allocated(dtpawuj(iuj)%vsh))  then
883      ABI_DEALLOCATE(dtpawuj(iuj)%vsh)
884    end if
885    ABI_ALLOCATE(dtpawuj(iuj)%vsh,(nspden,nnat))
886    if(allocated(dtpawuj(iuj)%occ))  then
887      ABI_DEALLOCATE(dtpawuj(iuj)%occ)
888    end if
889    ABI_ALLOCATE(dtpawuj(iuj)%occ,(nspden,nnat))
890    if(allocated(dtpawuj(iuj)%xred))  then
891      ABI_DEALLOCATE(dtpawuj(iuj)%xred)
892    end if
893    ABI_ALLOCATE(dtpawuj(iuj)%xred,(3,nnat))
894 
895    if (iuj==1) then
896      ABI_ALLOCATE(dtpawuj(0)%vsh,(nspden,nnat))
897      ABI_ALLOCATE(dtpawuj(0)%occ,(nspden,nnat))
898      ABI_ALLOCATE(dtpawuj(0)%xred,(3,nnat))
899      dtpawuj(0)%vsh=0
900      dtpawuj(0)%occ=0
901      dtpawuj(0)%xred=0
902    end if
903 
904    rrtab=(/0.75_dp,0.815_dp,1.0_dp/)*pawtab(typawujat)%rpaw
905    wftab=pawtab(typawujat)%phi(pawtab(typawujat)%mesh_size,pawtab(typawujat)%lnproju(1))
906 
907 !  DEBUG
908 !  write(std_out,*)' pawuj_red: rrtab ',rrtab
909 !  write(std_out,*)' pawuj_red: wftab ',wftab
910 !  END DEBUG
911 
912    do im1=1,ncoeff
913 !    DEBUG
914 !    write(std_out,*)' pawuj_red: ncoeff ',ncoeff,' im1 ',im1
915 !    END DEBUG
916      if (pawrad(typawujat)%mesh_type==1) then
917        im2=nint(rrtab(im1)/pawrad(typawujat)%rstep+1)
918      else if (pawrad(typawujat)%mesh_type==2) then
919        im2=nint(log(rrtab(im1)/pawrad(typawujat)%rstep+1)/pawrad(typawujat)%lstep+1)
920      else if (pawrad(typawujat)%mesh_type==3) then
921        im2=nint(log(rrtab(im1)/pawrad(typawujat)%rstep)/pawrad(typawujat)%lstep+1)
922      else if (pawrad(typawujat)%mesh_type==4) then
923        im2=nint(pawtab(typawujat)%mesh_size*(1-exp((-one)*rrtab(im1)/pawrad(typawujat)%rstep))+1)
924      end if
925 
926 !    DEBUG
927 !    write(std_out,*)' pawuj_red: im2 ',im2
928 !    END DEBUG
929 
930      rrtab(im1)=pawrad(typawujat)%rad(im2)
931      wftab(im1)=pawtab(typawujat)%phi(im2,pawtab(typawujat)%lnproju(1))
932    end do
933    write(message,fmt='(a,i3,a,10f10.5)')' pawuj_red: mesh_type',pawrad(typawujat)%mesh_type,' rrtab:', rrtab
934    call wrtout(std_out,message,'COLL')
935    write(message,fmt='(a,10f10.5)')' pawuj_red: wftab', wftab
936    call wrtout(std_out,message,'COLL')
937    a=reshape((/ (( rrtab(im2)**(im1-1), im1=1,3)  ,im2=1,3 )/),(/ncoeff,ncoeff/))
938    write(messg,fmt='(a)')'A'
939    call linvmat(a,b,ncoeff,messg,2,0.0_dp,3) ! linvmat(inmat,oumat,nat,nam,option,gam,prtvol)
940    write(std_out,*) 'pawuj_red: a,b ', a,b
941    wftab=matmul(wftab,b)
942    write(std_out,*) 'pawuj_red: wftab ', wftab
943    dtpawuj(iuj)%wfchr(4:6)=wftab
944 
945    dtpawuj(iuj)%nat=nnat
946    write(std_out,*) 'pawuj_red: m1'
947    dtpawuj(iuj)%vsh=reshape(pack(atvshift,atvshmusk),(/ nspden,nnat /))
948 !  factor in next line to compensate nocctot contains just occ of 1 spin channel for nspden=1
949    write(std_out,*) 'pawuj_red: m2'
950    dtpawuj(iuj)%occ=reshape(pack(nnocctot,dmusk),(/nspden,nnat/))*(3-nspden)
951    write(std_out,*) 'pawuj_red: m3'
952 !  dtpawuj(iuj)%occ=dtpawuj(iuj)%occ/pawtab(typawujat)%ph0phiint(1)
953 
954    write(std_out,*) 'pawuj_red: occ ', dtpawuj(iuj)%occ
955 
956    dtpawuj(iuj)%xred=reshape(pack(dtset%xred_orig(:,:,1),musk),(/3,nnat/))
957    dtpawuj(iuj)%ph0phiint=pawtab(typawujat)%ph0phiint(1)
958    dtpawuj(iuj)%wfchr(1:3)=(/ pawtab(typawujat)%zioneff(1)*(dtset%lpawu(typawujat)+2),&
959 &   one*(dtset%lpawu(typawujat)+1),one*(dtset%lpawu(typawujat))/)
960    dtpawuj(iuj)%pawrad=pawtab(typawujat)%rpaw
961 
962    write(std_out,*) 'pawuj_red: wfchr ',dtpawuj(iuj)%wfchr
963 
964 
965    write (hstr,'(I0)') iuj
966    write(message,'(a,a,I3,I3,a)') ch10, '---------- MARK ------ ',iuj,maxval(dtpawuj(:)%iuj) ,ch10
967    call wrtout(std_out,message,'COLL')
968    write(message,fmt='(a)') 'vsh'//trim(hstr)
969    call wrtout(std_out,message,'COLL')
970    call prmat(dtpawuj(iuj)%vsh(:,:),1,nnat*nspden,1)
971    write(message,fmt='(a)') 'occ'//trim(hstr)
972    call wrtout(std_out,message,'COLL')
973    call prmat(dtpawuj(iuj)%occ(:,:),1,nnat*nspden,1)
974    write(message, '(3a)' )'---------- MARK ---------- ',ch10
975    call wrtout(std_out,message,'COLL')
976  end if !usepawu
977 
978  ABI_DEALLOCATE(nnocctot)
979  ABI_DEALLOCATE(dmusk)
980  ABI_DEALLOCATE(atvshift)
981  ABI_DEALLOCATE(atvshmusk)
982 
983 !Destroy atom table used for parallelism
984  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
985 
986 end subroutine pawuj_red