flexible-atomic-code / fac Goto Github PK

Hello,
at the moment i am trying to install PFAC.
The command:
make
is working properly. It is also possible to install SFAC, but in the case of PFAC the command
make pfac
is not working. My C-compiler is gcc and the error message is

gcc: error: : No such file or directory
gcc: error: : No such file or directory
error: command 'gcc' failed with exit status 1
Makefile:42: recipe for target 'pfac' failed
make: *** [pfac] Error 1

For the command
make install-pfac
i get the error message

python/fac.c:23:10: fatal error: init.h: No such file or directory
#include "init.h"
^~~~~~~~
compilation terminated.
error: command 'gcc' failed with exit status 1
Makefile:56: recipe for target 'install-pfac' failed
make: *** [install-pfac] Error 1

Is there a possibility to fix this? Is there somewhere a list of packages that are required for installing the programm or is it a problem of the setup.py file?

Best regards

Dear all:
When I was using the CR Model of FAC, I found some problems. I want to call the AddIon function multiple times to make more than one ion in the model.And i know the order must be such that the number of electrons are in consecutive increasing order.But when I use it, there is an error like this:ERROR: Ion charge state does not match 0 0 11 10.I don't know how to solve it, hope someone can help me.
And my script is:
z = 13
nele = 11
nele1 = 10
temp = 10
dens = 1.0
# a is the atomic symbol corresponding to z
a = fac.ATOMICSYMBOL[z]
p = '%s%02d'%(a, nele)
add the ion with k number of electrons
AddIon(nele1, 0.0, '%sb'%p)
AddIon(nele, 0.0, '%sb'%p)
SetBlocks(-1)
SetEleDist(0, temp, -1, -1)

Hi,

I am using FAC to calculate the energy level. When l is 20, it works well. However, once l is more than 20, for example, l = 21, there is always the error 'Segmentation fault (core dumped)'. Is there any method to avoid such error?

The script is:
fac.SetAtom('Xe')
fac.Config('1s', group = 'i')
fac.Config('1s1 24[20]1', group = 'd') #This works. But fac.Config('1s1 24[21]1', group = 'd') will make mistake.

fac.SetBreit(-1)
fac.ConfigEnergy(0)
fac.OptimizeRadial(['i'])
fac.ConfigEnergy(1)
fac.Structure('Xe.lev.b', [ 'd'])
fac.MemENTable('Xe.lev.b')
fac.PrintTable('Xe.lev.b', 'Xe.lev', 1)

The files necessary to reproduce the issue can be found at: https://mcloud.mpi-hd.mpg.de/index.php/s/BYgbNNrxPskriCt

Thanks to the previous bug fix, I was able to get a CE file for the magnetic sublevels that I'm interested in. I load this using the pfac.pol functions and try to run the population iteration, no errors are given. However, in the populations file, most excitations above 15 eV appear to remain unpopulated. I've tried several different combinations of settings for e-density and SetMIteration, but nothing seems to help.

It seems that some highly excited states are actually populated, like ILEV 258 and 262, which belong to the same configuration. Is the a buffer overflow due to strange configurations or high J's again? As before, help is greatly appreciated!

Dear all,

Just wanted to highlight that there seems to be an overall confusion in the code, readme's etc. when referring to the type of parallelization - in a few instances I've come across, OpenMP (shared memory multiprocessing) is used to explain what's really implemented in MPI (message passing interface). Like this bit in the main readme:

2-4. Parallel computation
Some of the functions have a parallel version using MPI. You can build with MPI enabled using the option --with-mpi=***, where *** is the MPI implementation installed on your machine. It has been tested with lammpi and openmp. If you are using openmp, specify --with-mpi=omp.

where OpenMP and MPI is clearly mixed up. Or in codes such as mpiutil.c

