You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
I noticed that rwfnestimate and rmcdhf was producing a lot of NaNs when I was trying a simple helium calculation. I traced it back to rnucleus not handling finite nuclei very well for small masses, specifically, it is not able to reliably find a lower limit for the bracketing algorithm used to find the value of CPARM which is written to isodata.
Build steps
To aid in my debugging, I enabled a debug flag (-ffpe-trap=invalid) that halts execution when NaNs result from a calculation (so called signaling NaNs). I additionally added some printouts that let me follow the progress of the iterative algorithm:
diff --git a/configure.sh b/configure.sh
index 44d43309..062010e0 100755
--- a/configure.sh+++ b/configure.sh@@ -35,7 +35,7 @@ fi
# cmake ..
#
mkdir "${build_abspath}" && cd "${build_abspath}" \
- && cmake ${cmake_args} "${GRASP}" \+ && cmake ${cmake_args} -DCMAKE_Fortran_FLAGS="-fallow-argument-mismatch" -DCMAKE_Fortran_FLAGS_DEBUG="-g -O0 -fcheck=all -fbacktrace -ffpe-trap=invalid" "${GRASP}" \
|| exit
# Note: we need to use spaces, not tabs, to indent in the heredoc.
diff --git a/src/appl/rnucleus90/estrms.f90 b/src/appl/rnucleus90/estrms.f90
index 75e1eb77..c3d7cab7 100644
--- a/src/appl/rnucleus90/estrms.f90+++ b/src/appl/rnucleus90/estrms.f90@@ -44,6 +44,7 @@
DNUMER = 1.0D00 + (10.0D00/3.0D00)*PABC**2 + (7.0D00/3.0D00)*PABC**4 - &
120.0D00*ABC**5*SKFUN(5,CBAM)
DDENOM = 1.0D00 + PABC**2 - 6.0D00*ABC**3*SKFUN(3,CBAM)
+ write(*,*) "CPARM:",CPARM,"SQTBF:",SQTBF,"DNUMER:",DNUMER,"DDENOM:",DDENOM
ESTRMS = CPARM*SQTBF*SQRT(DNUMER/DDENOM)
!
RETURN
diff --git a/src/lib/lib9290/nucpot.f90 b/src/lib/lib9290/nucpot.f90
index 0c27fc07..4f772499 100644
--- a/src/lib/lib9290/nucpot.f90+++ b/src/lib/lib9290/nucpot.f90@@ -62,7 +62,9 @@
SABC3 = 6.0D00*ABC3
DMSAS = -SABC3*S3MCBA
EN = 1.0D00 + ABC2*PI2 + DMSAS
+ write(*,*) "ABC2:", ABC2, "PI2:", PI2, "DMSAS:", DMSAS
ZBN = Z/EN
+ write(*,*) "EN: ", EN, "Z: ", Z, "ZBN: ", ZBN
!
SET = .FALSE.
DO I = 1, N
@@ -72,6 +74,8 @@
RBC = RI/C
IF (RBC <= 1.0D00) THEN
CALL ES (RMCBA, S2RCBA, S3RCBA)
+ write (*,*) ZBN, DMSAS, SABC3, S3RCBA, RBC, H3PHP, &+ THABC2, S2RCBA, RBC
ZZ(I) = ZBN*(DMSAS + SABC3*S3RCBA + RBC*(H3PHP - THABC2*S2RCBA&
- 0.5D00*RBC*RBC))
ELSE
Then, I built Grasp using ./configure --debug; cd build-debug; make install.
Example calculations
✅ He point charge
#!/bin/bash
rnucleus <<EOF2 ! Atomic number0 ! Mass number0 ! Mass of nucleus0 ! Nuclear spin (I) (in units of h / 2 pi)0 ! Nuclear dipole moment (in nuclear magnetons)0 ! Nuclear quadrupole moment (in barns)EOF
cat isodata
Output:
Enter the atomic number:
Enter the mass number (0 if the nucleus is to be modelled as a point source:
Enter the mass of the neutral atom (in amu) (0 if the nucleus is to be static):
Enter the nuclear spin quantum number (I) (in units of h / 2 pi):
Enter the nuclear dipole moment (in nuclear magnetons):
Enter the nuclear quadrupole moment (in barns):
Atomic number:
2.0000000000000000
Mass number (integer) :
0.0000000000000000
Fermi distribution parameter a:
0.0000000000000000
Fermi distribution parameter c:
0.0000000000000000
Mass of nucleus (in amu):
0.0000000000000000
Nuclear spin (I) (in units of h / 2 pi):
0.0000000000000000
Nuclear dipole moment (in nuclear magnetons):
0.0000000000000000
Nuclear quadrupole moment (in barns):
0.0000000000000000
❌ He finite nucleus
#!/bin/bash
rnucleus <<EOF2 ! Atomic number4 ! Mass numbern ! Don't revise anything4.00260325413 ! Mass of nucleus0 ! Nuclear spin (I) (in units of h / 2 pi)0 ! Nuclear dipole moment (in nuclear magnetons)0 ! Nuclear quadrupole moment (in barns)n ! Don't revise anythingEOF
cat isodata
Here, we see the NaN occurs because we are trying to divide Infinity by Infinity. Running through a debugger, we see this happens when trying to find a lower limit for CPARM for the bracketing algorithm:
* thread #1, queue = 'com.apple.main-thread', stop reason = EXC_BAD_INSTRUCTION (code=1, subcode=0x1e601840)
frame #0: 0x00000001000050f4 rnucleus`estrms_ at estrms.f90:48:46
45 120.0D00*ABC**5*SKFUN(5,CBAM)
46 DDENOM = 1.0D00 + PABC**2 - 6.0D00*ABC**3*SKFUN(3,CBAM)
47 write(*,*) "CPARM:",CPARM,"SQTBF:",SQTBF,"DNUMER:",DNUMER,"DDENOM:",DDENOM
-> 48 ESTRMS = CPARM*SQTBF*SQRT(DNUMER/DDENOM)
49 !
50 RETURN
51 END FUNCTION ESTRMS
Target 0: (rnucleus) stopped.
warning: This version of LLDB has no plugin for the language "fortran90". Inspection of frame variables will be limited.
(lldb) bt
* thread #1, queue = 'com.apple.main-thread', stop reason = EXC_BAD_INSTRUCTION (code=1, subcode=0x1e601840)
* frame #0: 0x00000001000050f4 rnucleus`estrms_ at estrms.f90:48:46
frame #1: 0x00000001000051c8 rnucleus`getcpr_ at getcpr.f90:50:42
frame #2: 0x0000000100005e04 rnucleus`MAIN__ at geniso.f90:67:10
frame #3: 0x0000000100006618 rnucleus`main at geniso.f90:13:10
frame #4: 0x000000018abc20e0 dyld`start + 2360
✅ Xe finite nucleus
rnucleus <<EOF54 ! Atomic number131 ! Mass numbern ! Don't revise anything131.170428596076 ! Mass of nucleus0 ! Nuclear spin (I) (in units of h / 2 pi)0 ! Nuclear dipole moment (in nuclear magnetons)0 ! Nuclear quadrupole moment (in barns)n ! Don't revise anythingEOF
cat isodata
Output:
Enter the atomic number:
Enter the mass number (0 if the nucleus is to be modelled as a point source:
The default root mean squared radius is 4.7807998657226562 fm; (Angeli)
the default nuclear skin thickness is 2.2999999999999998 fm;
Revise these values?
CPARM: 2.3903999328613281 SQTBF: 0.77459666924148340 DNUMER: 3.1001987503262654 DDENOM: 1.4738104922203277
CPARM: 9.5615997314453125 SQTBF: 0.77459666924148340 DNUMER: 1.1006149549187731 DDENOM: 1.0295723209700669
CPARM: 5.9759998321533203 SQTBF: 0.77459666924148340 DNUMER: 1.2657234387500615 DDENOM: 1.0757051859788953
CPARM: 4.1831998825073242 SQTBF: 0.77459666924148340 DNUMER: 1.5706996655253853 DDENOM: 1.1545042606209786
CPARM: 5.0795998573303223 SQTBF: 0.77459666924148340 DNUMER: 1.3748924850426643 DDENOM: 1.1047826031094967
CPARM: 5.5277998447418213 SQTBF: 0.77459666924148340 DNUMER: 1.3131979032671883 DDENOM: 1.0884794793050161
CPARM: 5.7518998384475708 SQTBF: 0.77459666924148340 DNUMER: 1.2879792662424656 DDENOM: 1.0817192374492073
CPARM: 5.6398498415946960 SQTBF: 0.77459666924148340 DNUMER: 1.3001862243622659 DDENOM: 1.0849986365390318
CPARM: 5.5838248431682587 SQTBF: 0.77459666924148340 DNUMER: 1.3065870868286245 DDENOM: 1.0867128643059558
CPARM: 5.6118373423814774 SQTBF: 0.77459666924148340 DNUMER: 1.3033609784070448 DDENOM: 1.0858493326169478
CPARM: 5.6258435919880867 SQTBF: 0.77459666924148340 DNUMER: 1.3017672512855101 DDENOM: 1.0854223960961389
CPARM: 5.6328467167913914 SQTBF: 0.77459666924148340 DNUMER: 1.3009751588461611 DDENOM: 1.0852101211716068
CPARM: 5.6363482791930437 SQTBF: 0.77459666924148340 DNUMER: 1.3005802979218253 DDENOM: 1.0851042803142934
CPARM: 5.6380990603938699 SQTBF: 0.77459666924148340 DNUMER: 1.3003831628538103 DDENOM: 1.0850514338220065
CPARM: 5.6372236697934568 SQTBF: 0.77459666924148340 DNUMER: 1.3004817057992313 DDENOM: 1.0850778509131647
CPARM: 5.6376613650936633 SQTBF: 0.77459666924148340 DNUMER: 1.3004324281814412 DDENOM: 1.0850646408293170
CPARM: 5.6378802127437666 SQTBF: 0.77459666924148340 DNUMER: 1.3004077939816143 DDENOM: 1.0850580369411544
CPARM: 5.6377707889187150 SQTBF: 0.77459666924148340 DNUMER: 1.3004201106974924 DDENOM: 1.0850613387891015
CPARM: 5.6378255008312408 SQTBF: 0.77459666924148340 DNUMER: 1.3004139522435485 DDENOM: 1.0850596878410954
CPARM: 5.6377981448749779 SQTBF: 0.77459666924148340 DNUMER: 1.3004170314465189 DDENOM: 1.0850605133090900
CPARM: 5.6378118228531093 SQTBF: 0.77459666924148340 DNUMER: 1.3004154918390334 DDENOM: 1.0850601005735907
CPARM: 5.6378186618421751 SQTBF: 0.77459666924148340 DNUMER: 1.3004147220397910 DDENOM: 1.0850598942069674
CPARM: 5.6378152423476422 SQTBF: 0.77459666924148340 DNUMER: 1.3004151069390371 DDENOM: 1.0850599973901853
CPARM: 5.6378135326003758 SQTBF: 0.77459666924148340 DNUMER: 1.3004152993889417 DDENOM: 1.0850600489818645
CPARM: 5.6378143874740090 SQTBF: 0.77459666924148340 DNUMER: 1.3004152031639657 DDENOM: 1.0850600231860190
CPARM: 5.6378148149108256 SQTBF: 0.77459666924148340 DNUMER: 1.3004151550514957 DDENOM: 1.0850600102881005
CPARM: 5.6378146011924173 SQTBF: 0.77459666924148340 DNUMER: 1.3004151791077294 DDENOM: 1.0850600167370594
CPARM: 5.6378147080516214 SQTBF: 0.77459666924148340 DNUMER: 1.3004151670796122 DDENOM: 1.0850600135125799
CPARM: 5.6378147614812235 SQTBF: 0.77459666924148340 DNUMER: 1.3004151610655537 DDENOM: 1.0850600119003402
CPARM: 5.6378147881960246 SQTBF: 0.77459666924148340 DNUMER: 1.3004151580585246 DDENOM: 1.0850600110942203
CPARM: 5.6378148015534251 SQTBF: 0.77459666924148340 DNUMER: 1.3004151565550102 DDENOM: 1.0850600106911605
CPARM: 5.6378148082321253 SQTBF: 0.77459666924148340 DNUMER: 1.3004151558032528 DDENOM: 1.0850600104896306
CPARM: 5.6378148115714755 SQTBF: 0.77459666924148340 DNUMER: 1.3004151554273744 DDENOM: 1.0850600103888657
CPARM: 5.6378148132411505 SQTBF: 0.77459666924148340 DNUMER: 1.3004151552394347 DDENOM: 1.0850600103384831
CPARM: 5.6378148124063134 SQTBF: 0.77459666924148340 DNUMER: 1.3004151553334045 DDENOM: 1.0850600103636743
CPARM: 5.6378148119888944 SQTBF: 0.77459666924148340 DNUMER: 1.3004151553803895 DDENOM: 1.0850600103762698
CPARM: 5.6378148117801850 SQTBF: 0.77459666924148340 DNUMER: 1.3004151554038816 DDENOM: 1.0850600103825678
CPARM: 5.6378148116758302 SQTBF: 0.77459666924148340 DNUMER: 1.3004151554156280 DDENOM: 1.0850600103857166
CPARM: 5.6378148117280080 SQTBF: 0.77459666924148340 DNUMER: 1.3004151554097549 DDENOM: 1.0850600103841421
CPARM: 5.6378148117019187 SQTBF: 0.77459666924148340 DNUMER: 1.3004151554126915 DDENOM: 1.0850600103849295
CPARM: 5.6378148116888749 SQTBF: 0.77459666924148340 DNUMER: 1.3004151554141599 DDENOM: 1.0850600103853232
Enter the mass of the neutral atom (in amu) (0 if the nucleus is to be static):
Enter the nuclear spin quantum number (I) (in units of h / 2 pi):
Enter the nuclear dipole moment (in nuclear magnetons):
Enter the nuclear quadrupole moment (in barns):
Atomic number:
54.000000000000000
Mass number (integer) :
131.00000000000000
Fermi distribution parameter a:
0.52338755531043146
Fermi distribution parameter c:
5.6378148116888749
Mass of nucleus (in amu):
131.14080528131399
Nuclear spin (I) (in units of h / 2 pi):
0.0000000000000000
Nuclear dipole moment (in nuclear magnetons):
0.0000000000000000
Nuclear quadrupole moment (in barns):
0.0000000000000000
The text was updated successfully, but these errors were encountered:
Summary
I noticed that
rwfnestimateandrmcdhfwas producing a lot of NaNs when I was trying a simple helium calculation. I traced it back tornucleusnot handling finite nuclei very well for small masses, specifically, it is not able to reliably find a lower limit for the bracketing algorithm used to find the value ofCPARMwhich is written toisodata.Build steps
To aid in my debugging, I enabled a debug flag (
-ffpe-trap=invalid) that halts execution when NaNs result from a calculation (so called signaling NaNs). I additionally added some printouts that let me follow the progress of the iterative algorithm:Then, I built Grasp using
./configure --debug; cd build-debug; make install.Example calculations
✅ He point charge
Output:
❌ He finite nucleus
Output:
Here, we see the NaN occurs because we are trying to divide Infinity by Infinity. Running through a debugger, we see this happens when trying to find a lower limit for
CPARMfor the bracketing algorithm:✅ Xe finite nucleus
Output:
The text was updated successfully, but these errors were encountered: