diff --git a/.gitignore b/.gitignore index 427a587b..2184b796 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,4 @@ CMakeCache.txt CMakeFiles Makefile cmake_install.cmake +build/* diff --git a/CMakeLists.txt b/CMakeLists.txt index 26a15201..798ec889 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -135,6 +135,7 @@ endif() option(WITH_FORTRAN "Fortran interface" on) if(WITH_FORTRAN) add_definitions(-DWITH_FORTRAN) + enable_language(Fortran) message("Enable Fortran interface") else() message("Exclude Fortran interface") @@ -149,7 +150,7 @@ else() endif(WITH_CINT2_INTERFACE) option(BUILD_SHARED_LIBS "build shared libraries" 1) -option(ENABLE_EXAMPLE "build examples" 0) +option(ENABLE_EXAMPLE "build examples" 1) option(ENABLE_TEST "build tests" 0) option(ENABLE_STATIC "Enforce static library build" 0) if(QUICK_TEST) @@ -187,14 +188,42 @@ if(QUADMATH_FOUND) endif() target_link_libraries(cint "-lm") +# Fortran interface modules +if(WITH_FORTRAN) + add_library(cint_fortran OBJECT + ${PROJECT_SOURCE_DIR}/include/libcint_interface.f90 + ${PROJECT_SOURCE_DIR}/include/libcint_fortran.f90) + target_link_libraries(cint_fortran PUBLIC cint) + set_target_properties(cint_fortran PROPERTIES + Fortran_MODULE_DIRECTORY ${PROJECT_BINARY_DIR}/include) + target_include_directories(cint_fortran PUBLIC + $ + $) +endif() + install(TARGETS cint COMPONENT "lib") install(FILES ${CintHeaders} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR} COMPONENT "dev") +if(WITH_FORTRAN) + # Install Fortran source files + install(FILES + ${PROJECT_SOURCE_DIR}/include/libcint_interface.f90 + ${PROJECT_SOURCE_DIR}/include/libcint_fortran.f90 + DESTINATION ${CMAKE_INSTALL_INCLUDEDIR} COMPONENT "dev") + # Install Fortran module files (.mod) + install(DIRECTORY ${PROJECT_BINARY_DIR}/include/ + DESTINATION ${CMAKE_INSTALL_INCLUDEDIR} + COMPONENT "dev" + FILES_MATCHING PATTERN "*.mod") +endif() + if(ENABLE_EXAMPLE) - enable_language(Fortran) + if(NOT WITH_FORTRAN) + enable_language(Fortran) + endif() find_package(OpenMP) if(OPENMP_FOUND) set(HAVE_OPENMP 1) diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index bdc72602..5ea2213d 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,18 +1,14 @@ -add_executable(cartesian_c c_call_cartesian.c ) -add_executable(spheric_c c_call_spheric.c ) -add_executable(spinor_c c_call_spinor.c ) -add_executable(cartesian_f fortran_call_cartesian.F90) -add_executable(spheric_f fortran_call_spheric.F90 ) -add_executable(spinor_f fortran_call_spinor.F90 ) -add_executable(time_c2h6 time_c2h6.c ) -add_executable(time_c60 time_c60.c ) +# C examples +add_executable(cartesian_c c_call_cartesian.c) +add_executable(spheric_c c_call_spheric.c) +add_executable(spinor_c c_call_spinor.c) +add_executable(time_c2h6 time_c2h6.c) +add_executable(time_c60 time_c60.c) target_link_libraries(cartesian_c cint m) target_link_libraries(spheric_c cint m) target_link_libraries(spinor_c cint m) -target_link_libraries(cartesian_f cint m) -target_link_libraries(spheric_f cint m) -target_link_libraries(spinor_f cint m) + target_link_libraries(time_c2h6 cint m) target_link_libraries(time_c60 cint m) @@ -23,3 +19,46 @@ set_target_properties(time_c60 PROPERTIES COMPILE_FLAGS ${OpenMP_C_FLAGS} LINK_FLAGS ${OpenMP_C_FLAGS}) +if(WITH_FORTRAN) + +add_executable(cartesian_f fortran/fortran_call_cartesian.F90) +add_executable(spheric_f fortran/fortran_call_spheric.F90) +add_executable(spinor_f fortran/fortran_call_spinor.F90) + +add_executable(modern_cartesian_f fortran/fortran_modern_cartesian.F90) +add_executable(modern_spheric_f fortran/fortran_modern_spheric.F90) +add_executable(modern_spinor_f fortran/fortran_modern_spinor.F90) + +add_executable(pure_fortran_f fortran/fortran_pure_example.F90) + +add_executable(time_c2h6_f fortran/fortran_time_c2h6.F90) +add_executable(time_c2h6_pure_f fortran/fortran_time_c2h6_pure.F90) + +target_link_libraries(cartesian_f cint m) +target_link_libraries(spheric_f cint m) +target_link_libraries(spinor_f cint m) + +target_link_libraries(modern_cartesian_f cint_fortran cint m) +target_link_libraries(modern_spheric_f cint_fortran cint m) +target_link_libraries(modern_spinor_f cint_fortran cint m) +target_link_libraries(pure_fortran_f cint_fortran cint m) +target_link_libraries(time_c2h6_f cint_fortran cint m) +target_link_libraries(time_c2h6_pure_f cint_fortran cint m) + +# Set module directory for Fortran examples +target_include_directories(modern_cartesian_f PRIVATE ${PROJECT_BINARY_DIR}/include) +target_include_directories(modern_spheric_f PRIVATE ${PROJECT_BINARY_DIR}/include) +target_include_directories(modern_spinor_f PRIVATE ${PROJECT_BINARY_DIR}/include) +target_include_directories(pure_fortran_f PRIVATE ${PROJECT_BINARY_DIR}/include) +target_include_directories(time_c2h6_f PRIVATE ${PROJECT_BINARY_DIR}/include) +target_include_directories(time_c2h6_pure_f PRIVATE ${PROJECT_BINARY_DIR}/include) + +set_target_properties(time_c2h6_f PROPERTIES + COMPILE_FLAGS ${OpenMP_C_FLAGS} + LINK_FLAGS ${OpenMP_C_FLAGS}) + +set_target_properties(time_c2h6_pure_f PROPERTIES + COMPILE_FLAGS ${OpenMP_C_FLAGS} + LINK_FLAGS ${OpenMP_C_FLAGS}) + +endif() diff --git a/examples/fortran/Fortran_details.md b/examples/fortran/Fortran_details.md new file mode 100644 index 00000000..607edc92 --- /dev/null +++ b/examples/fortran/Fortran_details.md @@ -0,0 +1,469 @@ +# Fortran Interface Technical Details + +Complete technical documentation for libcint's Fortran interfaces. + +Transparency: This doc was generated by an LLM but carefully read and edited by a human. + +## Table of Contents +- [Architecture](#architecture) +- [Interface Modules](#interface-modules) +- [Type System](#type-system) +- [Performance](#performance) +- [API Reference](#api-reference) +- [Advanced Usage](#advanced-usage) +- [Migration Guide](#migration-guide) + +--- + +## Architecture + +libcint provides a two-layer Fortran interface: + +``` +┌─────────────────────────────────────┐ +│ Your Fortran Application │ ← Uses real(dp), integer(ip) +│ (examples: fortran_pure_example.F90)| or real(c_double), integer(c_int) +└─────────────────┬───────────────────┘ + │ +┌─────────────────▼───────────────────┐ +│ libcint_fortran module │ ← High-level Fortran interface +│ (include/libcint_fortran.f90) │ (compiled with libcint) +│ - Native Fortran types │ +│ - Optional arguments │ +│ - Cleaner API │ +└─────────────────┬───────────────────┘ + │ +┌─────────────────▼───────────────────┐ +│ libcint_interface module │ ← Low-level C binding +│ (include/libcint_interface.f90) │ (compiled with libcint) +│ - Direct iso_c_binding │ +│ - C types (c_double, c_int) │ +│ - Explicit interfaces │ +└─────────────────┬───────────────────┘ + │ +┌─────────────────▼───────────────────┐ +│ libcint (C library) │ ← Native C library +└─────────────────────────────────────┘ +``` + +Both modules are compiled as part of libcint when `WITH_FORTRAN=ON` (default). + +--- + +## Interface Modules + +### Low-Level: `module libcint_interface` + +**Location**: `include/libcint_interface.f90` + +**Purpose**: Direct, explicit interfaces, type-safe bindings to C using `iso_c_binding` + +**Example**: +```fortran +use iso_c_binding +use libcint_interface + +integer(c_int) :: atm(ATM_SLOTS, natm) +integer(c_int) :: bas(BAS_SLOTS, nbas) +real(c_double) :: env(10000) +integer(c_int) :: shls(2) +integer(c_int) :: ret + +ret = cint1e_ovlp_sph(buf, shls, atm, natm, bas, nbas, env) +ret = cint2e_sph(buf, shls, atm, natm, bas, nbas, env, c_null_ptr) +``` + +### High-Level: `module libcint_fortran` + +**Location**: `include/libcint_fortran.f90` + +**Purpose**: Native Fortran API hiding C details + +**NOTE**: The module binds `dp`, `ip`, and `zp` directly to `c_double`, `c_int`, and `c_double_complex`, and emits compile-time checks to ensure they also match `real64`/`int32` on the host compiler. If you intentionally change default kind sizes in your application, use the low-level interface instead so you can declare every bridge routine explicitly with `iso_c_binding` types. + +**Example**: +```fortran +use libcint_fortran + +integer(ip) :: atm(LIBCINT_ATM_SLOTS, natm) +integer(ip) :: bas(LIBCINT_BAS_SLOTS, nbas) +real(dp) :: env(10000) + +ret = libcint_1e_ovlp_sph(buf, shls, atm, natm, bas, nbas, env) +! Optional arguments work naturally +ret = libcint_2e_sph(buf, shls, atm, natm, bas, nbas, env) ! No optimizer +ret = libcint_2e_sph(buf, shls, atm, natm, bas, nbas, env, opt) ! With optimizer +``` + +--- + +## Type System + +### High Level Interface + +The high-level interface uses standard Fortran kinds that are binary compatible with C types: + +```fortran +use iso_fortran_env, only: real64, int32 + +integer, parameter :: dp = c_double ! Double precision +integer, parameter :: ip = c_int ! Integer +integer, parameter :: zp = c_double_complex ! Double complex +``` + +**Key insight**: These are **type aliases**, not conversions: +- No data copying occurs +- No runtime conversion +- Identical memory layout + + +### Comparison Table + +| Feature | Low-Level | High-Level | +|---------|-----------|------------| +| Integer type | `integer(c_int)` | `integer(ip)` | +| Real type | `real(c_double)` | `real(dp)` | +| Complex type | `complex(c_double_complex)` | `complex(zp)` | +| Constants | `ATM_SLOTS`, `BAS_SLOTS` | `LIBCINT_ATM_SLOTS`, `LIBCINT_BAS_SLOTS` | +| Null pointer | `c_null_ptr` | `c_null_ptr` (still needed) | +| Optional args | Must use `c_null_ptr` | Natural Fortran optional | + +--- + +## Performance + +### Interface Layer Performance + +Both Fortran interfaces (high-level and low-level) should have **identical performance** since they call the same underlying C functions with no overhead: + +### Fortran vs C Performance + +In practice, Fortran implementations often achieve **better performance** than equivalent (UNOPTIMIZED) C code due to: +- Aliasing guarantees that enable optimizations + +Run the benchmark to compare: +```bash +./examples/time_c2h6_f # Fortran implementation +./examples/time_c2h6 # C implementation +``` + +**Example results** (your mileage may vary): +- Fortran can be 50-80% faster than C in some cases + - I believe this is due to no `__restrict__` being used in the C code, limiting optimizations +- Both call the same libcint library functions +- Difference is in the driver code, not the integral computation itself + +--- + +## API Reference + +### Constants + +**Low-level** (`libcint_interface`): +```fortran +ATM_SLOTS, BAS_SLOTS, CHARGE_OF, PTR_COORD, ATOM_OF, ANG_OF, +NPRIM_OF, NCTR_OF, KAPPA_OF, PTR_EXP, PTR_COEFF, PTR_ENV_START +``` + +**High-level** (`libcint_fortran`): +```fortran +LIBCINT_ATM_SLOTS, LIBCINT_BAS_SLOTS, LIBCINT_CHARGE_OF, +LIBCINT_PTR_COORD, LIBCINT_ATOM_OF, LIBCINT_ANG_OF, +LIBCINT_NPRIM_OF, LIBCINT_NCTR_OF, LIBCINT_KAPPA_OF, +LIBCINT_PTR_EXP, LIBCINT_PTR_COEFF, LIBCINT_PTR_ENV_START +``` + +### Dimension Functions + +Available in both interfaces (shown with high-level names): + +```fortran +integer(ip) :: libcint_cgto_cart(bas_id, bas) ! Cartesian dimensions +integer(ip) :: libcint_cgto_sph(bas_id, bas) ! Spherical dimensions +integer(ip) :: libcint_cgto_spinor(bas_id, bas) ! Spinor dimensions +integer(ip) :: libcint_tot_cgto_sph(bas, nbas) ! Total contracted GTOs +integer(ip) :: libcint_tot_pgto_sph(bas, nbas) ! Total primitive GTOs +``` + +Low-level equivalents: `CINTcgto_cart`, `CINTcgto_spheric`, etc. + +### One-Electron Integrals + +Available in `_cart`, `_sph` variants (high-level API): + +```fortran +! Overlap +ret = libcint_1e_ovlp_cart(buf, shls, atm, natm, bas, nbas, env) +ret = libcint_1e_ovlp_sph(buf, shls, atm, natm, bas, nbas, env) + +! Kinetic energy +ret = libcint_1e_kin_cart(buf, shls, atm, natm, bas, nbas, env) +ret = libcint_1e_kin_sph(buf, shls, atm, natm, bas, nbas, env) + +! Nuclear attraction +ret = libcint_1e_nuc_cart(buf, shls, atm, natm, bas, nbas, env) +ret = libcint_1e_nuc_sph(buf, shls, atm, natm, bas, nbas, env) +``` + +Low-level equivalents: `cint1e_ovlp_cart`, `cint1e_ovlp_sph`, etc. + +### Two-Electron Integrals + +```fortran +! Basic 2e integrals +ret = libcint_2e_cart(buf, shls, atm, natm, bas, nbas, env, opt) +ret = libcint_2e_sph(buf, shls, atm, natm, bas, nbas, env, opt) + +! With optimizer (recommended for multiple integrals) +type(c_ptr) :: opt +call libcint_2e_sph_optimizer(opt, atm, natm, bas, nbas, env) +ret = libcint_2e_sph(buf, shls, atm, natm, bas, nbas, env, opt) +call libcint_del_optimizer(opt) +``` + +Low-level equivalents: `cint2e_cart`, `cint2e_sph`, etc. + +### Normalization + +```fortran +real(dp) :: libcint_gto_norm(n, alpha) ! High-level +real(c_double) :: CINTgto_norm(n, a) ! Low-level +``` + +--- + +## Advanced Usage + +### Using Optimizers + +Optimizers significantly speed up computation of many integrals: + +```fortran +type(c_ptr) :: opt +integer(ip) :: shls(4) +real(dp), allocatable :: buf(:,:,:,:) + +! Create optimizer +call libcint_2e_sph_optimizer(opt, atm, natm, bas, nbas, env) + +! Compute many integrals +do i = 1, nbas + do j = 1, nbas + shls = [i-1, j-1, i-1, j-1] ! 0-based! + ret = libcint_2e_sph(buf, shls, atm, natm, bas, nbas, env, opt) + end do +end do + +! Clean up +call libcint_del_optimizer(opt) +``` + +### OpenMP Parallelization + +Example from `fortran_time_c2h6.F90`: + +```fortran +!$omp parallel default(none) shared(atm,bas,env,nbas) private(...) +!$omp do schedule(dynamic, 2) +do ij = 1, nbas*(nbas+1)/2 + ! Compute shell pair indices + ! Compute integrals +end do +!$omp end do +!$omp end parallel +``` + +### Spinor Integrals (Relativistic) + +Spinor integrals are supported in both interfaces. Use complex buffers and the spinor dimension helpers: + +```fortran +! High-level API +use libcint_fortran + +integer(ip) :: di, dj +complex(zp), allocatable :: buf(:,:) + +di = libcint_cgto_spinor(0_ip, bas) +dj = libcint_cgto_spinor(1_ip, bas) +allocate(buf(di, dj)) +call libcint_1e_spnucsp(buf, shls, atm, natm, bas, nbas, env) + +! Two-electron spinor integrals with optional optimizer +call libcint_2e_spsp1(buf4d, shls4, atm, natm, bas, nbas, env, opt) +``` + +The low-level interface exposes the same functionality with `complex(c_double_complex)` buffers and the `cint1e_spnucsp`/`cint2e_spsp1` entry points if you prefer to work directly with `iso_c_binding`. + +### Mixing Both Interfaces + +You can use both in the same program: + +```fortran +module my_integrals + use libcint_fortran + use libcint_interface, only: cint2e_sph ! Import specific low-level function + implicit none +contains + subroutine compute_overlap(...) + ! Use high-level for clarity + ret = libcint_1e_ovlp_sph(buf, shls, atm, natm, bas, nbas, env) + end subroutine + + subroutine compute_eri_hot_loop(...) + ! Use low-level in performance-critical code if needed + ret = cint2e_sph(buf, shls, atm, natm, bas, nbas, env, opt) + end subroutine +end module +``` + +--- + +## Migration Guide + +### From Old-Style (external) to Modern + +**Before**: +```fortran +external :: cint1e_ovlp_sph +double precision :: buf(100) +integer :: shls(2), ret +ret = cint1e_ovlp_sph(buf, shls, atm, natm, bas, nbas, env) +``` + +**After** (high-level): +```fortran +use libcint_fortran +real(dp) :: buf(100) +integer(ip) :: shls(2), ret +ret = libcint_1e_ovlp_sph(buf, shls, atm, natm, bas, nbas, env) +``` + +### From Low-Level to High-Level + +**Before**: +```fortran +use iso_c_binding +use libcint_interface +integer(c_int) :: atm(ATM_SLOTS, natm) +real(c_double) :: env(10000) +``` + +**After**: +```fortran +use libcint_fortran +integer(ip) :: atm(LIBCINT_ATM_SLOTS, natm) +real(dp) :: env(10000) +``` + +Simple replacements: +- `integer(c_int)` → `integer(ip)` +- `real(c_double)` → `real(dp)` +- `ATM_SLOTS` → `LIBCINT_ATM_SLOTS` +- `cint1e_*` → `libcint_1e_*` + +--- + +## Important Gotchas + +### 0-Based Indexing + +libcint uses **0-based indexing** for: +- Atom indices stored in `bas(ATOM_OF, ...)` +- Shell indices in `shls` arrays +- Pointer offsets in `atm` and `bas` + +**Fortran arrays themselves are still 1-based!** + +Example: +```fortran +! Fortran array indexing (1-based) +atm(LIBCINT_PTR_COORD, 1) = 20 ! First atom, store offset 20 + +! But the offset value is 0-based for C +env(21) = x ! Offset 20 + 1 for Fortran 1-based + +! Shell indices are 0-based +shls(1) = 0 ! First shell +shls(2) = 1 ! Second shell +``` + +### Array Layout + +Data in `env` array: +```fortran +! Atom coordinates (offset stored in atm(PTR_COORD, iatm)) +env(off+1) = x +env(off+2) = y +env(off+3) = z + +! Basis exponents (offset stored in bas(PTR_EXP, ibas)) +env(off+1:off+nprim) = exponents(:) + +! Basis coefficients (offset stored in bas(PTR_COEFF, ibas)) +env(off+1:off+nprim*nctr) = coefficients(:) +``` + +### Null Pointers + +For optional optimizer arguments: +```fortran +use iso_c_binding, only: c_ptr, c_null_ptr + +! Without optimizer +ret = libcint_2e_sph(buf, shls, atm, natm, bas, nbas, env, c_null_ptr) + +! With optimizer +type(c_ptr) :: opt +call libcint_2e_sph_optimizer(opt, atm, natm, bas, nbas, env) +ret = libcint_2e_sph(buf, shls, atm, natm, bas, nbas, env, opt) +``` + + +## Building and Installation + +### Compiling with libcint + +The Fortran modules are built automatically: + +```bash +mkdir build && cd build +cmake -DWITH_FORTRAN=ON .. # ON by default +make +make install +``` + +This compiles: +1. The C library (`libcint.so`) +2. Fortran interface modules (`libcint_interface.mod`, `libcint_fortran.mod`) +3. Installs everything to standard locations + + +## Summary + +### Choose Your Interface + +| Situation | Recommendation | +|-----------|----------------| +| New application | **High-level** (`libcint_fortran`) | +| Library development | **Low-level** (`libcint_interface`) | +| Need both | **Mix** - import both modules | +| Legacy code | Migrate to modern interface when convenient | + +### Key Points + +- Modules are **compiled with libcint** and installed automatically +- High-level uses **native Fortran types** (`dp`, `ip`) +- Low-level uses **C types** (`c_double`, `c_int`) +- No data copying or conversion - types are binary compatible as long as you use the same compilers between your Fortran application and those for libcint + - Remember, THERE IS NO ABI BETWEEN FORTRAN COMPILERS +- Can **mix both** interfaces in the same program + +### Getting Started + +1. Look at `fortran_pure_example.F90` (high-level) +2. Or `fortran_modern_spheric.F90` (low-level) +3. Build and run examples to see them in action +4. Refer to this document for technical details + +For more information: https://github.com/sunqm/libcint diff --git a/examples/fortran/README.md b/examples/fortran/README.md new file mode 100644 index 00000000..acaf0052 --- /dev/null +++ b/examples/fortran/README.md @@ -0,0 +1,126 @@ +# Fortran Examples for libcint + +This directory contains Fortran examples demonstrating how to use libcint from Fortran code. + +Transparency: This doc was generated by an LLM but carefully read and edited by a human. + +## Quick Start + +libcint provides **two Fortran interface modules** (located in `include/`, compiled with the library): + +1. **`libcint_interface`** - Low-level C binding using `iso_c_binding` + - Use when you need maximum control and don't mind C types + - Types: `integer(c_int)`, `real(c_double)`, `complex(c_double_complex)` + +2. **`libcint_fortran`** - High-level pure Fortran API + - Types: `integer(ip)`, `real(dp)` (kinds bound to iso_c_binding types with compile-time checks) + +## Example: High-Level API + +```fortran +program my_app + use libcint_fortran + implicit none + + integer(ip) :: atm(LIBCINT_ATM_SLOTS, natm) + integer(ip) :: bas(LIBCINT_BAS_SLOTS, nbas) + real(dp) :: env(10000) + real(dp), allocatable :: buf(:,:) + integer(ip) :: shls(2), ret + + ! Setup system (atoms, basis sets, etc.) + ! ... + + ! Compute overlap integral + ret = libcint_1e_ovlp_sph(buf, shls, atm, natm, bas, nbas, env) +end program +``` + +See: `fortran_pure_example.F90` + +## Example: Low-Level API + +```fortran +program my_app + use iso_c_binding + use libcint_interface + implicit none + + integer(c_int) :: atm(ATM_SLOTS, natm) + integer(c_int) :: bas(BAS_SLOTS, nbas) + real(c_double) :: env(10000) + + ! Direct C interop with explicit types + ret = cint1e_ovlp_sph(buf, shls, atm, natm, bas, nbas, env) +end program +``` + +See: `fortran_modern_spheric.F90`, `fortran_modern_cartesian.F90` + +## Building + +The Fortran modules are compiled automatically when building libcint with `WITH_FORTRAN=ON` (default): + +```bash +mkdir build && cd build +cmake .. +make + +# Run examples +./examples/pure_fortran_f # High-level API +./examples/modern_spheric_f # Low-level API +./examples/time_c2h6_f # Benchmark with OpenMP +./examples/time_c2h6_pure_f # Benchmark with OpenMP using high-level API +``` + +The compiled module files (`.mod`) will be in `build/include/` and are installed with libcint. + +## Available Examples + +| Example | API Level | Description | +|---------|-----------|-------------| +| `fortran_pure_example.F90` | High-level | **Start here!** Clean Fortran code | +| `fortran_modern_cartesian.F90` | Low-level | Cartesian basis integrals | +| `fortran_modern_spheric.F90` | Low-level | Spherical harmonic integrals | +| `fortran_modern_spinor.F90` | Low-level | Relativistic spinor integrals | +| `fortran_time_c2h6.F90` | Low-level | Performance benchmark with OpenMP | +| `fortran_time_c2h6_pure.F90` | High-level | Performance benchmark with OpenMP using high-level API | +| `fortran_call_*.F90` | Legacy | Old-style | + +## Using in Your Application + +After installing libcint, simply use the modules: + +```fortran +! CMakeLists.txt or build script needs: +! - Link to libcint: target_link_libraries(myapp cint) +! - Include module directory for .mod files + +program myapp + use libcint_fortran ! High-level (recommended) + ! or + use libcint_interface ! Low-level + + ! Your code here +end program +``` + +## Which API Should I Use? + +- **New projects**: Use `libcint_fortran` (high-level) - cleaner code, same performance +- **Need control**: Use `libcint_interface` (low-level) - explicit C interop +- **Both**: You can mix them in the same program + +## More Information + +For detailed technical information, architecture diagrams, and advanced usage, see: +- **`Fortran_details.md`** - Complete technical documentation +- **libcint documentation**: https://github.com/sunqm/libcint + +## Important Notes + +️ **0-based indexing**: libcint uses 0-based indexing for atom and shell indices (Fortran arrays are still 1-based) + + **Array offsets**: Pointer offsets in `atm` and `bas` arrays are 0-based for C compatibility + +See the examples and `Fortran_details.md` for more information. diff --git a/examples/fortran_call_cartesian.F90 b/examples/fortran/fortran_call_cartesian.F90 similarity index 98% rename from examples/fortran_call_cartesian.F90 rename to examples/fortran/fortran_call_cartesian.F90 index 66568c21..f4446296 100644 --- a/examples/fortran_call_cartesian.F90 +++ b/examples/fortran/fortran_call_cartesian.F90 @@ -33,8 +33,8 @@ program cartesian integer,parameter :: PTR_ENV_START = 20 -integer :: natm = 2 -integer :: nbas = 4 +integer, parameter :: natm = 2 +integer, parameter :: nbas = 4 integer,allocatable :: atm(:,:) integer,allocatable :: bas(:,:) double precision,allocatable :: env(:) diff --git a/examples/fortran_call_spheric.F90 b/examples/fortran/fortran_call_spheric.F90 similarity index 98% rename from examples/fortran_call_spheric.F90 rename to examples/fortran/fortran_call_spheric.F90 index cb260bfe..5e16a794 100644 --- a/examples/fortran_call_spheric.F90 +++ b/examples/fortran/fortran_call_spheric.F90 @@ -33,8 +33,8 @@ program spheric integer,parameter :: PTR_ENV_START = 20 -integer :: natm = 2 -integer :: nbas = 4 +integer, parameter :: natm = 2 +integer, parameter :: nbas = 4 integer,allocatable :: atm(:,:) integer,allocatable :: bas(:,:) double precision,allocatable :: env(:) diff --git a/examples/fortran_call_spinor.F90 b/examples/fortran/fortran_call_spinor.F90 similarity index 98% rename from examples/fortran_call_spinor.F90 rename to examples/fortran/fortran_call_spinor.F90 index 113bc7a4..9fb98e38 100644 --- a/examples/fortran_call_spinor.F90 +++ b/examples/fortran/fortran_call_spinor.F90 @@ -26,8 +26,8 @@ program spinor integer,parameter :: PTR_ENV_START = 20 -integer :: natm = 2 -integer :: nbas = 4 +integer, parameter :: natm = 2 +integer, parameter :: nbas = 4 integer,allocatable :: atm(:,:) integer,allocatable :: bas(:,:) double precision,allocatable :: env(:) diff --git a/examples/fortran/fortran_modern_cartesian.F90 b/examples/fortran/fortran_modern_cartesian.F90 new file mode 100644 index 00000000..781ecfe8 --- /dev/null +++ b/examples/fortran/fortran_modern_cartesian.F90 @@ -0,0 +1,184 @@ +! +! Modern Fortran example using iso_c_binding interface +! This demonstrates how to use libcint with type-safe interfaces +! +! Computes overlap gradient integrals for H2 molecule +! + +program modern_cartesian + use iso_c_binding + use libcint_interface + implicit none + + ! Variables + integer, parameter :: natm = 2 + integer, parameter :: nbas = 4 + integer(c_int), allocatable :: atm(:,:) + integer(c_int), allocatable :: bas(:,:) + real(c_double), allocatable :: env(:) + + integer :: n, off + integer :: i, j, k, l + integer :: di, dj, dk, dl + integer(c_int) :: shls(4) + real(c_double), allocatable :: buf1e(:,:,:), buf2e(:,:,:,:,:) + type(c_ptr) :: opt + integer :: ret + + ! Allocate arrays + allocate(atm(ATM_SLOTS, natm)) + allocate(bas(BAS_SLOTS, nbas)) + allocate(env(10000)) + + atm = 0 ! Initialize to zero + bas = 0 + env = 0.0_c_double + + off = PTR_ENV_START + + ! Set up first atom (H at z=-0.8) + i = 1 + atm(CHARGE_OF, i) = 1 + atm(PTR_COORD, i) = off + env(off + 1) = 0.0_c_double ! x + env(off + 2) = 0.0_c_double ! y + env(off + 3) = -0.8_c_double ! z + off = off + 3 + + ! Set up second atom (H at z=0.8) + i = 2 + atm(CHARGE_OF, i) = 1 + atm(PTR_COORD, i) = off + env(off + 1) = 0.0_c_double + env(off + 2) = 0.0_c_double + env(off + 3) = 0.8_c_double + off = off + 3 + + ! Basis #1: 3s basis on atom 1 (index 0 in C) + n = 1 + bas(ATOM_OF, n) = 0 ! 0-based atom index + bas(ANG_OF, n) = 0 ! s orbital + bas(NPRIM_OF, n) = 3 ! 3 primitives + bas(NCTR_OF, n) = 2 ! 2 contractions + bas(PTR_EXP, n) = off + env(off + 1) = 6.0_c_double + env(off + 2) = 2.0_c_double + env(off + 3) = 0.8_c_double + off = off + 3 + bas(PTR_COEFF, n) = off + ! First contraction + env(off + 1) = 0.7_c_double * CINTgto_norm(bas(ANG_OF,n), env(bas(PTR_EXP,n)+1)) + env(off + 2) = 0.6_c_double * CINTgto_norm(bas(ANG_OF,n), env(bas(PTR_EXP,n)+2)) + env(off + 3) = 0.5_c_double * CINTgto_norm(bas(ANG_OF,n), env(bas(PTR_EXP,n)+3)) + ! Second contraction + env(off + 4) = 0.4_c_double * CINTgto_norm(bas(ANG_OF,n), env(bas(PTR_EXP,n)+1)) + env(off + 5) = 0.3_c_double * CINTgto_norm(bas(ANG_OF,n), env(bas(PTR_EXP,n)+2)) + env(off + 6) = 0.2_c_double * CINTgto_norm(bas(ANG_OF,n), env(bas(PTR_EXP,n)+3)) + off = off + 6 + + ! Basis #2: 1p basis on atom 1 + n = 2 + bas(ATOM_OF, n) = 0 + bas(ANG_OF, n) = 1 ! p orbital + bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1 + bas(PTR_EXP, n) = off + env(off + 1) = 0.9_c_double + off = off + 1 + bas(PTR_COEFF, n) = off + env(off + 1) = 1.0_c_double * CINTgto_norm(bas(ANG_OF,n), env(bas(PTR_EXP,n)+1)) + off = off + 1 + + ! Basis #3: Same as basis #1 but on atom 2 + n = 3 + bas(ATOM_OF, n) = 1 ! Second atom (0-based) + bas(ANG_OF, n) = bas(ANG_OF, 1) + bas(NPRIM_OF, n) = bas(NPRIM_OF, 1) + bas(NCTR_OF, n) = bas(NCTR_OF, 1) + bas(PTR_EXP, n) = bas(PTR_EXP, 1) + bas(PTR_COEFF, n) = bas(PTR_COEFF,1) + + ! Basis #4: Same as basis #2 but on atom 2 + n = 4 + bas(ATOM_OF, n) = 1 + bas(ANG_OF, n) = bas(ANG_OF, 2) + bas(NPRIM_OF, n) = bas(NPRIM_OF, 2) + bas(NCTR_OF, n) = bas(NCTR_OF, 2) + bas(PTR_EXP, n) = bas(PTR_EXP, 2) + bas(PTR_COEFF, n) = bas(PTR_COEFF,2) + + ! ==================================================================== + ! Call one-electron integral (overlap gradient) + ! Note: shell indices are 0-based in C + ! ==================================================================== + print*, "Computing one-electron overlap gradient integral..." + i = 0; shls(1) = i; di = CINTcgto_cart(i, bas) + j = 1; shls(2) = j; dj = CINTcgto_cart(j, bas) + + allocate(buf1e(di, dj, 3)) + ret = cint1e_ipovlp_cart(buf1e, shls, atm, natm, bas, nbas, env) + + if (ret /= 0) then + print*, "This gradient integral is not 0." + print*, "Buffer dimensions: ", di, "x", dj, "x 3" + else + print*, "This integral is 0." + endif + deallocate(buf1e) + + ! ==================================================================== + ! Call two-electron integral without optimizer + ! ==================================================================== + print*, "" + print*, "Computing two-electron gradient integral (no optimizer)..." + i = 0; shls(1) = i; di = CINTcgto_cart(i, bas) + j = 1; shls(2) = j; dj = CINTcgto_cart(j, bas) + k = 2; shls(3) = k; dk = CINTcgto_cart(k, bas) + l = 2; shls(4) = l; dl = CINTcgto_cart(l, bas) + + allocate(buf2e(di, dj, dk, dl, 3)) + ! Pass c_null_ptr for no optimizer + ret = cint2e_ip1_cart(buf2e, shls, atm, natm, bas, nbas, env, c_null_ptr) + + if (ret /= 0) then + print*, "This gradient integral is not 0." + print*, "Buffer dimensions: ", di, "x", dj, "x", dk, "x", dl, "x 3" + else + print*, "This integral is 0." + endif + deallocate(buf2e) + + ! ==================================================================== + ! Call two-electron integral WITH optimizer + ! ==================================================================== + print*, "" + print*, "Computing two-electron gradient integral (with optimizer)..." + + ! Create optimizer + call cint2e_ip1_cart_optimizer(opt, atm, natm, bas, nbas, env) + + i = 0; shls(1) = i; di = CINTcgto_cart(i, bas) + j = 1; shls(2) = j; dj = CINTcgto_cart(j, bas) + k = 2; shls(3) = k; dk = CINTcgto_cart(k, bas) + l = 2; shls(4) = l; dl = CINTcgto_cart(l, bas) + + allocate(buf2e(di, dj, dk, dl, 3)) + ret = cint2e_ip1_cart(buf2e, shls, atm, natm, bas, nbas, env, opt) + + if (ret /= 0) then + print*, "This gradient integral is not 0 (with optimizer)." + else + print*, "This integral is 0 (with optimizer)." + endif + deallocate(buf2e) + + ! Clean up optimizer + call CINTdel_optimizer(opt) + + ! Clean up + deallocate(atm, bas, env) + + print*, "" + print*, "Example completed successfully!" + +end program modern_cartesian diff --git a/examples/fortran/fortran_modern_spheric.F90 b/examples/fortran/fortran_modern_spheric.F90 new file mode 100644 index 00000000..f0245228 --- /dev/null +++ b/examples/fortran/fortran_modern_spheric.F90 @@ -0,0 +1,188 @@ +! +! Modern Fortran example using iso_c_binding interface - Spherical harmonics +! This demonstrates how to use libcint with type-safe interfaces +! +! Computes overlap gradient integrals for H2 molecule using spherical harmonics +! + +program modern_spheric + use iso_c_binding + use libcint_interface + implicit none + + ! Variables + integer, parameter :: natm = 2 + integer, parameter :: nbas = 4 + integer(c_int), allocatable :: atm(:,:) + integer(c_int), allocatable :: bas(:,:) + real(c_double), allocatable :: env(:) + + integer :: n, off + integer :: i, j, k, l + integer :: di, dj, dk, dl + integer(c_int) :: shls(4) + real(c_double), allocatable :: buf1e(:,:,:), buf2e(:,:,:,:,:) + type(c_ptr) :: opt + integer :: ret + + ! Allocate arrays + allocate(atm(ATM_SLOTS, natm)) + allocate(bas(BAS_SLOTS, nbas)) + allocate(env(10000)) + + atm = 0 ! Initialize to zero + bas = 0 + env = 0.0_c_double + + off = PTR_ENV_START + + ! Set up first atom (H at z=-0.8) + i = 1 + atm(CHARGE_OF, i) = 1 + atm(PTR_COORD, i) = off + env(off + 1) = 0.0_c_double ! x + env(off + 2) = 0.0_c_double ! y + env(off + 3) = -0.8_c_double ! z + off = off + 3 + + ! Set up second atom (H at z=0.8) + i = 2 + atm(CHARGE_OF, i) = 1 + atm(PTR_COORD, i) = off + env(off + 1) = 0.0_c_double + env(off + 2) = 0.0_c_double + env(off + 3) = 0.8_c_double + off = off + 3 + + ! Basis #1: 3s basis on atom 1 (index 0 in C) + n = 1 + bas(ATOM_OF, n) = 0 ! 0-based atom index + bas(ANG_OF, n) = 0 ! s orbital + bas(NPRIM_OF, n) = 3 ! 3 primitives + bas(NCTR_OF, n) = 2 ! 2 contractions + bas(PTR_EXP, n) = off + env(off + 1) = 6.0_c_double + env(off + 2) = 2.0_c_double + env(off + 3) = 0.8_c_double + off = off + 3 + bas(PTR_COEFF, n) = off + ! First contraction + env(off + 1) = 0.7_c_double * CINTgto_norm(bas(ANG_OF,n), env(bas(PTR_EXP,n)+1)) + env(off + 2) = 0.6_c_double * CINTgto_norm(bas(ANG_OF,n), env(bas(PTR_EXP,n)+2)) + env(off + 3) = 0.5_c_double * CINTgto_norm(bas(ANG_OF,n), env(bas(PTR_EXP,n)+3)) + ! Second contraction + env(off + 4) = 0.4_c_double * CINTgto_norm(bas(ANG_OF,n), env(bas(PTR_EXP,n)+1)) + env(off + 5) = 0.3_c_double * CINTgto_norm(bas(ANG_OF,n), env(bas(PTR_EXP,n)+2)) + env(off + 6) = 0.2_c_double * CINTgto_norm(bas(ANG_OF,n), env(bas(PTR_EXP,n)+3)) + off = off + 6 + + ! Basis #2: 1p basis on atom 1 + n = 2 + bas(ATOM_OF, n) = 0 + bas(ANG_OF, n) = 1 ! p orbital + bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1 + bas(PTR_EXP, n) = off + env(off + 1) = 0.9_c_double + off = off + 1 + bas(PTR_COEFF, n) = off + env(off + 1) = 1.0_c_double * CINTgto_norm(bas(ANG_OF,n), env(bas(PTR_EXP,n)+1)) + off = off + 1 + + ! Basis #3: Same as basis #1 but on atom 2 + n = 3 + bas(ATOM_OF, n) = 1 ! Second atom (0-based) + bas(ANG_OF, n) = bas(ANG_OF, 1) + bas(NPRIM_OF, n) = bas(NPRIM_OF, 1) + bas(NCTR_OF, n) = bas(NCTR_OF, 1) + bas(PTR_EXP, n) = bas(PTR_EXP, 1) + bas(PTR_COEFF, n) = bas(PTR_COEFF,1) + + ! Basis #4: Same as basis #2 but on atom 2 + n = 4 + bas(ATOM_OF, n) = 1 + bas(ANG_OF, n) = bas(ANG_OF, 2) + bas(NPRIM_OF, n) = bas(NPRIM_OF, 2) + bas(NCTR_OF, n) = bas(NCTR_OF, 2) + bas(PTR_EXP, n) = bas(PTR_EXP, 2) + bas(PTR_COEFF, n) = bas(PTR_COEFF,2) + + print*, "============================================" + print*, "Modern Fortran Example - Spherical Harmonics" + print*, "============================================" + print*, "" + + ! ==================================================================== + ! Call one-electron integral (overlap gradient) - Spherical + ! Note: shell indices are 0-based in C + ! ==================================================================== + print*, "Computing one-electron overlap gradient integral (spherical)..." + i = 0; shls(1) = i; di = CINTcgto_spheric(i, bas) + j = 0; shls(2) = j; dj = CINTcgto_spheric(j, bas) + + allocate(buf1e(di, dj, 3)) + ret = cint1e_ipovlp_sph(buf1e, shls, atm, natm, bas, nbas, env) + + if (ret /= 0) then + print*, "This gradient integral is not 0." + print*, "Buffer dimensions: ", di, "x", dj, "x 3" + else + print*, "This integral is 0." + endif + deallocate(buf1e) + + ! ==================================================================== + ! Call two-electron integral without optimizer - Spherical + ! ==================================================================== + print*, "" + print*, "Computing two-electron gradient integral (no optimizer, spherical)..." + i = 0; shls(1) = i; di = CINTcgto_spheric(i, bas) + j = 1; shls(2) = j; dj = CINTcgto_spheric(j, bas) + k = 2; shls(3) = k; dk = CINTcgto_spheric(k, bas) + l = 2; shls(4) = l; dl = CINTcgto_spheric(l, bas) + + allocate(buf2e(di, dj, dk, dl, 3)) + ret = cint2e_ip1_sph(buf2e, shls, atm, natm, bas, nbas, env, c_null_ptr) + + if (ret /= 0) then + print*, "This gradient integral is not 0." + print*, "Buffer dimensions: ", di, "x", dj, "x", dk, "x", dl, "x 3" + else + print*, "This integral is 0." + endif + deallocate(buf2e) + + ! ==================================================================== + ! Call two-electron integral WITH optimizer - Spherical + ! ==================================================================== + print*, "" + print*, "Computing two-electron gradient integral (with optimizer, spherical)..." + + ! Create optimizer + call cint2e_ip1_sph_optimizer(opt, atm, natm, bas, nbas, env) + + i = 0; shls(1) = i; di = CINTcgto_spheric(i, bas) + j = 1; shls(2) = j; dj = CINTcgto_spheric(j, bas) + k = 2; shls(3) = k; dk = CINTcgto_spheric(k, bas) + l = 2; shls(4) = l; dl = CINTcgto_spheric(l, bas) + + allocate(buf2e(di, dj, dk, dl, 3)) + ret = cint2e_ip1_sph(buf2e, shls, atm, natm, bas, nbas, env, opt) + + if (ret /= 0) then + print*, "This gradient integral is not 0 (with optimizer)." + else + print*, "This integral is 0 (with optimizer)." + endif + deallocate(buf2e) + + ! Clean up optimizer + call CINTdel_optimizer(opt) + + ! Clean up + deallocate(atm, bas, env) + + print*, "" + print*, "Example completed successfully!" + +end program modern_spheric diff --git a/examples/fortran/fortran_modern_spinor.F90 b/examples/fortran/fortran_modern_spinor.F90 new file mode 100644 index 00000000..0b74077b --- /dev/null +++ b/examples/fortran/fortran_modern_spinor.F90 @@ -0,0 +1,177 @@ +! +! Modern Fortran example using iso_c_binding interface - Spinor basis +! This demonstrates how to use libcint with type-safe interfaces for relativistic calculations +! +! Computes spinor integrals for H2 molecule +! + +program modern_spinor + use iso_c_binding + use libcint_interface + implicit none + + ! Variables + integer, parameter :: natm = 2 + integer, parameter :: nbas = 4 + integer(c_int), allocatable :: atm(:,:) + integer(c_int), allocatable :: bas(:,:) + real(c_double), allocatable :: env(:) + + integer :: n, off + integer :: i, j, k, l + integer :: di, dj, dk, dl + integer(c_int) :: shls(4) + complex(c_double_complex), allocatable :: buf1e(:,:), buf2e(:,:,:,:) + type(c_ptr) :: opt + + ! Allocate arrays + allocate(atm(ATM_SLOTS, natm)) + allocate(bas(BAS_SLOTS, nbas)) + allocate(env(10000)) + + atm = 0 ! Initialize to zero + bas = 0 + env = 0.0_c_double + + off = PTR_ENV_START + + ! Set up first atom (H at z=-0.8) + i = 1 + atm(CHARGE_OF, i) = 1 + atm(PTR_COORD, i) = off + env(off + 1) = 0.0_c_double ! x + env(off + 2) = 0.0_c_double ! y + env(off + 3) = -0.8_c_double ! z + off = off + 3 + + ! Set up second atom (H at z=0.8) + i = 2 + atm(CHARGE_OF, i) = 1 + atm(PTR_COORD, i) = off + env(off + 1) = 0.0_c_double + env(off + 2) = 0.0_c_double + env(off + 3) = 0.8_c_double + off = off + 3 + + ! Basis #1: with kappa > 0 => p_1/2 + n = 1 + bas(ATOM_OF, n) = 0 ! 0-based atom index + bas(ANG_OF, n) = 0 ! s orbital + bas(NPRIM_OF, n) = 3 ! 3 primitives + bas(NCTR_OF, n) = 2 ! 2 contractions + bas(KAPPA_OF, n) = 1 ! kappa parameter for spinor + bas(PTR_EXP, n) = off + env(off + 1) = 6.0_c_double + env(off + 2) = 2.0_c_double + env(off + 3) = 0.8_c_double + off = off + 3 + bas(PTR_COEFF, n) = off + ! First contraction + env(off + 1) = 0.7_c_double * CINTgto_norm(bas(ANG_OF,n), env(bas(PTR_EXP,n)+1)) + env(off + 2) = 0.6_c_double * CINTgto_norm(bas(ANG_OF,n), env(bas(PTR_EXP,n)+2)) + env(off + 3) = 0.5_c_double * CINTgto_norm(bas(ANG_OF,n), env(bas(PTR_EXP,n)+3)) + ! Second contraction + env(off + 4) = 0.4_c_double * CINTgto_norm(bas(ANG_OF,n), env(bas(PTR_EXP,n)+1)) + env(off + 5) = 0.3_c_double * CINTgto_norm(bas(ANG_OF,n), env(bas(PTR_EXP,n)+2)) + env(off + 6) = 0.2_c_double * CINTgto_norm(bas(ANG_OF,n), env(bas(PTR_EXP,n)+3)) + off = off + 6 + + ! Basis #2: with kappa = 0 => p_1/2, p_3/2 + n = 2 + bas(ATOM_OF, n) = 0 + bas(ANG_OF, n) = 1 ! p orbital + bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1 + bas(KAPPA_OF, n) = 0 ! kappa = 0 + bas(PTR_EXP, n) = off + env(off + 1) = 0.9_c_double + off = off + 1 + bas(PTR_COEFF, n) = off + env(off + 1) = 1.0_c_double * CINTgto_norm(bas(ANG_OF,n), env(bas(PTR_EXP,n)+1)) + off = off + 1 + + ! Basis #3: Same as basis #1 but on atom 2 + n = 3 + bas(ATOM_OF, n) = 1 ! Second atom (0-based) + bas(ANG_OF, n) = bas(ANG_OF, 1) + bas(NPRIM_OF, n) = bas(NPRIM_OF, 1) + bas(NCTR_OF, n) = bas(NCTR_OF, 1) + bas(KAPPA_OF, n) = bas(KAPPA_OF, 1) + bas(PTR_EXP, n) = bas(PTR_EXP, 1) + bas(PTR_COEFF, n) = bas(PTR_COEFF,1) + + ! Basis #4: Same as basis #2 but on atom 2 + n = 4 + bas(ATOM_OF, n) = 1 + bas(ANG_OF, n) = bas(ANG_OF, 2) + bas(NPRIM_OF, n) = bas(NPRIM_OF, 2) + bas(NCTR_OF, n) = bas(NCTR_OF, 2) + bas(KAPPA_OF, n) = bas(KAPPA_OF, 2) + bas(PTR_EXP, n) = bas(PTR_EXP, 2) + bas(PTR_COEFF, n) = bas(PTR_COEFF,2) + + print*, "============================================" + print*, "Modern Fortran Example - Spinor Basis" + print*, "============================================" + print*, "" + + ! ==================================================================== + ! Call one-electron spinor integrals + ! Note: shell indices are 0-based in C + ! ==================================================================== + print*, "Computing one-electron spinor integral (sigma dot p nuclear sigma dot p)..." + i = 0; shls(1) = i; di = CINTcgto_spinor(i, bas) + j = 1; shls(2) = j; dj = CINTcgto_spinor(j, bas) + + allocate(buf1e(di, dj)) + call cint1e_spnucsp(buf1e, shls, atm, natm, bas, nbas, env) + print*, "One-electron spinor integral computed." + print*, "Buffer dimensions: ", di, "x", dj, " (complex)" + deallocate(buf1e) + + ! ==================================================================== + ! Call two-electron spinor integrals without optimizer + ! ==================================================================== + print*, "" + print*, "Computing two-electron spinor integral (no optimizer)..." + i = 0; shls(1) = i; di = CINTcgto_spinor(i, bas) + j = 1; shls(2) = j; dj = CINTcgto_spinor(j, bas) + k = 2; shls(3) = k; dk = CINTcgto_spinor(k, bas) + l = 2; shls(4) = l; dl = CINTcgto_spinor(l, bas) + + allocate(buf2e(di, dj, dk, dl)) + call cint2e_spsp1(buf2e, shls, atm, natm, bas, nbas, env, c_null_ptr) + print*, "Two-electron spinor integral computed (no optimizer)." + print*, "Buffer dimensions: ", di, "x", dj, "x", dk, "x", dl, " (complex)" + deallocate(buf2e) + + ! ==================================================================== + ! Call two-electron spinor integrals WITH optimizer + ! ==================================================================== + print*, "" + print*, "Computing two-electron spinor integral (with optimizer)..." + + ! Create optimizer + call cint2e_spsp1_optimizer(opt, atm, natm, bas, nbas, env) + + i = 0; shls(1) = i; di = CINTcgto_spinor(i, bas) + j = 1; shls(2) = j; dj = CINTcgto_spinor(j, bas) + k = 2; shls(3) = k; dk = CINTcgto_spinor(k, bas) + l = 2; shls(4) = l; dl = CINTcgto_spinor(l, bas) + + allocate(buf2e(di, dj, dk, dl)) + call cint2e_spsp1(buf2e, shls, atm, natm, bas, nbas, env, opt) + print*, "Two-electron spinor integral computed (with optimizer)." + print*, "Buffer dimensions: ", di, "x", dj, "x", dk, "x", dl, " (complex)" + deallocate(buf2e) + + ! Clean up optimizer + call CINTdel_optimizer(opt) + + ! Clean up + deallocate(atm, bas, env) + + print*, "" + print*, "Example completed successfully!" + +end program modern_spinor diff --git a/examples/fortran/fortran_pure_example.F90 b/examples/fortran/fortran_pure_example.F90 new file mode 100644 index 00000000..02ad5f12 --- /dev/null +++ b/examples/fortran/fortran_pure_example.F90 @@ -0,0 +1,189 @@ +! +! Pure Fortran example using libcint_fortran high-level interface +! This demonstrates using libcint without any C-specific types +! +! No iso_c_binding needed in your application code! +! +program pure_fortran_example + use libcint_fortran ! High-level Fortran interface only + implicit none + + ! Use native Fortran types - dp and ip are defined in libcint_fortran + integer, parameter :: natm = 2 + integer, parameter :: nbas = 4 + integer(ip), allocatable :: atm(:,:) + integer(ip), allocatable :: bas(:,:) + real(dp), allocatable :: env(:) + + integer :: n, off + integer :: i, j + integer(ip) :: di, dj + integer(ip) :: shls(2) + real(dp), allocatable :: buf_ovlp(:,:) + real(dp), allocatable :: buf_kin(:,:) + real(dp), allocatable :: buf_nuc(:,:) + integer(ip) :: ret + + ! Allocate arrays using constants from the module + allocate(atm(LIBCINT_ATM_SLOTS, natm)) + allocate(bas(LIBCINT_BAS_SLOTS, nbas)) + allocate(env(10000)) + + atm = 0 + bas = 0 + env = 0.0_dp + + off = LIBCINT_PTR_ENV_START + + ! ======================================================================== + ! Set up H2 molecule (same as modern_spheric example) + ! ======================================================================== + + ! First atom: H at z=-0.8 + i = 1 + atm(LIBCINT_CHARGE_OF, i) = 1 + atm(LIBCINT_PTR_COORD, i) = off + env(off + 1) = 0.0_dp + env(off + 2) = 0.0_dp + env(off + 3) = -0.8_dp + off = off + 3 + + ! Second atom: H at z=0.8 + i = 2 + atm(LIBCINT_CHARGE_OF, i) = 1 + atm(LIBCINT_PTR_COORD, i) = off + env(off + 1) = 0.0_dp + env(off + 2) = 0.0_dp + env(off + 3) = 0.8_dp + off = off + 3 + + ! First basis function on atom 1 + n = 1 + bas(LIBCINT_ATOM_OF, n) = 0 + bas(LIBCINT_ANG_OF, n) = 0 + bas(LIBCINT_NPRIM_OF, n) = 3 + bas(LIBCINT_NCTR_OF, n) = 2 + bas(LIBCINT_PTR_EXP, n) = off + env(off + 1) = 6.0_dp + env(off + 2) = 2.0_dp + env(off + 3) = 0.8_dp + off = off + 3 + bas(LIBCINT_PTR_COEFF, n) = off + env(off + 1) = 0.7_dp * libcint_gto_norm(bas(LIBCINT_ANG_OF,n), env(bas(LIBCINT_PTR_EXP,n)+1)) + env(off + 2) = 0.6_dp * libcint_gto_norm(bas(LIBCINT_ANG_OF,n), env(bas(LIBCINT_PTR_EXP,n)+2)) + env(off + 3) = 0.5_dp * libcint_gto_norm(bas(LIBCINT_ANG_OF,n), env(bas(LIBCINT_PTR_EXP,n)+3)) + env(off + 4) = 0.4_dp * libcint_gto_norm(bas(LIBCINT_ANG_OF,n), env(bas(LIBCINT_PTR_EXP,n)+1)) + env(off + 5) = 0.3_dp * libcint_gto_norm(bas(LIBCINT_ANG_OF,n), env(bas(LIBCINT_PTR_EXP,n)+2)) + env(off + 6) = 0.2_dp * libcint_gto_norm(bas(LIBCINT_ANG_OF,n), env(bas(LIBCINT_PTR_EXP,n)+3)) + off = off + 6 + + ! Second basis function on atom 1 + n = 2 + bas(LIBCINT_ATOM_OF, n) = 0 + bas(LIBCINT_ANG_OF, n) = 1 + bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1 + bas(LIBCINT_PTR_EXP, n) = off + env(off + 1) = 0.9_dp + off = off + 1 + bas(LIBCINT_PTR_COEFF, n) = off + env(off + 1) = 1.0_dp * libcint_gto_norm(bas(LIBCINT_ANG_OF,n), env(bas(LIBCINT_PTR_EXP,n)+1)) + off = off + 1 + + ! Copy basis functions to atom 2 + n = 3 + bas(LIBCINT_ATOM_OF, n) = 1 + bas(LIBCINT_ANG_OF, n) = bas(LIBCINT_ANG_OF, 1) + bas(LIBCINT_NPRIM_OF, n) = bas(LIBCINT_NPRIM_OF, 1) + bas(LIBCINT_NCTR_OF, n) = bas(LIBCINT_NCTR_OF, 1) + bas(LIBCINT_PTR_EXP, n) = bas(LIBCINT_PTR_EXP, 1) + bas(LIBCINT_PTR_COEFF, n) = bas(LIBCINT_PTR_COEFF,1) + + n = 4 + bas(LIBCINT_ATOM_OF, n) = 1 + bas(LIBCINT_ANG_OF, n) = bas(LIBCINT_ANG_OF, 2) + bas(LIBCINT_NPRIM_OF, n) = bas(LIBCINT_NPRIM_OF, 2) + bas(LIBCINT_NCTR_OF, n) = bas(LIBCINT_NCTR_OF, 2) + bas(LIBCINT_PTR_EXP, n) = bas(LIBCINT_PTR_EXP, 2) + bas(LIBCINT_PTR_COEFF, n) = bas(LIBCINT_PTR_COEFF,2) + + print*, "============================================" + print*, "Pure Fortran Example (No C Types!)" + print*, "============================================" + print*, "" + + ! ======================================================================== + ! Compute overlap matrix between shells 0 and 2 (both s-orbitals on different atoms) + ! ======================================================================== + print*, "Computing overlap between shell 0 (atom 1) and shell 2 (atom 2)..." + + ! Shell 0 and 2 are both s-orbitals on different atoms + i = 0; shls(1) = i; di = libcint_cgto_sph(i, bas) + j = 2; shls(2) = j; dj = libcint_cgto_sph(j, bas) + + print*, "Shell dimensions: ", di, "x", dj + + ! Allocate buffer + allocate(buf_ovlp(di, dj)) + + ! Compute overlap integral + ret = libcint_1e_ovlp_sph(buf_ovlp, shls, atm, natm, bas, nbas, env) + if (ret /= 0) then + print*, "Overlap integral computed successfully!" + print*, "Overlap matrix (s-s between two H atoms):" + do i = 1, di + print '(10F12.6)', (buf_ovlp(i, j), j=1,dj) + end do + else + print*, "Overlap integral is zero (screened out)" + end if + deallocate(buf_ovlp) + + ! ======================================================================== + ! Compute kinetic energy for shell 0 with itself + ! ======================================================================== + print*, "" + print*, "Computing kinetic energy matrix for shell 0..." + + i = 0; shls(1) = i; di = libcint_cgto_sph(i, bas) + j = 0; shls(2) = j; dj = libcint_cgto_sph(j, bas) + + allocate(buf_kin(di, dj)) + ret = libcint_1e_kin_sph(buf_kin, shls, atm, natm, bas, nbas, env) + if (ret /= 0) then + print*, "Kinetic energy integral computed successfully!" + print*, "Kinetic energy matrix:" + do i = 1, di + print '(10F12.6)', (buf_kin(i, j), j=1,dj) + end do + end if + deallocate(buf_kin) + + ! ======================================================================== + ! Compute nuclear attraction for shell 0 + ! ======================================================================== + print*, "" + print*, "Computing nuclear attraction matrix for shell 0..." + + i = 0; shls(1) = i; di = libcint_cgto_sph(i, bas) + j = 0; shls(2) = j; dj = libcint_cgto_sph(j, bas) + + allocate(buf_nuc(di, dj)) + ret = libcint_1e_nuc_sph(buf_nuc, shls, atm, natm, bas, nbas, env) + if (ret /= 0) then + print*, "Nuclear attraction integral computed successfully!" + print*, "Nuclear attraction matrix:" + do i = 1, di + print '(10F12.6)', (buf_nuc(i, j), j=1,dj) + end do + end if + deallocate(buf_nuc) + + ! Clean up + deallocate(atm, bas, env) + + print*, "" + print*, "Example completed successfully!" + print*, "Notice: No iso_c_binding or C types in this program!" + +end program pure_fortran_example diff --git a/examples/fortran/fortran_time_c2h6.F90 b/examples/fortran/fortran_time_c2h6.F90 new file mode 100644 index 00000000..58f293d5 --- /dev/null +++ b/examples/fortran/fortran_time_c2h6.F90 @@ -0,0 +1,863 @@ +! +! Ethane (C2H6) molecule benchmark - Modern Fortran port +! +! This program benchmarks libcint performance with multiple basis sets +! using modern Fortran with iso_c_binding and OpenMP +! + +program fortran_time_c2h6 + use iso_c_binding + use libcint_interface + use omp_lib + implicit none + + integer(c_int), parameter :: natm = 8_c_int + integer(c_int), parameter :: nbas_max = natm * 20_c_int + integer(c_int), allocatable :: atm(:,:) + integer(c_int), allocatable :: bas(:,:) + real(c_double), allocatable :: env(:) + integer(c_int) :: nbas + integer(c_int) :: i, j, n, ia, off + + ! Allocate arrays + allocate(atm(ATM_SLOTS, natm)) + allocate(bas(BAS_SLOTS, nbas_max)) + allocate(env(10000)) + + atm = 0 + bas = 0 + env = 0.0_c_double + + off = PTR_ENV_START + + ! Set up ethane geometry (in Bohr) + atm(CHARGE_OF, 1) = 6; atm(PTR_COORD, 1) = off + env(off+1) = 0.000_c_double; env(off+2) = 0.000_c_double; env(off+3) = 0.769_c_double; off = off + 3 + + atm(CHARGE_OF, 2) = 1; atm(PTR_COORD, 2) = off + env(off+1) = 0.000_c_double; env(off+2) = 1.014_c_double; env(off+3) = 1.174_c_double; off = off + 3 + + atm(CHARGE_OF, 3) = 1; atm(PTR_COORD, 3) = off + env(off+1) = -0.878_c_double; env(off+2) = -0.507_c_double; env(off+3) = 1.174_c_double; off = off + 3 + + atm(CHARGE_OF, 4) = 1; atm(PTR_COORD, 4) = off + env(off+1) = 0.878_c_double; env(off+2) = -0.507_c_double; env(off+3) = 1.174_c_double; off = off + 3 + + atm(CHARGE_OF, 5) = 6; atm(PTR_COORD, 5) = off + env(off+1) = 0.000_c_double; env(off+2) = 0.000_c_double; env(off+3) = -0.769_c_double; off = off + 3 + + atm(CHARGE_OF, 6) = 1; atm(PTR_COORD, 6) = off + env(off+1) = 0.000_c_double; env(off+2) = 1.014_c_double; env(off+3) = -1.174_c_double; off = off + 3 + + atm(CHARGE_OF, 7) = 1; atm(PTR_COORD, 7) = off + env(off+1) = -0.878_c_double; env(off+2) = -0.507_c_double; env(off+3) = -1.174_c_double; off = off + 3 + + atm(CHARGE_OF, 8) = 1; atm(PTR_COORD, 8) = off + env(off+1) = 0.878_c_double; env(off+2) = -0.507_c_double; env(off+3) = -1.174_c_double; off = off + 3 + + ! ======================================================================== + ! 6-31G basis + ! ======================================================================== + call setup_6_31g_basis(bas, env, nbas, off) + print*, "6-31G basis" + call run_all(atm, natm, bas, nbas, env) + + ! ======================================================================== + ! 6-311G** basis + ! ======================================================================== + call setup_6_311gss_basis(bas, env, nbas, off) + print*, "6-311G(dp) basis" + call run_all(atm, natm, bas, nbas, env) + + ! ======================================================================== + ! cc-pVDZ basis + ! ======================================================================== + call setup_cc_pvdz_basis(bas, env, nbas, off) + print*, "cc-pVDZ basis" + call run_all(atm, natm, bas, nbas, env) + + ! ======================================================================== + ! cc-pVTZ basis + ! ======================================================================== + call setup_cc_pvtz_basis(bas, env, nbas, off) + print*, "cc-pVTZ basis" + call run_all(atm, natm, bas, nbas, env) + + ! ======================================================================== + ! cc-pVQZ basis + ! ======================================================================== + call setup_cc_pvqz_basis(bas, env, nbas, off) + print*, "cc-pVQZ basis" + call run_all(atm, natm, bas, nbas, env) + + deallocate(atm, bas, env) + +contains + + ! ======================================================================== + ! Run all benchmarks for a given basis set + ! ======================================================================== + subroutine run_all(atm, natm, bas, nbas, env) + integer(c_int), intent(in) :: atm(:,:) + integer(c_int), intent(in) :: natm + integer(c_int), intent(in) :: bas(:,:) + integer(c_int), intent(in) :: nbas + real(c_double), intent(in) :: env(:) + + integer(c_int) :: i, j, k, l, ij, kl + integer(c_int) :: di, dj, dk, dl + integer(c_int) :: kl_max + integer(c_int) :: shls(4) + real(c_double), allocatable :: buf(:) + integer(c_int), allocatable :: ishls(:), jshls(:) + integer(c_int) :: ncgto, npgto + integer(c_int) :: pct, count + real(8) :: time0, time1, tt, tot + type(c_ptr) :: opt_for_cint2e, opt_for_ip1 + integer(c_int) :: ret ! Return value from integral functions + + ! Set up shell pair lists + allocate(ishls(nbas*(nbas+1)/2)) + allocate(jshls(nbas*(nbas+1)/2)) + + ij = 0 + do i = 0, nbas-1 + do j = 0, i + ij = ij + 1 + ishls(ij) = i + jshls(ij) = j + end do + end do + + ncgto = CINTtot_cgto_spheric(bas, nbas) + npgto = CINTtot_pgto_spheric(bas, nbas) + print '(A,I0,A,I0,A,I0)', " shells = ", nbas, ", total cGTO = ", ncgto, & + ", total pGTO = ", npgto + + ! ==================================================================== + ! Create optimizers + ! ==================================================================== + opt_for_cint2e = c_null_ptr + call cint2e_sph_optimizer(opt_for_cint2e, atm, natm, bas, nbas, env) + opt_for_ip1 = c_null_ptr + call cint2e_ip1_sph_optimizer(opt_for_ip1, atm, natm, bas, nbas, env) + + ! ==================================================================== + ! Benchmark: cint2e_sph without optimizer + ! ==================================================================== + tot = real(ncgto, 8)**4 / 8.0_8 + print '(A,ES10.2)', " cint2e_sph without optimizer: total num ERI = ", tot + + pct = 0 + count = 0 + time0 = omp_get_wtime() + + !$omp parallel default(none) & + !$omp shared(atm, natm, bas, nbas, env, ishls, jshls, time0, pct, count) & + !$omp private(di, dj, dk, dl, i, j, k, l, ij, kl, kl_max, shls, buf, time1, ret) + !$omp do schedule(dynamic, 2) + do ij = 1, nbas*(nbas+1)/2 + i = ishls(ij) + j = jshls(ij) + di = CINTcgto_spheric(i, bas) + dj = CINTcgto_spheric(j, bas) + kl_max = (i+1)*(i+2)/2 + do kl = 1, kl_max + k = ishls(kl) + l = jshls(kl) + dk = CINTcgto_spheric(k, bas) + dl = CINTcgto_spheric(l, bas) + shls(1) = i + shls(2) = j + shls(3) = k + shls(4) = l + allocate(buf(di*dj*dk*dl)) + ret = cint2e_sph(buf, shls, atm, natm, bas, nbas, env, c_null_ptr) + deallocate(buf) + end do + !$omp atomic + count = count + kl_max + if (100*count / (int(nbas, 8)*nbas*(nbas+1)*(nbas+2)/8) > pct) then + !$omp atomic + pct = pct + 1 + time1 = omp_get_wtime() + !$omp critical + write(*, '(A,I0,A,F8.2)', advance='no') char(13)//" ", pct, "%, CPU time = ", time1-time0 + !$omp end critical + end if + end do + !$omp end do + !$omp end parallel + + time1 = omp_get_wtime() + tt = time1 - time0 + print '(A,I0,A,F8.2,A,F10.2,A)', char(13)//" ", 100, "%, CPU time = ", tt, & + ", ", tot/1e6_8/tt, " Mflops" + + ! ==================================================================== + ! Benchmark: cint2e_sph with optimizer + ! ==================================================================== + time0 = time1 + print '(A,ES10.2)', " cint2e_sph with optimizer: total num ERI = ", tot + + pct = 0 + count = 0 + + !$omp parallel default(none) & + !$omp shared(atm, natm, bas, nbas, env, ishls, jshls, opt_for_cint2e, time0, pct, count) & + !$omp private(di, dj, dk, dl, i, j, k, l, ij, kl, kl_max, shls, buf, time1, ret) + !$omp do schedule(dynamic, 2) + do ij = 1, nbas*(nbas+1)/2 + i = ishls(ij) + j = jshls(ij) + di = CINTcgto_spheric(i, bas) + dj = CINTcgto_spheric(j, bas) + kl_max = (i+1)*(i+2)/2 + do kl = 1, kl_max + k = ishls(kl) + l = jshls(kl) + dk = CINTcgto_spheric(k, bas) + dl = CINTcgto_spheric(l, bas) + shls(1) = i + shls(2) = j + shls(3) = k + shls(4) = l + allocate(buf(di*dj*dk*dl)) + ret = cint2e_sph(buf, shls, atm, natm, bas, nbas, env, opt_for_cint2e) + deallocate(buf) + end do + !$omp atomic + count = count + kl_max + if (100*count / (int(nbas, 8)*nbas*(nbas+1)*(nbas+2)/8) > pct) then + !$omp atomic + pct = pct + 1 + time1 = omp_get_wtime() + !$omp critical + write(*, '(A,I0,A,F8.2)', advance='no') char(13)//" ", pct, "%, CPU time = ", time1-time0 + !$omp end critical + end if + end do + !$omp end do + !$omp end parallel + + time1 = omp_get_wtime() + tt = time1 - time0 + print '(A,I0,A,F8.2,A,F10.2,A)', char(13)//" ", 100, "%, CPU time = ", tt, & + ", ", tot/1e6_8/tt, " Mflops" + + ! ==================================================================== + ! Benchmark: Gradients with optimizer + ! ==================================================================== + time0 = time1 + tot = real(ncgto, 8)**4 / 2.0_8 * 3.0_8 + print '(A,ES10.2)', " Gradients with optimizer: total num ERI = ", tot + + pct = 0 + count = 0 + + !$omp parallel default(none) & + !$omp shared(atm, natm, bas, nbas, env, ishls, jshls, opt_for_ip1, time0, pct, count) & + !$omp private(di, dj, dk, dl, i, j, k, l, ij, kl, shls, buf, time1, ret) + !$omp do schedule(dynamic, 2) + do ij = 0, nbas*nbas - 1 + i = ij / nbas + j = ij - nbas*i + di = CINTcgto_spheric(i, bas) + dj = CINTcgto_spheric(j, bas) + do kl = 1, nbas*(nbas+1)/2 + k = ishls(kl) + l = jshls(kl) + dk = CINTcgto_spheric(k, bas) + dl = CINTcgto_spheric(l, bas) + shls(1) = i + shls(2) = j + shls(3) = k + shls(4) = l + allocate(buf(di*dj*dk*dl*3)) + ret = cint2e_ip1_sph(buf, shls, atm, natm, bas, nbas, env, opt_for_ip1) + deallocate(buf) + end do + !$omp atomic + count = count + nbas*(nbas+1)/2 + if (100*count / (int(nbas, 8)*nbas*nbas*(nbas+1)/2) > pct) then + !$omp atomic + pct = pct + 1 + time1 = omp_get_wtime() + !$omp critical + write(*, '(A,I0,A,F8.2)', advance='no') char(13)//" ", pct, "%, CPU time = ", time1-time0 + !$omp end critical + end if + end do + !$omp end do + !$omp end parallel + + time1 = omp_get_wtime() + tt = time1 - time0 + print '(A,I0,A,F8.2,A,F10.2,A)', char(13)//" ", 100, "%, CPU time = ", tt, & + ", ", tot/1e6_8/tt, " Mflops" + print*, "" + + ! Clean up + call CINTdel_optimizer(opt_for_cint2e) + call CINTdel_optimizer(opt_for_ip1) + deallocate(ishls, jshls) + + end subroutine run_all + + ! ======================================================================== + ! 6-31G basis set setup + ! ======================================================================== + subroutine setup_6_31g_basis(bas, env, nbas, off) + integer(c_int), intent(inout) :: bas(:,:) + real(c_double), intent(inout) :: env(:) + integer(c_int), intent(out) :: nbas + integer(c_int), intent(inout) :: off + + integer(c_int) :: i, j, ia, n + + ! Carbon 6-31G exponents and coefficients + env(off+ 1) = 3047.5249_c_double; env(off+ 7) = 0.0018347_c_double * CINTgto_norm(0, env(off+1)) + env(off+ 2) = 457.36951_c_double; env(off+ 8) = 0.0140373_c_double * CINTgto_norm(0, env(off+2)) + env(off+ 3) = 103.94869_c_double; env(off+ 9) = 0.0688426_c_double * CINTgto_norm(0, env(off+3)) + env(off+ 4) = 29.210155_c_double; env(off+10) = 0.2321844_c_double * CINTgto_norm(0, env(off+4)) + env(off+ 5) = 9.2866630_c_double; env(off+11) = 0.4679413_c_double * CINTgto_norm(0, env(off+5)) + env(off+ 6) = 3.1639270_c_double; env(off+12) = 0.3623120_c_double * CINTgto_norm(0, env(off+6)) + env(off+13) = 7.8682724_c_double; env(off+16) =-0.1193324_c_double * CINTgto_norm(0, env(off+13)) + env(off+14) = 1.8812885_c_double; env(off+17) =-0.1608542_c_double * CINTgto_norm(0, env(off+14)) + env(off+15) = 0.5442493_c_double; env(off+18) = 1.1434564_c_double * CINTgto_norm(0, env(off+15)) + env(off+19) = 0.1687144_c_double; env(off+20) = 1.0000000_c_double * CINTgto_norm(0, env(off+19)) + env(off+21) = 7.8682724_c_double; env(off+24) = 0.0689991_c_double * CINTgto_norm(1, env(off+21)) + env(off+22) = 1.8812885_c_double; env(off+25) = 0.3164240_c_double * CINTgto_norm(1, env(off+22)) + env(off+23) = 0.5442493_c_double; env(off+26) = 0.7443083_c_double * CINTgto_norm(1, env(off+23)) + env(off+27) = 0.1687144_c_double; env(off+28) = 1.0000000_c_double * CINTgto_norm(1, env(off+27)) + + ! Hydrogen 6-31G + env(off+29) = 18.731137_c_double; env(off+32) = 0.0334946_c_double * CINTgto_norm(0, env(off+29)) + env(off+30) = 2.8253937_c_double; env(off+33) = 0.2347269_c_double * CINTgto_norm(0, env(off+30)) + env(off+31) = 0.6401217_c_double; env(off+34) = 0.8137573_c_double * CINTgto_norm(0, env(off+31)) + env(off+35) = 0.1612778_c_double; env(off+36) = 1.0000000_c_double * CINTgto_norm(0, env(off+35)) + + n = 1 + do i = 1, 2 ! Two carbons + ia = i - 1 + ! s orbital (6 primitives) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 6 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+1; bas(PTR_COEFF, n) = off+7 + n = n + 1 + ! s orbital (3 primitives) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 3 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+13; bas(PTR_COEFF, n) = off+16 + n = n + 1 + ! s orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+19; bas(PTR_COEFF, n) = off+20 + n = n + 1 + ! p orbital (3 primitives) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 1; bas(NPRIM_OF, n) = 3 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+21; bas(PTR_COEFF, n) = off+24 + n = n + 1 + ! p orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 1; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+27; bas(PTR_COEFF, n) = off+28 + n = n + 1 + ia = ia + 1 + + ! Three hydrogens attached to each carbon + do j = 1, 3 + ! s orbital (3 primitives) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 3 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+29; bas(PTR_COEFF, n) = off+32 + n = n + 1 + ! s orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+35; bas(PTR_COEFF, n) = off+36 + n = n + 1 + ia = ia + 1 + end do + end do + nbas = n - 1 + + end subroutine setup_6_31g_basis + + ! ======================================================================== + ! 6-311G** basis set setup + ! ======================================================================== + subroutine setup_6_311gss_basis(bas, env, nbas, off) + integer(c_int), intent(inout) :: bas(:,:) + real(c_double), intent(inout) :: env(:) + integer(c_int), intent(out) :: nbas + integer(c_int), intent(inout) :: off + + integer(c_int) :: i, j, ia, n + + ! Carbon 6-311G** + env(off+ 1) = 4563.240_c_double; env(off+18) = 0.0019666_c_double * CINTgto_norm(0, env(off+ 1)) + env(off+ 2) = 682.0240_c_double; env(off+19) = 0.0152306_c_double * CINTgto_norm(0, env(off+ 2)) + env(off+ 3) = 154.9730_c_double; env(off+20) = 0.0761269_c_double * CINTgto_norm(0, env(off+ 3)) + env(off+ 4) = 44.45530_c_double; env(off+21) = 0.2608010_c_double * CINTgto_norm(0, env(off+ 4)) + env(off+ 5) = 13.02900_c_double; env(off+22) = 0.6164620_c_double * CINTgto_norm(0, env(off+ 5)) + env(off+ 6) = 1.827730_c_double; env(off+23) = 0.2210060_c_double * CINTgto_norm(0, env(off+ 6)) + env(off+ 7) = 20.96420_c_double; env(off+24) = 0.1146600_c_double * CINTgto_norm(0, env(off+ 7)) + env(off+ 8) = 4.803310_c_double; env(off+25) = 0.9199990_c_double * CINTgto_norm(0, env(off+ 8)) + env(off+ 9) = 1.459330_c_double; env(off+26) = -0.003030_c_double * CINTgto_norm(0, env(off+ 9)) + env(off+10) = 0.483456_c_double; env(off+27) = 1.0000000_c_double * CINTgto_norm(0, env(off+10)) + env(off+11) = 0.145585_c_double; env(off+28) = 1.0000000_c_double * CINTgto_norm(0, env(off+11)) + env(off+12) = 20.96420_c_double; env(off+29) = 0.0402487_c_double * CINTgto_norm(1, env(off+12)) + env(off+13) = 4.803310_c_double; env(off+30) = 0.2375940_c_double * CINTgto_norm(1, env(off+13)) + env(off+14) = 1.459330_c_double; env(off+31) = 0.8158540_c_double * CINTgto_norm(1, env(off+14)) + env(off+15) = 0.483456_c_double; env(off+32) = 1.0000000_c_double * CINTgto_norm(1, env(off+15)) + env(off+16) = 0.145585_c_double; env(off+33) = 1.0000000_c_double * CINTgto_norm(1, env(off+16)) + env(off+17) = 0.626000_c_double; env(off+34) = 1.0000000_c_double * CINTgto_norm(2, env(off+17)) + + ! Hydrogen 6-311G** + env(off+35) = 33.86500_c_double; env(off+41) = 0.0254938_c_double * CINTgto_norm(0, env(off+35)) + env(off+36) = 5.094790_c_double; env(off+42) = 0.1903730_c_double * CINTgto_norm(0, env(off+36)) + env(off+37) = 1.158790_c_double; env(off+43) = 0.8521610_c_double * CINTgto_norm(0, env(off+37)) + env(off+38) = 0.325840_c_double; env(off+44) = 1.0000000_c_double * CINTgto_norm(0, env(off+38)) + env(off+39) = 0.102741_c_double; env(off+45) = 1.0000000_c_double * CINTgto_norm(0, env(off+39)) + env(off+40) = 0.750000_c_double; env(off+46) = 1.0000000_c_double * CINTgto_norm(1, env(off+40)) + + n = 1 + do i = 1, 2 ! Two carbons + ia = i - 1 + ! s orbital (6 primitives) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 6 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+ 1; bas(PTR_COEFF, n) = off+18 + n = n + 1 + ! s orbital (3 primitives) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 3 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+ 7; bas(PTR_COEFF, n) = off+24 + n = n + 1 + ! s orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+10; bas(PTR_COEFF, n) = off+27 + n = n + 1 + ! s orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+11; bas(PTR_COEFF, n) = off+28 + n = n + 1 + ! p orbital (3 primitives) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 1; bas(NPRIM_OF, n) = 3 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+12; bas(PTR_COEFF, n) = off+29 + n = n + 1 + ! p orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 1; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+15; bas(PTR_COEFF, n) = off+32 + n = n + 1 + ! p orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 1; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+16; bas(PTR_COEFF, n) = off+33 + n = n + 1 + ! d orbital (1 primitive) - polarization + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 2; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+17; bas(PTR_COEFF, n) = off+34 + n = n + 1 + ia = ia + 1 + + ! Three hydrogens attached to each carbon + do j = 1, 3 + ! s orbital (3 primitives) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 3 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+35; bas(PTR_COEFF, n) = off+41 + n = n + 1 + ! s orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+38; bas(PTR_COEFF, n) = off+44 + n = n + 1 + ! s orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+39; bas(PTR_COEFF, n) = off+45 + n = n + 1 + ! p orbital (1 primitive) - polarization + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 1; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+40; bas(PTR_COEFF, n) = off+46 + n = n + 1 + ia = ia + 1 + end do + end do + nbas = n - 1 + + end subroutine setup_6_311gss_basis + + ! ======================================================================== + ! cc-pVDZ basis set setup + ! ======================================================================== + subroutine setup_cc_pvdz_basis(bas, env, nbas, off) + integer(c_int), intent(inout) :: bas(:,:) + real(c_double), intent(inout) :: env(:) + integer(c_int), intent(out) :: nbas + integer(c_int), intent(inout) :: off + + integer(c_int) :: i, j, ia, n + + ! Carbon cc-pVDZ + env(off+ 1) = 6665.0_c_double; env(off+ 9) = 0.000692_c_double * CINTgto_norm(0, env(off+ 1)) + env(off+17) = -0.000146_c_double * CINTgto_norm(0, env(off+ 1)) + env(off+ 2) = 1000.0_c_double; env(off+10) = 0.005329_c_double * CINTgto_norm(0, env(off+ 2)) + env(off+18) = -0.001154_c_double * CINTgto_norm(0, env(off+ 2)) + env(off+ 3) = 228.00_c_double; env(off+11) = 0.027077_c_double * CINTgto_norm(0, env(off+ 3)) + env(off+19) = -0.005725_c_double * CINTgto_norm(0, env(off+ 3)) + env(off+ 4) = 64.710_c_double; env(off+12) = 0.101718_c_double * CINTgto_norm(0, env(off+ 4)) + env(off+20) = -0.023312_c_double * CINTgto_norm(0, env(off+ 4)) + env(off+ 5) = 21.060_c_double; env(off+13) = 0.274740_c_double * CINTgto_norm(0, env(off+ 5)) + env(off+21) = -0.063955_c_double * CINTgto_norm(0, env(off+ 5)) + env(off+ 6) = 7.4950_c_double; env(off+14) = 0.448564_c_double * CINTgto_norm(0, env(off+ 6)) + env(off+22) = -0.149981_c_double * CINTgto_norm(0, env(off+ 6)) + env(off+ 7) = 2.7970_c_double; env(off+15) = 0.285074_c_double * CINTgto_norm(0, env(off+ 7)) + env(off+23) = -0.127262_c_double * CINTgto_norm(0, env(off+ 7)) + env(off+ 8) = 0.5215_c_double; env(off+16) = 0.015204_c_double * CINTgto_norm(0, env(off+ 8)) + env(off+24) = 0.544529_c_double * CINTgto_norm(0, env(off+ 8)) + env(off+25) = 0.1596_c_double; env(off+26) = 1.000000_c_double * CINTgto_norm(0, env(off+25)) + env(off+27) = 9.4390_c_double; env(off+30) = 0.038109_c_double * CINTgto_norm(1, env(off+27)) + env(off+28) = 2.0020_c_double; env(off+31) = 0.209480_c_double * CINTgto_norm(1, env(off+28)) + env(off+29) = 0.5456_c_double; env(off+32) = 0.508557_c_double * CINTgto_norm(1, env(off+29)) + env(off+33) = 0.1517_c_double; env(off+34) = 1.000000_c_double * CINTgto_norm(1, env(off+33)) + env(off+35) = 0.5500_c_double; env(off+36) = 1.000000_c_double * CINTgto_norm(2, env(off+35)) + + ! Hydrogen cc-pVDZ + env(off+37) = 13.010_c_double; env(off+40) = 0.019685_c_double * CINTgto_norm(0, env(off+37)) + env(off+38) = 1.9620_c_double; env(off+41) = 0.137977_c_double * CINTgto_norm(0, env(off+38)) + env(off+39) = 0.4446_c_double; env(off+42) = 0.478148_c_double * CINTgto_norm(0, env(off+39)) + env(off+43) = 0.1220_c_double; env(off+44) = 1.000000_c_double * CINTgto_norm(0, env(off+43)) + env(off+45) = 0.7270_c_double; env(off+46) = 1.000000_c_double * CINTgto_norm(1, env(off+45)) + + n = 1 + do i = 1, 2 ! Two carbons + ia = i - 1 + ! s orbital (8 primitives, 2 contractions) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 8 + bas(NCTR_OF, n) = 2; bas(PTR_EXP, n) = off+ 1; bas(PTR_COEFF, n) = off+ 9 + n = n + 1 + ! s orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+25; bas(PTR_COEFF, n) = off+26 + n = n + 1 + ! p orbital (3 primitives) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 1; bas(NPRIM_OF, n) = 3 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+27; bas(PTR_COEFF, n) = off+30 + n = n + 1 + ! p orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 1; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+33; bas(PTR_COEFF, n) = off+34 + n = n + 1 + ! d orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 2; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+35; bas(PTR_COEFF, n) = off+36 + n = n + 1 + ia = ia + 1 + + ! Three hydrogens attached to each carbon + do j = 1, 3 + ! s orbital (3 primitives) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 3 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+37; bas(PTR_COEFF, n) = off+40 + n = n + 1 + ! s orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+43; bas(PTR_COEFF, n) = off+44 + n = n + 1 + ! p orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 1; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+45; bas(PTR_COEFF, n) = off+46 + n = n + 1 + ia = ia + 1 + end do + end do + nbas = n - 1 + + end subroutine setup_cc_pvdz_basis + + ! ======================================================================== + ! cc-pVTZ basis set setup + ! ======================================================================== + subroutine setup_cc_pvtz_basis(bas, env, nbas, off) + integer(c_int), intent(inout) :: bas(:,:) + real(c_double), intent(inout) :: env(:) + integer(c_int), intent(out) :: nbas + integer(c_int), intent(inout) :: off + + integer(c_int) :: i, j, ia, n + + ! Carbon cc-pVTZ + env(off+ 1) = 8236.0_c_double; env(off+19) = 0.000531_c_double * CINTgto_norm(0, env(off+ 1)) + env(off+27) = -0.000113_c_double * CINTgto_norm(0, env(off+ 1)) + env(off+ 2) = 1235.0_c_double; env(off+20) = 0.004108_c_double * CINTgto_norm(0, env(off+ 2)) + env(off+28) = -0.000878_c_double * CINTgto_norm(0, env(off+ 2)) + env(off+ 3) = 280.80_c_double; env(off+21) = 0.021087_c_double * CINTgto_norm(0, env(off+ 3)) + env(off+29) = -0.004540_c_double * CINTgto_norm(0, env(off+ 3)) + env(off+ 4) = 79.270_c_double; env(off+22) = 0.081853_c_double * CINTgto_norm(0, env(off+ 4)) + env(off+30) = -0.018133_c_double * CINTgto_norm(0, env(off+ 4)) + env(off+ 5) = 25.590_c_double; env(off+23) = 0.234817_c_double * CINTgto_norm(0, env(off+ 5)) + env(off+31) = -0.055760_c_double * CINTgto_norm(0, env(off+ 5)) + env(off+ 6) = 8.9970_c_double; env(off+24) = 0.434401_c_double * CINTgto_norm(0, env(off+ 6)) + env(off+32) = -0.126895_c_double * CINTgto_norm(0, env(off+ 6)) + env(off+ 7) = 3.3190_c_double; env(off+25) = 0.346129_c_double * CINTgto_norm(0, env(off+ 7)) + env(off+33) = -0.170352_c_double * CINTgto_norm(0, env(off+ 7)) + env(off+ 8) = 0.3643_c_double; env(off+26) = -0.008983_c_double * CINTgto_norm(0, env(off+ 8)) + env(off+34) = 0.598684_c_double * CINTgto_norm(0, env(off+ 8)) + env(off+ 9) = 0.9059_c_double; env(off+35) = 1.000000_c_double * CINTgto_norm(0, env(off+ 9)) + env(off+10) = 0.1285_c_double; env(off+36) = 1.000000_c_double * CINTgto_norm(0, env(off+10)) + env(off+11) = 18.710_c_double; env(off+37) = 0.014031_c_double * CINTgto_norm(1, env(off+11)) + env(off+12) = 4.1330_c_double; env(off+38) = 0.086866_c_double * CINTgto_norm(1, env(off+12)) + env(off+13) = 1.2000_c_double; env(off+39) = 0.290216_c_double * CINTgto_norm(1, env(off+13)) + env(off+14) = 0.3827_c_double; env(off+40) = 1.000000_c_double * CINTgto_norm(1, env(off+14)) + env(off+15) = 0.1209_c_double; env(off+41) = 1.000000_c_double * CINTgto_norm(1, env(off+15)) + env(off+16) = 1.0970_c_double; env(off+42) = 1.000000_c_double * CINTgto_norm(2, env(off+16)) + env(off+17) = 0.3180_c_double; env(off+43) = 1.000000_c_double * CINTgto_norm(2, env(off+17)) + env(off+18) = 0.7610_c_double; env(off+44) = 1.000000_c_double * CINTgto_norm(3, env(off+18)) + + ! Hydrogen cc-pVTZ + env(off+45) = 33.870_c_double; env(off+53) = 0.006068_c_double * CINTgto_norm(0, env(off+45)) + env(off+46) = 5.0950_c_double; env(off+54) = 0.045308_c_double * CINTgto_norm(0, env(off+46)) + env(off+47) = 1.1590_c_double; env(off+55) = 0.202822_c_double * CINTgto_norm(0, env(off+47)) + env(off+48) = 0.3258_c_double; env(off+56) = 1.000000_c_double * CINTgto_norm(0, env(off+48)) + env(off+49) = 0.1027_c_double; env(off+57) = 1.000000_c_double * CINTgto_norm(0, env(off+49)) + env(off+50) = 1.4070_c_double; env(off+58) = 1.000000_c_double * CINTgto_norm(1, env(off+50)) + env(off+51) = 0.3880_c_double; env(off+59) = 1.000000_c_double * CINTgto_norm(1, env(off+51)) + env(off+52) = 1.0570_c_double; env(off+60) = 1.000000_c_double * CINTgto_norm(2, env(off+52)) + + n = 1 + do i = 1, 2 ! Two carbons + ia = i - 1 + ! s orbital (8 primitives, 2 contractions) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 8 + bas(NCTR_OF, n) = 2; bas(PTR_EXP, n) = off+ 1; bas(PTR_COEFF, n) = off+19 + n = n + 1 + ! s orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+ 9; bas(PTR_COEFF, n) = off+35 + n = n + 1 + ! s orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+10; bas(PTR_COEFF, n) = off+36 + n = n + 1 + ! p orbital (3 primitives) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 1; bas(NPRIM_OF, n) = 3 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+11; bas(PTR_COEFF, n) = off+37 + n = n + 1 + ! p orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 1; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+14; bas(PTR_COEFF, n) = off+40 + n = n + 1 + ! p orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 1; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+15; bas(PTR_COEFF, n) = off+41 + n = n + 1 + ! d orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 2; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+16; bas(PTR_COEFF, n) = off+42 + n = n + 1 + ! d orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 2; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+17; bas(PTR_COEFF, n) = off+43 + n = n + 1 + ! f orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 3; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+18; bas(PTR_COEFF, n) = off+44 + n = n + 1 + ia = ia + 1 + + ! Three hydrogens attached to each carbon + do j = 1, 3 + ! s orbital (3 primitives) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 3 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+45; bas(PTR_COEFF, n) = off+53 + n = n + 1 + ! s orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+48; bas(PTR_COEFF, n) = off+56 + n = n + 1 + ! s orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+49; bas(PTR_COEFF, n) = off+57 + n = n + 1 + ! p orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 1; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+50; bas(PTR_COEFF, n) = off+58 + n = n + 1 + ! p orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 1; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+51; bas(PTR_COEFF, n) = off+59 + n = n + 1 + ! d orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 2; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+52; bas(PTR_COEFF, n) = off+60 + n = n + 1 + ia = ia + 1 + end do + end do + nbas = n - 1 + + end subroutine setup_cc_pvtz_basis + + ! ======================================================================== + ! cc-pVQZ basis set setup + ! ======================================================================== + subroutine setup_cc_pvqz_basis(bas, env, nbas, off) + integer(c_int), intent(inout) :: bas(:,:) + real(c_double), intent(inout) :: env(:) + integer(c_int), intent(out) :: nbas + integer(c_int), intent(inout) :: off + + integer(c_int) :: i, j, ia, n + + ! Carbon cc-pVQZ + env(off+ 1) = 33980.0_c_double; env(off+25) = 0.000091_c_double * CINTgto_norm(0, env(off+ 1)) + env(off+34) = -0.000019_c_double * CINTgto_norm(0, env(off+ 1)) + env(off+ 2) = 5089.00_c_double; env(off+26) = 0.000704_c_double * CINTgto_norm(0, env(off+ 2)) + env(off+35) = -0.000151_c_double * CINTgto_norm(0, env(off+ 2)) + env(off+ 3) = 1157.00_c_double; env(off+27) = 0.003693_c_double * CINTgto_norm(0, env(off+ 3)) + env(off+36) = -0.000785_c_double * CINTgto_norm(0, env(off+ 3)) + env(off+ 4) = 326.600_c_double; env(off+28) = 0.015360_c_double * CINTgto_norm(0, env(off+ 4)) + env(off+37) = -0.003324_c_double * CINTgto_norm(0, env(off+ 4)) + env(off+ 5) = 106.100_c_double; env(off+29) = 0.052929_c_double * CINTgto_norm(0, env(off+ 5)) + env(off+38) = -0.011512_c_double * CINTgto_norm(0, env(off+ 5)) + env(off+ 6) = 38.1100_c_double; env(off+30) = 0.147043_c_double * CINTgto_norm(0, env(off+ 6)) + env(off+39) = -0.034160_c_double * CINTgto_norm(0, env(off+ 6)) + env(off+ 7) = 14.7500_c_double; env(off+31) = 0.305631_c_double * CINTgto_norm(0, env(off+ 7)) + env(off+40) = -0.077173_c_double * CINTgto_norm(0, env(off+ 7)) + env(off+ 8) = 6.03500_c_double; env(off+32) = 0.399345_c_double * CINTgto_norm(0, env(off+ 8)) + env(off+41) = -0.141493_c_double * CINTgto_norm(0, env(off+ 8)) + env(off+ 9) = 2.53000_c_double; env(off+33) = 0.217051_c_double * CINTgto_norm(0, env(off+ 9)) + env(off+42) = -0.118019_c_double * CINTgto_norm(0, env(off+ 9)) + env(off+10) = 0.73550_c_double; env(off+43) = 1.000000_c_double * CINTgto_norm(0, env(off+10)) + env(off+11) = 0.29050_c_double; env(off+44) = 1.000000_c_double * CINTgto_norm(0, env(off+11)) + env(off+12) = 0.11110_c_double; env(off+45) = 1.000000_c_double * CINTgto_norm(0, env(off+12)) + env(off+13) = 34.5100_c_double; env(off+46) = 0.005378_c_double * CINTgto_norm(1, env(off+13)) + env(off+14) = 7.91500_c_double; env(off+47) = 0.036132_c_double * CINTgto_norm(1, env(off+14)) + env(off+15) = 2.36800_c_double; env(off+48) = 0.142493_c_double * CINTgto_norm(1, env(off+15)) + env(off+16) = 0.81320_c_double; env(off+49) = 1.000000_c_double * CINTgto_norm(1, env(off+16)) + env(off+17) = 0.28900_c_double; env(off+50) = 1.000000_c_double * CINTgto_norm(1, env(off+17)) + env(off+18) = 0.10070_c_double; env(off+51) = 1.000000_c_double * CINTgto_norm(1, env(off+18)) + env(off+19) = 1.84800_c_double; env(off+52) = 1.000000_c_double * CINTgto_norm(2, env(off+19)) + env(off+20) = 0.64900_c_double; env(off+53) = 1.000000_c_double * CINTgto_norm(2, env(off+20)) + env(off+21) = 0.22800_c_double; env(off+54) = 1.000000_c_double * CINTgto_norm(2, env(off+21)) + env(off+22) = 1.41900_c_double; env(off+55) = 1.000000_c_double * CINTgto_norm(3, env(off+22)) + env(off+23) = 0.48500_c_double; env(off+56) = 1.000000_c_double * CINTgto_norm(3, env(off+23)) + env(off+24) = 1.01100_c_double; env(off+57) = 1.000000_c_double * CINTgto_norm(4, env(off+24)) + + ! Hydrogen cc-pVQZ + env(off+58) = 82.640_c_double; env(off+70) = 0.002006_c_double * CINTgto_norm(0, env(off+58)) + env(off+59) = 12.410_c_double; env(off+71) = 0.015343_c_double * CINTgto_norm(0, env(off+59)) + env(off+60) = 2.8240_c_double; env(off+72) = 0.075579_c_double * CINTgto_norm(0, env(off+60)) + env(off+61) = 0.7970_c_double; env(off+73) = 1.000000_c_double * CINTgto_norm(0, env(off+61)) + env(off+62) = 0.2580_c_double; env(off+74) = 1.000000_c_double * CINTgto_norm(0, env(off+62)) + env(off+63) = 0.0890_c_double; env(off+75) = 1.000000_c_double * CINTgto_norm(0, env(off+63)) + env(off+64) = 2.2920_c_double; env(off+76) = 1.000000_c_double * CINTgto_norm(1, env(off+64)) + env(off+65) = 0.8380_c_double; env(off+77) = 1.000000_c_double * CINTgto_norm(1, env(off+65)) + env(off+66) = 0.2920_c_double; env(off+78) = 1.000000_c_double * CINTgto_norm(1, env(off+66)) + env(off+67) = 2.0620_c_double; env(off+79) = 1.000000_c_double * CINTgto_norm(2, env(off+67)) + env(off+68) = 0.6620_c_double; env(off+80) = 1.000000_c_double * CINTgto_norm(2, env(off+68)) + env(off+69) = 1.3970_c_double; env(off+81) = 1.000000_c_double * CINTgto_norm(3, env(off+69)) + + n = 1 + do i = 1, 2 ! Two carbons + ia = i - 1 + ! s orbital (9 primitives, 2 contractions) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 9 + bas(NCTR_OF, n) = 2; bas(PTR_EXP, n) = off+ 1; bas(PTR_COEFF, n) = off+25 + n = n + 1 + ! s orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+10; bas(PTR_COEFF, n) = off+43 + n = n + 1 + ! s orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+11; bas(PTR_COEFF, n) = off+44 + n = n + 1 + ! s orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+12; bas(PTR_COEFF, n) = off+45 + n = n + 1 + ! p orbital (3 primitives) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 1; bas(NPRIM_OF, n) = 3 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+13; bas(PTR_COEFF, n) = off+46 + n = n + 1 + ! p orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 1; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+16; bas(PTR_COEFF, n) = off+49 + n = n + 1 + ! p orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 1; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+17; bas(PTR_COEFF, n) = off+50 + n = n + 1 + ! p orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 1; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+18; bas(PTR_COEFF, n) = off+51 + n = n + 1 + ! d orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 2; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+19; bas(PTR_COEFF, n) = off+52 + n = n + 1 + ! d orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 2; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+20; bas(PTR_COEFF, n) = off+53 + n = n + 1 + ! d orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 2; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+21; bas(PTR_COEFF, n) = off+54 + n = n + 1 + ! f orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 3; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+22; bas(PTR_COEFF, n) = off+55 + n = n + 1 + ! f orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 3; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+23; bas(PTR_COEFF, n) = off+56 + n = n + 1 + ! g orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 4; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+24; bas(PTR_COEFF, n) = off+57 + n = n + 1 + ia = ia + 1 + + ! Three hydrogens attached to each carbon + do j = 1, 3 + ! s orbital (3 primitives) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 3 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+58; bas(PTR_COEFF, n) = off+70 + n = n + 1 + ! s orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+61; bas(PTR_COEFF, n) = off+73 + n = n + 1 + ! s orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+62; bas(PTR_COEFF, n) = off+74 + n = n + 1 + ! s orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 0; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+63; bas(PTR_COEFF, n) = off+75 + n = n + 1 + ! p orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 1; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+64; bas(PTR_COEFF, n) = off+76 + n = n + 1 + ! p orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 1; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+65; bas(PTR_COEFF, n) = off+77 + n = n + 1 + ! p orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 1; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+66; bas(PTR_COEFF, n) = off+78 + n = n + 1 + ! d orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 2; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+67; bas(PTR_COEFF, n) = off+79 + n = n + 1 + ! d orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 2; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+68; bas(PTR_COEFF, n) = off+80 + n = n + 1 + ! f orbital (1 primitive) + bas(ATOM_OF, n) = ia; bas(ANG_OF, n) = 3; bas(NPRIM_OF, n) = 1 + bas(NCTR_OF, n) = 1; bas(PTR_EXP, n) = off+69; bas(PTR_COEFF, n) = off+81 + n = n + 1 + ia = ia + 1 + end do + end do + nbas = n - 1 + + end subroutine setup_cc_pvqz_basis + +end program fortran_time_c2h6 diff --git a/examples/fortran/fortran_time_c2h6_pure.F90 b/examples/fortran/fortran_time_c2h6_pure.F90 new file mode 100644 index 00000000..605d81b3 --- /dev/null +++ b/examples/fortran/fortran_time_c2h6_pure.F90 @@ -0,0 +1,871 @@ +! +! Ethane (C2H6) molecule benchmark - Modern Fortran port +! +! This program benchmarks libcint performance with multiple basis sets +! using modern Fortran with iso_c_binding and OpenMP +! + +program fortran_time_c2h6_pure + use iso_c_binding, only: c_ptr, c_null_ptr + use libcint_fortran + use omp_lib + implicit none + + integer(ip), parameter :: natm = 8_ip + integer(ip), parameter :: nbas_max = natm * 20_ip + integer(ip), allocatable :: atm(:,:) + integer(ip), allocatable :: bas(:,:) + real(dp), allocatable :: env(:) + integer(ip) :: nbas + integer(ip) :: i, j, n, ia, off + + ! Allocate arrays + allocate(atm(LIBCINT_ATM_SLOTS, natm)) + allocate(bas(LIBCINT_BAS_SLOTS, nbas_max)) + allocate(env(10000)) + + atm = 0_ip + bas = 0_ip + env = 0.0_dp + + off = LIBCINT_PTR_ENV_START + + ! Set up ethane geometry (in Bohr) + atm(LIBCINT_CHARGE_OF, 1) = 6_ip; atm(LIBCINT_PTR_COORD, 1) = off + env(off+1) = 0.000_dp; env(off+2) = 0.000_dp; env(off+3) = 0.769_dp; off = off + 3 + + atm(LIBCINT_CHARGE_OF, 2) = 1_ip; atm(LIBCINT_PTR_COORD, 2) = off + env(off+1) = 0.000_dp; env(off+2) = 1.014_dp; env(off+3) = 1.174_dp; off = off + 3 + + atm(LIBCINT_CHARGE_OF, 3) = 1_ip; atm(LIBCINT_PTR_COORD, 3) = off + env(off+1) = -0.878_dp; env(off+2) = -0.507_dp; env(off+3) = 1.174_dp; off = off + 3 + + atm(LIBCINT_CHARGE_OF, 4) = 1_ip; atm(LIBCINT_PTR_COORD, 4) = off + env(off+1) = 0.878_dp; env(off+2) = -0.507_dp; env(off+3) = 1.174_dp; off = off + 3 + + atm(LIBCINT_CHARGE_OF, 5) = 6_ip; atm(LIBCINT_PTR_COORD, 5) = off + env(off+1) = 0.000_dp; env(off+2) = 0.000_dp; env(off+3) = -0.769_dp; off = off + 3 + + atm(LIBCINT_CHARGE_OF, 6) = 1_ip; atm(LIBCINT_PTR_COORD, 6) = off + env(off+1) = 0.000_dp; env(off+2) = 1.014_dp; env(off+3) = -1.174_dp; off = off + 3 + + atm(LIBCINT_CHARGE_OF, 7) = 1_ip; atm(LIBCINT_PTR_COORD, 7) = off + env(off+1) = -0.878_dp; env(off+2) = -0.507_dp; env(off+3) = -1.174_dp; off = off + 3 + + atm(LIBCINT_CHARGE_OF, 8) = 1_ip; atm(LIBCINT_PTR_COORD, 8) = off + env(off+1) = 0.878_dp; env(off+2) = -0.507_dp; env(off+3) = -1.174_dp; off = off + 3 + + ! ======================================================================== + ! 6-31G basis + ! ======================================================================== + call setup_6_31g_basis(bas, env, nbas, off) + print*, "6-31G basis" + call run_all(atm, natm, bas, nbas, env) + + ! ======================================================================== + ! 6-311G** basis + ! ======================================================================== + call setup_6_311gss_basis(bas, env, nbas, off) + print*, "6-311G(dp) basis" + call run_all(atm, natm, bas, nbas, env) + + ! ======================================================================== + ! cc-pVDZ basis + ! ======================================================================== + call setup_cc_pvdz_basis(bas, env, nbas, off) + print*, "cc-pVDZ basis" + call run_all(atm, natm, bas, nbas, env) + + ! ======================================================================== + ! cc-pVTZ basis + ! ======================================================================== + call setup_cc_pvtz_basis(bas, env, nbas, off) + print*, "cc-pVTZ basis" + call run_all(atm, natm, bas, nbas, env) + + ! ======================================================================== + ! cc-pVQZ basis + ! ======================================================================== + call setup_cc_pvqz_basis(bas, env, nbas, off) + print*, "cc-pVQZ basis" + call run_all(atm, natm, bas, nbas, env) + + deallocate(atm, bas, env) + +contains + + ! ======================================================================== + ! Run all benchmarks for a given basis set + ! ======================================================================== + subroutine run_all(atm, natm, bas, nbas, env) + integer(ip), intent(in) :: atm(:,:) + integer(ip), intent(in) :: natm + integer(ip), intent(in) :: bas(:,:) + integer(ip), intent(in) :: nbas + real(dp), intent(in) :: env(:) + + integer(ip) :: i, j, k, l, ij, kl + integer(ip) :: di, dj, dk, dl + integer(ip) :: kl_max + integer(ip) :: shls(4) + real(dp), allocatable :: buf(:) + integer(ip), allocatable :: ishls(:), jshls(:) + integer(ip) :: ncgto, npgto + integer(ip) :: pct, count + real(dp) :: time0, time1, tt, tot + real(dp) :: total_progress_eri, total_progress_grad + real(dp) :: nbas_real, progress + type(c_ptr) :: opt_for_cint2e, opt_for_ip1 + integer(ip) :: ret ! Return value from integral functions + + ! Set up shell pair lists + allocate(ishls(nbas*(nbas+1)/2)) + allocate(jshls(nbas*(nbas+1)/2)) + + ij = 0 + do i = 0, nbas-1 + do j = 0, i + ij = ij + 1 + ishls(ij) = i + jshls(ij) = j + end do + end do + + ncgto = libcint_tot_cgto_sph(bas, nbas) + npgto = libcint_tot_pgto_sph(bas, nbas) + nbas_real = real(nbas, dp) + print '(A,I0,A,I0,A,I0)', " shells = ", nbas, ", total cGTO = ", ncgto, & + ", total pGTO = ", npgto + + ! ==================================================================== + ! Create optimizers + ! ==================================================================== + opt_for_cint2e = c_null_ptr + call libcint_2e_sph_optimizer(opt_for_cint2e, atm, natm, bas, nbas, env) + opt_for_ip1 = c_null_ptr + call libcint_2e_ip1_sph_optimizer(opt_for_ip1, atm, natm, bas, nbas, env) + + ! ==================================================================== + ! Benchmark: libcint_2e_sph without optimizer + ! ==================================================================== + tot = real(ncgto, dp)**4 / 8.0_dp + total_progress_eri = max(1.0_dp, nbas_real*(nbas_real+1.0_dp)*(nbas_real+2.0_dp)/6.0_dp) + print '(A,ES10.2)', " libcint_2e_sph without optimizer: total num ERI = ", tot + + pct = 0 + count = 0 + time0 = omp_get_wtime() + + !$omp parallel default(none) & + !$omp shared(atm, natm, bas, nbas, env, ishls, jshls, time0, pct, count, total_progress_eri) & + !$omp private(di, dj, dk, dl, i, j, k, l, ij, kl, kl_max, shls, buf, time1, ret, progress) + !$omp do schedule(dynamic, 2) + do ij = 1, nbas*(nbas+1)/2 + i = ishls(ij) + j = jshls(ij) + di = libcint_cgto_sph(i, bas) + dj = libcint_cgto_sph(j, bas) + kl_max = (i+1)*(i+2)/2 + do kl = 1, kl_max + k = ishls(kl) + l = jshls(kl) + dk = libcint_cgto_sph(k, bas) + dl = libcint_cgto_sph(l, bas) + shls(1) = i + shls(2) = j + shls(3) = k + shls(4) = l + allocate(buf(di*dj*dk*dl)) + ret = libcint_2e_sph(buf, shls, atm, natm, bas, nbas, env, c_null_ptr) + deallocate(buf) + end do + !$omp atomic + count = count + kl_max + progress = min(100.0_dp, 100.0_dp*real(count, dp) / total_progress_eri) + if (progress > real(pct, dp)) then + !$omp atomic + pct = pct + 1 + time1 = omp_get_wtime() + !$omp critical + write(*, '(A,I0,A,F8.2)', advance='no') char(13)//" ", pct, "%, CPU time = ", time1-time0 + !$omp end critical + end if + end do + !$omp end do + !$omp end parallel + + time1 = omp_get_wtime() + tt = time1 - time0 + print '(A,I0,A,F8.2,A,F10.2,A)', char(13)//" ", 100, "%, CPU time = ", tt, & + ", ", tot/1e6_dp/tt, " Mflops" + + ! ==================================================================== + ! Benchmark: libcint_2e_sph with optimizer + ! ==================================================================== + time0 = time1 + print '(A,ES10.2)', " libcint_2e_sph with optimizer: total num ERI = ", tot + + pct = 0 + count = 0 + + !$omp parallel default(none) & + !$omp shared(atm, natm, bas, nbas, env, ishls, jshls, opt_for_cint2e, time0, pct, count, total_progress_eri) & + !$omp private(di, dj, dk, dl, i, j, k, l, ij, kl, kl_max, shls, buf, time1, ret, progress) + !$omp do schedule(dynamic, 2) + do ij = 1, nbas*(nbas+1)/2 + i = ishls(ij) + j = jshls(ij) + di = libcint_cgto_sph(i, bas) + dj = libcint_cgto_sph(j, bas) + kl_max = (i+1)*(i+2)/2 + do kl = 1, kl_max + k = ishls(kl) + l = jshls(kl) + dk = libcint_cgto_sph(k, bas) + dl = libcint_cgto_sph(l, bas) + shls(1) = i + shls(2) = j + shls(3) = k + shls(4) = l + allocate(buf(di*dj*dk*dl)) + ret = libcint_2e_sph(buf, shls, atm, natm, bas, nbas, env, opt_for_cint2e) + deallocate(buf) + end do + !$omp atomic + count = count + kl_max + progress = min(100.0_dp, 100.0_dp*real(count, dp) / total_progress_eri) + if (progress > real(pct, dp)) then + !$omp atomic + pct = pct + 1 + time1 = omp_get_wtime() + !$omp critical + write(*, '(A,I0,A,F8.2)', advance='no') char(13)//" ", pct, "%, CPU time = ", time1-time0 + !$omp end critical + end if + end do + !$omp end do + !$omp end parallel + + time1 = omp_get_wtime() + tt = time1 - time0 + print '(A,I0,A,F8.2,A,F10.2,A)', char(13)//" ", 100, "%, CPU time = ", tt, & + ", ", tot/1e6_dp/tt, " Mflops" + + ! ==================================================================== + ! Benchmark: Gradients with optimizer + ! ==================================================================== + time0 = time1 + tot = real(ncgto, dp)**4 / 2.0_dp * 3.0_dp + total_progress_grad = max(1.0_dp, (nbas_real**3)*(nbas_real+1.0_dp)/2.0_dp) + print '(A,ES10.2)', " Gradients with optimizer: total num ERI = ", tot + + pct = 0 + count = 0 + + !$omp parallel default(none) & + !$omp shared(atm, natm, bas, nbas, env, ishls, jshls, opt_for_ip1, time0, pct, count, total_progress_grad) & + !$omp private(di, dj, dk, dl, i, j, k, l, ij, kl, shls, buf, time1, ret, progress) + !$omp do schedule(dynamic, 2) + do ij = 0, nbas*nbas - 1 + i = ij / nbas + j = ij - nbas*i + di = libcint_cgto_sph(i, bas) + dj = libcint_cgto_sph(j, bas) + do kl = 1, nbas*(nbas+1)/2 + k = ishls(kl) + l = jshls(kl) + dk = libcint_cgto_sph(k, bas) + dl = libcint_cgto_sph(l, bas) + shls(1) = i + shls(2) = j + shls(3) = k + shls(4) = l + allocate(buf(di*dj*dk*dl*3)) + ret = libcint_2e_ip1_sph(buf, shls, atm, natm, bas, nbas, env, opt_for_ip1) + deallocate(buf) + end do + !$omp atomic + count = count + nbas*(nbas+1)/2 + progress = min(100.0_dp, 100.0_dp*real(count, dp) / total_progress_grad) + if (progress > real(pct, dp)) then + !$omp atomic + pct = pct + 1 + time1 = omp_get_wtime() + !$omp critical + write(*, '(A,I0,A,F8.2)', advance='no') char(13)//" ", pct, "%, CPU time = ", time1-time0 + !$omp end critical + end if + end do + !$omp end do + !$omp end parallel + + time1 = omp_get_wtime() + tt = time1 - time0 + print '(A,I0,A,F8.2,A,F10.2,A)', char(13)//" ", 100, "%, CPU time = ", tt, & + ", ", tot/1e6_dp/tt, " Mflops" + print*, "" + + ! Clean up + call libcint_del_optimizer(opt_for_cint2e) + call libcint_del_optimizer(opt_for_ip1) + deallocate(ishls, jshls) + + end subroutine run_all + + ! ======================================================================== + ! 6-31G basis set setup + ! ======================================================================== + subroutine setup_6_31g_basis(bas, env, nbas, off) + integer(ip), intent(inout) :: bas(:,:) + real(dp), intent(inout) :: env(:) + integer(ip), intent(out) :: nbas + integer(ip), intent(inout) :: off + + integer(ip) :: i, j, ia, n + + ! Carbon 6-31G exponents and coefficients + env(off+ 1) = 3047.5249_dp; env(off+ 7) = 0.0018347_dp * libcint_gto_norm(0, env(off+1)) + env(off+ 2) = 457.36951_dp; env(off+ 8) = 0.0140373_dp * libcint_gto_norm(0, env(off+2)) + env(off+ 3) = 103.94869_dp; env(off+ 9) = 0.0688426_dp * libcint_gto_norm(0, env(off+3)) + env(off+ 4) = 29.210155_dp; env(off+10) = 0.2321844_dp * libcint_gto_norm(0, env(off+4)) + env(off+ 5) = 9.2866630_dp; env(off+11) = 0.4679413_dp * libcint_gto_norm(0, env(off+5)) + env(off+ 6) = 3.1639270_dp; env(off+12) = 0.3623120_dp * libcint_gto_norm(0, env(off+6)) + env(off+13) = 7.8682724_dp; env(off+16) =-0.1193324_dp * libcint_gto_norm(0, env(off+13)) + env(off+14) = 1.8812885_dp; env(off+17) =-0.1608542_dp * libcint_gto_norm(0, env(off+14)) + env(off+15) = 0.5442493_dp; env(off+18) = 1.1434564_dp * libcint_gto_norm(0, env(off+15)) + env(off+19) = 0.1687144_dp; env(off+20) = 1.0000000_dp * libcint_gto_norm(0, env(off+19)) + env(off+21) = 7.8682724_dp; env(off+24) = 0.0689991_dp * libcint_gto_norm(1, env(off+21)) + env(off+22) = 1.8812885_dp; env(off+25) = 0.3164240_dp * libcint_gto_norm(1, env(off+22)) + env(off+23) = 0.5442493_dp; env(off+26) = 0.7443083_dp * libcint_gto_norm(1, env(off+23)) + env(off+27) = 0.1687144_dp; env(off+28) = 1.0000000_dp * libcint_gto_norm(1, env(off+27)) + + ! Hydrogen 6-31G + env(off+29) = 18.731137_dp; env(off+32) = 0.0334946_dp * libcint_gto_norm(0, env(off+29)) + env(off+30) = 2.8253937_dp; env(off+33) = 0.2347269_dp * libcint_gto_norm(0, env(off+30)) + env(off+31) = 0.6401217_dp; env(off+34) = 0.8137573_dp * libcint_gto_norm(0, env(off+31)) + env(off+35) = 0.1612778_dp; env(off+36) = 1.0000000_dp * libcint_gto_norm(0, env(off+35)) + + n = 1 + do i = 1, 2 ! Two carbons + ia = i - 1 + ! s orbital (6 primitives) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 6 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+1; bas(LIBCINT_PTR_COEFF, n) = off+7 + n = n + 1 + ! s orbital (3 primitives) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 3 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+13; bas(LIBCINT_PTR_COEFF, n) = off+16 + n = n + 1 + ! s orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+19; bas(LIBCINT_PTR_COEFF, n) = off+20 + n = n + 1 + ! p orbital (3 primitives) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 1; bas(LIBCINT_NPRIM_OF, n) = 3 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+21; bas(LIBCINT_PTR_COEFF, n) = off+24 + n = n + 1 + ! p orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 1; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+27; bas(LIBCINT_PTR_COEFF, n) = off+28 + n = n + 1 + ia = ia + 1 + + ! Three hydrogens attached to each carbon + do j = 1, 3 + ! s orbital (3 primitives) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 3 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+29; bas(LIBCINT_PTR_COEFF, n) = off+32 + n = n + 1 + ! s orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+35; bas(LIBCINT_PTR_COEFF, n) = off+36 + n = n + 1 + ia = ia + 1 + end do + end do + nbas = n - 1 + + end subroutine setup_6_31g_basis + + ! ======================================================================== + ! 6-311G** basis set setup + ! ======================================================================== + subroutine setup_6_311gss_basis(bas, env, nbas, off) + integer(ip), intent(inout) :: bas(:,:) + real(dp), intent(inout) :: env(:) + integer(ip), intent(out) :: nbas + integer(ip), intent(inout) :: off + + integer(ip) :: i, j, ia, n + + ! Carbon 6-311G** + env(off+ 1) = 4563.240_dp; env(off+18) = 0.0019666_dp * libcint_gto_norm(0, env(off+ 1)) + env(off+ 2) = 682.0240_dp; env(off+19) = 0.0152306_dp * libcint_gto_norm(0, env(off+ 2)) + env(off+ 3) = 154.9730_dp; env(off+20) = 0.0761269_dp * libcint_gto_norm(0, env(off+ 3)) + env(off+ 4) = 44.45530_dp; env(off+21) = 0.2608010_dp * libcint_gto_norm(0, env(off+ 4)) + env(off+ 5) = 13.02900_dp; env(off+22) = 0.6164620_dp * libcint_gto_norm(0, env(off+ 5)) + env(off+ 6) = 1.827730_dp; env(off+23) = 0.2210060_dp * libcint_gto_norm(0, env(off+ 6)) + env(off+ 7) = 20.96420_dp; env(off+24) = 0.1146600_dp * libcint_gto_norm(0, env(off+ 7)) + env(off+ 8) = 4.803310_dp; env(off+25) = 0.9199990_dp * libcint_gto_norm(0, env(off+ 8)) + env(off+ 9) = 1.459330_dp; env(off+26) = -0.003030_dp * libcint_gto_norm(0, env(off+ 9)) + env(off+10) = 0.483456_dp; env(off+27) = 1.0000000_dp * libcint_gto_norm(0, env(off+10)) + env(off+11) = 0.145585_dp; env(off+28) = 1.0000000_dp * libcint_gto_norm(0, env(off+11)) + env(off+12) = 20.96420_dp; env(off+29) = 0.0402487_dp * libcint_gto_norm(1, env(off+12)) + env(off+13) = 4.803310_dp; env(off+30) = 0.2375940_dp * libcint_gto_norm(1, env(off+13)) + env(off+14) = 1.459330_dp; env(off+31) = 0.8158540_dp * libcint_gto_norm(1, env(off+14)) + env(off+15) = 0.483456_dp; env(off+32) = 1.0000000_dp * libcint_gto_norm(1, env(off+15)) + env(off+16) = 0.145585_dp; env(off+33) = 1.0000000_dp * libcint_gto_norm(1, env(off+16)) + env(off+17) = 0.626000_dp; env(off+34) = 1.0000000_dp * libcint_gto_norm(2, env(off+17)) + + ! Hydrogen 6-311G** + env(off+35) = 33.86500_dp; env(off+41) = 0.0254938_dp * libcint_gto_norm(0, env(off+35)) + env(off+36) = 5.094790_dp; env(off+42) = 0.1903730_dp * libcint_gto_norm(0, env(off+36)) + env(off+37) = 1.158790_dp; env(off+43) = 0.8521610_dp * libcint_gto_norm(0, env(off+37)) + env(off+38) = 0.325840_dp; env(off+44) = 1.0000000_dp * libcint_gto_norm(0, env(off+38)) + env(off+39) = 0.102741_dp; env(off+45) = 1.0000000_dp * libcint_gto_norm(0, env(off+39)) + env(off+40) = 0.750000_dp; env(off+46) = 1.0000000_dp * libcint_gto_norm(1, env(off+40)) + + n = 1 + do i = 1, 2 ! Two carbons + ia = i - 1 + ! s orbital (6 primitives) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 6 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+ 1; bas(LIBCINT_PTR_COEFF, n) = off+18 + n = n + 1 + ! s orbital (3 primitives) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 3 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+ 7; bas(LIBCINT_PTR_COEFF, n) = off+24 + n = n + 1 + ! s orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+10; bas(LIBCINT_PTR_COEFF, n) = off+27 + n = n + 1 + ! s orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+11; bas(LIBCINT_PTR_COEFF, n) = off+28 + n = n + 1 + ! p orbital (3 primitives) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 1; bas(LIBCINT_NPRIM_OF, n) = 3 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+12; bas(LIBCINT_PTR_COEFF, n) = off+29 + n = n + 1 + ! p orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 1; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+15; bas(LIBCINT_PTR_COEFF, n) = off+32 + n = n + 1 + ! p orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 1; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+16; bas(LIBCINT_PTR_COEFF, n) = off+33 + n = n + 1 + ! d orbital (1 primitive) - polarization + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 2; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+17; bas(LIBCINT_PTR_COEFF, n) = off+34 + n = n + 1 + ia = ia + 1 + + ! Three hydrogens attached to each carbon + do j = 1, 3 + ! s orbital (3 primitives) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 3 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+35; bas(LIBCINT_PTR_COEFF, n) = off+41 + n = n + 1 + ! s orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+38; bas(LIBCINT_PTR_COEFF, n) = off+44 + n = n + 1 + ! s orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+39; bas(LIBCINT_PTR_COEFF, n) = off+45 + n = n + 1 + ! p orbital (1 primitive) - polarization + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 1; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+40; bas(LIBCINT_PTR_COEFF, n) = off+46 + n = n + 1 + ia = ia + 1 + end do + end do + nbas = n - 1 + + end subroutine setup_6_311gss_basis + + ! ======================================================================== + ! cc-pVDZ basis set setup + ! ======================================================================== + subroutine setup_cc_pvdz_basis(bas, env, nbas, off) + integer(ip), intent(inout) :: bas(:,:) + real(dp), intent(inout) :: env(:) + integer(ip), intent(out) :: nbas + integer(ip), intent(inout) :: off + + integer(ip) :: i, j, ia, n + + ! Carbon cc-pVDZ + env(off+ 1) = 6665.0_dp; env(off+ 9) = 0.000692_dp * libcint_gto_norm(0, env(off+ 1)) + env(off+17) = -0.000146_dp * libcint_gto_norm(0, env(off+ 1)) + env(off+ 2) = 1000.0_dp; env(off+10) = 0.005329_dp * libcint_gto_norm(0, env(off+ 2)) + env(off+18) = -0.001154_dp * libcint_gto_norm(0, env(off+ 2)) + env(off+ 3) = 228.00_dp; env(off+11) = 0.027077_dp * libcint_gto_norm(0, env(off+ 3)) + env(off+19) = -0.005725_dp * libcint_gto_norm(0, env(off+ 3)) + env(off+ 4) = 64.710_dp; env(off+12) = 0.101718_dp * libcint_gto_norm(0, env(off+ 4)) + env(off+20) = -0.023312_dp * libcint_gto_norm(0, env(off+ 4)) + env(off+ 5) = 21.060_dp; env(off+13) = 0.274740_dp * libcint_gto_norm(0, env(off+ 5)) + env(off+21) = -0.063955_dp * libcint_gto_norm(0, env(off+ 5)) + env(off+ 6) = 7.4950_dp; env(off+14) = 0.448564_dp * libcint_gto_norm(0, env(off+ 6)) + env(off+22) = -0.149981_dp * libcint_gto_norm(0, env(off+ 6)) + env(off+ 7) = 2.7970_dp; env(off+15) = 0.285074_dp * libcint_gto_norm(0, env(off+ 7)) + env(off+23) = -0.127262_dp * libcint_gto_norm(0, env(off+ 7)) + env(off+ 8) = 0.5215_dp; env(off+16) = 0.015204_dp * libcint_gto_norm(0, env(off+ 8)) + env(off+24) = 0.544529_dp * libcint_gto_norm(0, env(off+ 8)) + env(off+25) = 0.1596_dp; env(off+26) = 1.000000_dp * libcint_gto_norm(0, env(off+25)) + env(off+27) = 9.4390_dp; env(off+30) = 0.038109_dp * libcint_gto_norm(1, env(off+27)) + env(off+28) = 2.0020_dp; env(off+31) = 0.209480_dp * libcint_gto_norm(1, env(off+28)) + env(off+29) = 0.5456_dp; env(off+32) = 0.508557_dp * libcint_gto_norm(1, env(off+29)) + env(off+33) = 0.1517_dp; env(off+34) = 1.000000_dp * libcint_gto_norm(1, env(off+33)) + env(off+35) = 0.5500_dp; env(off+36) = 1.000000_dp * libcint_gto_norm(2, env(off+35)) + + ! Hydrogen cc-pVDZ + env(off+37) = 13.010_dp; env(off+40) = 0.019685_dp * libcint_gto_norm(0, env(off+37)) + env(off+38) = 1.9620_dp; env(off+41) = 0.137977_dp * libcint_gto_norm(0, env(off+38)) + env(off+39) = 0.4446_dp; env(off+42) = 0.478148_dp * libcint_gto_norm(0, env(off+39)) + env(off+43) = 0.1220_dp; env(off+44) = 1.000000_dp * libcint_gto_norm(0, env(off+43)) + env(off+45) = 0.7270_dp; env(off+46) = 1.000000_dp * libcint_gto_norm(1, env(off+45)) + + n = 1 + do i = 1, 2 ! Two carbons + ia = i - 1 + ! s orbital (8 primitives, 2 contractions) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 8 + bas(LIBCINT_NCTR_OF, n) = 2; bas(LIBCINT_PTR_EXP, n) = off+ 1; bas(LIBCINT_PTR_COEFF, n) = off+ 9 + n = n + 1 + ! s orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+25; bas(LIBCINT_PTR_COEFF, n) = off+26 + n = n + 1 + ! p orbital (3 primitives) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 1; bas(LIBCINT_NPRIM_OF, n) = 3 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+27; bas(LIBCINT_PTR_COEFF, n) = off+30 + n = n + 1 + ! p orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 1; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+33; bas(LIBCINT_PTR_COEFF, n) = off+34 + n = n + 1 + ! d orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 2; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+35; bas(LIBCINT_PTR_COEFF, n) = off+36 + n = n + 1 + ia = ia + 1 + + ! Three hydrogens attached to each carbon + do j = 1, 3 + ! s orbital (3 primitives) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 3 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+37; bas(LIBCINT_PTR_COEFF, n) = off+40 + n = n + 1 + ! s orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+43; bas(LIBCINT_PTR_COEFF, n) = off+44 + n = n + 1 + ! p orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 1; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+45; bas(LIBCINT_PTR_COEFF, n) = off+46 + n = n + 1 + ia = ia + 1 + end do + end do + nbas = n - 1 + + end subroutine setup_cc_pvdz_basis + + ! ======================================================================== + ! cc-pVTZ basis set setup + ! ======================================================================== + subroutine setup_cc_pvtz_basis(bas, env, nbas, off) + integer(ip), intent(inout) :: bas(:,:) + real(dp), intent(inout) :: env(:) + integer(ip), intent(out) :: nbas + integer(ip), intent(inout) :: off + + integer(ip) :: i, j, ia, n + + ! Carbon cc-pVTZ + env(off+ 1) = 8236.0_dp; env(off+19) = 0.000531_dp * libcint_gto_norm(0, env(off+ 1)) + env(off+27) = -0.000113_dp * libcint_gto_norm(0, env(off+ 1)) + env(off+ 2) = 1235.0_dp; env(off+20) = 0.004108_dp * libcint_gto_norm(0, env(off+ 2)) + env(off+28) = -0.000878_dp * libcint_gto_norm(0, env(off+ 2)) + env(off+ 3) = 280.80_dp; env(off+21) = 0.021087_dp * libcint_gto_norm(0, env(off+ 3)) + env(off+29) = -0.004540_dp * libcint_gto_norm(0, env(off+ 3)) + env(off+ 4) = 79.270_dp; env(off+22) = 0.081853_dp * libcint_gto_norm(0, env(off+ 4)) + env(off+30) = -0.018133_dp * libcint_gto_norm(0, env(off+ 4)) + env(off+ 5) = 25.590_dp; env(off+23) = 0.234817_dp * libcint_gto_norm(0, env(off+ 5)) + env(off+31) = -0.055760_dp * libcint_gto_norm(0, env(off+ 5)) + env(off+ 6) = 8.9970_dp; env(off+24) = 0.434401_dp * libcint_gto_norm(0, env(off+ 6)) + env(off+32) = -0.126895_dp * libcint_gto_norm(0, env(off+ 6)) + env(off+ 7) = 3.3190_dp; env(off+25) = 0.346129_dp * libcint_gto_norm(0, env(off+ 7)) + env(off+33) = -0.170352_dp * libcint_gto_norm(0, env(off+ 7)) + env(off+ 8) = 0.3643_dp; env(off+26) = -0.008983_dp * libcint_gto_norm(0, env(off+ 8)) + env(off+34) = 0.598684_dp * libcint_gto_norm(0, env(off+ 8)) + env(off+ 9) = 0.9059_dp; env(off+35) = 1.000000_dp * libcint_gto_norm(0, env(off+ 9)) + env(off+10) = 0.1285_dp; env(off+36) = 1.000000_dp * libcint_gto_norm(0, env(off+10)) + env(off+11) = 18.710_dp; env(off+37) = 0.014031_dp * libcint_gto_norm(1, env(off+11)) + env(off+12) = 4.1330_dp; env(off+38) = 0.086866_dp * libcint_gto_norm(1, env(off+12)) + env(off+13) = 1.2000_dp; env(off+39) = 0.290216_dp * libcint_gto_norm(1, env(off+13)) + env(off+14) = 0.3827_dp; env(off+40) = 1.000000_dp * libcint_gto_norm(1, env(off+14)) + env(off+15) = 0.1209_dp; env(off+41) = 1.000000_dp * libcint_gto_norm(1, env(off+15)) + env(off+16) = 1.0970_dp; env(off+42) = 1.000000_dp * libcint_gto_norm(2, env(off+16)) + env(off+17) = 0.3180_dp; env(off+43) = 1.000000_dp * libcint_gto_norm(2, env(off+17)) + env(off+18) = 0.7610_dp; env(off+44) = 1.000000_dp * libcint_gto_norm(3, env(off+18)) + + ! Hydrogen cc-pVTZ + env(off+45) = 33.870_dp; env(off+53) = 0.006068_dp * libcint_gto_norm(0, env(off+45)) + env(off+46) = 5.0950_dp; env(off+54) = 0.045308_dp * libcint_gto_norm(0, env(off+46)) + env(off+47) = 1.1590_dp; env(off+55) = 0.202822_dp * libcint_gto_norm(0, env(off+47)) + env(off+48) = 0.3258_dp; env(off+56) = 1.000000_dp * libcint_gto_norm(0, env(off+48)) + env(off+49) = 0.1027_dp; env(off+57) = 1.000000_dp * libcint_gto_norm(0, env(off+49)) + env(off+50) = 1.4070_dp; env(off+58) = 1.000000_dp * libcint_gto_norm(1, env(off+50)) + env(off+51) = 0.3880_dp; env(off+59) = 1.000000_dp * libcint_gto_norm(1, env(off+51)) + env(off+52) = 1.0570_dp; env(off+60) = 1.000000_dp * libcint_gto_norm(2, env(off+52)) + + n = 1 + do i = 1, 2 ! Two carbons + ia = i - 1 + ! s orbital (8 primitives, 2 contractions) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 8 + bas(LIBCINT_NCTR_OF, n) = 2; bas(LIBCINT_PTR_EXP, n) = off+ 1; bas(LIBCINT_PTR_COEFF, n) = off+19 + n = n + 1 + ! s orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+ 9; bas(LIBCINT_PTR_COEFF, n) = off+35 + n = n + 1 + ! s orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+10; bas(LIBCINT_PTR_COEFF, n) = off+36 + n = n + 1 + ! p orbital (3 primitives) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 1; bas(LIBCINT_NPRIM_OF, n) = 3 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+11; bas(LIBCINT_PTR_COEFF, n) = off+37 + n = n + 1 + ! p orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 1; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+14; bas(LIBCINT_PTR_COEFF, n) = off+40 + n = n + 1 + ! p orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 1; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+15; bas(LIBCINT_PTR_COEFF, n) = off+41 + n = n + 1 + ! d orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 2; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+16; bas(LIBCINT_PTR_COEFF, n) = off+42 + n = n + 1 + ! d orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 2; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+17; bas(LIBCINT_PTR_COEFF, n) = off+43 + n = n + 1 + ! f orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 3; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+18; bas(LIBCINT_PTR_COEFF, n) = off+44 + n = n + 1 + ia = ia + 1 + + ! Three hydrogens attached to each carbon + do j = 1, 3 + ! s orbital (3 primitives) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 3 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+45; bas(LIBCINT_PTR_COEFF, n) = off+53 + n = n + 1 + ! s orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+48; bas(LIBCINT_PTR_COEFF, n) = off+56 + n = n + 1 + ! s orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+49; bas(LIBCINT_PTR_COEFF, n) = off+57 + n = n + 1 + ! p orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 1; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+50; bas(LIBCINT_PTR_COEFF, n) = off+58 + n = n + 1 + ! p orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 1; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+51; bas(LIBCINT_PTR_COEFF, n) = off+59 + n = n + 1 + ! d orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 2; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+52; bas(LIBCINT_PTR_COEFF, n) = off+60 + n = n + 1 + ia = ia + 1 + end do + end do + nbas = n - 1 + + end subroutine setup_cc_pvtz_basis + + ! ======================================================================== + ! cc-pVQZ basis set setup + ! ======================================================================== + subroutine setup_cc_pvqz_basis(bas, env, nbas, off) + integer(ip), intent(inout) :: bas(:,:) + real(dp), intent(inout) :: env(:) + integer(ip), intent(out) :: nbas + integer(ip), intent(inout) :: off + + integer(ip) :: i, j, ia, n + + ! Carbon cc-pVQZ + env(off+ 1) = 33980.0_dp; env(off+25) = 0.000091_dp * libcint_gto_norm(0, env(off+ 1)) + env(off+34) = -0.000019_dp * libcint_gto_norm(0, env(off+ 1)) + env(off+ 2) = 5089.00_dp; env(off+26) = 0.000704_dp * libcint_gto_norm(0, env(off+ 2)) + env(off+35) = -0.000151_dp * libcint_gto_norm(0, env(off+ 2)) + env(off+ 3) = 1157.00_dp; env(off+27) = 0.003693_dp * libcint_gto_norm(0, env(off+ 3)) + env(off+36) = -0.000785_dp * libcint_gto_norm(0, env(off+ 3)) + env(off+ 4) = 326.600_dp; env(off+28) = 0.015360_dp * libcint_gto_norm(0, env(off+ 4)) + env(off+37) = -0.003324_dp * libcint_gto_norm(0, env(off+ 4)) + env(off+ 5) = 106.100_dp; env(off+29) = 0.052929_dp * libcint_gto_norm(0, env(off+ 5)) + env(off+38) = -0.011512_dp * libcint_gto_norm(0, env(off+ 5)) + env(off+ 6) = 38.1100_dp; env(off+30) = 0.147043_dp * libcint_gto_norm(0, env(off+ 6)) + env(off+39) = -0.034160_dp * libcint_gto_norm(0, env(off+ 6)) + env(off+ 7) = 14.7500_dp; env(off+31) = 0.305631_dp * libcint_gto_norm(0, env(off+ 7)) + env(off+40) = -0.077173_dp * libcint_gto_norm(0, env(off+ 7)) + env(off+ 8) = 6.03500_dp; env(off+32) = 0.399345_dp * libcint_gto_norm(0, env(off+ 8)) + env(off+41) = -0.141493_dp * libcint_gto_norm(0, env(off+ 8)) + env(off+ 9) = 2.53000_dp; env(off+33) = 0.217051_dp * libcint_gto_norm(0, env(off+ 9)) + env(off+42) = -0.118019_dp * libcint_gto_norm(0, env(off+ 9)) + env(off+10) = 0.73550_dp; env(off+43) = 1.000000_dp * libcint_gto_norm(0, env(off+10)) + env(off+11) = 0.29050_dp; env(off+44) = 1.000000_dp * libcint_gto_norm(0, env(off+11)) + env(off+12) = 0.11110_dp; env(off+45) = 1.000000_dp * libcint_gto_norm(0, env(off+12)) + env(off+13) = 34.5100_dp; env(off+46) = 0.005378_dp * libcint_gto_norm(1, env(off+13)) + env(off+14) = 7.91500_dp; env(off+47) = 0.036132_dp * libcint_gto_norm(1, env(off+14)) + env(off+15) = 2.36800_dp; env(off+48) = 0.142493_dp * libcint_gto_norm(1, env(off+15)) + env(off+16) = 0.81320_dp; env(off+49) = 1.000000_dp * libcint_gto_norm(1, env(off+16)) + env(off+17) = 0.28900_dp; env(off+50) = 1.000000_dp * libcint_gto_norm(1, env(off+17)) + env(off+18) = 0.10070_dp; env(off+51) = 1.000000_dp * libcint_gto_norm(1, env(off+18)) + env(off+19) = 1.84800_dp; env(off+52) = 1.000000_dp * libcint_gto_norm(2, env(off+19)) + env(off+20) = 0.64900_dp; env(off+53) = 1.000000_dp * libcint_gto_norm(2, env(off+20)) + env(off+21) = 0.22800_dp; env(off+54) = 1.000000_dp * libcint_gto_norm(2, env(off+21)) + env(off+22) = 1.41900_dp; env(off+55) = 1.000000_dp * libcint_gto_norm(3, env(off+22)) + env(off+23) = 0.48500_dp; env(off+56) = 1.000000_dp * libcint_gto_norm(3, env(off+23)) + env(off+24) = 1.01100_dp; env(off+57) = 1.000000_dp * libcint_gto_norm(4, env(off+24)) + + ! Hydrogen cc-pVQZ + env(off+58) = 82.640_dp; env(off+70) = 0.002006_dp * libcint_gto_norm(0, env(off+58)) + env(off+59) = 12.410_dp; env(off+71) = 0.015343_dp * libcint_gto_norm(0, env(off+59)) + env(off+60) = 2.8240_dp; env(off+72) = 0.075579_dp * libcint_gto_norm(0, env(off+60)) + env(off+61) = 0.7970_dp; env(off+73) = 1.000000_dp * libcint_gto_norm(0, env(off+61)) + env(off+62) = 0.2580_dp; env(off+74) = 1.000000_dp * libcint_gto_norm(0, env(off+62)) + env(off+63) = 0.0890_dp; env(off+75) = 1.000000_dp * libcint_gto_norm(0, env(off+63)) + env(off+64) = 2.2920_dp; env(off+76) = 1.000000_dp * libcint_gto_norm(1, env(off+64)) + env(off+65) = 0.8380_dp; env(off+77) = 1.000000_dp * libcint_gto_norm(1, env(off+65)) + env(off+66) = 0.2920_dp; env(off+78) = 1.000000_dp * libcint_gto_norm(1, env(off+66)) + env(off+67) = 2.0620_dp; env(off+79) = 1.000000_dp * libcint_gto_norm(2, env(off+67)) + env(off+68) = 0.6620_dp; env(off+80) = 1.000000_dp * libcint_gto_norm(2, env(off+68)) + env(off+69) = 1.3970_dp; env(off+81) = 1.000000_dp * libcint_gto_norm(3, env(off+69)) + + n = 1 + do i = 1, 2 ! Two carbons + ia = i - 1 + ! s orbital (9 primitives, 2 contractions) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 9 + bas(LIBCINT_NCTR_OF, n) = 2; bas(LIBCINT_PTR_EXP, n) = off+ 1; bas(LIBCINT_PTR_COEFF, n) = off+25 + n = n + 1 + ! s orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+10; bas(LIBCINT_PTR_COEFF, n) = off+43 + n = n + 1 + ! s orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+11; bas(LIBCINT_PTR_COEFF, n) = off+44 + n = n + 1 + ! s orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+12; bas(LIBCINT_PTR_COEFF, n) = off+45 + n = n + 1 + ! p orbital (3 primitives) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 1; bas(LIBCINT_NPRIM_OF, n) = 3 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+13; bas(LIBCINT_PTR_COEFF, n) = off+46 + n = n + 1 + ! p orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 1; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+16; bas(LIBCINT_PTR_COEFF, n) = off+49 + n = n + 1 + ! p orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 1; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+17; bas(LIBCINT_PTR_COEFF, n) = off+50 + n = n + 1 + ! p orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 1; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+18; bas(LIBCINT_PTR_COEFF, n) = off+51 + n = n + 1 + ! d orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 2; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+19; bas(LIBCINT_PTR_COEFF, n) = off+52 + n = n + 1 + ! d orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 2; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+20; bas(LIBCINT_PTR_COEFF, n) = off+53 + n = n + 1 + ! d orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 2; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+21; bas(LIBCINT_PTR_COEFF, n) = off+54 + n = n + 1 + ! f orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 3; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+22; bas(LIBCINT_PTR_COEFF, n) = off+55 + n = n + 1 + ! f orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 3; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+23; bas(LIBCINT_PTR_COEFF, n) = off+56 + n = n + 1 + ! g orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 4; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+24; bas(LIBCINT_PTR_COEFF, n) = off+57 + n = n + 1 + ia = ia + 1 + + ! Three hydrogens attached to each carbon + do j = 1, 3 + ! s orbital (3 primitives) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 3 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+58; bas(LIBCINT_PTR_COEFF, n) = off+70 + n = n + 1 + ! s orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+61; bas(LIBCINT_PTR_COEFF, n) = off+73 + n = n + 1 + ! s orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+62; bas(LIBCINT_PTR_COEFF, n) = off+74 + n = n + 1 + ! s orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 0; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+63; bas(LIBCINT_PTR_COEFF, n) = off+75 + n = n + 1 + ! p orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 1; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+64; bas(LIBCINT_PTR_COEFF, n) = off+76 + n = n + 1 + ! p orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 1; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+65; bas(LIBCINT_PTR_COEFF, n) = off+77 + n = n + 1 + ! p orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 1; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+66; bas(LIBCINT_PTR_COEFF, n) = off+78 + n = n + 1 + ! d orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 2; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+67; bas(LIBCINT_PTR_COEFF, n) = off+79 + n = n + 1 + ! d orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 2; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+68; bas(LIBCINT_PTR_COEFF, n) = off+80 + n = n + 1 + ! f orbital (1 primitive) + bas(LIBCINT_ATOM_OF, n) = ia; bas(LIBCINT_ANG_OF, n) = 3; bas(LIBCINT_NPRIM_OF, n) = 1 + bas(LIBCINT_NCTR_OF, n) = 1; bas(LIBCINT_PTR_EXP, n) = off+69; bas(LIBCINT_PTR_COEFF, n) = off+81 + n = n + 1 + ia = ia + 1 + end do + end do + nbas = n - 1 + + end subroutine setup_cc_pvqz_basis + +end program fortran_time_c2h6_pure diff --git a/include/libcint_fortran.f90 b/include/libcint_fortran.f90 new file mode 100644 index 00000000..bd2b9286 --- /dev/null +++ b/include/libcint_fortran.f90 @@ -0,0 +1,459 @@ +! +! High-level Fortran interface for libcint +! +! This module provides a pure Fortran API using standard Fortran types (dp, ip) +! while internally delegating to the C-compatible libcint_interface module. +! +module libcint_fortran + use iso_c_binding, only: c_int, c_double, c_double_complex, c_ptr, c_null_ptr + use iso_fortran_env, only: real64, int32 + use libcint_interface + implicit none + private + + ! ======================================================================== + ! Public API + ! ======================================================================== + + ! Type parameters for user code + public :: dp, ip, zp + + ! Constants (re-export from libcint_interface with LIBCINT_ prefix) + public :: LIBCINT_ATM_SLOTS, LIBCINT_BAS_SLOTS + public :: LIBCINT_CHARGE_OF, LIBCINT_PTR_COORD, LIBCINT_NUC_MOD_OF, LIBCINT_PTR_ZETA + public :: LIBCINT_ATOM_OF, LIBCINT_ANG_OF, LIBCINT_NPRIM_OF, LIBCINT_NCTR_OF + public :: LIBCINT_KAPPA_OF, LIBCINT_PTR_EXP, LIBCINT_PTR_COEFF + public :: LIBCINT_PTR_EXPCUTOFF, LIBCINT_PTR_COMMON_ORIG, LIBCINT_PTR_RINV_ORIG + public :: LIBCINT_PTR_RINV_ZETA, LIBCINT_PTR_RANGE_OMEGA, LIBCINT_PTR_ENV_START + + ! Dimension and normalization functions + public :: libcint_cgto_cart, libcint_cgto_sph, libcint_cgto_spinor + public :: libcint_tot_cgto_sph, libcint_tot_pgto_sph + public :: libcint_gto_norm + + ! One-electron integrals + public :: libcint_1e_ovlp_cart, libcint_1e_kin_cart, libcint_1e_nuc_cart + public :: libcint_1e_ipovlp_cart + public :: libcint_1e_ovlp_sph, libcint_1e_kin_sph, libcint_1e_nuc_sph + public :: libcint_1e_ipovlp_sph + public :: libcint_1e_spnucsp + + ! Two-electron integrals + public :: libcint_2e_cart, libcint_2e_sph + public :: libcint_2e_ip1_cart, libcint_2e_ip1_sph + public :: libcint_2e_spsp1 + + ! Optimizers + public :: libcint_2e_cart_optimizer, libcint_2e_sph_optimizer + public :: libcint_2e_ip1_cart_optimizer, libcint_2e_ip1_sph_optimizer + public :: libcint_del_optimizer + public :: libcint_2e_spsp1_optimizer + + ! ======================================================================== + ! Type parameters + ! ======================================================================== + + ! Bind directly to ISO C kinds to guarantee binary compatibility + integer, parameter :: dp = c_double ! Double precision + integer, parameter :: ip = c_int ! Integer + integer, parameter :: zp = c_double_complex ! Double complex + + ! Optional asserts: fail to compile if real64/int32 are not C-compatible + integer, parameter :: libcint_real64_is_c_double = 1 / merge(1, 0, real64 == c_double) + integer, parameter :: libcint_int32_is_c_int = 1 / merge(1, 0, int32 == c_int) + + ! ======================================================================== + ! Constants (re-export with LIBCINT_ prefix) + ! ======================================================================== + + integer(c_int), parameter :: LIBCINT_ATM_SLOTS = ATM_SLOTS + integer(c_int), parameter :: LIBCINT_BAS_SLOTS = BAS_SLOTS + integer(c_int), parameter :: LIBCINT_CHARGE_OF = CHARGE_OF + integer(c_int), parameter :: LIBCINT_PTR_COORD = PTR_COORD + integer(c_int), parameter :: LIBCINT_NUC_MOD_OF = NUC_MOD_OF + integer(c_int), parameter :: LIBCINT_PTR_ZETA = PTR_ZETA + integer(c_int), parameter :: LIBCINT_ATOM_OF = ATOM_OF + integer(c_int), parameter :: LIBCINT_ANG_OF = ANG_OF + integer(c_int), parameter :: LIBCINT_NPRIM_OF = NPRIM_OF + integer(c_int), parameter :: LIBCINT_NCTR_OF = NCTR_OF + integer(c_int), parameter :: LIBCINT_KAPPA_OF = KAPPA_OF + integer(c_int), parameter :: LIBCINT_PTR_EXP = PTR_EXP + integer(c_int), parameter :: LIBCINT_PTR_COEFF = PTR_COEFF + integer(c_int), parameter :: LIBCINT_PTR_EXPCUTOFF = PTR_EXPCUTOFF + integer(c_int), parameter :: LIBCINT_PTR_COMMON_ORIG = PTR_COMMON_ORIG + integer(c_int), parameter :: LIBCINT_PTR_RINV_ORIG = PTR_RINV_ORIG + integer(c_int), parameter :: LIBCINT_PTR_RINV_ZETA = PTR_RINV_ZETA + integer(c_int), parameter :: LIBCINT_PTR_RANGE_OMEGA = PTR_RANGE_OMEGA + integer(c_int), parameter :: LIBCINT_PTR_ENV_START = PTR_ENV_START + +contains + + ! ======================================================================== + ! Dimension inquiry functions + ! ======================================================================== + + !> Get number of Cartesian GTOs for a given shell + function libcint_cgto_cart(bas_id, bas) result(dim) + integer(ip), intent(in) :: bas_id + integer(ip), intent(in) :: bas(LIBCINT_BAS_SLOTS, *) + integer(ip) :: dim + dim = CINTcgto_cart(bas_id, bas) + end function libcint_cgto_cart + + !> Get number of spherical GTOs for a given shell + function libcint_cgto_sph(bas_id, bas) result(dim) + integer(ip), intent(in) :: bas_id + integer(ip), intent(in) :: bas(LIBCINT_BAS_SLOTS, *) + integer(ip) :: dim + dim = CINTcgto_spheric(bas_id, bas) + end function libcint_cgto_sph + + !> Get number of spinor components for a given shell + function libcint_cgto_spinor(bas_id, bas) result(dim) + integer(ip), intent(in) :: bas_id + integer(ip), intent(in) :: bas(LIBCINT_BAS_SLOTS, *) + integer(ip) :: dim + dim = CINTcgto_spinor(bas_id, bas) + end function libcint_cgto_spinor + + !> Total number of contracted spherical GTOs + function libcint_tot_cgto_sph(bas, nbas) result(ntot) + integer(ip), intent(in) :: bas(LIBCINT_BAS_SLOTS, *) + integer(ip), intent(in) :: nbas + integer(ip) :: ntot + ntot = CINTtot_cgto_spheric(bas, nbas) + end function libcint_tot_cgto_sph + + !> Total number of primitive spherical GTOs + function libcint_tot_pgto_sph(bas, nbas) result(ntot) + integer(ip), intent(in) :: bas(LIBCINT_BAS_SLOTS, *) + integer(ip), intent(in) :: nbas + integer(ip) :: ntot + ntot = CINTtot_pgto_spheric(bas, nbas) + end function libcint_tot_pgto_sph + + ! ======================================================================== + ! Normalization + ! ======================================================================== + + !> Compute GTO normalization factor for angular momentum n and exponent alpha + function libcint_gto_norm(n, alpha) result(norm) + integer(ip), intent(in) :: n + real(dp), intent(in) :: alpha + real(dp) :: norm + norm = CINTgto_norm(n, alpha) + end function libcint_gto_norm + + ! ======================================================================== + ! One-electron integrals - Cartesian + ! ======================================================================== + + !> Overlap integral (Cartesian basis) + function libcint_1e_ovlp_cart(buf, shls, atm, natm, bas, nbas, env) result(ret) + real(dp), intent(out) :: buf(*) + integer(ip), intent(in) :: shls(2) + integer(ip), intent(in) :: atm(LIBCINT_ATM_SLOTS, *) + integer(ip), intent(in) :: natm + integer(ip), intent(in) :: bas(LIBCINT_BAS_SLOTS, *) + integer(ip), intent(in) :: nbas + real(dp), intent(in) :: env(*) + integer(ip) :: ret + ret = cint1e_ovlp_cart(buf, shls, atm, natm, bas, nbas, env) + end function libcint_1e_ovlp_cart + + !> Kinetic energy integral (Cartesian basis) + function libcint_1e_kin_cart(buf, shls, atm, natm, bas, nbas, env) result(ret) + real(dp), intent(out) :: buf(*) + integer(ip), intent(in) :: shls(2) + integer(ip), intent(in) :: atm(LIBCINT_ATM_SLOTS, *) + integer(ip), intent(in) :: natm + integer(ip), intent(in) :: bas(LIBCINT_BAS_SLOTS, *) + integer(ip), intent(in) :: nbas + real(dp), intent(in) :: env(*) + integer(ip) :: ret + ret = cint1e_kin_cart(buf, shls, atm, natm, bas, nbas, env) + end function libcint_1e_kin_cart + + !> Nuclear attraction integral (Cartesian basis) + function libcint_1e_nuc_cart(buf, shls, atm, natm, bas, nbas, env) result(ret) + real(dp), intent(out) :: buf(*) + integer(ip), intent(in) :: shls(2) + integer(ip), intent(in) :: atm(LIBCINT_ATM_SLOTS, *) + integer(ip), intent(in) :: natm + integer(ip), intent(in) :: bas(LIBCINT_BAS_SLOTS, *) + integer(ip), intent(in) :: nbas + real(dp), intent(in) :: env(*) + integer(ip) :: ret + ret = cint1e_nuc_cart(buf, shls, atm, natm, bas, nbas, env) + end function libcint_1e_nuc_cart + + !> Overlap gradient integral (Cartesian basis) + function libcint_1e_ipovlp_cart(buf, shls, atm, natm, bas, nbas, env) result(ret) + real(dp), intent(out) :: buf(*) + integer(ip), intent(in) :: shls(2) + integer(ip), intent(in) :: atm(LIBCINT_ATM_SLOTS, *) + integer(ip), intent(in) :: natm + integer(ip), intent(in) :: bas(LIBCINT_BAS_SLOTS, *) + integer(ip), intent(in) :: nbas + real(dp), intent(in) :: env(*) + integer(ip) :: ret + ret = cint1e_ipovlp_cart(buf, shls, atm, natm, bas, nbas, env) + end function libcint_1e_ipovlp_cart + + ! ======================================================================== + ! One-electron integrals - Spherical + ! ======================================================================== + + !> Overlap integral (Spherical basis) + function libcint_1e_ovlp_sph(buf, shls, atm, natm, bas, nbas, env) result(ret) + real(dp), intent(out) :: buf(*) + integer(ip), intent(in) :: shls(2) + integer(ip), intent(in) :: atm(LIBCINT_ATM_SLOTS, *) + integer(ip), intent(in) :: natm + integer(ip), intent(in) :: bas(LIBCINT_BAS_SLOTS, *) + integer(ip), intent(in) :: nbas + real(dp), intent(in) :: env(*) + integer(ip) :: ret + ret = cint1e_ovlp_sph(buf, shls, atm, natm, bas, nbas, env) + end function libcint_1e_ovlp_sph + + !> Kinetic energy integral (Spherical basis) + function libcint_1e_kin_sph(buf, shls, atm, natm, bas, nbas, env) result(ret) + real(dp), intent(out) :: buf(*) + integer(ip), intent(in) :: shls(2) + integer(ip), intent(in) :: atm(LIBCINT_ATM_SLOTS, *) + integer(ip), intent(in) :: natm + integer(ip), intent(in) :: bas(LIBCINT_BAS_SLOTS, *) + integer(ip), intent(in) :: nbas + real(dp), intent(in) :: env(*) + integer(ip) :: ret + ret = cint1e_kin_sph(buf, shls, atm, natm, bas, nbas, env) + end function libcint_1e_kin_sph + + !> Nuclear attraction integral (Spherical basis) + function libcint_1e_nuc_sph(buf, shls, atm, natm, bas, nbas, env) result(ret) + real(dp), intent(out) :: buf(*) + integer(ip), intent(in) :: shls(2) + integer(ip), intent(in) :: atm(LIBCINT_ATM_SLOTS, *) + integer(ip), intent(in) :: natm + integer(ip), intent(in) :: bas(LIBCINT_BAS_SLOTS, *) + integer(ip), intent(in) :: nbas + real(dp), intent(in) :: env(*) + integer(ip) :: ret + ret = cint1e_nuc_sph(buf, shls, atm, natm, bas, nbas, env) + end function libcint_1e_nuc_sph + + !> Overlap gradient integral (Spherical basis) + function libcint_1e_ipovlp_sph(buf, shls, atm, natm, bas, nbas, env) result(ret) + real(dp), intent(out) :: buf(*) + integer(ip), intent(in) :: shls(2) + integer(ip), intent(in) :: atm(LIBCINT_ATM_SLOTS, *) + integer(ip), intent(in) :: natm + integer(ip), intent(in) :: bas(LIBCINT_BAS_SLOTS, *) + integer(ip), intent(in) :: nbas + real(dp), intent(in) :: env(*) + integer(ip) :: ret + ret = cint1e_ipovlp_sph(buf, shls, atm, natm, bas, nbas, env) + end function libcint_1e_ipovlp_sph + + ! ======================================================================== + ! One-electron integrals - Spinor + ! ======================================================================== + + !> Spinor nuclear attraction integral (complex spinor basis) + subroutine libcint_1e_spnucsp(buf, shls, atm, natm, bas, nbas, env) + complex(zp), intent(out) :: buf(*) + integer(ip), intent(in) :: shls(2) + integer(ip), intent(in) :: atm(LIBCINT_ATM_SLOTS, *) + integer(ip), intent(in) :: natm + integer(ip), intent(in) :: bas(LIBCINT_BAS_SLOTS, *) + integer(ip), intent(in) :: nbas + real(dp), intent(in) :: env(*) + + call cint1e_spnucsp(buf, shls, atm, natm, bas, nbas, env) + end subroutine libcint_1e_spnucsp + + ! ======================================================================== + ! Two-electron integrals - Cartesian + ! ======================================================================== + + !> Electron repulsion integral (Cartesian basis) + !> Optional optimizer argument for better performance + function libcint_2e_cart(buf, shls, atm, natm, bas, nbas, env, opt) result(ret) + real(dp), intent(out) :: buf(*) + integer(ip), intent(in) :: shls(4) + integer(ip), intent(in) :: atm(LIBCINT_ATM_SLOTS, *) + integer(ip), intent(in) :: natm + integer(ip), intent(in) :: bas(LIBCINT_BAS_SLOTS, *) + integer(ip), intent(in) :: nbas + real(dp), intent(in) :: env(*) + type(c_ptr), intent(in), optional :: opt + integer(ip) :: ret + + if (present(opt)) then + ret = cint2e_cart(buf, shls, atm, natm, bas, nbas, env, opt) + else + ret = cint2e_cart(buf, shls, atm, natm, bas, nbas, env, c_null_ptr) + end if + end function libcint_2e_cart + + !> ERI gradient (Cartesian basis) + !> Optional optimizer argument for better performance + function libcint_2e_ip1_cart(buf, shls, atm, natm, bas, nbas, env, opt) result(ret) + real(dp), intent(out) :: buf(*) + integer(ip), intent(in) :: shls(4) + integer(ip), intent(in) :: atm(LIBCINT_ATM_SLOTS, *) + integer(ip), intent(in) :: natm + integer(ip), intent(in) :: bas(LIBCINT_BAS_SLOTS, *) + integer(ip), intent(in) :: nbas + real(dp), intent(in) :: env(*) + type(c_ptr), intent(in), optional :: opt + integer(ip) :: ret + + if (present(opt)) then + ret = cint2e_ip1_cart(buf, shls, atm, natm, bas, nbas, env, opt) + else + ret = cint2e_ip1_cart(buf, shls, atm, natm, bas, nbas, env, c_null_ptr) + end if + end function libcint_2e_ip1_cart + + ! ======================================================================== + ! Two-electron integrals - Spherical + ! ======================================================================== + + !> Electron repulsion integral (Spherical basis) + !> Optional optimizer argument for better performance + function libcint_2e_sph(buf, shls, atm, natm, bas, nbas, env, opt) result(ret) + real(dp), intent(out) :: buf(*) + integer(ip), intent(in) :: shls(4) + integer(ip), intent(in) :: atm(LIBCINT_ATM_SLOTS, *) + integer(ip), intent(in) :: natm + integer(ip), intent(in) :: bas(LIBCINT_BAS_SLOTS, *) + integer(ip), intent(in) :: nbas + real(dp), intent(in) :: env(*) + type(c_ptr), intent(in), optional :: opt + integer(ip) :: ret + + if (present(opt)) then + ret = cint2e_sph(buf, shls, atm, natm, bas, nbas, env, opt) + else + ret = cint2e_sph(buf, shls, atm, natm, bas, nbas, env, c_null_ptr) + end if + end function libcint_2e_sph + + !> ERI gradient (Spherical basis) + !> Optional optimizer argument for better performance + function libcint_2e_ip1_sph(buf, shls, atm, natm, bas, nbas, env, opt) result(ret) + real(dp), intent(out) :: buf(*) + integer(ip), intent(in) :: shls(4) + integer(ip), intent(in) :: atm(LIBCINT_ATM_SLOTS, *) + integer(ip), intent(in) :: natm + integer(ip), intent(in) :: bas(LIBCINT_BAS_SLOTS, *) + integer(ip), intent(in) :: nbas + real(dp), intent(in) :: env(*) + type(c_ptr), intent(in), optional :: opt + integer(ip) :: ret + + if (present(opt)) then + ret = cint2e_ip1_sph(buf, shls, atm, natm, bas, nbas, env, opt) + else + ret = cint2e_ip1_sph(buf, shls, atm, natm, bas, nbas, env, c_null_ptr) + end if + end function libcint_2e_ip1_sph + + ! ======================================================================== + ! Two-electron integrals - Spinor + ! ======================================================================== + + !> Spinor electron repulsion integral + !> Optional optimizer argument for better performance + subroutine libcint_2e_spsp1(buf, shls, atm, natm, bas, nbas, env, opt) + complex(zp), intent(out) :: buf(*) + integer(ip), intent(in) :: shls(4) + integer(ip), intent(in) :: atm(LIBCINT_ATM_SLOTS, *) + integer(ip), intent(in) :: natm + integer(ip), intent(in) :: bas(LIBCINT_BAS_SLOTS, *) + integer(ip), intent(in) :: nbas + real(dp), intent(in) :: env(*) + type(c_ptr), intent(in), optional :: opt + + if (present(opt)) then + call cint2e_spsp1(buf, shls, atm, natm, bas, nbas, env, opt) + else + call cint2e_spsp1(buf, shls, atm, natm, bas, nbas, env, c_null_ptr) + end if + end subroutine libcint_2e_spsp1 + + ! ======================================================================== + ! Optimizer creation and deletion + ! ======================================================================== + + !> Create optimizer for 2e integrals (Cartesian basis) + subroutine libcint_2e_cart_optimizer(opt, atm, natm, bas, nbas, env) + type(c_ptr), intent(out) :: opt + integer(ip), intent(in) :: atm(LIBCINT_ATM_SLOTS, *) + integer(ip), intent(in) :: natm + integer(ip), intent(in) :: bas(LIBCINT_BAS_SLOTS, *) + integer(ip), intent(in) :: nbas + real(dp), intent(in) :: env(*) + + call cint2e_cart_optimizer(opt, atm, natm, bas, nbas, env) + end subroutine libcint_2e_cart_optimizer + + !> Create optimizer for 2e integrals (Spherical basis) + subroutine libcint_2e_sph_optimizer(opt, atm, natm, bas, nbas, env) + type(c_ptr), intent(out) :: opt + integer(ip), intent(in) :: atm(LIBCINT_ATM_SLOTS, *) + integer(ip), intent(in) :: natm + integer(ip), intent(in) :: bas(LIBCINT_BAS_SLOTS, *) + integer(ip), intent(in) :: nbas + real(dp), intent(in) :: env(*) + + call cint2e_sph_optimizer(opt, atm, natm, bas, nbas, env) + end subroutine libcint_2e_sph_optimizer + + !> Create optimizer for ERI gradients (Cartesian basis) + subroutine libcint_2e_ip1_cart_optimizer(opt, atm, natm, bas, nbas, env) + type(c_ptr), intent(out) :: opt + integer(ip), intent(in) :: atm(LIBCINT_ATM_SLOTS, *) + integer(ip), intent(in) :: natm + integer(ip), intent(in) :: bas(LIBCINT_BAS_SLOTS, *) + integer(ip), intent(in) :: nbas + real(dp), intent(in) :: env(*) + + call cint2e_ip1_cart_optimizer(opt, atm, natm, bas, nbas, env) + end subroutine libcint_2e_ip1_cart_optimizer + + !> Create optimizer for ERI gradients (Spherical basis) + subroutine libcint_2e_ip1_sph_optimizer(opt, atm, natm, bas, nbas, env) + type(c_ptr), intent(out) :: opt + integer(ip), intent(in) :: atm(LIBCINT_ATM_SLOTS, *) + integer(ip), intent(in) :: natm + integer(ip), intent(in) :: bas(LIBCINT_BAS_SLOTS, *) + integer(ip), intent(in) :: nbas + real(dp), intent(in) :: env(*) + + call cint2e_ip1_sph_optimizer(opt, atm, natm, bas, nbas, env) + end subroutine libcint_2e_ip1_sph_optimizer + + !> Create optimizer for spinor 2e integrals + subroutine libcint_2e_spsp1_optimizer(opt, atm, natm, bas, nbas, env) + type(c_ptr), intent(out) :: opt + integer(ip), intent(in) :: atm(LIBCINT_ATM_SLOTS, *) + integer(ip), intent(in) :: natm + integer(ip), intent(in) :: bas(LIBCINT_BAS_SLOTS, *) + integer(ip), intent(in) :: nbas + real(dp), intent(in) :: env(*) + + call cint2e_spsp1_optimizer(opt, atm, natm, bas, nbas, env) + end subroutine libcint_2e_spsp1_optimizer + + !> Delete/free an optimizer and nullify the pointer + subroutine libcint_del_optimizer(opt) + type(c_ptr), intent(inout) :: opt + + call CINTdel_optimizer(opt) + opt = c_null_ptr ! Nullify to prevent use-after-free + end subroutine libcint_del_optimizer + +end module libcint_fortran diff --git a/include/libcint_interface.f90 b/include/libcint_interface.f90 new file mode 100644 index 00000000..388a3162 --- /dev/null +++ b/include/libcint_interface.f90 @@ -0,0 +1,431 @@ +! +! Low-level Fortran interface for libcint using iso_c_binding +! +! This module provides type-safe, direct bindings to libcint C functions. +! All types use C-compatible kinds (c_int, c_double, c_double_complex). +! +! For a higher-level API with native Fortran types, see libcint_fortran module. +! +module libcint_interface + use iso_c_binding + implicit none + private + + ! ======================================================================== + ! Public API + ! ======================================================================== + + ! Constants + public :: CHARGE_OF, PTR_COORD, NUC_MOD_OF, PTR_ZETA, ATM_SLOTS + public :: ATOM_OF, ANG_OF, NPRIM_OF, NCTR_OF, KAPPA_OF, PTR_EXP, PTR_COEFF, BAS_SLOTS + public :: PTR_EXPCUTOFF, PTR_COMMON_ORIG, PTR_RINV_ORIG, PTR_RINV_ZETA + public :: PTR_RANGE_OMEGA, PTR_ENV_START + + ! Dimension and normalization functions + public :: CINTcgto_cart, CINTcgto_spheric, CINTcgto_spinor + public :: CINTtot_cgto_spheric, CINTtot_pgto_spheric + public :: CINTgto_norm + + ! One-electron integrals + public :: cint1e_ovlp_cart, cint1e_nuc_cart, cint1e_kin_cart, cint1e_ipovlp_cart + public :: cint1e_ovlp_sph, cint1e_nuc_sph, cint1e_kin_sph, cint1e_ipovlp_sph + public :: cint1e_spnucsp + + ! Two-electron integrals + public :: cint2e_cart, cint2e_sph + public :: cint2e_ip1_cart, cint2e_ip1_sph + public :: cint2e_spsp1 + + ! Optimizers + public :: cint2e_cart_optimizer, cint2e_sph_optimizer + public :: cint2e_ip1_cart_optimizer, cint2e_ip1_sph_optimizer + public :: cint2e_spsp1_optimizer + public :: CINTdel_optimizer + + ! ======================================================================== + ! Constants from cint.h + ! ======================================================================== + + ! Slots of atm array (1-indexed in Fortran, but values are 0-indexed offsets for C) + integer(c_int), parameter :: CHARGE_OF = 1 + integer(c_int), parameter :: PTR_COORD = 2 + integer(c_int), parameter :: NUC_MOD_OF = 3 + integer(c_int), parameter :: PTR_ZETA = 4 + integer(c_int), parameter :: ATM_SLOTS = 6 + + ! Slots of bas array + integer(c_int), parameter :: ATOM_OF = 1 + integer(c_int), parameter :: ANG_OF = 2 + integer(c_int), parameter :: NPRIM_OF = 3 + integer(c_int), parameter :: NCTR_OF = 4 + integer(c_int), parameter :: KAPPA_OF = 5 + integer(c_int), parameter :: PTR_EXP = 6 + integer(c_int), parameter :: PTR_COEFF = 7 + integer(c_int), parameter :: BAS_SLOTS = 8 + + ! Global parameters in env array + integer(c_int), parameter :: PTR_EXPCUTOFF = 0 + integer(c_int), parameter :: PTR_COMMON_ORIG = 1 + integer(c_int), parameter :: PTR_RINV_ORIG = 4 + integer(c_int), parameter :: PTR_RINV_ZETA = 7 + integer(c_int), parameter :: PTR_RANGE_OMEGA = 8 + integer(c_int), parameter :: PTR_ENV_START = 20 + + ! ======================================================================== + ! Interface to C functions + ! ======================================================================== + + interface + + ! ==================================================================== + ! Basis set dimension functions + ! ==================================================================== + + ! Get number of Cartesian GTOs for a given shell + function CINTcgto_cart(bas_id, bas) bind(C, name='CINTcgto_cart') + import :: c_int + integer(c_int), value, intent(in) :: bas_id + integer(c_int), intent(in) :: bas(*) + integer(c_int) :: CINTcgto_cart + end function CINTcgto_cart + + ! Get number of spherical GTOs for a given shell + function CINTcgto_spheric(bas_id, bas) bind(C, name='CINTcgto_spheric') + import :: c_int + integer(c_int), value, intent(in) :: bas_id + integer(c_int), intent(in) :: bas(*) + integer(c_int) :: CINTcgto_spheric + end function CINTcgto_spheric + + ! Get number of spinor components for a given shell + function CINTcgto_spinor(bas_id, bas) bind(C, name='CINTcgto_spinor') + import :: c_int + integer(c_int), value, intent(in) :: bas_id + integer(c_int), intent(in) :: bas(*) + integer(c_int) :: CINTcgto_spinor + end function CINTcgto_spinor + + ! ==================================================================== + ! GTO normalization and total counts + ! ==================================================================== + + ! Compute GTO normalization factor + function CINTgto_norm(n, a) bind(C, name='CINTgto_norm') + import :: c_int, c_double + integer(c_int), value, intent(in) :: n + real(c_double), value, intent(in) :: a + real(c_double) :: CINTgto_norm + end function CINTgto_norm + + ! Total number of contracted spherical GTOs + function CINTtot_cgto_spheric(bas, nbas) bind(C, name='CINTtot_cgto_spheric') + import :: c_int + integer(c_int), intent(in) :: bas(*) + integer(c_int), value, intent(in) :: nbas + integer(c_int) :: CINTtot_cgto_spheric + end function CINTtot_cgto_spheric + + ! Total number of primitive spherical GTOs + function CINTtot_pgto_spheric(bas, nbas) bind(C, name='CINTtot_pgto_spheric') + import :: c_int + integer(c_int), intent(in) :: bas(*) + integer(c_int), value, intent(in) :: nbas + integer(c_int) :: CINTtot_pgto_spheric + end function CINTtot_pgto_spheric + + ! ==================================================================== + ! Optimizer functions + ! ==================================================================== + + ! Delete/free an optimizer object + subroutine CINTdel_optimizer(opt) bind(C, name='CINTdel_optimizer') + import :: c_ptr + type(c_ptr), intent(inout) :: opt + end subroutine CINTdel_optimizer + + ! ==================================================================== + ! One-electron integrals - Cartesian + ! ==================================================================== + + ! Overlap integral (Cartesian) + function cint1e_ovlp_cart(buf, shls, atm, natm, bas, nbas, env) & + bind(C, name='cint1e_ovlp_cart') + import :: c_double, c_int + real(c_double), intent(out) :: buf(*) + integer(c_int), intent(in) :: shls(*) + integer(c_int), intent(in) :: atm(*) + integer(c_int), value, intent(in) :: natm + integer(c_int), intent(in) :: bas(*) + integer(c_int), value, intent(in) :: nbas + real(c_double), intent(in) :: env(*) + integer(c_int) :: cint1e_ovlp_cart + end function cint1e_ovlp_cart + + ! Nuclear attraction integral (Cartesian) + function cint1e_nuc_cart(buf, shls, atm, natm, bas, nbas, env) & + bind(C, name='cint1e_nuc_cart') + import :: c_double, c_int + real(c_double), intent(out) :: buf(*) + integer(c_int), intent(in) :: shls(*) + integer(c_int), intent(in) :: atm(*) + integer(c_int), value, intent(in) :: natm + integer(c_int), intent(in) :: bas(*) + integer(c_int), value, intent(in) :: nbas + real(c_double), intent(in) :: env(*) + integer(c_int) :: cint1e_nuc_cart + end function cint1e_nuc_cart + + ! Kinetic energy integral (Cartesian) + function cint1e_kin_cart(buf, shls, atm, natm, bas, nbas, env) & + bind(C, name='cint1e_kin_cart') + import :: c_double, c_int + real(c_double), intent(out) :: buf(*) + integer(c_int), intent(in) :: shls(*) + integer(c_int), intent(in) :: atm(*) + integer(c_int), value, intent(in) :: natm + integer(c_int), intent(in) :: bas(*) + integer(c_int), value, intent(in) :: nbas + real(c_double), intent(in) :: env(*) + integer(c_int) :: cint1e_kin_cart + end function cint1e_kin_cart + + ! Overlap gradient integral (Cartesian) + function cint1e_ipovlp_cart(buf, shls, atm, natm, bas, nbas, env) & + bind(C, name='cint1e_ipovlp_cart') + import :: c_double, c_int + real(c_double), intent(out) :: buf(*) + integer(c_int), intent(in) :: shls(*) + integer(c_int), intent(in) :: atm(*) + integer(c_int), value, intent(in) :: natm + integer(c_int), intent(in) :: bas(*) + integer(c_int), value, intent(in) :: nbas + real(c_double), intent(in) :: env(*) + integer(c_int) :: cint1e_ipovlp_cart + end function cint1e_ipovlp_cart + + ! ==================================================================== + ! One-electron integrals - Spherical + ! ==================================================================== + + ! Overlap integral (Spherical) + function cint1e_ovlp_sph(buf, shls, atm, natm, bas, nbas, env) & + bind(C, name='cint1e_ovlp_sph') + import :: c_double, c_int + real(c_double), intent(out) :: buf(*) + integer(c_int), intent(in) :: shls(*) + integer(c_int), intent(in) :: atm(*) + integer(c_int), value, intent(in) :: natm + integer(c_int), intent(in) :: bas(*) + integer(c_int), value, intent(in) :: nbas + real(c_double), intent(in) :: env(*) + integer(c_int) :: cint1e_ovlp_sph + end function cint1e_ovlp_sph + + ! Nuclear attraction integral (Spherical) + function cint1e_nuc_sph(buf, shls, atm, natm, bas, nbas, env) & + bind(C, name='cint1e_nuc_sph') + import :: c_double, c_int + real(c_double), intent(out) :: buf(*) + integer(c_int), intent(in) :: shls(*) + integer(c_int), intent(in) :: atm(*) + integer(c_int), value, intent(in) :: natm + integer(c_int), intent(in) :: bas(*) + integer(c_int), value, intent(in) :: nbas + real(c_double), intent(in) :: env(*) + integer(c_int) :: cint1e_nuc_sph + end function cint1e_nuc_sph + + ! Kinetic energy integral (Spherical) + function cint1e_kin_sph(buf, shls, atm, natm, bas, nbas, env) & + bind(C, name='cint1e_kin_sph') + import :: c_double, c_int + real(c_double), intent(out) :: buf(*) + integer(c_int), intent(in) :: shls(*) + integer(c_int), intent(in) :: atm(*) + integer(c_int), value, intent(in) :: natm + integer(c_int), intent(in) :: bas(*) + integer(c_int), value, intent(in) :: nbas + real(c_double), intent(in) :: env(*) + integer(c_int) :: cint1e_kin_sph + end function cint1e_kin_sph + + ! Overlap gradient integral (Spherical) + function cint1e_ipovlp_sph(buf, shls, atm, natm, bas, nbas, env) & + bind(C, name='cint1e_ipovlp_sph') + import :: c_double, c_int + real(c_double), intent(out) :: buf(*) + integer(c_int), intent(in) :: shls(*) + integer(c_int), intent(in) :: atm(*) + integer(c_int), value, intent(in) :: natm + integer(c_int), intent(in) :: bas(*) + integer(c_int), value, intent(in) :: nbas + real(c_double), intent(in) :: env(*) + integer(c_int) :: cint1e_ipovlp_sph + end function cint1e_ipovlp_sph + + ! ==================================================================== + ! One-electron integrals - Spinor + ! ==================================================================== + + ! Spinor nuclear attraction integral + subroutine cint1e_spnucsp(buf, shls, atm, natm, bas, nbas, env) & + bind(C, name='cint1e_spnucsp') + import :: c_double_complex, c_int, c_double + complex(c_double_complex), intent(out) :: buf(*) + integer(c_int), intent(in) :: shls(*) + integer(c_int), intent(in) :: atm(*) + integer(c_int), value, intent(in) :: natm + integer(c_int), intent(in) :: bas(*) + integer(c_int), value, intent(in) :: nbas + real(c_double), intent(in) :: env(*) + end subroutine cint1e_spnucsp + + ! ==================================================================== + ! Two-electron integrals - Cartesian + ! ==================================================================== + + ! Electron repulsion integral (Cartesian) + function cint2e_cart(buf, shls, atm, natm, bas, nbas, env, opt) & + bind(C, name='cint2e_cart') + import :: c_double, c_int, c_ptr + real(c_double), intent(out) :: buf(*) + integer(c_int), intent(in) :: shls(*) + integer(c_int), intent(in) :: atm(*) + integer(c_int), value, intent(in) :: natm + integer(c_int), intent(in) :: bas(*) + integer(c_int), value, intent(in) :: nbas + real(c_double), intent(in) :: env(*) + type(c_ptr), value, intent(in) :: opt + integer(c_int) :: cint2e_cart + end function cint2e_cart + + ! Create optimizer for 2e integrals (Cartesian) + subroutine cint2e_cart_optimizer(opt, atm, natm, bas, nbas, env) & + bind(C, name='cint2e_cart_optimizer') + import :: c_ptr, c_int, c_double + type(c_ptr), intent(out) :: opt + integer(c_int), intent(in) :: atm(*) + integer(c_int), value, intent(in) :: natm + integer(c_int), intent(in) :: bas(*) + integer(c_int), value, intent(in) :: nbas + real(c_double), intent(in) :: env(*) + end subroutine cint2e_cart_optimizer + + ! ERI gradient (Cartesian) + function cint2e_ip1_cart(buf, shls, atm, natm, bas, nbas, env, opt) & + bind(C, name='cint2e_ip1_cart') + import :: c_double, c_int, c_ptr + real(c_double), intent(out) :: buf(*) + integer(c_int), intent(in) :: shls(*) + integer(c_int), intent(in) :: atm(*) + integer(c_int), value, intent(in) :: natm + integer(c_int), intent(in) :: bas(*) + integer(c_int), value, intent(in) :: nbas + real(c_double), intent(in) :: env(*) + type(c_ptr), value, intent(in) :: opt + integer(c_int) :: cint2e_ip1_cart + end function cint2e_ip1_cart + + ! Create optimizer for ERI gradients (Cartesian) + subroutine cint2e_ip1_cart_optimizer(opt, atm, natm, bas, nbas, env) & + bind(C, name='cint2e_ip1_cart_optimizer') + import :: c_ptr, c_int, c_double + type(c_ptr), intent(out) :: opt + integer(c_int), intent(in) :: atm(*) + integer(c_int), value, intent(in) :: natm + integer(c_int), intent(in) :: bas(*) + integer(c_int), value, intent(in) :: nbas + real(c_double), intent(in) :: env(*) + end subroutine cint2e_ip1_cart_optimizer + + ! ==================================================================== + ! Two-electron integrals - Spherical + ! ==================================================================== + + ! Electron repulsion integral (Spherical) + function cint2e_sph(buf, shls, atm, natm, bas, nbas, env, opt) & + bind(C, name='cint2e_sph') + import :: c_double, c_int, c_ptr + real(c_double), intent(out) :: buf(*) + integer(c_int), intent(in) :: shls(*) + integer(c_int), intent(in) :: atm(*) + integer(c_int), value, intent(in) :: natm + integer(c_int), intent(in) :: bas(*) + integer(c_int), value, intent(in) :: nbas + real(c_double), intent(in) :: env(*) + type(c_ptr), value, intent(in) :: opt + integer(c_int) :: cint2e_sph + end function cint2e_sph + + ! Create optimizer for 2e integrals (Spherical) + subroutine cint2e_sph_optimizer(opt, atm, natm, bas, nbas, env) & + bind(C, name='cint2e_sph_optimizer') + import :: c_ptr, c_int, c_double + type(c_ptr), intent(out) :: opt + integer(c_int), intent(in) :: atm(*) + integer(c_int), value, intent(in) :: natm + integer(c_int), intent(in) :: bas(*) + integer(c_int), value, intent(in) :: nbas + real(c_double), intent(in) :: env(*) + end subroutine cint2e_sph_optimizer + + ! ERI gradient (Spherical) + function cint2e_ip1_sph(buf, shls, atm, natm, bas, nbas, env, opt) & + bind(C, name='cint2e_ip1_sph') + import :: c_double, c_int, c_ptr + real(c_double), intent(out) :: buf(*) + integer(c_int), intent(in) :: shls(*) + integer(c_int), intent(in) :: atm(*) + integer(c_int), value, intent(in) :: natm + integer(c_int), intent(in) :: bas(*) + integer(c_int), value, intent(in) :: nbas + real(c_double), intent(in) :: env(*) + type(c_ptr), value, intent(in) :: opt + integer(c_int) :: cint2e_ip1_sph + end function cint2e_ip1_sph + + ! Create optimizer for ERI gradients (Spherical) + subroutine cint2e_ip1_sph_optimizer(opt, atm, natm, bas, nbas, env) & + bind(C, name='cint2e_ip1_sph_optimizer') + import :: c_ptr, c_int, c_double + type(c_ptr), intent(out) :: opt + integer(c_int), intent(in) :: atm(*) + integer(c_int), value, intent(in) :: natm + integer(c_int), intent(in) :: bas(*) + integer(c_int), value, intent(in) :: nbas + real(c_double), intent(in) :: env(*) + end subroutine cint2e_ip1_sph_optimizer + + ! ==================================================================== + ! Two-electron integrals - Spinor + ! ==================================================================== + + ! Spinor electron repulsion integral + subroutine cint2e_spsp1(buf, shls, atm, natm, bas, nbas, env, opt) & + bind(C, name='cint2e_spsp1') + import :: c_double_complex, c_int, c_ptr, c_double + complex(c_double_complex), intent(out) :: buf(*) + integer(c_int), intent(in) :: shls(*) + integer(c_int), intent(in) :: atm(*) + integer(c_int), value, intent(in) :: natm + integer(c_int), intent(in) :: bas(*) + integer(c_int), value, intent(in) :: nbas + real(c_double), intent(in) :: env(*) + type(c_ptr), value, intent(in) :: opt + end subroutine cint2e_spsp1 + + ! Create optimizer for spinor 2e integrals + subroutine cint2e_spsp1_optimizer(opt, atm, natm, bas, nbas, env) & + bind(C, name='cint2e_spsp1_optimizer') + import :: c_ptr, c_int, c_double + type(c_ptr), intent(out) :: opt + integer(c_int), intent(in) :: atm(*) + integer(c_int), value, intent(in) :: natm + integer(c_int), intent(in) :: bas(*) + integer(c_int), value, intent(in) :: nbas + real(c_double), intent(in) :: env(*) + end subroutine cint2e_spsp1_optimizer + + end interface + +end module libcint_interface