178 void ResetWidMPI(void) {
179 #if USE_MPI == 2
180   if (!MPIReady()) {
181     printf("openmp not initialized\n");
182     Abort(1);
183   }
184 #endif

where openmp is given in the error message, while it's actually MPI that is tested?

This is a classic confusion for all of us I guess, since the most common implementation of MPI is called open-mpi. The code seems to use both MPI (open-mpi) and OpenMP, which is really good of course for HPC computations if used correctly (OpenMP for local machine vecorization and MPI for node communication). I guess the important bit to understand is wether InitializeMPI() defines no. of OpenMP threads or MPI processes? It seems to be the latter.

Maybe this doesn't matter too much, and I'm sure you're aware of it :) Just thought it deserved to be mentioned as I was a bit confused there for a while.

Thanks for a good job, looking forward to test the code out! Cheers!
Jon Grumer
The CompAS collaboration (grasp, atsp etc)
Theoretical Astrophysics, Uppsala University

The documentation(generated from .tex to .pdf) describes the format of DB_RR(readable form) as follows:

2.2.4 DB RR
...
#bound 2J free 2J Delta E L 
7 1 0 0 2.0465E+03 0 
# the parameters in the fitting formula 
3.8124E-02 4.9724E+00 1.2195E+00 2.1768E+03 
#user egrid RR cross sec. PI cross sec. gf 
1.0069E+02 1.7567E-01 1.9604E+00 3.0537E-02 
1.2357E+03 1.4375E-02 8.4255E-01 1.3124E-02 
2.9799E+03 5.5123E-03 3.3223E-01 5.1751E-03 
5.6604E+03 2.4847E-03 1.2100E-01 1.8848E-03 
9.7797E+03 1.1476E-03 4.1004E-02 6.3872E-04 
1.6110E+04 5.2175E-04 1.3029E-02 2.0295E-04 
8 1 0 0 1.9975E+03 1 
3.4223E-02 5.3145E+00 1.2206E+00 2.1537E+03 
...

So the 2nd, 3rd, 4th columns are supposed to be RR cross section, PI cross section and something called gf respectively.
If I use PrintTable, everything goes fine, But when I use InterpCross for interpolation(and set m to 0), it becomes five columns as follows:

 1.0140E+02  1.1500E+02  1.9520E-02  3.1983E-04  1.2531E+00
 1.0240E+02  1.1600E+02  1.9001E-02  3.1368E-04  1.2198E+00
 1.0340E+02  1.1700E+02  1.8500E-02  3.0769E-04  1.1877E+00
 1.0440E+02  1.1800E+02  1.8016E-02  3.0187E-04  1.1566E+00
 1.0540E+02  1.1900E+02  1.7549E-02  2.9620E-04  1.1266E+00
 1.0640E+02  1.2000E+02  1.7097E-02  2.9069E-04  1.0976E+00
 1.0740E+02  1.2100E+02  1.6660E-02  2.8532E-04  1.0695E+00
 1.0840E+02  1.2200E+02  1.6238E-02  2.8010E-04  1.0424E+00
 1.0940E+02  1.2300E+02  1.5829E-02  2.7501E-04  1.0162E+00
 1.1040E+02  1.2400E+02  1.5434E-02  2.7006E-04  9.9084E-01

Turns out they are electron energy, photon energy, gf factor, RR cross section, and PI cross section respectively. The order has been changed.(checked with TotalPICross, TotalRRCross, PrintTable, and compared with another HFS-based program)

It might not even be an issue, but still. And if this is not repeatable, then just ignore it.

I noticed that if we specify multipole (m) in TransitionTable (even I set the default value -1), the value of strength in the binary file changes.
The ascii file does not change with or without multipole argument.

I suppose there were some updates in the binary file format?
Here is the relevant page of the document.

Hi,

I am trying to get the Hamiltonian of the system using StructureMBPT of the first type as in the manual and couldn't get the output in ACSI format which is named hamilton in my example.

Here is the example
##########################
from pfac import fac
fac.SetAtom('C')
fac.Config('1s2 2s2','1s2 2s1 2p1','1s2 2p2', group = 'n3')
fac.ConfigEnergy(0)
fac.OptimizeRadial(['n3'])
fac.ConfigEnergy(1)
n11 = [1,2]
n22 = [3,4,5,6]
fac.StructureMBPT('c34.lev.b',"hamilton",['n3'],3,n11,n22,1)
fac.MemENTable('c34.lev.b')
fac.PrintTable('c34.lev.b', 'c34.lev2', 1)
#############################

Dear Prof. Gu and Fuji-san,

Are there some instructions on how to use parallelized FAC code?
I would like to try to use openmp for the calculation.

If no, is it possible to add them to the manual?

By the way, the manual now is missing in the master branch.

I just noticed SetCILevel(-1) makes core dumped with the latest master.

The working example is as follows.
I will look inside the code later

SetAtom('O')
Closed('1s')
Config('g1', '2*6')
Config('g3', '2*5 3*1')
Config('g4', '2*5 4*1')
Config('ion', '2*5')
ConfigEnergy(0)
OptimizeRadial('g1')
ConfigEnergy(1)
SetCILevel(-1)
Structure('O_like/O_woCI.en', ['g1'])
Structure('O_like/O_woCI.en', ['g3'])
Structure('O_like/O_woCI.en', ['g4'])
Structure('O_like/O_woCI.en', ['ion'])

Hi,
I used FAC 1.1.5 and FAC 1.1.1 to calculate the collision strengths between levels of Beryllium-like Oxygen OV and it works well, but the results obtained by FAC 1.1.5 are different by more than 20 percent compared to the FAC 1.1.1 versions.

FAC 1.1.5
    2  2     13  2  6.2193E+01 1
 5.8390E+00  3.8615E-01  
     2  2     14  4  6.2202E+01 1
 5.8390E+00  6.7959E-02 
     3  4     13  2  6.2156E+01 1
 5.8690E+00  6.8104E-02  
     3  4     14  4  6.2165E+01 1
 5.8590E+00  6.8839E-01

FAC 1.1.1
     2  2     13  2  6.362E+01 1
 5.8390E+00  4.762E-01  
     2  2     14  4  6.2369E+01 1
  5.8390E+00  7.517E-02 
     3  4     13  2  6.2327E+01 1
 5.8690E+00  7.531E-02  
     3  4     14  4  6.2334E+01 1
 5.8590E+00  8.420E-01

The calculation of fe17_excitation.py script is the same comparted to ref result and comparted to the FAC 1.1.1 versions result.
Could you tell me is it normal or not?

Sincerely

the new setup.py no longer respects the CC env variables to override the c compiler. How can we specify a different compiler than the default used by python?

On a fresh install (and also on a month old install) I get a segmentation fault when running:

from pfac import fac

fac.SetAtom('Pr')
fac.Closed('1s', '2s', '2p', '3s', '3p', '3d', '4s', '4p', '4d')

fac.Config('n1', '5s2 5p2')
fac.Config('n2', '4f1 5s2 5p1')
fac.Config('n3', '5s2 5p1 5d1')
fac.Config('n4', '4f2 5s2')
fac.Config('n5', '4f1 5s2 5d1')
fac.Config('n6', '4f1 5s1 5p2')
fac.Config('n7', '4f2 5s1 5p1')
fac.Config('n8', '5s1 5p3')

fac.ConfigEnergy(0)
fac.OptimizeRadial('n1')
fac.ConfigEnergy(1)

print 'Structure'

fac.Structure('9plus_bin.en', ['n1', 'n2', 'n3', 'n4', 'n5', 'n6', 'n7', 'n8'])
fac.MemENTable('9plus_bin.en')
fac.PrintTable('9plus_bin.en', '9plus.en')

print 'CETable'

fac.CETableMSub('9plus_MSub_bin.ce', ['n7'],['n2'])

It occurs for this specific set of configurations (n2 to n7), other sets I tested run fine though. Any idea what causes this? Help is much appreciated!

What is the output of MaxwellRate?

This is a part of the output of MaxwellRate('Ne03b.ce', Ne03b.cer', -1, -1, 10.0) in the demo directory,

#     1	 1	    2	 3	 2.0010E-01	    1	1
 1.0000E+01	 2.5776E-01	 3.1999E+01


#     0	 1	    1	 1	 1.5998E+01	    1	1
 1.0000E+01	 1.1336E+00	 2.8994E+01

From the source code below,

Hi all,

Feel bad for asking this, but really tried to work things out myself and feel I don't get much further without a lot of effort. Could I please ask someone to clarify the current table output from the various RECORD's for me, the RT RECORD in particular.

I'm trying to learn how to dig out the various level in- and out-flows of a CRM-model, limited to say electron collision excitation/de-excitation and radiative decays. The RT record seems helpful here but can't make sense of the table produced by RateTable().

Take the NeL demo case as an example:

              NB          TR          CE          RR          AI          CI
     0    0   1.0000E+00  2.0000E+00  0.0000E+00  3.0000E+00  1.0000E+00 -1.0000E+00 1*2.2*1
     1    1   8.2916E-08  2.0000E+00  1.5998E+01  3.0000E+00  1.0000E+00 -1.0000E+00 1*2.2*1
     2    2   1.5890E-07  4.0000E+00  1.6198E+01  3.0000E+00  1.0000E+00 -1.0000E+00 1*2.2*1
     3    3   3.0589E-11  2.0000E+00  1.3631E+02  3.0000E+00  1.0000E+00 -1.0000E+00 1*2.3*1
     4    4   5.6260E-12  2.0000E+00  1.4073E+02  3.0000E+00  1.0000E+00 -1.0000E+00 1*2.3*1

In column order:

• dir should give the direction (0:in, 1:out)) of the transition but seems to be just level index?
• iblock should be index of the superlevel block from which the rates originate?
• NB gives level populations - correct?
• TR should be total rad trans rates - but seems to give 2J+1?
• CE should be col exc rates - but seems to give the excitation energy?
• RR should be rad rec/photoionization rates - but seems to give the number of electrons?
• AI should be autoionization rates - gives 1?
• CI should be col ionization - gives -1?
• icomplex is correct I guess?

I guess AI and CI could be 1 and -1 to show that these contributions are not included.

I have a similar issue, although not as severe, with understanding other RECORDS and HEADERS fully, most notably the SP record, but let's start with RT.

I know that the manual is outdated, and I'm a total FAC noob learning, but these outputs just seems to be way to off?

Even a short answer would be much appreciated :) Maybe @fujiisoup knows?

