diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index 05739e3d..abd7784f 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -47,6 +47,7 @@ jobs: - name: Configure MOPAC with CMake run: | cmake -B build \ + -DCMAKE_BUILD_TYPE=Debug \ -DTHREADS_KEYWORD=OFF \ -DSTATIC_BUILD=ON @@ -54,6 +55,12 @@ jobs: run: | cmake --build build -- -j$NUM_CORES +# - name: Run Valgrind on MOPAC API tester +# run: | +# sudo apt update && sudo apt install -y valgrind +# cd build/tests +# valgrind -s --track-origins=yes --leak-check=yes ./mopac-api-test + - name: Test MOPAC with CTest run: | cd build @@ -92,6 +99,7 @@ jobs: - name: Configure MOPAC with CMake run: | cmake -B build \ + -DCMAKE_BUILD_TYPE=Debug \ -DTHREADS_KEYWORD=OFF \ -DENABLE_COVERAGE=ON @@ -99,6 +107,12 @@ jobs: run: | cmake --build build -- -j$NUM_CORES +# - name: Run Valgrind on MOPAC API tester +# run: | +# sudo apt update && sudo apt install -y valgrind +# cd build/tests +# valgrind -s --track-origins=yes --leak-check=yes ./mopac-api-test + - name: Test MOPAC with CTest run: | cd build diff --git a/CMakeLists.txt b/CMakeLists.txt index 2b1e5afa..9a5625a1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -134,6 +134,11 @@ if(BUILD_BZ) endif() endif() +# Use malloc in API w/ Intel Fortran compiler because it can't deallocate Fortran-allocated memory passed through a C interface +if(CMAKE_Fortran_COMPILER_ID STREQUAL "Intel") + target_compile_definitions(mopac-core PRIVATE MOPAC_API_MALLOC) +endif() + # MOPAC shared library ABI compatibility version set_target_properties(mopac-core PROPERTIES SOVERSION 2) @@ -151,6 +156,7 @@ if(APPLE) target_compile_definitions(mopac-core PRIVATE MOPAC_OS="MacOS") elseif(WIN32) target_compile_definitions(mopac-core PRIVATE MOPAC_OS="Windows") + target_compile_definitions(mopac-core PRIVATE WIN32) elseif(UNIX) target_compile_definitions(mopac-core PRIVATE MOPAC_OS="Linux") endif() @@ -262,12 +268,10 @@ if(TESTS) endif() message(FATAL_ERROR "Python3 and Numpy are required for MOPAC testing (testing can be disabled with -DTESTS=OFF)") endif() + enable_testing() add_executable(mopac-api-test) target_link_libraries(mopac-api-test mopac-core) - if (STATIC_BUILD) - target_link_options(mopac-api-test PUBLIC "-static") - endif() - enable_testing() + add_subdirectory(include) add_subdirectory(tests) endif() @@ -320,6 +324,10 @@ install(FILES "${CMAKE_SOURCE_DIR}/.github/mopac.ico" DESTINATION "." COMPONENT install(FILES "${CMAKE_SOURCE_DIR}/CITATION.cff" "${CMAKE_SOURCE_DIR}/COPYING" "${CMAKE_SOURCE_DIR}/COPYING.lesser" DESTINATION "." COMPONENT redist EXCLUDE_FROM_ALL) +install(FILES "${CMAKE_SOURCE_DIR}/include/mopac_api.h" DESTINATION ${CMAKE_INSTALL_INCLUDEDIR} COMPONENT redist EXCLUDE_FROM_ALL) +install(FILES "${CMAKE_SOURCE_DIR}/include/mopac_api_internal.F90" DESTINATION ${CMAKE_INSTALL_INCLUDEDIR} COMPONENT redist EXCLUDE_FROM_ALL) +install(FILES "${CMAKE_SOURCE_DIR}/include/mopac_api_f.F90" DESTINATION ${CMAKE_INSTALL_INCLUDEDIR} COMPONENT redist EXCLUDE_FROM_ALL) +install(FILES "${CMAKE_SOURCE_DIR}/include/LICENSE" DESTINATION ${CMAKE_INSTALL_INCLUDEDIR} COMPONENT redist EXCLUDE_FROM_ALL) # Install the executables and library install(TARGETS mopac RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR} COMPONENT main) diff --git a/include/CMakeLists.txt b/include/CMakeLists.txt new file mode 100644 index 00000000..1ab8df15 --- /dev/null +++ b/include/CMakeLists.txt @@ -0,0 +1,18 @@ +# Molecular Orbital PACkage (MOPAC) +# Copyright (C) 2021, Virginia Polytechnic Institute and State University +# +# MOPAC is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# MOPAC is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with this program. If not, see . + +target_sources(mopac-api-test PRIVATE mopac_api_f.F90) +target_sources(mopac-api-test PRIVATE mopac_api_internal.F90) diff --git a/include/LICENSE b/include/LICENSE new file mode 100644 index 00000000..1a36f5fe --- /dev/null +++ b/include/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2021, Virginia Polytechnic Institute and State University + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/include/mopac_api.h b/include/mopac_api.h new file mode 100644 index 00000000..0bb453ed --- /dev/null +++ b/include/mopac_api.h @@ -0,0 +1,176 @@ +/* Molecular Orbital PACkage (MOPAC) + * Copyright (C) 2021, Virginia Polytechnic Institute and State University + */ + +/* Diskless/stateless Application Programming Interface (API) to core MOPAC operations */ +#ifndef MOPAC_API_H +#define MOPAC_API_H + +/* data that defines the atomistic system and MOPAC job options */ +struct mopac_system { + /* number of atoms */ + int natom; + /* number of atoms that are allowed to move (first natom_move atoms in array) */ + int natom_move; + /* net charge */ + int charge; + /* number of spin excitations, floor[(number of alpha electrons)/2 - (number of beta electrons)/2] */ + int spin; + /* semiempirical model: PM7 = 0, PM6-D3H4 = 1, PM6-ORG = 2, PM6 = 3, AM1 = 4, RM1 = 5 */ + int model; + /* dielectric constant for COSMO implicit solvent, must be 1 (no solvent) for nlattice > 0 */ + double epsilon; + /* atomic number of each atom */ + int *atom; /* [natom] */ + /* (x,y,z) coordinates of each atom (Angstroms) */ + double *coord; /* [3*natom] */ + /* number of lattice vectors / translation vectors / periodic dimensions */ + int nlattice; + /* number of lattice vectors that are allowed to move (first nlattice_move vectors in array) */ + int nlattice_move; + /* external hydrostatic pressure (Gigapascals) */ + double pressure; + /* (x,y,z) coordinates of each lattice vectors (Angstroms) */ + double *lattice; /* [3*nlattice] */ + /* relative numerical tolerances (equivalent to GNORM and RELSCF keyword values) */ + double tolerance; + /* time limit for a MOPAC calculation (seconds) */ + int max_time; +}; + +/* calculated ground-state properties of an atomistic system and MOPAC job info */ +struct mopac_properties { + /* heat of formation (kcal/mol) */ + double heat; + /* dipole moment vector (Debye) */ + double dipole[3]; + /* atomic partial charges */ + double *charge; /* [natom] */ + /* (x,y,z) coordinates of each moveable atom (Angstroms) */ + double *coord_update; /* [3*natom_move] */ + /* (x,y,z) heat gradients for each moveable atom (kcal/mol/Angstrom) */ + double *coord_deriv; /* [3*natom_move] */ + /* vibrational frequencies of normal modes (1/cm) */ + double *freq; /* [3*natom_move], NULL if unavailable */ + /* (x,y,z) displacement vectors of normal modes */ + double *disp; /* [3*natom_move,3*natom_move], NULL if unavailable */ + /* bond-order matrix in compressed sparse column (CSC) matrix format (0-based indexing) + * with insignificant bond orders (<0.01) truncated + * diagonal matrix entries are atomic valencies */ + /* > first index of each atom in CSC bond-order matrix */ + int *bond_index; /* [natom+1] */ + /* > list of atoms bonded to each atom in CSC format */ + int *bond_atom; /* [bond_index[natom]] */ + /* > bond order of atoms bonded to each atom in CSC format */ + double *bond_order; /* [bond_index[natom]] */ + /* (x,y,z) coordinates of each moveable lattice vectors (Angstroms) */ + double *lattice_update; /* [3*nlattice_move] */ + /* (x,y,z) heat gradients for each moveable lattice vector (kcal/mol/Angstrom) */ + double *lattice_deriv; /* [3*nlattice_move] */ + /* stress tensor (Gigapascals) in Voigt form (xx, yy, zz, yz, xz, xy) */ + double stress[6]; /* 0 if unavailable */ + /* number of MOPAC error messages (negative value indicates that allocation of error_msg failed) */ + int nerror; + /* text of MOPAC error messages */ + char **error_msg; /* [nerror][*] */ +}; + +/* data that describes the electronic state using standard molecular orbitals */ +struct mopac_state { + /* MOPAC data format is adapted from molkst_C and Common_arrays_C modules */ + /* > number of matrix elements in packed lower triangle matrix format */ + int mpack; /* 0 if state is unavailable */ + /* > flag for unrestricted Hartree-Fock ground state (0 == restricted, 1 == unrestricted) */ + int uhf; + /* > alpha density matrix */ + double *pa; /* [mpack] */ + /* > beta density matrix */ + double *pb; /* [mpack], NULL if uhf == 0 */ +}; + +/* data that describes the electronic state using localized molecular orbitals */ +struct mozyme_state { + /* MOZYME data format is adapted from molkst_C, Common_arrays_C, and MOZYME_C modules */ + /* > number of real atoms */ + int numat; /* 0 if state is unavailable */ + /* > number of Lewis bonds per real atom */ + int *nbonds; /* [numat] */ + /* > list of Lewis-bonded real atoms for each real atom */ + int *ibonds; /* [9,numat] */ + /* > number of orbitals per real atom */ + int *iorbs; /* [numat] */ + /* > number of occupied molecular orbitals */ + int noccupied; + /* > number of atoms in each occupied LMO */ + int *ncf; /* [noccupied] */ + /* > number of virtual molecular orbitals */ + int nvirtual; + /* > number of atoms in each virtual LMO */ + int *nce; /* [nvirtual] */ + /* > size of array icocc */ + int icocc_dim; + /* > index of each real atom in the occupied LMOs */ + int *icocc; /* [iccoc_dim] */ + /* > size of array icvir */ + int icvir_dim; + /* > index of each real atom in the virtual LMOs */ + int *icvir; /* [icvir_dim] */ + /* > size of array cocc */ + int cocc_dim; + /* > atomic orbital coefficients of the occupied LMOs */ + double *cocc; /* [cocc_dim] */ + /* > size of array cvir */ + int cvir_dim; + /* > atomic orbital coefficients of the virtual LMOs */ + double *cvir; /* [cvir_dim] */ +}; + +/* MOPAC electronic ground state calculation */ +void mopac_scf(struct mopac_system *system, + struct mopac_state *state, + struct mopac_properties *properties); + +/* MOPAC geometry relaxation */ +void mopac_relax(struct mopac_system *system, + struct mopac_state *state, + struct mopac_properties *properties); + +/* MOPAC vibrational calculation */ +void mopac_vibe(struct mopac_system *system, + struct mopac_state *state, + struct mopac_properties *properties); + +/* MOZYME electronic ground state calculation */ +void mozyme_scf(struct mopac_system *system, + struct mozyme_state *state, + struct mopac_properties *properties); + +/* MOZYME geometry relaxation */ +void mozyme_relax(struct mopac_system *system, + struct mozyme_state *state, + struct mopac_properties *properties); + +/* MOZYME vibrational calculation */ +void mozyme_vibe(struct mopac_system *system, + struct mozyme_state *state, + struct mopac_properties *properties); + +/* allocate memory for mopac_state */ +void create_mopac_state(struct mopac_state *state); + +/* allocate memory for mozyme_state */ +void create_mozyme_state(struct mozyme_state *state); + +/* deallocate memory in mopac_properties */ +void destroy_mopac_properties(struct mopac_properties *properties); + +/* deallocate memory in mopac_state */ +void destroy_mopac_state(struct mopac_state *state); + +/* deallocate memory in mozyme_state */ +void destroy_mozyme_state(struct mozyme_state *state); + +/* run MOPAC conventionally from an input file */ +int run_mopac_from_input(char *path_to_file); + +#endif \ No newline at end of file diff --git a/include/mopac_api_f.F90 b/include/mopac_api_f.F90 new file mode 100644 index 00000000..d05d95de --- /dev/null +++ b/include/mopac_api_f.F90 @@ -0,0 +1,188 @@ +! Molecular Orbital PACkage (MOPAC) +! Copyright (C) 2021, Virginia Polytechnic Institute and State University + +! Fortran wrapper for diskless/stateless Application Programming Interface (API) to core MOPAC operations +module mopac_api_f + implicit none + + private + ! public derived types of the MOPAC API + public :: mopac_system_f, mopac_properties_f, mopac_state_f, mozyme_state_f + ! public subroutines of the MOPAC API + public :: mopac_scf_f, mopac_relax_f, mopac_vibe_f, mozyme_scf_f, mozyme_relax_f, mozyme_vibe_f + ! public subroutine of the simple, disk-based MOPAC API + public :: run_mopac_from_input_f + + ! data that defines the atomistic system and MOPAC job options (Fortran wrapper) + type :: mopac_system_f + ! number of atoms + integer :: natom = 0 + ! number of atoms that are allowed to move (first natom_move atoms in array) + integer :: natom_move = 0 + ! net charge + integer :: charge = 0 + ! number of spin excitations, floor[(number of alpha electrons)/2 - (number of beta electrons)/2] + integer :: spin = 0 + ! semiempirical model: PM7 = 0, PM6-D3H4 = 1, PM6-ORG = 2, PM6 = 3, AM1 = 4, RM1 = 5 + integer :: model = 0 + ! dielectric constant for COSMO implicit solvent, must be 1 (no solvent) for nlattice > 0 + double precision :: epsilon = 1.d0 + ! atomic number of each atom [natom] + integer, dimension (:), allocatable :: atom + ! (x,y,z) coordinates of each atom (Angstroms) [3*natom] + double precision, dimension (:), allocatable :: coord + ! number of lattice vectors / translation vectors / periodic dimensions + integer :: nlattice = 0 + ! number of lattice vectors that are allowed to move (first nlattice_move vectors in array) + integer :: nlattice_move = 0 + ! external hydrostatic pressure (Gigapascals) + double precision :: pressure = 0.d0 + ! (x,y,z) coordinates of each lattice vectors (Angstroms) [3*nlattice] + double precision, dimension (:), allocatable :: lattice + ! relative numerical tolerances (equivalent to GNORM and RELSCF keyword values) + double precision :: tolerance = 1.d0 + ! time limit for a MOPAC calculation (seconds) + integer :: max_time = 3600 + end type + + ! calculated ground-state properties of an atomistic system and MOPAC job info (Fortran wrapper) + type :: mopac_properties_f + ! heat of formation (kcal/mol) + double precision :: heat + ! dipole moment vector (Debye) + double precision, dimension (3) :: dipole + ! atomic partial charges [natom] + double precision, dimension (:), allocatable :: charge + ! (x,y,z) coordinates of each moveable atom (Angstroms) [3*natom_move] + double precision, dimension (:), allocatable :: coord_update + ! (x,y,z) heat gradients for each moveable atom (kcal/mol/Angstrom) [3*natom_move] + double precision, dimension (:), allocatable :: coord_deriv + ! vibrational frequencies of normal modes (1/cm) [3*natom_move], if available + double precision, dimension (:), allocatable :: freq + ! (x,y,z) displacement vectors of normal modes [3*natom_move,3*natom_move], if available + double precision, dimension (:,:), allocatable :: disp + ! bond-order matrix in compressed sparse column (CSC) matrix format + ! with insignificant bond orders (<0.01) truncated + ! diagonal matrix entries are atomic valencies + ! > first index of each atom in CSC bond-order matrix [natom+1] + integer, dimension (:), allocatable :: bond_index + ! > list of atoms bonded to each atom in CSC format [bond_index(natom+1)-1] + integer, dimension (:), allocatable :: bond_atom + ! > bond order of atoms bonded to each atom in CSC format [bond_index(natom+1)-1] + double precision, dimension (:), allocatable :: bond_order + ! (x,y,z) coordinates of each moveable lattice vectors (Angstroms) [3*nlattice_move] + double precision, dimension (:), allocatable :: lattice_update + ! (x,y,z) heat gradients for each moveable lattice vector (kcal/mol/Angstrom) [3*nlattice_move] + double precision, dimension (:), allocatable :: lattice_deriv + ! stress tensor (Gigapascals) in Voigt form (xx, yy, zz, yz, xz, xy), 0 if available + double precision, dimension (6) :: stress + ! number of MOPAC error messages (negative value indicates that allocation of error_msg failed) + integer :: nerror + ! text of MOPAC error messages [nerror,120] + character(len=120), dimension (:), allocatable :: error_msg + end type + + ! data that describes the electronic state using standard molecular orbitals (Fortran wrapper) + type :: mopac_state_f + ! MOPAC data format is adapted from molkst_C and Common_arrays_C modules + ! > number of matrix elements in packed lower triangle matrix format + integer :: mpack = 0 ! 0 if state is unavailable + ! > flag for unrestricted Hartree-Fock ground state (0 == restricted, 1 == unrestricted) + logical :: uhf + ! > alpha density matrix [mpack] + double precision, dimension (:), allocatable :: pa + ! > beta density matrix [mpack], empty if uhf == 0 + double precision, dimension (:), allocatable :: pb + end type + + ! data that describes the electronic state using localized molecular orbitals (Fortran wrapper) + type :: mozyme_state_f + ! MOZYME data format is adapted from molkst_C, Common_arrays_C, and MOZYME_C modules + ! > number of real atoms + integer :: numat = 0 ! 0 if state is unavailable + ! > number of Lewis bonds per real atom [numat] + integer, dimension (:), allocatable :: nbonds + ! > list of Lewis-bonded real atoms for each real atom [9,numat] + integer, dimension (:,:), allocatable :: ibonds + ! > number of orbitals per real atom [numat] + integer, dimension (:), allocatable :: iorbs + ! > number of occupied molecular orbitals + integer :: noccupied + ! > number of atoms in each occupied LMO [noccupied] + integer, dimension (:), allocatable :: ncf + ! > number of virtual molecular orbitals + integer :: nvirtual + ! > number of atoms in each virtual LMO [nvirtual] + integer, dimension (:), allocatable :: nce + ! > size of array icocc + integer :: icocc_dim + ! > index of each real atom in the occupied LMOs [iccoc_dim] + integer, dimension (:), allocatable :: icocc + ! > size of array icvir + integer :: icvir_dim + ! > index of each real atom in the virtual LMOs [icvir_dim] + integer, dimension (:), allocatable :: icvir + ! > size of array cocc + integer :: cocc_dim + ! > atomic orbital coefficients of the occupied LMOs [cocc_dim] + double precision, dimension (:), allocatable :: cocc + ! > size of array cvir + integer :: cvir_dim + ! > atomic orbital coefficients of the virtual LMOs [cvir_dim] + double precision, dimension (:), allocatable :: cvir + end type + + ! library subroutines + interface + + ! MOPAC electronic ground state calculation + module subroutine mopac_scf_f(system, state, properties) + type(mopac_system_f), intent(in) :: system + type(mopac_state_f), intent(inout) :: state + type(mopac_properties_f), intent(out) :: properties + end subroutine mopac_scf_f + + ! MOPAC geometry relaxation + module subroutine mopac_relax_f(system, state, properties) + type(mopac_system_f), intent(in) :: system + type(mopac_state_f), intent(inout) :: state + type(mopac_properties_f), intent(out) :: properties + end subroutine mopac_relax_f + + ! MOPAC vibrational calculation + module subroutine mopac_vibe_f(system, state, properties) + type(mopac_system_f), intent(in) :: system + type(mopac_state_f), intent(inout) :: state + type(mopac_properties_f), intent(out) :: properties + end subroutine mopac_vibe_f + + ! MOZYME electronic ground state calculation + module subroutine mozyme_scf_f(system, state, properties) + type(mopac_system_f), intent(in) :: system + type(mozyme_state_f), intent(inout) :: state + type(mopac_properties_f), intent(out) :: properties + end subroutine mozyme_scf_f + + ! MOZYME geometry relaxation + module subroutine mozyme_relax_f(system, state, properties) + type(mopac_system_f), intent(in) :: system + type(mozyme_state_f), intent(inout) :: state + type(mopac_properties_f), intent(out) :: properties + end subroutine mozyme_relax_f + + ! MOZYME vibrational calculation + module subroutine mozyme_vibe_f(system, state, properties) + type(mopac_system_f), intent(in) :: system + type(mozyme_state_f), intent(inout) :: state + type(mopac_properties_f), intent(out) :: properties + end subroutine mozyme_vibe_f + + ! run MOPAC conventionally from an input file + module function run_mopac_from_input_f(path_to_file) + logical :: run_mopac_from_input_f + character(len=240), intent(in) :: path_to_file + end function run_mopac_from_input_f + + end interface + +end module mopac_api_f diff --git a/include/mopac_api_internal.F90 b/include/mopac_api_internal.F90 new file mode 100644 index 00000000..9d2814b0 --- /dev/null +++ b/include/mopac_api_internal.F90 @@ -0,0 +1,610 @@ +! Molecular Orbital PACkage (MOPAC) +! Copyright (C) 2021, Virginia Polytechnic Institute and State University + +! C-bound types for diskless/stateless Application Programming Interface (API) to core MOPAC operations +module mopac_api_c + use, intrinsic :: iso_c_binding + implicit none + + ! data that defines the atomistic system and MOPAC job options + type, bind(c) :: mopac_system + ! number of atoms + integer(c_int) :: natom + ! number of atoms that are allowed to move (first natom_move atoms in array) + integer(c_int) :: natom_move + ! net charge + integer(c_int) :: charge + ! number of spin excitations, floor[(number of alpha electrons)/2 - (number of beta electrons)/2] + integer(c_int) :: spin + ! semiempirical model: PM7 = 0, PM6-D3H4 = 1, PM6-ORG = 2, PM6 = 3, AM1 = 4, RM1 = 5 + integer(c_int) :: model + ! dielectric constant for COSMO implicit solvent, must be 1 (no solvent) for nlattice > 0 + real(c_double) :: epsilon + ! atomic number of each atom + type(c_ptr) :: atom ! integer(c_int)[natom] + ! (x,y,z) coordinates of each atom (Angstroms) + type(c_ptr) :: coord ! real(c_double)[3*natom] + ! number of lattice vectors / translation vectors / periodic dimensions + integer(c_int) :: nlattice + ! number of lattice vectors that are allowed to move (first nlattice_move vectors in array) + integer(c_int) :: nlattice_move + ! external hydrostatic pressure (Gigapascals) + real(c_double) :: pressure + ! (x,y,z) coordinates of each lattice vectors (Angstroms) + type(c_ptr) :: lattice ! real(c_double)[3*nlattice] + ! relative numerical tolerances (equivalent to GNORM and RELSCF keyword values) + real(c_double) :: tolerance + ! time limit for a MOPAC calculation (seconds) + integer(c_int) :: max_time + end type + + ! calculated ground-state properties of an atomistic system and MOPAC job info + type, bind(c) :: mopac_properties + ! heat of formation (kcal/mol) + real(c_double) :: heat + ! dipole moment vector (Debye) + real(c_double), dimension (3) :: dipole + ! atomic partial charges + type(c_ptr) :: charge ! real(c_double)[natom] + ! (x,y,z) coordinates of each moveable atom (Angstroms) + type(c_ptr) :: coord_update ! real(c_double)[3*natom_move] + ! (x,y,z) heat gradients for each moveable atom (kcal/mol/Angstrom) + type(c_ptr) :: coord_deriv ! real(c_double)[3*natom_move] + ! vibrational frequencies of normal modes (1/cm) + type(c_ptr) :: freq ! real(c_double)[3*natom_move], NULL if unavailable + ! (x,y,z) displacement vectors of normal modes + type(c_ptr) :: disp ! real(c_double)[3*natom_move,3*natom_move], NULL if unavailable + ! bond-order matrix in compressed sparse column (CSC) matrix format (0-based indexing) + ! with insignificant bond orders (<0.01) truncated + ! diagonal matrix entries are atomic valencies + ! > first index of each atom in CSC bond-order matrix + type(c_ptr) :: bond_index ! integer(c_int)[natom+1] + ! > list of atoms bonded to each atom in CSC format + type(c_ptr) :: bond_atom ! integer(c_int)[bond_index(natom+1)] + ! > bond order of atoms bonded to each atom in CSC format + type(c_ptr) :: bond_order ! real(c_double)[bond_index(natom+1)] + ! (x,y,z) coordinates of each moveable lattice vectors (Angstroms) + type(c_ptr) :: lattice_update ! real(c_double)[3*nlattice_move] + ! (x,y,z) heat gradients for each moveable lattice vector (kcal/mol/Angstrom) + type(c_ptr) :: lattice_deriv ! real(c_double)[3*nlattice_move] + ! stress tensor (Gigapascals) in Voigt form (xx, yy, zz, yz, xz, xy) + real(c_double), dimension (6) :: stress ! 0 if unavailable + ! number of MOPAC error messages (negative value indicates that allocation of error_msg failed) + integer(c_int) :: nerror + ! text of MOPAC error messages + type(c_ptr) :: error_msg ! type(c_ptr)[nerror] -> character(kind=c_char)[*] + end type + + ! data that describes the electronic state using standard molecular orbitals + type, bind(c) :: mopac_state + ! MOPAC data format is adapted from molkst_C and Common_arrays_C modules + ! > number of matrix elements in packed lower triangle matrix format + integer(c_int) :: mpack ! 0 if state is unavailable + ! > flag for unrestricted Hartree-Fock ground state (0 == restricted, 1 == unrestricted) + integer(c_int) :: uhf + ! > alpha density matrix + type(c_ptr) :: pa ! real(c_double)[mpack] + ! > beta density matrix + type(c_ptr) :: pb ! real(c_double)[mpack], NULL if uhf == 0 + end type + + ! data that describes the electronic state using localized molecular orbitals + type, bind(c) :: mozyme_state + ! MOZYME data format is adapted from molkst_C, Common_arrays_C, and MOZYME_C modules + ! > number of real atoms + integer(c_int) :: numat ! 0 if state is unavailable + ! > number of Lewis bonds per real atom + type(c_ptr) :: nbonds ! integer(c_int)[numat] + ! > list of Lewis-bonded real atoms for each real atom + type(c_ptr) :: ibonds ! integer(c_int)[9,numat] + ! > number of orbitals per real atom + type(c_ptr) :: iorbs ! integer(c_int)[numat] + ! > number of occupied molecular orbitals + integer(c_int) :: noccupied + ! > number of atoms in each occupied LMO + type(c_ptr) :: ncf ! integer(c_int)[noccupied] + ! > number of virtual molecular orbitals + integer(c_int) :: nvirtual + ! > number of atoms in each virtual LMO + type(c_ptr) :: nce ! integer(c_int)[nvirtual] + ! > size of array icocc + integer(c_int) :: icocc_dim + ! > index of each real atom in the occupied LMOs + type(c_ptr) :: icocc ! integer(c_int)[iccoc_dim] + ! > size of array icvir + integer(c_int) :: icvir_dim + ! > index of each real atom in the virtual LMOs + type(c_ptr) :: icvir ! integer(c_int)[icvir_dim] + ! > size of array cocc + integer(c_int) :: cocc_dim + ! > atomic orbital coefficients of the occupied LMOs + type(c_ptr) :: cocc ! real(c_double)[cocc_dim] + ! > size of array cvir + integer(c_int) :: cvir_dim + ! > atomic orbital coefficients of the virtual LMOs + type(c_ptr) :: cvir ! real(c_double)[cvir_dim] + end type + +end module mopac_api_c + +! Implementation of Fortran API wrapper +submodule (mopac_api_f) mopac_api_f_internal + use, intrinsic :: iso_c_binding + use mopac_api_c + implicit none + + ! C-bound library subroutines + interface + + ! MOPAC electronic ground state calculation + subroutine mopac_scf(system, state, properties) bind(c) + use mopac_api_c + type(mopac_system), intent(in) :: system + type(mopac_state), intent(inout) :: state + type(mopac_properties), intent(out) :: properties + end subroutine mopac_scf + + ! MOPAC geometry relaxation + subroutine mopac_relax(system, state, properties) bind(c) + use mopac_api_c + type(mopac_system), intent(in) :: system + type(mopac_state), intent(inout) :: state + type(mopac_properties), intent(out) :: properties + end subroutine mopac_relax + + ! MOPAC vibrational calculation + subroutine mopac_vibe(system, state, properties) bind(c) + use mopac_api_c + type(mopac_system), intent(in) :: system + type(mopac_state), intent(inout) :: state + type(mopac_properties), intent(out) :: properties + end subroutine mopac_vibe + + ! MOZYME electronic ground state calculation + subroutine mozyme_scf(system, state, properties) bind(c) + use mopac_api_c + type(mopac_system), intent(in) :: system + type(mozyme_state), intent(inout) :: state + type(mopac_properties), intent(out) :: properties + end subroutine mozyme_scf + + ! MOZYME geometry relaxation + subroutine mozyme_relax(system, state, properties) bind(c) + use mopac_api_c + type(mopac_system), intent(in) :: system + type(mozyme_state), intent(inout) :: state + type(mopac_properties), intent(out) :: properties + end subroutine mozyme_relax + + ! MOZYME vibrational calculation + subroutine mozyme_vibe(system, state, properties) bind(c) + use mopac_api_c + type(mopac_system), intent(in) :: system + type(mozyme_state), intent(inout) :: state + type(mopac_properties), intent(out) :: properties + end subroutine mozyme_vibe + + ! allocate memory for mopac_state + subroutine create_mopac_state(state) bind(c) + use mopac_api_c + type(mopac_state), intent(inout) :: state + end subroutine create_mopac_state + + ! allocate memory for mozyme_state + subroutine create_mozyme_state(state) bind(c) + use mopac_api_c + type(mozyme_state), intent(inout) :: state + end subroutine create_mozyme_state + + ! deallocate memory in mopac_properties (associated system is needed for size info) + subroutine destroy_mopac_properties(properties) bind(c) + use mopac_api_c + type(mopac_properties), intent(in) :: properties + end subroutine destroy_mopac_properties + + ! deallocate memory in mopac_state + subroutine destroy_mopac_state(state) bind(c) + use mopac_api_c + type(mopac_state), intent(in) :: state + end subroutine destroy_mopac_state + + ! deallocate memory in mozyme_state + subroutine destroy_mozyme_state(state) bind(c) + use mopac_api_c + type(mozyme_state), intent(in) :: state + end subroutine destroy_mozyme_state + + ! run MOPAC conventionally from an input file + function run_mopac_from_input(path_to_file) bind(c) + use iso_c_binding + integer(c_int) :: run_mopac_from_input + character(kind=c_char), dimension(*), intent(in) :: path_to_file + end function run_mopac_from_input + + end interface + +contains + + module subroutine mopac_scf_f(system, state, properties) + type(mopac_system_f), intent(in) :: system + type(mopac_state_f), intent(inout) :: state + type(mopac_properties_f), intent(out) :: properties + type(mopac_system) :: system_c + type(mopac_state) :: state_c + type(mopac_properties) :: properties_c + integer(c_int), pointer :: iwork(:) + real(c_double), pointer :: rwork(:) + call mopac_system_f2c(system, system_c, iwork, rwork) + call mopac_state_f2c(state, state_c) + call mopac_scf(system_c, state_c, properties_c) + call mopac_state_c2f(state_c, state) + call mopac_properties_c2f(system_c, properties_c, properties) + deallocate(iwork) + deallocate(rwork) + end subroutine mopac_scf_f + + module subroutine mopac_relax_f(system, state, properties) + type(mopac_system_f), intent(in) :: system + type(mopac_state_f), intent(inout) :: state + type(mopac_properties_f), intent(out) :: properties + type(mopac_system) :: system_c + type(mopac_state) :: state_c + type(mopac_properties) :: properties_c + integer(c_int), pointer :: iwork(:) + real(c_double), pointer :: rwork(:) + call mopac_system_f2c(system, system_c, iwork, rwork) + call mopac_state_f2c(state, state_c) + call mopac_relax(system_c, state_c, properties_c) + call mopac_state_c2f(state_c, state) + call mopac_properties_c2f(system_c, properties_c, properties) + deallocate(iwork) + deallocate(rwork) + end subroutine mopac_relax_f + + module subroutine mopac_vibe_f(system, state, properties) + type(mopac_system_f), intent(in) :: system + type(mopac_state_f), intent(inout) :: state + type(mopac_properties_f), intent(out) :: properties + type(mopac_system) :: system_c + type(mopac_state) :: state_c + type(mopac_properties) :: properties_c + integer(c_int), pointer :: iwork(:) + real(c_double), pointer :: rwork(:) + call mopac_system_f2c(system, system_c, iwork, rwork) + call mopac_state_f2c(state, state_c) + call mopac_vibe(system_c, state_c, properties_c) + call mopac_state_c2f(state_c, state) + call mopac_properties_c2f(system_c, properties_c, properties) + deallocate(iwork) + deallocate(rwork) + end subroutine mopac_vibe_f + + module subroutine mozyme_scf_f(system, state, properties) + type(mopac_system_f), intent(in) :: system + type(mozyme_state_f), intent(inout) :: state + type(mopac_properties_f), intent(out) :: properties + type(mopac_system) :: system_c + type(mozyme_state) :: state_c + type(mopac_properties) :: properties_c + integer(c_int), pointer :: iwork(:) + real(c_double), pointer :: rwork(:) + call mopac_system_f2c(system, system_c, iwork, rwork) + call mozyme_state_f2c(state, state_c) + call mozyme_scf(system_c, state_c, properties_c) + call mozyme_state_c2f(state_c, state) + call mopac_properties_c2f(system_c, properties_c, properties) + deallocate(iwork) + deallocate(rwork) + end subroutine mozyme_scf_f + + module subroutine mozyme_relax_f(system, state, properties) + type(mopac_system_f), intent(in) :: system + type(mozyme_state_f), intent(inout) :: state + type(mopac_properties_f), intent(out) :: properties + type(mopac_system) :: system_c + type(mozyme_state) :: state_c + type(mopac_properties) :: properties_c + integer(c_int), pointer :: iwork(:) + real(c_double), pointer :: rwork(:) + call mopac_system_f2c(system, system_c, iwork, rwork) + call mozyme_state_f2c(state, state_c) + call mozyme_relax(system_c, state_c, properties_c) + call mozyme_state_c2f(state_c, state) + call mopac_properties_c2f(system_c, properties_c, properties) + deallocate(iwork) + deallocate(rwork) + end subroutine mozyme_relax_f + + module subroutine mozyme_vibe_f(system, state, properties) + type(mopac_system_f), intent(in) :: system + type(mozyme_state_f), intent(inout) :: state + type(mopac_properties_f), intent(out) :: properties + type(mopac_system) :: system_c + type(mozyme_state) :: state_c + type(mopac_properties) :: properties_c + integer(c_int), pointer :: iwork(:) + real(c_double), pointer :: rwork(:) + call mopac_system_f2c(system, system_c, iwork, rwork) + call mozyme_state_f2c(state, state_c) + call mozyme_vibe(system_c, state_c, properties_c) + call mozyme_state_c2f(state_c, state) + call mopac_properties_c2f(system_c, properties_c, properties) + deallocate(iwork) + deallocate(rwork) + end subroutine mozyme_vibe_f + + module function run_mopac_from_input_f(path_to_file) + logical :: run_mopac_from_input_f + character(len=240), intent(in) :: path_to_file + character(kind=c_char), allocatable :: path_to_file_c(:) + integer :: i, size, status + size = len_trim(path_to_file) + allocate(path_to_file_c(size+1), stat=status) + if (status /= 0) then + write(*,*) "ERROR: Failed to allocate memory in MOPAC API wrapper" + stop 1 + end if + do i=1, size + path_to_file_c(i) = path_to_file(i:i) + end do + path_to_file_c(size+1) = c_null_char + status = run_mopac_from_input(path_to_file_c) + if (status == 0) then + run_mopac_from_input_f = .true. + else + run_mopac_from_input_f = .false. + end if + end function run_mopac_from_input_f + + subroutine mopac_system_f2c(system_f, system_c, iwork, rwork) + type(mopac_system_f), intent(in) :: system_f + type(mopac_system), intent(out) :: system_c + integer(c_int), pointer :: iwork(:) + real(c_double), pointer :: rwork(:) + integer :: status + system_c%natom = system_f%natom + system_c%natom_move = system_f%natom_move + system_c%charge = system_f%charge + system_c%spin = system_f%spin + system_c%model = system_f%model + system_c%epsilon = system_f%epsilon + system_c%nlattice = system_f%nlattice + system_c%nlattice_move = system_f%nlattice_move + system_c%pressure = system_f%pressure + system_c%tolerance = system_f%tolerance + system_c%max_time = system_f%max_time + allocate(iwork(system_f%natom), stat=status) + if (status /= 0) then + write(*,*) "ERROR: Failed to allocate memory in MOPAC API wrapper" + stop 1 + end if + iwork = system_f%atom + system_c%atom = c_loc(iwork) + allocate(rwork(3*system_f%natom+3*system_f%nlattice), stat=status) + if (status /= 0) then + write(*,*) "ERROR: Failed to allocate memory in MOPAC API wrapper" + stop 1 + end if + rwork(:3*system_f%natom) = system_f%coord + system_c%coord = c_loc(rwork) + if (system_f%nlattice > 0) then + rwork(3*system_f%natom+1:) = system_f%lattice + system_c%lattice = c_loc(rwork(3*system_f%natom+1)) + end if + end subroutine mopac_system_f2c + + subroutine mopac_properties_c2f(system_c, properties_c, properties_f) + type(mopac_system), intent(in) :: system_c + type(mopac_properties), intent(in) :: properties_c + type(mopac_properties_f), intent(out) :: properties_f + character(kind=c_char), pointer :: cptr(:) + type(c_ptr), pointer :: pptr(:) + integer :: i, j, size, status + properties_f%nerror = properties_c%nerror + if (properties_c%nerror == 0) then + properties_f%heat = properties_c%heat + properties_f%dipole = properties_c%dipole + properties_f%stress = properties_c%stress + call copy_real(properties_f%charge, properties_c%charge, [system_c%natom]) + call copy_real(properties_f%coord_update, properties_c%coord_update, [3*system_c%natom_move]) + call copy_real(properties_f%coord_deriv, properties_c%coord_deriv, [3*system_c%natom_move]) + if (c_associated(properties_c%freq)) then + call copy_real(properties_f%freq, properties_c%freq, [3*system_c%natom_move]) + call copy_real2(properties_f%disp, properties_c%disp, [3*system_c%natom_move,3*system_c%natom_move]) + end if + call copy_int(properties_f%bond_index, properties_c%bond_index, [system_c%natom+1]) + do i=1, system_c%natom+1 + properties_f%bond_index(i) = properties_f%bond_index(i) + 1 + end do + call copy_int(properties_f%bond_atom, properties_c%bond_atom, & + [properties_f%bond_index(system_c%natom+1)-1]) + call copy_real(properties_f%bond_order, properties_c%bond_order, & + [properties_f%bond_index(system_c%natom+1)-1]) + call copy_real(properties_f%lattice_update, properties_c%lattice_update, [3*system_c%nlattice_move]) + call copy_real(properties_f%lattice_deriv, properties_c%lattice_deriv, [3*system_c%nlattice_move]) + else if (properties_c%nerror > 0) then + allocate(properties_f%error_msg(properties_c%nerror), stat=status) + if (status /= 0) then + write(*,*) "ERROR: Failed to allocate memory in MOPAC API wrapper" + stop 1 + end if + call c_f_pointer(properties_c%error_msg, pptr, [properties_c%nerror]) + do i=1, properties_c%nerror + properties_f%error_msg(i) = ' ' + call c_f_pointer(pptr(i), cptr, [120]) + do size=1, 120 + if (cptr(size) == c_null_char) exit + end do + size = size - 1 + do j=1, size + properties_f%error_msg(i)(j:j) = cptr(j) + end do + end do + end if + call destroy_mopac_properties(properties_c) + end subroutine mopac_properties_c2f + + subroutine mopac_state_f2c(state_f, state_c) + type(mopac_state_f), intent(in) :: state_f + type(mopac_state), intent(out) :: state_c + real(c_double), pointer :: rptr(:) + state_c%mpack = state_f%mpack + if (state_f%uhf) then + state_c%uhf = 1 + else + state_c%uhf = 0 + end if + call create_mopac_state(state_c) + if (state_f%mpack /= 0) then + call c_f_pointer(state_c%pa, rptr, [state_c%mpack]) + rptr = state_f%pa + if (state_f%uhf) then + call c_f_pointer(state_c%pb, rptr, [state_c%mpack]) + rptr = state_f%pb + end if + end if + end subroutine mopac_state_f2c + + subroutine mopac_state_c2f(state_c, state_f) + type(mopac_state), intent(in) :: state_c + type(mopac_state_f), intent(out) :: state_f + state_f%mpack = state_c%mpack + state_f%uhf = (state_c%uhf /= 0) + if (allocated(state_f%pa)) deallocate(state_f%pa) + if (allocated(state_f%pb)) deallocate(state_f%pb) + if (state_f%mpack /= 0) then + call copy_real(state_f%pa, state_c%pa, [state_c%mpack]) + if (state_f%uhf) call copy_real(state_f%pb, state_c%pb, [state_c%mpack]) + end if + call destroy_mopac_state(state_c) + end subroutine mopac_state_c2f + + subroutine mozyme_state_f2c(state_f, state_c) + type(mozyme_state_f), intent(in) :: state_f + type(mozyme_state), intent(out) :: state_c + integer(c_int), pointer :: iptr(:), iptr2(:,:) + real(c_double), pointer :: rptr(:) + state_c%numat = state_f%numat + state_c%noccupied = state_f%noccupied + state_c%nvirtual = state_f%nvirtual + state_c%icocc_dim = state_f%icocc_dim + state_c%icvir_dim = state_f%icvir_dim + state_c%cocc_dim = state_f%cocc_dim + state_c%cvir_dim = state_f%cvir_dim + call create_mozyme_state(state_c) + if (state_f%numat /= 0) then + call c_f_pointer(state_c%nbonds, iptr, [state_c%numat]) + iptr = state_f%nbonds + call c_f_pointer(state_c%ibonds, iptr2, [9,state_c%numat]) + iptr2 = state_f%ibonds + call c_f_pointer(state_c%iorbs, iptr, [state_c%numat]) + iptr = state_f%iorbs + call c_f_pointer(state_c%ncf, iptr, [state_c%noccupied]) + iptr = state_f%ncf + call c_f_pointer(state_c%nce, iptr, [state_c%nvirtual]) + iptr = state_f%nce + call c_f_pointer(state_c%icocc, iptr, [state_c%icocc_dim]) + iptr = state_f%icocc + call c_f_pointer(state_c%icvir, iptr, [state_c%icvir_dim]) + iptr = state_f%icvir + call c_f_pointer(state_c%cocc, rptr, [state_c%cocc_dim]) + rptr = state_f%cocc + call c_f_pointer(state_c%cvir, rptr, [state_c%cvir_dim]) + rptr = state_f%cvir + end if + end subroutine mozyme_state_f2c + + subroutine mozyme_state_c2f(state_c, state_f) + type(mozyme_state), intent(in) :: state_c + type(mozyme_state_f), intent(out) :: state_f + state_f%numat = state_c%numat + state_f%noccupied = state_c%noccupied + state_f%nvirtual = state_c%nvirtual + state_f%icocc_dim = state_c%icocc_dim + state_f%icvir_dim = state_c%icvir_dim + state_f%cocc_dim = state_c%cocc_dim + state_f%cvir_dim = state_c%cvir_dim + if (allocated(state_f%nbonds)) deallocate(state_f%nbonds) + if (allocated(state_f%ibonds)) deallocate(state_f%ibonds) + if (allocated(state_f%iorbs)) deallocate(state_f%iorbs) + if (allocated(state_f%ncf)) deallocate(state_f%ncf) + if (allocated(state_f%nce)) deallocate(state_f%nce) + if (allocated(state_f%icocc)) deallocate(state_f%icocc) + if (allocated(state_f%icvir)) deallocate(state_f%icvir) + if (allocated(state_f%cocc)) deallocate(state_f%cocc) + if (allocated(state_f%cvir)) deallocate(state_f%cvir) + if (state_f%numat /= 0) then + call copy_int(state_f%nbonds, state_c%nbonds, [state_c%numat]) + call copy_int2(state_f%ibonds, state_c%ibonds, [9,state_c%numat]) + call copy_int(state_f%iorbs, state_c%iorbs, [state_c%numat]) + call copy_int(state_f%ncf, state_c%ncf, [state_c%noccupied]) + call copy_int(state_f%nce, state_c%nce, [state_c%nvirtual]) + call copy_int(state_f%icocc, state_c%icocc, [state_c%icocc_dim]) + call copy_int(state_f%icvir, state_c%icvir, [state_c%icvir_dim]) + call copy_real(state_f%cocc, state_c%cocc, [state_c%cocc_dim]) + call copy_real(state_f%cvir, state_c%cvir, [state_c%cvir_dim]) + end if + call destroy_mozyme_state(state_c) + end subroutine mozyme_state_c2f + + subroutine copy_int(data, ptr, size) + integer, allocatable, intent(out) :: data(:) + type(c_ptr), intent(in) :: ptr + integer, intent(in) :: size(1) + integer(c_int), pointer :: data_c(:) + integer :: status + call c_f_pointer(ptr, data_c, size) + allocate(data(size(1)), stat=status) + if (status /= 0) then + write(*,*) "ERROR: Failed to allocate memory in MOPAC API wrapper" + stop 1 + end if + data = data_c + end subroutine copy_int + + subroutine copy_int2(data, ptr, size) + integer, allocatable, intent(out) :: data(:,:) + type(c_ptr), intent(in) :: ptr + integer, intent(in) :: size(2) + integer(c_int), pointer :: data_c(:,:) + integer :: status + call c_f_pointer(ptr, data_c, size) + allocate(data(size(1),size(2)), stat=status) + if (status /= 0) then + write(*,*) "ERROR: Failed to allocate memory in MOPAC API wrapper" + stop 1 + end if + data = data_c + end subroutine copy_int2 + + subroutine copy_real(data, ptr, size) + double precision, allocatable, intent(out) :: data(:) + type(c_ptr), intent(in) :: ptr + integer, intent(in) :: size(1) + real(c_double), pointer :: data_c(:) + integer :: status + call c_f_pointer(ptr, data_c, size) + allocate(data(size(1)), stat=status) + if (status /= 0) then + write(*,*) "ERROR: Failed to allocate memory in MOPAC API wrapper" + stop 1 + end if + data = data_c + end subroutine copy_real + + subroutine copy_real2(data, ptr, size) + double precision, allocatable, intent(out) :: data(:,:) + type(c_ptr), intent(in) :: ptr + integer, intent(in) :: size(2) + real(c_double), pointer :: data_c(:,:) + integer :: status + call c_f_pointer(ptr, data_c, size) + allocate(data(size(1),size(2)), stat=status) + if (status /= 0) then + write(*,*) "ERROR: Failed to allocate memory in MOPAC API wrapper" + stop 1 + end if + data = data_c + end subroutine copy_real2 + +end submodule mopac_api_f_internal diff --git a/src/MOZYME/iter_for_MOZYME.F90 b/src/MOZYME/iter_for_MOZYME.F90 index 056dc1ac..68c7b26e 100644 --- a/src/MOZYME/iter_for_MOZYME.F90 +++ b/src/MOZYME/iter_for_MOZYME.F90 @@ -263,6 +263,7 @@ subroutine iter_for_MOZYME (ee) return end if icocc(:icocc_dim) = iwork(:icocc_dim) + icocc(icocc_dim+1:) = 0 icocc_dim = Nint(icocc_dim*1.6) deallocate (iwork) @@ -279,6 +280,7 @@ subroutine iter_for_MOZYME (ee) return end if cocc(:cocc_dim) = rwork(:cocc_dim) + cocc(cocc_dim+1:) = 0.0d0 cocc_dim = Nint(cocc_dim*1.6) deallocate (rwork) @@ -295,6 +297,7 @@ subroutine iter_for_MOZYME (ee) return end if icvir(:icvir_dim) = iwork(:icvir_dim) + icvir(icvir_dim+1:) = 0 icvir_dim = Nint(icvir_dim*1.6) deallocate (iwork) @@ -311,6 +314,7 @@ subroutine iter_for_MOZYME (ee) return end if cvir(:cvir_dim) = rwork(:cvir_dim) + cvir(cvir_dim+1:) = 0.0d0 cvir_dim = Nint(cvir_dim*1.6) deallocate (rwork) end if @@ -361,6 +365,7 @@ subroutine iter_for_MOZYME (ee) return end if icocc(:icocc_dim) = iwork(:icocc_dim) + icocc(icocc_dim+1:) = 0 icocc_dim = Nint(icocc_dim*1.6) deallocate (iwork) @@ -377,6 +382,7 @@ subroutine iter_for_MOZYME (ee) return end if cocc(:cocc_dim) = rwork(:cocc_dim) + cocc(cocc_dim+1:) = 0.0d0 cocc_dim = Nint(cocc_dim*1.6) deallocate (rwork) @@ -393,6 +399,7 @@ subroutine iter_for_MOZYME (ee) return end if icvir(:icvir_dim) = iwork(:icvir_dim) + icvir(icvir_dim+1:) = 0 icvir_dim = Nint(icvir_dim*1.6) deallocate (iwork) @@ -409,6 +416,7 @@ subroutine iter_for_MOZYME (ee) return end if cvir(:cvir_dim) = rwork(:cvir_dim) + cvir(cvir_dim+1:) = 0.0d0 cvir_dim = Nint(cvir_dim*1.6) deallocate (rwork) end if diff --git a/src/chanel_C.F90 b/src/chanel_C.F90 index 274990c5..c18c3827 100644 --- a/src/chanel_C.F90 +++ b/src/chanel_C.F90 @@ -15,7 +15,9 @@ ! along with this program. If not, see . module chanel_C +#ifdef WIN32 !dec$ attributes dllexport :: iend, end_fn +#endif integer :: & & iw0 = -1, & ! Abbrev. output, for GUI (By default, not used) & ifiles_1 = 1, & diff --git a/src/interface/CMakeLists.txt b/src/interface/CMakeLists.txt index 19007113..4f258a4c 100644 --- a/src/interface/CMakeLists.txt +++ b/src/interface/CMakeLists.txt @@ -18,7 +18,7 @@ # Create a list of source files (src_list) with the .F90 extension set(src_list mopac_api mopac_api_initialize mopac_api_finalize - mopac_api_saveload mopac_api_operations + mopac_api_saveload mopac_api_createdestroy mopac_api_operations ) #----------------------------------------------- # Add a list of source files to the target diff --git a/src/interface/mopac_api.F90 b/src/interface/mopac_api.F90 index 19316b4c..e222a16a 100644 --- a/src/interface/mopac_api.F90 +++ b/src/interface/mopac_api.F90 @@ -15,8 +15,8 @@ ! along with this program. If not, see . ! Diskless/stateless Application Programming Interface (API) to core MOPAC operations -! NOTE: include this file or compile it to generate a .mod file to use the API module mopac_api + use, intrinsic :: iso_c_binding implicit none private @@ -24,188 +24,305 @@ module mopac_api public :: mopac_system, mopac_properties, mopac_state, mozyme_state ! public subroutines of the MOPAC API public :: mopac_scf, mopac_relax, mopac_vibe, mozyme_scf, mozyme_relax, mozyme_vibe + ! public subroutines to create derived types of the API that it can destroy + public :: create_mopac_state, create_mozyme_state + ! public subroutines to destroy derived types created by the API + public :: destroy_mopac_properties, destroy_mopac_state, destroy_mozyme_state ! public subroutine of the simple, disk-based MOPAC API public :: run_mopac_from_input ! data that defines the atomistic system and MOPAC job options - type :: mopac_system + type, bind(c) :: mopac_system ! number of atoms - integer :: natom = 0 + integer(c_int) :: natom ! number of atoms that are allowed to move (first natom_move atoms in array) - integer :: natom_move = 0 + integer(c_int) :: natom_move ! net charge - integer :: charge = 0 + integer(c_int) :: charge ! number of spin excitations, floor[(number of alpha electrons)/2 - (number of beta electrons)/2] - integer :: spin = 0 - ! dielectric constant for COSMO implicit solvent, must be 1 (no solvent) for nlattice > 0 - double precision :: epsilon = 1.d0 + integer(c_int) :: spin ! semiempirical model: PM7 = 0, PM6-D3H4 = 1, PM6-ORG = 2, PM6 = 3, AM1 = 4, RM1 = 5 - integer :: model = 0 - ! atomic number of each atom [natom] - integer, dimension (:), allocatable :: atom - ! (x,y,z) coordinates of each atom (Angstroms) [3*natom] - double precision, dimension (:), allocatable :: coord + integer(c_int) :: model + ! dielectric constant for COSMO implicit solvent, must be 1 (no solvent) for nlattice > 0 + real(c_double) :: epsilon + ! atomic number of each atom + type(c_ptr) :: atom ! integer(c_int)[natom] + ! (x,y,z) coordinates of each atom (Angstroms) + type(c_ptr) :: coord ! real(c_double)[3*natom] ! number of lattice vectors / translation vectors / periodic dimensions - integer :: nlattice = 0 + integer(c_int) :: nlattice ! number of lattice vectors that are allowed to move (first nlattice_move vectors in array) - integer :: nlattice_move = 0 + integer(c_int) :: nlattice_move ! external hydrostatic pressure (Gigapascals) - double precision :: pressure = 0.d0 - ! (x,y,z) coordinates of each lattice vectors (Angstroms) [3*nlattice] - double precision, dimension (:), allocatable :: lattice + real(c_double) :: pressure + ! (x,y,z) coordinates of each lattice vectors (Angstroms) + type(c_ptr) :: lattice ! real(c_double)[3*nlattice] ! relative numerical tolerances (equivalent to GNORM and RELSCF keyword values) - double precision :: tolerance = 1.d0 + real(c_double) :: tolerance ! time limit for a MOPAC calculation (seconds) - integer :: max_time = 3600 + integer(c_int) :: max_time end type ! calculated ground-state properties of an atomistic system and MOPAC job info - type :: mopac_properties + type, bind(c) :: mopac_properties ! heat of formation (kcal/mol) - double precision :: heat + real(c_double) :: heat ! dipole moment vector (Debye) - double precision, dimension (3) :: dipole - ! atomic partial charges [natom] - double precision, dimension (:), allocatable :: charge - ! (x,y,z) coordinates of each moveable atom (Angstroms) [3*natom_move] - double precision, dimension (:), allocatable :: coord_update - ! (x,y,z) heat gradients for each moveable atom (kcal/mol/Angstrom) [3*natom_move] - double precision, dimension (:), allocatable :: coord_deriv - ! availability of vibrational properties - logical :: calc_vibe - ! vibrational frequencies of normal modes (1/cm) [3*natom_move] - double precision, dimension (:), allocatable :: freq - ! (x,y,z) displacement vectors of normal modes [3*natom_move,3*natom_move] - double precision, dimension (:,:), allocatable :: disp - ! bond-order matrix in compressed sparse column (CSC) matrix format + real(c_double), dimension (3) :: dipole + ! atomic partial charges + type(c_ptr) :: charge ! real(c_double)[natom] + ! (x,y,z) coordinates of each moveable atom (Angstroms) + type(c_ptr) :: coord_update ! real(c_double)[3*natom_move] + ! (x,y,z) heat gradients for each moveable atom (kcal/mol/Angstrom) + type(c_ptr) :: coord_deriv ! real(c_double)[3*natom_move] + ! vibrational frequencies of normal modes (1/cm) + type(c_ptr) :: freq ! real(c_double)[3*natom_move], NULL if unavailable + ! (x,y,z) displacement vectors of normal modes + type(c_ptr) :: disp ! real(c_double)[3*natom_move,3*natom_move], NULL if unavailable + ! bond-order matrix in compressed sparse column (CSC) matrix format (0-based indexing) ! with insignificant bond orders (<0.01) truncated ! diagonal matrix entries are atomic valencies - ! > first index of each atom in CSC bond-order matrix [natom+1] - integer, dimension (:), allocatable :: bond_index - ! > list of atoms bonded to each atom in CSC format [bond_index(natom+1)-1] - integer, dimension (:), allocatable :: bond_atom - ! > bond order of atoms bonded to each atom in CSC format [bond_index(natom+1)-1] - double precision, dimension (:), allocatable :: bond_order - ! (x,y,z) coordinates of each moveable lattice vectors (Angstroms) [3*nlattice_move] - double precision, dimension (:), allocatable :: lattice_update - ! (x,y,z) heat gradients for each moveable lattice vector (kcal/mol/Angstrom) [3*nlattice_move] - double precision, dimension (:), allocatable :: lattice_deriv - ! stress tensor (Gigapascals) in Voigt form (xx, yy, zz, yz, xz, xy), if available - double precision, dimension (6) :: stress + ! > first index of each atom in CSC bond-order matrix + type(c_ptr) :: bond_index ! integer(c_int)[natom+1] + ! > list of atoms bonded to each atom in CSC format + type(c_ptr) :: bond_atom ! integer(c_int)[bond_index(natom+1)] + ! > bond order of atoms bonded to each atom in CSC format + type(c_ptr) :: bond_order ! real(c_double)[bond_index(natom+1)] + ! (x,y,z) coordinates of each moveable lattice vectors (Angstroms) + type(c_ptr) :: lattice_update ! real(c_double)[3*nlattice_move] + ! (x,y,z) heat gradients for each moveable lattice vector (kcal/mol/Angstrom) + type(c_ptr) :: lattice_deriv ! real(c_double)[3*nlattice_move] + ! stress tensor (Gigapascals) in Voigt form (xx, yy, zz, yz, xz, xy) + real(c_double), dimension (6) :: stress ! 0 if unavailable ! number of MOPAC error messages (negative value indicates that allocation of error_msg failed) - integer :: nerror - ! text of MOPAC error messages [nerror,120] - character*120, dimension (:), allocatable :: error_msg + integer(c_int) :: nerror + ! text of MOPAC error messages + type(c_ptr) :: error_msg ! type(c_ptr)[nerror] -> character(kind=c_char)[*] end type ! data that describes the electronic state using standard molecular orbitals - type :: mopac_state - ! availability of a saved state - logical :: save_state = .false. - ! MOPAC data format is copied from molkst_C and Common_arrays_C modules + type, bind(c) :: mopac_state + ! MOPAC data format is adapted from molkst_C and Common_arrays_C modules ! > number of matrix elements in packed lower triangle matrix format - integer :: mpack - ! > alpha density matrix [mpack] - double precision, dimension (:), allocatable :: pa - ! > beta density matrix [mpack] - double precision, dimension (:), allocatable :: pb + integer(c_int) :: mpack ! 0 if state is unavailable + ! > flag for unrestricted Hartree-Fock ground state (0 == restricted, 1 == unrestricted) + integer(c_int) :: uhf + ! > alpha density matrix + type(c_ptr) :: pa ! real(c_double)[mpack] + ! > beta density matrix + type(c_ptr) :: pb ! real(c_double)[mpack], NULL if uhf == 0 end type ! data that describes the electronic state using localized molecular orbitals - type :: mozyme_state - ! availability of a saved state - logical :: save_state = .false. - ! MOZYME data format is copied from molkst_C, Common_arrays_C, and MOZYME_C modules + type, bind(c) :: mozyme_state + ! MOZYME data format is adapted from molkst_C, Common_arrays_C, and MOZYME_C modules ! > number of real atoms - integer :: numat - ! > number of Lewis bonds per real atom [numat] - integer, dimension (:), allocatable :: nbonds - ! > list of Lewis-bonded real atoms for each real atom [9,numat] - integer, dimension (:,:), allocatable :: ibonds - ! > number of orbitals per real atom [numat] - integer, dimension (:), allocatable :: iorbs + integer(c_int) :: numat ! 0 if state is unavailable + ! > number of Lewis bonds per real atom + type(c_ptr) :: nbonds ! integer(c_int)[numat] + ! > list of Lewis-bonded real atoms for each real atom + type(c_ptr) :: ibonds ! integer(c_int)[9,numat] + ! > number of orbitals per real atom + type(c_ptr) :: iorbs ! integer(c_int)[numat] ! > number of occupied molecular orbitals - integer :: noccupied - ! > number of atoms in each occupied LMO [noccupied] - integer, dimension (:), allocatable :: ncf + integer(c_int) :: noccupied + ! > number of atoms in each occupied LMO + type(c_ptr) :: ncf ! integer(c_int)[noccupied] ! > number of virtual molecular orbitals - integer :: nvirtual - ! > number of atoms in each virtual LMO [nvirtual] - integer, dimension (:), allocatable :: nce + integer(c_int) :: nvirtual + ! > number of atoms in each virtual LMO + type(c_ptr) :: nce ! integer(c_int)[nvirtual] ! > size of array icocc - integer :: icocc_dim - ! > index of each real atom in the occupied LMOs [iccoc_dim] - integer, dimension (:), allocatable :: icocc + integer(c_int) :: icocc_dim + ! > index of each real atom in the occupied LMOs + type(c_ptr) :: icocc ! integer(c_int)[iccoc_dim] ! > size of array icvir - integer :: icvir_dim - ! > index of each real atom in the virtual LMOs [icvir_dim] - integer, dimension (:), allocatable :: icvir + integer(c_int) :: icvir_dim + ! > index of each real atom in the virtual LMOs + type(c_ptr) :: icvir ! integer(c_int)[icvir_dim] ! > size of array cocc - integer :: cocc_dim - ! > atomic orbital coefficients of the occupied LMOs [cocc_dim] - double precision, dimension (:), allocatable :: cocc + integer(c_int) :: cocc_dim + ! > atomic orbital coefficients of the occupied LMOs + type(c_ptr) :: cocc ! real(c_double)[cocc_dim] ! > size of array cvir - integer :: cvir_dim - ! > atomic orbital coefficients of the virtual LMOs [cvir_dim] - double precision, dimension (:), allocatable :: cvir + integer(c_int) :: cvir_dim + ! > atomic orbital coefficients of the virtual LMOs + type(c_ptr) :: cvir ! real(c_double)[cvir_dim] end type interface ! MOPAC electronic ground state calculation - module subroutine mopac_scf(system, state, properties) - !dec$ attributes dllexport :: mopac_scf + module subroutine mopac_scf(system, state, properties) bind(c) +#ifdef WIN32 +!dec$ attributes dllexport :: mopac_scf +#endif type(mopac_system), intent(in) :: system type(mopac_state), intent(inout) :: state type(mopac_properties), intent(out) :: properties end subroutine mopac_scf ! MOPAC geometry relaxation - module subroutine mopac_relax(system, state, properties) - !dec$ attributes dllexport :: mopac_relax + module subroutine mopac_relax(system, state, properties) bind(c) +#ifdef WIN32 +!dec$ attributes dllexport :: mopac_relax +#endif type(mopac_system), intent(in) :: system type(mopac_state), intent(inout) :: state type(mopac_properties), intent(out) :: properties end subroutine mopac_relax ! MOPAC vibrational calculation - module subroutine mopac_vibe(system, state, properties) - !dec$ attributes dllexport :: mopac_vibe + module subroutine mopac_vibe(system, state, properties) bind(c) +#ifdef WIN32 +!dec$ attributes dllexport :: mopac_vibe +#endif type(mopac_system), intent(in) :: system type(mopac_state), intent(inout) :: state type(mopac_properties), intent(out) :: properties end subroutine mopac_vibe ! MOZYME electronic ground state calculation - module subroutine mozyme_scf(system, state, properties) - !dec$ attributes dllexport :: mozyme_scf + module subroutine mozyme_scf(system, state, properties) bind(c) +#ifdef WIN32 +!dec$ attributes dllexport :: mozyme_scf +#endif type(mopac_system), intent(in) :: system type(mozyme_state), intent(inout) :: state type(mopac_properties), intent(out) :: properties end subroutine mozyme_scf ! MOZYME geometry relaxation - module subroutine mozyme_relax(system, state, properties) - !dec$ attributes dllexport :: mozyme_relax + module subroutine mozyme_relax(system, state, properties) bind(c) +#ifdef WIN32 +!dec$ attributes dllexport :: mozyme_relax +#endif type(mopac_system), intent(in) :: system type(mozyme_state), intent(inout) :: state type(mopac_properties), intent(out) :: properties end subroutine mozyme_relax ! MOZYME vibrational calculation - module subroutine mozyme_vibe(system, state, properties) - !dec$ attributes dllexport :: mozyme_vibe + module subroutine mozyme_vibe(system, state, properties) bind(c) +#ifdef WIN32 +!dec$ attributes dllexport :: mozyme_vibe +#endif type(mopac_system), intent(in) :: system type(mozyme_state), intent(inout) :: state type(mopac_properties), intent(out) :: properties end subroutine mozyme_vibe - ! Run MOPAC conventionally from an input file - module subroutine run_mopac_from_input(path_to_file) - !dec$ attributes dllexport :: run_mopac_from_input - character(len=240), intent(in) :: path_to_file - end subroutine run_mopac_from_input - + ! allocate memory for mopac_state + module subroutine create_mopac_state(state) bind(c) +#ifdef WIN32 +!dec$ attributes dllexport :: create_mopac_state +#endif + type(mopac_state), intent(inout) :: state + end subroutine create_mopac_state + + ! allocate memory for mozyme_state + module subroutine create_mozyme_state(state) bind(c) +#ifdef WIN32 +!dec$ attributes dllexport :: create_mozyme_state +#endif + type(mozyme_state), intent(inout) :: state + end subroutine create_mozyme_state + + ! deallocate memory in mopac_properties + module subroutine destroy_mopac_properties(properties) bind(c) +#ifdef WIN32 +!dec$ attributes dllexport :: destroy_mopac_properties +#endif + type(mopac_properties), intent(in) :: properties + end subroutine destroy_mopac_properties + + ! deallocate memory in mopac_state + module subroutine destroy_mopac_state(state) bind(c) +#ifdef WIN32 +!dec$ attributes dllexport :: destroy_mopac_state +#endif + type(mopac_state), intent(in) :: state + end subroutine destroy_mopac_state + + ! deallocate memory in mozyme_state + module subroutine destroy_mozyme_state(state) bind(c) +#ifdef WIN32 +!dec$ attributes dllexport :: destroy_mozyme_state +#endif + type(mozyme_state), intent(in) :: state + end subroutine destroy_mozyme_state + + ! run MOPAC conventionally from an input file + module function run_mopac_from_input(path_to_file) bind(c) +#ifdef WIN32 +!dec$ attributes dllexport :: run_mopac_from_input +#endif + integer(c_int) :: run_mopac_from_input + character(kind=c_char), dimension(*), intent(in) :: path_to_file + end function run_mopac_from_input + + end interface + + ! allocate memory (C or Fortran memory manager, depending on compiler) + interface + module function create_int(size) + type(c_ptr) :: create_int + integer, intent(in) :: size + end function create_int + module function create_int2(size, size2) + type(c_ptr) :: create_int2 + integer, intent(in) :: size + integer, intent(in) :: size2 + end function create_int2 + module function create_real(size) + type(c_ptr) :: create_real + integer, intent(in) :: size + end function create_real + end interface + + ! allocate memory & copy data (C or Fortran memory manager, depending on compiler) + interface create_copy + module function create_copy_int(array, size) + type(c_ptr) :: create_copy_int + integer, intent(in) :: array(:) + integer, intent(in) :: size(1) + end function create_copy_int + module function create_copy_int2(array, size) + type(c_ptr) :: create_copy_int2 + integer, intent(in) :: array(:,:) + integer, intent(in) :: size(2) + end function create_copy_int2 + module function create_copy_real(array, size) + type(c_ptr) :: create_copy_real + double precision, intent(in) :: array(:) + integer, intent(in) :: size(1) + end function create_copy_real + module function create_copy_char(array, size) + type(c_ptr) :: create_copy_char + character(len=*), intent(in) :: array + integer, intent(in) :: size(1) + end function create_copy_char + module function create_copy_ptr(array, size) + type(c_ptr) :: create_copy_ptr + type(c_ptr), intent(in) :: array(:) + integer, intent(in) :: size(1) + end function create_copy_ptr + end interface + + ! deallocate memory (C or Fortran memory manager, depending on compiler) + interface + module subroutine destroy_int(copy) + type(c_ptr), intent(in) :: copy + end subroutine destroy_int + module subroutine destroy_real(copy) + type(c_ptr), intent(in) :: copy + end subroutine destroy_real + module subroutine destroy_char(copy) + type(c_ptr), intent(in) :: copy + end subroutine destroy_char + module subroutine destroy_ptr(copy) + type(c_ptr), intent(in) :: copy + end subroutine destroy_ptr end interface end module mopac_api diff --git a/src/interface/mopac_api_createdestroy.F90 b/src/interface/mopac_api_createdestroy.F90 new file mode 100644 index 00000000..c4863b4f --- /dev/null +++ b/src/interface/mopac_api_createdestroy.F90 @@ -0,0 +1,373 @@ +! Molecular Orbital PACkage (MOPAC) +! Copyright (C) 2021, Virginia Polytechnic Institute and State University +! +! MOPAC is free software: you can redistribute it and/or modify it under +! the terms of the GNU Lesser General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! MOPAC is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU Lesser General Public License for more details. +! +! You should have received a copy of the GNU Lesser General Public License +! along with this program. If not, see . + +! Diskless/stateless Application Programming Interface (API) to core MOPAC operations +submodule (mopac_api) mopac_api_createdestroy + implicit none + +contains + + ! allocate memory for mopac_state + module subroutine create_mopac_state(state) bind(c) +#ifdef WIN32 +!dec$ attributes dllexport :: create_mopac_state +#endif + type(mopac_state), intent(inout) :: state + if (state%mpack /= 0) then + state%pa = create_real(state%mpack) + if (state%uhf == 1) state%pb = create_real(state%mpack) + end if + end subroutine create_mopac_state + + ! allocate memory for mozyme_state + module subroutine create_mozyme_state(state) bind(c) +#ifdef WIN32 +!dec$ attributes dllexport :: create_mozyme_state +#endif + type(mozyme_state), intent(inout) :: state + if (state%numat /= 0) then + state%nbonds = create_int(state%numat) + state%ibonds = create_int2(9,state%numat) + state%iorbs = create_int(state%numat) + state%ncf = create_int(state%noccupied) + state%nce = create_int(state%nvirtual) + state%icocc = create_int(state%icocc_dim) + state%icvir = create_int(state%icvir_dim) + state%cocc = create_real(state%cocc_dim) + state%cvir = create_real(state%cvir_dim) + end if + end subroutine create_mozyme_state + + ! deallocate memory in mopac_properties + module subroutine destroy_mopac_properties(properties) bind(c) +#ifdef WIN32 +!dec$ attributes dllexport :: destroy_mopac_properties +#endif + type(mopac_properties), intent(in) :: properties + integer :: i + type(c_ptr), pointer :: pptr(:) + call destroy_real(properties%charge) + call destroy_int(properties%bond_index) + call destroy_int(properties%bond_atom) + call destroy_real(properties%bond_order) + call destroy_real(properties%coord_update) + call destroy_real(properties%coord_deriv) + call destroy_real(properties%freq) + call destroy_real(properties%disp) + call destroy_real(properties%lattice_update) + call destroy_real(properties%lattice_deriv) + if (properties%nerror > 0) then + call c_f_pointer(properties%error_msg, pptr, [properties%nerror]) + do i=1, properties%nerror + call destroy_char(pptr(i)) + end do + call destroy_ptr(properties%error_msg) + end if + end subroutine destroy_mopac_properties + + ! deallocate memory in mopac_state + module subroutine destroy_mopac_state(state) bind(c) +#ifdef WIN32 +!dec$ attributes dllexport :: destroy_mopac_state +#endif + type(mopac_state), intent(in) :: state + if (state%mpack /= 0) then + call destroy_real(state%pa) + if (state%uhf == 1) call destroy_real(state%pb) + end if + end subroutine destroy_mopac_state + + ! deallocate memory in mozyme_state + module subroutine destroy_mozyme_state(state) bind(c) +#ifdef WIN32 +!dec$ attributes dllexport :: destroy_mozyme_state +#endif + type(mozyme_state), intent(in) :: state + if (state%numat /= 0) then + call destroy_int(state%nbonds) + call destroy_int(state%ibonds) + call destroy_int(state%iorbs) + call destroy_int(state%ncf) + call destroy_int(state%nce) + call destroy_int(state%icocc) + call destroy_int(state%icvir) + call destroy_real(state%cocc) + call destroy_real(state%cvir) + end if + end subroutine destroy_mozyme_state + + ! Unfortunately, the Intel Fortran compiler does not adhere to the Fortran standard for pointers. + ! Memory allocated to a Fortran pointer erroneously cannot be deallocated if its memory is passed + ! through a C pointer in a C-bound interface and then reassigned to a Fortran pointer because hidden + ! information about memory allocation is contained within the original Fortran pointer and is not + ! retained by the C pointer. To get around this, MOPAC's API will support both C and Fortran memory + ! managers through the presence/absence of the MOPAC_API_MALLOC preprocessor variable. + + ! allocate memory (C or Fortran memory manager, depending on compiler) + module function create_int(size) + integer, intent(in) :: size + type(c_ptr) :: create_int +#ifndef MOPAC_API_MALLOC + integer(c_int), pointer :: ptr(:) + integer :: status + allocate(ptr(size), stat=status) + if (status /= 0) then + write(*,*) "ERROR: Failed to allocate memory in MOPAC API I/O" + stop 1 + end if + create_int = c_loc(ptr) +#else + integer(c_intptr_t) :: dummy + integer(c_int) :: mold + dummy = malloc(c_sizeof(mold)*size) + create_int = transfer(dummy, create_int) +#endif + end function create_int + module function create_int2(size, size2) + integer, intent(in) :: size + integer, intent(in) :: size2 + type(c_ptr) :: create_int2 +#ifndef MOPAC_API_MALLOC + integer(c_int), pointer :: ptr(:,:) + integer :: status + allocate(ptr(size,size2), stat=status) + if (status /= 0) then + write(*,*) "ERROR: Failed to allocate memory in MOPAC API I/O" + stop 1 + end if + create_int2 = c_loc(ptr) +#else + integer(c_intptr_t) :: dummy + integer(c_int) :: mold + dummy = malloc(c_sizeof(mold)*size*size2) + create_int2 = transfer(dummy, create_int2) +#endif + end function create_int2 + module function create_real(size) + integer, intent(in) :: size + type(c_ptr) :: create_real +#ifndef MOPAC_API_MALLOC + real(c_double), pointer :: ptr(:) + integer :: status + allocate(ptr(size), stat=status) + if (status /= 0) then + write(*,*) "ERROR: Failed to allocate memory in MOPAC API I/O" + stop 1 + end if + create_real = c_loc(ptr) +#else + integer(c_intptr_t) :: dummy + real(c_double) :: mold + dummy = malloc(c_sizeof(mold)*size) + create_real = transfer(dummy, create_real) +#endif + end function create_real + + ! allocate memory & copy data (C or Fortran memory manager, depending on compiler) + module function create_copy_int(array, size) + integer, intent(in) :: array(:) + integer, intent(in) :: size(1) + type(c_ptr) :: create_copy_int + integer(c_int), pointer :: ptr(:) +#ifndef MOPAC_API_MALLOC + integer :: status + allocate(ptr(size(1)), stat=status) + if (status /= 0) then + write(*,*) "ERROR: Failed to allocate memory in MOPAC API I/O" + stop 1 + end if + create_copy_int = c_loc(ptr) +#else + integer(c_intptr_t) :: dummy + integer(c_int) :: mold + dummy = malloc(c_sizeof(mold)*size(1)) + create_copy_int = transfer(dummy, create_copy_int) + call c_f_pointer(create_copy_int, ptr, size) +#endif + ptr = array(:size(1)) + end function create_copy_int + module function create_copy_int2(array, size) + integer, intent(in) :: array(:,:) + integer, intent(in) :: size(2) + type(c_ptr) :: create_copy_int2 + integer(c_int), pointer :: ptr(:,:) +#ifndef MOPAC_API_MALLOC + integer :: status + allocate(ptr(size(1),size(2)), stat=status) + if (status /= 0) then + write(*,*) "ERROR: Failed to allocate memory in MOPAC API I/O" + stop 1 + end if + create_copy_int2 = c_loc(ptr) +#else + integer(c_intptr_t) :: dummy + integer(c_int) :: mold + dummy = malloc(c_sizeof(mold)*size(1)*size(2)) + create_copy_int2 = transfer(dummy, create_copy_int2) + call c_f_pointer(create_copy_int2, ptr, size) +#endif + ptr = array(:size(1),:size(2)) + end function create_copy_int2 + module function create_copy_real(array, size) + double precision, intent(in) :: array(:) + integer, intent(in) :: size(1) + type(c_ptr) :: create_copy_real + real(c_double), pointer :: ptr(:) +#ifndef MOPAC_API_MALLOC + integer :: status + allocate(ptr(size(1)), stat=status) + if (status /= 0) then + write(*,*) "ERROR: Failed to allocate memory in MOPAC API I/O" + stop 1 + end if + create_copy_real = c_loc(ptr) +#else + integer(c_intptr_t) :: dummy + real(c_double) :: mold + dummy = malloc(c_sizeof(mold)*size(1)) + create_copy_real = transfer(dummy, create_copy_real) + call c_f_pointer(create_copy_real, ptr, size) +#endif + ptr = array(:size(1)) + end function create_copy_real + module function create_copy_char(array, size) + character(len=*), intent(in) :: array + integer, intent(in) :: size(1) + type(c_ptr) :: create_copy_char + character(kind=c_char), pointer :: ptr(:) + integer :: i +#ifndef MOPAC_API_MALLOC + integer :: status + allocate(ptr(size(1)), stat=status) + if (status /= 0) then + write(*,*) "ERROR: Failed to allocate memory in MOPAC API I/O" + stop 1 + end if + create_copy_char = c_loc(ptr) +#else + integer(c_intptr_t) :: dummy + character(kind=c_char) :: mold + dummy = malloc(c_sizeof(mold)*size(1)) + create_copy_char = transfer(dummy, create_copy_char) + call c_f_pointer(create_copy_char, ptr, size) +#endif + do i=1, size(1)-1 + ptr(i) = array(i:i) + end do + ptr(size(1)) = c_null_char + end function create_copy_char + module function create_copy_ptr(array, size) + type(c_ptr), intent(in) :: array(:) + integer, intent(in) :: size(1) + type(c_ptr) :: create_copy_ptr + type(c_ptr), pointer :: ptr(:) +#ifndef MOPAC_API_MALLOC + integer :: status + allocate(ptr(size(1)), stat=status) + if (status /= 0) then + write(*,*) "ERROR: Failed to allocate memory in MOPAC API I/O" + stop 1 + end if + create_copy_ptr = c_loc(ptr) +#else + integer(c_intptr_t) :: dummy + type(c_ptr) :: mold + dummy = malloc(c_sizeof(mold)*size(1)) + create_copy_ptr = transfer(dummy, create_copy_ptr) + call c_f_pointer(create_copy_ptr, ptr, size) +#endif + ptr = array(:size(1)) + end function create_copy_ptr + + ! deallocate memory (C or Fortran memory manager, depending on compiler) + module subroutine destroy_int(copy) + type(c_ptr), intent(in) :: copy +#ifndef MOPAC_API_MALLOC + integer(c_int), pointer :: ptr + integer :: status + if (c_associated(copy)) then + call c_f_pointer(copy, ptr) + deallocate(ptr, stat=status) + if (status /= 0) then + write(*,*) "ERROR: Failed to deallocate memory in MOPAC API I/O" + stop 1 + end if + end if +#else + integer(c_intptr_t) :: copy2 + copy2 = transfer(copy, copy2) + call free(copy2) +#endif + end subroutine destroy_int + module subroutine destroy_real(copy) + type(c_ptr), intent(in) :: copy +#ifndef MOPAC_API_MALLOC + real(c_double), pointer :: ptr + integer :: status + if (c_associated(copy)) then + call c_f_pointer(copy, ptr) + deallocate(ptr, stat=status) + if (status /= 0) then + write(*,*) "ERROR: Failed to deallocate memory in MOPAC API I/O" + stop 1 + end if + end if +#else + integer(c_intptr_t) :: copy2 + copy2 = transfer(copy, copy2) + call free(copy2) +#endif + end subroutine destroy_real + module subroutine destroy_char(copy) + type(c_ptr), intent(in) :: copy +#ifndef MOPAC_API_MALLOC + character(kind=c_char), pointer :: ptr + integer :: status + if (c_associated(copy)) then + call c_f_pointer(copy, ptr) + deallocate(ptr, stat=status) + if (status /= 0) then + write(*,*) "ERROR: Failed to deallocate memory in MOPAC API I/O" + stop 1 + end if + end if +#else + integer(c_intptr_t) :: copy2 + copy2 = transfer(copy, copy2) + call free(copy2) +#endif + end subroutine destroy_char + module subroutine destroy_ptr(copy) + type(c_ptr), intent(in) :: copy +#ifndef MOPAC_API_MALLOC + type(c_ptr), pointer :: ptr + integer :: status + if (c_associated(copy)) then + call c_f_pointer(copy, ptr) + deallocate(ptr, stat=status) + if (status /= 0) then + write(*,*) "ERROR: Failed to deallocate memory in MOPAC API I/O" + stop 1 + end if + end if +#else + integer(c_intptr_t) :: copy2 + copy2 = transfer(copy, copy2) + call free(copy2) +#endif + end subroutine destroy_ptr + +end submodule mopac_api_createdestroy diff --git a/src/interface/mopac_api_finalize.F90 b/src/interface/mopac_api_finalize.F90 index 6e88aebf..61e5a538 100644 --- a/src/interface/mopac_api_finalize.F90 +++ b/src/interface/mopac_api_finalize.F90 @@ -47,39 +47,42 @@ ! save properties and clean up after a MOPAC/MOZYME calculation module subroutine mopac_finalize(properties) type(mopac_properties), intent(out) :: properties - integer :: status, i + integer :: status, i, j, size + type(c_ptr), allocatable :: pptr(:) ! close dummy output file to free up /dev/null close(iw) - ! deallocate any prior arrays - if (allocated(properties%coord_update)) deallocate(properties%coord_update) - if (allocated(properties%coord_deriv)) deallocate(properties%coord_deriv) - if (allocated(properties%freq)) deallocate(properties%freq) - if (allocated(properties%disp)) deallocate(properties%disp) - if (allocated(properties%charge)) deallocate(properties%charge) - if (allocated(properties%bond_index)) deallocate(properties%bond_index) - if (allocated(properties%bond_atom)) deallocate(properties%bond_atom) - if (allocated(properties%bond_order)) deallocate(properties%bond_order) - if (allocated(properties%lattice_update)) deallocate(properties%lattice_update) - if (allocated(properties%lattice_deriv)) deallocate(properties%lattice_deriv) - ! record properties if (.not. moperr) call mopac_record(properties) - ! collect error messages + ! collect error messages & assign NULL pointers for memory safety if (moperr) then + properties%charge = c_null_ptr + properties%coord_update = c_null_ptr + properties%coord_deriv = c_null_ptr + properties%lattice_update = c_null_ptr + properties%lattice_deriv = c_null_ptr + properties%freq = c_null_ptr + properties%disp = c_null_ptr + properties%bond_index = c_null_ptr + properties%bond_atom = c_null_ptr + properties%bond_order = c_null_ptr call summary("",0) properties%nerror = dummy - allocate(properties%error_msg(properties%nerror), stat=status) + allocate(pptr(properties%nerror), stat=status) if (status /= 0) then properties%nerror = -properties%nerror + properties%error_msg = c_null_ptr + return else do i=1, properties%nerror call summary("",-i) - properties%error_msg(i) = trim(errtxt) + size = len_trim(errtxt) + 1 + pptr(i) = create_copy(errtxt, [size]) end do end if + properties%error_msg = create_copy(pptr, [properties%nerror]) call summary("",-abs(properties%nerror)-1) else properties%nerror = 0 @@ -98,6 +101,8 @@ subroutine mopac_record(properties) double precision, external :: dipole, dipole_for_MOZYME integer :: status, i, j, k, kk, kl, ku, io, jo, natom_move, nlattice_move double precision :: valenc, sum, dumy(3) + integer(c_int), pointer :: bond_index(:), bond_atom(:) + real(c_double), pointer :: bond_order(:) ! trigger charge & dipole calculation call chrge (p, q) @@ -115,12 +120,7 @@ subroutine mopac_record(properties) end if ! save basic properties properties%heat = escf - allocate(properties%charge(numat), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOPAC_FINALIZE") - return - end if - properties%charge = q(:numat) + properties%charge = create_copy(q, [numat]) properties%stress = voigt natom_move = nvar/3 nlattice_move = 0 @@ -131,62 +131,32 @@ subroutine mopac_record(properties) end if end do ! save properties of moveable coordinates - allocate(properties%coord_update(3*natom_move), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOPAC_FINALIZE") - return - end if - properties%coord_update = xparam(:3*natom_move) - allocate(properties%coord_deriv(3*natom_move), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOPAC_FINALIZE") - return - end if - properties%coord_deriv = grad(:3*natom_move) + properties%coord_update = create_copy(xparam, [3*natom_move]) + properties%coord_deriv = create_copy(grad, [3*natom_move]) if (nlattice_move > 0) then - allocate(properties%lattice_update(3*nlattice_move), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOPAC_FINALIZE") - return - end if - properties%lattice_update = xparam(3*natom_move+1:) - allocate(properties%lattice_deriv(3*nlattice_move), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOPAC_FINALIZE") - return - end if - properties%lattice_deriv = grad(3*natom_move+1:) + properties%lattice_update = create_copy(xparam(3*natom_move+1:), [3*nlattice_move]) + properties%lattice_deriv = create_copy(grad(3*natom_move+1:), [3*nlattice_move]) + else + properties%lattice_update = c_null_ptr + properties%lattice_deriv = c_null_ptr end if ! save vibrational properties if available if (index(keywrd, " FORCETS") /= 0) then - properties%calc_vibe = .true. - allocate(properties%freq(nvar), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOPAC_FINALIZE") - return - end if - allocate(properties%disp(nvar,nvar), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOPAC_FINALIZE") - return - end if - properties%freq = freq - properties%disp = reshape(cnorml,[nvar, nvar]) + properties%freq = create_copy(freq, [nvar]) + properties%disp = create_copy(cnorml, [nvar*nvar]) else - properties%calc_vibe = .false. + properties%freq = c_null_ptr + properties%disp = c_null_ptr end if ! prune bond order matrix - allocate(properties%bond_index(numat+1), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOPAC_FINALIZE") - return - end if + properties%bond_index = create_int(numat+1) + call c_f_pointer(properties%bond_index, bond_index, [numat+1]) if (mozyme) then ! 1st pass to populate bond_index - properties%bond_index(1) = 1 + bond_index(1) = 0 do i = 1, numat io = iorbs(i) - properties%bond_index(i+1) = properties%bond_index(i) + bond_index(i+1) = bond_index(i) valenc = 0.d0 if (io == 1) then kk = ijbo (i, i) + 1 @@ -203,7 +173,7 @@ subroutine mopac_record(properties) end do end if if (valenc > 0.01d0) then - properties%bond_index(i+1) = properties%bond_index(i+1) + 1 + bond_index(i+1) = bond_index(i+1) + 1 end if do j = 1, numat jo = iorbs(j) @@ -215,22 +185,16 @@ subroutine mopac_record(properties) sum = sum + p(k) ** 2 end do if (sum > 0.01d0) then - properties%bond_index(i+1) = properties%bond_index(i+1) + 1 + bond_index(i+1) = bond_index(i+1) + 1 end if end if end do end do ! 2nd pass to populate bond_atom and bond_order - allocate(properties%bond_atom(properties%bond_index(numat+1)), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOPAC_FINALIZE") - return - end if - allocate(properties%bond_order(properties%bond_index(numat+1)), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOPAC_FINALIZE") - return - end if + properties%bond_atom = create_int(bond_index(numat+1)) + call c_f_pointer(properties%bond_atom, bond_atom, [bond_index(numat+1)]) + properties%bond_order = create_real(bond_index(numat+1)) + call c_f_pointer(properties%bond_order, bond_order, [bond_index(numat+1)]) do i = 1, numat io = iorbs(i) valenc = 0.d0 @@ -248,7 +212,7 @@ subroutine mopac_record(properties) valenc = valenc + 2.d0 * p(kk) end do end if - kk = properties%bond_index(i) + kk = bond_index(i) + 1 do j = 1, numat jo = iorbs(j) if (i /= j .and. ijbo (i, j) >= 0) then @@ -259,13 +223,13 @@ subroutine mopac_record(properties) sum = sum + p(k) ** 2 end do if (sum > 0.01d0) then - properties%bond_atom(kk) = j - properties%bond_order(kk) = sum + bond_atom(kk) = j + bond_order(kk) = sum kk = kk + 1 end if else if (valenc > 0.01d0) then - properties%bond_atom(kk) = j - properties%bond_order(kk) = valenc + bond_atom(kk) = j + bond_order(kk) = valenc kk = kk + 1 end if end do @@ -273,51 +237,45 @@ subroutine mopac_record(properties) else call bonds() ! 1st pass to populate bond_index - properties%bond_index(1) = 1 + bond_index(1) = 0 do i = 1, numat - properties%bond_index(i+1) = properties%bond_index(i) + bond_index(i+1) = bond_index(i) ku = i*(i-1)/2 + 1 kl = (i+1)*(i+2)/2 - 1 do j = 1, i if (bondab(ku) > 0.01d0) then - properties%bond_index(i+1) = properties%bond_index(i+1) + 1 + bond_index(i+1) = bond_index(i+1) + 1 end if ku = ku + 1 end do do j = i+1, numat if (bondab(kl) > 0.01d0) then - properties%bond_index(i+1) = properties%bond_index(i+1) + 1 + bond_index(i+1) = bond_index(i+1) + 1 end if kl = kl + j end do end do ! 2nd pass to populate bond_atom and bond_order - allocate(properties%bond_atom(properties%bond_index(numat+1)), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOPAC_FINALIZE") - return - end if - allocate(properties%bond_order(properties%bond_index(numat+1)), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOPAC_FINALIZE") - return - end if + properties%bond_atom = create_int(bond_index(numat+1)) + call c_f_pointer(properties%bond_atom, bond_atom, [bond_index(numat+1)]) + properties%bond_order = create_real(bond_index(numat+1)) + call c_f_pointer(properties%bond_order, bond_order, [bond_index(numat+1)]) do i = 1, numat ku = i*(i-1)/2 + 1 kl = (i+1)*(i+2)/2 - 1 - kk = properties%bond_index(i) + kk = bond_index(i) + 1 do j = 1, i if (bondab(ku) > 0.01d0) then - properties%bond_atom(kk) = j - properties%bond_order(kk) = bondab(ku) + bond_atom(kk) = j + bond_order(kk) = bondab(ku) kk = kk + 1 end if ku = ku + 1 end do do j = i+1, numat if (bondab(kl) > 0.01d0) then - properties%bond_atom(kk) = j - properties%bond_order(kk) = bondab(kl) + bond_atom(kk) = j + bond_order(kk) = bondab(kl) kk = kk + 1 end if kl = kl + j diff --git a/src/interface/mopac_api_initialize.F90 b/src/interface/mopac_api_initialize.F90 index 34e9f5bf..1b44bf69 100644 --- a/src/interface/mopac_api_initialize.F90 +++ b/src/interface/mopac_api_initialize.F90 @@ -78,6 +78,9 @@ module subroutine mopac_initialize(system) character(100) :: num2str integer :: i, j, nelectron, status double precision :: eat + integer(c_int), pointer :: atom(:) + real(c_double), pointer :: coord(:) + real(c_double), pointer :: lattice(:) ! validity checks of data in system if (system%natom <= 0) call mopend("Invalid number of atoms") @@ -85,20 +88,26 @@ module subroutine mopac_initialize(system) call mopend("Invalid number of moveable atoms") if (system%epsilon /= 1.d0 .and. system%nlattice > 0) & call mopend("COSMO solvent not available for periodic systems") - if (size(system%atom) < system%natom) call mopend("List of atomic numbers is too small") - if (size(system%coord) < 3*system%natom) call mopend("List of atomic coordinates is too small") + if (system%natom > 0 .and. .not. c_associated(system%atom)) & + call mopend("Array of atomic numbers not found") + if (system%natom > 0 .and. .not. c_associated(system%coord)) & + call mopend("Array of atomic coordinates is too small") if (system%nlattice < 0 .or. system%nlattice > 3) call mopend("Invalid number of lattice vectors") if (system%nlattice_move < 0 .or. system%nlattice_move > system%nlattice) & call mopend("Invalid number of moveable lattice vectors") - if (system%nlattice > 0) then - if(size(system%lattice) < 3*system%nlattice) call mopend("List of lattice vectors is too small") - end if + if (system%nlattice > 0 .and. .not. c_associated(system%lattice)) & + call mopend("List of lattice vectors is too small") if (system%tolerance <= 0.d0) call mopend("Relative numerical tolerance must be a positive number") if (system%max_time <= 0) call mopend("Time limit must be a positive number") if (system%nlattice_move > 0 .and. index(keywrd," FORCETS") /= 0) & call mopend("Lattice vectors cannot move during vibrational calculations") if (moperr) return + ! convert C pointers + call c_f_pointer(system%atom, atom, [system%natom]) + call c_f_pointer(system%coord, coord, [3*system%natom]) + if (system%nlattice > 0) call c_f_pointer(system%lattice, lattice, [3*system%nlattice]) + use_disk = .false. ! load ios, iop, and iod data into tore tore = ios + iop + iod @@ -126,9 +135,10 @@ module subroutine mopac_initialize(system) write(num2str,'(i10)') system%max_time keywrd = trim(keywrd) // " T=" // adjustl(num2str) ! open dummy output file - output_fn = '/dev/null' -#ifdef MOPAC_OS - if (MOPAC_OS == "Windows") output_fn = 'NUL' +#ifdef WIN32 + output_fn = 'NUL' +#else + output_fn = '/dev/null' #endif close(iw) open(unit=iw, file=output_fn, status='UNKNOWN', position='asis', iostat=status) @@ -193,20 +203,20 @@ module subroutine mopac_initialize(system) end if ! insert geometry information into MOPAC data structures id = system%nlattice - nat(:system%natom) = system%atom(:) - labels(:system%natom) = system%atom(:) - geo(:,:system%natom) = reshape(system%coord,[3, system%natom]) + nat(:system%natom) = atom(:system%natom) + labels(:system%natom) = atom(:system%natom) + geo(:,:system%natom) = reshape(coord,[3, system%natom]) if (id > 0) then nat(system%natom+1:system%natom+id) = 107 labels(system%natom+1:system%natom+id) = 107 - geo(:,system%natom+1:system%natom+id) = reshape(system%lattice,[3, id]) + geo(:,system%natom+1:system%natom+id) = reshape(lattice,[3, id]) end if atmass(:natoms) = ams(labels(:)) ! set optimization flags & xparam nvar = 3*system%natom_move + 3*system%nlattice_move lopt(:,:system%natom_move) = 1 lopt(:,system%natom_move+1:) = 0 - xparam(:3*system%natom_move) = system%coord(:3*system%natom_move) + xparam(:3*system%natom_move) = coord(:3*system%natom_move) do i=1, system%natom_move do j=1, 3 loc(1,3*(i-1)+j) = i @@ -215,7 +225,7 @@ module subroutine mopac_initialize(system) end do if (system%nlattice_move > 0) then lopt(:,numat+1:numat+system%nlattice_move) = 1 - xparam(3*system%natom_move+1:) = system%lattice(:3*system%nlattice_move) + xparam(3*system%natom_move+1:) = lattice(:3*system%nlattice_move) do i=1, system%nlattice_move do j=1, 3 loc(1,3*system%natom_move+3*(i-1)+j) = numat+i diff --git a/src/interface/mopac_api_operations.F90 b/src/interface/mopac_api_operations.F90 index 08da0530..27352002 100644 --- a/src/interface/mopac_api_operations.F90 +++ b/src/interface/mopac_api_operations.F90 @@ -16,7 +16,7 @@ submodule (mopac_api) mopac_api_operations use Common_arrays_C, only: xparam, grad, lopt - use molkst_C, only: keywrd, escf, moperr, nvar, gui, jobnam + use molkst_C, only: keywrd, escf, moperr, nvar, gui, jobnam, run implicit none interface @@ -33,7 +33,7 @@ end subroutine mopac_finalize ! save MOPAC density matrices module subroutine mopac_save(state) - type(mopac_state), intent(out) :: state + type(mopac_state), intent(inout) :: state end subroutine mopac_save ! load MOPAC density matrices, or construct initial guesses @@ -43,7 +43,7 @@ end subroutine mopac_load ! save MOZYME density matrix module subroutine mozyme_save(state) - type(mozyme_state), intent(out) :: state + type(mozyme_state), intent(inout) :: state end subroutine mozyme_save ! load MOZYME density matrix, or construct initial guess @@ -56,8 +56,10 @@ end subroutine mozyme_load contains ! MOPAC electronic ground state calculation - module subroutine mopac_scf(system, state, properties) - !dec$ attributes dllexport :: mopac_scf + module subroutine mopac_scf(system, state, properties) bind(c) +#ifdef WIN32 +!dec$ attributes dllexport :: mopac_scf +#endif type(mopac_system), intent(in) :: system type(mopac_state), intent(inout) :: state type(mopac_properties), intent(out) :: properties @@ -72,8 +74,10 @@ module subroutine mopac_scf(system, state, properties) end subroutine mopac_scf ! MOPAC geometry relaxation - module subroutine mopac_relax(system, state, properties) - !dec$ attributes dllexport :: mopac_relax + module subroutine mopac_relax(system, state, properties) bind(c) +#ifdef WIN32 +!dec$ attributes dllexport :: mopac_relax +#endif type(mopac_system), intent(in) :: system type(mopac_state), intent(inout) :: state type(mopac_properties), intent(out) :: properties @@ -88,8 +92,10 @@ module subroutine mopac_relax(system, state, properties) end subroutine mopac_relax ! MOPAC vibrational calculation - module subroutine mopac_vibe(system, state, properties) - !dec$ attributes dllexport :: mopac_vibe + module subroutine mopac_vibe(system, state, properties) bind(c) +#ifdef WIN32 +!dec$ attributes dllexport :: mopac_vibe +#endif type(mopac_system), intent(in) :: system type(mopac_state), intent(inout) :: state type(mopac_properties), intent(out) :: properties @@ -114,8 +120,10 @@ module subroutine mopac_vibe(system, state, properties) end subroutine mopac_vibe ! MOZYME electronic ground state calculation - module subroutine mozyme_scf(system, state, properties) - !dec$ attributes dllexport :: mozyme_scf + module subroutine mozyme_scf(system, state, properties) bind(c) +#ifdef WIN32 +!dec$ attributes dllexport :: mozyme_scf +#endif type(mopac_system), intent(in) :: system type(mozyme_state), intent(inout) :: state type(mopac_properties), intent(out) :: properties @@ -130,8 +138,10 @@ module subroutine mozyme_scf(system, state, properties) end subroutine mozyme_scf ! MOZYME geometry relaxation - module subroutine mozyme_relax(system, state, properties) - !dec$ attributes dllexport :: mozyme_relax + module subroutine mozyme_relax(system, state, properties) bind(c) +#ifdef WIN32 +!dec$ attributes dllexport :: mozyme_relax +#endif type(mopac_system), intent(in) :: system type(mozyme_state), intent(inout) :: state type(mopac_properties), intent(out) :: properties @@ -146,8 +156,10 @@ module subroutine mozyme_relax(system, state, properties) end subroutine mozyme_relax ! MOZYME vibrational calculation - module subroutine mozyme_vibe(system, state, properties) - !dec$ attributes dllexport :: mozyme_vibe + module subroutine mozyme_vibe(system, state, properties) bind(c) +#ifdef WIN32 +!dec$ attributes dllexport :: mozyme_vibe +#endif type(mopac_system), intent(in) :: system type(mozyme_state), intent(inout) :: state type(mopac_properties), intent(out) :: properties @@ -172,14 +184,30 @@ module subroutine mozyme_vibe(system, state, properties) end subroutine mozyme_vibe ! Run MOPAC conventionally from an input file - module subroutine run_mopac_from_input(path_to_file) - !dec$ attributes dllexport :: run_mopac_from_input - character(len=240), intent(in) :: path_to_file - jobnam = trim(path_to_file) + module function run_mopac_from_input(path_to_file) bind(c) +#ifdef WIN32 +!dec$ attributes dllexport :: run_mopac_from_input +#endif + integer(c_int) :: run_mopac_from_input + character(kind=c_char), dimension(*), intent(in) :: path_to_file + integer :: i + i = 1 + do + if(path_to_file(i) == ' ' .or. path_to_file(i) == c_null_char) exit + jobnam(i:i) = path_to_file(i) + i = i + 1 + end do gui = .false. + run = 2 call run_mopac + if (moperr) then + run_mopac_from_input = 1 + else + run_mopac_from_input = 0 + end if + run = 1 gui = .true. jobnam = ' ' - end subroutine run_mopac_from_input + end function run_mopac_from_input end submodule mopac_api_operations diff --git a/src/interface/mopac_api_saveload.F90 b/src/interface/mopac_api_saveload.F90 index cb87131f..2bd4ec4e 100644 --- a/src/interface/mopac_api_saveload.F90 +++ b/src/interface/mopac_api_saveload.F90 @@ -25,37 +25,35 @@ ! save MOPAC density matrices module subroutine mopac_save(state) - type(mopac_state), intent(out) :: state - integer :: status + type(mopac_state), intent(inout) :: state - state%save_state = .true. - state%mpack = mpack + if (state%mpack > 0) call destroy_mopac_state(state) - if (allocated(state%pa)) deallocate(state%pa) - allocate(state%pa(mpack), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOPAC_SAVE") - return - end if - state%pa = pa - if (uhf) then - if (allocated(state%pb)) deallocate(state%pb) - allocate(state%pb(mpack), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOPAC_SAVE") - return + state%mpack = mpack + if (mpack > 0) then + state%pa = create_copy(pa, [mpack]) + if (uhf) then + state%uhf = 1 + state%pb = create_copy(pb, [mpack]) + else + state%uhf = 0 + state%pb = c_null_ptr end if - state%pb = pb end if + end subroutine mopac_save ! load MOPAC density matrices, or construct initial guesses module subroutine mopac_load(state) type(mopac_state), intent(in) :: state integer :: status + real(c_double), pointer :: rptr(:) - if(state%save_state) then - ! TO DO: compatibility tests + if(state%mpack > 0) then + if (state%mpack /= mpack .or. (state%uhf /= 0 .neqv. uhf)) then + call mopend("Attempting to load incompatible MOPAC state") + return + end if keywrd = trim(keywrd) // " OLDENS" mpack = state%mpack if (allocated(pa)) deallocate(pa) @@ -63,8 +61,9 @@ module subroutine mopac_load(state) if (status /= 0) then call mopend("Failed to allocate memory in MOPAC_LOAD") return - end if - pa = state%pa + end if + call c_f_pointer(state%pa, rptr, [mpack]) + pa = rptr if(uhf) then if (allocated(pb)) deallocate(pb) allocate(pb(mpack), stat=status) @@ -72,96 +71,51 @@ module subroutine mopac_load(state) call mopend("Failed to allocate memory in MOPAC_LOAD") return end if - pb = state%pb - end if + call c_f_pointer(state%pb, rptr, [mpack]) + pb = rptr + end if end if end subroutine mopac_load ! save MOZYME density matrix module subroutine mozyme_save(state) - type(mozyme_state), intent(out) :: state - integer :: status + type(mozyme_state), intent(inout) :: state + + if (state%numat > 0) call destroy_mozyme_state(state) - state%save_state = .true. state%numat = numat - state%noccupied = noccupied - state%nvirtual = nvirtual - state%icocc_dim = icocc_dim - state%icvir_dim = icvir_dim - state%cocc_dim = cocc_dim - state%cvir_dim = cvir_dim - if (allocated(state%nbonds)) deallocate(state%nbonds) - if (allocated(state%ibonds)) deallocate(state%ibonds) - if (allocated(state%iorbs)) deallocate(state%iorbs) - if (allocated(state%ncf)) deallocate(state%ncf) - if (allocated(state%nce)) deallocate(state%nce) - if (allocated(state%icocc)) deallocate(state%icocc) - if (allocated(state%icvir)) deallocate(state%icvir) - if (allocated(state%cocc)) deallocate(state%cocc) - if (allocated(state%cvir)) deallocate(state%cvir) - allocate(state%nbonds(numat), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOZYME_SAVE") - return - end if - allocate(state%ibonds(9,numat), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOZYME_SAVE") - return - end if - allocate(state%iorbs(numat), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOZYME_SAVE") - return - end if - allocate(state%ncf(noccupied), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOZYME_SAVE") - return - end if - allocate(state%nce(nvirtual), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOZYME_SAVE") - return - end if - allocate(state%icocc(icocc_dim), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOZYME_SAVE") - return - end if - allocate(state%icvir(icvir_dim), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOZYME_SAVE") - return - end if - allocate(state%cocc(cocc_dim), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOZYME_SAVE") - return - end if - allocate(state%cvir(cvir_dim), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOZYME_SAVE") - return + if (numat > 0) then + state%noccupied = noccupied + state%nvirtual = nvirtual + state%icocc_dim = icocc_dim + state%icvir_dim = icvir_dim + state%cocc_dim = cocc_dim + state%cvir_dim = cvir_dim + + state%nbonds = create_copy(nbonds, [numat]) + state%ibonds = create_copy(ibonds, [9,numat]) + state%iorbs = create_copy(iorbs, [numat]) + state%ncf = create_copy(ncf, [noccupied]) + state%nce = create_copy(nce, [nvirtual]) + state%icocc = create_copy(icocc, [icocc_dim]) + state%icvir = create_copy(icvir, [icvir_dim]) + state%cocc = create_copy(cocc, [cocc_dim]) + state%cvir = create_copy(cvir, [cvir_dim]) end if - state%nbonds = nbonds - state%ibonds = ibonds - state%iorbs = iorbs - state%ncf = ncf - state%nce = nce - state%icocc = icocc - state%icvir = icvir - state%cocc = cocc - state%cvir = cvir end subroutine mozyme_save ! load MOZYME density matrix, or construct initial guess module subroutine mozyme_load(state) type(mozyme_state), intent(in) :: state integer :: status, i, j, k, l + integer(c_int), pointer :: iptr(:), iptr2(:,:) + real(c_double), pointer :: rptr(:) - if(state%save_state) then - ! TO DO: compatibility tests + if(state%numat > 0) then + if (state%numat /= numat) then + call mopend("Attempting to load incompatible MOZYME state") + return + end if keywrd = trim(keywrd) // " OLDENS" numat = state%numat noccupied = state%noccupied @@ -248,15 +202,24 @@ module subroutine mozyme_load(state) call mopend("Failed to allocate memory in MOZYME_LOAD") return end if - nbonds = state%nbonds - ibonds = state%ibonds - iorbs = state%iorbs - ncf = state%ncf - nce = state%nce - icocc = state%icocc - icvir = state%icvir - cocc = state%cocc - cvir = state%cvir + call c_f_pointer(state%nbonds, iptr, [numat]) + nbonds = iptr + call c_f_pointer(state%ibonds, iptr2, [9,numat]) + ibonds = iptr2 + call c_f_pointer(state%iorbs, iptr, [numat]) + iorbs = iptr + call c_f_pointer(state%ncf, iptr, [noccupied]) + ncf = iptr + call c_f_pointer(state%nce, iptr, [nvirtual]) + nce = iptr + call c_f_pointer(state%icocc, iptr, [icocc_dim]) + icocc = iptr + call c_f_pointer(state%icvir, iptr, [icvir_dim]) + icvir = iptr + call c_f_pointer(state%cocc, rptr, [cocc_dim]) + cocc = rptr + call c_f_pointer(state%cvir, rptr, [cvir_dim]) + cvir = rptr ! reconstruct nncf, nnce, ncocc, & ncvir j = 0 do i = 1, noccupied diff --git a/src/molkst_C.F90 b/src/molkst_C.F90 index dbe05625..0f11f27d 100644 --- a/src/molkst_C.F90 +++ b/src/molkst_C.F90 @@ -15,7 +15,9 @@ ! along with this program. If not, see . module molkst_C +#ifdef WIN32 !dec$ attributes dllexport :: gui +#endif ! ! This module contains all the scalars relating to the system being calculated. ! Every entry is unique, and has the same meaning in every subroutine. diff --git a/src/run_mopac.F90 b/src/run_mopac.F90 index 68e35824..00efd586 100644 --- a/src/run_mopac.F90 +++ b/src/run_mopac.F90 @@ -15,7 +15,9 @@ ! along with this program. If not, see . subroutine run_mopac +#ifdef WIN32 !dec$ attributes dllexport :: run_mopac +#endif !----------------------------------------------- ! M o d u l e s !----------------------------------------------- @@ -32,7 +34,7 @@ subroutine run_mopac sparkle, itemp_1, maxtxt, koment, sz, ss2, keywrd_quoted, & nl_atoms, use_ref_geo, prt_coords, pdb_label, step, & density, norbs, method_indo, nclose, nopen, backslash, gui, os, git_hash, verson, & - use_disk + use_disk, run ! USE parameters_C, only : tore, ios, iop, iod, eisol, eheat, zs, eheat_sparkles, gss ! @@ -1030,8 +1032,11 @@ subroutine run_mopac write (iw, '(3/,'' TOTAL JOB TIME: '',F16.2,'' SECONDS'')') tim write (iw, '(/,'' == MOPAC DONE =='')') call fdate (line) - write(*,'(//10x,a,/)')"MOPAC Job: """//trim(job_fn)//""" ended normally on "// & - line(5:10)//", "//trim(line(21:))//", at"//line(11:16)//"." + ! suppress stdout when running as an API call + if (run /= 2) then + write(*,'(//10x,a,/)')"MOPAC Job: """//trim(job_fn)//""" ended normally on "// & + line(5:10)//", "//trim(line(21:))//", at"//line(11:16)//"." + end if ! ! Delete files that are definitely not wanted ! diff --git a/src/run_param.F90 b/src/run_param.F90 index 0c3d22b2..840b8708 100644 --- a/src/run_param.F90 +++ b/src/run_param.F90 @@ -15,7 +15,9 @@ ! along with this program. If not, see . subroutine run_param +#ifdef WIN32 !dec$ attributes dllexport :: run_param +#endif use param_global_C, only : numvar, large, power, ifiles_8, & & contrl, fnsnew, nfns, xparamp, diffns, fns, penalty, & save_parameters diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 8ecb9334..04b48ae2 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -31,17 +31,7 @@ macro(add_mopac_test _name _files) ${CMAKE_BINARY_DIR}/mopac${CMAKE_EXECUTABLE_SUFFIX} ${HOF_ERROR} ${_files}) endmacro() -#=============================================== -# Create a list of source files (src_list) with the .F90 extension -set(src_list - mopac_api_test - ) -#----------------------------------------------- -# Add a list of source files to the target -foreach(idx IN LISTS src_list) - target_sources(mopac-api-test PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/${idx}.F90) -endforeach() -#=============================================== +target_sources(mopac-api-test PRIVATE mopac_api_test.F90) add_subdirectory(INDO-dev) add_subdirectory(keywords) diff --git a/tests/mopac_api_test.F90 b/tests/mopac_api_test.F90 index 2238cb18..7ac27ab2 100644 --- a/tests/mopac_api_test.F90 +++ b/tests/mopac_api_test.F90 @@ -15,32 +15,42 @@ ! along with this program. If not, see . program mopac_api_test - use mopac_api + use mopac_api_f implicit none integer :: nfail nfail = 0 call test_mopac_scf1(nfail) + write(*,*) "Passed API test 1" call test_mopac_relax1(nfail) + write(*,*) "Passed API test 2" call test_mopac_vibe1(nfail) + write(*,*) "Passed API test 3" call test_mozyme_scf1(nfail) + write(*,*) "Passed API test 4" call test_mozyme_relax1(nfail) + write(*,*) "Passed API test 5" call test_mozyme_vibe1(nfail) + write(*,*) "Passed API test 6" call test_mopac_restart1(nfail) + write(*,*) "Passed API test 7" call test_mozyme_restart1(nfail) + write(*,*) "Passed API test 8" call test_cosmo1(nfail) + write(*,*) "Passed API test 9" call test_crystal1(nfail) + write(*,*) "Passed API test 10" call exit(nfail) end program mopac_api_test subroutine test_mopac_scf1(nfail) - use mopac_api + use mopac_api_f implicit none integer, intent(inout) :: nfail - type(mopac_system) :: test_in - type(mopac_state) :: test_restore - type(mopac_properties) :: test_target - type(mopac_properties) :: test_out + type(mopac_system_f) :: test_in + type(mopac_state_f) :: test_restore + type(mopac_properties_f) :: test_target + type(mopac_properties_f) :: test_out character(50) :: test_name integer :: i test_name = 'H2O SCF' @@ -104,21 +114,20 @@ subroutine test_mopac_scf1(nfail) test_target%bond_order(5) = 0.895d0 test_target%bond_order(6) = 0.895d0 test_target%bond_order(7) = 1.791d0 - test_target%calc_vibe = .false. test_target%nerror = 0 - call mopac_scf(test_in, test_restore, test_out) + call mopac_scf_f(test_in, test_restore, test_out) call test_output(test_name, test_in, test_target, test_out, nfail) end subroutine test_mopac_scf1 subroutine test_mopac_relax1(nfail) - use mopac_api + use mopac_api_f implicit none integer, intent(inout) :: nfail - type(mopac_system) :: test_in - type(mopac_state) :: test_restore - type(mopac_properties) :: test_target - type(mopac_properties) :: test_out + type(mopac_system_f) :: test_in + type(mopac_state_f) :: test_restore + type(mopac_properties_f) :: test_target + type(mopac_properties_f) :: test_out character(50) :: test_name integer :: i test_name = 'H2O relax' @@ -190,21 +199,20 @@ subroutine test_mopac_relax1(nfail) test_target%bond_order(5) = 0.894d0 test_target%bond_order(6) = 0.894d0 test_target%bond_order(7) = 1.788d0 - test_target%calc_vibe = .false. test_target%nerror = 0 - call mopac_relax(test_in, test_restore, test_out) + call mopac_relax_f(test_in, test_restore, test_out) call test_output(test_name, test_in, test_target, test_out, nfail) end subroutine test_mopac_relax1 subroutine test_mopac_vibe1(nfail) - use mopac_api + use mopac_api_f implicit none integer, intent(inout) :: nfail - type(mopac_system) :: test_in - type(mopac_state) :: test_restore - type(mopac_properties) :: test_target - type(mopac_properties) :: test_out + type(mopac_system_f) :: test_in + type(mopac_state_f) :: test_restore + type(mopac_properties_f) :: test_target + type(mopac_properties_f) :: test_out character(50) :: test_name test_name = 'H2O vibe' @@ -267,7 +275,6 @@ subroutine test_mopac_vibe1(nfail) test_target%bond_order(5) = 0.894d0 test_target%bond_order(6) = 0.894d0 test_target%bond_order(7) = 1.788d0 - test_target%calc_vibe = .true. test_target%nerror = 0 allocate(test_target%freq(9)) test_target%freq(1) = -18.0d0 @@ -362,18 +369,18 @@ subroutine test_mopac_vibe1(nfail) test_target%disp(8,9) = 0.1483d0 test_target%disp(9,9) = 0.0000d0 - call mopac_vibe(test_in, test_restore, test_out) + call mopac_vibe_f(test_in, test_restore, test_out) call test_output(test_name, test_in, test_target, test_out, nfail) end subroutine test_mopac_vibe1 subroutine test_mozyme_scf1(nfail) - use mopac_api + use mopac_api_f implicit none integer, intent(inout) :: nfail - type(mopac_system) :: test_in - type(mozyme_state) :: test_restore - type(mopac_properties) :: test_target - type(mopac_properties) :: test_out + type(mopac_system_f) :: test_in + type(mozyme_state_f) :: test_restore + type(mopac_properties_f) :: test_target + type(mopac_properties_f) :: test_out character(50) :: test_name integer :: i test_name = 'H2O SCF MOZYME' @@ -437,21 +444,20 @@ subroutine test_mozyme_scf1(nfail) test_target%bond_order(5) = 0.896d0 test_target%bond_order(6) = 0.896d0 test_target%bond_order(7) = 1.793d0 - test_target%calc_vibe = .false. test_target%nerror = 0 - call mozyme_scf(test_in, test_restore, test_out) + call mozyme_scf_f(test_in, test_restore, test_out) call test_output(test_name, test_in, test_target, test_out, nfail) end subroutine test_mozyme_scf1 subroutine test_mozyme_relax1(nfail) - use mopac_api + use mopac_api_f implicit none integer, intent(inout) :: nfail - type(mopac_system) :: test_in - type(mozyme_state) :: test_restore - type(mopac_properties) :: test_target - type(mopac_properties) :: test_out + type(mopac_system_f) :: test_in + type(mozyme_state_f) :: test_restore + type(mopac_properties_f) :: test_target + type(mopac_properties_f) :: test_out character(50) :: test_name integer :: i test_name = 'H2O relax MOZYME' @@ -523,22 +529,21 @@ subroutine test_mozyme_relax1(nfail) test_target%bond_order(5) = 0.894d0 test_target%bond_order(6) = 0.894d0 test_target%bond_order(7) = 1.788d0 - test_target%calc_vibe = .false. test_target%nerror = 0 - call mozyme_relax(test_in, test_restore, test_out) + call mozyme_relax_f(test_in, test_restore, test_out) call test_output(test_name, test_in, test_target, test_out, nfail) end subroutine test_mozyme_relax1 subroutine test_mozyme_vibe1(nfail) - use mopac_api + use mopac_api_f implicit none integer, intent(inout) :: nfail - type(mopac_system) :: test_in - type(mozyme_state) :: test_restore - type(mopac_properties) :: test_target - type(mopac_properties) :: test_out + type(mopac_system_f) :: test_in + type(mozyme_state_f) :: test_restore + type(mopac_properties_f) :: test_target + type(mopac_properties_f) :: test_out character(50) :: test_name test_name = 'H2O vibe MOZYME' @@ -602,7 +607,6 @@ subroutine test_mozyme_vibe1(nfail) test_target%bond_order(5) = 0.894d0 test_target%bond_order(6) = 0.894d0 test_target%bond_order(7) = 1.788d0 - test_target%calc_vibe = .true. test_target%nerror = 0 allocate(test_target%freq(9)) test_target%freq(1) = -18.0d0 @@ -697,18 +701,18 @@ subroutine test_mozyme_vibe1(nfail) test_target%disp(8,9) = 0.1483d0 test_target%disp(9,9) = 0.0000d0 - call mozyme_vibe(test_in, test_restore, test_out) + call mozyme_vibe_f(test_in, test_restore, test_out) call test_output(test_name, test_in, test_target, test_out, nfail) end subroutine test_mozyme_vibe1 subroutine test_mopac_restart1(nfail) - use mopac_api + use mopac_api_f implicit none integer, intent(inout) :: nfail - type(mopac_system) :: test_in - type(mopac_state) :: test_restore - type(mopac_properties) :: test_target - type(mopac_properties) :: test_out + type(mopac_system_f) :: test_in + type(mopac_state_f) :: test_restore + type(mopac_properties_f) :: test_target + type(mopac_properties_f) :: test_out character(50) :: test_name integer :: i test_name = 'H2O SCF restart' @@ -772,22 +776,21 @@ subroutine test_mopac_restart1(nfail) test_target%bond_order(5) = 0.895d0 test_target%bond_order(6) = 0.895d0 test_target%bond_order(7) = 1.791d0 - test_target%calc_vibe = .false. test_target%nerror = 0 - call mopac_scf(test_in, test_restore, test_out) - call mopac_scf(test_in, test_restore, test_out) + call mopac_scf_f(test_in, test_restore, test_out) + call mopac_scf_f(test_in, test_restore, test_out) call test_output(test_name, test_in, test_target, test_out, nfail) end subroutine test_mopac_restart1 subroutine test_mozyme_restart1(nfail) - use mopac_api + use mopac_api_f implicit none integer, intent(inout) :: nfail - type(mopac_system) :: test_in - type(mozyme_state) :: test_restore - type(mopac_properties) :: test_target - type(mopac_properties) :: test_out + type(mopac_system_f) :: test_in + type(mozyme_state_f) :: test_restore + type(mopac_properties_f) :: test_target + type(mopac_properties_f) :: test_out character(50) :: test_name integer :: i test_name = 'H2O SCF MOZYME restart' @@ -854,22 +857,21 @@ subroutine test_mozyme_restart1(nfail) test_target%bond_order(5) = 0.895d0 test_target%bond_order(6) = 0.895d0 test_target%bond_order(7) = 1.791d0 - test_target%calc_vibe = .false. test_target%nerror = 0 - call mozyme_scf(test_in, test_restore, test_out) - call mozyme_scf(test_in, test_restore, test_out) + call mozyme_scf_f(test_in, test_restore, test_out) + call mozyme_scf_f(test_in, test_restore, test_out) call test_output(test_name, test_in, test_target, test_out, nfail) end subroutine test_mozyme_restart1 subroutine test_cosmo1(nfail) - use mopac_api + use mopac_api_f implicit none integer, intent(inout) :: nfail - type(mopac_system) :: test_in - type(mopac_state) :: test_restore - type(mopac_properties) :: test_target - type(mopac_properties) :: test_out + type(mopac_system_f) :: test_in + type(mopac_state_f) :: test_restore + type(mopac_properties_f) :: test_target + type(mopac_properties_f) :: test_out character(50) :: test_name integer :: i test_name = 'H2O SCF COSMO' @@ -934,22 +936,20 @@ subroutine test_cosmo1(nfail) test_target%bond_order(5) = 0.859d0 test_target%bond_order(6) = 0.859d0 test_target%bond_order(7) = 1.718d0 - test_target%calc_vibe = .false. test_target%nerror = 0 - call mopac_scf(test_in, test_restore, test_out) + call mopac_scf_f(test_in, test_restore, test_out) call test_output(test_name, test_in, test_target, test_out, nfail) end subroutine test_cosmo1 - subroutine test_crystal1(nfail) - use mopac_api + use mopac_api_f implicit none integer, intent(inout) :: nfail - type(mopac_system) :: test_in - type(mopac_state) :: test_restore - type(mopac_properties) :: test_target - type(mopac_properties) :: test_out + type(mopac_system_f) :: test_in + type(mopac_state_f) :: test_restore + type(mopac_properties_f) :: test_target + type(mopac_properties_f) :: test_out character(50) :: test_name integer :: i test_name = 'LiF SCF' @@ -2189,19 +2189,18 @@ subroutine test_crystal1(nfail) test_target%bond_order(376) = 0.037d0 test_target%bond_order(377) = 0.037d0 test_target%bond_order(378) = 0.283d0 - test_target%calc_vibe = .false. test_target%nerror = 0 - call mopac_scf(test_in, test_restore, test_out) + call mopac_scf_f(test_in, test_restore, test_out) call test_output(test_name, test_in, test_target, test_out, nfail) end subroutine test_crystal1 subroutine test_output(name, input, target, output, nfail) - use mopac_api + use mopac_api_f implicit none character(50), intent(in) :: name - type(mopac_system), intent(in) :: input - type(mopac_properties), intent(in) :: target, output + type(mopac_system_f), intent(in) :: input + type(mopac_properties_f), intent(in) :: target, output integer, intent(inout) :: nfail integer :: i, j double precision :: heat_tol, coord_tol, deriv_tol, freq_tol, charge_tol, stress_tol @@ -2250,11 +2249,11 @@ subroutine test_output(name, input, target, output, nfail) end if end do ! compare vibrational properties - if(target%calc_vibe .neqv. output%calc_vibe) then + if(allocated(target%freq) .neqv. allocated(output%freq)) then nfail = nfail + 1 - write(*,*) "calc_vibe mistmatch in test '", trim(name), "':", target%calc_vibe, & - "vs", output%calc_vibe - else if(target%calc_vibe .eqv. .true.) then + write(*,*) "freq allocation mistmatch in test '", trim(name), "':", allocated(target%freq), & + "vs", allocated(output%freq) + else if(allocated(target%freq) .eqv. .true.) then do i = 1, 3*input%natom_move if(abs(target%freq(i) - output%freq(i)) > freq_tol) then nfail = nfail + 1