From 40151d3b2cd9268cc476245818a1167b2aefa68a Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Wed, 23 Oct 2024 14:22:49 -0400 Subject: [PATCH 01/30] Switch MOPAC API to strict C binding --- src/interface/mopac_api.F90 | 233 ++++++++++++++----------- src/interface/mopac_api_operations.F90 | 22 +-- 2 files changed, 144 insertions(+), 111 deletions(-) diff --git a/src/interface/mopac_api.F90 b/src/interface/mopac_api.F90 index 19316b4c..8bc5112a 100644 --- a/src/interface/mopac_api.F90 +++ b/src/interface/mopac_api.F90 @@ -15,195 +15,226 @@ ! 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 iso_c_binding implicit none private ! public derived types of the MOPAC API public :: mopac_system, mopac_properties, mopac_state, mozyme_state + ! public subroutines to create or destroy derived types + public :: create_mopac_system, destroy_mopac_system, destroy_mopac_properties, & + destroy_mopac_state, destroy_mozyme_state ! public subroutines of the MOPAC API public :: mopac_scf, mopac_relax, mopac_vibe, mozyme_scf, mozyme_relax, mozyme_vibe ! 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 + integer(c_int) :: spin ! dielectric constant for COSMO implicit solvent, must be 1 (no solvent) for nlattice > 0 - double precision :: epsilon = 1.d0 + real(c_double) :: epsilon ! 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 + ! 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 + 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 ! 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)-1] + ! > bond order of atoms bonded to each atom in CSC format + type(c_ptr) :: bond_order ! real(c_double)[bond_index(natom+1)-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 + ! > alpha density matrix + type(c_ptr) :: pa ! real(c_double)[mpack] + ! > beta density matrix + type(c_ptr) :: pb ! real(c_double)[mpack] 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 + 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 [numat] - integer, dimension (:), allocatable :: iorbs + 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 + ! allocate memory & initialize mopac_system + module subroutine create_mopac_system(natom, natom_move, charge, spin, epsilon, model, & + atom, coord, nlattice, nlattice_move, pressure, lattice, tolerance, max_time, & + system) bind(c) + integer(c_int), intent(in) :: natom + integer(c_int), intent(in) :: natom_move + integer(c_int), intent(in) :: charge + integer(c_int), intent(in) :: spin + real(c_double), intent(in) :: epsilon + integer(c_int), intent(in) :: model + integer(c_int), dimension(natom), intent(in) :: atom + real(c_double), dimension(3*natom), intent(in) :: coord + integer(c_int), intent(in) :: nlattice + integer(c_int), intent(in) :: nlattice_move + real(c_double), intent(in) :: pressure + real(c_double), dimension(3*nlattice), intent(in) :: lattice + real(c_double), intent(in) :: tolerance + integer(c_int), intent(in) :: max_time + type(mopac_system), intent(out) :: system + end subroutine create_mopac_system + + ! deallocate memory in mopac_system + module subroutine destroy_mopac_system(system) bind(c) + type(mopac_system), intent(in) :: system + end subroutine destroy_mopac_system + + ! deallocate memory in mopac_properties + module subroutine destroy_mopac_properties(properties) bind(c) + type(mopac_properties), intent(in) :: properties + end subroutine destroy_mopac_properties + + ! deallocate memory in mopac_state + module subroutine destroy_mopac_state(state) bind(c) + type(mopac_state), intent(in) :: state + end subroutine destroy_mopac_state + + ! deallocate memory in mozyme_state + module subroutine destroy_mozyme_state(state) bind(c) + type(mozyme_state), intent(in) :: state + end subroutine destroy_mozyme_state + ! 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) 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) 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) 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) 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) 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) 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 + module subroutine run_mopac_from_input(path_to_file) bind(c) + character(kind=c_char,len=*), intent(in) :: path_to_file end subroutine run_mopac_from_input end interface diff --git a/src/interface/mopac_api_operations.F90 b/src/interface/mopac_api_operations.F90 index 08da0530..55d9ffbc 100644 --- a/src/interface/mopac_api_operations.F90 +++ b/src/interface/mopac_api_operations.F90 @@ -15,8 +15,9 @@ ! along with this program. If not, see . submodule (mopac_api) mopac_api_operations + use iso_c_binding 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 @@ -57,7 +58,6 @@ end subroutine mozyme_load ! MOPAC electronic ground state calculation module subroutine mopac_scf(system, state, properties) - !dec$ attributes dllexport :: mopac_scf type(mopac_system), intent(in) :: system type(mopac_state), intent(inout) :: state type(mopac_properties), intent(out) :: properties @@ -73,7 +73,6 @@ end subroutine mopac_scf ! MOPAC geometry relaxation module subroutine mopac_relax(system, state, properties) - !dec$ attributes dllexport :: mopac_relax type(mopac_system), intent(in) :: system type(mopac_state), intent(inout) :: state type(mopac_properties), intent(out) :: properties @@ -89,7 +88,6 @@ end subroutine mopac_relax ! MOPAC vibrational calculation module subroutine mopac_vibe(system, state, properties) - !dec$ attributes dllexport :: mopac_vibe type(mopac_system), intent(in) :: system type(mopac_state), intent(inout) :: state type(mopac_properties), intent(out) :: properties @@ -115,7 +113,6 @@ end subroutine mopac_vibe ! MOZYME electronic ground state calculation module subroutine mozyme_scf(system, state, properties) - !dec$ attributes dllexport :: mozyme_scf type(mopac_system), intent(in) :: system type(mozyme_state), intent(inout) :: state type(mopac_properties), intent(out) :: properties @@ -131,7 +128,6 @@ end subroutine mozyme_scf ! MOZYME geometry relaxation module subroutine mozyme_relax(system, state, properties) - !dec$ attributes dllexport :: mozyme_relax type(mopac_system), intent(in) :: system type(mozyme_state), intent(inout) :: state type(mopac_properties), intent(out) :: properties @@ -147,7 +143,6 @@ end subroutine mozyme_relax ! MOZYME vibrational calculation module subroutine mozyme_vibe(system, state, properties) - !dec$ attributes dllexport :: mozyme_vibe type(mopac_system), intent(in) :: system type(mozyme_state), intent(inout) :: state type(mopac_properties), intent(out) :: properties @@ -173,11 +168,18 @@ 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) + character(kind=c_char,len=*), 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) = path_to_file(i) + i = i + 1 + end do gui = .false. + run = 2 call run_mopac + run = 1 gui = .true. jobnam = ' ' end subroutine run_mopac_from_input From 55d02bbc3ed8fa6316d34836791c4a28fe73c5c4 Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Wed, 23 Oct 2024 15:32:38 -0400 Subject: [PATCH 02/30] Create C header file to specify API --- include/mopac_api.h | 190 ++++++++++++++++++++++++++++++++++++ src/interface/mopac_api.F90 | 17 ++-- 2 files changed, 199 insertions(+), 8 deletions(-) create mode 100644 include/mopac_api.h diff --git a/include/mopac_api.h b/include/mopac_api.h new file mode 100644 index 00000000..c55ee886 --- /dev/null +++ b/include/mopac_api.h @@ -0,0 +1,190 @@ +/* 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 */ +#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 + * 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]-1] */ + /* > bond order of atoms bonded to each atom in CSC format */ + double *bond_order; /* [bond_index[natom]-1] */ + /* (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 */ + /* > alpha density matrix */ + double *pa; /* [mpack] */ + /* > beta density matrix */ + double *pb; /* [mpack] */ +}; + +/* 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] */ +}; + +/* allocate memory & initialize mopac_system */ +void create_mopac_system(int *natom, int *natom_move, int *charge, int *spin, int *model, + double *epsilon, int *atom, double *coord, int *nlattice, int *nlattice_move, + double *pressure, double *lattice, double *tolerance, int *max_time, + struct mopac_system *system); + +/* deallocate memory in mopac_system */ +void destroy_mopac_system(struct mopac_system *system); + +/* 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); + +/* 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); + +/* run MOPAC conventionally from an input file */ +void run_mopac_from_input(char *path_to_file); + +#endif \ No newline at end of file diff --git a/src/interface/mopac_api.F90 b/src/interface/mopac_api.F90 index 8bc5112a..336d0e16 100644 --- a/src/interface/mopac_api.F90 +++ b/src/interface/mopac_api.F90 @@ -40,10 +40,10 @@ module mopac_api integer(c_int) :: charge ! number of spin excitations, floor[(number of alpha electrons)/2 - (number of beta electrons)/2] integer(c_int) :: spin - ! dielectric constant for COSMO implicit solvent, must be 1 (no solvent) for nlattice > 0 - real(c_double) :: epsilon ! 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) @@ -119,7 +119,7 @@ module mopac_api 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 [numat] + ! > number of orbitals per real atom type(c_ptr) :: iorbs ! integer(c_int)[numat] ! > number of occupied molecular orbitals integer(c_int) :: noccupied @@ -150,15 +150,16 @@ module mopac_api interface ! allocate memory & initialize mopac_system - module subroutine create_mopac_system(natom, natom_move, charge, spin, epsilon, model, & - atom, coord, nlattice, nlattice_move, pressure, lattice, tolerance, max_time, & - system) bind(c) + module subroutine create_mopac_system(natom, natom_move, charge, spin, model, & + epsilon, atom, coord, nlattice, nlattice_move, & + pressure, lattice, tolerance, max_time, & + system) bind(c) integer(c_int), intent(in) :: natom integer(c_int), intent(in) :: natom_move integer(c_int), intent(in) :: charge integer(c_int), intent(in) :: spin - real(c_double), intent(in) :: epsilon integer(c_int), intent(in) :: model + real(c_double), intent(in) :: epsilon integer(c_int), dimension(natom), intent(in) :: atom real(c_double), dimension(3*natom), intent(in) :: coord integer(c_int), intent(in) :: nlattice @@ -232,7 +233,7 @@ module subroutine mozyme_vibe(system, state, properties) bind(c) type(mopac_properties), intent(out) :: properties end subroutine mozyme_vibe - ! Run MOPAC conventionally from an input file + ! run MOPAC conventionally from an input file module subroutine run_mopac_from_input(path_to_file) bind(c) character(kind=c_char,len=*), intent(in) :: path_to_file end subroutine run_mopac_from_input From 189c891164361f0af9529adc84cfc23117cf1185 Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Wed, 23 Oct 2024 16:30:58 -0400 Subject: [PATCH 03/30] Package API header in installer --- CMakeLists.txt | 1 + src/interface/mopac_api.F90 | 7 ++++--- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index d69859ec..b5e1c475 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -319,6 +319,7 @@ 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(PUBLIC_HEADER "${CMAKE_SOURCE_DIR}/include/mopac_api.h" COMPONENT redist EXCLUDE_FROM_ALL) # Install the executables and library install(TARGETS mopac RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR} COMPONENT main) diff --git a/src/interface/mopac_api.F90 b/src/interface/mopac_api.F90 index 336d0e16..4ba7a84e 100644 --- a/src/interface/mopac_api.F90 +++ b/src/interface/mopac_api.F90 @@ -22,9 +22,10 @@ module mopac_api private ! public derived types of the MOPAC API public :: mopac_system, mopac_properties, mopac_state, mozyme_state - ! public subroutines to create or destroy derived types - public :: create_mopac_system, destroy_mopac_system, destroy_mopac_properties, & - destroy_mopac_state, destroy_mozyme_state + ! public subroutines to create or destroy a MOPAC system + public :: create_mopac_system, destroy_mopac_system + ! public subroutines to destroy the other derived types of the API + public :: destroy_mopac_properties, destroy_mopac_state, destroy_mozyme_state ! public subroutines of the MOPAC API public :: mopac_scf, mopac_relax, mopac_vibe, mozyme_scf, mozyme_relax, mozyme_vibe ! public subroutine of the simple, disk-based MOPAC API From 98b4bc9e68e3c2836b69cf9b986a79b1211cac66 Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Wed, 23 Oct 2024 19:23:01 -0400 Subject: [PATCH 04/30] Implement create/destroy API subroutines --- src/interface/CMakeLists.txt | 2 +- src/interface/mopac_api.F90 | 2 +- src/interface/mopac_api_createdestroy.F90 | 186 ++++++++++++++++++++++ 3 files changed, 188 insertions(+), 2 deletions(-) create mode 100644 src/interface/mopac_api_createdestroy.F90 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 4ba7a84e..9f4c7cd2 100644 --- a/src/interface/mopac_api.F90 +++ b/src/interface/mopac_api.F90 @@ -16,7 +16,7 @@ ! Diskless/stateless Application Programming Interface (API) to core MOPAC operations module mopac_api - use iso_c_binding + use iso_c_binding, only: c_int, c_double, c_char, c_ptr implicit none private diff --git a/src/interface/mopac_api_createdestroy.F90 b/src/interface/mopac_api_createdestroy.F90 new file mode 100644 index 00000000..51a18865 --- /dev/null +++ b/src/interface/mopac_api_createdestroy.F90 @@ -0,0 +1,186 @@ +! 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 + use iso_c_binding, only: c_int, c_double, c_char, c_ptr, c_loc, c_f_pointer, c_associated + implicit none + +contains + + ! allocate memory & initialize mopac_system + module subroutine create_mopac_system(natom, natom_move, charge, spin, model, & + epsilon, atom, coord, nlattice, nlattice_move, & + pressure, lattice, tolerance, max_time, & + system) bind(c) + integer(c_int), intent(in) :: natom + integer(c_int), intent(in) :: natom_move + integer(c_int), intent(in) :: charge + integer(c_int), intent(in) :: spin + integer(c_int), intent(in) :: model + real(c_double), intent(in) :: epsilon + integer(c_int), dimension(natom), intent(in) :: atom + real(c_double), dimension(3*natom), intent(in) :: coord + integer(c_int), intent(in) :: nlattice + integer(c_int), intent(in) :: nlattice_move + real(c_double), intent(in) :: pressure + real(c_double), dimension(3*nlattice), intent(in) :: lattice + real(c_double), intent(in) :: tolerance + integer(c_int), intent(in) :: max_time + type(mopac_system), intent(out) :: system + integer(c_int), pointer :: iptr(:) + real(c_double), pointer :: rptr(:) + system%natom = natom + system%natom_move = natom_move + system%charge = charge + system%spin = spin + system%model = model + system%epsilon = epsilon + if (natom > 0) then + allocate(iptr(natom)) + iptr(1:natom) = atom(1:natom) + system%atom = call c_loc(iptr) + allocate(rptr(3*natom)) + rptr(1:3*natom) = coord(1:3*natom) + system%coord = call c_loc(rptr) + end if + system%nlattice = nlattice + system%nlattice_move = nlattice_move + system%pressure = pressure + if (nlattice > 0) then + allocate(rptr(3*nlattice)) + rptr(1:3*nlattice) = lattice(1:3*nlattice) + system%lattice = call c_loc(rptr) + end if + system%tolerance = tolerance + system%max_time = max_time + end subroutine create_mopac_system + + ! deallocate memory in mopac_system + module subroutine destroy_mopac_system(system) bind(c) + type(mopac_system), intent(in) :: system + integer(c_int), pointer :: iptr(:) + real(c_double), pointer :: rptr(:) + if (system%natom > 0) then + call c_f_pointer(system%atom, iptr) + deallocate(iptr) + call c_f_pointer(system%coord, rptr) + deallocate(rptr) + end if + if (system%nlattice > 0) then + call c_f_pointer(system%lattice, rptr) + deallocate(rptr) + end if + end subroutine destroy_mopac_system + + ! deallocate memory in mopac_properties + module subroutine destroy_mopac_properties(properties) bind(c) + type(mopac_properties), intent(in) :: properties + integer :: i + integer(c_int), pointer :: iptr(:) + real(c_double), pointer :: rptr(:) + character(kind=c_char), pointer :: cptr(:) + type(c_ptr), pointer :: pptr(:) + if (.not. c_associated(properties%charge)) then + call c_f_pointer(properties%charge, rptr) + deallocate(rptr) + end if + if (.not. c_associated(properties%coord_update)) then + call c_f_pointer(properties%coord_update, rptr) + deallocate(rptr) + end if + if (.not. c_associated(properties%coord_deriv)) then + call c_f_pointer(properties%coord_deriv, rptr) + deallocate(rptr) + end if + if (.not. c_associated(properties%freq)) then + call c_f_pointer(properties%freq, rptr) + deallocate(rptr) + end if + if (.not. c_associated(properties%disp)) then + call c_f_pointer(properties%disp, rptr) + deallocate(rptr) + end if + if (.not. c_associated(properties%bond_index)) then + call c_f_pointer(properties%bond_index, iptr) + deallocate(iptr) + end if + if (.not. c_associated(properties%bond_atom)) then + call c_f_pointer(properties%bond_atom, iptr) + deallocate(iptr) + end if + if (.not. c_associated(properties%bond_order)) then + call c_f_pointer(properties%bond_order, rptr) + deallocate(rptr) + end if + if (.not. c_associated(properties%lattice_update)) then + call c_f_pointer(properties%lattice_update, rptr) + deallocate(rptr) + end if + if (.not. c_associated(properties%lattice_deriv)) then + call c_f_pointer(properties%lattice_deriv, rptr) + deallocate(rptr) + end if + if (properties%nerror > 0) then + call c_f_pointer(properties%error_msg, pptr) + do i=1, properties%nerror + call c_f_pointer(pptr(i), cptr) + deallocate(cptr) + end do + deallocate(pptr) + end if + end subroutine destroy_mopac_properties + + ! deallocate memory in mopac_state + module subroutine destroy_mopac_state(state) bind(c) + type(mopac_state), intent(in) :: state + real(c_double), pointer :: rptr(:) + if (state%mpack /= 0) then + call c_f_pointer(state%pa, rptr) + deallocate(rptr) + call c_f_pointer(state%pb, rptr) + deallocate(rptr) + end if + end subroutine destroy_mopac_state + + ! deallocate memory in mozyme_state + module subroutine destroy_mozyme_state(state) bind(c) + type(mozyme_state), intent(in) :: state + integer(c_int), pointer :: iptr(:) + real(c_double), pointer :: rptr(:) + if (state%numat /= 0) then + call c_f_pointer(state%nbonds, iptr) + deallocate(iptr) + call c_f_pointer(state%ibonds, iptr) + deallocate(iptr) + call c_f_pointer(state%iorbs, iptr) + deallocate(iptr) + call c_f_pointer(state%ncf, iptr) + deallocate(iptr) + call c_f_pointer(state%nce, iptr) + deallocate(iptr) + call c_f_pointer(state%icocc, iptr) + deallocate(iptr) + call c_f_pointer(state%icvir, iptr) + deallocate(iptr) + call c_f_pointer(state%cocc, rptr) + deallocate(rptr) + call c_f_pointer(state%cvir, rptr) + deallocate(rptr) + end if + end subroutine destroy_mozyme_state + +end submodule mopac_api_createdestroy From 4bf8f5cc70b135bcc00a46d728208e2fa8f93567 Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Thu, 24 Oct 2024 16:18:47 -0400 Subject: [PATCH 05/30] Create a Fortran header file --- include/mopac_api.H90 | 250 ++++++++++++++++++++++++++++++++++++++++++ include/mopac_api.h | 11 +- 2 files changed, 256 insertions(+), 5 deletions(-) create mode 100644 include/mopac_api.H90 diff --git a/include/mopac_api.H90 b/include/mopac_api.H90 new file mode 100644 index 00000000..2ff1264a --- /dev/null +++ b/include/mopac_api.H90 @@ -0,0 +1,250 @@ +! 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 + +! data that defines the atomistic system and MOPAC job options +type, bind(c) :: mopac_system + use iso_c_binding, only: c_int, c_double, c_ptr + implicit none + ! 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 + use iso_c_binding, only: c_int, c_double, c_ptr + implicit none + ! 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 + use iso_c_binding, only: c_int, c_ptr + implicit none + ! 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 + ! > alpha density matrix + type(c_ptr) :: pa ! real(c_double)[mpack] + ! > beta density matrix + type(c_ptr) :: pb ! real(c_double)[mpack] +end type + +! data that describes the electronic state using localized molecular orbitals +type, bind(c) :: mozyme_state + use iso_c_binding, only: c_int, c_ptr + implicit none + ! 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 + +interface + + ! allocate memory & initialize mopac_system + subroutine create_mopac_system(natom, natom_move, charge, spin, model, & + epsilon, atom, coord, nlattice, nlattice_move, & + pressure, lattice, tolerance, max_time, & + system) bind(c) + use iso_c_binding, only: c_int, c_double + implicit none + integer(c_int), intent(in) :: natom + integer(c_int), intent(in) :: natom_move + integer(c_int), intent(in) :: charge + integer(c_int), intent(in) :: spin + integer(c_int), intent(in) :: model + real(c_double), intent(in) :: epsilon + integer(c_int), dimension(natom), intent(in) :: atom + real(c_double), dimension(3*natom), intent(in) :: coord + integer(c_int), intent(in) :: nlattice + integer(c_int), intent(in) :: nlattice_move + real(c_double), intent(in) :: pressure + real(c_double), dimension(3*nlattice), intent(in) :: lattice + real(c_double), intent(in) :: tolerance + integer(c_int), intent(in) :: max_time + type(mopac_system), intent(out) :: system + end subroutine create_mopac_system + + ! deallocate memory in mopac_system + subroutine destroy_mopac_system(system) bind(c) + implicit none + type(mopac_system), intent(in) :: system + end subroutine destroy_mopac_system + + ! deallocate memory in mopac_properties (associated system is needed for size info) + subroutine destroy_mopac_properties(system, properties) bind(c) + implict none + type(mopac_system), intent(in) :: system + type(mopac_properties), intent(in) :: properties + end subroutine destroy_mopac_properties + + ! deallocate memory in mopac_state + subroutine destroy_mopac_state(state) bind(c) + implicit none + type(mopac_state), intent(in) :: state + end subroutine destroy_mopac_state + + ! deallocate memory in mozyme_state + subroutine destroy_mozyme_state(state) bind(c) + implicit none + type(mozyme_state), intent(in) :: state + end subroutine destroy_mozyme_state + + ! MOPAC electronic ground state calculation + subroutine mopac_scf(system, state, properties) bind(c) + implicit none + 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) + implicit none + 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) + implicit none + 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) + implicit none + 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) + implicit none + 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) + implicit none + 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 + subroutine run_mopac_from_input(path_to_file) bind(c) + use iso_c_binding, only: c_char + implicit none + character(kind=c_char,len=*), intent(in) :: path_to_file + end subroutine run_mopac_from_input + +end interface diff --git a/include/mopac_api.h b/include/mopac_api.h index c55ee886..68bb7bc9 100644 --- a/include/mopac_api.h +++ b/include/mopac_api.h @@ -67,15 +67,15 @@ struct mopac_properties { 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 + /* 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]-1] */ + int *bond_atom; /* [bond_index[natom]] */ /* > bond order of atoms bonded to each atom in CSC format */ - double *bond_order; /* [bond_index[natom]-1] */ + 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) */ @@ -145,8 +145,9 @@ void create_mopac_system(int *natom, int *natom_move, int *charge, int *spin, in /* deallocate memory in mopac_system */ void destroy_mopac_system(struct mopac_system *system); -/* deallocate memory in mopac_properties */ -void destroy_mopac_properties(struct mopac_properties *properties); +/* deallocate memory in mopac_properties (associated system is needed for size info) */ +void destroy_mopac_properties(struct mopac_system *system, + struct mopac_properties *properties); /* deallocate memory in mopac_state */ void destroy_mopac_state(struct mopac_state *state); From 3440aecb142a98addac36736dd8f452b857415fa Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Thu, 24 Oct 2024 16:19:43 -0400 Subject: [PATCH 06/30] Switch to memory management for C binding --- src/interface/mopac_api_createdestroy.F90 | 105 +++++++------- src/interface/mopac_api_finalize.F90 | 122 +++++++++------- src/interface/mopac_api_initialize.F90 | 33 +++-- src/interface/mopac_api_saveload.F90 | 167 +++++++++++++++------- 4 files changed, 263 insertions(+), 164 deletions(-) diff --git a/src/interface/mopac_api_createdestroy.F90 b/src/interface/mopac_api_createdestroy.F90 index 51a18865..cc8484bb 100644 --- a/src/interface/mopac_api_createdestroy.F90 +++ b/src/interface/mopac_api_createdestroy.F90 @@ -16,7 +16,7 @@ ! Diskless/stateless Application Programming Interface (API) to core MOPAC operations submodule (mopac_api) mopac_api_createdestroy - use iso_c_binding, only: c_int, c_double, c_char, c_ptr, c_loc, c_f_pointer, c_associated + use iso_c_binding, only: c_int, c_double, c_char, c_ptr, c_loc, c_f_pointer, c_null_char, c_associated implicit none contains @@ -41,6 +41,7 @@ module subroutine create_mopac_system(natom, natom_move, charge, spin, model, & real(c_double), intent(in) :: tolerance integer(c_int), intent(in) :: max_time type(mopac_system), intent(out) :: system + integer :: status integer(c_int), pointer :: iptr(:) real(c_double), pointer :: rptr(:) system%natom = natom @@ -50,10 +51,18 @@ module subroutine create_mopac_system(natom, natom_move, charge, spin, model, & system%model = model system%epsilon = epsilon if (natom > 0) then - allocate(iptr(natom)) + allocate(iptr(natom), stat=status) + if (status /= 0) then + write(*,*) "ERROR: Failed to allocate memory in CREATE_MOPAC_SYSTEM" + stop 1 + end if iptr(1:natom) = atom(1:natom) system%atom = call c_loc(iptr) - allocate(rptr(3*natom)) + allocate(rptr(3*natom), stat=status) + if (status /= 0) then + write(*,*) "ERROR: Failed to allocate memory in CREATE_MOPAC_SYSTEM" + stop 1 + end if rptr(1:3*natom) = coord(1:3*natom) system%coord = call c_loc(rptr) end if @@ -61,7 +70,11 @@ module subroutine create_mopac_system(natom, natom_move, charge, spin, model, & system%nlattice_move = nlattice_move system%pressure = pressure if (nlattice > 0) then - allocate(rptr(3*nlattice)) + allocate(rptr(3*nlattice), stat=status) + if (status /= 0) then + write(*,*) "ERROR: Failed to allocate memory in CREATE_MOPAC_SYSTEM" + stop 1 + end if rptr(1:3*nlattice) = lattice(1:3*nlattice) system%lattice = call c_loc(rptr) end if @@ -75,69 +88,65 @@ module subroutine destroy_mopac_system(system) bind(c) integer(c_int), pointer :: iptr(:) real(c_double), pointer :: rptr(:) if (system%natom > 0) then - call c_f_pointer(system%atom, iptr) + call c_f_pointer(system%atom, iptr, [system%natom]) deallocate(iptr) - call c_f_pointer(system%coord, rptr) + call c_f_pointer(system%coord, rptr, [3*system%natom]) deallocate(rptr) end if if (system%nlattice > 0) then - call c_f_pointer(system%lattice, rptr) + call c_f_pointer(system%lattice, rptr, [3*system%nlattice]) deallocate(rptr) end if end subroutine destroy_mopac_system ! deallocate memory in mopac_properties - module subroutine destroy_mopac_properties(properties) bind(c) + module subroutine destroy_mopac_properties(system, properties) bind(c) + type(mopac_system), intent(in) :: system type(mopac_properties), intent(in) :: properties - integer :: i + integer :: i, j integer(c_int), pointer :: iptr(:) real(c_double), pointer :: rptr(:) character(kind=c_char), pointer :: cptr(:) type(c_ptr), pointer :: pptr(:) - if (.not. c_associated(properties%charge)) then - call c_f_pointer(properties%charge, rptr) + if (system%natom > 0) then + call c_f_pointer(properties%charge, rptr, [system%natom]) deallocate(rptr) - end if - if (.not. c_associated(properties%coord_update)) then - call c_f_pointer(properties%coord_update, rptr) + call c_f_pointer(properties%bond_index, iptr, [system%natom+1]) + i = iptr(system%natom+1) + deallocate(iptr) + call c_f_pointer(properties%bond_atom, iptr, [i]) + deallocate(iptr) + call c_f_pointer(properties%bond_order, rptr, [i]) deallocate(rptr) end if - if (.not. c_associated(properties%coord_deriv)) then - call c_f_pointer(properties%coord_deriv, rptr) + if (system%natom_move > 0) then + call c_f_pointer(properties%coord_update, rptr, [3*system%natom_move]) deallocate(rptr) - end if - if (.not. c_associated(properties%freq)) then - call c_f_pointer(properties%freq, rptr) + call c_f_pointer(properties%coord_deriv, rptr, [3*system%natom_move]) deallocate(rptr) end if - if (.not. c_associated(properties%disp)) then - call c_f_pointer(properties%disp, rptr) + if (c_associated(properties%freq)) then + call c_f_pointer(properties%freq, rptr, [3*system%natom_move]) deallocate(rptr) end if - if (.not. c_associated(properties%bond_index)) then - call c_f_pointer(properties%bond_index, iptr) - deallocate(iptr) - end if - if (.not. c_associated(properties%bond_atom)) then - call c_f_pointer(properties%bond_atom, iptr) - deallocate(iptr) - end if - if (.not. c_associated(properties%bond_order)) then - call c_f_pointer(properties%bond_order, rptr) + if (c_associated(properties%disp)) then + call c_f_pointer(properties%disp, rptr, [3*system%natom_move,3*system%natom_move]) deallocate(rptr) end if - if (.not. c_associated(properties%lattice_update)) then - call c_f_pointer(properties%lattice_update, rptr) + if (system%nlattice_move > 0) then + call c_f_pointer(properties%lattice_update, rptr, [3*system%nlattice_move]) deallocate(rptr) - end if - if (.not. c_associated(properties%lattice_deriv)) then - call c_f_pointer(properties%lattice_deriv, rptr) + call c_f_pointer(properties%lattice_deriv, rptr, [3*system%nlattice_move]) deallocate(rptr) end if if (properties%nerror > 0) then - call c_f_pointer(properties%error_msg, pptr) + call c_f_pointer(properties%error_msg, pptr, [properties%nerror]) do i=1, properties%nerror call c_f_pointer(pptr(i), cptr) + do j=1, 120 ! using the hard-coded max size of MOPAC's error messages + if (cptr(j) == c_null_char) exit + end do + call c_f_pointer(pptr(i), cptr, [j]) deallocate(cptr) end do deallocate(pptr) @@ -149,9 +158,9 @@ module subroutine destroy_mopac_state(state) bind(c) type(mopac_state), intent(in) :: state real(c_double), pointer :: rptr(:) if (state%mpack /= 0) then - call c_f_pointer(state%pa, rptr) + call c_f_pointer(state%pa, rptr, [state%mpack]) deallocate(rptr) - call c_f_pointer(state%pb, rptr) + call c_f_pointer(state%pb, rptr, [state%mpack]) deallocate(rptr) end if end subroutine destroy_mopac_state @@ -162,23 +171,23 @@ module subroutine destroy_mozyme_state(state) bind(c) integer(c_int), pointer :: iptr(:) real(c_double), pointer :: rptr(:) if (state%numat /= 0) then - call c_f_pointer(state%nbonds, iptr) + call c_f_pointer(state%nbonds, iptr, [state%numat]) deallocate(iptr) - call c_f_pointer(state%ibonds, iptr) + call c_f_pointer(state%ibonds, iptr, [9,state%numat]) deallocate(iptr) - call c_f_pointer(state%iorbs, iptr) + call c_f_pointer(state%iorbs, iptr, [state%numat]) deallocate(iptr) - call c_f_pointer(state%ncf, iptr) + call c_f_pointer(state%ncf, iptr, [state%noccupied]) deallocate(iptr) - call c_f_pointer(state%nce, iptr) + call c_f_pointer(state%nce, iptr, [state%nvirtual]) deallocate(iptr) - call c_f_pointer(state%icocc, iptr) + call c_f_pointer(state%icocc, iptr, [state%icocc_dim]) deallocate(iptr) - call c_f_pointer(state%icvir, iptr) + call c_f_pointer(state%icvir, iptr, [state%icvir_dim]) deallocate(iptr) - call c_f_pointer(state%cocc, rptr) + call c_f_pointer(state%cocc, rptr, [state%cocc_dim]) deallocate(rptr) - call c_f_pointer(state%cvir, rptr) + call c_f_pointer(state%cvir, rptr, [state%cvir_dim]) deallocate(rptr) end if end subroutine destroy_mozyme_state diff --git a/src/interface/mopac_api_finalize.F90 b/src/interface/mopac_api_finalize.F90 index 6e88aebf..fd634f51 100644 --- a/src/interface/mopac_api_finalize.F90 +++ b/src/interface/mopac_api_finalize.F90 @@ -15,6 +15,8 @@ ! along with this program. If not, see . submodule (mopac_api:mopac_api_operations) mopac_api_finalize + use iso_c_binding, only: c_int, c_double, c_char, c_ptr, c_loc, c_null_char, c_null_ptr + use chanel_C, only : iw ! file handle for main output file use Common_arrays_C, only : xparam, & ! values of coordinates undergoing optimization geo, & ! raw coordinates of atoms @@ -47,23 +49,13 @@ ! 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 + character(kind=c_int), pointer :: cptr(:) + type(c_ptr), pointer :: 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) @@ -71,13 +63,25 @@ module subroutine mopac_finalize(properties) if (moperr) then call summary("",0) properties%nerror = dummy - allocate(properties%error_msg(properties%nerror), stat=status) + allocate(pptr(properties%nerror), stat=status) + properties%error_msg = c_loc(pptr) if (status /= 0) then properties%nerror = -properties%nerror + return else do i=1, properties%nerror call summary("",-i) - properties%error_msg(i) = trim(errtxt) + size = len_trim(errtxt) + allocate(cptr(size+1), stat=status) + pptr(i) = c_loc(cptr) + if (status /= 0) then + properties%nerror = -properties%nerror + return + end if + do j=1, size + cptr(j) = errtxt(j) + end do + cptr(size+1) = c_null_char end do end if call summary("",-abs(properties%nerror)-1) @@ -98,6 +102,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 :: rptr(:), bond_order(:) ! trigger charge & dipole calculation call chrge (p, q) @@ -115,12 +121,13 @@ subroutine mopac_record(properties) end if ! save basic properties properties%heat = escf - allocate(properties%charge(numat), stat=status) + allocate(rptr(numat), stat=status) if (status /= 0) then call mopend("Failed to allocate memory in MOPAC_FINALIZE") return end if - properties%charge = q(:numat) + rptr = q(:numat) + properties%charge = c_loc(rptr) properties%stress = voigt natom_move = nvar/3 nlattice_move = 0 @@ -131,62 +138,69 @@ subroutine mopac_record(properties) end if end do ! save properties of moveable coordinates - allocate(properties%coord_update(3*natom_move), stat=status) + allocate(rptr(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) + rptr = xparam(:3*natom_move) + properties%coord_update = c_loc(rptr) + allocate(rptr(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) + rptr = grad(:3*natom_move) + properties%coord_deriv = c_loc(rptr) if (nlattice_move > 0) then - allocate(properties%lattice_update(3*nlattice_move), stat=status) + allocate(rptr(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) + rptr = xparam(3*natom_move+1:) + properties%lattice_update = c_loc(rptr) + allocate(rptr(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:) + rptr = grad(3*natom_move+1:) + properties%lattice_deriv = c_loc(rptr) end if ! save vibrational properties if available if (index(keywrd, " FORCETS") /= 0) then - properties%calc_vibe = .true. - allocate(properties%freq(nvar), stat=status) + allocate(rptr(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) + rptr = freq + properties%freq = c_loc(rptr) + allocate(rptr(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]) + rptr = cnorml + properties%disp = c_loc(rptr) 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) + allocate(bond_index(numat+1), stat=status) + properties%bond_index = c_loc(bond_index) if (status /= 0) then call mopend("Failed to allocate memory in MOPAC_FINALIZE") return end if if (mozyme) then ! 1st pass to populate bond_index - properties%bond_index(1) = 1 + bond_index(1) = 1 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 +217,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,7 +229,7 @@ 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 @@ -226,11 +240,13 @@ subroutine mopac_record(properties) call mopend("Failed to allocate memory in MOPAC_FINALIZE") return end if + properties%bond_atom = c_loc(bond_atom) 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_order = c_loc(bond_order) do i = 1, numat io = iorbs(i) valenc = 0.d0 @@ -248,7 +264,7 @@ subroutine mopac_record(properties) valenc = valenc + 2.d0 * p(kk) end do end if - kk = properties%bond_index(i) + kk = bond_index(i) do j = 1, numat jo = iorbs(j) if (i /= j .and. ijbo (i, j) >= 0) then @@ -259,13 +275,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 +289,53 @@ subroutine mopac_record(properties) else call bonds() ! 1st pass to populate bond_index - properties%bond_index(1) = 1 + bond_index(1) = 1 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) + allocate(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) + properties%bond_atom = c_loc(bond_atom) + allocate(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_order = c_loc(bond_order) 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) 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..f827d318 100644 --- a/src/interface/mopac_api_initialize.F90 +++ b/src/interface/mopac_api_initialize.F90 @@ -15,6 +15,8 @@ ! along with this program. If not, see . submodule (mopac_api:mopac_api_operations) mopac_api_initialize + use iso_c_binding, only: c_int, c_double, c_f_pointer, c_associated + use chanel_C, only : output_fn, & ! name of output file iw ! file handle for main output file use conref_C, only : fpcref ! values of fundamental constants (data) @@ -78,6 +80,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 +90,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) + call c_f_pointer(system%coord, coord) + call c_f_pointer(system%lattice, lattice) + use_disk = .false. ! load ios, iop, and iod data into tore tore = ios + iop + iod @@ -193,20 +204,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 +226,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_saveload.F90 b/src/interface/mopac_api_saveload.F90 index cb87131f..2606afdf 100644 --- a/src/interface/mopac_api_saveload.F90 +++ b/src/interface/mopac_api_saveload.F90 @@ -15,6 +15,8 @@ ! along with this program. If not, see . submodule (mopac_api:mopac_api_operations) mopac_api_saveload + use iso_c_binding, only: c_int, c_double, c_loc, c_f_pointer + use Common_arrays_C, only: pa, pb, nbonds, ibonds use molkst_C, only: keywrd, uhf, mpack, numat use MOZYME_C, only: iorbs, noccupied, ncf, nvirtual, nce, icocc_dim, & @@ -27,25 +29,33 @@ module subroutine mopac_save(state) type(mopac_state), intent(out) :: state integer :: status + real(c_double), pointer :: rptr(:) - state%save_state = .true. - state%mpack = mpack + if (state%mpack > 0) then + call c_f_pointer(state%pa, rptr, [state%mpack]) + deallocate(rptr) + if (uhf) then + call c_f_pointer(state%pb, rptr, [state%mpack]) + deallocate(rptr) + end if + end if - if (allocated(state%pa)) deallocate(state%pa) - allocate(state%pa(mpack), stat=status) + state%mpack = mpack + allocate(rptr(mpack), stat=status) if (status /= 0) then call mopend("Failed to allocate memory in MOPAC_SAVE") return end if - state%pa = pa + rptr = pa + state%pa = c_loc(rptr) if (uhf) then - if (allocated(state%pb)) deallocate(state%pb) - allocate(state%pb(mpack), stat=status) + allocate(rptr(mpack), stat=status) if (status /= 0) then call mopend("Failed to allocate memory in MOPAC_SAVE") return end if - state%pb = pb + rptr = pb + state%pb = c_loc(rptr) end if end subroutine mopac_save @@ -53,9 +63,13 @@ end subroutine mopac_save 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) 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 +77,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,8 +87,9 @@ 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 @@ -81,8 +97,30 @@ end subroutine mopac_load module subroutine mozyme_save(state) type(mozyme_state), intent(out) :: state integer :: status + integer(c_int), pointer :: iptr(:) + real(c_double), pointer :: rptr(:) + + if (state%numat > 0) then + call c_f_pointer(state%nbonds, iptr, [state%numat]) + deallocate(iptr) + call c_f_pointer(state%ibonds, iptr, [9,state%numat]) + deallocate(iptr) + call c_f_pointer(state%iorbs, iptr, [state%numat]) + deallocate(iptr) + call c_f_pointer(state%ncf, iptr, [state%noccupied]) + deallocate(iptr) + call c_f_pointer(state%nce, iptr, [state%nvirtual]) + deallocate(iptr) + call c_f_pointer(state%icocc, iptr, [state%icocc_dim]) + deallocate(iptr) + call c_f_pointer(state%icvir, iptr, [state%icvir_dim]) + deallocate(iptr) + call c_f_pointer(state%cocc, rptr, [state%cocc_dim]) + deallocate(rptr) + call c_f_pointer(state%cvir, rptr, [state%cvir_dim]) + deallocate(rptr) + end if - state%save_state = .true. state%numat = numat state%noccupied = noccupied state%nvirtual = nvirtual @@ -90,78 +128,92 @@ module subroutine mozyme_save(state) 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) + + allocate(iptr(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) + iptr = nbonds + state%nbonds = c_loc(iptr) + + allocate(iptr(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) + iptr = ibonds + state%ibonds = c_loc(iptr) + + allocate(iptr(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) + iptr = iorbs + state%iorbs = c_loc(iptr) + + allocate(iptr(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) + iptr = ncf + state%ncf = c_loc(iptr) + + allocate(iptr(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) + iptr = nce + state%nce = c_loc(iptr) + + allocate(iptr(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) + iptr = icocc + state%icocc = c_loc(iptr) + + allocate(iptr(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) + iptr = icvir + state%icvir = c_loc(iptr) + + allocate(rptr(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) + rptr = cocc + state%cocc = c_loc(rptr) + + allocate(rptr(cvir_dim), stat=status) if (status /= 0) then call mopend("Failed to allocate memory in MOZYME_SAVE") return 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 + rptr = cvir + state%cvir = c_loc(rptr) 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(:) + 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 +300,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, iptr, [9,numat]) + ibonds = iptr + 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 From b61aa6b544c93be6125dd9c9cb1a9c8595bc5b6b Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Thu, 24 Oct 2024 16:28:15 -0400 Subject: [PATCH 07/30] Switch bond order matrix to 0-based indexing --- src/interface/mopac_api.F90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/interface/mopac_api.F90 b/src/interface/mopac_api.F90 index 9f4c7cd2..6944926b 100644 --- a/src/interface/mopac_api.F90 +++ b/src/interface/mopac_api.F90 @@ -177,8 +177,9 @@ module subroutine destroy_mopac_system(system) bind(c) type(mopac_system), intent(in) :: system end subroutine destroy_mopac_system - ! deallocate memory in mopac_properties - module subroutine destroy_mopac_properties(properties) bind(c) + ! deallocate memory in mopac_properties (associated system is needed for size info) + module subroutine destroy_mopac_properties(system, properties) bind(c) + type(mopac_system), intent(in) :: system type(mopac_properties), intent(in) :: properties end subroutine destroy_mopac_properties From 66cb16f3dee59f537b2ba9fcebce7896d9bdd21a Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Fri, 25 Oct 2024 12:21:46 -0400 Subject: [PATCH 08/30] Clean source up for compilability --- CMakeLists.txt | 2 +- src/interface/mopac_api.F90 | 10 +++++----- src/interface/mopac_api_createdestroy.F90 | 13 ++++++------- src/interface/mopac_api_finalize.F90 | 20 +++++++++----------- src/interface/mopac_api_initialize.F90 | 8 +++----- src/interface/mopac_api_operations.F90 | 5 ++--- src/interface/mopac_api_saveload.F90 | 20 +++++++++----------- 7 files changed, 35 insertions(+), 43 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index b5e1c475..0a5201f4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -319,7 +319,7 @@ 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(PUBLIC_HEADER "${CMAKE_SOURCE_DIR}/include/mopac_api.h" COMPONENT redist EXCLUDE_FROM_ALL) +install(FILES "${CMAKE_SOURCE_DIR}/include/mopac_api.h" 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/src/interface/mopac_api.F90 b/src/interface/mopac_api.F90 index 6944926b..1baa63fb 100644 --- a/src/interface/mopac_api.F90 +++ b/src/interface/mopac_api.F90 @@ -16,7 +16,7 @@ ! Diskless/stateless Application Programming Interface (API) to core MOPAC operations module mopac_api - use iso_c_binding, only: c_int, c_double, c_char, c_ptr + use, intrinsic :: iso_c_binding implicit none private @@ -79,15 +79,15 @@ module mopac_api 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 + ! 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)-1] + 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)-1] + 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) @@ -237,7 +237,7 @@ end subroutine mozyme_vibe ! run MOPAC conventionally from an input file module subroutine run_mopac_from_input(path_to_file) bind(c) - character(kind=c_char,len=*), intent(in) :: path_to_file + character(kind=c_char), dimension(*), intent(in) :: path_to_file end subroutine run_mopac_from_input end interface diff --git a/src/interface/mopac_api_createdestroy.F90 b/src/interface/mopac_api_createdestroy.F90 index cc8484bb..9c5551ab 100644 --- a/src/interface/mopac_api_createdestroy.F90 +++ b/src/interface/mopac_api_createdestroy.F90 @@ -16,7 +16,6 @@ ! Diskless/stateless Application Programming Interface (API) to core MOPAC operations submodule (mopac_api) mopac_api_createdestroy - use iso_c_binding, only: c_int, c_double, c_char, c_ptr, c_loc, c_f_pointer, c_null_char, c_associated implicit none contains @@ -57,14 +56,14 @@ module subroutine create_mopac_system(natom, natom_move, charge, spin, model, & stop 1 end if iptr(1:natom) = atom(1:natom) - system%atom = call c_loc(iptr) + system%atom = c_loc(iptr) allocate(rptr(3*natom), stat=status) if (status /= 0) then write(*,*) "ERROR: Failed to allocate memory in CREATE_MOPAC_SYSTEM" stop 1 end if rptr(1:3*natom) = coord(1:3*natom) - system%coord = call c_loc(rptr) + system%coord = c_loc(rptr) end if system%nlattice = nlattice system%nlattice_move = nlattice_move @@ -76,7 +75,7 @@ module subroutine create_mopac_system(natom, natom_move, charge, spin, model, & stop 1 end if rptr(1:3*nlattice) = lattice(1:3*nlattice) - system%lattice = call c_loc(rptr) + system%lattice = c_loc(rptr) end if system%tolerance = tolerance system%max_time = max_time @@ -130,7 +129,7 @@ module subroutine destroy_mopac_properties(system, properties) bind(c) deallocate(rptr) end if if (c_associated(properties%disp)) then - call c_f_pointer(properties%disp, rptr, [3*system%natom_move,3*system%natom_move]) + call c_f_pointer(properties%disp, rptr, [3*system%natom_move*3*system%natom_move]) deallocate(rptr) end if if (system%nlattice_move > 0) then @@ -142,7 +141,7 @@ module subroutine destroy_mopac_properties(system, properties) bind(c) if (properties%nerror > 0) then call c_f_pointer(properties%error_msg, pptr, [properties%nerror]) do i=1, properties%nerror - call c_f_pointer(pptr(i), cptr) + call c_f_pointer(pptr(i), cptr, [120]) do j=1, 120 ! using the hard-coded max size of MOPAC's error messages if (cptr(j) == c_null_char) exit end do @@ -173,7 +172,7 @@ module subroutine destroy_mozyme_state(state) bind(c) if (state%numat /= 0) then call c_f_pointer(state%nbonds, iptr, [state%numat]) deallocate(iptr) - call c_f_pointer(state%ibonds, iptr, [9,state%numat]) + call c_f_pointer(state%ibonds, iptr, [9*state%numat]) deallocate(iptr) call c_f_pointer(state%iorbs, iptr, [state%numat]) deallocate(iptr) diff --git a/src/interface/mopac_api_finalize.F90 b/src/interface/mopac_api_finalize.F90 index fd634f51..27020184 100644 --- a/src/interface/mopac_api_finalize.F90 +++ b/src/interface/mopac_api_finalize.F90 @@ -15,8 +15,6 @@ ! along with this program. If not, see . submodule (mopac_api:mopac_api_operations) mopac_api_finalize - use iso_c_binding, only: c_int, c_double, c_char, c_ptr, c_loc, c_null_char, c_null_ptr - use chanel_C, only : iw ! file handle for main output file use Common_arrays_C, only : xparam, & ! values of coordinates undergoing optimization geo, & ! raw coordinates of atoms @@ -79,7 +77,7 @@ module subroutine mopac_finalize(properties) return end if do j=1, size - cptr(j) = errtxt(j) + cptr(j) = errtxt(j:j) end do cptr(size+1) = c_null_char end do @@ -197,7 +195,7 @@ subroutine mopac_record(properties) end if if (mozyme) then ! 1st pass to populate bond_index - bond_index(1) = 1 + bond_index(1) = 0 do i = 1, numat io = iorbs(i) bond_index(i+1) = bond_index(i) @@ -235,13 +233,13 @@ subroutine mopac_record(properties) end do end do ! 2nd pass to populate bond_atom and bond_order - allocate(properties%bond_atom(properties%bond_index(numat+1)), stat=status) + allocate(bond_atom(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 = c_loc(bond_atom) - allocate(properties%bond_order(properties%bond_index(numat+1)), stat=status) + allocate(bond_order(bond_index(numat+1)), stat=status) if (status /= 0) then call mopend("Failed to allocate memory in MOPAC_FINALIZE") return @@ -264,7 +262,7 @@ subroutine mopac_record(properties) valenc = valenc + 2.d0 * p(kk) end do end if - kk = bond_index(i) + kk = bond_index(i) + 1 do j = 1, numat jo = iorbs(j) if (i /= j .and. ijbo (i, j) >= 0) then @@ -289,7 +287,7 @@ subroutine mopac_record(properties) else call bonds() ! 1st pass to populate bond_index - bond_index(1) = 1 + bond_index(1) = 0 do i = 1, numat bond_index(i+1) = bond_index(i) ku = i*(i-1)/2 + 1 @@ -308,13 +306,13 @@ subroutine mopac_record(properties) end do end do ! 2nd pass to populate bond_atom and bond_order - allocate(bond_atom(properties%bond_index(numat+1)), stat=status) + allocate(bond_atom(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 = c_loc(bond_atom) - allocate(bond_order(properties%bond_index(numat+1)), stat=status) + allocate(bond_order(bond_index(numat+1)), stat=status) if (status /= 0) then call mopend("Failed to allocate memory in MOPAC_FINALIZE") return @@ -323,7 +321,7 @@ subroutine mopac_record(properties) do i = 1, numat ku = i*(i-1)/2 + 1 kl = (i+1)*(i+2)/2 - 1 - kk = bond_index(i) + kk = bond_index(i) + 1 do j = 1, i if (bondab(ku) > 0.01d0) then bond_atom(kk) = j diff --git a/src/interface/mopac_api_initialize.F90 b/src/interface/mopac_api_initialize.F90 index f827d318..24c0618e 100644 --- a/src/interface/mopac_api_initialize.F90 +++ b/src/interface/mopac_api_initialize.F90 @@ -15,8 +15,6 @@ ! along with this program. If not, see . submodule (mopac_api:mopac_api_operations) mopac_api_initialize - use iso_c_binding, only: c_int, c_double, c_f_pointer, c_associated - use chanel_C, only : output_fn, & ! name of output file iw ! file handle for main output file use conref_C, only : fpcref ! values of fundamental constants (data) @@ -106,9 +104,9 @@ module subroutine mopac_initialize(system) if (moperr) return ! convert C pointers - call c_f_pointer(system%atom, atom) - call c_f_pointer(system%coord, coord) - call c_f_pointer(system%lattice, lattice) + 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 diff --git a/src/interface/mopac_api_operations.F90 b/src/interface/mopac_api_operations.F90 index 55d9ffbc..7dae759e 100644 --- a/src/interface/mopac_api_operations.F90 +++ b/src/interface/mopac_api_operations.F90 @@ -15,7 +15,6 @@ ! along with this program. If not, see . submodule (mopac_api) mopac_api_operations - use iso_c_binding use Common_arrays_C, only: xparam, grad, lopt use molkst_C, only: keywrd, escf, moperr, nvar, gui, jobnam, run implicit none @@ -168,12 +167,12 @@ end subroutine mozyme_vibe ! Run MOPAC conventionally from an input file module subroutine run_mopac_from_input(path_to_file) - character(kind=c_char,len=*), intent(in) :: path_to_file + 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) = path_to_file(i) + jobnam(i:i) = path_to_file(i) i = i + 1 end do gui = .false. diff --git a/src/interface/mopac_api_saveload.F90 b/src/interface/mopac_api_saveload.F90 index 2606afdf..c32fc43e 100644 --- a/src/interface/mopac_api_saveload.F90 +++ b/src/interface/mopac_api_saveload.F90 @@ -15,8 +15,6 @@ ! along with this program. If not, see . submodule (mopac_api:mopac_api_operations) mopac_api_saveload - use iso_c_binding, only: c_int, c_double, c_loc, c_f_pointer - use Common_arrays_C, only: pa, pb, nbonds, ibonds use molkst_C, only: keywrd, uhf, mpack, numat use MOZYME_C, only: iorbs, noccupied, ncf, nvirtual, nce, icocc_dim, & @@ -97,14 +95,14 @@ end subroutine mopac_load module subroutine mozyme_save(state) type(mozyme_state), intent(out) :: state integer :: status - integer(c_int), pointer :: iptr(:) + integer(c_int), pointer :: iptr(:), iptr2(:,:) real(c_double), pointer :: rptr(:) if (state%numat > 0) then call c_f_pointer(state%nbonds, iptr, [state%numat]) deallocate(iptr) - call c_f_pointer(state%ibonds, iptr, [9,state%numat]) - deallocate(iptr) + call c_f_pointer(state%ibonds, iptr2, [9,state%numat]) + deallocate(iptr2) call c_f_pointer(state%iorbs, iptr, [state%numat]) deallocate(iptr) call c_f_pointer(state%ncf, iptr, [state%noccupied]) @@ -137,13 +135,13 @@ module subroutine mozyme_save(state) iptr = nbonds state%nbonds = c_loc(iptr) - allocate(iptr(9,numat), stat=status) + allocate(iptr2(9,numat), stat=status) if (status /= 0) then call mopend("Failed to allocate memory in MOZYME_SAVE") return end if - iptr = ibonds - state%ibonds = c_loc(iptr) + iptr2 = ibonds + state%ibonds = c_loc(iptr2) allocate(iptr(numat), stat=status) if (status /= 0) then @@ -206,7 +204,7 @@ end subroutine mozyme_save module subroutine mozyme_load(state) type(mozyme_state), intent(in) :: state integer :: status, i, j, k, l - integer(c_int), pointer :: iptr(:) + integer(c_int), pointer :: iptr(:), iptr2(:,:) real(c_double), pointer :: rptr(:) if(state%numat > 0) then @@ -302,8 +300,8 @@ module subroutine mozyme_load(state) end if call c_f_pointer(state%nbonds, iptr, [numat]) nbonds = iptr - call c_f_pointer(state%ibonds, iptr, [9,numat]) - ibonds = 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]) From 4cd2dbd0bd3ff5da8a62b5fb38753368292f2da5 Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Wed, 6 Nov 2024 12:36:21 -0500 Subject: [PATCH 09/30] Finish switch of API to C binding w/ Fortran wrapper --- CMakeLists.txt | 13 +- include/mopac_api.H90 | 250 ------------- include/mopac_api.h | 51 +-- src/interface/mopac_api_createdestroy.F90 | 423 ++++++++++++++-------- src/interface/mopac_api_finalize.F90 | 132 ++----- src/interface/mopac_api_operations.F90 | 14 +- src/interface/mopac_api_saveload.F90 | 138 +------ tests/CMakeLists.txt | 7 + tests/mopac_api_test.F90 | 161 ++++---- 9 files changed, 455 insertions(+), 734 deletions(-) delete mode 100644 include/mopac_api.H90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 0a5201f4..28249488 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) @@ -261,11 +266,6 @@ if(TESTS) endif() message(FATAL_ERROR "Python3 and Numpy are required for MOPAC testing (testing can be disabled with -DTESTS=OFF)") endif() - 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(tests) endif() @@ -320,6 +320,9 @@ 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/mopac_api.H90 b/include/mopac_api.H90 deleted file mode 100644 index 2ff1264a..00000000 --- a/include/mopac_api.H90 +++ /dev/null @@ -1,250 +0,0 @@ -! 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 - -! data that defines the atomistic system and MOPAC job options -type, bind(c) :: mopac_system - use iso_c_binding, only: c_int, c_double, c_ptr - implicit none - ! 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 - use iso_c_binding, only: c_int, c_double, c_ptr - implicit none - ! 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 - use iso_c_binding, only: c_int, c_ptr - implicit none - ! 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 - ! > alpha density matrix - type(c_ptr) :: pa ! real(c_double)[mpack] - ! > beta density matrix - type(c_ptr) :: pb ! real(c_double)[mpack] -end type - -! data that describes the electronic state using localized molecular orbitals -type, bind(c) :: mozyme_state - use iso_c_binding, only: c_int, c_ptr - implicit none - ! 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 - -interface - - ! allocate memory & initialize mopac_system - subroutine create_mopac_system(natom, natom_move, charge, spin, model, & - epsilon, atom, coord, nlattice, nlattice_move, & - pressure, lattice, tolerance, max_time, & - system) bind(c) - use iso_c_binding, only: c_int, c_double - implicit none - integer(c_int), intent(in) :: natom - integer(c_int), intent(in) :: natom_move - integer(c_int), intent(in) :: charge - integer(c_int), intent(in) :: spin - integer(c_int), intent(in) :: model - real(c_double), intent(in) :: epsilon - integer(c_int), dimension(natom), intent(in) :: atom - real(c_double), dimension(3*natom), intent(in) :: coord - integer(c_int), intent(in) :: nlattice - integer(c_int), intent(in) :: nlattice_move - real(c_double), intent(in) :: pressure - real(c_double), dimension(3*nlattice), intent(in) :: lattice - real(c_double), intent(in) :: tolerance - integer(c_int), intent(in) :: max_time - type(mopac_system), intent(out) :: system - end subroutine create_mopac_system - - ! deallocate memory in mopac_system - subroutine destroy_mopac_system(system) bind(c) - implicit none - type(mopac_system), intent(in) :: system - end subroutine destroy_mopac_system - - ! deallocate memory in mopac_properties (associated system is needed for size info) - subroutine destroy_mopac_properties(system, properties) bind(c) - implict none - type(mopac_system), intent(in) :: system - type(mopac_properties), intent(in) :: properties - end subroutine destroy_mopac_properties - - ! deallocate memory in mopac_state - subroutine destroy_mopac_state(state) bind(c) - implicit none - type(mopac_state), intent(in) :: state - end subroutine destroy_mopac_state - - ! deallocate memory in mozyme_state - subroutine destroy_mozyme_state(state) bind(c) - implicit none - type(mozyme_state), intent(in) :: state - end subroutine destroy_mozyme_state - - ! MOPAC electronic ground state calculation - subroutine mopac_scf(system, state, properties) bind(c) - implicit none - 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) - implicit none - 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) - implicit none - 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) - implicit none - 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) - implicit none - 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) - implicit none - 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 - subroutine run_mopac_from_input(path_to_file) bind(c) - use iso_c_binding, only: c_char - implicit none - character(kind=c_char,len=*), intent(in) :: path_to_file - end subroutine run_mopac_from_input - -end interface diff --git a/include/mopac_api.h b/include/mopac_api.h index 68bb7bc9..177e5219 100644 --- a/include/mopac_api.h +++ b/include/mopac_api.h @@ -1,18 +1,5 @@ /* 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 */ @@ -93,10 +80,12 @@ 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] */ + double *pb; /* [mpack], NULL if uhf == 0 */ }; /* data that describes the electronic state using localized molecular orbitals */ @@ -136,25 +125,6 @@ struct mozyme_state { double *cvir; /* [cvir_dim] */ }; -/* allocate memory & initialize mopac_system */ -void create_mopac_system(int *natom, int *natom_move, int *charge, int *spin, int *model, - double *epsilon, int *atom, double *coord, int *nlattice, int *nlattice_move, - double *pressure, double *lattice, double *tolerance, int *max_time, - struct mopac_system *system); - -/* deallocate memory in mopac_system */ -void destroy_mopac_system(struct mopac_system *system); - -/* deallocate memory in mopac_properties (associated system is needed for size info) */ -void destroy_mopac_properties(struct mopac_system *system, - 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); - /* MOPAC electronic ground state calculation */ void mopac_scf(struct mopac_system *system, struct mopac_state *state, @@ -185,6 +155,21 @@ 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 */ void run_mopac_from_input(char *path_to_file); diff --git a/src/interface/mopac_api_createdestroy.F90 b/src/interface/mopac_api_createdestroy.F90 index 9c5551ab..99a687d7 100644 --- a/src/interface/mopac_api_createdestroy.F90 +++ b/src/interface/mopac_api_createdestroy.F90 @@ -20,175 +20,310 @@ contains - ! allocate memory & initialize mopac_system - module subroutine create_mopac_system(natom, natom_move, charge, spin, model, & - epsilon, atom, coord, nlattice, nlattice_move, & - pressure, lattice, tolerance, max_time, & - system) bind(c) - integer(c_int), intent(in) :: natom - integer(c_int), intent(in) :: natom_move - integer(c_int), intent(in) :: charge - integer(c_int), intent(in) :: spin - integer(c_int), intent(in) :: model - real(c_double), intent(in) :: epsilon - integer(c_int), dimension(natom), intent(in) :: atom - real(c_double), dimension(3*natom), intent(in) :: coord - integer(c_int), intent(in) :: nlattice - integer(c_int), intent(in) :: nlattice_move - real(c_double), intent(in) :: pressure - real(c_double), dimension(3*nlattice), intent(in) :: lattice - real(c_double), intent(in) :: tolerance - integer(c_int), intent(in) :: max_time - type(mopac_system), intent(out) :: system - integer :: status - integer(c_int), pointer :: iptr(:) - real(c_double), pointer :: rptr(:) - system%natom = natom - system%natom_move = natom_move - system%charge = charge - system%spin = spin - system%model = model - system%epsilon = epsilon - if (natom > 0) then - allocate(iptr(natom), stat=status) - if (status /= 0) then - write(*,*) "ERROR: Failed to allocate memory in CREATE_MOPAC_SYSTEM" - stop 1 - end if - iptr(1:natom) = atom(1:natom) - system%atom = c_loc(iptr) - allocate(rptr(3*natom), stat=status) - if (status /= 0) then - write(*,*) "ERROR: Failed to allocate memory in CREATE_MOPAC_SYSTEM" - stop 1 - end if - rptr(1:3*natom) = coord(1:3*natom) - system%coord = c_loc(rptr) - end if - system%nlattice = nlattice - system%nlattice_move = nlattice_move - system%pressure = pressure - if (nlattice > 0) then - allocate(rptr(3*nlattice), stat=status) - if (status /= 0) then - write(*,*) "ERROR: Failed to allocate memory in CREATE_MOPAC_SYSTEM" - stop 1 - end if - rptr(1:3*nlattice) = lattice(1:3*nlattice) - system%lattice = c_loc(rptr) + ! allocate memory for mopac_state + module subroutine create_mopac_state(state) bind(c) + 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 - system%tolerance = tolerance - system%max_time = max_time - end subroutine create_mopac_system + end subroutine create_mopac_state - ! deallocate memory in mopac_system - module subroutine destroy_mopac_system(system) bind(c) - type(mopac_system), intent(in) :: system - integer(c_int), pointer :: iptr(:) - real(c_double), pointer :: rptr(:) - if (system%natom > 0) then - call c_f_pointer(system%atom, iptr, [system%natom]) - deallocate(iptr) - call c_f_pointer(system%coord, rptr, [3*system%natom]) - deallocate(rptr) - end if - if (system%nlattice > 0) then - call c_f_pointer(system%lattice, rptr, [3*system%nlattice]) - deallocate(rptr) - end if - end subroutine destroy_mopac_system + ! allocate memory for mozyme_state + module subroutine create_mozyme_state(state) bind(c) + 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(system, properties) bind(c) - type(mopac_system), intent(in) :: system + module subroutine destroy_mopac_properties(properties) bind(c) type(mopac_properties), intent(in) :: properties - integer :: i, j - integer(c_int), pointer :: iptr(:) - real(c_double), pointer :: rptr(:) - character(kind=c_char), pointer :: cptr(:) + integer :: i type(c_ptr), pointer :: pptr(:) - if (system%natom > 0) then - call c_f_pointer(properties%charge, rptr, [system%natom]) - deallocate(rptr) - call c_f_pointer(properties%bond_index, iptr, [system%natom+1]) - i = iptr(system%natom+1) - deallocate(iptr) - call c_f_pointer(properties%bond_atom, iptr, [i]) - deallocate(iptr) - call c_f_pointer(properties%bond_order, rptr, [i]) - deallocate(rptr) - end if - if (system%natom_move > 0) then - call c_f_pointer(properties%coord_update, rptr, [3*system%natom_move]) - deallocate(rptr) - call c_f_pointer(properties%coord_deriv, rptr, [3*system%natom_move]) - deallocate(rptr) - end if - if (c_associated(properties%freq)) then - call c_f_pointer(properties%freq, rptr, [3*system%natom_move]) - deallocate(rptr) - end if - if (c_associated(properties%disp)) then - call c_f_pointer(properties%disp, rptr, [3*system%natom_move*3*system%natom_move]) - deallocate(rptr) - end if - if (system%nlattice_move > 0) then - call c_f_pointer(properties%lattice_update, rptr, [3*system%nlattice_move]) - deallocate(rptr) - call c_f_pointer(properties%lattice_deriv, rptr, [3*system%nlattice_move]) - deallocate(rptr) - end if + 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 c_f_pointer(pptr(i), cptr, [120]) - do j=1, 120 ! using the hard-coded max size of MOPAC's error messages - if (cptr(j) == c_null_char) exit - end do - call c_f_pointer(pptr(i), cptr, [j]) - deallocate(cptr) + call destroy_char(pptr(i)) end do - deallocate(pptr) + 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) type(mopac_state), intent(in) :: state - real(c_double), pointer :: rptr(:) if (state%mpack /= 0) then - call c_f_pointer(state%pa, rptr, [state%mpack]) - deallocate(rptr) - call c_f_pointer(state%pb, rptr, [state%mpack]) - deallocate(rptr) + 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) type(mozyme_state), intent(in) :: state - integer(c_int), pointer :: iptr(:) - real(c_double), pointer :: rptr(:) if (state%numat /= 0) then - call c_f_pointer(state%nbonds, iptr, [state%numat]) - deallocate(iptr) - call c_f_pointer(state%ibonds, iptr, [9*state%numat]) - deallocate(iptr) - call c_f_pointer(state%iorbs, iptr, [state%numat]) - deallocate(iptr) - call c_f_pointer(state%ncf, iptr, [state%noccupied]) - deallocate(iptr) - call c_f_pointer(state%nce, iptr, [state%nvirtual]) - deallocate(iptr) - call c_f_pointer(state%icocc, iptr, [state%icocc_dim]) - deallocate(iptr) - call c_f_pointer(state%icvir, iptr, [state%icvir_dim]) - deallocate(iptr) - call c_f_pointer(state%cocc, rptr, [state%cocc_dim]) - deallocate(rptr) - call c_f_pointer(state%cvir, rptr, [state%cvir_dim]) - deallocate(rptr) + 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 + integer(c_int), pointer :: ptr(:) +#ifndef MOPAC_API_MALLOC + 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 + create_int = malloc(c_sizeof(ptr)) + call c_f_pointer(create_int, ptr, [size]) +#endif + end function create_int + module function create_int2(size, size2) + integer, intent(in) :: size + integer, intent(in) :: size2 + type(c_ptr) :: create_int2 + integer(c_int), pointer :: ptr(:,:) +#ifndef MOPAC_API_MALLOC + 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 + create_int2 = malloc(c_sizeof(ptr)) + call c_f_pointer(create_int2, ptr, [size,size2]) +#endif + end function create_int2 + module function create_real(size) + integer, intent(in) :: size + type(c_ptr) :: create_real + real(c_double), pointer :: ptr(:) +#ifndef MOPAC_API_MALLOC + 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 + create_real = malloc(c_sizeof(ptr)) + call c_f_pointer(create_real, ptr, [size]) +#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 + create_copy_int = malloc(c_sizeof(ptr)) + call c_f_pointer(create_copy_int, ptr, size) +#endif + ptr = array + 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 + create_copy_int2 = malloc(c_sizeof(ptr)) + call c_f_pointer(create_copy_int2, ptr, size) +#endif + ptr = array + 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 + create_copy_real = malloc(c_sizeof(ptr)) + call c_f_pointer(create_copy_real, ptr, size) +#endif + ptr = array + 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 + create_copy_char = malloc(c_sizeof(ptr)) + 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 + create_copy_ptr = malloc(c_sizeof(ptr)) + call c_f_pointer(create_copy_ptr, ptr, size) +#endif + ptr = array + 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 + call free(copy) +#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 + call free(copy) +#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 + call free(copy) +#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 + call free(copy) +#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 27020184..61e5a538 100644 --- a/src/interface/mopac_api_finalize.F90 +++ b/src/interface/mopac_api_finalize.F90 @@ -48,8 +48,7 @@ module subroutine mopac_finalize(properties) type(mopac_properties), intent(out) :: properties integer :: status, i, j, size - character(kind=c_int), pointer :: cptr(:) - type(c_ptr), pointer :: pptr(:) + type(c_ptr), allocatable :: pptr(:) ! close dummy output file to free up /dev/null close(iw) @@ -57,31 +56,33 @@ module subroutine mopac_finalize(properties) ! 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(pptr(properties%nerror), stat=status) - properties%error_msg = c_loc(pptr) if (status /= 0) then properties%nerror = -properties%nerror + properties%error_msg = c_null_ptr return else do i=1, properties%nerror call summary("",-i) - size = len_trim(errtxt) - allocate(cptr(size+1), stat=status) - pptr(i) = c_loc(cptr) - if (status /= 0) then - properties%nerror = -properties%nerror - return - end if - do j=1, size - cptr(j) = errtxt(j:j) - end do - cptr(size+1) = c_null_char + 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 @@ -101,7 +102,7 @@ subroutine mopac_record(properties) 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 :: rptr(:), bond_order(:) + real(c_double), pointer :: bond_order(:) ! trigger charge & dipole calculation call chrge (p, q) @@ -119,13 +120,7 @@ subroutine mopac_record(properties) end if ! save basic properties properties%heat = escf - allocate(rptr(numat), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOPAC_FINALIZE") - return - end if - rptr = q(:numat) - properties%charge = c_loc(rptr) + properties%charge = create_copy(q, [numat]) properties%stress = voigt natom_move = nvar/3 nlattice_move = 0 @@ -136,63 +131,26 @@ subroutine mopac_record(properties) end if end do ! save properties of moveable coordinates - allocate(rptr(3*natom_move), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOPAC_FINALIZE") - return - end if - rptr = xparam(:3*natom_move) - properties%coord_update = c_loc(rptr) - allocate(rptr(3*natom_move), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOPAC_FINALIZE") - return - end if - rptr = grad(:3*natom_move) - properties%coord_deriv = c_loc(rptr) + properties%coord_update = create_copy(xparam, [3*natom_move]) + properties%coord_deriv = create_copy(grad, [3*natom_move]) if (nlattice_move > 0) then - allocate(rptr(3*nlattice_move), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOPAC_FINALIZE") - return - end if - rptr = xparam(3*natom_move+1:) - properties%lattice_update = c_loc(rptr) - allocate(rptr(3*nlattice_move), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOPAC_FINALIZE") - return - end if - rptr = grad(3*natom_move+1:) - properties%lattice_deriv = c_loc(rptr) + 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 - allocate(rptr(nvar), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOPAC_FINALIZE") - return - end if - rptr = freq - properties%freq = c_loc(rptr) - allocate(rptr(nvar*nvar), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOPAC_FINALIZE") - return - end if - rptr = cnorml - properties%disp = c_loc(rptr) + properties%freq = create_copy(freq, [nvar]) + properties%disp = create_copy(cnorml, [nvar*nvar]) else properties%freq = c_null_ptr properties%disp = c_null_ptr end if ! prune bond order matrix - allocate(bond_index(numat+1), stat=status) - properties%bond_index = c_loc(bond_index) - 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 bond_index(1) = 0 @@ -233,18 +191,10 @@ subroutine mopac_record(properties) end do end do ! 2nd pass to populate bond_atom and bond_order - allocate(bond_atom(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 = c_loc(bond_atom) - allocate(bond_order(bond_index(numat+1)), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOPAC_FINALIZE") - return - end if - properties%bond_order = c_loc(bond_order) + 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 @@ -306,18 +256,10 @@ subroutine mopac_record(properties) end do end do ! 2nd pass to populate bond_atom and bond_order - allocate(bond_atom(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 = c_loc(bond_atom) - allocate(bond_order(bond_index(numat+1)), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOPAC_FINALIZE") - return - end if - properties%bond_order = c_loc(bond_order) + 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 diff --git a/src/interface/mopac_api_operations.F90 b/src/interface/mopac_api_operations.F90 index 7dae759e..27a5844c 100644 --- a/src/interface/mopac_api_operations.F90 +++ b/src/interface/mopac_api_operations.F90 @@ -56,7 +56,7 @@ end subroutine mozyme_load contains ! MOPAC electronic ground state calculation - module subroutine mopac_scf(system, state, properties) + module subroutine mopac_scf(system, state, properties) bind(c) type(mopac_system), intent(in) :: system type(mopac_state), intent(inout) :: state type(mopac_properties), intent(out) :: properties @@ -71,7 +71,7 @@ module subroutine mopac_scf(system, state, properties) end subroutine mopac_scf ! MOPAC geometry relaxation - module subroutine mopac_relax(system, state, properties) + module subroutine mopac_relax(system, state, properties) bind(c) type(mopac_system), intent(in) :: system type(mopac_state), intent(inout) :: state type(mopac_properties), intent(out) :: properties @@ -86,7 +86,7 @@ module subroutine mopac_relax(system, state, properties) end subroutine mopac_relax ! MOPAC vibrational calculation - module subroutine mopac_vibe(system, state, properties) + module subroutine mopac_vibe(system, state, properties) bind(c) type(mopac_system), intent(in) :: system type(mopac_state), intent(inout) :: state type(mopac_properties), intent(out) :: properties @@ -111,7 +111,7 @@ module subroutine mopac_vibe(system, state, properties) end subroutine mopac_vibe ! MOZYME electronic ground state calculation - module subroutine mozyme_scf(system, state, properties) + module subroutine mozyme_scf(system, state, properties) bind(c) type(mopac_system), intent(in) :: system type(mozyme_state), intent(inout) :: state type(mopac_properties), intent(out) :: properties @@ -126,7 +126,7 @@ module subroutine mozyme_scf(system, state, properties) end subroutine mozyme_scf ! MOZYME geometry relaxation - module subroutine mozyme_relax(system, state, properties) + module subroutine mozyme_relax(system, state, properties) bind(c) type(mopac_system), intent(in) :: system type(mozyme_state), intent(inout) :: state type(mopac_properties), intent(out) :: properties @@ -141,7 +141,7 @@ module subroutine mozyme_relax(system, state, properties) end subroutine mozyme_relax ! MOZYME vibrational calculation - module subroutine mozyme_vibe(system, state, properties) + module subroutine mozyme_vibe(system, state, properties) bind(c) type(mopac_system), intent(in) :: system type(mozyme_state), intent(inout) :: state type(mopac_properties), intent(out) :: properties @@ -166,7 +166,7 @@ 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) + module subroutine run_mopac_from_input(path_to_file) bind(c) character(kind=c_char), dimension(*), intent(in) :: path_to_file integer :: i i = 1 diff --git a/src/interface/mopac_api_saveload.F90 b/src/interface/mopac_api_saveload.F90 index c32fc43e..4817d1af 100644 --- a/src/interface/mopac_api_saveload.F90 +++ b/src/interface/mopac_api_saveload.F90 @@ -26,35 +26,19 @@ ! save MOPAC density matrices module subroutine mopac_save(state) type(mopac_state), intent(out) :: state - integer :: status - real(c_double), pointer :: rptr(:) - if (state%mpack > 0) then - call c_f_pointer(state%pa, rptr, [state%mpack]) - deallocate(rptr) - if (uhf) then - call c_f_pointer(state%pb, rptr, [state%mpack]) - deallocate(rptr) - end if - end if + if (state%mpack > 0) call destroy_mopac_state(state) state%mpack = mpack - allocate(rptr(mpack), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOPAC_SAVE") - return - end if - rptr = pa - state%pa = c_loc(rptr) + state%pa = create_copy(pa, [mpack]) if (uhf) then - allocate(rptr(mpack), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOPAC_SAVE") - return - end if - rptr = pb - state%pb = c_loc(rptr) + state%uhf = 1 + state%pb = create_copy(pb, [mpack]) + else + state%uhf = 0 + state%pb = c_null_ptr end if + end subroutine mopac_save ! load MOPAC density matrices, or construct initial guesses @@ -64,7 +48,7 @@ module subroutine mopac_load(state) real(c_double), pointer :: rptr(:) if(state%mpack > 0) then - if (state%mpack /= mpack) then + if (state%mpack /= mpack .or. (state%uhf /= 0 .neqv. uhf)) then call mopend("Attempting to load incompatible MOPAC state") return end if @@ -94,30 +78,8 @@ end subroutine mopac_load ! save MOZYME density matrix module subroutine mozyme_save(state) type(mozyme_state), intent(out) :: state - integer :: status - integer(c_int), pointer :: iptr(:), iptr2(:,:) - real(c_double), pointer :: rptr(:) - if (state%numat > 0) then - call c_f_pointer(state%nbonds, iptr, [state%numat]) - deallocate(iptr) - call c_f_pointer(state%ibonds, iptr2, [9,state%numat]) - deallocate(iptr2) - call c_f_pointer(state%iorbs, iptr, [state%numat]) - deallocate(iptr) - call c_f_pointer(state%ncf, iptr, [state%noccupied]) - deallocate(iptr) - call c_f_pointer(state%nce, iptr, [state%nvirtual]) - deallocate(iptr) - call c_f_pointer(state%icocc, iptr, [state%icocc_dim]) - deallocate(iptr) - call c_f_pointer(state%icvir, iptr, [state%icvir_dim]) - deallocate(iptr) - call c_f_pointer(state%cocc, rptr, [state%cocc_dim]) - deallocate(rptr) - call c_f_pointer(state%cvir, rptr, [state%cvir_dim]) - deallocate(rptr) - end if + if (state%numat > 0) call destroy_mozyme_state(state) state%numat = numat state%noccupied = noccupied @@ -127,77 +89,15 @@ module subroutine mozyme_save(state) state%cocc_dim = cocc_dim state%cvir_dim = cvir_dim - allocate(iptr(numat), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOZYME_SAVE") - return - end if - iptr = nbonds - state%nbonds = c_loc(iptr) - - allocate(iptr2(9,numat), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOZYME_SAVE") - return - end if - iptr2 = ibonds - state%ibonds = c_loc(iptr2) - - allocate(iptr(numat), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOZYME_SAVE") - return - end if - iptr = iorbs - state%iorbs = c_loc(iptr) - - allocate(iptr(noccupied), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOZYME_SAVE") - return - end if - iptr = ncf - state%ncf = c_loc(iptr) - - allocate(iptr(nvirtual), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOZYME_SAVE") - return - end if - iptr = nce - state%nce = c_loc(iptr) - - allocate(iptr(icocc_dim), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOZYME_SAVE") - return - end if - iptr = icocc - state%icocc = c_loc(iptr) - - allocate(iptr(icvir_dim), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOZYME_SAVE") - return - end if - iptr = icvir - state%icvir = c_loc(iptr) - - allocate(rptr(cocc_dim), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOZYME_SAVE") - return - end if - rptr = cocc - state%cocc = c_loc(rptr) - - allocate(rptr(cvir_dim), stat=status) - if (status /= 0) then - call mopend("Failed to allocate memory in MOZYME_SAVE") - return - end if - rptr = cvir - state%cvir = c_loc(rptr) + 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 subroutine mozyme_save ! load MOZYME density matrix, or construct initial guess diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 8ecb9334..08100971 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -32,6 +32,11 @@ macro(add_mopac_test _name _files) endmacro() #=============================================== +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() # Create a list of source files (src_list) with the .F90 extension set(src_list mopac_api_test @@ -41,6 +46,8 @@ set(src_list 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 ${CMAKE_SOURCE_DIR}/include/mopac_api_internal.F90) +target_sources(mopac-api-test PRIVATE ${CMAKE_SOURCE_DIR}/include/mopac_api_f.F90) #=============================================== add_subdirectory(INDO-dev) 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 From 9f3453e17cd5ca00957341b4afd138127bdf2ef3 Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Wed, 6 Nov 2024 13:04:14 -0500 Subject: [PATCH 10/30] Add missing header/wrapper files for API --- include/LICENSE | 21 ++ include/mopac_api_f.F90 | 187 ++++++++++ include/mopac_api_internal.F90 | 603 +++++++++++++++++++++++++++++++++ src/interface/mopac_api.F90 | 144 +++++--- 4 files changed, 906 insertions(+), 49 deletions(-) create mode 100644 include/LICENSE create mode 100644 include/mopac_api_f.F90 create mode 100644 include/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_f.F90 b/include/mopac_api_f.F90 new file mode 100644 index 00000000..14e7c0fb --- /dev/null +++ b/include/mopac_api_f.F90 @@ -0,0 +1,187 @@ +! 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 subroutine run_mopac_from_input_f(path_to_file) + character(len=240), intent(in) :: path_to_file + end subroutine 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..dc914047 --- /dev/null +++ b/include/mopac_api_internal.F90 @@ -0,0 +1,603 @@ +! 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 + module 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 + module 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 + subroutine run_mopac_from_input(path_to_file) bind(c) + use iso_c_binding + character(kind=c_char), dimension(*), intent(in) :: path_to_file + end subroutine run_mopac_from_input + + end interface + +contains + + 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 + + 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 + + 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 + + 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 + + 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 + + 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 + + subroutine run_mopac_from_input_f(path_to_file) + 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 + call run_mopac_from_input(path_to_file_c) + end subroutine 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(out) :: 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/interface/mopac_api.F90 b/src/interface/mopac_api.F90 index 1baa63fb..294dac4d 100644 --- a/src/interface/mopac_api.F90 +++ b/src/interface/mopac_api.F90 @@ -22,12 +22,12 @@ module mopac_api private ! public derived types of the MOPAC API public :: mopac_system, mopac_properties, mopac_state, mozyme_state - ! public subroutines to create or destroy a MOPAC system - public :: create_mopac_system, destroy_mopac_system - ! public subroutines to destroy the other derived types of the API - public :: destroy_mopac_properties, destroy_mopac_state, destroy_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 @@ -105,10 +105,12 @@ module mopac_api ! 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] + type(c_ptr) :: pb ! real(c_double)[mpack], NULL if uhf == 0 end type ! data that describes the electronic state using localized molecular orbitals @@ -150,49 +152,6 @@ module mopac_api interface - ! allocate memory & initialize mopac_system - module subroutine create_mopac_system(natom, natom_move, charge, spin, model, & - epsilon, atom, coord, nlattice, nlattice_move, & - pressure, lattice, tolerance, max_time, & - system) bind(c) - integer(c_int), intent(in) :: natom - integer(c_int), intent(in) :: natom_move - integer(c_int), intent(in) :: charge - integer(c_int), intent(in) :: spin - integer(c_int), intent(in) :: model - real(c_double), intent(in) :: epsilon - integer(c_int), dimension(natom), intent(in) :: atom - real(c_double), dimension(3*natom), intent(in) :: coord - integer(c_int), intent(in) :: nlattice - integer(c_int), intent(in) :: nlattice_move - real(c_double), intent(in) :: pressure - real(c_double), dimension(3*nlattice), intent(in) :: lattice - real(c_double), intent(in) :: tolerance - integer(c_int), intent(in) :: max_time - type(mopac_system), intent(out) :: system - end subroutine create_mopac_system - - ! deallocate memory in mopac_system - module subroutine destroy_mopac_system(system) bind(c) - type(mopac_system), intent(in) :: system - end subroutine destroy_mopac_system - - ! deallocate memory in mopac_properties (associated system is needed for size info) - module subroutine destroy_mopac_properties(system, properties) bind(c) - type(mopac_system), intent(in) :: system - type(mopac_properties), intent(in) :: properties - end subroutine destroy_mopac_properties - - ! deallocate memory in mopac_state - module subroutine destroy_mopac_state(state) bind(c) - type(mopac_state), intent(in) :: state - end subroutine destroy_mopac_state - - ! deallocate memory in mozyme_state - module subroutine destroy_mozyme_state(state) bind(c) - type(mozyme_state), intent(in) :: state - end subroutine destroy_mozyme_state - ! MOPAC electronic ground state calculation module subroutine mopac_scf(system, state, properties) bind(c) type(mopac_system), intent(in) :: system @@ -235,11 +194,98 @@ module subroutine mozyme_vibe(system, state, properties) bind(c) type(mopac_properties), intent(out) :: properties end subroutine mozyme_vibe + ! allocate memory for mopac_state + module subroutine create_mopac_state(state) bind(c) + type(mopac_state), intent(inout) :: state + end subroutine create_mopac_state + + ! allocate memory for mozyme_state + module subroutine create_mozyme_state(state) bind(c) + type(mozyme_state), intent(inout) :: state + end subroutine create_mozyme_state + + ! deallocate memory in mopac_properties + module subroutine destroy_mopac_properties(properties) bind(c) + type(mopac_properties), intent(in) :: properties + end subroutine destroy_mopac_properties + + ! deallocate memory in mopac_state + module subroutine destroy_mopac_state(state) bind(c) + type(mopac_state), intent(in) :: state + end subroutine destroy_mopac_state + + ! deallocate memory in mozyme_state + module subroutine destroy_mozyme_state(state) bind(c) + type(mozyme_state), intent(in) :: state + end subroutine destroy_mozyme_state + ! run MOPAC conventionally from an input file module subroutine run_mopac_from_input(path_to_file) bind(c) character(kind=c_char), dimension(*), intent(in) :: path_to_file end subroutine 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 From ca9c1c85249d2e2a7d7e301fa7c5997cae113f72 Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Wed, 6 Nov 2024 16:41:52 -0500 Subject: [PATCH 11/30] Fix syntax for ifort's malloc/free --- src/interface/mopac_api_createdestroy.F90 | 65 ++++++++++++++++------- 1 file changed, 47 insertions(+), 18 deletions(-) diff --git a/src/interface/mopac_api_createdestroy.F90 b/src/interface/mopac_api_createdestroy.F90 index 99a687d7..71339158 100644 --- a/src/interface/mopac_api_createdestroy.F90 +++ b/src/interface/mopac_api_createdestroy.F90 @@ -105,8 +105,8 @@ end subroutine destroy_mozyme_state module function create_int(size) integer, intent(in) :: size type(c_ptr) :: create_int - integer(c_int), pointer :: ptr(:) #ifndef MOPAC_API_MALLOC + integer(c_int), pointer :: ptr(:) integer :: status allocate(ptr(size), stat=status) if (status /= 0) then @@ -115,16 +115,18 @@ module function create_int(size) end if create_int = c_loc(ptr) #else - create_int = malloc(c_sizeof(ptr)) - call c_f_pointer(create_int, ptr, [size]) + integer(c_intptr_t) :: dummy + integer(c_int) :: mold + dummy = malloc(c_sizeof(mold)*size(1)) + 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 - integer(c_int), pointer :: ptr(:,:) #ifndef MOPAC_API_MALLOC + integer(c_int), pointer :: ptr(:,:) integer :: status allocate(ptr(size,size2), stat=status) if (status /= 0) then @@ -133,15 +135,17 @@ module function create_int2(size, size2) end if create_int2 = c_loc(ptr) #else - create_int2 = malloc(c_sizeof(ptr)) - call c_f_pointer(create_int2, ptr, [size,size2]) + integer(c_intptr_t) :: dummy + integer(c_int) :: mold + dummy = malloc(c_sizeof(mold)*size(1)*size(2)) + 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 - real(c_double), pointer :: ptr(:) #ifndef MOPAC_API_MALLOC + real(c_double), pointer :: ptr(:) integer :: status allocate(ptr(size), stat=status) if (status /= 0) then @@ -150,8 +154,10 @@ module function create_real(size) end if create_real = c_loc(ptr) #else - create_real = malloc(c_sizeof(ptr)) - call c_f_pointer(create_real, ptr, [size]) + integer(c_intptr_t) :: dummy + real(c_double) :: mold + dummy = malloc(c_sizeof(mold)*size(1)) + create_real = transfer(dummy, create_real) #endif end function create_real @@ -170,7 +176,10 @@ module function create_copy_int(array, size) end if create_copy_int = c_loc(ptr) #else - create_copy_int = malloc(c_sizeof(ptr)) + 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 @@ -189,7 +198,10 @@ module function create_copy_int2(array, size) end if create_copy_int2 = c_loc(ptr) #else - create_copy_int2 = malloc(c_sizeof(ptr)) + 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 @@ -208,7 +220,10 @@ module function create_copy_real(array, size) end if create_copy_real = c_loc(ptr) #else - create_copy_real = malloc(c_sizeof(ptr)) + 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 @@ -228,7 +243,10 @@ module function create_copy_char(array, size) end if create_copy_char = c_loc(ptr) #else - create_copy_char = malloc(c_sizeof(ptr)) + 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 @@ -250,7 +268,10 @@ module function create_copy_ptr(array, size) end if create_copy_ptr = c_loc(ptr) #else - create_copy_ptr = malloc(c_sizeof(ptr)) + 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 @@ -271,7 +292,9 @@ module subroutine destroy_int(copy) end if end if #else - call free(copy) + integer(c_intptr_t) :: copy2 + copy2 = transfer(copy, copy2) + call free(copy2) #endif end subroutine destroy_int module subroutine destroy_real(copy) @@ -288,7 +311,9 @@ module subroutine destroy_real(copy) end if end if #else - call free(copy) + integer(c_intptr_t) :: copy2 + copy2 = transfer(copy, copy2) + call free(copy2) #endif end subroutine destroy_real module subroutine destroy_char(copy) @@ -305,7 +330,9 @@ module subroutine destroy_char(copy) end if end if #else - call free(copy) + integer(c_intptr_t) :: copy2 + copy2 = transfer(copy, copy2) + call free(copy2) #endif end subroutine destroy_char module subroutine destroy_ptr(copy) @@ -322,7 +349,9 @@ module subroutine destroy_ptr(copy) end if end if #else - call free(copy) + integer(c_intptr_t) :: copy2 + copy2 = transfer(copy, copy2) + call free(copy2) #endif end subroutine destroy_ptr From 4e0bf5c60b2ba46dba315c3172897a6731865520 Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Wed, 6 Nov 2024 22:28:16 -0500 Subject: [PATCH 12/30] Fix more syntax for ifort's malloc/free --- src/interface/mopac_api_createdestroy.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/interface/mopac_api_createdestroy.F90 b/src/interface/mopac_api_createdestroy.F90 index 71339158..10a892f1 100644 --- a/src/interface/mopac_api_createdestroy.F90 +++ b/src/interface/mopac_api_createdestroy.F90 @@ -117,7 +117,7 @@ module function create_int(size) #else integer(c_intptr_t) :: dummy integer(c_int) :: mold - dummy = malloc(c_sizeof(mold)*size(1)) + dummy = malloc(c_sizeof(mold)*size) create_int = transfer(dummy, create_int) #endif end function create_int @@ -137,7 +137,7 @@ module function create_int2(size, size2) #else integer(c_intptr_t) :: dummy integer(c_int) :: mold - dummy = malloc(c_sizeof(mold)*size(1)*size(2)) + dummy = malloc(c_sizeof(mold)*size*size2) create_int2 = transfer(dummy, create_int2) #endif end function create_int2 @@ -156,7 +156,7 @@ module function create_real(size) #else integer(c_intptr_t) :: dummy real(c_double) :: mold - dummy = malloc(c_sizeof(mold)*size(1)) + dummy = malloc(c_sizeof(mold)*size) create_real = transfer(dummy, create_real) #endif end function create_real From e9fe88630448e57a03960f0190f8d7d40b80600a Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Wed, 6 Nov 2024 22:35:58 -0500 Subject: [PATCH 13/30] Fix submodule syntax for Fortran API wrapper --- include/mopac_api_internal.F90 | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/include/mopac_api_internal.F90 b/include/mopac_api_internal.F90 index dc914047..143f38b2 100644 --- a/include/mopac_api_internal.F90 +++ b/include/mopac_api_internal.F90 @@ -185,13 +185,13 @@ subroutine mozyme_vibe(system, state, properties) bind(c) end subroutine mozyme_vibe ! allocate memory for mopac_state - module subroutine create_mopac_state(state) bind(c) + 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 - module subroutine create_mozyme_state(state) bind(c) + subroutine create_mozyme_state(state) bind(c) use mopac_api_c type(mozyme_state), intent(inout) :: state end subroutine create_mozyme_state @@ -224,7 +224,7 @@ end subroutine run_mopac_from_input contains - subroutine mopac_scf_f(system, state, properties) + 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 @@ -242,7 +242,7 @@ subroutine mopac_scf_f(system, state, properties) deallocate(rwork) end subroutine mopac_scf_f - subroutine mopac_relax_f(system, state, properties) + 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 @@ -260,7 +260,7 @@ subroutine mopac_relax_f(system, state, properties) deallocate(rwork) end subroutine mopac_relax_f - subroutine mopac_vibe_f(system, state, properties) + 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 @@ -278,7 +278,7 @@ subroutine mopac_vibe_f(system, state, properties) deallocate(rwork) end subroutine mopac_vibe_f - subroutine mozyme_scf_f(system, state, properties) + 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 @@ -296,7 +296,7 @@ subroutine mozyme_scf_f(system, state, properties) deallocate(rwork) end subroutine mozyme_scf_f - subroutine mozyme_relax_f(system, state, properties) + 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 @@ -314,7 +314,7 @@ subroutine mozyme_relax_f(system, state, properties) deallocate(rwork) end subroutine mozyme_relax_f - subroutine mozyme_vibe_f(system, state, properties) + 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 @@ -332,7 +332,7 @@ subroutine mozyme_vibe_f(system, state, properties) deallocate(rwork) end subroutine mozyme_vibe_f - subroutine run_mopac_from_input_f(path_to_file) + module subroutine run_mopac_from_input_f(path_to_file) character(len=240), intent(in) :: path_to_file character(kind=c_char), allocatable :: path_to_file_c(:) integer :: i, size, status @@ -387,7 +387,7 @@ subroutine mopac_system_f2c(system_f, system_c, iwork, rwork) end subroutine mopac_system_f2c subroutine mopac_properties_c2f(system_c, properties_c, properties_f) - type(mopac_system), intent(out) :: system_c + 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(:) From d99326ad47b8d1e092a15a331f670d80f9e3c9c2 Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Thu, 7 Nov 2024 11:22:06 -0500 Subject: [PATCH 14/30] Add dllexport commands w/ preprocessor barriers to API --- CMakeLists.txt | 1 + src/chanel_C.F90 | 2 ++ src/interface/mopac_api.F90 | 3 +++ src/molkst_C.F90 | 2 ++ src/run_mopac.F90 | 2 ++ src/run_param.F90 | 2 ++ 6 files changed, 12 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 28249488..8be586a1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -156,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() 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/mopac_api.F90 b/src/interface/mopac_api.F90 index 294dac4d..20583252 100644 --- a/src/interface/mopac_api.F90 +++ b/src/interface/mopac_api.F90 @@ -16,6 +16,9 @@ ! Diskless/stateless Application Programming Interface (API) to core MOPAC operations module mopac_api +#ifdef WIN32 +!dec$ attributes dllexport :: mopac_scf, mopac_relax, mopac_vibe, mozyme_scf, mozyme_relax, mozyme_vibe +#endif use, intrinsic :: iso_c_binding implicit none 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 e842d1aa..85bfd70b 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 !----------------------------------------------- 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 From 41beccd92fc25efd976dda3b49e1b049e987372f Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Thu, 7 Nov 2024 11:30:17 -0500 Subject: [PATCH 15/30] Add API dllexport commands to submodule --- src/interface/mopac_api_operations.F90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/interface/mopac_api_operations.F90 b/src/interface/mopac_api_operations.F90 index 27a5844c..dd4e1255 100644 --- a/src/interface/mopac_api_operations.F90 +++ b/src/interface/mopac_api_operations.F90 @@ -15,6 +15,9 @@ ! along with this program. If not, see . submodule (mopac_api) mopac_api_operations +#ifdef WIN32 +!dec$ attributes dllexport :: mopac_scf, mopac_relax, mopac_vibe, mozyme_scf, mozyme_relax, mozyme_vibe +#endif use Common_arrays_C, only: xparam, grad, lopt use molkst_C, only: keywrd, escf, moperr, nvar, gui, jobnam, run implicit none From 56662c2afb7edc2545c875e9705854dcc8adbdaf Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Thu, 7 Nov 2024 11:51:27 -0500 Subject: [PATCH 16/30] Remove dllexport commands in scope of subroutine interfaces --- src/interface/mopac_api.F90 | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/interface/mopac_api.F90 b/src/interface/mopac_api.F90 index 20583252..294dac4d 100644 --- a/src/interface/mopac_api.F90 +++ b/src/interface/mopac_api.F90 @@ -16,9 +16,6 @@ ! Diskless/stateless Application Programming Interface (API) to core MOPAC operations module mopac_api -#ifdef WIN32 -!dec$ attributes dllexport :: mopac_scf, mopac_relax, mopac_vibe, mozyme_scf, mozyme_relax, mozyme_vibe -#endif use, intrinsic :: iso_c_binding implicit none From d073ba8570ba3dc41c6d61e04f5ed2ff7f3c6c64 Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Thu, 7 Nov 2024 12:08:33 -0500 Subject: [PATCH 17/30] Move dllexport commands to parent module interface block --- src/interface/mopac_api.F90 | 36 ++++++++++++++++++++++++++ src/interface/mopac_api_operations.F90 | 3 --- 2 files changed, 36 insertions(+), 3 deletions(-) diff --git a/src/interface/mopac_api.F90 b/src/interface/mopac_api.F90 index 294dac4d..4d969449 100644 --- a/src/interface/mopac_api.F90 +++ b/src/interface/mopac_api.F90 @@ -154,6 +154,9 @@ module mopac_api ! MOPAC electronic ground state calculation 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 @@ -161,6 +164,9 @@ end subroutine mopac_scf ! MOPAC geometry relaxation 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 @@ -168,6 +174,9 @@ end subroutine mopac_relax ! MOPAC vibrational calculation 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 @@ -175,6 +184,9 @@ end subroutine mopac_vibe ! MOZYME electronic ground state calculation 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 @@ -182,6 +194,9 @@ end subroutine mozyme_scf ! MOZYME geometry relaxation 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 @@ -189,6 +204,9 @@ end subroutine mozyme_relax ! MOZYME vibrational calculation 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 @@ -196,31 +214,49 @@ end subroutine mozyme_vibe ! 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 subroutine run_mopac_from_input(path_to_file) bind(c) +#ifdef WIN32 +!dec$ attributes dllexport :: run_mopac_from_input +#endif character(kind=c_char), dimension(*), intent(in) :: path_to_file end subroutine run_mopac_from_input diff --git a/src/interface/mopac_api_operations.F90 b/src/interface/mopac_api_operations.F90 index dd4e1255..27a5844c 100644 --- a/src/interface/mopac_api_operations.F90 +++ b/src/interface/mopac_api_operations.F90 @@ -15,9 +15,6 @@ ! along with this program. If not, see . submodule (mopac_api) mopac_api_operations -#ifdef WIN32 -!dec$ attributes dllexport :: mopac_scf, mopac_relax, mopac_vibe, mozyme_scf, mozyme_relax, mozyme_vibe -#endif use Common_arrays_C, only: xparam, grad, lopt use molkst_C, only: keywrd, escf, moperr, nvar, gui, jobnam, run implicit none From 6a67b289b0b812e6d189a2d11985273aa7561cd0 Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Thu, 7 Nov 2024 12:26:00 -0500 Subject: [PATCH 18/30] Add dllexport commands to subroutine bodies in submodules --- src/interface/mopac_api_createdestroy.F90 | 15 +++++++++++++++ src/interface/mopac_api_operations.F90 | 21 +++++++++++++++++++++ 2 files changed, 36 insertions(+) diff --git a/src/interface/mopac_api_createdestroy.F90 b/src/interface/mopac_api_createdestroy.F90 index 10a892f1..fae1222a 100644 --- a/src/interface/mopac_api_createdestroy.F90 +++ b/src/interface/mopac_api_createdestroy.F90 @@ -22,6 +22,9 @@ ! 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) @@ -31,6 +34,9 @@ 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) @@ -47,6 +53,9 @@ 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(:) @@ -71,6 +80,9 @@ 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) @@ -80,6 +92,9 @@ 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) diff --git a/src/interface/mopac_api_operations.F90 b/src/interface/mopac_api_operations.F90 index 27a5844c..8672ec41 100644 --- a/src/interface/mopac_api_operations.F90 +++ b/src/interface/mopac_api_operations.F90 @@ -57,6 +57,9 @@ end subroutine mozyme_load ! MOPAC electronic ground state calculation 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,6 +75,9 @@ end subroutine mopac_scf ! MOPAC geometry relaxation 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 @@ -87,6 +93,9 @@ end subroutine mopac_relax ! MOPAC vibrational calculation 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 @@ -112,6 +121,9 @@ end subroutine mopac_vibe ! MOZYME electronic ground state calculation 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 @@ -127,6 +139,9 @@ end subroutine mozyme_scf ! MOZYME geometry relaxation 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 @@ -142,6 +157,9 @@ end subroutine mozyme_relax ! MOZYME vibrational calculation 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 @@ -167,6 +185,9 @@ end subroutine mozyme_vibe ! Run MOPAC conventionally from an input file module subroutine run_mopac_from_input(path_to_file) bind(c) +#ifdef WIN32 +!dec$ attributes dllexport :: run_mopac_from_input +#endif character(kind=c_char), dimension(*), intent(in) :: path_to_file integer :: i i = 1 From 8caa8b5986763312a05292aeee2f5afc93b83315 Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Thu, 7 Nov 2024 14:13:43 -0500 Subject: [PATCH 19/30] Suppress stdout output for MOPAC when run as API call --- src/run_mopac.F90 | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/run_mopac.F90 b/src/run_mopac.F90 index 85bfd70b..f63f53dd 100644 --- a/src/run_mopac.F90 +++ b/src/run_mopac.F90 @@ -34,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 ! @@ -1013,8 +1013,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 ! From 8f12de86f2f0e9fd35038999ebedb32f4f024027 Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Thu, 7 Nov 2024 14:23:14 -0500 Subject: [PATCH 20/30] Add return status to MOPAC API call --- include/mopac_api.h | 2 +- include/mopac_api_f.F90 | 5 +++-- include/mopac_api_internal.F90 | 17 ++++++++++++----- src/interface/mopac_api.F90 | 5 +++-- src/interface/mopac_api_initialize.F90 | 7 ++++--- src/interface/mopac_api_operations.F90 | 10 ++++++++-- 6 files changed, 31 insertions(+), 15 deletions(-) diff --git a/include/mopac_api.h b/include/mopac_api.h index 177e5219..0bb453ed 100644 --- a/include/mopac_api.h +++ b/include/mopac_api.h @@ -171,6 +171,6 @@ void destroy_mopac_state(struct mopac_state *state); void destroy_mozyme_state(struct mozyme_state *state); /* run MOPAC conventionally from an input file */ -void run_mopac_from_input(char *path_to_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 index 14e7c0fb..d05d95de 100644 --- a/include/mopac_api_f.F90 +++ b/include/mopac_api_f.F90 @@ -178,9 +178,10 @@ module subroutine mozyme_vibe_f(system, state, properties) end subroutine mozyme_vibe_f ! run MOPAC conventionally from an input file - module subroutine run_mopac_from_input_f(path_to_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 subroutine run_mopac_from_input_f + end function run_mopac_from_input_f end interface diff --git a/include/mopac_api_internal.F90 b/include/mopac_api_internal.F90 index 143f38b2..9d2814b0 100644 --- a/include/mopac_api_internal.F90 +++ b/include/mopac_api_internal.F90 @@ -215,10 +215,11 @@ subroutine destroy_mozyme_state(state) bind(c) end subroutine destroy_mozyme_state ! run MOPAC conventionally from an input file - subroutine run_mopac_from_input(path_to_file) bind(c) + 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 subroutine run_mopac_from_input + end function run_mopac_from_input end interface @@ -332,7 +333,8 @@ module subroutine mozyme_vibe_f(system, state, properties) deallocate(rwork) end subroutine mozyme_vibe_f - module subroutine run_mopac_from_input_f(path_to_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 character(kind=c_char), allocatable :: path_to_file_c(:) integer :: i, size, status @@ -346,8 +348,13 @@ module subroutine run_mopac_from_input_f(path_to_file) path_to_file_c(i) = path_to_file(i:i) end do path_to_file_c(size+1) = c_null_char - call run_mopac_from_input(path_to_file_c) - end subroutine run_mopac_from_input_f + 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 diff --git a/src/interface/mopac_api.F90 b/src/interface/mopac_api.F90 index 4d969449..e222a16a 100644 --- a/src/interface/mopac_api.F90 +++ b/src/interface/mopac_api.F90 @@ -253,12 +253,13 @@ module subroutine destroy_mozyme_state(state) bind(c) end subroutine destroy_mozyme_state ! run MOPAC conventionally from an input file - module subroutine run_mopac_from_input(path_to_file) bind(c) + 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 subroutine run_mopac_from_input + end function run_mopac_from_input end interface diff --git a/src/interface/mopac_api_initialize.F90 b/src/interface/mopac_api_initialize.F90 index 24c0618e..1b44bf69 100644 --- a/src/interface/mopac_api_initialize.F90 +++ b/src/interface/mopac_api_initialize.F90 @@ -135,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) diff --git a/src/interface/mopac_api_operations.F90 b/src/interface/mopac_api_operations.F90 index 8672ec41..2ce535ab 100644 --- a/src/interface/mopac_api_operations.F90 +++ b/src/interface/mopac_api_operations.F90 @@ -184,10 +184,11 @@ module subroutine mozyme_vibe(system, state, properties) bind(c) end subroutine mozyme_vibe ! Run MOPAC conventionally from an input file - module subroutine run_mopac_from_input(path_to_file) bind(c) + 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 @@ -199,9 +200,14 @@ module subroutine run_mopac_from_input(path_to_file) bind(c) 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 From 59980c4310c14c29bc48a1df8ac4f6305e179c7b Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Thu, 7 Nov 2024 15:36:19 -0500 Subject: [PATCH 21/30] Remove state allocations for empty states --- src/interface/mopac_api_saveload.F90 | 48 +++++++++++++++------------- 1 file changed, 26 insertions(+), 22 deletions(-) diff --git a/src/interface/mopac_api_saveload.F90 b/src/interface/mopac_api_saveload.F90 index 4817d1af..e80cce80 100644 --- a/src/interface/mopac_api_saveload.F90 +++ b/src/interface/mopac_api_saveload.F90 @@ -30,13 +30,15 @@ module subroutine mopac_save(state) if (state%mpack > 0) call destroy_mopac_state(state) state%mpack = mpack - 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 + 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 end if end subroutine mopac_save @@ -82,22 +84,24 @@ module subroutine mozyme_save(state) if (state%numat > 0) call destroy_mozyme_state(state) 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 (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]) + 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 end subroutine mozyme_save ! load MOZYME density matrix, or construct initial guess From d05b447fa8bde4858a40c6d890a5564c628833c4 Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Thu, 7 Nov 2024 15:36:41 -0500 Subject: [PATCH 22/30] Temporary debug printing --- include/mopac_api_internal.F90 | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/include/mopac_api_internal.F90 b/include/mopac_api_internal.F90 index 9d2814b0..42234e1d 100644 --- a/include/mopac_api_internal.F90 +++ b/include/mopac_api_internal.F90 @@ -534,14 +534,23 @@ subroutine mozyme_state_c2f(state_c, state_f) if (allocated(state_f%cocc)) deallocate(state_f%cocc) if (allocated(state_f%cvir)) deallocate(state_f%cvir) if (state_f%numat /= 0) then +write(*,*) "nbonds", state_c%numat call copy_int(state_f%nbonds, state_c%nbonds, [state_c%numat]) +write(*,*) "ibonds" call copy_int2(state_f%ibonds, state_c%ibonds, [9,state_c%numat]) +write(*,*) "iorbs" call copy_int(state_f%iorbs, state_c%iorbs, [state_c%numat]) +write(*,*) "ncf", state_c%noccupied call copy_int(state_f%ncf, state_c%ncf, [state_c%noccupied]) +write(*,*) "nce", state_c%nvirtual call copy_int(state_f%nce, state_c%nce, [state_c%nvirtual]) +write(*,*) "icocc", state_c%icocc_dim call copy_int(state_f%icocc, state_c%icocc, [state_c%icocc_dim]) +write(*,*) "icvir", state_c%icvir_dim call copy_int(state_f%icvir, state_c%icvir, [state_c%icvir_dim]) +write(*,*) "cocc", state_c%cocc_dim call copy_real(state_f%cocc, state_c%cocc, [state_c%cocc_dim]) +write(*,*) "cvir", state_c%cvir_dim call copy_real(state_f%cvir, state_c%cvir, [state_c%cvir_dim]) end if call destroy_mozyme_state(state_c) From 4972293c27061c8dbb3fed4998e39efd33aa0deb Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Thu, 7 Nov 2024 15:49:54 -0500 Subject: [PATCH 23/30] Adjust intent of API save operations --- src/interface/mopac_api_operations.F90 | 4 ++-- src/interface/mopac_api_saveload.F90 | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/interface/mopac_api_operations.F90 b/src/interface/mopac_api_operations.F90 index 2ce535ab..27352002 100644 --- a/src/interface/mopac_api_operations.F90 +++ b/src/interface/mopac_api_operations.F90 @@ -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 diff --git a/src/interface/mopac_api_saveload.F90 b/src/interface/mopac_api_saveload.F90 index e80cce80..2bd4ec4e 100644 --- a/src/interface/mopac_api_saveload.F90 +++ b/src/interface/mopac_api_saveload.F90 @@ -25,7 +25,7 @@ ! save MOPAC density matrices module subroutine mopac_save(state) - type(mopac_state), intent(out) :: state + type(mopac_state), intent(inout) :: state if (state%mpack > 0) call destroy_mopac_state(state) @@ -79,7 +79,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 if (state%numat > 0) call destroy_mozyme_state(state) From 2167de94eedb0d6ac481c61bd0a9c821c63a712c Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Thu, 7 Nov 2024 15:58:53 -0500 Subject: [PATCH 24/30] Zero-pad MOZYME state variables when their memory footprint is expanded --- src/MOZYME/iter_for_MOZYME.F90 | 8 ++++++++ 1 file changed, 8 insertions(+) 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 From 1938c33bef6feddf1357b66fc93f66d487eb6338 Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Fri, 8 Nov 2024 15:25:52 -0500 Subject: [PATCH 25/30] Add valgrind testing to API tester on Linux --- .github/workflows/CI.yaml | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index 05739e3d..cacf3f00 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -54,6 +54,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 --leak-check=yes ./mopac-api-test + - name: Test MOPAC with CTest run: | cd build @@ -99,6 +105,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 --leak-check=yes ./mopac-api-test + - name: Test MOPAC with CTest run: | cd build From 6cd4e1ade436a26425a4aee9f50a0163c4e89d0f Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Fri, 8 Nov 2024 15:33:35 -0500 Subject: [PATCH 26/30] Switch non-production Linux builds to debug --- .github/workflows/CI.yaml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index cacf3f00..89e0df2a 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 @@ -98,6 +99,7 @@ jobs: - name: Configure MOPAC with CMake run: | cmake -B build \ + -DCMAKE_BUILD_TYPE=Debug \ -DTHREADS_KEYWORD=OFF \ -DENABLE_COVERAGE=ON From 7451b891a38274d7921d60b0f0f32d9b7df1629d Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Fri, 8 Nov 2024 15:44:36 -0500 Subject: [PATCH 27/30] Request more info from valgrind --- .github/workflows/CI.yaml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index 89e0df2a..00f88e15 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -59,7 +59,7 @@ jobs: run: | sudo apt update && sudo apt install -y valgrind cd build/tests - valgrind --leak-check=yes ./mopac-api-test + valgrind -s --track-origins=yes --leak-check=yes ./mopac-api-test - name: Test MOPAC with CTest run: | @@ -111,7 +111,7 @@ jobs: run: | sudo apt update && sudo apt install -y valgrind cd build/tests - valgrind --leak-check=yes ./mopac-api-test + valgrind -s --track-origins=yes --leak-check=yes ./mopac-api-test - name: Test MOPAC with CTest run: | From 644661182b48d06187646012225410421839cb62 Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Fri, 8 Nov 2024 16:51:36 -0500 Subject: [PATCH 28/30] Don't trust stride when copying internal MOPAC arrays to API I/O --- src/interface/mopac_api_createdestroy.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/interface/mopac_api_createdestroy.F90 b/src/interface/mopac_api_createdestroy.F90 index fae1222a..c4863b4f 100644 --- a/src/interface/mopac_api_createdestroy.F90 +++ b/src/interface/mopac_api_createdestroy.F90 @@ -197,7 +197,7 @@ module function create_copy_int(array, size) create_copy_int = transfer(dummy, create_copy_int) call c_f_pointer(create_copy_int, ptr, size) #endif - ptr = array + ptr = array(:size(1)) end function create_copy_int module function create_copy_int2(array, size) integer, intent(in) :: array(:,:) @@ -219,7 +219,7 @@ module function create_copy_int2(array, size) create_copy_int2 = transfer(dummy, create_copy_int2) call c_f_pointer(create_copy_int2, ptr, size) #endif - ptr = array + ptr = array(:size(1),:size(2)) end function create_copy_int2 module function create_copy_real(array, size) double precision, intent(in) :: array(:) @@ -241,7 +241,7 @@ module function create_copy_real(array, size) create_copy_real = transfer(dummy, create_copy_real) call c_f_pointer(create_copy_real, ptr, size) #endif - ptr = array + ptr = array(:size(1)) end function create_copy_real module function create_copy_char(array, size) character(len=*), intent(in) :: array @@ -289,7 +289,7 @@ module function create_copy_ptr(array, size) create_copy_ptr = transfer(dummy, create_copy_ptr) call c_f_pointer(create_copy_ptr, ptr, size) #endif - ptr = array + ptr = array(:size(1)) end function create_copy_ptr ! deallocate memory (C or Fortran memory manager, depending on compiler) From 6ccf65dffe568146f7574c138cff1a29c5f454aa Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Fri, 8 Nov 2024 18:41:18 -0500 Subject: [PATCH 29/30] Deactivate valgrind tests --- .github/workflows/CI.yaml | 20 ++++++++++---------- include/mopac_api_internal.F90 | 9 --------- 2 files changed, 10 insertions(+), 19 deletions(-) diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index 00f88e15..abd7784f 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -55,11 +55,11 @@ 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: 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: | @@ -107,11 +107,11 @@ 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: 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: | diff --git a/include/mopac_api_internal.F90 b/include/mopac_api_internal.F90 index 42234e1d..9d2814b0 100644 --- a/include/mopac_api_internal.F90 +++ b/include/mopac_api_internal.F90 @@ -534,23 +534,14 @@ subroutine mozyme_state_c2f(state_c, state_f) if (allocated(state_f%cocc)) deallocate(state_f%cocc) if (allocated(state_f%cvir)) deallocate(state_f%cvir) if (state_f%numat /= 0) then -write(*,*) "nbonds", state_c%numat call copy_int(state_f%nbonds, state_c%nbonds, [state_c%numat]) -write(*,*) "ibonds" call copy_int2(state_f%ibonds, state_c%ibonds, [9,state_c%numat]) -write(*,*) "iorbs" call copy_int(state_f%iorbs, state_c%iorbs, [state_c%numat]) -write(*,*) "ncf", state_c%noccupied call copy_int(state_f%ncf, state_c%ncf, [state_c%noccupied]) -write(*,*) "nce", state_c%nvirtual call copy_int(state_f%nce, state_c%nce, [state_c%nvirtual]) -write(*,*) "icocc", state_c%icocc_dim call copy_int(state_f%icocc, state_c%icocc, [state_c%icocc_dim]) -write(*,*) "icvir", state_c%icvir_dim call copy_int(state_f%icvir, state_c%icvir, [state_c%icvir_dim]) -write(*,*) "cocc", state_c%cocc_dim call copy_real(state_f%cocc, state_c%cocc, [state_c%cocc_dim]) -write(*,*) "cvir", state_c%cvir_dim call copy_real(state_f%cvir, state_c%cvir, [state_c%cvir_dim]) end if call destroy_mozyme_state(state_c) From cd0d325bdc1e3e8c375646966cfb7701b596b5ae Mon Sep 17 00:00:00 2001 From: Jonathan Moussa Date: Fri, 8 Nov 2024 18:50:39 -0500 Subject: [PATCH 30/30] Localize CMake source file inclusion for the API tester --- CMakeLists.txt | 3 +++ include/CMakeLists.txt | 18 ++++++++++++++++++ tests/CMakeLists.txt | 19 +------------------ 3 files changed, 22 insertions(+), 18 deletions(-) create mode 100644 include/CMakeLists.txt diff --git a/CMakeLists.txt b/CMakeLists.txt index 8be586a1..bf6a2d0d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -268,6 +268,9 @@ if(TESTS) 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) + add_subdirectory(include) add_subdirectory(tests) endif() 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/tests/CMakeLists.txt b/tests/CMakeLists.txt index 08100971..04b48ae2 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -31,24 +31,7 @@ macro(add_mopac_test _name _files) ${CMAKE_BINARY_DIR}/mopac${CMAKE_EXECUTABLE_SUFFIX} ${HOF_ERROR} ${_files}) endmacro() -#=============================================== -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() -# 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 ${CMAKE_SOURCE_DIR}/include/mopac_api_internal.F90) -target_sources(mopac-api-test PRIVATE ${CMAKE_SOURCE_DIR}/include/mopac_api_f.F90) -#=============================================== +target_sources(mopac-api-test PRIVATE mopac_api_test.F90) add_subdirectory(INDO-dev) add_subdirectory(keywords)