All the best,
Jon

Hi,

I have been working a bit with FAC lately and I find that it would be very convenient to be able to print the full electron configuration in the verbose ASCII level file (DB_EN). Eliminating all empty and full shells makes it unnecessarily hard to reconstruct the actual electron configuration of each level (programmatically).

When I look at the file in person I manage to figure out what the distribution should be, but automating this is very difficult.

Is there a specific reason why they aren't printed fully? (after all they are somewhat accessible with the GetConfigNR() and ListConfig()) methods in pfac, but even with those lists it is hard to match them automatically, since different configurations get mapped onto the same identifier when dropping the open and full orbitals.)

Are you ware of any workaround, and do you think that this could be a convenience feature to add in the future?

Cheers and thanks for this very versatile code!

Dear all:
When I was using the CR Model of FAC, I found some problems. I want to call the AddIon function multiple times to make more than one ion in the model.And i know the order must be such that the number of electrons are in consecutive increasing order.But when I use it, there is an error like this:ERROR: Ion charge state does not match 0 0 11 10.I don't know how to solve it, hope someone can help me.
And my script is:
z = 13
nele = 11
nele1 = 10
temp = 10
dens = 1.0
# a is the atomic symbol corresponding to z
a = fac.ATOMICSYMBOL[z]
p = '%s%02d'%(a, nele)
add the ion with k number of electrons
AddIon(nele1, 0.0, '%sb'%p)
AddIon(nele, 0.0, '%sb'%p)
SetBlocks(-1)
SetEleDist(0, temp, -1, -1)

Hello,
I haven't found a description in the manual of the functions output. The output are 9 tables!
Thank you

check still fails. it's comparing with an old copy of ref/* data. where is it getting those old data from?

scrm crashes with segfault after #343 for some scripts.
I will send reproducible examples soon.

Hello All,
I am trying to do a mix mode calculation which is done in a new release.
I am giving UTA(0) for few groups and UTA(1) for few groups. But it is taking only the UTA type which I have given later. please suggest me weather mix mode is possible in 1.1.5 latest release.

hello,
I am tying InitialMPI(2) for the excitation calculation.
but the run time without using InitialMPI is the same as without using InitialMPI(2).
and also when I am seeing the processors in use while running is only one showing only one.
please tell whether MPI is implemented in excitation or structure calculation program.
I am using the same script
use_openmp = False
if len(sys.argv) == 2 and sys.argv[1] == 'openmp':
use_openmp = True
if use_openmp:
# enable openmp with 2 cores
InitializeMPI(2)
.
.
.
if use_openmp:
FinalizeMPI()

Thanks!!!

I'm sorry, but I happen to notice that Structure acts weird when giving different configurations their reletivistic names. Here's an example for Si and Si+'s configurations.(ASCII form)

FAC 1.1.5
Endian	= 0
TSess	= 1568247079
Type	= 1
Verbose	= 1
Si Z	=  14.0
NBlocks	= 2
E0	= 0, -7.87413423E+03

NELE	= 14
NLEV	= 5
  ILEV  IBASE    ENERGY       P   VNL   2J
     0     -1  0.00000000E+00 0   301    0 1*2.2*8.3*4                      3p2                                              3p-2(0)0
     1     -1  1.00789669E-02 0   301    2 1*2.2*8.3*4                      3p2                                              3p-1(1)1.3p+1(3)2
     2     -1  2.95059453E-02 0   301    4 1*2.2*8.3*4                      3p2                                              3p+2(4)4
     3     -1  1.01392715E+00 0   301    4 1*2.2*8.3*4                      3p2                                              3p-1(1)1.3p+1(3)4
     4     -1  2.50507337E+00 0   301    0 1*2.2*8.3*4                      3p2                                              3p+2(0)0

NELE	= 13
NLEV	= 47
  ILEV  IBASE    ENERGY       P   VNL   2J
     5     -1  1.14765135E+01 1   301    1 1*2.2*8.3*3                      3p1                                              3p-1(1)1
     6     -1  1.15032211E+01 1   301    3 1*2.2*8.3*3                      3p1                                              3p+1(3)3
     7     -1  1.56048383E+01 0   301    3 1*2.2*8.3*3                      3s1.3p2                                          3s+1(1)1.3p-1(1)0.3p+1(3)3
     8     -1  1.56458258E+01 0   301    5 1*2.2*8.3*3                      3s1.3p2                                          3s+1(1)1.3p+2(4)5
     9     -1  1.56532582E+01 0   301    1 1*2.2*8.3*3                      3s1.3p2                                          3s+1(1)1
    10     -1  1.86042573E+01 0   301    3 1*2.2*8.3*3                      3s1.3p2                                          3s+1(1)1.3p-1(1)0.3p+1(3)3
    11     -1  1.86120359E+01 0   301    5 1*2.2*8.3*3                      3s1.3p2                                          3s+1(1)1.3p-1(1)2.3p+1(3)5
    12     -1  2.01549603E+01 0   301    1 1*2.2*8.3*3                      3s1.3p2                                          3s+1(1)1.3p+2(0)1
    13     -1  2.16084508E+01 0   301    1 1*2.2*8.3*3                      3s1.3p2                                          3s+1(1)1.3p-1(1)2.3p+1(3)1
    14     -1  2.16169392E+01 0   301    3 1*2.2*8.3*3                      3s1.3p2                                          3s+1(1)1.3p+2(4)3
    15     -1  1.07595671E+02 1   301    5 1*2.2*7.3*4                      2p5.3p2                                          2p+3(3)3.3p-1(1)2.3p+1(3)5
  ...

