TABLE OF CONTENTS


ABINIT/chiscwrt [ Functions ]

[ Top ] [ Functions ]

NAME

  chiscwrt

FUNCTION

  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 =

COPYRIGHT

  Copyright (C) 1998-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 .
  For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

OUTPUT

  chi_sc = first line of reponse matrix in super-cell

SIDE EFFECTS

PARENTS

      pawuj_det

CHILDREN

      prmat,wrtout

SOURCE

606 subroutine chiscwrt(chi_org,disv_org,nat_org,sdisv_org,smult_org,nsh_org,chi_sc,&
607 & disv_sc,nat_sc,smult_sc,nsh_sc,opt,prtvol)
608 
609  use defs_basis
610  use m_profiling_abi
611 
612  use m_pptools,    only : prmat
613 
614 !This section has been created automatically by the script Abilint (TD).
615 !Do not modify the following lines by hand.
616 #undef ABI_FUNC
617 #define ABI_FUNC 'chiscwrt'
618  use interfaces_14_hidewrite
619 !End of the abilint section
620 
621  implicit none
622 
623 !Arguments ------------------------------------
624 !scalars
625  integer,intent(in)              :: nat_org,nsh_org,nsh_sc
626  integer,intent(in)              :: nat_sc
627  integer,intent(in),optional     :: prtvol
628  integer,intent(in),optional     :: opt
629 !arrays
630  real(dp),intent(in)             :: chi_org(nat_org),disv_org(nat_org),sdisv_org(nsh_org)
631  integer,intent(in)              :: smult_org(nsh_org), smult_sc(nsh_sc)
632  real(dp),intent(in)             :: disv_sc(nat_sc)
633  real(dp),intent(out)            :: chi_sc(nat_sc)
634 
635 !Local variables-------------------------------
636 !scalars
637  integer                      :: iatom,jatom,jsh,optt
638  character(len=500)           :: message
639 !arrays
640  real(dp)                     :: chi_orgl(nat_org)
641 
642 ! *************************************************************************
643 
644  if (present(opt)) then
645    optt=opt
646  else
647    optt=1
648  end if
649 
650 
651  do iatom=1,nat_org
652    do jsh=1,nsh_org
653      if (disv_org(iatom)==sdisv_org(jsh)) then
654        if (opt==2) then
655          chi_orgl(iatom)=chi_org(iatom)
656        else if (opt==1) then
657          chi_orgl(iatom)=chi_org(iatom)*smult_org(jsh)/smult_sc(jsh)
658        end if
659        exit
660      end if
661    end do !iatom
662  end do  !jsh
663 
664  if (prtvol>1) then
665    write(message,fmt='(a,150f10.5)')' chiscwrt: chi at input ',chi_org
666    call wrtout(std_out,message,'COLL')
667    write(message,fmt='(a,150f10.5)')' chiscwrt: chi after division ',chi_orgl
668    call wrtout(std_out,message,'COLL')
669  end if
670 
671  do iatom=1,nat_sc
672    do jatom=1,nat_org
673      if (disv_org(jatom)==disv_sc(iatom)) then
674        chi_sc(iatom)=chi_orgl(jatom)
675        exit
676      else if (jatom==nat_org) then
677        chi_sc(iatom)=0_dp
678      end if
679    end do
680  end do
681 
682  if (prtvol>1) then
683    write(message,'(a)')' chiscwrt, chi_sc '
684    call wrtout(std_out,message,'COLL')
685    call prmat(chi_sc,1,nat_sc,1,std_out)
686  end if
687 
688 end subroutine chiscwrt

ABINIT/pawuj_det [ Functions ]

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

