diff --git a/CHANGES b/CHANGES index acaa5ae6a..3cbd46742 100644 --- a/CHANGES +++ b/CHANGES @@ -4,6 +4,9 @@ Following is a brief summary of changes made in each release of Libint. - unreleased: + - PR #447: TennoGmEval (Yukawa/Slater core integrals) fixed for U < 1e-7 at small T (#444) and no + longer returns NaN for large U and T; the slow, thread-unsafe Gaussian-fit fallback + is replaced by series expansions - PR #449: generated and generic RR functions declare their target LIBINT2_RESTRICT (new config.h macro), up to ~1.2x faster high-am and derivative ERIs - PR #446: a value of 0 for LIBINT2_OPT_AM (or for the per-class LIBINT2__MAX_AM / diff --git a/export/tests/unit/test-core-ints.cc b/export/tests/unit/test-core-ints.cc index 995beb09f..4c9f8e05f 100644 --- a/export/tests/unit/test-core-ints.cc +++ b/export/tests/unit/test-core-ints.cc @@ -421,7 +421,8 @@ TEST_CASE("Slater/Yukawa core integral values", "[core-ints]") { } } - // eval at T=0, U=10^-8 => outside of supported U range, use Gaussian fit + // eval at T=0, U=10^-8 => U below the interpolation range, use the + // small-U expansion { const std::vector ref_values{ 0.9998227746131371278708362646069674307491, @@ -463,7 +464,7 @@ TEST_CASE("Slater/Yukawa core integral values", "[core-ints]") { 0.01369863013312753231830929991270401008582}; tenno_eval->eval_yukawa(yukawa, 4e-8, 0, mmax, zeta); for (int m = 0; m <= mmax; ++m) { - CHECK(std::abs(yukawa[m] - ref_values.at(m)) <= 1.4e-11); + CHECK(std::abs(yukawa[m] - ref_values.at(m)) <= 1.4e-14); } } { @@ -507,7 +508,714 @@ TEST_CASE("Slater/Yukawa core integral values", "[core-ints]") { 7.717538100870460006421265480944800459593e-12}; tenno_eval->eval_slater(stg, 4e-8, 0, mmax, zeta); for (int m = 0; m <= mmax; ++m) { - CHECK(std::abs(stg[m] - ref_values.at(m)) <= 1.4e-11); + CHECK(std::abs(stg[m] - ref_values.at(m)) <= 1.4e-14); + } + } + + // eval at T=10^-3, U=10^-8 => U below the interpolation range at small T + // (issue #444): the upward recursion is unstable here, use the small-U + // expansion about the Boys function + { + const std::vector ref_values{ + 0.9994895412626637181927989321139209904782, + 0.3331333980792379141818065719644582162965, + 0.1998571970647880653193612908486582327537, + 0.1427460766166482567546903558833141777502, + 0.1110202401533734552692202013845432609636, + 0.09083220095381904394579519352024581804025, + 0.07685643951964740158043846972602571565582, + 0.06660786934262397701281409329109961410404, + 0.05877092155672695855203690380381568940329, + 0.05258398157034528452577387674630548779668, + 0.04757558930192272809373602644318475546167, + 0.04343827934096485365344432735071245940245, + 0.03996298016421419980936688672556020878952, + 0.03700257037279519270694202011789602499991, + 0.03445051567740735183159110143721740320497, + 0.03222777574480945199521932524456029979332, + 0.03027447236416686594538102276767381745039, + 0.02854441434354971285371703498325302592057, + 0.02700139817720274768909307691019393546523, + 0.02561664700747984969137279193016858807316, + 0.02436699918355491976924094539024917746628, + 0.0232336023548295102101988300080700092745, + 0.0222009558169655207963642820866659529741, + 0.02125619737274563246377046814425464667853, + 0.02038856484442523423353076979903046084935, + 0.01958898429271661493250613012298490049158, + 0.01884975147183338913098188759467317432731, + 0.01816428278715906415849200007243355243716, + 0.01752691868428324283621471791358061733428, + 0.01693276702775286968242353813960740012941, + 0.0163775772913466402991248967374956021728, + 0.01585763871270356681036873402326031529813, + 0.01536969725063160605384614012091191487579, + 0.01491088741608804300754344949658634471926, + 0.01447867595891742499394009027596195331115, + 0.01407081507254061778586601229389684931815, + 0.01368530329119515339000778891106620267421}; + tenno_eval->eval_yukawa(yukawa, 4e-8, 1e-3, mmax, zeta); + for (int m = 0; m <= mmax; ++m) { + CHECK(std::abs(yukawa[m] - ref_values.at(m)) <= 1.4e-14); + } + } + // eval at T=1, U=10^-8 + { + const std::vector ref_values{ + 0.7466469135079010595229990656238259201744, + 0.1894723402481559596145888890545509681242, + 0.1002687972529819137031320437310866747016, + 0.06673227444145702594162782039287514357657, + 0.04962324096206640252775462332565300610595, + 0.0393648644109003949917041793009862976826, + 0.03256703417046342127727512635224043678795, + 0.02774600191593972161353038612587419840864, + 0.02415529410949709300835022363607731786431, + 0.02138027962246414893261223196123062938245, + 0.01917293606924119515702524863446278177116, + 0.01737610835511418457564471490725109448663, + 0.01588552568982132251456428792290963921795, + 0.01462935071080645418543355971208720557398, + 0.01355651416902122760380439617808972075487, + 0.01262973501138014656546640135590611434358, + 0.01182117222623625265717961197385830038347, + 0.01110962127347435815950317814259553616346, + 0.01047865181829182925590625898648756887622, + 0.009915338163773893171247487764321371913136, + 0.009409373712656274224482418882598908448131, + 0.008952440622885842441866633725022067106999, + 0.00853775290041818882893348225206819790357, + 0.008159719763212494082099890009470356357639, + 0.007813693935149979135767418431157763202282, + 0.007495780907050525660664807303633664827702, + 0.007202692622204182900690492519602203695079, + 0.006931633977647495141041423295376037233017, + 0.006680213871611881802919084549015515381071, + 0.006446374824533810361906975976625987591383, + 0.006228336804828383594612724258927045713049, + 0.006024552026008287083264308435614220291351, + 0.005833668295823250373347666587244748043931, + 0.005654499088779996596120149837366962060983, + 0.005485998946745408130495638204539461192467, + 0.00532724313354041059213759893708247579087, + 0.005177410709823404690576958490653838898368}; + tenno_eval->eval_yukawa(yukawa, 4e-8, 1., mmax, zeta); + for (int m = 0; m <= mmax; ++m) { + CHECK(std::abs(yukawa[m] - ref_values.at(m)) <= 1.4e-14); + } + } + // eval at T=1024, U=10^-8 => at the edge of the T range of the small-U + // expansion + { + const std::vector ref_values{ + 0.02751791228765790220176742068804349405662, + 0.00001352247408510689100383727157248387325126, + 1.980858037771798152810344241497333038265e-8, + 4.836092399319902821207875430678618241734e-11, + 1.652963203730472421114169151447029181582e-13, + 7.264003176640161342691974989635802730212e-16, + 3.901565820425090301711921248504729442037e-18, + 2.476580575856080429813362383521067481352e-20, + 1.813902169969377309442433553547778770983e-22, + 1.505680754035183885384874513527197565862e-24, + 1.396871970433243835989611066498235392917e-26, + 1.43233956984666155401121109210019386258e-28, + 1.608584614268415850036402690308915616983e-30, + 1.963604405341285934777010511725833721922e-32, + 2.588736433661325623818895697724750542467e-34, + 3.665691430829455750650383046290956503087e-36, + 5.548654274081347441055328226211476906504e-38, + 8.940703045706412461608268201452384301589e-40, + 1.527952234848793868745040376765087416888e-41, + 2.760460667848924916252391907322357163834e-43, + 5.256736772559206124581592865022103658188e-45, + 1.052374087870355404324683300638135081203e-46, + 2.209574552234992036025517912292889721492e-48, + 4.855022312271621941497027085085285094894e-50, + 1.114189712382399100196993872314669062724e-51, + 2.665785980358496990428568752091025443873e-53, + 6.63843199348939312874293001712060753256e-55, + 1.717953618098192365919771834787742692421e-56, + 4.613645035306981387864554524507414624378e-58, + 1.284071144767603395042306346369464370923e-59, + 3.699228440625859734809196875425258925095e-61, + 1.101820983624608073039881833425554136571e-62, + 3.389390757236132652670141722053187777784e-64, + 1.07573438195307905113659842779893004105e-65, + 3.519248255309299212567315683425618025439e-67, + 1.185684237147511765494642990378732906998e-68, + 4.110526442447411445897436271013338529393e-70}; + tenno_eval->eval_yukawa(yukawa, 4e-8, 1024., mmax, zeta); + for (int m = 0; m <= mmax; ++m) { + CHECK(std::abs(yukawa[m] - ref_values.at(m)) <= 1.4e-14); + } + } + // same for the Slater-type geminal + { + const std::vector ref_values{ + 0.0001772053770789240503456743416852583782928, + 1.332712286366851608021984720298925548364e-8, + 2.665524020288996977248905622315999670855e-9, + 1.142222408962796171293418699306881100072e-9, + 6.345167292654960297094030899754183357317e-10, + 4.037607839910882264685001572859488584664e-10, + 2.795152286834328473071344758844020476887e-10, + 2.049714035404684913524875286985220310355e-10, + 1.56738955717940369215543789745678494015e-10, + 1.237387997276334805252605411502040321322e-10, + 1.001678453684511286407570060624146467003e-10, + 8.274619921915748880583398184944592118448e-11, + 6.950598353501307688154881250304501225842e-11, + 5.920819582838014204849733215328367579229e-11, + 5.104109390775681750701837361357243589887e-11, + 4.445479865195799672743552385314206823293e-11, + 3.906606761285172099676604953772964685852e-11, + 3.460116041234306183327975568841583059655e-11, + 3.086032332693930329247916146118180910676e-11, + 2.769502339445795995440569960050694784128e-11, + 2.499295647849859844263693079838821213766e-11, + 2.266793657450819118084230764358336383571e-11, + 2.0652930757279788276690958428081126008e-11, + 1.889516888439776665187627884822612591139e-11, + 1.735265056640796460479396690448371658349e-11, + 1.599161103417238602049279352091120715551e-11, + 1.478465641766451603048485056623452328529e-11, + 1.370937369348649944979775044479243780313e-11, + 1.274728205751642644554564317705870205752e-11, + 1.188303313060746307582359547946434409743e-11, + 1.110379472812458766597282804223595913216e-11, + 1.03987715728614697751232542847057374934e-11, + 9.758829241439215130451878046968008446728e-12, + 9.176196690871260926053812486511403130578e-12, + 8.64422914341236027206718441248782816232e-12, + 8.157217727536144161481559641302079860027e-12, + 7.710235626909287917164467656612932878741e-12}; + tenno_eval->eval_slater(stg, 4e-8, 1e-3, mmax, zeta); + for (int m = 0; m <= mmax; ++m) { + CHECK(std::abs(stg[m] - ref_values.at(m)) <= 1.4e-14); + } + } + + { + const std::vector ref_values{ + 0.0001771932269149111436812921867555710432325, + 1.1143491465194901998168203531385499041e-8, + 1.784070859903480918229136906469285868451e-9, + 6.707304562304977552300844667642306225004e-10, + 3.421806695878124682774639413444427494124e-10, + 2.051675310233201507210088804933341684668e-10, + 1.359566048087394742885810589749172178931e-10, + 9.64206450904739932748948045273247675862e-11, + 7.181415612885257210360324979593761088667e-11, + 5.550028974065888151475983349693376963711e-11, + 4.414687106445907551173966653535695222576e-11, + 3.593655428254021162761067454423374569072e-11, + 2.981165330585724122160853968682910537358e-11, + 2.512349958029736658261456421644867287947e-11, + 2.145673083570453163258327067994969638213e-11, + 1.853558315282162076675989644367212822581e-11, + 1.617125570287787816573578764095627920213e-11, + 1.423101905523788995352867662525528440027e-11, + 1.261938910365057807193838312215934574477e-11, + 1.126627309035872169317542444332393926169e-11, + 1.01192890223523789353013776344492693001e-11, + 9.138661795408635652315703151536826822638e-12, + 8.293754449353072258663029459077384068589e-12, + 7.560662744113894936671844851956830918628e-12, + 6.920516561250298926649431566251863107134e-12, + 6.358260561989069502052222550481967491592e-12, + 5.86176569692685519948629568062922265247e-12, + 5.421172891133755192981384484523329241243e-12, + 5.028402120712266762446774927210437038907e-12, + 4.676780941561428820242171447790555793763e-12, + 4.360760394108535345885034353978837566679e-12, + 4.075695576401930226968316466256508433956e-12, + 3.817674603700734198332836967389444948413e-12, + 3.583384140865075544550334997555719658956e-12, + 3.37000284069176931249023265655001737032e-12, + 3.17511626409995076716078534913970803194e-12, + 2.996648474340118031212808928572737850036e-12}; + tenno_eval->eval_slater(stg, 4e-8, 1., mmax, zeta); + for (int m = 0; m <= mmax; ++m) { + CHECK(std::abs(stg[m] - ref_values.at(m)) <= 1.4e-14); + } + } + // eval at U=0 (zeta=0): Yukawa reduces to Coulomb, i.e. G_m(T,0) = F_m(T) + { + auto boys = new scalar_type[mmax + 1]; + auto fm_eval = libint2::FmEval_Chebyshev7::instance(mmax); + for (const scalar_type T : {0., 1e-3, 1., 1024., 2000.}) { + fm_eval->eval(boys, T, mmax); + tenno_eval->eval_yukawa(yukawa, one_over_rho, T, mmax, 0.); + for (int m = 0; m <= mmax; ++m) { + CHECK(std::abs(yukawa[m] - boys[m]) <= + std::numeric_limits::epsilon()); + } + } + delete[] boys; + } + // U + 2 sqrt(U T) > 709 or U > Umax: used to be computed via the Gaussian + // fit of the kernel, now by overflow-safe closed forms for G_{-1} and G_0, + // and by the asymptotic expansion in 1/(U-T) for U > Umax. For the Slater + // the reference values are G_{m-1} - G_m, i.e. stg[m] / (zeta/(2 rho)), with + // zeta = 1, rho = 1/(4 U) => zeta/(2 rho) = 2 U + // T=0.3, U=800.: U + 2 sqrt(U T) > 709 inside the interpolation range + { + const scalar_type T = 0.3, U = 800.; + const std::vector ref_yukawa{ + 0.000462318558174057428145724240231768592938, + 0.0004617429589755465886073842825157757786696, + 0.0004611687895010978868452631676003397575046, + 0.0004605960444369418986123142423767874610295, + 0.000460024718495576959722480898940420087213, + 0.0004594548064156072505522276258447437498464, + 0.0004588863029615820752836047852457476780284, + 0.0004583192029238363256288070699464123105504, + 0.0004577535011183321188766385409593288692063, + 0.000457189192386501600200646321252359124746, + 0.0004566262715950908992669438681735698634124, + 0.0004560647336360052312769265293378242382611, + 0.0004555045734261551326761998911081112626973, + 0.0004549457859073038218561081673411510048647, + 0.0004543883660459156752682782889537087747552, + 0.0004538323088330058094655980128038167819245, + 0.0004532776092839907596750356750591461096202, + 0.0004527242624385402455986974087363341906237, + 0.0004521722633604300152295168043143331365264, + 0.0004516216071373957575569940332469314451553, + 0.0004510722888809880751264581362441210350506, + 0.0004505243037264285075024291054706708389617, + 0.0004499776468324665967728170133593832900546, + 0.0004494323133812379863159250607145303481223, + 0.0004488882985781235441365331846321971113595, + 0.0004483455976516095021607397873489054499456, + 0.000447804205853148602961742080821139900264, + 0.0004472641184570222454703512065871040745985, + 0.0004467253307602036213047772638294872383705, + 0.0004461878380822218334340920997173662707672, + 0.0004456516357650269889687944878121986430935, + 0.0004451167191728562579500733108554637898592, + 0.0004445830836921008900866996093252925474445, + 0.0004440507247311741814649877617820119607014, + 0.0004435196377203803833329594034100084933132, + 0.0004429898181117845451347306144880849642982, + 0.0004424612613790832840452329447512493066665}; + const std::vector ref_slater{ + 5.770342627732822242085632976655127223561e-7, + 5.755991985108395383399577159928142683746e-7, + 5.741694744487017621211149154360211649835e-7, + 5.727450641559882329489252235522964750813e-7, + 5.713259413649388898333434363673738165526e-7, + 5.699120799697091702532730956763373665916e-7, + 5.685034540251752686228405989960718179326e-7, + 5.671000377457496547977152993353674780802e-7, + 5.657018055042067521685289870834413440412e-7, + 5.643087318305186759922197069697444603613e-7, + 5.629207914107009337024530787892613335515e-7, + 5.615379590856679900173388357456251513302e-7, + 5.601602098500986007266382297129755637366e-7, + 5.587875188513108200917237669602578326584e-7, + 5.574198613881465878298783874422301094623e-7, + 5.560572129098658026802761498919928307539e-7, + 5.546995490150497905623377446706723042623e-7, + 5.533468454505140763382663228119189965442e-7, + 5.519990781102303691806044220010540972781e-7, + 5.50656223034257672522771067401691371043e-7, + 5.493182564076824305358970028104101047638e-7, + 5.479851545595676240290307734501960889147e-7, + 5.466568939619107296120921112875489070923e-7, + 5.453334512286104568919526448529419323234e-7, + 5.440148031144421793918760823332367627635e-7, + 5.427009265140419757933972832916614138515e-7, + 5.413917984608991989977065277655496816646e-7, + 5.400873961263574913908742340358256654781e-7, + 5.387876968186241655739427576168362279668e-7, + 5.37492677981787870685164112120967603363e-7, + 5.362023171948444652976119051676276736736e-7, + 5.349165921707310187211769567348532342552e-7, + 5.336354807553678633737015301712424147072e-7, + 5.323589609267086217118475432805867431407e-7, + 5.31087010793798132028358372003467388212e-7, + 5.29819608595838198228788921923529015015e-7, + 5.285567327012610894976697368356576316509e-7}; + tenno_eval->eval_yukawa(yukawa, 4 * U, T, mmax, 1.); + tenno_eval->eval_slater(stg, 4 * U, T, mmax, 1.); + for (int m = 0; m <= mmax; ++m) { + // interpolation range: absolute precision of the table + CHECK(std::abs(yukawa[m] - ref_yukawa.at(m)) <= 1.4e-14); + CHECK(std::abs(stg[m] / (2 * U) - ref_slater.at(m)) <= 1.4e-14); + } + } + // T=5., U=1000.: same, at the U edge of the interpolation range + { + const scalar_type T = 5., U = 1000.; + const std::vector ref_yukawa{ + 3.380777475199456148771987531212095238007e-6, + 3.377391669675422760741969607009080056501e-6, + 3.374012632247146919015254809679346733329e-6, + 3.370640342661415944296706491813259781938e-6, + 3.367274780745665300453814165296203628222e-6, + 3.363915926407577966144888796576114768059e-6, + 3.360563759634686189917368420632027004353e-6, + 3.357218260493975612265495947202163832595e-6, + 3.353879409131491738267085732366221734621e-6, + 3.350547185771948744548392870612918582955e-6, + 3.347221570718340604454250612566467346867e-6, + 3.34390254435155451542765809413087313452e-6, + 3.340590087129986612728893814952052697905e-6, + 3.33728417958915995374901104871233376377e-6, + 3.333984802341344757296250507091415858244e-6, + 3.330691936075180882356493898189433858026e-6, + 3.32740556155530253095039018782799172378e-6, + 3.324125659621965159830224942876719408781e-6, + 3.320852211190674585878982550824437801745e-6, + 3.317585197251818270192381698551876737662e-6, + 3.314324598870298765941956474397454740951e-6, + 3.311070397185169315233518940562487085415e-6, + 3.307822573409271580290584009067217772534e-6, + 3.304581108828875494406573838457472174503e-6, + 3.301345984803321218222856539351248842098e-6, + 3.298117182764663187001922419473129341931e-6, + 3.294894684217316234676269894720303178463e-6, + 3.291678470737703780563872021801050346997e-6, + 3.288468523973908064750432749123987716158e-6, + 3.285264825645322418247028715374374496477e-6, + 3.282067357542305554139176930663927875948e-6, + 3.278876101525837866049880037082434453791e-6, + 3.275691039527179720344788051562487363456e-6, + 3.272512153547531728612287436800633735739e-6, + 3.269339425657696987055093824219293835749e-6, + 3.266172837997745269532792432397449729417e-6, + 3.26301237277667916109674879903823534377e-6}; + const std::vector ref_slater{ + 3.392593954054785275360534631556239081609e-9, + 3.385805524033388030017924203015181505526e-9, + 3.37903742827584172671479732973332317204e-9, + 3.372289585730974718548317866086951391031e-9, + 3.365561915750643842892326517056153715857e-9, + 3.358854338087334308925368720088860162761e-9, + 3.352166772891776227520375944087763706039e-9, + 3.345499140710577651872473429863171758143e-9, + 3.338851362483873998410214835942097974236e-9, + 3.332223359542993718692861753303151665769e-9, + 3.325615053608140094142258046451236087945e-9, + 3.31902636678608902659251843559421234751e-9, + 3.312457221567902698764279178820436614717e-9, + 3.305907540826658979882766239718934134732e-9, + 3.299377247815196452760541620917905526098e-9, + 3.292866266163874939756608901982000217955e-9, + 3.286374519878351406103710361442134246475e-9, + 3.279901933337371120165244951272314999034e-9, + 3.273448431290573951242392052281607035688e-9, + 3.267013938856315686600852272561064082759e-9, + 3.260598381519504250425224154421996711336e-9, + 3.254201685129450708437533834967655535754e-9, + 3.247823775897734942934931495269312880822e-9, + 3.241464580396085884010170609745598031413e-9, + 3.235124025554276183717299106223332405166e-9, + 3.228802038658031220934119878119500167379e-9, + 3.222498547346952325652524752826163467325e-9, + 3.216213479612454112397872919252831466513e-9, + 3.209946763795715813439272677062630838214e-9, + 3.203698328585646503404033749613219681104e-9, + 3.197468103016864107851784710446620529692e-9, + 3.191256016467688089296893581493422156497e-9, + 3.185061998658145705091985519947090334857e-9, + 3.178885979647991732500614761853627717771e-9, + 3.172727889834741557193612581339899989197e-9, + 3.166587659951717522301391821844106332e-9, + 3.160465221066108436043633359214385646836e-9}; + tenno_eval->eval_yukawa(yukawa, 4 * U, T, mmax, 1.); + tenno_eval->eval_slater(stg, 4 * U, T, mmax, 1.); + for (int m = 0; m <= mmax; ++m) { + // interpolation range: absolute precision of the table + CHECK(std::abs(yukawa[m] - ref_yukawa.at(m)) <= 1.4e-14); + CHECK(std::abs(stg[m] / (2 * U) - ref_slater.at(m)) <= 1.4e-14); + } + } + // T=0.3, U=1e4: U > Umax, asymptotic expansion in 1/(U-T) + { + const scalar_type T = 0.3, U = 1e4; + const std::vector ref_yukawa{ + 0.00003703646699129072631570988572577535441631, + 0.0000370327645292388267471979459347530421356, + 0.00003702906280729454594254839249079211734362, + 0.0000370253618252360113298022328298851940676, + 0.00003702166158284143900447687972445174758839, + 0.00003701796207988913368528199567845816343258, + 0.00003701426331615748866986187280186563878238, + 0.0000370105652914249857905643296177496226342, + 0.00003700686800547019537023610627024634153011, + 0.00003700317145807177617804473961628051303718, + 0.00003699947564900847538532689969881361263646, + 0.00003699578057805912852146316911412404456714, + 0.00003699208624500265942977924680038928670166, + 0.00003698839264961808022347355778958554693593, + 0.00003698469979168449124157125047945269708564, + 0.00003698100767098108100490456299699125308819, + 0.00003697731628728712617211954023966396061043, + 0.00003697362564038199149570908319516713613065, + 0.00003696993573004512977807231215531631835912, + 0.0000369662465560560818276002254542570166332, + 0.0000369625581181944764147876353758644148003, + 0.00003695887041624003022837136288983581420512, + 0.00003695518344997254783149467289060639082758, + 0.00003695149721917192161789793162783251246501, + 0.00003694781172361813176813546803178742519119, + 0.00003694412696309124620581862065160158721434, + 0.00003694044293737142055388495193885431574472, + 0.00003693675964623889809089361162358473159915, + 0.0000369330770894740097073468309443382500338, + 0.00003692939526685717386203752950840008771009, + 0.00003692571417816889653842301657288944775081, + 0.00003692203382318977120102476855189722050921, + 0.00003691835420170047875185426456934620991304, + 0.000036914675313481787486864861891736077005, + 0.00003691099715831455305242969308940539651237, + 0.00003690731973597971840184556678940046186023, + 0.00003690364304625831375186285389748476189788}; + const std::vector ref_slater{ + 3.703202381479616120433684207243016343284e-9, + 3.702462051899568511939791022312280713845e-9, + 3.701721944280804649553443960924791982133e-9, + 3.700982058534612746159660906923276015789e-9, + 3.700242394572325325353105433446479212652e-9, + 3.699502952305319194884045993584155805074e-9, + 3.698763731645015420122876592524650200226e-9, + 3.698024732502879297543184116016148179726e-9, + 3.697285954790420328223347503281104096753e-9, + 3.696547398419192191366653965828492925042e-9, + 3.695809063300792717839917466900400719895e-9, + 3.695070949346863863730584689568069317834e-9, + 3.694333056469091683922313734757865481938e-9, + 3.693595384579206305689010803739765730601e-9, + 3.692857933588981902307310132849850286842e-9, + 3.692120703410236666687482461443997454918e-9, + 3.691383693954832785022757327292477754561e-9, + 3.690646905134676410457044496824479787521e-9, + 3.689910336861717636771039850817771525348e-9, + 3.689173989047950472086701059301725923726e-9, + 3.688437861605412812590078392601832901737e-9, + 3.687701954446186416272486028600595171904e-9, + 3.686966267482396876689999229423377548945e-9, + 3.686230800626213596741262773878362564942e-9, + 3.685495553789849762463596045087273822237e-9, + 3.684760526885562316847380185837976851305e-9, + 3.684025719825651933668712747271469621119e-9, + 3.683291132522462991340315269584145568387e-9, + 3.682556764888383546780679246481565347207e-9, + 3.681822616835845309301435938162323712447e-9, + 3.681088688277323614512935510639959281924e-9, + 3.680354979125337398248020992227241600628e-9, + 3.679621489292449170503982551010596164148e-9, + 3.678888218691264989402677610132908040473e-9, + 3.67815516723443443516880233068049263504e-9, + 3.677422334834650584126300004934652131576e-9, + 3.676689721404649982712891915699962352805e-9}; + tenno_eval->eval_yukawa(yukawa, 4 * U, T, mmax, 1.); + tenno_eval->eval_slater(stg, 4 * U, T, mmax, 1.); + for (int m = 0; m <= mmax; ++m) { + // asymptotic expansion: full relative precision + CHECK(std::abs(yukawa[m] - ref_yukawa.at(m)) <= + 1e-14 * std::abs(ref_yukawa.at(m))); + CHECK(std::abs(stg[m] / (2 * U) - ref_slater.at(m)) <= + 1e-14 * std::abs(ref_slater.at(m))); + } + } + // T=10., U=2000.: U > Umax + { + const scalar_type T = 10., U = 2000.; + const std::vector ref_yukawa{ + 1.139837266419768490071553164005174761447e-8, + 1.13926549453786311281203703555364039419e-8, + 1.138694295710119803274860553614947158904e-8, + 1.138123669075894835268544446378764438228e-8, + 1.137553613776266169358523476468948144513e-8, + 1.136984128954029151962647037410863121099e-8, + 1.136415213753692227326573361612552429644e-8, + 1.135846867321472662334102364467732109482e-8, + 1.135279088805292284107671293108233821321e-8, + 1.134711877354773230354415689935369454916e-8, + 1.134145232121233712413375724332314056714e-8, + 1.133579152257683790959604694869769553257e-8, + 1.13301363691882116432111246280999513944e-8, + 1.132448685261026969364751749713353386112e-8, + 1.131884296442361594907329621359003769732e-8, + 1.131320469622560507608400090888181498881e-8, + 1.130757203963030090301366609924384080132e-8, + 1.130194498626843492719695281258482318738e-8, + 1.129632352778736494575210924326108894543e-8, + 1.12907076558510338094561865894671401299e-8, + 1.128509736213992829928563447414820044408e-8, + 1.127949263835103812519709053790132499446e-8, + 1.127389347619781504672486145859742565867e-8, + 1.126829986741013211497326783457869472825e-8, + 1.126271180373424303558369310321455245005e-8, + 1.125712927693274165225783699108408725742e-8, + 1.125155227878452155042032694264440062029e-8, + 1.124598080108473578060548658729460123209e-8, + 1.124041483564475670115469861640976555075e-8, + 1.123485437429213593981243048815756634092e-8, + 1.122929940887056447381061519448741895991e-8, + 1.122374993123983282803269594716938411489e-8, + 1.121820593327579139085025310353684004791e-8, + 1.121266740687031084722673399284104123717e-8, + 1.12071343439312427286844015558549858312e-8, + 1.120160673638238007973220590837743665506e-8, + 1.119608457616341824035386411820655446857e-8}; + const std::vector ref_slater{ + 5.723457984371717597770218953574027817154e-12, + 5.717718819053772595161284515343672563633e-12, + 5.711988277433095371764819386932352868457e-12, + 5.706266342249680063161072361827206753257e-12, + 5.700552996286659100209699098162937151428e-12, + 5.694848222370173958764390580850234141311e-12, + 5.689152003369246360736757983106914550205e-12, + 5.683464322195649924709971448203201620765e-12, + 5.677785161803782264310713594982881612533e-12, + 5.672114505190537532556031728643664041794e-12, + 5.666452335395179410399656030553982025629e-12, + 5.660798635499214537710294625445034569173e-12, + 5.655153388626266384922320597744138166499e-12, + 5.649516577941949563607130966417533285631e-12, + 5.643888186653744574221283543496163795095e-12, + 5.638268198010872989295304708222708517289e-12, + 5.632656595304173070334809637974187488494e-12, + 5.627053361865975816713286659017613933318e-12, + 5.621458481069981444843569323734241950354e-12, + 5.615871936331136295922653793948815533429e-12, + 5.610293711105510170552115318939685817547e-12, + 5.604723788890174088543936246875449621712e-12, + 5.599162153223078472229079303899335787933e-12, + 5.593608787682931751593624018730930419259e-12, + 5.588063675889079389574731364142278204665e-12, + 5.582526801501383325856112130465192627782e-12, + 5.576998148220101837510048439686637128988e-12, + 5.571477699785769814840355349799388206191e-12, + 5.565965439979079450787970884835681335334e-12, + 5.560461352620761342268128252199209835235e-12, + 5.554965421571466001815293670147381005714e-12, + 5.54947763073164577791924731803484502285e-12, + 5.543997964041437182442843632544066973814e-12, + 5.538526405480543623519110695798810746265e-12, + 5.53306293906811854233243698605540597159e-12, + 5.527607548862648952195647477549176138317e-12, + 5.522160218961839378341790170882186489509e-12}; + tenno_eval->eval_yukawa(yukawa, 4 * U, T, mmax, 1.); + tenno_eval->eval_slater(stg, 4 * U, T, mmax, 1.); + for (int m = 0; m <= mmax; ++m) { + // asymptotic expansion: full relative precision + CHECK(std::abs(yukawa[m] - ref_yukawa.at(m)) <= + 1e-14 * std::abs(ref_yukawa.at(m))); + CHECK(std::abs(stg[m] / (2 * U) - ref_slater.at(m)) <= + 1e-14 * std::abs(ref_slater.at(m))); + } + } + // T=0., U=1e6: U > Umax, T = 0 + { + const scalar_type T = 0., U = 1e6; + const std::vector ref_yukawa{ + 4.999992500018749934375295310875791807343e-7, + 4.999987500043749803126082805461771536661e-7, + 4.999982500078749566877815291385335724467e-7, + 4.999977500123749195631032761221580894494e-7, + 4.99997250017874865938639520424566791935e-7, + 4.999967500243747928144682605832845419901e-7, + 4.99996250031874697190679494685847676447e-7, + 4.999957500403745760673752203098070667819e-7, + 4.999952500498744264446694344627315389893e-7, + 4.999947500603742453226881335222116534296e-7, + 4.999942500718740297015693131758638446484e-7, + 4.999937500843737765814629683613349211634e-7, + 4.999932500978734829625310932063069252182e-7, + 4.999927501123731458449476809685023524989e-7, + 4.999922501278727622288987239756897318117e-7, + 4.999917501443723291145822135656895647197e-7, + 4.999912501618718435022081400263806251342e-7, + 4.999907501803713023919984925357066188607e-7, + 4.999902501998707027841872591016832030954e-7, + 4.999897502203700416790204265024053658691e-7, + 4.999892502418693160767559802260551654387e-7, + 4.999887502643685229776639044109098296202e-7, + 4.999882502878676593820261817853502150638e-7, + 4.999877503123667222901367936078696264666e-7, + 4.999872503378657087023017196070829957213e-7, + 4.999867503643646156188389379217364209978e-7, + 4.99986250391863440040078425040717065756e-7, + 4.999857504203621789663621557430634176866e-7, + 4.999852504498608293980441030379759075773e-7, + 4.999847504803593883354902381048278881023e-7, + 4.999842505118578527790785302331769725329e-7, + 4.999837505443562197291989467627767333653e-7, + 4.999832505778544861862534530235887608648e-7, + 4.999827506123526491506560122757950815227e-7, + 4.999822506478507056228325856498109364239e-7, + 4.999817506843486526032211320862979195221e-7, + 4.999812507218464870922716082761774758215e-7}; + const std::vector ref_slater{ + 4.99998500005624973750147655275476140132e-13, + 4.999975000131249212505414020270682006176e-13, + 4.999965000236248267514076435812193819563e-13, + 4.999955000371246782530163754829972506346e-13, + 4.99994500053624463755697591297514417238e-13, + 4.999935000731241712598412822499449159599e-13, + 4.999925000956237887658974368655431039589e-13, + 4.999915001211233042743760406096650805446e-13, + 4.999905001496227057858470755277926261725e-13, + 4.999895001811219813009405198855596612268e-13, + 4.999885002156211188203463478087812245712e-13, + 4.999875002531201063448145289234849718485e-13, + 4.99986500293618931875155027995945193507e-13, + 4.999855003371175834122378045727193525353e-13, + 4.999845003836160489569928126206871418848e-13, + 4.999835004331143165104100001670920615587e-13, + 4.999825004856123740735393089395855153492e-13, + 4.999815005411102096474906740062734272013e-13, + 4.999805005996078112334340234157653771834e-13, + 4.99979500661105166832599277837226257045e-13, + 4.999785007256022644462763502004304453409e-13, + 4.999775007930990920758151453358185021015e-13, + 4.999765008635956377226255596145563830302e-13, + 4.999755009370918893881774805885971732058e-13, + 4.999745010135878350740007866307453402718e-13, + 4.999735010930834627816853465747235070915e-13, + 4.999725011755787605128810193552417438481e-13, + 4.999715012610737162692976536480693795706e-13, + 4.999705013495683180527050875101093330659e-13, + 4.999695014410625538649331480194749632345e-13, + 4.999685015355564117078716509155694387526e-13, + 4.999675016330498795834704002391676270989e-13, + 4.999665017335429454937391879725005029052e-13, + 4.999655018370355974407477936793420756128e-13, + 4.999645019435278234266259841450988364126e-13, + 4.999635020530196114535635130169017244498e-13, + 4.999625021655109495238101204437006122724e-13}; + tenno_eval->eval_yukawa(yukawa, 4 * U, T, mmax, 1.); + tenno_eval->eval_slater(stg, 4 * U, T, mmax, 1.); + for (int m = 0; m <= mmax; ++m) { + // asymptotic expansion: full relative precision + CHECK(std::abs(yukawa[m] - ref_yukawa.at(m)) <= + 1e-14 * std::abs(ref_yukawa.at(m))); + CHECK(std::abs(stg[m] / (2 * U) - ref_slater.at(m)) <= + 1e-14 * std::abs(ref_slater.at(m))); + } + } + // eval at T > Tmax with U + 2 sqrt(U T) > 709 => the closed-form G_{-1} + // and G_0 that seed the upward recursion used to overflow (and return NaN); + // the exact values are below exp(U - 2 sqrt(U T)), i.e. negligible + { + for (const auto& TU : + {std::make_pair(1100., 1000.), std::make_pair(1280., 125.)}) { + const auto T = TU.first; + const auto U = TU.second; + tenno_eval->eval_yukawa(yukawa, 4 * U, T, mmax, 1.); + tenno_eval->eval_slater(stg, 4 * U, T, mmax, 1.); + for (int m = 0; m <= mmax; ++m) { + CHECK(std::isfinite(yukawa[m])); + CHECK(std::abs(yukawa[m]) <= 1.4e-14); + CHECK(std::isfinite(stg[m])); + CHECK(std::abs(stg[m]) <= 1.4e-14); + } } } @@ -544,6 +1252,44 @@ TEST_CASE("Slater/Yukawa core integral values", "[core-ints]") { } } +// the TennoGmEval singleton is shared by all threads; this exercises every +// evaluation path, including the Gaussian-fit fallback (large U), which must +// not use state shared between threads +TEST_CASE("TennoGmEval shared across threads", "[core-ints][threads]") { + using scalar_type = libint2::scalar_type; + const int mmax = 20; + const auto tenno_eval = libint2::TennoGmEval::instance(mmax); + + // evaluates Yukawa and Slater Gm on a grid of (T, U) values + auto compute = [&]() { + std::vector result; + std::vector Gm(mmax + 1); + auto append = [&]() { result.insert(result.end(), Gm.begin(), Gm.end()); }; + for (scalar_type U : {1e-8, 1e-3, 1., 800., 1e4}) { + for (scalar_type T : {0., 0.3, 10., 1023.5, 1100.}) { + const scalar_type zeta = 1., one_over_rho = 4 * U; + tenno_eval->eval_yukawa(Gm.data(), one_over_rho, T, mmax, zeta); + append(); + tenno_eval->eval_slater(Gm.data(), one_over_rho, T, mmax, zeta); + append(); + } + } + return result; + }; + + const auto reference = compute(); + for (const auto v : reference) REQUIRE(std::isfinite(v)); + const int nthreads = 4; + std::vector> results(nthreads); + std::vector threads; + for (int t = 0; t != nthreads; ++t) + threads.emplace_back([&, t]() { + for (int rep = 0; rep != 10; ++rep) results[t] = compute(); + }); + for (auto& thread : threads) thread.join(); + for (int t = 0; t != nthreads; ++t) REQUIRE(results[t] == reference); +} + // core evaluators wrapped by GenericGmEval are shared by all threads (see // #431), so evaluating them concurrently must reproduce serial results exactly TEST_CASE("GenericGmEval shared across threads", "[core-ints][threads]") { diff --git a/include/libint2/boys.h b/include/libint2/boys.h index 18e70be2f..9587502a9 100644 --- a/include/libint2/boys.h +++ b/include/libint2/boys.h @@ -900,10 +900,11 @@ constexpr double pow10(long exp) { * @tparam Real real type * @warning only @p Real = double is supported * @warning guarantees absolute precision of only about 1e-14 - * @warning for vast majority of practical uses the integral is computed - * by interpolation or closed-form evaluation. But rarely the - * the "analytic" route will cause overflow, so will represent kernel - * as a linear combination of Gaussians to compute without overflow. + * @note depending on \f$ T \f$ and \f$ U \f$ the core integrals are computed + * by Chebyshev interpolation ( \f$ 10^{-7} \leq U \leq 10^3 \f$, + * \f$ T \leq 2^{10} \f$ ), by expansion in \f$ U \f$ about the Boys function + * ( \f$ U < 10^{-7} \f$ ), by asymptotic expansion in \f$ 1/(U-T) \f$ + * ( \f$ U > 10^3 \f$ ), or by upward recursion ( \f$ T > 2^{10} \f$ ). */ template struct TennoGmEval { @@ -926,58 +927,6 @@ struct TennoGmEval { static constexpr std::size_t ORDERp1 = interpolation_order + 1; static constexpr Real maxabsprecision = 1.4e-14; //!< guaranteed abs precision of the interpolation table for m>0 - // map of zeta to existing gaussian fit to slater to avoid - // recomputing - - // to avoid recomputing Gaussian fits of the exponential cache - // them for every value of zeta - using gaussian_fit_t = std::vector>; - using gaussian_fits_t = std::map; // zeta -> fit - // TODO replace by atomic shared_ptr to avoid leaking memory at the end - using gaussian_fits_aptr_t = std::atomic; - static gaussian_fits_aptr_t& gaussian_fits_aptr_accessor() { - static gaussian_fits_aptr_t gaussian_fits_aptr_; - return gaussian_fits_aptr_; - } - static gaussian_fit_t& get_gaussian_fit(double zeta) { - while (gaussian_fits_aptr_accessor().load() == nullptr || - gaussian_fits_aptr_accessor().load()->count(zeta) == 0) { - // remember ptr to the cache we will be adding new fit to - auto* gaussian_fits_ptr = gaussian_fits_aptr_accessor().load(); - - // use existing cache, or create a new one - auto new_gaussian_fits_ptr = new gaussian_fits_t( - gaussian_fits_ptr == nullptr ? gaussian_fits_t() - : *gaussian_fits_ptr); - // constexpr double low_end_range = 1e-5; - // constexpr double high_end_range = 5; - // constexpr double eps = 1e-12; - std::vector> geminal = - fit_slater(zeta /* , low_end_range, high_end_range, eps */); - gaussian_fits_t::iterator it; - bool inserted; - std::tie(it, inserted) = new_gaussian_fits_ptr->emplace(zeta, geminal); - - // inserted may be false if someone inserted between count query and - // grabbing gaussian_fits_ptr ... if inserted new fit try to replace the - // value in gaussian_fits_aptr_ - if (inserted) { - auto* expected_gaussian_fits_ptr = gaussian_fits_ptr; - if (gaussian_fits_aptr_accessor().compare_exchange_strong( - expected_gaussian_fits_ptr, new_gaussian_fits_ptr)) { - // can't delete since someone might be using it ... need atomic - // shared_ptr - // delete gaussian_fits_ptr; - break; - } else { - delete new_gaussian_fits_ptr; - } - } - } - return gaussian_fits_aptr_accessor().load()->at(zeta); - } - mutable std::unique_ptr> yukawa_gaussian_fit_eval_; - mutable std::unique_ptr> slater_gaussian_fit_eval_; public: /// \param m_max maximum value of the Gm function index @@ -1018,6 +967,7 @@ struct TennoGmEval { "TennoGmEval::init() : requested mmax exceeds the " "hard-coded mmax"); init_table(); + fm_eval_ = FmEval_Chebyshev7::instance(mmax_); } ~TennoGmEval() { @@ -1057,27 +1007,23 @@ struct TennoGmEval { /// @param[in] T \f$ T \f$ /// @param[in] mmax \f$ m_\mathrm{max} \f$ /// @param[in] zeta \f$ \zeta \f$ - /// @throw std::invalid_argument if \f$ \zeta^2/(4 \rho) \f$ exceeds \c Umax void eval_yukawa(Real* Gm, Real one_over_rho, Real T, size_t mmax, Real zeta) const { assert(mmax <= mmax_); assert(T >= 0); const auto U = 0.25 * zeta * zeta * one_over_rho; assert(U >= 0); - bool success = false; // did interpolation or recursion work? - if (T > Tmax || U < Umin) { - success = eval_urr(Gm, T, U, /* zeta_over_two_rho = */ 0, mmax); + if (U < Umin && (T <= Tmax || U == 0)) { + // U is below the interpolation range: the upward recursion is unstable + // for small T (#444), expand in U about the Boys function instead; for + // U = 0 the expansion is exact for any T + eval_small_U(Gm, T, U, /* zeta_over_two_rho = */ 0, mmax); + } else if (T > Tmax) { + eval_urr(Gm, T, U, /* zeta_over_two_rho = */ 0, mmax); + } else if (U > Umax) { + eval_large_U(Gm, T, U, /* zeta_over_two_rho = */ 0, mmax); } else { - success = - interpolate_Gm(Gm, T, U, /* zeta_over_two_rho = */ 0, mmax); - } - if (!success) { // fallback is to use the Gaussian fit - auto& gaussian_fit = get_gaussian_fit(zeta); - if (!yukawa_gaussian_fit_eval_) - yukawa_gaussian_fit_eval_ = std::unique_ptr>( - new GaussianGmEval(mmax_, precision_)); - yukawa_gaussian_fit_eval_->eval(Gm, 1.0 / one_over_rho, T, mmax, - gaussian_fit); + interpolate_Gm(Gm, T, U, /* zeta_over_two_rho = */ 0, mmax); } } /// @param[in] Gm pointer to array of @c mmax+1 @c Real elements, on @@ -1091,7 +1037,6 @@ struct TennoGmEval { /// @param[in] T \f$ T \f$ /// @param[in] mmax \f$ m_\mathrm{max} \f$ /// @param[in] zeta \f$ \zeta \f$ - /// @throw std::invalid_argument if \f$ \zeta^2/(4 \rho) \f$ exceeds \c Umax void eval_slater(Real* Gm, Real one_over_rho, Real T, size_t mmax, Real zeta) const { assert(mmax <= mmax_); @@ -1099,24 +1044,24 @@ struct TennoGmEval { const auto U = 0.25 * zeta * zeta * one_over_rho; assert(U > 0); // integral factorizes into 2 overlaps for U = 0 const auto zeta_over_two_rho = 0.5 * zeta * one_over_rho; - bool success = false; // did interpolation or recursion work? - if (T > Tmax) { - success = eval_urr(Gm, T, U, zeta_over_two_rho, mmax); + if (U < Umin && T <= Tmax) { + // U is below the interpolation range: expand in U about the Boys + // function + eval_small_U(Gm, T, U, zeta_over_two_rho, mmax); + } else if (T > Tmax) { + eval_urr(Gm, T, U, zeta_over_two_rho, mmax); + } else if (U > Umax) { + eval_large_U(Gm, T, U, zeta_over_two_rho, mmax); } else { - success = interpolate_Gm(Gm, T, U, zeta_over_two_rho, mmax); - } - if (!success) { // fallback is to use the Gaussian fit - auto& gaussian_fit = get_gaussian_fit(zeta); - if (!slater_gaussian_fit_eval_) - slater_gaussian_fit_eval_ = std::unique_ptr>( - new GaussianGmEval(mmax_, precision_)); - slater_gaussian_fit_eval_->eval(Gm, 1.0 / one_over_rho, T, mmax, - gaussian_fit); + interpolate_Gm(Gm, T, U, zeta_over_two_rho, mmax); } } /// Fits a slater with an undetermined number of gaussians, borrowed from /// https://github.com/m-a-d-n-e-s-s/madness/blob/master/src/madness/mra/gfit.h#L140 + /// @note not used by TennoGmEval itself + /// (eval_yukawa() and eval_slater() no longer need a Gaussian fit of the + /// kernel), kept as a utility /// @param[in] zeta exponent of Slater /// @param[in] lo smallest lengthscale to be represented /// @param[in] hi largest lengthscale to be represented @@ -1181,17 +1126,41 @@ struct TennoGmEval { } private: + /// @param[in] x the argument of erfc + /// @param[in] T the value of \f$ T \f$ + /// @return \f$ \exp(x^2 - T) \mathrm{erfc}(x) \f$ , evaluated without + /// overflow for large \f$ x \f$ (where \f$ \exp(x^2) \f$ overflows and + /// \f$ \mathrm{erfc}(x) \f$ underflows) as \f$ \exp(-T) + /// \mathrm{erfcx}(x) \f$ + static inline Real exp_x2_mT_erfc(Real x, Real T) { + // for x <= 25 exp(x^2 - T) is finite (x^2 - T <= 0 for x < 0 since then + // T > U in all uses) and erfc(x) > 1e-273 is a normal number. + // N.B. x^2 must be computed from x (not as U -/+ 2 sqrt(U T)) so that the + // rounding errors of exp(x^2 - T) erfc(x) for x = sqrt(U) -/+ sqrt(T) are + // correlated: G_0 is their difference divided by sqrt(T) + if (x <= 25) return exp(x * x - T) * erfc(x); + // asymptotic expansion: erfcx(x) = 1/(x sqrt(pi)) sum_n (-1)^n + // (2n-1)!!/(2x^2)^n ; for x > 25 the terms decrease by at least (2n-1)/1250 + // per order, 10 terms give full double precision + const Real oo_sqrtPi(0.56418958354775628694807945156077258584405062932900); + const Real oo_two_x2 = 1 / (2 * x * x); + Real term = 1; + Real sum = 1; + for (int n = 1; n <= 10; ++n) { + term *= -(2 * n - 1) * oo_two_x2; + sum += term; + } + return exp(-T) * oo_sqrtPi / x * sum; + } + /// Computes G_0 and (optionally) G_{-1} /// @tparam compute_Gm1, if true, compute G_0 and G_{-1}, otherwise G_0 only /// @param[in] T the value of \f$ T \f$ - /// @param[in] U the value of \f$ U \f$ - /// @return if {T,U} can be handled without overflow return {G_0,G_{-1}} if - // @c compute_Gm1==true , else {G_0,0}; - // if {T,U} cannot be handled without overflow - // return nullopt + /// @param[in] U the value of \f$ U \f$, must not exceed \c Umax if + /// \f$ T = 0 \f$ + /// @return {G_0,G_{-1}} if @c compute_Gm1==true , else {G_0,0} template - static inline optional> eval_G0_and_maybe_Gm1(Real T, - Real U) { + static inline std::tuple eval_G0_and_maybe_Gm1(Real T, Real U) { assert(T >= 0); assert(U >= 0); const Real sqrtPi(1.7724538509055160272981674833411451827975494561224); @@ -1217,17 +1186,9 @@ struct TennoGmEval { const Real lambda = sqrt_U + sqrt_T; const Real sqrtPi_over_4 = sqrtPi * 0.25; const Real pfac = sqrtPi_over_4; - const Real erfc_k = exp(kappa * kappa - T) * erfc(kappa); - // TODO when lambda * lambda - T exceeds 709 exp overflows, so need to - // handle this more carefully - if (lambda * lambda - T > 709) { - libint2::verbose_stream() - << "TennoGmEval::eval_G0_and_maybe_Gm1(): for large T and/or U " - << "current implementation overflows, will represent " - << "Slater function as a linear combination of Gaussians\n"; - return libint2::nullopt; - } - const Real erfc_l = exp(lambda * lambda - T) * erfc(lambda); + // exp(kappa^2 - T) erfc(kappa) and exp(lambda^2 - T) erfc(lambda) + const Real erfc_k = exp_x2_mT_erfc(kappa, T); + const Real erfc_l = exp_x2_mT_erfc(lambda, T); G_m1 = compute_Gm1 ? pfac * (erfc_k + erfc_l) * oo_sqrt_U : 0.0; G_0 = pfac * (erfc_k - erfc_l) * oo_sqrt_T; @@ -1273,8 +1234,9 @@ struct TennoGmEval { } /// Computes core integrals by upward recursion from \f$ G_{-1} (T,U) \f$ and - /// \f$ G_0 (T,U) \f$ this is unstable for small T, so use when interpolation - /// does not apply. + /// \f$ G_0 (T,U) \f$ this is unstable for small T, so use only for + /// \f$ T > \f$ \c Tmax (for \f$ U < \f$ \c Umin at smaller T use + /// eval_small_U() instead). /// @tparam[in] Exp if true, will compute core integrals of the exponential /// operator, otherwise will compute Yukawa /// @param[out] Gm_vec pointer to an array with at least @c mmax+1 elements, @@ -1288,30 +1250,15 @@ struct TennoGmEval { /// only used if @c Exp==true /// @param[in] mmax the maximum value of \f$ m \f$ template - static bool eval_urr(Real* Gm_vec, Real T, Real U, Real zeta_over_two_rho, + static void eval_urr(Real* Gm_vec, Real T, Real U, Real zeta_over_two_rho, size_t mmax) { - assert(T >= 0); + assert(T > 0); assert(U > 0); - if (T == 0) return false; - - const Real sqrt_U = sqrt(U); - const Real sqrt_T = sqrt(T); - const Real oo_sqrt_T = 1 / sqrt_T; - const Real oo_sqrt_U = 1 / sqrt_U; - const Real kappa = sqrt_U - sqrt_T; - const Real lambda = sqrt_U + sqrt_T; - const Real sqrtPi_over_4( - 0.44311346272637900682454187083528629569938736403060); - const Real pfac = sqrtPi_over_4; - const Real erfc_k = exp(kappa * kappa - T) * erfc(kappa); - const Real erfc_l = exp(lambda * lambda - T) * erfc(lambda); - - Real Gmm1; // will contain G[m-1] - Real Gm; // will contain G[m] - Real Gmp1; // will contain G[m+1] - Gmm1 = pfac * (erfc_k + erfc_l) * oo_sqrt_U; // G_{-1} - Gm = pfac * (erfc_k - erfc_l) * oo_sqrt_T; // G_{0} + Real Gmm1; // will contain G[m-1] + Real Gm; // will contain G[m] + Real Gmp1; // will contain G[m+1] + std::tie(Gm, Gmm1) = eval_G0_and_maybe_Gm1(T, U); if #if __cplusplus >= 201703L constexpr @@ -1344,8 +1291,194 @@ struct TennoGmEval { Gm = Gmp1; } } + } + + /// Computes core integrals for large \f$ U \f$ by the asymptotic (Watson) + /// expansion in powers of \f$ 1/V \f$ , \f$ V = U - T \f$ . Substituting + /// \f$ t^{-2} = 1 + s \f$ , + /// \f[ + /// G_m(T,U) = \frac{e^{-T}}{2} \int_0^\infty e^{-V s} h_m(s) ds , \quad + /// h_m(s) = (1+s)^{-m-3/2} e^{-T s^2/(1+s)} = \sum_n b_n s^n , + /// \f] + /// hence \f$ G_m = \frac{e^{-T}}{2} \sum_n b_n n!/V^{n+1} \f$ and + /// \f$ G_{m-1} - G_m = \frac{e^{-T}}{2} \sum_n b_n (n+1)!/V^{n+2} \f$ + /// (evaluated directly, without cancellation). The coefficients follow + /// from \f$ (1+s)^2 h_m' = -[a(1+s) + T(2s+s^2)] h_m \f$, \f$ a = m+3/2 \f$: + /// \f$ (n+1) b_{n+1} = -(2n+a) b_n - (n-1+a+2T) b_{n-1} - T b_{n-2} \f$ . + /// For \f$ U > 10^3 \f$, \f$ T \leq 2^{10} \f$ the values are below + /// \f$ 10^{-300} \f$ unless \f$ V > 250 \f$, where the series converges to + /// double precision within a few dozen terms. See issue #444. + /// @tparam[in] Exp if true, will compute core integrals of the exponential + /// operator, otherwise will compute Yukawa + /// @param[out] Gm_vec pointer to an array with at least @c mmax+1 elements, + /// on output contains: + /// - if `Exp==false`, `Gm_vec[m]` contains \f$ G_m(T,U), m=0..mmax\f$ ; + /// - if `Exp==true`, `Gm_vec[m]` contains \f$ \frac{\zeta}{2 \rho} \left( + /// G_{m-1}(T,U) - G_m(T,U) \right), m=0..mmax \f$; + /// @param[in] T the value of \f$ T \f$, must not exceed \c Tmax + /// @param[in] U the value of \f$ U \f$, must exceed \c Umax + /// @param[in] zeta_over_two_rho the value of \f$ \frac{\zeta}{2 \rho} \f$, + /// only used if @c Exp==true + /// @param[in] mmax the maximum value of \f$ m \f$ + template + static void eval_large_U(Real* Gm_vec, Real T, Real U, Real zeta_over_two_rho, + size_t mmax) { + assert(T >= 0 && T <= Tmax); + assert(U > Umax); + // values below are bounded by max(e^{-T}, e^{U - 2 sqrt(U T)}), which + // underflows unless U - T is large + static_assert(cheb_table_umaxlog10 >= 3 && cheb_table_tmaxlog2 <= 10, + "TennoGmEval::eval_large_U(): convergence analysis assumes " + "Umax >= 1e3 and Tmax <= 1024"); + const Real half_exp_mT = 0.5 * exp(-T); + const Real V = U - T; + if (half_exp_mT == 0 || V <= 0) { // e^{-T} underflows, so do the values + std::fill(Gm_vec, Gm_vec + mmax + 1, Real(0)); + return; + } + const Real oo_V = 1 / V; + // the series is asymptotic; for the values that do not underflow + // (V > 255, T < 745) the terms decrease monotonically well beyond nmax + constexpr int nmax = 60; + const Real eps = std::numeric_limits::epsilon() / 4; + + for (int m = 0; m <= static_cast(mmax); ++m) { + const Real a = m + 1.5; + // c_n = b_n n! / V^n ; the n-th term of G_m (V/(e^{-T}/2) times) is c_n, + // of G_{m-1} - G_m (V^2/(e^{-T}/2) times) is (n+1) c_n + Real c_nm2 = 0, c_nm1 = 0, c_n = 1; + Real sum = 1; // the n = 0 term + int nsmall = 0; // number of consecutive negligible terms + for (int n = 0; n < nmax && nsmall < 2; ++n) { + // c_{n+1} = -[(2n+a) c_n/V + (n-1+a+2T) n c_{n-1}/V^2 + // + T n (n-1) c_{n-2}/V^3] + const Real c_np1 = + -oo_V * ((2 * n + a) * c_n + oo_V * n * + ((n - 1 + a + 2 * T) * c_nm1 + + oo_V * T * (n - 1) * c_nm2)); + const Real term = Exp ? (n + 2) * c_np1 : c_np1; + sum += term; + // N.B. require 2 consecutive negligible terms since a coefficient + // can be accidentally small + nsmall = (std::abs(term) <= eps * std::abs(sum)) ? nsmall + 1 : 0; + c_nm2 = c_nm1; + c_nm1 = c_n; + c_n = c_np1; + } + if +#if __cplusplus >= 201703L + constexpr +#endif + (!Exp) { + Gm_vec[m] = half_exp_mT * oo_V * sum; + } else { + Gm_vec[m] = half_exp_mT * oo_V * oo_V * sum * zeta_over_two_rho; + } + } + } - return true; + /// Computes core integrals for small \f$ U \f$ by expanding in powers of + /// \f$ U \f$ about the Boys function: + /// \f[ + /// G_m(T,U) = e^U \left[ \sum_{k=0}^{K} \frac{(-U)^k}{k!} F_{m-k}(T) + + /// \frac{U^{m+1/2}}{2} \sum_{j=0}^{J} \Gamma(-m-j-\tfrac{1}{2}) + /// \frac{(-TU)^j}{j!} \right], + /// \f] + /// where for \f$ n<0 \f$ \f$ F_n(T) \f$ is the Hadamard finite part of the + /// (divergent) Boys integral, obtained from \f$ F_0(T) \f$ by the downward + /// recursion \f$ F_n = (2T F_{n+1} + e^{-T})/(2n+1) \f$. The first sum is + /// the Taylor series of \f$ e^{-U/t^2} \f$ integrated term by term, the + /// second collects the terms nonanalytic in \f$ U \f$ produced by the + /// \f$ t \to 0 \f$ region; the latter are only significant for + /// \f$ m \leq 2 \f$. For \f$ U \leq 10^{-7} \f$ and \f$ T \leq 1024 \f$ the + /// truncation error is below \f$ 10^{-20} \f$ (verified against 50-digit + /// quadrature). Unlike eval_urr() this is stable for small (even zero) T, + /// and unlike interpolate_Gm() it is exact for U = 0. See issue #444. + /// @tparam[in] Exp if true, will compute core integrals of the exponential + /// operator, otherwise will compute Yukawa + /// @param[out] Gm_vec pointer to an array with at least @c mmax+1 elements, + /// on output contains: + /// - if `Exp==false`, `Gm_vec[m]` contains \f$ G_m(T,U), m=0..mmax\f$ ; + /// - if `Exp==true`, `Gm_vec[m]` contains \f$ \frac{\zeta}{2 \rho} \left( + /// G_{m-1}(T,U) - G_m(T,U) \right), m=0..mmax \f$; + /// @param[in] T the value of \f$ T \f$, must not exceed \c Tmax unless + /// \f$ U = 0 \f$ + /// @param[in] U the value of \f$ U \f$, must be smaller than \c Umin + /// (and nonzero if @c Exp==true) + /// @param[in] zeta_over_two_rho the value of \f$ \frac{\zeta}{2 \rho} \f$, + /// only used if @c Exp==true + /// @param[in] mmax the maximum value of \f$ m \f$ + template + void eval_small_U(Real* Gm_vec, Real T, Real U, Real zeta_over_two_rho, + size_t mmax) const { + assert(T >= 0 && (T <= Tmax || U == 0)); + assert(U >= 0 && U < Umin); + assert(!Exp || U > 0); + assert(mmax <= mmax_); + // the truncation orders below assume U < 1e-7 and T <= 1024 + static_assert(cheb_table_uminlog10 <= -7 && cheb_table_tmaxlog2 <= 10, + "TennoGmEval::eval_small_U(): truncation orders assume " + "Umin <= 1e-7 and Tmax <= 1024"); + constexpr int K = 4; // order of the Taylor series in U + constexpr int J = 4; // order of the series in TU of the nonanalytic + // part + constexpr int mmin = -1; // need G_{-1} if Exp + constexpr int mmax_nonanalytic = 2; // U^{m+1/2} is negligible for m > 2 + constexpr int nmin = mmin - K; // lowest-order F_n needed + const Real sqrtPi(1.7724538509055160272981674833411451827975494561224); + + // F0[n] = F_n(T), n = nmin .. mmax + Real F[static_cast(cheb_table_mmax) + 1 - nmin]; + Real* const F0 = F - nmin; + fm_eval_->eval(F0, T, static_cast(mmax)); + const Real exp_mT = exp(-T); + const Real two_T = 2 * T; + for (int n = -1; n >= nmin; --n) + F0[n] = (two_T * F0[n + 1] + exp_mT) / (2 * n + 1); + + // (-U)^k / k! + Real Uk_over_kfac[K + 1]; + Uk_over_kfac[0] = 1; + for (int k = 1; k <= K; ++k) + Uk_over_kfac[k] = Uk_over_kfac[k - 1] * (-U) / k; + + // gamma_mnmh[n - mmin] = Gamma(-n-1/2), n = mmin .. mmax_nonanalytic + J + constexpr int ngamma = mmax_nonanalytic + J - mmin + 1; + Real gamma_mnmh[ngamma]; + gamma_mnmh[0] = sqrtPi; // Gamma(1/2) + for (int n = mmin + 1; n <= mmax_nonanalytic + J; ++n) + gamma_mnmh[n - mmin] = gamma_mnmh[n - mmin - 1] * (-2) / (2 * n + 1); + + const Real exp_U = exp(U); + const Real minus_TU = -T * U; + Real U_pow_m_half = Exp ? 1 / sqrt(U) : sqrt(U); // U^{m+1/2} + + Real Gmm1 = 0; // will contain G[m-1] + for (int m = Exp ? mmin : 0; m <= static_cast(mmax); ++m) { + Real Gm = 0; + for (int k = 0; k <= K; ++k) Gm += Uk_over_kfac[k] * F0[m - k]; + if (m <= mmax_nonanalytic) { + Real sum = 0; + Real TUj_over_jfac = 1; // (-TU)^j / j! + for (int j = 0; j <= J; ++j) { + sum += gamma_mnmh[m + j - mmin] * TUj_over_jfac; + TUj_over_jfac *= minus_TU / (j + 1); + } + Gm += 0.5 * U_pow_m_half * sum; + U_pow_m_half *= U; + } + Gm *= exp_U; + if +#if __cplusplus >= 201703L + constexpr +#endif + (!Exp) { + Gm_vec[m] = Gm; + } else { + if (m != mmin) Gm_vec[m] = (Gmm1 - Gm) * zeta_over_two_rho; + Gmm1 = Gm; + } + } } /// compute core integrals by Chebyshev interpolation @@ -1359,26 +1492,12 @@ struct TennoGmEval { /// @param[in] zeta_over_two_rho the value of \f$ \frac{\zeta}{2 \rho} \f$, /// only used if @c Exp==true /// @param[in] mmax the maximum value of \f$ m \f$ - /// @throw std::invalid_argument if \p U exceeds \c Umax - /// @return true if {T,U} can be handled without overflow; false return - /// indicates that Gaussian fit of the kernel should be used to - /// compute Gm + /// @pre \f$ T \leq \f$ \c Tmax and \c Umin \f$ \leq U \leq \f$ \c Umax template - bool interpolate_Gm(Real* Gm_vec, Real T, Real U, Real zeta_over_two_rho, + void interpolate_Gm(Real* Gm_vec, Real T, Real U, Real zeta_over_two_rho, long mmax) const { - assert(T >= 0); - if (U > Umax) { - libint2::verbose_stream() - << "TennoGmEval::eval_{yukawa,slater}() : arguments out of range, " - << "zeta*zeta*one_over_rho/4=" << std::to_string(U) - << " cannot exceed " << std::to_string(Umax) - << " will represent Slater function as a linear combination of" - << " Gaussian functions " << std::endl; - return false; - } - if (U < Umin || U > Umax) { - return false; - } + assert(T >= 0 && T <= Tmax); + assert(U >= Umin && U <= Umax); // maps x in [0,2] to [-1/2,1/2] in a linear fashion auto linear_map_02 = [](Real x) { @@ -1476,11 +1595,7 @@ struct TennoGmEval { // precision of interpolation for m=-1,0 can be insufficient, just // evaluate explicitly Real G0; - const auto eval_G0_Gm1_result = eval_G0_and_maybe_Gm1(T, U); - if (!eval_G0_Gm1_result.has_value()) { - return false; - } - std::tie(Exp ? G0 : Gm_vec[0], Gmm1) = eval_G0_Gm1_result.value(); + std::tie(Exp ? G0 : Gm_vec[0], Gmm1) = eval_G0_and_maybe_Gm1(T, U); if #if __cplusplus >= 201703L constexpr @@ -1661,8 +1776,6 @@ struct TennoGmEval { Gmm1 = Gm; } } - - return true; } private: @@ -1670,6 +1783,8 @@ struct TennoGmEval { Real precision_; ExpensiveNumbers numbers_; Real* c_ = nullptr; + std::shared_ptr> + fm_eval_; // Boys function evaluator used by eval_small_U() void init_table() { // get memory