Fortran wrapper

[KBGrammer]

I'm working on a fortran wrapper so that I can call NCrystal in MCNP. I would like to be able to call ncrystal_get_flatcompos so that I can get the composition with isotope fractions so that I can make the MCNP composition of a material match that of the NCrystal material without doing so by hand. However, I haven't figured out a way to poll the NCrystal AtomDB for the element fractions, so I haven't found the function pointer to give it for natelemprovider_raw.

Is there an internal function that I can give it for that parameter?

I've written one that works based on Nubase 2020, which works alright. I was wondering if there's a built in one, though.

[tkittel]

That is a very good question indeed! The short answer is that currently, NCrystal does in fact not provide such numbers. Currently, the ncrystal_get_flatcompos functionality has been needed in the context of Geant4 and OpenMC hooks, and in both of those cases we ended up deciding that it might be better to use the abundance fractions already provided by Geant4 or OpenMC itself (hence the decision to let ncrystal_get_flatcompos accept that data via a function pointer). That way, we won't end up assigning a non-zero fraction for some ultra-rare/obscure isotope for which the application itself (Geant4/OpenMC) might not provide support (i.e. some datafiles might be missing). However, the argument against this is of course that one would get a more uniform experience when using NCrystal as backend in different applications, if NCrystal provided its own database of abundances.

However, I think it is reasonable to suggest that at the very least NCrystal should provide such numbers, so that people with use-cases like your own might make their own choice (+ we might need such data anyway eventually for our internal modelling). In fact, if you look at the file where NCrystal's (currently hardcoded) database of atomic data, you can see that we could actually easily provide these numbers (there would mainly be some API coding work to do):

https://github.com/mctools/ncrystal/blob/v3.7.1/ncrystal_core/src/NCAtomDB.cc#L245

In particular note the comment in lines 271-272, and the commented-out abundance values below (e.g. in line 277 for H1).

Based on this I am happy to open a feature request to make NCrystal also provide this information, and to let the ncrystal_get_flatcompos function default to using those.

Btw., note that the abundances are only actually strictly needed/used if the prefernatelem parameter is 0, or if you are dealing with a rather exotic material where somehow a particular element appears both as a natural element ("O") and as a particular isotope ("O16").

The next NCrystal release might be a while off still though, since I am currently deeply engrossed in some other work, so I hope it is correctly understood that you have a work-around with the Nubase 2020 DB for now?

Finally, it sounds like an interesting wrapper you are working on. If you get it to work, it might be interesting to share more widely. Well, as widely as one can with MCNP-adjacent code... the MCNP developers themselves might at least be interested I think.

Cheers,
Thomas

[tkittel]

Great, this will be very useful indeed!

For reference/curiosity I wanted to double-check that your result above would be compatible with the fractions I got if I fed the same cfg-string into our Geant4 material converter. So I did that and got the printout below:

===================================================================================
Material defined by string "NCrystal::cfg=[gasmix::0.72xCO2+0.28xAr/massfractions/3bar/250K]":

 Material: NCrystalRel::gasmix::0.72xCO2+0.28xAr/massfractions/3bar/250K    density:  6.176 mg/cm3  RadL:  47.323 m    Nucl.Int.Length: 150.937 m  
                       Imean: 108.986 eV   temperature: 250.00 K  pressure:   1.00 atm

   --->  Element: C (C)   Z =  6.0   N =    12   A = 12.011 g/mole
         --->  Isotope:   C12   Z =  6   N =  12   A =  12.00 g/mole   abundance: 98.930 %
         --->  Isotope:   C13   Z =  6   N =  13   A =  13.00 g/mole   abundance:  1.070 %
          ElmMassFraction:  19.65 %  ElmAbundance  29.17 % 

   --->  Element: O (O)   Z =  8.0   N =    16   A = 15.999 g/mole
         --->  Isotope:   O16   Z =  8   N =  16   A =  15.99 g/mole   abundance: 99.757 %
         --->  Isotope:   O17   Z =  8   N =  17   A =  17.00 g/mole   abundance:  0.038 %
         --->  Isotope:   O18   Z =  8   N =  18   A =  18.00 g/mole   abundance:  0.205 %
          ElmMassFraction:  52.35 %  ElmAbundance  58.34 % 

   --->  Element: Ar (Ar)   Z = 18.0   N =    40   A = 39.948 g/mole
         --->  Isotope:  Ar36   Z = 18   N =  36   A =  35.97 g/mole   abundance:  0.337 %
         --->  Isotope:  Ar38   Z = 18   N =  38   A =  37.96 g/mole   abundance:  0.063 %
         --->  Isotope:  Ar40   Z = 18   N =  40   A =  39.96 g/mole   abundance: 99.600 %
          ElmMassFraction:  28.00 %  ElmAbundance  12.50 % 