It can be seen that ILEV=7 and ILEV=10 have the same reletivistic configuration name, but they have different energy, suggesting Structure name them incorrectly.
Maybe their names are essentially the same? However, when BasisTable is used, it gives a different result.

============Basis Table===========================
#    0    0  0       2
     0    0  0       0   0     0     0   1*2.2*8.3*4                      3p2                                              3p-2(0)0
     0    0  0       1   0     2     0   1*2.2*8.3*4                      3p2                                              3p+2(0)0

#    2    0  1       9
     2    0  1       0   1     0     0   1*1.2*8.3*4                      1s1.3p2                                          1s+1(1)1
     2    0  1       1   1     1     7   1*1.2*8.3*4                      1s1.3p2                                          1s+1(1)1.3p-1(1)2.3p+1(3)1
     2    0  1       2   1     2     0   1*1.2*8.3*4                      1s1.3p2                                          1s+1(1)1.3p+2(0)1
     2    0  1       3   1     3     0   1*2.2*7.3*4                      2s1.3p2                                          2s+1(1)1
     2    0  1       4   1     4     7   1*2.2*7.3*4                      2s1.3p2                                          2s+1(1)1.3p-1(1)2.3p+1(3)1
     2    0  1       5   1     5     0   1*2.2*7.3*4                      2s1.3p2                                          2s+1(1)1.3p+2(0)1
     2    0  1       6   1    12     0   1*2.2*8.3*3                      3s1.3p2                                          3s+1(1)1
     2    0  1       7   1    13     7   1*2.2*8.3*3                      3s1.3p2                                          3s+1(1)1.3p-1(1)2.3p+1(3)1
     2    0  1       8   1    14     0   1*2.2*8.3*3                      3s1.3p2                                          3s+1(1)1.3p+2(0)1

#    3    1  1       7
     3    1  1       0   1     7     0   1*2.2*7.3*4                      2p5.3p2                                          2p-1(1)1
     3    1  1       1   1     8     0   1*2.2*7.3*4                      2p5.3p2                                          2p+3(3)3.3p-1(1)2.3p+1(3)1
     3    1  1       2   1     8    21   1*2.2*7.3*4                      2p5.3p2                                          2p+3(3)3.3p-1(1)4.3p+1(3)1
     3    1  1       3   1     9     7   1*2.2*7.3*4                      2p5.3p2                                          2p-1(1)1.3p-1(1)2.3p+1(3)1
     3    1  1       4   1    10     6   1*2.2*7.3*4                      2p5.3p2                                          2p+3(3)3.3p+2(4)1
     3    1  1       5   1    11     0   1*2.2*7.3*4                      2p5.3p2                                          2p-1(1)1.3p+2(0)1
     3    1  1       6   1    16     0   1*2.2*8.3*3                      3p1                                              3p-1(1)1

#    4    0  2       1
     4    0  2       0   0     1     0   1*2.2*8.3*4                      3p2                                              3p-1(1)1.3p+1(3)2

#    6    0  3       9
     6    0  3       0   1     1     0   1*1.2*8.3*4                      1s1.3p2                                          1s+1(1)1.3p-1(1)0.3p+1(3)3
     6    0  3       1   1     1    14   1*1.2*8.3*4                      1s1.3p2                                          1s+1(1)1.3p-1(1)2.3p+1(3)3
     6    0  3       2   1     2     6   1*1.2*8.3*4                      1s1.3p2                                          1s+1(1)1.3p+2(4)3
     6    0  3       3   1     4     0   1*2.2*7.3*4                      2s1.3p2                                          2s+1(1)1.3p-1(1)0.3p+1(3)3
     6    0  3       4   1     4    14   1*2.2*7.3*4                      2s1.3p2                                          2s+1(1)1.3p-1(1)2.3p+1(3)3
     6    0  3       5   1     5     6   1*2.2*7.3*4                      2s1.3p2                                          2s+1(1)1.3p+2(4)3
     6    0  3       6   1    13     0   1*2.2*8.3*3                      3s1.3p2                                          3s+1(1)1.3p-1(1)0.3p+1(3)3
     6    0  3       7   1    13    14   1*2.2*8.3*3                      3s1.3p2                                          3s+1(1)1.3p-1(1)2.3p+1(3)3
     6    0  3       8   1    14     6   1*2.2*8.3*3                      3s1.3p2                                          3s+1(1)1.3p+2(4)3

#    7    1  3       9
     7    1  3       0   1     6     0   1*2.2*7.3*4                      2p5.3p2                                          2p+3(3)3
     7    1  3       1   1     8     7   1*2.2*7.3*4                      2p5.3p2                                          2p+3(3)3.3p-1(1)2.3p+1(3)3
     7    1  3       2   1     8    28   1*2.2*7.3*4                      2p5.3p2                                          2p+3(3)3.3p-1(1)4.3p+1(3)3
     7    1  3       3   1     9     0   1*2.2*7.3*4                      2p5.3p2                                          2p-1(1)1.3p-1(1)0.3p+1(3)3
     7    1  3       4   1     9    14   1*2.2*7.3*4                      2p5.3p2                                          2p-1(1)1.3p-1(1)2.3p+1(3)3
     7    1  3       5   1    10     0   1*2.2*7.3*4                      2p5.3p2                                          2p+3(3)3.3p+2(0)3
     7    1  3       6   1    10    12   1*2.2*7.3*4                      2p5.3p2                                          2p+3(3)3.3p+2(4)3
     7    1  3       7   1    11     6   1*2.2*7.3*4                      2p5.3p2                                          2p-1(1)1.3p+2(4)3
     7    1  3       8   1    15     0   1*2.2*8.3*3                      3p1                                              3p+1(3)3

#    8    0  4       2
     8    0  4       0   0     1     7   1*2.2*8.3*4                      3p2                                              3p-1(1)1.3p+1(3)4
     8    0  4       1   0     2     6   1*2.2*8.3*4                      3p2                                              3p+2(4)4

#   10    0  5       6
    10    0  5       0   1     1    21   1*1.2*8.3*4                      1s1.3p2                                          1s+1(1)1.3p-1(1)2.3p+1(3)5
    10    0  5       1   1     2    12   1*1.2*8.3*4                      1s1.3p2                                          1s+1(1)1.3p+2(4)5
    10    0  5       2   1     4    21   1*2.2*7.3*4                      2s1.3p2                                          2s+1(1)1.3p-1(1)2.3p+1(3)5
    10    0  5       3   1     5    12   1*2.2*7.3*4                      2s1.3p2                                          2s+1(1)1.3p+2(4)5
    10    0  5       4   1    13    21   1*2.2*8.3*3                      3s1.3p2                                          3s+1(1)1.3p-1(1)2.3p+1(3)5
    10    0  5       5   1    14    12   1*2.2*8.3*3                      3s1.3p2                                          3s+1(1)1.3p+2(4)5

