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