===================================================================================

Which if I manually massage it into a format like in your listing, is:

6   12 0.28857881
6   13 0.00312119
8   16 0.581982338
8   17 0.000221692
8   18 0.00119597
18  36 0.00042125
18  38 7.875e-05
18  40 0.1245

Or, with relative differences (a.k.a. abs(x/xref-1)) to your output above calculated as:

 6   12 0.28857881      # relative difference wrt KBGrammer impl of 2.74639e-05
 6   13 0.00312119      # relative difference wrt KBGrammer impl of 0.00950829
 8   16 0.581982338     # relative difference wrt KBGrammer impl of 7.21467e-05
 8   17 0.000221692     # relative difference wrt KBGrammer impl of 0.00905495
 8   18 0.00119597      # relative difference wrt KBGrammer impl of 0.0025173
18   36 0.00042125      # relative difference wrt KBGrammer impl of 0.0104792
18   38 7.875e-05       # relative difference wrt KBGrammer impl of 0.00187471
18   40 0.1245          # relative difference wrt KBGrammer impl of 0.000249261

So I think it is safe to conclude that this is basically consistent :-)

[marquezj]

Hi Kyle. This is a bit tough because it should reflect whatever library is used for the high energy part. Maybe the wrapper could take a parameter for the library used ("endf70", "endf71", "endf80") and provide the abundances that correspond. For instance, for ENDF/B-VII.1 and earlier, carbon is natural. This is probably very close to NUBASE, but not exactly it.

Also, probably the wrapper would have to map the temperature to whatever extension the MCNP library has for it. Or, take as a parameter which suffix should be used for all the nuclides.

[tkittel]

@marquezj I would argue that it is not really a given that you absolutely need to reflect the library used for the high energy part. Conversely, one might equally well argue that the low energy part is the most important, and that the abundances used there (in the name of reproducibility of NCrystal physics) should be the same in all applications where NCrystal is integrated - and therefore that the high energy should just accept whatever NCrystal decides.

I can see arguments either way, IMO it is absolutely not a clear cut case. However, the only reason it is important at all is of course that there is a significant variation across codes and databases. In the example above though, the variations are all so small as to be dwarfed by other effects (e.g. the introduction of sampling artifacts, etc.).

[KBGrammer]

The wrapper I've written does take a temperature as an option, and it can be given to the MCNP material. I don't think it's right to go fishing through the endf-vii cross section library for 600C, though.

I'm not sure that saying "endf70" to get endf-vii is the right way because the library names could be all over the place rather than just the stock ones, in principle.

I had thought about carbon specifically, given that it's 6000.XX in earlier libraries. I wonder if the best way would be to have the user specify the isotopic composition and libraries that they want as normal as an MCNP material and then just update the fractions from the NCrystal information. In that case, the wrapper could look at the MCNP material, see that 6000 is given, and add the fractions for C-12 and C-13 (or pull a natural and isotopic composition from NCrystal and switch between them). It requires the user to pay closer attention to the material definition, and that's not a bad thing.