#   11    1  5       5
    11    1  5       0   1     8    14   1*2.2*7.3*4                      2p5.3p2                                          2p+3(3)3.3p-1(1)2.3p+1(3)5
    11    1  5       1   1     8    35   1*2.2*7.3*4                      2p5.3p2                                          2p+3(3)3.3p-1(1)4.3p+1(3)5
    11    1  5       2   1     9    21   1*2.2*7.3*4                      2p5.3p2                                          2p-1(1)1.3p-1(1)2.3p+1(3)5
    11    1  5       3   1    10    18   1*2.2*7.3*4                      2p5.3p2                                          2p+3(3)3.3p+2(4)5
    11    1  5       4   1    11    12   1*2.2*7.3*4                      2p5.3p2                                          2p-1(1)1.3p+2(4)5

#   15    1  7       2
    15    1  7       0   1     8    42   1*2.2*7.3*4                      2p5.3p2                                          2p+3(3)3.3p-1(1)4.3p+1(3)7
    15    1  7       1   1    10    24   1*2.2*7.3*4                      2p5.3p2                                          2p+3(3)3.3p+2(4)7

In the file generated by BasisTable, there are no longer same reletivistic configuration names. And the total configuration number is the same as Structure's. Therefore, these configurations should be the same ones in Structure's output file. In short, Struture is indeed naming configurations incorrectly.

Now I want to find out the real names of these configurations via BasisTable. But because the "Basis Table" part doesn't have ILEV or energy information, this can't be done. Actually the "Mixing Coefficients" part of BasisTable's output file(not shown above) does contain ILEV and energy information, but it doesn't contain configuration information, or I can't deduce them from those columns yet.

Structure's issues mentioned above is quite common for multi-electron atoms(ions). Hope these issues be fixed soon. Also an explanation of BasisTable's output file's structure would be really helpful.

Thank you very much!

Probably the documentation for AI_header is outdated.
Did we remove the entry channel and new entry emin?

I think PR#9 breaks some parts of the code.
@mfgu , can you investigate this?

To reproduce this, do

cd demo/structure
make

It works until PR#9 was merged.

Sorry,I want to change UTYPE from 1 to 0 in CI_Table, but I can't find any state about how to change UTYPE

version: FAC 1.1.5

Hello, developer：

I have some problems when I calculated a spectral model involving more than one ion using function spectrum(). In fact, the ions with k electrons and k-1 electrons are introduced by AddIon() and SetBlock() respectively.
NormalizeMode was set to 2.
EN, CE, TR, CI, RR, AI data of ion k and EN, CE, TR data of ion k-1 has already calculated.

1. segfault & bus error
   However, if the CE data of ion k-1 is larger than hundreds of KBs (maybe 500K), errors occur when SetCERates() is called, mainly segfault and sometimes bus error.
   I find that the segfault in SetBlock() has been fixed last year. I'm not sure whether that was just the BUG I encountered.

2. mismatched population data in different parts
   2.1 two ions
   To avoid the segfualt or bus error, I reduced the number of the configurations involved in CE. Though the code run successfully, new problem occurred. For an example:
   EN data of ion k :

NELE    = 68
NLEV    = 3192
...
NELE    = 67
NLEV    = 493 
 ILEV  IBASE    ENERGY       P   VNL   2J
  3192     -1  1.19279240E+02 1   501    7 4*31.5*8                         4f13                                             4f+7(7)7
  3193     -1  1.21438905E+02 1   501    5 4*31.5*8                         4f13                                             4f-5(5)5
...
  3673     -1  2.24462514E+02 1   503    5 4*31.5*8                    4f13.5p5.5f1                                     4f-5(5)5.5p-1(1)4.5f+1(7)5
...

And EN data of ion k-1

NELE    = 67
NLEV    = 493
  ILEV  IBASE    ENERGY       P   VNL   2J
     0     -1  0.00000000E+00 1   501    7 4*31.5*8                         4f13                                             4f+7(7)7
     1     -1  2.15966489E+00 1   501    5 4*31.5*8                         4f13                                             4f-5(5)5
...
    481     -1  1.05183273E+02 1   503    5 4*31.5*8                    4f13.5p5.5f1                                     4f-5(5)5.5p-1(1)4.5f+1(7)5
...

After CRM, part of SP data is followed:

#population
NELE    = 67
NTRANS  = 1
TYPE    = 0000000
IBLK    = 0
ICOMP   = 4*31.5*8
FBLK    = 0
FCOMP   =
  3192      0  1.192792E+02  9.8041E-01  8.0000E+00  1.5485E-03
...
  3193      0  1.214389E+02  8.0384E-12  6.0000E+00  4.0073E+07

...
  3672     19  2.240807E+02  1.7817E-13  4.0000E+00  1.2970E+08
  3673     20  2.244625E+02  3.9569E-05  6.0000E+00  9.7364E+01
  3674     21  2.247520E+02  0.0000E+00  1.2000E+01  0.0000E+00

#line emissivity TR
NELE    = 67
NTRANS  = 3
TYPE    = 0000505
IBLK    = 138
ICOMP   = 4*31.5*8
FBLK    = 5
FCOMP   = 4*31.5*8
     1      0  2.159665E+00  3.2163E-03  8.1284E+01  9.7364E+01
#population = 3.9569E-05.

#line emissivity RR
NELE    = 68
NTRANS  = 1
TYPE    = 0000005
IBLK    = 1
ICOMP   = 4*31.5*8
FBLK    = 190
FCOMP   = 4*31.5*9
3193     50  4.350177E+01  8.1631E-24  1.0155E-12  4.0073E+07
#population = 8.0385E-12

According to the population part, the population of level 3193 belonging to ion k-1 should be 8.0384E-12(It’s too small. It should be comparable to 3192). But the TR 3193→3192 (1→0) in line emissivity part suggested that its population is 3.2163E-03/ 8.1284E +01 = 3.9569E-05 which was the population of level 3673 (3673’s population shouldn’t be so large, I think.). It's not a coincidence since the population is still matched level 3673 when I chose different temperature, though there might be a very little deviation. But the matched level would change if base data was different. And I checked some other levels of ion k-1, they all had such phenomenon.
However, the population calculated from the RR 3193→50 is 8.0385E-12. It was equal to the population part.
If I reduced the total number of levels of ion k-1 , the problem could be avoided. But no more than about 200~300 levels...

