start pyscf interface

.gitignore

@@ -1,2 +1,3 @@
 *.o
 *.
+*.so

src/python/pyscf_interface.f90

@@ -0,0 +1,29 @@
+module pyscf_interface
+
+    use, intrinsic :: iso_c_binding
+
+    implicit none
+
+    interface
+
+        subroutine call_mol_prop(xyz, input_basis, charge, multiplicity, unit_type, cartesian, &
+                                 natoms, nalpha, nbeta, Enuc) bind(C, name="call_mol_prop")
+
+            import c_char, c_int, c_double, c_ptr
+
+            character(kind=c_char), intent(in) :: xyz(*)
+            character(kind=c_char), intent(in) :: input_basis(*)
+            integer(c_int), intent(in), value  :: charge
+            integer(c_int), intent(in), value  :: multiplicity
+            character(kind=c_char), intent(in) :: unit_type(*)
+            integer(c_int), intent(in), value  :: cartesian
+            integer(c_int), intent(out)        :: natoms
+            integer(c_int), intent(out)        :: nalpha
+            integer(c_int), intent(out)        :: nbeta
+            real(c_double), intent(out)        :: Enuc
+        end subroutine call_mol_prop
+
+    end interface
+
+end module pyscf_interface
+

src/python/pyscf_module.py

@@ -0,0 +1,52 @@
+
+import numpy as np
+from pyscf import gto
+
+
+def read_xyz(xyz):
+
+    f = open(xyz + '.xyz', 'r')
+
+    lines = f.read().splitlines()
+
+    # nb of atoms
+    nbAt = int(lines.pop(0))
+    lines.pop(0)
+
+    # list of atoms positions
+    list_pos_atom = []
+    for line in lines:
+        tmp = line.split()
+        atom = tmp[0]
+        pos = (float(tmp[1]), float(tmp[2]), float(tmp[3]))
+        list_pos_atom.append([atom,pos])
+
+    f.close()
+
+    return list_pos_atom
+
+
+def mol_prop(xyz: str, input_basis: str, charge: int, multiplicity: int, unit: str, cartesian: bool):
+
+    list_pos_atom = read_xyz(xyz)
+
+    mol = gto.M(
+        atom = list_pos_atom,
+        basis = input_basis,
+        charge = charge,
+        spin = multiplicity - 1
+    )
+
+    mol.unit = unit
+    mol.cart = cartesian
+
+    mol.build()
+
+    natoms = mol.natm       # nb of atoms
+    nalpha = mol.nelec[0]   # nb of alpha-electrons
+    nbeta = mol.nelec[1]    # nb of beta-electrons
+    Enuc = mol.energy_nuc() # nuclear energy
+
+    return natoms, nalpha, nbeta, Enuc
+
+

src/python/pyscf_wrapper.c

@@ -0,0 +1,88 @@
+#include <Python.h>
+#include <stdio.h>
+
+
+void call_mol_prop(const char *xyz, const char *input_basis, int charge, int multiplicity, 
+                   const char *unit, int cartesian, 
+                   int *natoms, int *nalpha, int *nbeta, double *Enuc) {
+
+    PyObject *pName, *pModule, *pFunc, *pArgs, *pValue;
+    PyObject *py_xyz, *py_input_basis, *py_unit, *py_cartesian;
+    PyObject *py_result;
+    
+    // Initialize the Python interpreter
+    Py_Initialize();
+
+    PyRun_SimpleString("import sys");
+    PyRun_SimpleString("sys.path.append('.')");
+
+    // Import the Python module
+    pName = PyUnicode_DecodeFSDefault("pyscf_module");
+    pModule = PyImport_Import(pName);
+    Py_XDECREF(pName);
+
+    if (pModule != NULL) {
+
+        // Get the Python function
+        pFunc = PyObject_GetAttrString(pModule, "mol_prop");
+
+        if (pFunc && PyCallable_Check(pFunc)) {
+
+            // Convert C strings to Python strings
+            py_xyz = PyUnicode_FromString(xyz);
+            py_input_basis = PyUnicode_FromString(input_basis);
+            py_unit = PyUnicode_FromString(unit);
+
+            // Convert C int to Python boolean
+            py_cartesian = PyBool_FromLong(cartesian);
+
+            // Create a tuple to hold the arguments for the Python function
+            pArgs = PyTuple_Pack(6, 
+                                 py_xyz, 
+                                 py_input_basis, 
+                                 PyLong_FromLong(charge), 
+                                 PyLong_FromLong(multiplicity), 
+                                 py_unit, 
+                                 py_cartesian);
+            
+            // Call the Python function
+            pValue = PyObject_CallObject(pFunc, pArgs);
+            Py_XDECREF(pArgs);
+
+            if (pValue != NULL) {
+                // Unpack the result tuple
+                if (PyTuple_Check(pValue) && PyTuple_Size(pValue) == 4) {
+                    *natoms = (int)PyLong_AsLong(PyTuple_GetItem(pValue, 0));
+                    *nalpha = (int)PyLong_AsLong(PyTuple_GetItem(pValue, 1));
+                    *nbeta = (int)PyLong_AsLong(PyTuple_GetItem(pValue, 2));
+                    *Enuc = PyFloat_AsDouble(PyTuple_GetItem(pValue, 3));
+                } else {
+                    PyErr_Print();
+                }
+                Py_XDECREF(pValue);  // Clean up reference to result tuple
+            } else {
+                // Handle error if Python function call fails
+                PyErr_Print();
+            }
+
+            // Clean up references
+            Py_XDECREF(py_xyz);
+            Py_XDECREF(py_input_basis);
+            Py_XDECREF(py_unit);
+            Py_XDECREF(py_cartesian);
+        } else {
+            // Handle error if Python function is not callable
+            PyErr_Print();
+        }
+
+        Py_XDECREF(pFunc);  // Clean up reference to Python function object
+        Py_XDECREF(pModule);  // Clean up reference to Python module object
+    } else {
+        // Handle error if Python module could not be imported
+        PyErr_Print();
+    }
+
+    // Finalize the Python interpreter
+    Py_Finalize();
+}
+

src/python/test_python_call.f90

@@ -0,0 +1,28 @@
+program test_python_call
+
+    use pyscf_interface
+    implicit none
+
+    character(len = 100) :: xyz, input_basis, unit_type
+    integer :: charge, multiplicity, cartesian
+    integer :: natoms, nalpha, nbeta
+    double precision :: Enuc
+
+    xyz = "../../mol/H2O"
+    input_basis = "sto-3g"
+    charge = 0
+    multiplicity = 1
+    unit_type = "Angstrom"
+    cartesian = 0
+
+    call call_mol_prop(trim(adjustl(xyz)), trim(adjustl(input_basis)), charge, multiplicity, trim(adjustl(unit_type)), cartesian, &
+                       natoms, nalpha, nbeta, Enuc)
+
+    print *, "Number of atoms:", natoms
+    print *, "Number of alpha electrons:", nalpha
+    print *, "Number of beta electrons:", nbeta
+    print *, "Nuclear energy:", Enuc
+
+end program test_python_call
+
+