Take the CO2+Ar mix above. If the user gives the material as the following:

m1    6000.70c 1
      8016.70c 1
      8017.70c 1
      8018.70c 1
      18036.70c 1
      18038.70c 1
mn1   gasmix::0.72xCO2+0.28xAr/massfractions/3bar/250K;

The wrapper could throw a warning because 18040 is missing and the NCrystal material wants it to be present, and the fractions for C-12 and C-13 could be added together for C-natural because natural carbon is requested. It could also throw a warning if the temperature doesn't match.

Doing it with warnings also allows MCNP to read the libraries it wants to read at the time it wants to read them without having to read NCrystal earlier to update the material card. My implementation is going to work like a thermal treatment card, where you would say the following:

m1    6000.70c 1
      8016.70c 1
      8017.70c 1
      8018.70c 1
      18036.70c 1
      18038.70c 1
      18040.70c 1
mn1   gasmix::0.72xCO2+0.28xAr/massfractions/3bar/250K;

Using warnings also allows users to specify a material where the material card is something like aluminum 6061 with alloying components in it but they want to apply just an Al-27 powder treatment to it at low energies. It's not particularly rigorous, but throwing a warning and letting the user do the strange thing they want to do is probably better than a fatal error. Maybe a better example than 6061 is a material where you care about the NCrystal treatment at low energies but something that doesn't matter much at low energies matters at high energies.

I also could see validity in letting the Al-27 powder treatment only affect 13027.70c in a material and letting the alloying components just use their ENDF cross sections, which is how MCNP does it with a thermal treatment. However, maybe it's much easier to make an NCrystal material that has all of that than it is to make a new ENDF library with all of that.

[KBGrammer]

Looking closer at how to incorporate the wrapper, I think the idea that NCrystal could control all interactions below a certain energy isn't going to work. The mechanism for absorption in the analogous McStas component is to either terminate the track or reduce the weight, but this would need to create a particle on absorption in MCNP.

I suspect that I'd need to handle 1 NCrystal material as 2 separate materials with a scattering material and an absorption material. It looks like attaching "coh_elas=0;incoh_elas=0;sans=0;" will give me absorption only and "inelas=0;" would give me scattering only. I think that this means NCrystal could handle the scattering part and MCNP routines and data tables would then need to handle the absorption processes.

[marquezj]

Hi Kyle. Take a look at how we implemented it in OpenMC. The logic is very similar to MCNP. NCrystal only replaces elastic scattering below 5 eV. Above that energy and all other reactions (including capture) are handled by the regular processes.

…

On Fri, Sep 29, 2023, 18:04 Kyle Grammer ***@***.***> wrote: Looking closer at how to incorporate the wrapper, I think the idea that NCrystal could control all interactions below a certain energy isn't going to work. The mechanism for absorption in the analogous McStas component is to either terminate the track or reduce the weight, but this would need to create a particle on absorption in MCNP. I suspect that I'd need to handle 1 NCrystal material as 2 separate materials with a scattering material and an absorption material. It looks like attaching "coh_elas=0;incoh_elas=0;sans=0;" will give me absorption only and "inelas=0;" would give me scattering only. I think that this means NCrystal could handle the scattering part and MCNP routines and data tables would then need to handle the absorption processes. — Reply to this email directly, view it on GitHub <#144 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AFLRDEM4RYXMRZ5BE7G7JG3X43WRDANCNFSM6AAAAAA5GKFMQM> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[tkittel]

FYI the changes needed for OpenMC was in this PR: openmc-dev/openmc#2222