2.2 three ions
The numbers of levels and CE transitions of ion k-1 are limited. But the ion k-1, w7+, is my interested ion. I modified k from 68 to 67. Therefore, ion k → ion k+1, ion k-1 → ion k. They are both introduced by AddIon(). And another ion, ion k-1, that is w8+, was introduced by SetBlock(). Therefore, more CE transitions of W7+ were allowed. However, except the error happened on ion k-1, new problem occurred again. It's opposite to the former one. The population calculated from RR process between ion k and ion k+1 matched the wrong level of ion k. But the TR process was right.

NELE	= 67                                                     
NTRANS	= 1                                                    
TYPE	= 0000000                                                
IBLK	= 65                                                     
ICOMP	= 4*31.5*8                                               
FBLK	= 0                                                      
FCOMP	=                                                        
     0      0  1.192792E+02  5.0628E-01  8.0000E+00  2.4123E+02
…
     1      0  1.214389E+02  1.1523E-01  6.0000E+00  3.2322E+02
…
     5      0  1.616404E+02  2.3321E-07  8.0000E+00  6.5483E+06
#lev0: 4f13   lev5: 4f13.5p5.5d1

NELE	= 68                                                     
NTRANS	= 2416                                                 
TYPE	= 0000005                                                
IBLK	= 70                                                     
ICOMP	= 4*31.5*8                                               
FBLK	= 853                                                    
FCOMP	= 4*31.5*9                                               
3192    640  2.915478E+00  1.3869E-18  5.9468E-12  6.5483E+06
#population = 2.3322E-07

This phenomenon can’t be repeated in the demo script involving Fe23+, Fe22+, Fe21+(NELE=3, 4, 5). But I can’t avoid that in my own script even if I reduced the levels of ion k in the EN data of ion k+1 to only 5. And the matching result in that case is shown below.

     level number in the transition         level number whose (direct) population equals to the  calculated population 
                    3192 (0)               →                              6
                    3193 (1)               →                              5
                    3194 (2)               →                              0
                    3195 (3)               →                              1
                    3196 (4)               →                              7

There're different phenomenon in different situations. I'm not sure whether there're also similar phenomenon in CE, CI, AI processes since we can't see that intuitively.

3. NormalizeMode
   By the way, does the norm_mode=2 mean that the total population is fixed to the total number set by SetAbund() and the populations of different ions can be compared? And I can't understand the result when m=3. What's the meaning when m=3?

Hi,

I have seen that you started using Travis and that the release source code packages have disappeared (the link in the readme is broken at least). Does this mean that you are switching over to a more continuous release model and I should download the source code directly from the repository for buliding fac in the future?

Also I saw that the python 3.6 builds are passing on Travis - does this mean fac is ready to be used with python 3 now?

Keep up the good work guys, cheers!

Results of cr model change depending on the order of Structure executed.

For the example shown below, the calculated spectral energies and intensities are

They show a significant difference.
.en, .tr files look (effectively) the same.

The following is an example script that generates the above result

SetAtom('Si')
Closed('1s')
Closed('2s')
Closed('2p')

Config('lev1', '3s2 3p1')
Config('lev2', '3s2 3d1')
Config('lev3', '3s2 4s1')

ConfigEnergy(0)
OptimizeRadial(['lev1'])
ConfigEnergy(1)

Structure('text2.en',['lev1','lev2'])  # line 1
Structure('text2.en',['lev3'])  # line 2
MemENTable('text2.en')

TransitionTable('text2.tr', ['lev1', 'lev2', 'lev3'], ['lev1', 'lev2', 'lev3'])
CETable('text2.ce', ['lev1'], ['lev1', 'lev2', 'lev3'])

AddIon(13,1,text2)
SetBlocks(-1)
SetEleDist(0, 1000, -1, -1)
SetTRRates(0)
SetCERates(1)
SetAbund(13, 1.0)
SetEleDensity(1000)
InitBlocks()
SetIteration(1e-4, 0.5)
LevelPopulation()
SpecTable('text2.sp',0)

For example 2, we interchanged line 1 and line 2 and kept the rest same.

Hi,
I am PHD student and i'am working in EIE , one part of my work assigned to modified excitation programme (collision strength ) .
But i have problem , when i try to do a very small change , the fac work normaly also when i deleted excitation.c the fac work
Could you help me please ?.

Dear Dr.Gu,
I am currently using the wave function output from PFAC, but I am not quite sure about the units of the amplitude for the continuum radial wave functions.

