TABLE OF CONTENTS
ABINIT/chiscwrt [ 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 ]
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