As Ignacio says, NCrystal only modifies the scattering below 5eV. Absorption at any energy is handled by whatever was already the in OpenMC (similarly to what we do for Geant4). Technically, there is nothing special needed in the NCrystal cfg-string in order to separate out scattering processes from absorption processes: By calling NCrystal::createScatter you get the scattering processes, and by calling NCrystal::createAbsorption you get the absorption processes. Now, in general you have the right idea IMO: The material setup should be automated by asking NCrystal for the (flattened) composition. Because ideally (that's my guiding principle at least), users should simply be able to supply the same cfg-string when working with any application (MCNP/OpenMC/McStas/Geant4/...) and get the same physics (at least scattering physics below 5eV).

[KBGrammer]

A keyword "arbitrary" that just echoes everything that follows to the config string should handle that with no problems unless an issue crops up due to string length limits. I don't know how it will behave when such a string crosses input lines, for example.

[marquezj]

Well, maybe multiline reading can be handled with a terminator:

phases<0.9*Al_sg225.ncmat&
0.1*solid::U/19gcm3/U_is_0.20_U235_0.80_U238>;
density=6.1gcm3
#TERMINATOR

The terminator should be a combination of symbols that is not part of neither NCrystal or MCNP syntax. Or it could be wrapped between tags:

BEGIN_NCRYSTAL_CFG
phases<0.9*Al_sg225.ncmat&
0.1*solid::U/19gcm3/U_is_0.20_U235_0.80_U238>;
density=6.1gcm3
END_NCRYSTAL_CFG

And then the liens are concatenated before being sent to NCrystal.

[KBGrammer]

It turns out that MCNP has syntax that it ignores or uses as a delimter that NCrystal uses for important syntax, for instance ampersand and colon. This can work if the MCNP input line uses latex style commands, so "\amp" instead of "&" and "\colon" instead of ":", for instance. These commands get replaced with the character of choice. It will also happily break up such config string across input lines, so it can be arbitrarily long without issue. Take the following from the examples above:

mn15    arbitrary phases\langle0.9*Al_sg225.ncmat\amp0.1*solid\colon\colon
        U/19gcm3/U_is_0.20_U235_0.80_U238\rangle;
        density=6.1gcm3

which produces the following config string for NCrystal:
phases<0.9*Al_sg225.ncmat&0.1*solid::U/19gcm3/U_is_0.20_U235_0.80_U238>;density=6.1gcm3;

Are there other such special characters that might appear in an NCrystal config string besides equals, ampersand, colon, <, and >? Slashes are handled without issue. In particular, does NCrystal ever use parentheses, square brackets, commas, $, or #?

[tkittel]

Wow, that seems troublesome indeed! :-)

Is the support for $\LaTeX$ commands like \amp for &, etc. a feature you have implemented yourself, or a general MCNP feature? If the former, perhaps it would be better to find a more standard encoding that could be supported in tool-chains (e.g. in python code). What about https://en.wikipedia.org/wiki/Percent-encoding ?

As for your question: Yes we use commas (,) - see for instance the dir1 parameter in https://github.com/mctools/ncrystal/wiki/CfgRefDoc). In general, only "'|(){}[] are explicitly forbidden currently (see the list in this line: https://github.com/mctools/ncrystal/blob/master/ncrystal_core/include/NCrystal/internal/NCCfgTypes.hh#L43), although I could probably easily be persuaded to add $ and # to that list (and they are probably not used in the wild right now). (UPDATE: I have opened #153 for that),

Also be aware, that while in general only ASCII characters are supported in NCrystal cfg-strings, the exception is the filename which can be in UTF8 (because in theory an internaltional user might have a non-ASCII character in a folder name, and use an absolute path to indicate that file). Of course, for C/C++ code that is no trouble, since UTF8 has no null-bytes (which is part of the reason UTF8 is great). You could of course make the reasonable restriction of only supporting pure ASCII NCrystal-cfg strings, if this somehow gives problems.

[KBGrammer]

This was something I came up with. Whatever is used has to be passed down from MCNP, and delimiters get stripped. Percent encoding should work just fine.

I believe it's going to have to be ASCII only. I don't think I'll be changing how MCNP interprets text.

I'll take into account that list of forbidden characters. Most of them are delimiters in MCNP anyway. It's just colon and some of the multiphase characters that need some other treatment.