According to the review of "the flexible atomic code" the orthogonality condition is
Integral dr [P(E)P(E')+Q(E)Q(E')]=Pi*Delta(E-E').
Is that right? Does this mean that the Continuum wave function are dimensionless?
What about the bound state wave functions? Is the energy of the amplitudes in terms of Hartree or eV or something different?
Best Regards

hello all,

while running C-R model i found

SetBlocks{den, pref}
which setup the superlevel blocks for the spectral model.
The optional {den} is the abundance of the ion with one less
electron than the lowest charge state included in the model, and {pref} is
the base file name for that ion.

As I understand pref is the base file ion with one less
electron than the lowest charge state included in the model.

but when I an running this, for example I am running the example given in demo/FeL for fe(with 9,11 electrons), atomic data is generated successfully. however when I am running the spec.py neles = range(9,11) is is showing

ERROR: Ionized block 4 of ion 9 does not match a block in file data/Fe08b.en
nlevels = 290 VS -1

I also want to know that does the present CR model includes the ionization balance calculations when I consider more than one ion?

Hi @mfgu

Did you change the unit of the values in .en (binary) file?
Values changed depending on the version of fac by a factor of 1e28.

Hi, @mfgu
I often encountered segfault in scrm.
A minimal reproducing script is as follows;
Some scripts are working but some don't.
I have no idea about the cause.

FeII.sf

InitializeMPI(8)
SetAtom('Fe')
Config('3d5.4s1.4p1', '1s2 2s2 2p6 3s2 3p6 3d5 4s1 4p1')
ConfigEnergy(0)
OptimizeRadial('3d5.4s1.4p1')
ConfigEnergy(1)
Structure('FeII.en', ['3d5.4s1.4p1'])
MemENTable('FeII.en')
TransitionTable('FeII.tr', ['3d5.4s1.4p1'], ['3d5.4s1.4p1'])
CETable('FeII.ce', ['3d5.4s1.4p1'], ['3d5.4s1.4p1'])
FinalizeMPI(8)

FeII.crm

InitializeMPI(8)
AddIon(26, 1, 'FeII')
SetBlocks(-1)
SetEleDist(0, 1.0, -1, -1)
SetEleDensity(1.0)
SetTRRates(0)
SetCERates(1)
InitBlocks()
SetIteration(1e-6, 0.5, 1024)
LevelPopulation()
SpecTable('tmp.sp', 0)
FinalizeMPI()

Hello,

I have been trying to install the pfac interface and I'm facing very similar problem as described here (#264). The sfac interface perfectly works.

I'm using Python 2.7.16 (no anaconda), gcc as c compiler and f77 as fortran compiler.

gcc build/temp.linux-x86_64-2.7/python/fac.o -L/home/simone/Scrivania/fac-master -L/usr/lib/gcc/x86_64-linux-gnu/8 -L/usr/lib/gcc/x86_64-linux-gnu/8/../../../x86_64-linux-gnu -L/usr/lib/gcc/x86_64-linux-gnu/8/../../../../lib -L/lib/x86_64-linux-gnu -L/lib/../lib -L/usr/lib/x86_64-linux-gnu -L/usr/lib/../lib -L/usr/lib/gcc/x86_64-linux-gnu/8/../../.. -lfac -lmlapack -lmblas -lpthread -lgfortran -lm -lquadmath -o build/lib.linux-x86_64-2.7/pfac/fac.so -pthread -shared -fno-strict-aliasing -DNDEBUG -g -fwrapv -O2 -Wall -Wstrict-prototypes -Wdate-time -D_FORTIFY_SOURCE=2 -g -fdebug-prefix-map=/build/python2.7-S1JTsz/python2.7-2.7.16=. -fstack-protector-strong -Wformat -Werror=format-security /usr/bin/ld: /home/simone/Scrivania/fac-master/libfac.a(dcoul.o): relocation R_X86_64_PC32 against symbol ofcoul_' can not be used when making a shared object; recompile with -fPIC
/usr/bin/ld: link finale non riuscito: bad value
collect2: error: ld returned 1 exit status
error: command 'gcc' failed with exit status 1
make: *** [Makefile:42: pfac] Error 1
`
I didn't find any workaround. I tried to run the code till the 'build/lib.linux-x86_64-2.7/pfac/fac.so -pthread -shared' but I always obtain the Relocation error. Do you have any idea?

Many thanks.

Sometimes, TransitionTable does not compute the transition rate (I want an oscillator strength and its matrix element)

SetAtom('Fe')
Config('lev_3d6_4s2', '1s2 2s2 2p6 3s2 3p6 3d6 4s2')
Config('lev_3d7_4s1', '1s2 2s2 2p6 3s2 3p6 3d7 4s1')
Config('lev_3d8', '1s2 2s2 2p6 3s2 3p6 3d8')
ConfigEnergy(0)
OptimizeRadial('lev_3d6_4s2')
ConfigEnergy(1)
Structure('Fe_small.en', ['lev_3d6_4s2', 'lev_3d7_4s1', 'lev_3d8'])
MemENTable('Fe_small.en')
PrintTable('Fe_small.en', 'Fe_small.en.txt')
TransitionTable('Fe_small.tr', ['lev_3d6_4s2', 'lev_3d7_4s1', 'lev_3d8'], ['lev_3d6_4s2', 'lev_3d7_4s1', 'lev_3d8'], -1)
PrintTable('Fe_small.tr', 'Fe_small.tr.txt')

I got only the header and no valus in 'Fe_small.tr.txt'

FAC 1.1.5
Endian	= 0
TSess	= 1576472878
Type	= 2
Verbose	= 1
Fe Z	=  26.0
NBlocks	= 0

m=-2, 1, 3 only give values.
If I compute m=0, the result is similar to that computed from m=-2, 1, 3.
I don't think the electric dipole (m=-1) is that weak and electric quadrupole is that strong.
Or is there any bug in my script?

HI,
We still have problems with cross-section values between magnetic sub-levels of the OV ion, when I compared the values of the cross section between a magnetic sublevel to a magnetic sublevel for the transitions from (2s2p 3P2) and (2s3p 3P2) by the values calculated by Zhang, we found differences especially at the big energy levels for example t:

Lv. initial 2s2p 3P2 Ji = 2 Lvl. final 2s3p 3P2 Jf = 2
 Electron energy scattered = 61.0249 eV
                      

                    Zhang              FAC1.1.5

   Mi   Mf    Sigma Mi->Mf    Sigma Mi->Mf
   -2   -2        1.194E-19        1.453E-18
   -2   -1        1.398E-21        6.005E-20
   -2    0        2.684E-21        2.831E-20
   -2    1        4.900E-24        9.510E-21
   -2    2        1.838E-37        3.027E-31
   -1   -2        1.417E-21       4.939E-20
   -1   -1        8.604E-20       1.488E-18
   -1    0        2.435E-22       4.094E-20
   -1    1        4.020E-21       2.192E-20
   -1    2        4.900E-24       9.510E-01
    0   -2        2.682E-21       1.683E-20
    0   -1        2.682E-21       4.135E-20
    0    0        7.649E-20       1.508E-18

and thanks.

In #89, we made a backward incompatible change.
I think we can now consider releasing FAC 1.2.

In order to make a release, we may need

• fill up what's new entry
• complete documentation for the newly added functions

Any thoughts?

I noticed that OptimizeRadial does not update a, and thus a==0 is always used.
I checked both wih GetPotential and SavePotential functions.

Is it an intended behavior?
From the manual, it looks that both lambda and a should be updated.

The parameterized potential that is used in FAC is

By setting ip=1 in Structure(fn, [, hfn ], g [, p [, ip ] ]), fac ignores the non-diagonal element within p.

But how does it compute the eigenvalues from the perturbed hamiltonian?
I computed the eigenvalues from the hamiltonian but I found if ip=0, the eigenvalues of the output hamiltonian (computed by the direct eigenvalue decomposition) is consistent with the energy level in .en file in 1e-10 eV order level (the ionization energy is ~3000 eV),
but if ip==1, the eigenvalues from the perturbed hamiltonian are different from those in .en file in the order of 0.01 eV.

Here is the script I tested

SetAtom('W')
Closed('1s')
Closed('2s')
Closed('2p')

Config('g1_0', '3s2 3p4')
Config('g1_2', '3s2 3p3 3d1')
Config('g1_3', '3s2 3p3 4s1')
Config('g1_4', '3s2 3p3 4p1')
Config('g1_5', '3s2 3p3 4d1')

# configuration to be considered as a perturbation
Config('g2', '3s1 3p4 4p1')

ConfigEnergy(0)
OptimizeRadial('g1_0')
ConfigEnergy(1)

Structure('test_W.en', 'test_W.ham', ['g1_0', 'g1_2', 'g1_3', 'g1_4', 'g1_5'], ['g2'], 1)
# for full hamiltonian
# Structure('test_W.en', 'test_W.ham', ['g1_0', 'g1_2', 'g1_3', 'g1_4', 'g1_5'], ['g2'])  

Hi!

I wanted to extract the CE rates from the demo example on Fe 17 excitation python code but I encountered an issue when the file type to be read in the CERate function seems not to recognize the FAC binary format created just before by the CETable function.

I'm using the FAC version 1.1.5. I attach the piece of code from the demo folder where I added 2 lines to extract the CE rates:

temperature = [500, 700]
fac.CERate('ne.ce.b','ne.rce', -1, -1, temperature)

Just after the PrintTable of the CETable function and just before the function that cleans the parallel version of the run (which I don't use). The output is also in the zip file. Both the CE table and the energy level files match the reference files in the ref folder. The output error message I have when I run the fe17_excitation.py file modified with the CE rates is the following:

> python fe17_excitation.py
File ne.ce.b is not in FAC binary format

I'm surprised of this output since I used a similar piece of code using FAC version 1.1.0 years ago and it worked just fine, creating the RCE file with the rates inside. When I run the same old code on this new FAC version, I have the same error output, saying the binary format is not recognized as a FAC binary file. Here are my questions:

• Am I doing something wrong on the code used to create the CE rates?
• If not, could this be related to a problem in the way the header is being written on the binary file?
• or the way it's being read?

Thanks a lot for looking into this!

Eric

fe17_excitation.zip

I added a validation test with openmp in #16 , and I noticed that the transition order changes if we switch on openmp.

e.g. in demo/excitation/ne.ce, we got
without openmp, the transitions are orderd in the ascent manner

     0  0      1  4  7.2447E+02 1
-1.0000E+00  0.0000E+00  0.0000E+00
...
     0  0      2  2  7.2648E+02 1
 8.3122E-03 -6.1412E-03  9.4388E+04
...
     0  0      3  0  7.3714E+02 1
-1.0000E+00  0.0000E+00  0.0000E+00

but with openmp, it is randomly ordered

     0  0      1  4  7.2434E+02 1
-1.0000E+00  0.0000E+00  0.0000E+00
...
     0  0      3  0  7.3700E+02 1
-1.0000E+00  0.0000E+00  0.0000E+00
...
     0  0      4  2  7.3831E+02 1
 8.3792E-03 -6.1549E-03  9.4387E+04

My current test code does not take account the ordering, and it simply reports the results are different.
Is it an intended behavior?
If so, I will update the test code.

But I personally think it would be better for users if the transition order is kept same independent from how many cores are used.

The document says that a data block in .sp has 4 columns,
i.e. upper, lower, Delta E, emissivity

But the generated file by FAC==1.1.4 seems to have six columns, like

NELE	= 1
NTRANS	= 1
TYPE	= 0000201
IBLK	= 2
ICOMP	= 2*1 
FBLK	= 1
FCOMP	= 1*1 
     1      0  9.183927E+01  4.5097E-04  5.0776E+10  5.0776E+10

Can anyone give me a suggestion about the fifth and sixth columns?
I will update the documents when I have time.

I think by #21 , the ascii file is now sorted even when openmp is enabled.

But the new behavior looks somehow unstable.
The tests passed when the PR is sent, but failed after the merge commit.

In my local environment,

cd demo/NeL
make openmp

sometimes fails and causes a segfault (but sometimes succeeds).
I guess it is failing when writing to the file?

I think it is a good idea to list up what is missing in the current documentation.
This issue can be used for this purpose.

What I noticed is

• WaveFuncTable #133
• New option for Structure to save Hamiltonian
• RMBPT

Note

Please modify this comment and add items that should be added / modified to the documentation.

From #89, the automatic testings are failing.
@mfgu did you change something in the output?

I will update the testing code if necessary.

why do you want to import setuptool in the setup.py. this is not working on my system.
python setup.py install
does not install the module files into the site-package dir.

Hi, @mfgu
Thanks for keeping developing and maintaining this good library.

I would like to compute the oscillator strength between two states that are specified by given mixing coefficients.
These mixing coefficients may be different from those by FAC (but the basis functions should be the same).

I think if I have an overlap integral among all the basis functions, I may be able to compute it with the given mixing coefficients and simple algebraic operations.

Is it possible to implement such a function in FAC?

I looked inside the FAC code but it looks that it computes oscillator strength level by level, not basis by basis.
Is there any reason not to compute the overlap integral first?
I guess FAC computes the same overlap integrals for many times, but is this understanding correct?

I am writing my master thesis at the moment using the flexible atomic code and there are some questions concerning the output of the radial wave functions of the bound and continuum states:

-)Continuum wave function:
Since the essential output of the continuum wave function at order 1keV was only a fraction of 1 Bohr radius, I changed
some values in the file "consts.h" in the directory "faclib" in the following way:
--) I changed the values
#define MAXRP 3000 /* maximum radial mesh /
#define DMAXRP 1200 / default radial mesh points /
#define GRIDASYMP 36 / no. points in one wavelength near infinity /
to
#define MAXRP 3000 / maximum radial mesh /
#define DMAXRP 2500 / default radial mesh points /
#define GRIDASYMP 419 / no. points in one wavelength near infinity */
Is there a reason why these values are chosen to be the ones in the file, because the output looks fine to me?
-) Ionization: (I wrote a sequence of the code at the end)
--)If I want calculate the continuum states, does the Config() should include the electron configuration before or after the ionisation process?
--)Is there a main difference, if I use the function Closed() for the closed shells instead of Config() (?), because the function OptimizeRadial() only contains the groups I include in the argument as far as I understood?
--) Since in the manual it is written that one should NOT use the function OptimizeRadial() for levels bigger than n=4, does this has any influence on the WaveFuncTable() output?

Example: Ionization Xe-2p

fac.SetAtom('Xe')
fac.Config('1* 2', group='1s')
fac.Config('2* 2', group='2s')
fac.Config('2* 6', group='2p')
fac.Config('3* 2', group='3s')
fac.Config('3* 6', group='3p')
fac.Config('3* 10', group='3d')
fac.Closed('4s')
fac.Closed('4p')
fac.Closed('4d')
fac.Closed('5p')
fac.Closed('5s')
fac.ConfigEnergy(0)
fac.OptimizeRadial(['1s','2s','2p','3s','3p','3d'])
fac.ConfigEnergy(1)
fac.WaveFuncTable('FXe2pbd.txt', 2, -2)
fac.WaveFuncTable('FXe2p100.txt', 0, -2, 100)

Thank you very much in advance