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 : !!***
|