COPYRIGHT

 Copyright (C) 1998-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 .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

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

 34 #if defined HAVE_CONFIG_H
 35 #include "config.h"
 36 #endif
 37 
 38 #include "abi_common.h"
 39 
 40 subroutine pawuj_det(dtpawuj,ndtpawuj,ujdet_filename,ures)
 41 
 42  use defs_basis
 43  use defs_abitypes
 44  use defs_parameters
 45  use m_profiling_abi
 46  use m_errors
 47 
 48  use m_special_funcs,  only : iradfnh
 49  use m_pptools,        only : prmat
 50  use m_geometry,       only : shellstruct
 51 
 52 !This section has been created automatically by the script Abilint (TD).
 53 !Do not modify the following lines by hand.
 54 #undef ABI_FUNC
 55 #define ABI_FUNC 'pawuj_det'
 56  use interfaces_14_hidewrite
 57  use interfaces_41_geometry
 58  use interfaces_57_iovars
 59  use interfaces_65_paw, except_this_one => pawuj_det
 60 !End of the abilint section
 61 
 62  implicit none
 63 
 64 !Arguments ------------------------------------
 65 !scalars
 66 !arrays
 67  integer                        :: ndtpawuj
 68  type(macro_uj_type),intent(in) :: dtpawuj(0:ndtpawuj)
 69  real(dp),intent(out)           :: ures
 70  character(len=*),intent(in)    :: ujdet_filename
 71 
 72 !Local variables-------------------------------
 73 !scalars
 74  integer,parameter           :: natmax=2,nwfchr=6
 75  integer                     :: ii,jj,nat_org,jdtset,nspden,macro_uj,kdtset,marr
 76  integer                     :: im1,ndtuj,idtset, nsh_org, nsh_sc,nat_sc,maxnat
 77  integer                     :: pawujat,pawprtvol,pawujoption
 78  integer                     :: dmatpuopt,invopt
 79  real(dp)                    :: pawujga,ph0phiint,intg,fcorr,eyp
 80 
 81  character(len=500)          :: message
 82  character(len=2)            :: hstr
 83 !arrays
 84  integer                     :: ext(3)
 85  real(dp)                    :: rprimd_sc(3,3),vsh(ndtpawuj),a(5),b(5)
 86  integer,allocatable         :: narrm(:)
 87  integer,allocatable         :: idum2(:,:),jdtset_(:),smult_org(:),smult_sc(:)
 88  real(dp),allocatable        :: chih(:,:),dqarr(:,:),dqarrr(:,:),dparr(:,:),dparrr(:,:),xred_org(:,:),drarr(:,:)
 89  real(dp),allocatable        :: magv_org(:),magv_sc(:),chi_org(:),chi0_org(:),chi0_sc(:), chi_sc(:), xred_sc(:,:)
 90  real(dp),allocatable        :: sdistv_org(:),sdistv_sc(:),distv_org(:),distv_sc(:)
 91  integer                     :: ncid=0
 92 ! *********************************************************************
 93 
 94  DBG_ENTER("COLL")
 95 
 96 !write(std_out,*) 'pawuj 01'
 97 !###########################################################
 98 !### 01. Allocations
 99 
