Branch data Line data Source code
1 : : !!****p* ABINIT/lruj
2 : : !! NAME
3 : : !! lruj
4 : : !!
5 : : !! FUNCTION
6 : : !! Linear Response U and J:
7 : : !! Determines Hubbard U or Hund's J from series of *DS*_LRUJ.nc
8 : : !! files containing information regarding the perturbation applied to a particular
9 : : !! atom and the resulting occupations/magnetizations. The procedure implemented
10 : : !! is that of the SCF linear response for the Hubbard U (Phys. Rev. B 71,035105)
11 : : !! and the Hund's J (Phys. Rev. B 98, 235157) parameters.
12 : : !! This protocol was coded up in November 2022 by Lorien MacEnulty (macenulty.com),
13 : : !! doctoral researcher in the Quantum Theory of Materials group (theoryofmaterials.com)
14 : : !! at Trinity College Dublin, headed by Dr. David O'Regan.
15 : : !!
16 : : !! COPYRIGHT
17 : : !! Copyright (C) 1998-2026 ABINIT group (LMac)
18 : : !! This file is distributed under the terms of the
19 : : !! GNU General Public License, see ~abinit/COPYING
20 : : !! or http://www.gnu.org/copyleft/gpl.txt .
21 : : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
22 : : !!
23 : : !! INPUTS
24 : : !! Executed as ./lruj *LRUJ.nc FILE1 FILE2 FILE3 ... [--d 5] [--help] [--version]
25 : : !! --d <n> = Command line argument: highest degree of intended polynomial fits
26 : : !! *DS*_LRUJ.nc files = gives data from perturbative Abinit calculations
27 : : !!
28 : : !! OUTPUT
29 : : !! std_out = log file
30 : : !!
31 : : !! SOURCE
32 : :
33 : : #if defined HAVE_CONFIG_H
34 : : #include "config.h"
35 : : #endif
36 : :
37 : : #include "abi_common.h"
38 : :
39 : :
40 : 4 : program lruj
41 : :
42 : 4 : use defs_basis
43 : : use m_xmpi
44 : : use m_abicore
45 : : use m_build_info
46 : : use m_errors
47 : : use m_argparse
48 : : use m_crystal
49 : : use netcdf
50 : : use m_nctk
51 : : use m_yaml
52 : :
53 : : use m_fstrings, only : itoa, sjoin, ltoa
54 : : use m_specialmsg, only : specialmsg_getcount, herald
55 : : use m_numeric_tools, only : polynomial_regression
56 : : use m_sort, only : sort_dp
57 : : use m_common, only : crystal_from_file
58 : : use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq
59 : : use m_paw_uj, only : pawuj_ini,pawuj_free,pawuj_det, macro_uj_type
60 : :
61 : : implicit none
62 : :
63 : : !Local variables-------------------------------
64 : :
65 : : !scalars
66 : : integer,parameter :: master=0
67 : : integer :: nproc,my_rank,comm
68 : : logical :: iam_master
69 : :
70 : : integer :: ncid,nnat,natom,prtvol,nargs,nfiles,ndtpawuj,degarg
71 : : integer :: ndata,nspden,macro_uj,pawujat,dmatpuopt
72 : : integer :: degree,mdegree,ii,ipert
73 : : real(dp) :: diem,ph0phiint,signum !diemix, diemixmag,
74 : : type(yamldoc_t) :: ydoc
75 : : !type(crystal_t) :: cryst
76 : :
77 : : !arrays
78 : 4 : integer, allocatable :: iperm(:),pawujat_file(:),macrouj_file(:),dmatpuopt_file(:)
79 : 4 : real(dp), allocatable :: diem_file(:),ph0phiint_file(:),nspden_file(:)
80 : :
81 : 4 : real(dp), allocatable :: perts(:),occs0(:),occs(:)
82 : 4 : real(dp), allocatable :: uj_perts(:),luocc(:,:),luocc_nnat(:,:)
83 : 4 : real(dp), allocatable :: chi0coeffs(:),chicoeffs(:),chi0(:),chi(:),hubpar(:)
84 : 4 : real(dp), allocatable :: chi0err(:),chierr(:),hubparerr(:)
85 : :
86 : : !characters
87 : : character(len=1) :: parname
88 : : character(len=5) :: degreename
89 : : character(len=12) :: regname
90 : : character(len=14) :: occmag
91 : : character(len=24) :: codename
92 : : character(len=30) :: diem_token
93 : : character(len=500) :: message,arg,msg,pertname
94 : 4 : character(len=fnlen),allocatable :: file_paths(:)
95 : :
96 : :
97 : : !##########################################################################################################
98 : : !################################## Set up MPI architecture (unused) ####################################
99 : :
100 : : !Change communicator for I/O (mandatory!)
101 : 4 : call abi_io_redirect(new_io_comm=xmpi_world)
102 : :
103 : : !Initialize MPI (not used but necessary)
104 : 4 : call xmpi_init()
105 : 4 : comm = xmpi_world
106 : 4 : nproc = xmpi_comm_size(comm)
107 : 4 : my_rank = xmpi_comm_rank(comm)
108 : 4 : iam_master = (my_rank == master)
109 : :
110 : : !Initialize memory profiling if it is activated
111 : : !if a full abimem.mocc report is desired, set the argument of abimem_init to "2" instead of "0"
112 : : !note that abimem.mocc files can easily be multiple GB in size so don't use this option normally
113 : : #ifdef HAVE_MEM_PROFILING
114 : : call abimem_init(0)
115 : : #endif
116 : :
117 : : !No MPI functionality needed for main procedure.
118 [ - + ]: 4 : if (my_rank /= master) goto 100
119 : :
120 : : !##########################################################################################################
121 : : !###################################### Read command line options #######################################
122 : :
123 : : !Count arguments and number of files (= #perturbations)
124 : 4 : nargs = command_argument_count()
125 [ + - - + : 4 : ABI_MALLOC(file_paths, (nargs))
- + - + -
+ - + -
+ ]
126 : 4 : nfiles = 0
127 [ + + ]: 22 : do ii=1,nargs
128 : 18 : call get_command_argument(ii, arg)
129 [ - + ]: 18 : if (arg(1:1) == "-") exit
130 : 18 : nfiles = nfiles + 1
131 [ + - + - ]: 22 : file_paths(nfiles) = trim(arg)
132 : : end do
133 : :
134 : : !Assess options
135 [ + + ]: 22 : do ii=1,command_argument_count()
136 : 18 : call get_command_argument(ii, arg)
137 [ - + ]: 22 : if (arg == "--version") then
138 : 0 : write(std_out,"(a)") trim(abinit_version); goto 100
139 [ - + ]: 18 : else if (arg == "-h" .or. arg == "--help") then
140 : : !Document the options.
141 : 0 : call lruj_show_help()
142 : 0 : goto 100
143 : : end if
144 : : end do
145 : :
146 : : !If no files found, exit program.
147 [ - + ]: 4 : if (nfiles == 0) then
148 : 0 : write(std_out, *) "Empty file list!"
149 : 0 : goto 100
150 : : end if
151 : :
152 : : !Get other options from the CLI. e.g. --prtvol 0 -d 3.0
153 : : !Should be documented in lruj_show_help
154 [ - + ]: 4 : ABI_CHECK(get_arg("prtvol", prtvol, msg, default=0) == 0, msg)
155 [ - + ]: 4 : ABI_CHECK(get_arg("d", degarg, msg, default=1) == 0, msg)
156 : :
157 : : !Print header
158 : 4 : codename='LRUJ'//repeat(' ',18)
159 : 4 : call herald(codename, abinit_version, std_out)
160 : :
161 : : !##########################################################################################################
162 : : !###################################### Read *DSi*_LRUJ.nc files ########################################
163 : :
164 [ + - - + : 4 : ABI_MALLOC(uj_perts,(nfiles))
- + - + -
+ - + -
+ ]
165 [ + - - + : 4 : ABI_MALLOC(macrouj_file, (nfiles))
- + - + -
+ - + -
+ ]
166 : :
167 : : !Read perturbation strengths, parameter type, #spins and #atoms
168 : : !from each file.
169 [ + + ]: 22 : do ii=1,nfiles
170 [ - + ]: 18 : NCF_CHECK(nctk_open_read(ncid, file_paths(ii), xmpi_comm_self))
171 [ - + ]: 18 : NCF_CHECK(nf90_get_var(ncid, vid("uj_pert"), uj_perts(ii)))
172 [ - + ]: 18 : NCF_CHECK(nf90_get_var(ncid, vid("macro_uj"), macrouj_file(ii)))
173 : 18 : macro_uj=macrouj_file(1)
174 : : !Make sure ndtpawuj is always 4.
175 [ - + ]: 18 : NCF_CHECK(nctk_get_dim(ncid, "ndtpawuj", ndtpawuj))
176 [ - + ]: 18 : ABI_CHECK_IEQ(ndtpawuj, 4, "Wrong ndtpawuj")
177 [ - + ]: 18 : NCF_CHECK(nctk_get_dim(ncid, "nspden", nspden))
178 [ - + ]: 18 : NCF_CHECK(nctk_get_dim(ncid, "nnat", nnat))
179 [ - + ]: 18 : NCF_CHECK(nctk_get_dim(ncid, "natom", natom))
180 [ - + ]: 22 : NCF_CHECK(nf90_close(ncid))
181 : : end do
182 : :
183 : : !Sort files by perturbation magnitude.
184 [ + - - + : 4 : ABI_MALLOC(iperm, (nfiles))
- + - + -
+ - + -
+ ]
185 [ - + + + : 58 : iperm = [(ii, ii=1,nfiles)]
- + + + +
- + - - -
+ + ]
186 : 4 : call sort_dp(nfiles, uj_perts, iperm, tol12)
187 [ - + + + : 40 : file_paths(1:nfiles) = file_paths(iperm(:))
+ + ]
188 [ - + ]: 4 : ABI_FREE(iperm)
189 : :
190 : : !Allocate main data-holding arrays.
191 [ + - - + : 4 : ABI_MALLOC(luocc, (ndtpawuj, nfiles))
+ - - + -
+ - + - +
- + - + ]
192 [ + - - + : 4 : ABI_MALLOC(luocc_nnat, (ndtpawuj, nnat))
+ - - + -
+ - + - +
- + - + ]
193 [ + - - + : 4 : ABI_MALLOC(pawujat_file, (nfiles))
- + - + -
+ - + -
+ ]
194 [ + - - + : 4 : ABI_MALLOC(diem_file, (nfiles))
- + - + -
+ - + -
+ ]
195 [ + - - + : 4 : ABI_MALLOC(dmatpuopt_file, (nfiles))
- + - + -
+ - + -
+ ]
196 [ + - - + : 4 : ABI_MALLOC(ph0phiint_file, (nfiles))
- + - + -
+ - + -
+ ]
197 [ + - - + : 4 : ABI_MALLOC(nspden_file, (nfiles))
- + - + -
+ - + -
+ ]
198 : :
199 : : !Set macro_uj-specific variables, strings and constants.
200 [ + + ]: 4 : if (macro_uj==4) then !Calculation of the Hunds J parameter
201 : 2 : diem_token="diemixmag" !Unscreened response in Hund's J impacted by diemixmag
202 : 2 : pertname='beta ' !Hund's J perturbation: +beta to spin up, -beta to down
203 : 2 : parname='J'
204 : 2 : occmag='Magnetizations' !Magnetic moments are monitored.
205 : 2 : signum=-1.0d0 !Hund's J is -1*(1/chi0-1/chi)
206 : : else
207 : 2 : diem_token="diemix" !Unscreened response in Hubbard U impacted by diemix
208 : 2 : pertname='alpha' !Hubbard U perturbation; applied equally to spins up and down
209 : 2 : parname='U'
210 : 2 : occmag=' Occupations' !Total occupation is monitored.
211 : 2 : signum=1.0d0 !Hubbard U is 1*(1/chi0-1/chi)
212 : : end if
213 : :
214 : : !Allocate perturbation and occupation arrays. Set the first
215 : : !to the unperturbed case (i.e., when perturbation=0.0d0).
216 [ + - - + : 4 : ABI_MALLOC(perts,(0:nfiles))
- + - + -
+ - + -
+ ]
217 [ + - - + : 4 : ABI_MALLOC(occs0,(0:nfiles))
- + - + -
+ - + -
+ ]
218 [ + - - + : 4 : ABI_MALLOC(occs,(0:nfiles))
- + - + -
+ - + -
+ ]
219 : 4 : perts(0)=0.0d0 !Unperturbed case.
220 : :
221 : 4 : write(std_out,'(a,i2)') ' Number of perturbations detected: ',nfiles
222 : :
223 : : !Open _LRUJ.nc files and read in the relevant data.
224 [ + + ]: 22 : do ii=1,nfiles
225 [ - + ]: 18 : NCF_CHECK(nctk_open_read(ncid, file_paths(ii), xmpi_comm_self))
226 [ - + ]: 18 : NCF_CHECK(nf90_get_var(ncid, vid("pawujat"), pawujat_file(ii)))
227 : 18 : pawujat=pawujat_file(1)
228 [ - + ]: 18 : NCF_CHECK(nf90_get_var(ncid, vid("nspden"), nspden_file(ii)))
229 : 18 : nspden=nspden_file(1)
230 [ - + ]: 18 : NCF_CHECK(nf90_get_var(ncid, vid("luocc"), luocc_nnat))
231 [ + + ]: 90 : luocc(:,ii) = luocc_nnat(:,pawujat)
232 [ - + ]: 18 : NCF_CHECK(nf90_get_var(ncid, vid(diem_token), diem_file(ii)))
233 : 18 : diem=diem_file(1)
234 [ - + ]: 18 : NCF_CHECK(nf90_get_var(ncid, vid("dmatpuopt"), dmatpuopt_file(ii)))
235 : 18 : dmatpuopt=dmatpuopt_file(1)
236 [ - + ]: 18 : NCF_CHECK(nf90_get_var(ncid, vid("ph0phiint"), ph0phiint_file(ii)))
237 : 18 : ph0phiint=ph0phiint_file(1)
238 [ - + ]: 18 : NCF_CHECK(nf90_close(ncid))
239 : : !Testing if the unperturbed occupancies are equal across each run. If they
240 : : !aren't, then exit. If they are equal, then save them in appropriate arrays.
241 [ - + ]: 22 : if ((ii>1).and.((occs0(0)/=luocc(1,ii)).or.(occs(0)/=luocc(2,ii)))) then
242 : 0 : write(std_out,'(a)') "ERROR: Unperturbed ground state occupations across LRUJ datasets are not equal."
243 : 0 : write(std_out,'(a)') "Check the consistency of input variables in your perturbative calculations:"
244 : 0 : write(std_out,'(a)') " 1. Are they each reading in the same WFK file?"
245 : 0 : write(std_out,'(a)') " 2. Are macro-uj, pawujat, dmatpuopt, diemix(mag) consistent"
246 : 0 : write(std_out,'(a)') " across all perturbations?"
247 : 0 : write(std_out,'(a)') "If not, relaunch perturbative Abinit calculations, then"
248 : 0 : write(std_out,'(2a)') "reexecute lruj utility. Exiting.",ch10
249 : 0 : goto 100
250 : : else
251 : 18 : perts(ii)=uj_perts(ii)*Ha_eV
252 : 18 : occs0(0)=luocc(1,ii)
253 : 18 : occs(0)=luocc(2,ii)
254 : 18 : occs0(ii)=luocc(3,ii)
255 : 18 : occs(ii)=luocc(4,ii)
256 : : end if
257 : : end do
258 : :
259 : : !##########################################################################################################
260 : : !#################################### Tests on input information ########################################
261 : :
262 : : !Tests if we have enough data points (at least 3) to conduct distinct regression.
263 : 4 : ndata=nfiles+1
264 : 4 : write(std_out,'(a,i2,a)') ' Including unperturbed state, we have ',ndata,' data points.'
265 [ - + ]: 4 : if (ndata==0) then
266 : 0 : ABI_ERROR('No linear response data points found.')
267 [ - + ]: 4 : else if (ndata==1) then
268 : : msg = sjoin('Only one data point found. This utility needs',ch10,&
269 : : 'at least three (3) data points (two non-zero perturbations and one unperturbed) to compute',ch10,&
270 [ # # ]: 0 : 'the Hubbard parameter.')
271 : 0 : ABI_ERROR(msg)
272 [ - + ]: 4 : else if (ndata==2) then
273 : : msg = sjoin('Only two data points found. The scalar Hubbard Parameter from',ch10,&
274 : : 'the two-point linear regression scheme has already been printed in your .abo file. Try bashing',ch10,&
275 [ # # ]: 0 : '==> grep "two-point regression" <run_name.abo> ',ch10,'to find the result of this calculation.')
276 : 0 : ABI_ERROR(msg)
277 : : end if
278 : :
279 : : !pawujat consistency check.
280 [ + + - + : 22 : if (any(pawujat_file /= pawujat_file(1))) then
- + ]
281 : : msg = sjoin("Found different values of pawujat in files: ", ltoa(pawujat_file),ch10,&
282 [ # # ]: 0 : "Perturbed atom has to be consistent across perturbations to compute U or J.")
283 : 0 : ABI_ERROR(msg)
284 : : end if
285 : :
286 : : !dmatpuopt consistency check
287 [ + + - + : 22 : if (any(dmatpuopt_file /= dmatpuopt_file(1))) then
- + ]
288 : : msg = sjoin("Found different values of dmatpuopt in files: ", ltoa(dmatpuopt_file),ch10,&
289 [ # # ]: 0 : "PAW projector must be consistent across perturbations to compute U or J.")
290 : 0 : ABI_ERROR(msg)
291 : : end if
292 : :
293 : : !macro_uj consistency check
294 [ + + - + : 22 : if (any(macrouj_file /= macrouj_file(1))) then
- + ]
295 : : msg = sjoin("Found different values of macro_uj in files: ",ltoa(macrouj_file),ch10,&
296 [ # # ]: 0 : "Perturbation protocol and occupancy monitoring must be consistent to compute U or J.")
297 : 0 : ABI_ERROR(msg)
298 : : end if
299 : :
300 : : !diemix/diemixmag consistency check
301 [ + + - + : 22 : if (any(diem_file /= diem_file(1))) then
- + ]
302 : : msg = sjoin("Found different values of mixing constant in files: ",ltoa(diem_file),ch10,&
303 [ # # ]: 0 : "Unscreened response functions will factor into U (J) incorrectly.")
304 : 0 : ABI_ERROR(msg)
305 : : end if
306 : :
307 : : !Tests consistency of macro_uj, then writes message about macro_uj procedure selected.
308 : : !Also assigns Hubbard parameter-specific variables.
309 [ - + ]: 4 : if (nspden==1) then
310 : 0 : write(message,'(a)') ' Determination of U-parameter for unpolarized structure (non standard)'
311 [ + + ]: 4 : else if (macro_uj==1.and.nspden==2) then
312 : 2 : write(message,'(a)') ' Standard determination of the Hubbard U parameter.'
313 [ - + ]: 2 : else if (macro_uj==2.and.nspden==2) then
314 : 0 : write(message,'(a)') ' Determination of parameter on single spin channel (experimental)'
315 : 0 : pertname='Pert. '
316 [ - + ]: 2 : else if (macro_uj==3.and.nspden==2) then
317 : 0 : parname='J'
318 : 0 : pertname='Pert. '
319 : 0 : write(message,'(a)') ' Determination of (not Hunds) J-parameter on single spin channel (experimental)'
320 [ + - ]: 2 : else if (macro_uj==4.and.nspden==2) then
321 : 2 : write(message,'(a)') ' Hunds J determination, implemented by L. MacEnulty August 2021'
322 : : end if
323 : 4 : call wrtout(std_out,message)
324 : :
325 : : !Tests compatibility of nspden and macro_uj
326 [ - + ]: 4 : if (macro_uj>1.and.nspden==1) then
327 : : msg = sjoin('U on a single spin channel (or J) can only be determined for nspden=2 ,',ch10,&
328 [ # # ]: 0 : 'Cannot calculate the chosen Hubbard parameter.')
329 : 0 : ABI_ERROR(msg)
330 : : end if
331 : :
332 : : !Tests if perturbations are too small.
333 [ + - - + : 26 : if (maxval(abs(uj_perts))<0.00000001) then
- - + + +
+ - + ]
334 : : msg = sjoin('Perturbation magnitudes are too small.',ch10,&
335 [ # # ]: 0 : 'Rerun perturbative Abinit calculations with pawujv >> 1d-8.')
336 : 0 : ABI_ERROR(msg)
337 : : end if
338 : :
339 : : !##########################################################################################################
340 : : !############################### Calculation of the Response Functions ##################################
341 : :
342 : : !Test compatibility of polynomial degree (if present as an argument) with
343 : : !number of data points. Otherwise, default to the following:
344 : : !If we have 3 data points, conduct at maximum a linear regression.
345 : : !If 4 data points, conduct linear and quadratic regressions.
346 : : !If more than 5 data points, conduct linear, quadratic and cubic regressions.
347 [ - + ]: 4 : if (degarg/=1) then
348 [ # # ]: 0 : if (degarg>ndata-2) then
349 : 0 : write(std_out,'(4a,i2,3a,i2,2a)') ch10,' ERROR: Your chosen polynomial degree is too large. The resulting',ch10,&
350 : 0 : & ' parameters will certainly be overfitted. Either conduct ',degarg+1,' perturbations,',ch10,&
351 : 0 : & ' or execute this utility again with --d',ndata-2,' or smaller. Exiting program.',ch10
352 : 0 : goto 100
353 : : else
354 : 0 : mdegree=degarg
355 : : end if
356 : : else !Default max polynomial degree
357 [ - + ]: 4 : if (ndata<=4) then
358 : 0 : mdegree=ndata-2
359 : : else
360 : 4 : mdegree=3
361 : : end if
362 : : end if
363 : 4 : write(std_out,'(a,i2)') ' Maximum degree of polynomials analyzed: ',mdegree
364 : :
365 : : !Write warning about response matrices
366 : 4 : write(std_out,'(5a)') ' NOTE: Unlike the ujdet utility, lruj treats the ',ch10,&
367 : 4 : ' response functions as scalars, not matrices!',ch10,&
368 : 8 : ' See lruj tutorial for more information.'
369 : :
370 : : !Allocate the response and error arrays
371 [ + - - + : 4 : ABI_MALLOC(chi0err,(mdegree))
- + - + -
+ - + -
+ ]
372 [ + - - + : 4 : ABI_MALLOC(chierr,(mdegree))
- + - + -
+ - + -
+ ]
373 [ + - - + : 4 : ABI_MALLOC(chi0,(mdegree))
- + - + -
+ - + -
+ ]
374 [ + - - + : 4 : ABI_MALLOC(chi,(mdegree))
- + - + -
+ - + -
+ ]
375 [ + - - + : 4 : ABI_MALLOC(hubpar,(mdegree))
- + - + -
+ - + -
+ ]
376 [ + - - + : 4 : ABI_MALLOC(hubparerr,(mdegree))
- + - + -
+ - + -
+ ]
377 : :
378 : : !Start to write information in YAML format to plot with AbiPY
379 : 4 : ydoc = yamldoc_open('LRUJ_Abipy_Plots') !, width=11, real_fmt='(3f8.3)')
380 : 4 : call ydoc%add_int("natom",natom)
381 : 4 : call ydoc%add_int("ndata",ndata)
382 : 4 : call ydoc%add_int("pawujat",pawujat)
383 : 4 : call ydoc%add_int("macro_uj",macro_uj)
384 : 4 : call ydoc%add_string("diem_token",diem_token)
385 : 4 : call ydoc%add_real("diem",diem)
386 : :
387 : : !For all regressions, call subroutine to calculate polynomial fit for chi0 and chi.
388 [ + + ]: 16 : do degree=1,mdegree
389 [ + - + - : 12 : ABI_MALLOC(chi0coeffs,(degree+1))
- + - + -
+ - + - +
- + ]
390 [ + - + - : 12 : ABI_MALLOC(chicoeffs,(degree+1))
- + - + -
+ - + - +
- + ]
391 : 12 : call polynomial_regression(degree,ndata,perts,occs0,chi0coeffs,chi0err(degree))
392 : 12 : call polynomial_regression(degree,ndata,perts,occs,chicoeffs,chierr(degree))
393 : :
394 : : !YAML doc information on regression coefficients
395 : 12 : write(message, '(a,i0)' ) 'chi0_coefficients_degree',degree
396 : 12 : call ydoc%add_real1d(message,chi0coeffs)
397 : 12 : write(message, '(a,i0)' ) 'chi_coefficients_degree',degree
398 : 12 : call ydoc%add_real1d(message,chicoeffs)
399 : :
400 : 12 : chi0(degree)=chi0coeffs(2)/diem !The derivative of all polynomial regressions
401 : 12 : chi(degree)=chicoeffs(2) !at pert=0.0 is just the second coefficient.
402 : 12 : chi0err(degree)=chi0err(degree)/diem !Chi0 error divided by diem also.
403 : 12 : hubpar(degree)=signum*(1.0d0/chi0(degree)-1.0d0/chi(degree))
404 : 12 : hubparerr(degree)=sqrt((chi0err(degree)/chi0(degree)**2)**2+(chierr(degree)/chi(degree)**2)**2)
405 [ - + ]: 12 : ABI_FREE(chi0coeffs)
406 [ - + ]: 16 : ABI_FREE(chicoeffs)
407 : : end do
408 : :
409 : : !##########################################################################################################
410 : : !############################# Printing information on Hubbard Parameters ################################
411 : :
412 : : !Printing relevant information about the Hubbard parameter just calculated.
413 : 4 : write(message,'(3a)') ch10,ch10, &
414 : 8 : '*************************************************************************************************'
415 : 4 : call wrtout(std_out,message)
416 : : write(message,'(4a)') &
417 : 4 : '************************************** Linear Response ',parname,' **************************************',ch10
418 : 4 : call wrtout(std_out,message)
419 : 4 : write(message, '(a,i4)' ) ' Total number of atoms: ',natom
420 : 4 : call wrtout(std_out,message)
421 : 4 : write(message, '(a,i4)' ) ' Index of perturbed atom: ',pawujat
422 : 4 : call wrtout(std_out,message)
423 : 4 : write(message, '(a,i4)' ) ' Value of macro_uj: ',macro_uj
424 : 4 : call wrtout(std_out,message)
425 : 4 : write(message, '(a,i4)' ) ' Value of dmatpuopt: ',dmatpuopt
426 : 4 : call wrtout(std_out,message)
427 : 4 : write(message, '(a,f6.3)' ) ' Mixing constant factored out of Chi0: ',diem
428 : 4 : call wrtout(std_out,message)
429 : : write(message, '(a,f12.5,2a)' )&
430 : 4 : ' Percentage of AE orbital within the PAW sphere of perturbed subspace: ',ph0phiint*100.00,'%',ch10
431 : 4 : call wrtout(std_out,message)
432 : :
433 : 4 : write(message, fmt='(10a)')' Perturbations ',occmag,ch10,&
434 : 4 : ' --------------- -----------------------------',ch10,&
435 [ + - ]: 4 : ' ',trim(pertname),' [eV] Unscreened Screened',ch10,&
436 : 8 : ' --------------- -----------------------------'
437 : 4 : call wrtout(std_out,message)
438 [ + + ]: 22 : do ipert=1,nfiles
439 [ + + ]: 18 : if ((perts(ipert)>0.0d0).and.(perts(ipert-1)<0.0d0)) then
440 : 4 : write(message, fmt='(3f15.10)') perts(0),occs0(0),occs(0)
441 : 4 : call wrtout(std_out,message)
442 : : end if
443 : 18 : write(message, fmt='(3f15.10)') perts(ipert),occs0(ipert),occs(ipert)
444 : 22 : call wrtout(std_out,message)
445 : : end do
446 : :
447 : 4 : write(message, fmt='(11a)') ' RMS Errors',&
448 : 4 : ch10,' ---------------------------------------',ch10,&
449 : 4 : ' Regression Chi0 [eV^-1] Chi [eV^-1] ',parname,' [eV] | Chi0 [eV^-1] Chi [eV^-1] ',parname,&
450 : 8 : ' [eV]',ch10,'--------------------------------------------------------|---------------------------------------'
451 : 4 : call wrtout(std_out,message)
452 [ + + ]: 16 : do degree=1,mdegree
453 [ + + ]: 12 : if (degree==1) then
454 : 4 : regname=' Linear: '
455 [ + + ]: 8 : else if (degree==2) then
456 : 4 : regname=' Quadratic: '
457 [ + - ]: 4 : else if (degree==3) then
458 : 4 : regname=' Cubic: '
459 : : else
460 : 0 : write(degreename,'(i2)') degree
461 [ # # # # ]: 0 : regname=' Degree'//trim(degreename)//' : '
462 : : end if
463 : 12 : write(message,fmt='(a,3f14.7,a,3f13.7)') regname,chi0(degree),chi(degree),hubpar(degree),&
464 : 24 : ' |',chi0err(degree),chierr(degree),hubparerr(degree)
465 : 16 : call wrtout(std_out,message)
466 : : end do
467 : :
468 : : write(message,'(3a)') &
469 : 4 : '*************************************************************************************************',ch10,&
470 : 8 : '*************************************************************************************************'
471 : 4 : call wrtout(std_out,message)
472 : :
473 : :
474 : : !##########################################################################################################
475 : : !############################################ Deallocations ##############################################
476 : :
477 : :
478 [ - + ]: 4 : ABI_FREE(perts)
479 [ - + ]: 4 : ABI_FREE(occs0)
480 [ - + ]: 4 : ABI_FREE(occs)
481 [ - + ]: 4 : ABI_FREE(luocc_nnat)
482 [ - + ]: 4 : ABI_FREE(pawujat_file)
483 [ - + ]: 4 : ABI_FREE(diem_file)
484 [ - + ]: 4 : ABI_FREE(dmatpuopt_file)
485 [ - + ]: 4 : ABI_FREE(ph0phiint_file)
486 [ - + ]: 4 : ABI_FREE(macrouj_file)
487 [ - + ]: 4 : ABI_FREE(nspden_file)
488 [ - + ]: 4 : ABI_FREE(chi0err)
489 [ - + ]: 4 : ABI_FREE(chierr)
490 [ - + ]: 4 : ABI_FREE(chi0)
491 [ - + ]: 4 : ABI_FREE(chi)
492 [ - + ]: 4 : ABI_FREE(hubpar)
493 [ - + ]: 4 : ABI_FREE(hubparerr)
494 [ - + ]: 4 : ABI_FREE(luocc)
495 [ - + ]: 4 : ABI_FREE(uj_perts)
496 [ - + ]: 4 : ABI_FREE(file_paths)
497 : :
498 : : !Ending herald.
499 : 4 : write(std_out,*) ch10,'Linear Response UJ (LRUJ) program complete. Live long and prosper. ~LMac',ch10
500 : :
501 : 4 : call ydoc%write_and_free(std_out)
502 : :
503 : : ! Writes information on file about the memory before ending mpi module, if memory profiling is enabled
504 : 4 : call abinit_doctor("__lruj")
505 : :
506 : 4 : 100 call xmpi_end()
507 : :
508 : : !##########################################################################################################
509 : : !##################################### Subroutines and Functions ########################################
510 : :
511 : : contains
512 : :
513 : : ! Show command line help
514 : 0 : subroutine lruj_show_help()
515 : :
516 : 0 : write(std_out,"(a)")" "
517 : 0 : write(std_out,"(a)")" Linear Response Hubbard U and Hund's J (LRUJ) Utility"
518 : 0 : write(std_out,"(a)")"-----------------------------------------------------------------------------"
519 : 0 : write(std_out,"(a)")"To execute the LRUJ utility, execute: "
520 : 0 : write(std_out,"(a)")" ./lruj FILE1 FILE2 FILE3 ... [options]"
521 : 0 : write(std_out,"(2a)")" ^ input files must be _LRUJ.nc from Abinit run",ch10
522 : 0 : write(std_out,"(a)")" --version Show version number and exit."
523 : 0 : write(std_out,"(a)")" -h, --help Show this help and exit."
524 : 0 : write(std_out,"(a)")" --d <n> Set the maximum degree n polynomial calculated for"
525 : 0 : write(std_out,"(a)")" the response functions chi and chi0."
526 : 0 : write(std_out,"(a)")" (i.e., 1=linear, 2=quadratic, 3=cubic, etc.)"
527 : 0 : write(std_out,"(a)")" NOTE: a degree n polynomial will require at minimum"
528 : 0 : write(std_out,"(a)")" n+2 points (n+1 perturbations and the unperturbed"
529 : 0 : write(std_out,"(a)")" case) or more so as to avoid overfitting."
530 : :
531 : 4 : end subroutine lruj_show_help
532 : :
533 : : ! Function to simplify reading in of variables from netcdf files.
534 : 144 : integer function vid(vname)
535 : : character(len=*),intent(in) :: vname
536 : 288 : vid = nctk_idname(ncid, vname)
537 : 144 : end function vid
538 : :
539 : : end program lruj
540 : : !!***
|
