Merge pull request #65 from TREX-CoE/fix-basis-structure

FIX: data format in the basis group consistent with ECP

configure.ac

@@ -2,7 +2,7 @@
 # Process this file with autoconf to produce a configure script.
 
 AC_PREREQ([2.69])
-AC_INIT([trexio], [1.1.0], [https://github.com/TREX-CoE/trexio/issues])
+AC_INIT([trexio], [2.0.0], [https://github.com/TREX-CoE/trexio/issues])
 AM_INIT_AUTOMAKE([subdir-objects color-tests parallel-tests silent-rules 1.11])
 AM_MAINTAINER_MODE()
 LT_PREREQ([2.2])

python/pytrexio/_version.py

@@ -1 +1 @@
-__version__ = "0.1.0"
+__version__ = "1.0.0"

python/test/test_api.py

@@ -71,11 +71,15 @@ charges_np = np.array(charges, dtype=np.int32)
 # from the size of the list/array by SWIG using typemaps from numpy.i
 trexio.write_nucleus_charge(test_file, charges_np)
 
+basis_shell_num = 24
 # initialize arrays of nuclear indices as a list and convert it to numpy array
-indices = [i for i in range(nucleus_num)]
+indices = [i for i in range(basis_shell_num)]
 # type cast is important here because by default numpy transforms a list of integers into int64 array
 indices_np = np.array(indices, dtype=np.int64)
 
+# first write basis_shell_num because it is needed to check dimensions of basis_nucleus_index in TREXIO >= 2.0.0
+trexio.write_basis_shell_num(test_file, basis_shell_num)
+
 # function call below works with both lists and numpy arrays, dimension needed for memory-safety is derived 
 # from the size of the list/array by SWIG using typemacs from numpy.i
 trexio.write_basis_nucleus_index(test_file, indices_np)
@@ -160,20 +164,20 @@ except trexio.Error:
     print("Unsafe call to safe API: checked")
 
 # safe call to read array of int values (nuclear indices)
-rindices_np_16 = trexio.read_basis_nucleus_index(test_file2, dim=nucleus_num, dtype=np.int16)
+rindices_np_16 = trexio.read_basis_nucleus_index(test_file2, dim=basis_shell_num, dtype=np.int16)
 assert rindices_np_16.dtype is np.dtype(np.int16)
-for i in range(nucleus_num):
+for i in range(basis_shell_num):
     assert rindices_np_16[i]==indices_np[i]
 
-rindices_np_32 = trexio.read_basis_nucleus_index(test_file2, dim=nucleus_num, dtype=np.int32)
+rindices_np_32 = trexio.read_basis_nucleus_index(test_file2, dim=basis_shell_num, dtype=np.int32)
 assert rindices_np_32.dtype is np.dtype(np.int32)
-for i in range(nucleus_num):
+for i in range(basis_shell_num):
     assert rindices_np_32[i]==indices_np[i]
 
 rindices_np_64 = trexio.read_basis_nucleus_index(test_file2)
 assert rindices_np_64.dtype is np.dtype(np.int64)
-assert rindices_np_64.size==nucleus_num
+assert rindices_np_64.size==basis_shell_num
-for i in range(nucleus_num):
+for i in range(basis_shell_num):
     assert rindices_np_64[i]==indices_np[i]
 
 # read nuclear coordinates without providing optional argument dim

python/test/test_pytrexio.py

@@ -59,11 +59,16 @@ assert rc==0
 #        charges[i] = 1.
 #rc = trexio_write_nucleus_charge(test_file, charges)
 
+basis_num = 24
 # initialize arrays of nuclear indices as a list and convert it to numpy array
-indices = [i for i in range(nucleus_num)]
+indices = [i for i in range(basis_num)]
 # type cast is important here because by default numpy transforms a list of integers into int64 array
 indices_np = np.array(indices, dtype=np.int32)
 
+# first write basis_num because it is needed to check dimensions of basis_nucleus_index in TREXIO >= 2.0.0
+rc = trexio_write_basis_shell_num(test_file, basis_num)
+assert rc==0
+
 # function call below works with both lists and numpy arrays, dimension needed for memory-safety is derived 
 # from the size of the list/array by SWIG using typemacs from numpy.i
 rc = trexio_write_safe_basis_nucleus_index(test_file, indices_np)
@@ -124,8 +129,12 @@ assert rc==23
 #for i in range(nucleus_num):
 #    assert charges2[i]==charges[i]
 
+result_basis = trexio_read_basis_shell_num(test_file2)
+assert result[0]==0
+assert result[1]==basis_num
+
 # safe call to read_safe array of int values
-rc, rindices_np = trexio_read_safe_basis_nucleus_index(test_file2, nucleus_num)
+rc, rindices_np = trexio_read_safe_basis_nucleus_index(test_file2, basis_num)
 assert rc==0
 assert rindices_np.dtype is np.dtype(np.int32)
 for i in range(nucleus_num):

tests/io_dset_int_hdf5.c

@@ -25,7 +25,7 @@ static int test_write_dset (const char* file_name, const back_end_t backend) {
   assert (file != NULL);
 
   // write numerical attribute in an empty file
-  rc = trexio_write_nucleus_num(file, num);
+  rc = trexio_write_basis_shell_num(file, num);
   assert (rc == TREXIO_SUCCESS);
   
   // write numerical (integer) dataset in a file
@@ -91,7 +91,7 @@ static int test_read_dset (const char* file_name, const back_end_t backend) {
   assert (file != NULL);
 
   // read numerical attribute from the file
-  rc = trexio_read_nucleus_num(file, &num);
+  rc = trexio_read_basis_shell_num(file, &num);
   assert (rc == TREXIO_SUCCESS);
   assert (num == 12);
 

tests/io_dset_int_text.c

@@ -25,7 +25,7 @@ static int test_write_dset (const char* file_name, const back_end_t backend) {
   assert (file != NULL);
 
   // write numerical attribute in an empty file
-  rc = trexio_write_nucleus_num(file, num);
+  rc = trexio_write_basis_shell_num(file, num);
   assert (rc == TREXIO_SUCCESS);
   
   // write numerical (integer) dataset in a file
@@ -91,7 +91,7 @@ static int test_read_dset (const char* file_name, const back_end_t backend) {
   assert (file != NULL);
 
   // read numerical attribute from the file
-  rc = trexio_read_nucleus_num(file, &num);
+  rc = trexio_read_basis_shell_num(file, &num);
   assert (rc == TREXIO_SUCCESS);
   assert (num == 12);
 

tests/test_f.f90

@@ -44,9 +44,9 @@ subroutine test_write(file_name, back_end)
 
   integer :: rc = 1
 
-  integer :: num
+  integer :: num, basis_shell_num
 
-  integer :: basis_nucleus_index(12)
+  integer :: basis_nucleus_index(24)
   double precision :: charge(12)
   double precision :: coord(3,12)
 
@@ -70,7 +70,8 @@ subroutine test_write(file_name, back_end)
                        0.00000000d0,  2.47304151d0 ,  0.00000000d0 /), &
                        shape(coord) )
 
-  basis_nucleus_index = (/ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 /)
+  basis_shell_num = 24
+  basis_nucleus_index = (/ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24 /)
 
   label = [character(len=8) :: 'C', 'Na','C', 'C 66', 'C','C', 'H 99', 'Ru', 'H', 'H',  'H', 'H' ]
 
@@ -104,6 +105,8 @@ subroutine test_write(file_name, back_end)
   deallocate(sym_str)
   call trexio_assert(rc, TREXIO_SUCCESS, 'SUCCESS WRITE POINT GROUP')
 
+  rc = trexio_write_basis_shell_num(trex_file, basis_shell_num)
+  call trexio_assert(rc, TREXIO_SUCCESS, 'SUCCESS WRITE BASIS NUM')
 
   rc = trexio_write_basis_nucleus_index(trex_file, basis_nucleus_index)
   call trexio_assert(rc, TREXIO_SUCCESS, 'SUCCESS WRITE INDEX')
@@ -137,9 +140,9 @@ subroutine test_read(file_name, back_end)
 
   integer :: i, j, k, ind, offset, flag
   integer :: rc = 1
-  integer :: num, num_read
+  integer :: num, num_read, basis_shell_num
 
-  integer :: basis_nucleus_index(12)
+  integer :: basis_nucleus_index(24)
   double precision :: charge(12)
   double precision :: coord(3,12)
 
@@ -152,6 +155,7 @@ subroutine test_read(file_name, back_end)
   character*(128) :: str
 
   num = 12
+  basis_shell_num = 24
 
 ! ================= START OF TEST ===================== !
 

trex.org

@@ -215,20 +215,18 @@ arrays are 0-based. Hence, we introduce the ~index~ type which is an
   All the basis set parameters are stored in one-dimensional arrays:
 
   #+NAME: basis
-  | Variable            | Type    | Dimensions         | Description                                              |
+  | Variable        | Type    | Dimensions          | Description                                                  |
-  |---------------------+---------+--------------------+----------------------------------------------------------|
+  |-----------------+---------+---------------------+--------------------------------------------------------------|
-  | ~type~              | ~str~   |                    | Type of basis set: "Gaussian" or "Slater"                |
+  | ~type~          | ~str~   |                     | Type of basis set: "Gaussian" or "Slater"                    |
-  | ~num~               | ~dim~   |                    | Total Number of shells                                   |
+  | ~prim_num~      | ~dim~   |                     | Total number of primitives                                   |
-  | ~prim_num~          | ~dim~   |                    | Total number of primitives                               |
+  | ~shell_num~     | ~dim~   |                     | Total number of shells                                       |
-  | ~nucleus_index~     | ~index~ | ~(nucleus.num)~    | Index of the first shell of each nucleus ($A$)           |
+  | ~nucleus_index~ | ~index~ | ~(basis.shell_num)~ | One-to-one correspondence between shells and atomic indices  |
-  | ~nucleus_shell_num~ | ~int~   | ~(nucleus.num)~    | Number of shells for each nucleus                        |
+  | ~ang_mom~       | ~int~   | ~(basis.shell_num)~ | One-to-one correspondence between shells and angular momenta |
-  | ~shell_ang_mom~     | ~int~   | ~(basis.num)~      | Angular momentum ~0:S, 1:P, 2:D, ...~                    |
+  | ~shell_factor~  | ~float~ | ~(basis.shell_num)~ | Normalization factor of each shell ($\mathcal{N}_s$)         |
-  | ~shell_prim_num~    | ~int~   | ~(basis.num)~      | Number of primitives in the shell ($N_{\text{prim}}$)    |
+  | ~shell_index~   | ~index~ | ~(basis.prim_num)~  | One-to-one correspondence between primitives and shell index |
-  | ~shell_factor~      | ~float~ | ~(basis.num)~      | Normalization factor of the shell ($\mathcal{N}_s$)      |
+  | ~exponent~      | ~float~ | ~(basis.prim_num)~  | Exponents of the primitives ($\gamma_{ks}$)                  |
-  | ~shell_prim_index~  | ~index~ | ~(basis.num)~      | Index of the first primitive in the complete list        |
+  | ~coefficient~   | ~float~ | ~(basis.prim_num)~  | Coefficients of the primitives ($a_{ks}$)                    |
-  | ~exponent~          | ~float~ | ~(basis.prim_num)~ | Exponents of the primitives ($\gamma_{ks}$)              |
+  | ~prim_factor~   | ~float~ | ~(basis.prim_num)~  | Normalization coefficients for the primitives ($f_{ks}$)     |
-  | ~coefficient~       | ~float~ | ~(basis.prim_num)~ | Coefficients of the primitives ($a_{ks}$)                |
-  | ~prim_factor~       | ~float~ | ~(basis.prim_num)~ | Normalization coefficients for the primitives ($f_{ks}$) |
 
   #+CALL: json(data=basis, title="basis")
 
@@ -236,22 +234,22 @@ arrays are 0-based. Hence, we introduce the ~index~ type which is an
   :RESULTS:
   #+begin_src python :tangle trex.json
       "basis": {
-		       "type" : [ "str"  , []                   ]
+		   "type" : [ "str"  , []                    ]
-	,               "num" : [ "dim"  , []                   ]
+	,      "prim_num" : [ "dim"  , []                    ]
-	,          "prim_num" : [ "dim"  , []                   ]
+	,     "shell_num" : [ "dim"  , []                    ]
-	,     "nucleus_index" : [ "index", [ "nucleus.num" ]    ]
+	, "nucleus_index" : [ "index", [ "basis.shell_num" ] ]
-	, "nucleus_shell_num" : [ "int"  , [ "nucleus.num" ]    ]
+	,       "ang_mom" : [ "int"  , [ "basis.shell_num" ] ]
-	,     "shell_ang_mom" : [ "int"  , [ "basis.num" ]      ]
+	,  "shell_factor" : [ "float", [ "basis.shell_num" ] ]
-	,    "shell_prim_num" : [ "int"  , [ "basis.num" ]      ]
+	,   "shell_index" : [ "index", [ "basis.prim_num" ]  ]
-	,      "shell_factor" : [ "float", [ "basis.num" ]      ]
+	,      "exponent" : [ "float", [ "basis.prim_num" ]  ]
-	,  "shell_prim_index" : [ "index", [ "basis.num" ]      ]
+	,   "coefficient" : [ "float", [ "basis.prim_num" ]  ]
-	,          "exponent" : [ "float", [ "basis.prim_num" ] ]
+	,   "prim_factor" : [ "float", [ "basis.prim_num" ]  ]
-	,       "coefficient" : [ "float", [ "basis.prim_num" ] ]
-	,       "prim_factor" : [ "float", [ "basis.prim_num" ] ]
       } ,
   #+end_src
   :END:
 
+** Example
+
   For example, consider H_2 with the following basis set (in GAMESS
   format), where both the AOs and primitives are considered normalized:
 
@@ -272,31 +270,44 @@ P   1
 P   1
 1         3.880000E-01           1.000000E+00
 D   1
-1         1.057000E+00           1.0000000
+1         1.057000E+00           1.000000E+00
    #+END_EXAMPLE
 
-   we have:
+   In TREXIO representaion we have:
 
    #+BEGIN_EXAMPLE
-type     = "Gaussian"
+type  = "Gaussian"
-num      = 12
+prim_num   = 20
-prim_num = 20
+shell_num   = 12
 
-nucleus_index    = [0 , 6]
+# 6 shells per H atom
-shell_ang_mom    = [0 , 0 , 0 , 1 , 1 , 2 , 0 , 0 , 0 , 1 , 1 , 2 ]
+nucleus_index = 
-shell_prim_num   = [5 , 1 , 1 , 1 , 1 , 1 , 5 , 1 , 1 , 1 , 1 , 1 ]
+[ 0, 0, 0, 0, 0, 0,  
-shell_prim_index = [0 , 5 , 6 , 7 , 8 , 9 , 10, 15, 16, 17, 18, 19]
+  1, 1, 1, 1, 1, 1 ] 
-shell_factor     = [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.]
 
+# 3 shells in S (l=0), 2 in P (l=1), 1 in D (l=2)
+ang_mom =
+[ 0, 0, 0, 1, 1, 2, 
+  0, 0, 0, 1, 1, 2 ]
+
+# no need to renormalize shells 
+shell_factor = 
+[ 1., 1., 1., 1., 1., 1.,  
+  1., 1., 1., 1., 1., 1. ] 
+
+# 5 primitives for the first S shell and then 1 primitive per remaining shells in each H atom
+shell_index = 
+[ 0, 0, 0, 0, 0, 1, 2, 3, 4, 5, 
+  6, 6, 6, 6, 6, 7, 8, 9, 10, 11 ]
+
+# parameters of the primitives (10 per H atom)
 exponent =
-[ 33.87, 5.095, 1.159, 0.3258, 0.1027, 0.3258, 0.1027, 1.407,
+[ 33.87, 5.095, 1.159, 0.3258, 0.1027, 0.3258, 0.1027, 1.407, 0.388, 1.057, 
-  0.388, 1.057, 33.87, 5.095, 1.159, 0.3258, 0.1027, 0.3258, 0.1027, 1.407,
+  33.87, 5.095, 1.159, 0.3258, 0.1027, 0.3258, 0.1027, 1.407, 0.388, 1.057 ]
-  0.388, 1.057]
 
 coefficient =
-[ 0.006068, 0.045308, 0.202822, 0.503903, 0.383421, 1.0, 1.0,
+[ 0.006068, 0.045308, 0.202822, 0.503903, 0.383421, 1.0, 1.0, 1.0, 1.0, 1.0, 
-  1.0, 1.0, 1.0, 0.006068, 0.045308, 0.202822, 0.503903, 0.383421, 1.0, 1.0,
+  0.006068, 0.045308, 0.202822, 0.503903, 0.383421, 1.0, 1.0, 1.0, 1.0, 1.0 ]
-  1.0, 1.0, 1.0]
 
 prim_factor =
 [ 1.0006253235944540e+01, 2.4169531573445120e+00, 7.9610924849766440e-01