100 !Initializations
101  ndtuj=count(dtpawuj(:)%iuj/=-1)-1 ! number of datasets initialized by pawuj_red
102  ABI_ALLOCATE(jdtset_,(0:ndtuj))
103  jdtset_(0:ndtuj)=pack(dtpawuj(:)%iuj,dtpawuj(:)%iuj/=-1)
104  jdtset=maxval(dtpawuj(:)%iuj)
105 
106 !DEBUG
107 !write(message,'(10(a,i3))')'pawuj_det jdtset ',jdtset,&
108 !& ' ndtuj ', ndtuj,' ndtpawuj ',ndtpawuj
109 !call wrtout(std_out,message,'COLL')
110 !call flush_unit(6)
111 !END DEBUG
112 
113  nspden=dtpawuj(jdtset)%nspden
114  nat_org=dtpawuj(jdtset)%nat
115  macro_uj=dtpawuj(jdtset)%macro_uj
116  pawujat=dtpawuj(jdtset)%pawujat
117  pawprtvol=dtpawuj(jdtset)%pawprtvol
118  pawujga=dtpawuj(jdtset)%pawujga
119  pawujoption=dtpawuj(jdtset)%option
120  ph0phiint=dtpawuj(jdtset)%ph0phiint
121  dmatpuopt=dtpawuj(jdtset)%dmatpuopt
122  marr=maxval((/ 9, nspden*nat_org ,nat_org*3 /))
123  eyp=2.5_dp ! for dmatpuopt==2 and 3
124  if (dmatpuopt==1) eyp=eyp+3.0_dp
125  if (dmatpuopt>=3) eyp=(eyp+3.0_dp-dmatpuopt)
126 
127  ABI_ALLOCATE(chih,(ndtpawuj,nat_org))
128  ABI_ALLOCATE(idum2,(marr,0:ndtuj))
129  ABI_ALLOCATE(drarr,(marr,0:ndtuj))
130  ABI_ALLOCATE(magv_org,(nat_org))
131  ABI_ALLOCATE(xred_org,(3,nat_org))
132  ABI_ALLOCATE(chi0_org,(nat_org))
133  ABI_ALLOCATE(chi_org,(nat_org))
134  ABI_ALLOCATE(dparr,(marr,0:ndtuj))
135  ABI_ALLOCATE(dparrr,(marr,0:ndtuj))
136  ABI_ALLOCATE(dqarr,(marr,0:ndtuj))
137  ABI_ALLOCATE(dqarrr,(marr,0:ndtuj))
138  ABI_ALLOCATE(distv_org,(nat_org))
139  ABI_ALLOCATE(narrm,(0:ndtuj))
140  dparr=-one ;  dparrr=-one ;  dqarr=-one ;  dqarrr=-one
141 !DEBUG
142 !write(message,fmt='((a,i3,a))')'pawuj_det init sg'
143 !call wrtout(std_out,message,'COLL')
144 !END DEBUG
145  idum2=1 ; drarr=one
146  chih=zero
147 
148 !write(std_out,*) 'pawuj 02'
149 !###########################################################
150 !### 02. Create the file UJDET.nc
151 
152  if(.false.)write(std_out,*)ujdet_filename ! This is for the abirules
153 
154 !write(std_out,*) 'pawuj 03'
155 !###########################################################
156 !### 03. Write out the Input for UJDET
157 
158  write(message, '(3a)' ) ch10,&
159 & ' # input for ujdet, cut it using ''sed -n "/MARK/,/MARK/p" abi.out > ujdet.in ''------- ',ch10
160  call wrtout(ab_out,message,'COLL')
161 
162  if (ndtuj/=4) then
163    write (hstr,'(I0)') jdtset              ! convert integer to string
164  else
165    hstr=''
166  end if
167  if (ndtuj>=2.or.jdtset==1) then
168    idum2(1,1:ndtuj)=4 !dtpawuj(:)%ndtuj
169    call prttagm(dparr,idum2,ab_out,jdtset_,1,marr,1,narrm,ncid,ndtuj,'ndtset','INT',0)
170  end if
171  idum2(1,0:ndtuj)=pack(dtpawuj(:)%nat,dtpawuj(:)%iuj/=-1)
172  call prttagm(dparr,idum2,ab_out,jdtset_,1,marr,1,narrm,ncid,ndtuj,'nat'//trim(hstr),'INT',0)
173 
174  idum2(1,0:ndtuj)=pack(dtpawuj(:)%nspden,dtpawuj(:)%iuj/=-1)
175  call prttagm(dparr,idum2,ab_out,jdtset_,1,marr,1,narrm,ncid,ndtuj,'nspden'//trim(hstr),'INT',0)
176 
177  idum2(1,0:ndtuj)=pack(dtpawuj(:)%macro_uj,dtpawuj(:)%iuj/=-1)
178  call prttagm(dparr,idum2,ab_out,jdtset_,1,marr,1,narrm,ncid,ndtuj,'macro_uj'//trim(hstr),'INT',0)
179 
180  idum2(1,0:ndtuj)=pack(dtpawuj(:)%pawujat,dtpawuj(:)%iuj/=-1)
181  call prttagm(dparr,idum2,ab_out,jdtset_,1,marr,1,narrm,ncid,ndtuj,'pawujat'//trim(hstr),'INT',0)
182 
183  dparr(1,0:ndtuj)=pack(dtpawuj(:)%pawujga,dtpawuj(:)%iuj/=-1)
184  call prttagm(dparr,idum2,ab_out,jdtset_,1,marr,1,narrm,ncid,ndtuj,'pawujga'//trim(hstr),'DPR',0)
185 
186  dparr(1,0:ndtuj)=pack(dtpawuj(:)%pawujrad,dtpawuj(:)%iuj/=-1)
187  call prttagm(dparr,idum2,ab_out,jdtset_,1,marr,1,narrm,ncid,ndtuj,'pawujrad'//trim(hstr),'DPR',0)
188 
189  dparr(1,0:ndtuj)=pack(dtpawuj(:)%pawrad,dtpawuj(:)%iuj/=-1)
190  call prttagm(dparr,idum2,ab_out,jdtset_,1,marr,1,narrm,ncid,ndtuj,'pawrad'//trim(hstr),'DPR',0)
191 
192  dparr(1,0:ndtuj)=pack(dtpawuj(:)%ph0phiint,dtpawuj(:)%iuj/=-1)
193  call prttagm(dparr,idum2,ab_out,jdtset_,1,marr,1,narrm,ncid,ndtuj,'ph0phiint'//trim(hstr),'DPR',0)
194 
195  idum2(1,0:ndtuj)=pack(dtpawuj(:)%pawprtvol,dtpawuj(:)%iuj/=-1)
196  call prttagm(dparr,idum2,ab_out,jdtset_,1,marr,1,narrm,ncid,ndtuj,'pawprtvol'//trim(hstr),'INT',0)
197 
198  idum2(1,0:ndtuj)=pack(dtpawuj(:)%option,dtpawuj(:)%iuj/=-1)
199  call prttagm(dparr,idum2,ab_out,jdtset_,1,marr,1,narrm,ncid,ndtuj,'pawujopt'//trim(hstr),'INT',0)
200 
201  idum2(1,0:ndtuj)=pack(dtpawuj(:)%dmatpuopt,dtpawuj(:)%iuj/=-1)
202  call prttagm(dparr,idum2,ab_out,jdtset_,1,marr,1,narrm,ncid,ndtuj,'dmatpuopt'//trim(hstr),'INT',0)
203 
204  kdtset=0
205 
206  do idtset=0,ndtpawuj
207 !  DEBUG
208 !  write(message,fmt='((a,i3,a))')'pawuj_det m2, idtset ',idtset,ch10
209 !  call wrtout(std_out,message,'COLL')
210 !  call flush_unit(6)
211 !  END DEBUG
212    if (dtpawuj(idtset)%iuj/=-1) then
213      dparr(1:nspden*nat_org,kdtset)=reshape(dtpawuj(idtset)%vsh,(/nspden*nat_org/))
214      dparrr(1:nspden*nat_org,kdtset)=reshape(dtpawuj(idtset)%occ,(/nspden*nat_org/))
215 !    DEBUG
216 !    write(message,fmt='((a,i3,a))')'pawuj_det m3, idtset ',idtset,ch10
217 !    call wrtout(std_out,message,'COLL')
218 !    write(std_out,*)' marr,narrm,ncid,ndtuj,nat_org,kdtset,idtset=',marr,narrm,ncid,ndtuj,nat_org,kdtset,idtset
219 !    write(std_out,*)' dtpawuj(idtset)%xred=',dtpawuj(idtset)%xred
220 !    call flush_unit(6)
221 !    END DEBUG
222      dqarr(1:nat_org*3,kdtset)=reshape(dtpawuj(idtset)%xred,(/nat_org*3/))
223      dqarrr(1:3*3,kdtset)=reshape(dtpawuj(idtset)%rprimd,(/3*3/))
224      idum2(1:3,kdtset)=reshape(dtpawuj(idtset)%scdim,(/3/))
225      drarr(1:nwfchr,kdtset)=reshape(dtpawuj(idtset)%wfchr,(/nwfchr/))
226 !    DEBUG
227 !    write(message,fmt='((a,i3,a))')'pawuj_det m4, idtset ',idtset,ch10
228 !    call wrtout(std_out,message,'COLL')
229 !    call flush_unit(6)
230 !    END DEBUG
231      kdtset=kdtset+1
232    end if
233  end do
234 
235 !DEBUG
236 !write(message,fmt='((a,i3,a))')'pawuj_det m5'
237 !call wrtout(std_out,message,'COLL')
238 !END DEBUG
239  call prttagm(dparr,idum2,ab_out,jdtset_,2,marr,nspden*nat_org,narrm,ncid,ndtuj,'vsh'//trim(hstr),'DPR',0)
240  call prttagm(dparrr,idum2,ab_out,jdtset_,2,marr,nspden*nat_org,narrm,ncid,ndtuj,'occ'//trim(hstr),'DPR',0)
241 !DEBUG
242 !write(message,fmt='((a,i3,a))')'pawuj_det m6'
243 !call wrtout(std_out,message,'COLL')
244 !END DEBUG
245  call prttagm(dqarr,idum2,ab_out,jdtset_,2,marr,nat_org*3,narrm,ncid,ndtuj,'xred'//trim(hstr),'DPR',0)
246  call prttagm(dqarrr,idum2,ab_out,jdtset_,2,marr,3*3,narrm,ncid,ndtuj,'rprimd'//trim(hstr),'DPR',0)
247  call prttagm(dqarrr,idum2,ab_out,jdtset_,2,marr,3,narrm,ncid,ndtuj,'scdim'//trim(hstr),'INT',0)
248  call prttagm(drarr,idum2,ab_out,jdtset_,2,marr,nwfchr,narrm,ncid,ndtuj,'wfchr'//trim(hstr),'DPR',0)
249  ABI_DEALLOCATE(narrm)
250 
251  write(message, '( 15a )'  ) ch10,' # further possible options: ',ch10,&
252 & ' #    scdim    2 2 2 ',ch10,&
253 & ' #    mdist    10.0  ',ch10,&
254 & ' #  pawujga    2 '    ,ch10,&
255 & ' # pawujopt    2 '    ,ch10,&
256 & ' # pawujrad    3.0'   ,ch10,&
257 & ' # ------- end input for ujdet: end-MARK  -------- ',ch10
258  call wrtout(ab_out,message,'COLL')
259 
260 !write(std_out,*) 'pawuj 04'
261 !###########################################################
262 !### 04. Test
263 
264  if (ndtuj/=4)  return
265 
266  write(message, '(3a)' ) ch10,' ---------- calculate U, (J) start ---------- ',ch10
267  call wrtout(ab_out,message,'COLL')
268 
269  if (all(dtpawuj(1:ndtpawuj)%pawujat==pawujat)) then
270    write (message,fmt='(a,i3)') ' All pawujat  ok and equal to ',pawujat
271    call wrtout(ab_out,message,'COLL')
272  else
273    write (message,fmt='(a,4i3,2a)') ' Differing values of pawujat were found: ',dtpawuj(1:ndtuj)%pawujat,ch10,&
274 &   'No determination of U.'
275    call wrtout(ab_out,message,'COLL')
276    return
277  end if
278 
279  if (all(dtpawuj(1:ndtpawuj)%macro_uj==macro_uj)) then
280    if (nspden==1) then
281      write(message,fmt='(2a)') ' pawuj_det found nspden==1, determination',&
282 &     ' of U-parameter for unpol. struct. (non standard)'
283    else if (macro_uj==1.and.nspden==2) then
284      write(message,fmt='(2a)') ' pawuj_det: found macro_uj=1 and nspden=2:',&
285 &     ' standard determination of U-parameter'
286    else if (macro_uj==2.and.nspden==2) then
287      write(message,fmt='(2a)') ' pawuj_det: found macro_uj=2 and nspden=2:',&
288 &     ' determination of U on single spin channel (experimental)'
289 
290    else if (macro_uj==3.and.nspden==2) then
291      write(message,fmt='(2a)') ' pawuj_det: found macro_uj=3 and nspden=2,',&
292 &     ' determination of J-parameter on single spin channel (experimental)'
293    end if
294    write (message,fmt='(a,i3,a,a)') ' All macro_uj ok and equal to ',macro_uj,ch10,trim(message)
295    write (message,'(a,i3)') ' All macro_uj ok and equal to ',macro_uj
296    call wrtout(ab_out,message,'COLL')
297  else
298    write (message,fmt='(a,10i3)') ' Differing values of macro_uj were found: ',dtpawuj(:)%macro_uj
299    write (message,fmt='(3a)')trim(message),ch10,' No determination of U.'
300    call wrtout(ab_out,message,'COLL')
301    return
302  end if
303 
304  if (macro_uj>1.and.nspden==1) then
305    write (message,'(4a,2a)') ' U on a single spin channel (or J) can only be determined for nspden=2 ,',ch10,&
306 &   'No determination of U.'
307    call wrtout(ab_out,message,'COLL')
308    return
309  end if
310 
311 !Calculation of response matrix
312 
313  do jdtset=1,4
314    if (nspden==1) then
315      chih(jdtset,1:nat_org)=dtpawuj(jdtset)%occ(1,:)
316    else if (macro_uj==1.and.nspden==2) then
317      chih(jdtset,1:nat_org)=dtpawuj(jdtset)%occ(1,:)+dtpawuj(jdtset)%occ(2,:)
318    else if (macro_uj==2.and.nspden==2) then
319      chih(jdtset,1:nat_org)=dtpawuj(jdtset)%occ(1,:)
320    else if (macro_uj==3.and.nspden==2) then
321      chih(jdtset,1:nat_org)=dtpawuj(jdtset)%occ(2,:)
322    end if
323    vsh(jdtset)=dtpawuj(jdtset)%vsh(1,pawujat)
324    if (pawprtvol==3) then
325      write(message,fmt='(2a,i3,a,f15.12)') ch10,' Potential shift vsh(',jdtset,') =',vsh(jdtset)
326      call wrtout(std_out,message,'COLL')
327      write(message,fmt='( a,i3,a,120f15.9)') ' Occupations occ(',jdtset,') ',chih(jdtset,1:nat_org)
328      call wrtout(std_out,message,'COLL')
329    end if
330  end do
331 
332  if (any(abs((/(vsh(ii)-vsh(ii+2), ii=1,2) /))<0.00000001)) then
333    write(message, '(2a,18f10.7,a)' )  ch10,' vshift is too small: ',abs((/(vsh(ii)-vsh(ii+2), ii=1,2) /))
334    call wrtout(ab_out,message,'COLL')
335    return
336  end if
337 
338 !DEBUG
339 !write(message,fmt='(a)')'pawuj_det: after test vsh'
340 !call wrtout(std_out,message,'COLL')
341 !END DEBUG
342 
343  chi0_org=(chih(1,1:nat_org)-chih(3,1:nat_org))/(vsh(1)-vsh(3))/dtpawuj(1)%diemix
344  chi_org=(chih(2,1:nat_org)-chih(4,1:nat_org))/(vsh(2)-vsh(4))
345 
346  if (pawprtvol==3) then
347    write(message,fmt='(2a, 150f15.10)') ch10,' Chi_0n ',chi0_org
348    call wrtout(std_out,message,'COLL')
349    write(message,fmt='(a, 150f15.10)') ' Chi_n ',chi_org
350    call wrtout(std_out,message,'COLL')
351  end if
352 
353  write(message,fmt='(a)')': '
354  if (nspden==2) then
355    magv_org=dtpawuj(1)%occ(1,:)-dtpawuj(1)%occ(2,:)
356    if (all(abs(magv_org)<0.001)) then
357      magv_org=(/(1,im1=1,nat_org)/)
358    else
359      magv_org=abs(magv_org)/magv_org
360      if (magv_org(1).lt.0) magv_org=magv_org*(-1_dp)
361      if (all(magv_org(2:nat_org).lt.0)) then
362        magv_org=abs(magv_org)
363        write(message,'(a)')', (reset to fm): '
364      end if
365    end if
366  else
367    magv_org=(/(1,im1=1,nat_org)/)
368  end if
369 
370  if (pawprtvol==3) then
371    write(message,fmt='(3a, 150f4.1)') ch10,' Magnetization',trim(message),magv_org
372    call wrtout(std_out,message,'COLL')
373  end if
374 
375 !Case of extrapolation to larger r_paw: calculate intg
376 
377  if (all(dtpawuj(1)%wfchr(:)/=0).and.ph0phiint/=1) then
378    if (dtpawuj(1)%pawujrad<20.and.dtpawuj(1)%pawujrad>dtpawuj(1)%pawrad) then
379      fcorr=(1-ph0phiint)/(IRadFnH(dtpawuj(1)%pawrad,20.0_dp,nint(dtpawuj(1)%wfchr(2)),&
380 &     nint(dtpawuj(1)%wfchr(3)),dtpawuj(1)%wfchr(1)))
381      intg=ph0phiint/(1-fcorr*IRadFnH(dtpawuj(1)%pawujrad,20.0_dp,nint(dtpawuj(1)%wfchr(2)),&
382 &     nint(dtpawuj(1)%wfchr(3)),dtpawuj(1)%wfchr(1)))
383      write(message, fmt='(a,f12.5,a,f12.5)') ' pawuj_det: met2 extrapolation to ', dtpawuj(1)%pawujrad,' using factor ',intg
384      call wrtout(std_out,message,'COLL')
385    else if (dtpawuj(1)%pawujrad<dtpawuj(1)%pawrad) then
386      a=0 ; a(1:3)=dtpawuj(1)%wfchr(4:6)
387      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/)
388      b=(/(dtpawuj(1)%pawujrad**(im1)-dtpawuj(1)%pawrad**(im1),im1=1,5)/)
389      intg=dot_product(a,b)
390      intg=ph0phiint/(ph0phiint+intg)
391      write(message, fmt='(a,f12.5,a,f12.5)') ' pawuj_det: met1 extrapolation to ', dtpawuj(1)%pawujrad,' using factor ',intg
392      call wrtout(std_out,message,'COLL')
393    else if (dtpawuj(1)%pawujrad==dtpawuj(1)%pawrad) then
394      intg=one
395      write(message, fmt='(a,2i7,3f12.5)') ' pawuj_det: no extrapolation (pawujrad=pawrad)'
396      call wrtout(std_out,message,'COLL')
397    else
398      intg=ph0phiint
399      write(message, fmt='(a,3f12.5)') ' pawuj_det: extrapolation to r->\inf using factor ',intg
400      call wrtout(std_out,message,'COLL')
401    end if
402  else
403    write(message, fmt='(a,2i7,3f12.5)') ' pawuj_det: no extrapolation (ph0phiint=1 or wfchr=0)'
404    call wrtout(std_out,message,'COLL')
405    intg=one
406  end if
407 
408 
409 !Determine U in primitive cell
410 
411  write(message,fmt='(a)')' pawuj_det: determine U in primitive cell'
412  call wrtout(std_out,message,'COLL')
413 
414  call lcalcu(int(magv_org),nat_org,dtpawuj(1)%rprimd,dtpawuj(1)%xred,chi_org,chi0_org,pawujat,ures,pawprtvol,pawujga,pawujoption)
415 
416 !Begin calculate U in supercell
417 
418 !Analize shell structure of primitive cell
419 !and atomic distances in primitive cell
420  ABI_ALLOCATE(smult_org,(nat_org))
421  ABI_ALLOCATE(sdistv_org,(nat_org))
422  call shellstruct(dtpawuj(1)%xred,dtpawuj(1)%rprimd,nat_org,&
423 & int(magv_org),distv_org,smult_org,sdistv_org,nsh_org,pawujat,pawprtvol)
424 
425  ii=1
426  write(message, fmt='(8a)') ' URES ','     ii','    nat','       r_max','    U(J)[eV]','   U_ASA[eV]','   U_inf[eV]',ch10
427  write(message, fmt='(a,2i7,4f12.5)') trim(message)//' URES ',ii,nat_org,maxval(abs(distv_org)),ures,ures*exp(log(intg)*eyp),&
428 & ures*exp(log(ph0phiint)*eyp)
429  call wrtout(std_out,message,'COLL')
430  call wrtout(ab_out,message,'COLL')
431 
432  if (pawprtvol>1) then
433    write(message,fmt='(a, 150f10.5)')' pawuj_det: ionic distances in original cell (distv_org) ', distv_org
434    call wrtout(std_out,message,'COLL')
435  end if
436 
437 !Construct supercell, calculate limit dimensions of supercell
438  ii=1
439  maxnat=product(dtpawuj(1)%scdim(:))*nat_org
440  if (maxnat==0) then
441    maxnat=dtpawuj(1)%scdim(1)
442    ext=(/ii, ii, ii/)
443  else
444    jj=1
445    do while (jj<=3)
446      ext(jj)=minval( (/ii, dtpawuj(1)%scdim(jj) /) )
447      jj=jj+1
448    end do
449  end if
450  ext=ext+(/ 1,1,1 /)
451  ii=ii+1
452 
453 
454  nat_sc=product(ext)*nat_org
455 
456 !DEBUG
457 !write(message,fmt='(a,3i4)')'pawuj_det: ext ',ext
458 !call wrtout(std_out,message,'COLL')
459 !END DEBUG
460 
461  do while (nat_sc<=maxnat)
462    ABI_ALLOCATE(chi0_sc,(nat_sc))
463    ABI_ALLOCATE(chi_sc,(nat_sc))
464    ABI_ALLOCATE(distv_sc,(nat_sc))
465    ABI_ALLOCATE(magv_sc,(nat_sc))
466    ABI_ALLOCATE(sdistv_sc,(nat_sc))
467    ABI_ALLOCATE(smult_sc,(nat_sc))
468    ABI_ALLOCATE(xred_sc,(3,nat_sc))
469 
470 !  Determine positions=xred_sc and supercell dimensions=rpimd_sc
471    call mksupercell(dtpawuj(1)%xred,int(magv_org),dtpawuj(1)%rprimd,nat_org,&
472 &   nat_sc,xred_sc,magv_sc,rprimd_sc,ext,pawprtvol)
473 
474 !  Determine shell structure of supercell: multiplicities (smult_sc), radii of shells (sdistv_sc)
475 !  number of shells (nsh_sc) and atomic distances in supercell (distv_sc)
476 
477    write(message,fmt='(a,3i3,a)')' pawuj_det: determine shell structure of ',ext(1:3),' supercell'
478    call wrtout(std_out,message,'COLL')
479 
480    call shellstruct(xred_sc,rprimd_sc,nat_sc,int(magv_sc),distv_sc,smult_sc,sdistv_sc,nsh_sc,pawujat,pawprtvol)
481 
482    if (pawprtvol>=2) then
483      write(message,fmt='(a)')' pawuj_det: ionic distances in supercell (distv_sc) '
484      call wrtout(std_out,message,'COLL')
485      call prmat(distv_sc,1,nat_sc,1,std_out)
486    end if
487 
488 !  Determine chi and chi0 in supercell (chi0_sc, chi_sc)
489 !  DEBUG
490 !  write(message,fmt='(a)')'pawuj_det:  chi and chi0 in supercell'
491 !  call wrtout(std_out,message,'COLL')
492 !  END DEBUG
493 
494    if (pawujoption>2) then
495      invopt=2
496    else
497      invopt=1
498    end if
499 
500    call chiscwrt(chi0_org,distv_org,nat_org,sdistv_org,smult_org,nsh_org,&
501 &   chi0_sc,distv_sc,nat_sc,smult_sc,nsh_sc,invopt,pawprtvol)
502    call chiscwrt(chi_org,distv_org,nat_org,sdistv_org,smult_org,nsh_org,&
503 &   chi_sc,distv_sc,nat_sc,smult_sc,nsh_sc,invopt,pawprtvol)
504 
505 !  Calculate U in supercell
506 !  DEBUG
507 !  write(message,fmt='(a)')'pawuj_det:   U in supercell'
508 !  call wrtout(std_out,message,'COLL')
509 !  END DEBUG
510    call lcalcu(int(magv_sc),nat_sc,rprimd_sc,xred_sc,chi_sc,chi0_sc,pawujat,ures,pawprtvol,pawujga,pawujoption)
511 
512    write(message, fmt='(a,2i7,4f12.5)') ' URES ',ii,nat_sc,maxval(abs(distv_sc)),ures,ures*exp(log(intg)*eyp),&
513 &   ures*exp(log(ph0phiint)*eyp)
514    call wrtout(std_out,message,'COLL')
515    call wrtout(ab_out,message,'COLL')
516 
517    ABI_DEALLOCATE(chi0_sc)
518    ABI_DEALLOCATE(chi_sc)
519    ABI_DEALLOCATE(distv_sc)
520    ABI_DEALLOCATE(magv_sc)
521    ABI_DEALLOCATE(sdistv_sc)
522    ABI_DEALLOCATE(smult_sc)
523    ABI_DEALLOCATE(xred_sc)
524 
525    if (product(dtpawuj(1)%scdim(:))==0) then
526      ext=(/ii, ii, ii/)
527    else
528      jj=1
529      do while (jj<=3)
530        ext(jj)=minval( (/ii, dtpawuj(1)%scdim(jj) /) )
531        jj=jj+1
532      end do
533    end if
534    ext=ext+(/ 1,1,1 /)
535    ii=ii+1
536 
537    nat_sc=product(ext)*nat_org
538 
539  end do
540 
541  write(message, '(3a)' )ch10,' ---------- calculate U, (J) end -------------- ',ch10
542  call wrtout(ab_out,message,'COLL')
543 
544  ABI_DEALLOCATE(dparrr)
545  ABI_DEALLOCATE(dqarr)
546  ABI_DEALLOCATE(dqarrr)
547  ABI_DEALLOCATE(jdtset_)
548  ABI_DEALLOCATE(chi_org)
549  ABI_DEALLOCATE(chi0_org)
550  ABI_DEALLOCATE(smult_org)
551  ABI_DEALLOCATE(sdistv_org)
552  ABI_DEALLOCATE(chih)
553  ABI_DEALLOCATE(idum2)
554  ABI_DEALLOCATE(drarr)
555  ABI_DEALLOCATE(dparr)
556  ABI_DEALLOCATE(magv_org)
557  ABI_DEALLOCATE(xred_org)
558  ABI_DEALLOCATE(distv_org)
559 
560  DBG_EXIT("COLL")
561 
562 end subroutine pawuj_det