diff --git a/.gitignore b/.gitignore index 288e294..0664d48 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,3 @@ *.o *. +*.so diff --git a/src/python/pyscf_interface.f90 b/src/python/pyscf_interface.f90 new file mode 100644 index 0000000..46788a2 --- /dev/null +++ b/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 + diff --git a/src/python/pyscf_module.py b/src/python/pyscf_module.py new file mode 100644 index 0000000..a89c8e6 --- /dev/null +++ b/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 + + diff --git a/src/python/pyscf_wrapper.c b/src/python/pyscf_wrapper.c new file mode 100644 index 0000000..fe6359c --- /dev/null +++ b/src/python/pyscf_wrapper.c @@ -0,0 +1,88 @@ +#include +#include + + +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(); +} + diff --git a/src/python/test_python_call.f90 b/src/python/test_python_call.f90 new file mode 100644 index 0000000..e73b130 --- /dev/null +++ b/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 + +