diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index c318e9c4..20f8b0f4 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -20,6 +20,10 @@ jobs: cmargs: > -D BUILD_SHARED_LIBS=ON -D LIBEFP_ENABLE_OPENMP=ON + # uncomment the following two options to look for possible memory-related bugs in openmp + # however, the tests will take ~10 times longer + #-D CMAKE_C_FLAGS="-fsanitize=thread -g" + #-D CMAKE_EXE_LINKER_FLAGS="-fsanitize=thread" pytest-marker-expr: "test" # i.e., all - label: L-Gnu @@ -30,6 +34,11 @@ jobs: cmargs: > -D BUILD_SHARED_LIBS=ON -D LIBEFP_ENABLE_OPENMP=ON + # uncomment the following two options to look for possible memory-related bugs in openmp + # however, the tests will take ~10 times longer + # this workflow was crashing with these options in pytests due to conflict of TSan and libgomp + #-D CMAKE_C_FLAGS="-fsanitize=thread -g" + #-D CMAKE_EXE_LINKER_FLAGS="-fsanitize=thread" pytest-marker-expr: "test" - label: L-Intel @@ -48,18 +57,20 @@ jobs: - label: M-Clang # NaNs in tests on macos-latest (macos-12) - runs-on: macos-13 + #runs-on: macos-13 + runs-on: macos-latest python-version: "3.10" blas: OBL build_type: Release cmargs: > -D BUILD_SHARED_LIBS=ON -D LIBEFP_ENABLE_OPENMP=ON - pytest-marker-expr: "test" + pytest-marker-expr: "test" - label: M-Clang # NaNs in tests on macos-latest (macos-12) - runs-on: macos-13 + #runs-on: macos-13 + runs-on: macos-15-intel python-version: "3.10" blas: ACC build_type: Release @@ -124,11 +135,20 @@ jobs: : sed -E -i.bak "s;#${{ matrix.cfg.blas }};;g" export.yaml fi + # if [[ "${{ runner.os }}" == "Windows" ]]; then + # : + # sed -i "s;fortran-compiler;m2w64-gcc-fortran;g" export.yaml + # sed -i "s;#${{ matrix.cfg.blas }};;g" export.yaml + # sed -i "s;openmp;pthreads;g" export.yaml # W openblas is pthreads + # sed -i "s;#- libpython;- libpython;g" export.yaml + # fi if [[ "${{ runner.os }}" == "Windows" ]]; then - : - sed -i "s;fortran-compiler;m2w64-gcc-fortran;g" export.yaml + # Use the modern conda-forge compilers instead of the old m2w64 ones + sed -i "s;c-compiler;gcc;g" export.yaml + sed -i "s;cxx-compiler;gxx;g" export.yaml + sed -i "s;fortran-compiler;gfortran;g" export.yaml sed -i "s;#${{ matrix.cfg.blas }};;g" export.yaml - sed -i "s;openmp;pthreads;g" export.yaml # W openblas is pthreads + sed -i "s;openmp;pthreads;g" export.yaml sed -i "s;#- libpython;- libpython;g" export.yaml fi # model sed for L/W @@ -182,6 +202,17 @@ jobs: - name: Test (CTest) - unit tests run: ctest --output-on-failure --test-dir "${{ github.workspace }}/build" + # - name: Test (CTest) with Valgrind + # run: | + # sudo apt-get update && sudo apt-get install -y valgrind + # # Run valgrind directly on the ctest executable. The tests will take forever + # valgrind --tool=memcheck \ + # --leak-check=full \ + # --track-origins=yes \ + # --trace-children=yes \ + # --error-exitcode=1 \ + # ctest --output-on-failure --test-dir "${{ github.workspace }}/build" + - name: Test (find_package) - consume installation run: | mkdir test_installed_library && cd test_installed_library @@ -242,11 +273,19 @@ jobs: - name: Install Psi4 for QM/EFP testing if: ${{ matrix.cfg.blas == 'MKL' }} run: conda install psi4 -c conda-forge - + - name: Test (pytest) -- unit tests Python bindings run: | - PYTHONPATH="${{ github.workspace }}/installed/lib" \ - pytest --cache-clear -v -rws --color=yes \ - --durations=50 --durations-min=1 --strict-markers \ - -k "${{ matrix.cfg.pytest-marker-expr }}" \ - "${{ github.workspace }}/installed/" + # Setup PYTHONPATH + SEP=$(python -c "import os; print(os.pathsep)") + LIB_PATH="${{ github.workspace }}/installed/lib" + TEST_PATH="${{ github.workspace }}/installed/lib/pylibefp/tests" + export PYTHONPATH="${LIB_PATH}${SEP}${TEST_PATH}" + echo "PYTHONPATH is set to: $PYTHONPATH" + # PYTHONPATH="${{ github.workspace }}/installed/lib:${{ github.workspace }}/installed/lib/pylibefp/tests" \ + + # Run Pytest + pytest -s --cache-clear -v -rws --color=yes \ + --durations=50 --durations-min=1 --strict-markers \ + -k "${{ matrix.cfg.pytest-marker-expr }}" \ + "${{ github.workspace }}/installed/" diff --git a/CMakeLists.txt b/CMakeLists.txt index 87b1d09c..94a71967 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -32,7 +32,7 @@ project( LANGUAGES C CXX ) -set(CMAKE_INSTALL_PREFIX "$ENV{LIBEFP_DIR}" CACHE PATH "Installation directory" FORCE) +#set(CMAKE_INSTALL_PREFIX "$ENV{LIBEFP_DIR}" CACHE PATH "Installation directory" FORCE) #set(CMAKE_INSTALL_PREFIX "${LIBEFP_DIR}" CACHE PATH "Installation directory") #set(CMAKE_INSTALL_LIBDIR "lib") @@ -81,7 +81,7 @@ if((${BUILD_SHARED_LIBS}) AND NOT ${BUILD_FPIC}) endif() ### This option turns OPENMP ON and OFF! -option_with_print(LIBEFP_ENABLE_OPENMP "Enable OpenMP parallelization. Psi4 wants OFF" OFF) +option_with_print(LIBEFP_ENABLE_OPENMP "Enable OpenMP parallelization. Psi4 wants OFF" ON) option_with_print(ENABLE_GENERIC "Enable mostly static linking in shared library" OFF) include(xhost) # defines: option(ENABLE_XHOST "Enable processor-specific optimization" ON) diff --git a/efpmd/src/common.c b/efpmd/src/common.c index defdc5bf..fd899c8a 100644 --- a/efpmd/src/common.c +++ b/efpmd/src/common.c @@ -146,20 +146,24 @@ void print_geometry_pbc(struct efp *efp, int ligand) size_t n_frags; check_fail(efp_get_frag_count(efp, &n_frags)); + if (ligand < 0) + msg(" WARNING! Specify ligand tor printing PBC geometry\n\n"); + size_t lig = (size_t)ligand; + msg(" GEOMETRY IN PBC CELL (ANGSTROMS)\n\n"); - double bx[6]; + double bx[6] = {0.0}; check_fail(efp_get_periodic_box(efp, bx)); six_t box = {bx[0],bx[1],bx[2],bx[3],bx[4],bx[5]}; //six_t box = {10.66, 12.03, 10.872, 90.0, 115.83, 90.0}; double lig_com[6]; - check_fail(efp_get_frag_xyzabc(efp,ligand,lig_com)); + check_fail(efp_get_frag_xyzabc(efp,lig,lig_com)); for (size_t i = 0; i < n_frags; i++) { vec_t cell = {0.0,0.0,0.0}; - if (i != ligand) { + if (i != lig) { double frag_com[6]; check_fail(efp_get_frag_xyzabc(efp, i, frag_com)); vec_t dr = {frag_com[0] - lig_com[0], frag_com[1] - lig_com[1], frag_com[2] - lig_com[2]}; @@ -199,7 +203,7 @@ void print_geometry_pbc(struct efp *efp, int ligand) void print_energy(struct state *state) { // printf("Inside print_energy\n"); - struct efp_energy energy; + struct efp_energy energy = {0.0}; check_fail(efp_get_energy(state->efp, &energy)); @@ -346,13 +350,15 @@ void print_pair_energy(struct state *state) { check_fail(efp_get_coordinates(state->efp, coord)); struct efp_energy *energies; - energies = xmalloc(n_frags * sizeof(struct efp_energy)); + energies = xcalloc(n_frags, sizeof(struct efp_energy)); check_fail(efp_get_pairwise_energy(state->efp, energies)); char ligand[32]; size_t lig_atoms; - size_t ligand_index = cfg_get_int(state->cfg, "ligand"); + size_t ligand_index = (size_t)cfg_get_int(state->cfg, "ligand"); + if (ligand_index >= n_frags) + error("Ligand index is out of bound"); check_fail(efp_get_frag_name(state->efp, ligand_index, sizeof(ligand),ligand)); check_fail(efp_get_frag_atom_count(state->efp, ligand_index, &lig_atoms)); struct efp_atom latoms[lig_atoms]; @@ -360,11 +366,9 @@ void print_pair_energy(struct state *state) { char frag_name[32]; size_t frag_atoms; - double lattice_energy[6]; - for (size_t j=0; j<6; j++){ - lattice_energy[j]=0.0; - } - for (size_t i=0; i efp, i, sizeof(frag_name),frag_name)); check_fail(efp_get_frag_atom_count(state->efp, i, &frag_atoms)); diff --git a/efpmd/src/efield.c b/efpmd/src/efield.c index 54605ae9..57cbc194 100644 --- a/efpmd/src/efield.c +++ b/efpmd/src/efield.c @@ -63,12 +63,12 @@ sim_efield(struct state *state) msg("ELECTRIC FIELD IS IN ATOMIC UNITS\n\n"); for (size_t i = 0; i < n_frags; i++) { - double field[3]; + double field[3] = {0.0}; struct efp_atom *atoms; size_t n_atoms; check_fail(efp_get_frag_atom_count(state->efp, i, &n_atoms)); - atoms = xmalloc(n_atoms * sizeof(struct efp_atom)); + atoms = xcalloc(n_atoms, sizeof(struct efp_atom)); check_fail(efp_get_frag_atoms(state->efp, i, n_atoms, atoms)); for (size_t j = 0; j < n_atoms; j++) { @@ -116,7 +116,7 @@ sim_elpot(struct state *state) size_t n_atoms; check_fail(efp_get_frag_atom_count(state->efp, i, &n_atoms)); - atoms = xmalloc(n_atoms * sizeof(struct efp_atom)); + atoms = xcalloc(n_atoms, sizeof(struct efp_atom)); check_fail(efp_get_frag_atoms(state->efp, i, n_atoms, atoms)); msg("ELECTROSTATIC POTENTIAL ON FRAGMENT %zu\n", i); @@ -133,38 +133,13 @@ sim_elpot(struct state *state) msg("ELECTROSTATIC POTENTIAL JOB COMPLETED SUCCESSFULLY\n"); } -/* -void get_frag_elpot(struct state *state) { - size_t spec_frag; - spec_frag = cfg_get_int(state->cfg, "special_fragment"); - - double elpot; - struct efp_atom *atoms; - size_t n_atoms; - - - check_fail(efp_get_frag_atom_count(state->efp, spec_frag, &n_atoms)); // SKP - atoms = xmalloc(n_atoms * sizeof(struct efp_atom)); - check_fail(efp_get_frag_atoms(state->efp, spec_frag, n_atoms, atoms)); - //state->spec_elpot = malloc(n_atoms * sizeof(double)); - state->spec_elpot = xcalloc(n_atoms, sizeof(double)); - - - for (size_t j = 0; j < n_atoms; j++) { - check_fail(efp_get_elec_potential(state->efp, spec_frag, &atoms[j].x, &elpot)); - state->spec_elpot[j] = elpot; - print_elpot(atoms + j, elpot); - } - - free(atoms); -} -*/ - void sim_frag_elpot(struct state *state) { - // size_t chosen_frag = cfg_get_int(state->cfg, "frag_num"); + if (cfg_get_int(state->cfg, "special_fragment") < 0) + error("Special fragment is not defined in fragment elec potential calculation"); + size_t spec_frag; - spec_frag = cfg_get_int(state->cfg, "special_fragment"); + spec_frag = (size_t)cfg_get_int(state->cfg, "special_fragment"); msg("\n=============FRAG-ELECTROSTATIC POTENTIAL JOB======================\n\n"); @@ -172,13 +147,13 @@ void sim_frag_elpot(struct state *state) { msg("\nCOORDINATES IN ANGSTROMS, ELECTROSTATIC POTENTIAL IN ATOMIC UNITS\n"); msg(" ATOM X Y Z ELPOT \n\n"); - double elpot; + double elpot = 0.0; struct efp_atom *atoms; size_t n_atoms; - check_fail(efp_get_frag_atom_count(state->efp, spec_frag, &n_atoms)); // SKP - atoms = xmalloc(n_atoms * sizeof(struct efp_atom)); + check_fail(efp_get_frag_atom_count(state->efp, spec_frag, &n_atoms)); + atoms = xcalloc(n_atoms, sizeof(struct efp_atom)); check_fail(efp_get_frag_atoms(state->efp, spec_frag, n_atoms, atoms)); msg("ELECTROSTATIC POTENTIAL ON FRAGMENT %zu\n", spec_frag); diff --git a/efpmd/src/energy.c b/efpmd/src/energy.c index e61dcc89..9c6de4e8 100644 --- a/efpmd/src/energy.c +++ b/efpmd/src/energy.c @@ -37,8 +37,8 @@ void compute_energy(struct state *state, bool do_grad) { struct efp_atom *atoms; - struct efp_energy efp_energy; - double xyz[3], xyzabc[6], *grad; + struct efp_energy efp_energy = {0.0}; + double xyz[3] = {0.0}, xyzabc[6] = {0.0}, *grad; size_t ifrag, nfrag, iatom, natom, spec_frag, n_special_atoms; int itotal; @@ -94,13 +94,16 @@ void compute_energy(struct state *state, bool do_grad) if (cfg_get_bool(state->cfg, "enable_torch") && cfg_get_int(state->cfg, "opt_special_frag") > -1) { - spec_frag = cfg_get_int(state->cfg, "special_fragment"); - check_fail(efp_get_frag_atom_count(state->efp, spec_frag, &n_special_atoms)); // SKP + spec_frag = (size_t)cfg_get_int(state->cfg, "special_fragment"); + check_fail(efp_get_frag_atom_count(state->efp, spec_frag, &n_special_atoms)); + + if (n_special_atoms < 1) + error("ML special fragment does not have any atoms!"); if (cfg_get_bool(state->cfg, "enable_elpot")) { double *elpot; struct efp_atom *atoms; - atoms = xmalloc(n_special_atoms * sizeof(struct efp_atom)); + atoms = xcalloc(n_special_atoms, sizeof(struct efp_atom)); check_fail(efp_get_frag_atoms(state->efp, spec_frag, n_special_atoms, atoms)); elpot = xcalloc(n_special_atoms, sizeof(double)); for (size_t j = 0; j < n_special_atoms; j++) { @@ -189,7 +192,7 @@ void compute_energy(struct state *state, bool do_grad) for (ifrag = 0, itotal = 0; ifrag < nfrag; ifrag++) { check_fail(efp_get_frag_atom_count(state->efp, ifrag, &natom)); - atoms = xmalloc(natom * sizeof(struct efp_atom)); + atoms = xcalloc(natom, sizeof(struct efp_atom)); check_fail(efp_get_frag_atoms(state->efp, ifrag, natom, atoms)); for (iatom = 0; iatom < natom; iatom++, itotal++) @@ -204,7 +207,7 @@ void compute_energy(struct state *state, bool do_grad) for (ifrag = 0, itotal = 0, grad = state->grad; ifrag < nfrag; ifrag++, grad += 6) { check_fail(efp_get_frag_xyzabc(state->efp, ifrag, xyzabc)); check_fail(efp_get_frag_atom_count(state->efp, ifrag, &natom)); - atoms = xmalloc(natom * sizeof(struct efp_atom)); + atoms = xcalloc(natom, sizeof(struct efp_atom)); check_fail(efp_get_frag_atoms(state->efp, ifrag, natom, atoms)); for (iatom = 0; iatom < natom; iatom++, itotal++) { diff --git a/efpmd/src/gtest.c b/efpmd/src/gtest.c index 403b7246..fb102b4b 100644 --- a/efpmd/src/gtest.c +++ b/efpmd/src/gtest.c @@ -70,10 +70,10 @@ static void test_cgrad(struct state *state, const double *cgrad) check_fail(efp_get_point_charge_coordinates(state->efp, xyz)); for (size_t i = 0; i < n_charges; i++) { - double ngrad[3]; + double ngrad[3] = {0.0}; for (size_t j = 0; j < 3; j++) { - double e1, e2; + double e1 = 0.0, e2 = 0.0; double coord = xyz[3 * i + j]; xyz[3 * i + j] = coord - dstep; @@ -105,20 +105,23 @@ static void test_fgrad(struct state *state, const double *fgrad) size_t n_frags, spec_frag; check_fail(efp_get_frag_count(state->efp, &n_frags)); - spec_frag = n_frags + 1; // make in inactive if torch model is off + spec_frag = n_frags + 1; // make it inactive if torch model is off if (cfg_get_bool(state->cfg, "enable_torch") ) - spec_frag = cfg_get_int(state->cfg, "special_fragment"); + if (cfg_get_int(state->cfg, "special_fragment") < 0) + error("Special fragment is not defined in ML model"); + spec_frag = (size_t)cfg_get_int(state->cfg, "special_fragment"); + double xyzabc[6 * n_frags]; check_fail(efp_get_coordinates(state->efp, xyzabc)); for (size_t i = 0, k=0; i < n_frags; i++) { if (i == spec_frag) continue; - double deriv[3], ngrad[6]; + double deriv[3] = {0.0}, ngrad[6] = {0.0}; for (size_t j = 0; j < 6; j++) { - double e1, e2; + double e1 = 0.0, e2 = 0.0; double coord = xyzabc[6 * i + j]; double step = j < 3 ? dstep : astep; @@ -155,17 +158,19 @@ static void test_agrad(struct state *state, const double *agrad) size_t spec_frag, n_special_atoms; - spec_frag = cfg_get_int(state->cfg, "special_fragment"); + if (cfg_get_int(state->cfg, "special_fragment") < 0) + error("Special fragment is not defined in ML model!"); + spec_frag = (size_t)cfg_get_int(state->cfg, "special_fragment"); check_fail(efp_get_frag_atom_count(state->efp, spec_frag, &n_special_atoms)); double atom_coord[3 * n_special_atoms]; // = (double*)malloc(3 * n_special_atoms * sizeof(double)); check_fail(efp_get_frag_atom_coord(state->efp, spec_frag, atom_coord)); for (size_t i = 0; i < n_special_atoms; i++) { - double ngrad[3]; + double ngrad[3] = {0.0}; for (size_t j = 0; j < 3; j++) { - double e1, e2; + double e1 = 0.0, e2 = 0.0; double coord = atom_coord[3 * i + j]; atom_coord[3 * i + j] = coord - dstep; @@ -221,12 +226,12 @@ static void test_agrad(struct state *state, const double *agrad) spec_frag = cfg_get_int(state->cfg, "special_fragment"); check_fail(efp_get_frag_atom_count(state->efp, spec_frag, &n_special_atoms)); - double atom_coord[3 * n_special_atoms]; // = (double*)malloc(3 * n_special_atoms * sizeof(double)); + double atom_coord[3 * n_special_atoms] = {0.0}; // = (double*)malloc(3 * n_special_atoms * sizeof(double)); check_fail(efp_get_frag_atom_coord(state->efp, spec_frag, atom_coord)); for (size_t i = 0; i < n_special_atoms; i++) { - double ngrad[3]; - double tmp_agrad[3]; + double ngrad[3] = {0.0}; + double tmp_agrad[3] = {0.0}; for (size_t j = 0; j < 3; j++) { double e1, e2, e3, e4; @@ -300,7 +305,10 @@ static void test_grad(struct state *state) #ifdef TORCH_SWITCH // models with libtorch optimized fragment if (cfg_get_bool(state->cfg, "enable_torch")) { - spec_frag = cfg_get_int(state->cfg, "special_fragment"); + if (cfg_get_int(state->cfg, "special_fragment") < 0) + error("Special fragment is not defined in ML model"); + + spec_frag = (size_t)cfg_get_int(state->cfg, "special_fragment"); check_fail(efp_get_frag_atom_count(state->efp, spec_frag, &n_special_atoms)); // check efp and charges gradient on non-special fragment diff --git a/efpmd/src/main.c b/efpmd/src/main.c index 5a1ccd4f..f7e3f0c3 100644 --- a/efpmd/src/main.c +++ b/efpmd/src/main.c @@ -513,17 +513,17 @@ static void state_init(struct state *state, const struct cfg *cfg, const struct if (state->torch->nn_type != 3) load_ani_model(state->torch->ani_model, state->torch->nn_type, cfg_get_string(state->cfg, "ml_path")); if (state->torch->nn_type == 3) load_custom_ani_model(state->torch->ani_model, state->torch->aev, state->torch->custom_model, ml_location); - spec_frag = cfg_get_int(cfg, "special_fragment"); + spec_frag = (size_t)cfg_get_int(cfg, "special_fragment"); check_fail(efp_get_frag_atom_count(state->efp, spec_frag, &n_special_atoms)); torch_init(state->torch, n_special_atoms); state->torch_grad = xcalloc(n_special_atoms * 3, sizeof(double)); // special fragment atomic coordinates - double *atom_coord = (double*)malloc(3 * n_special_atoms * sizeof(double)); + double *atom_coord = (double*)calloc(3 * n_special_atoms, sizeof(double)); check_fail(efp_get_frag_atom_coord(state->efp, spec_frag, atom_coord)); - int *atom_znuc = (int*)malloc(3 * n_special_atoms * sizeof(int)); + int *atom_znuc = (int*)calloc(3 * n_special_atoms, sizeof(int)); check_fail(efp_get_frag_atom_znuc(state->efp, spec_frag, atom_znuc)); torch_set_coord(state->torch, atom_coord); diff --git a/efpmd/src/opt.c b/efpmd/src/opt.c index ad4704ad..08b6ced6 100644 --- a/efpmd/src/opt.c +++ b/efpmd/src/opt.c @@ -48,7 +48,7 @@ static double compute_efp(size_t n, const double *x, double *gx, void *data) if (cfg_get_bool(state->cfg, "enable_torch") && cfg_get_int(state->cfg, "opt_special_frag") > -1) { // prepare for optimization of atom coordinates of a special fragment // through forces provided externally - spec_frag = cfg_get_int(state->cfg, "special_fragment"); + spec_frag = (size_t)cfg_get_int(state->cfg, "special_fragment"); check_fail(efp_get_frag_atom_count(state->efp, spec_frag, &n_special_atoms)); switch(cfg_get_int(state->cfg, "opt_special_frag")) { @@ -238,8 +238,8 @@ static void print_status(struct state *state, double e_diff, double rms_grad, do void sim_opt(struct state *state) { - size_t n_frags, n_charge, n_coord, n_special_atoms, spec_frag; - double rms_grad, max_grad; + size_t n_frags=0, n_charge=0, n_coord=0, n_special_atoms=0, spec_frag; + double rms_grad = 0.0, max_grad = 0.0; int static algorithm_switch; @@ -275,7 +275,7 @@ void sim_opt(struct state *state) #ifdef TORCH_SWITCH if (cfg_get_bool(state->cfg, "enable_torch") && cfg_get_int(state->cfg, "opt_special_frag") > -1) { - spec_frag = cfg_get_int(state->cfg, "special_fragment"); + spec_frag = (size_t)cfg_get_int(state->cfg, "special_fragment"); check_fail(efp_get_frag_atom_count(state->efp, spec_frag, &n_special_atoms)); if (cfg_get_int(state->cfg, "opt_special_frag") == 0) diff --git a/efpmd/src/parse.c b/efpmd/src/parse.c index 348ef358..6b11b4e5 100644 --- a/efpmd/src/parse.c +++ b/efpmd/src/parse.c @@ -127,7 +127,7 @@ static void parse_frag(struct stream *stream, enum efp_coord_type coord_type, } } else { - frag->coord = (double*)malloc(n_rows * n_cols * sizeof(double)); + frag->coord = (double*)calloc(n_rows * n_cols, sizeof(double)); for (int i = 0, idx = 0; i < n_rows; i++) { for (int j = 0; j < n_cols; j++, idx++) { if (!efp_stream_parse_double(stream, frag->coord + idx)) @@ -203,6 +203,8 @@ struct sys *parse_input(struct cfg *cfg, const char *path) if (is_keyword(efp_stream_get_ptr(stream), "fragment")) { struct frag frag; + memset(&frag, 0, sizeof(frag)); + enum efp_coord_type coord_type; efp_stream_advance(stream, strlen("fragment")); diff --git a/efpmd/src/sp.c b/efpmd/src/sp.c index db540c33..0882c6de 100644 --- a/efpmd/src/sp.c +++ b/efpmd/src/sp.c @@ -33,11 +33,15 @@ void sim_sp(struct state *state) msg("SINGLE POINT ENERGY JOB\n\n\n"); print_geometry(state->efp); - if (cfg_get_bool(state->cfg, "print_pbc") && cfg_get_bool(state->cfg, "enable_pbc")) - print_geometry_pbc(state->efp,cfg_get_int(state->cfg, "ligand")); + if (cfg_get_bool(state->cfg, "print_pbc") && cfg_get_bool(state->cfg, "enable_pbc") ) { + if (cfg_get_int(state->cfg, "ligand") > -1) + print_geometry_pbc(state->efp,cfg_get_int(state->cfg, "ligand")); + else + error("Ligand is not specified for printing PBC geometry"); + } compute_energy(state, false); - if (cfg_get_bool(state->cfg, "enable_pairwise") && cfg_get_int(state->cfg, "ligand") != -1){ + if (cfg_get_bool(state->cfg, "enable_pairwise") && cfg_get_int(state->cfg, "ligand") > -1){ print_pair_energy(state); } print_energy(state); diff --git a/efpmd/src/torch.c b/efpmd/src/torch.c index 2fe3b6e6..7a592c41 100644 --- a/efpmd/src/torch.c +++ b/efpmd/src/torch.c @@ -53,10 +53,10 @@ void get_torch_type(struct torch *torch, const char *str) { void torch_init(struct torch *torch, size_t natom) { torch->natoms = natom; - torch->atom_coords = malloc(3*natom*sizeof(double)); - torch->atom_types = malloc(natom*sizeof(int)); - torch->grad = malloc(3*natom*sizeof(double)); - torch->elpot = malloc(natom*sizeof(double)); + torch->atom_coords = calloc(3*natom, sizeof(double)); + torch->atom_types = calloc(natom, sizeof(int)); + torch->grad = calloc(3*natom, sizeof(double)); + torch->elpot = calloc(natom, sizeof(double)); } void torch_get_atom_count(struct torch *torch , size_t natom) { @@ -101,13 +101,13 @@ void torch_custom_compute(struct torch *torch, int print) { float *frag_coord; float *elecpots_data; //float custom_energy; - double custom_energy; + double custom_energy = 0.0; - elecpots_data = malloc(n_atoms * sizeof(float)); - gradients = malloc(n_atoms * 3 * sizeof(float)); - forces = malloc(n_atoms * 3 * sizeof(float)); + elecpots_data = calloc(n_atoms, sizeof(float)); + gradients = calloc(n_atoms * 3, sizeof(float)); + forces = calloc(n_atoms * 3, sizeof(float)); - frag_coord = malloc(n_atoms*3* sizeof(float)); + frag_coord = calloc(n_atoms*3, sizeof(float)); for (size_t i=0; iatom_coords[i*3] * BOHR_RADIUS); @@ -208,12 +208,12 @@ void torch_compute(struct torch *torch, const char* nn_path, int print) { size_t n_atoms = torch->natoms; float *gradients, *forces, *frag_coord; - double ani_energy; + double ani_energy = 0.0; - gradients = malloc(n_atoms * 3 * sizeof(float)); - forces = malloc(n_atoms * 3 * sizeof(float)); + gradients = calloc(n_atoms * 3, sizeof(float)); + forces = calloc(n_atoms * 3, sizeof(float)); - frag_coord = malloc(n_atoms*3* sizeof(float)); + frag_coord = calloc(n_atoms*3, sizeof(float)); for (size_t i=0; iatom_coords[i*3] * BOHR_RADIUS); frag_coord[i*3+1] = (float)(torch->atom_coords[i*3+1] * BOHR_RADIUS); diff --git a/lib/pylibefp/__pycache__/__init__.cpython-312.pyc b/lib/pylibefp/__pycache__/__init__.cpython-312.pyc index 9cdbc471..b376489a 100644 Binary files a/lib/pylibefp/__pycache__/__init__.cpython-312.pyc and b/lib/pylibefp/__pycache__/__init__.cpython-312.pyc differ diff --git a/lib/pylibefp/__pycache__/wrapper.cpython-312.pyc b/lib/pylibefp/__pycache__/wrapper.cpython-312.pyc index d4eec368..871c3d7f 100644 Binary files a/lib/pylibefp/__pycache__/wrapper.cpython-312.pyc and b/lib/pylibefp/__pycache__/wrapper.cpython-312.pyc differ diff --git a/lib/pylibefp/tests/conftest.py b/lib/pylibefp/tests/conftest.py index c5875acb..9ef25c1e 100644 --- a/lib/pylibefp/tests/conftest.py +++ b/lib/pylibefp/tests/conftest.py @@ -1,4 +1,11 @@ import pytest +from libefp2py import read_libefp_input + + +@pytest.fixture +def pyjob_prepper(): + """Converts efpmd input into py format.""" + return read_libefp_input @pytest.fixture(scope="session", autouse=True) diff --git a/lib/pylibefp/tests/test_coverage.py b/lib/pylibefp/tests/test_coverage.py index 97c80b63..82c15541 100644 --- a/lib/pylibefp/tests/test_coverage.py +++ b/lib/pylibefp/tests/test_coverage.py @@ -6,12 +6,12 @@ from qcelemental.testing import compare -def test_grad_fail(): - asdf = system_1() - asdf.compute(do_gradient=False) - - with pytest.raises(pylibefp.Fatal) as e_info: - grad = asdf.get_gradient() +#def test_grad_fail(): +# asdf = system_1() +# asdf.compute(do_gradient=False) +# +# with pytest.raises(pylibefp.Fatal) as e_info: +# grad = asdf.get_gradient() #def test_frag_file_fail(): diff --git a/lib/pylibefp/tests/test_efpefp.py b/lib/pylibefp/tests/test_efpefp.py index ab87085c..fa276a89 100644 --- a/lib/pylibefp/tests/test_efpefp.py +++ b/lib/pylibefp/tests/test_efpefp.py @@ -118,7 +118,7 @@ def test_total_1a(): 'pol': True, # 'pol_damp': 'tt', 'disp': True, 'disp_damp': 'tt', - 'print': 2 + 'print': 1 }) asdf.compute() ene = asdf.get_energy() diff --git a/lib/pylibefp/tests/test_efpefp_new.py b/lib/pylibefp/tests/test_efpefp_new.py deleted file mode 100644 index 44284c98..00000000 --- a/lib/pylibefp/tests/test_efpefp_new.py +++ /dev/null @@ -1,57 +0,0 @@ -import libefp2py -import pylibefp -from qcelemental.testing import compare, compare_values -import pprint - - -b2a = 0.529177 -a2b = 1.0 / b2a - -def frag_setup(test_name): - # coordinates in Bohr - coord_type, frags, frag_coords, efp_options, if_gradient, ref_energy, periodic_box = libefp2py.read_libefp_input(test_name) - #print(frag_coords) - - efp = pylibefp.core.efp() - efp.add_potential(frags) - efp.add_fragment(frags) - for i in range(len(frags)): - efp.set_frag_coordinates(i, coord_type, frag_coords[i]) - efp.prepare() - - efp.set_opts(efp_options) - if periodic_box: - #print('box1', periodic_box) - efp.set_periodic_box(periodic_box) - #print('box2', efp.get_periodic_box()) - - #print(frag_coords) - #pprint.pprint(efp_options) - efp.compute(do_gradient = if_gradient) - ene = efp.get_energy() - - # print pairwise components - #if 'enable_pairwise' in efp_options.keys(): - # if efp_options['enable_pairwise'] in [True, 'true', 1]: - # efp.print_pairwise_energies() - - print(efp.energy_summary()) - if if_gradient: - print(efp.gradient_summary()) - if ref_energy != 0.0: - assert compare_values(ref_energy, ene['total'], atol=1.e-5, return_message=True), 'FAILED' - -##### -if __name__ == '__main__': - files = ['atom_coord.in', 'atom_coord_2.in', 'grad_1.in', 'lj_1.in', 'lj_2.in', - 'pairwise_0.in', 'pairwise_1.in', 'pairwise_2.in', 'pairwise_x.in', 'pbc_1.in', 'pbc_2.in', - 'reduced.in', 'spec_frag_0.in', 'spec_frag_1.in', 'spec_frag_2.in', 'spec_frag_base.in', 'spec_frag_ref.in', - 'symm_1.in', 'symm_2.in', 'symm_2full.in', 'symm_2pw.in'] - - # running for all tests in files list - for f in files: - print(f'\nComputing {f}...') - frag_setup('../'+f) - - # single test execution - frag_setup('../symm_2pw.in') diff --git a/lib/pylibefp/tests/test_lori.py b/lib/pylibefp/tests/test_lori.py deleted file mode 100644 index 27f98e97..00000000 --- a/lib/pylibefp/tests/test_lori.py +++ /dev/null @@ -1,26 +0,0 @@ -import pylibefp - - - -efp = pylibefp.core.efp() - -frags = ["h2o_l", "nh3_l"] -efp.add_potential(frags) -efp.add_fragment(frags) -efp.set_frag_coordinates(0, "xyzabc", - [0.0, 0.0, 0.0, 1.0, 2.0, 3.0]) -efp.set_frag_coordinates(1, "xyzabc", - [9.0, 0.0, 0.0, 5.0, 2.0, 8.0]) -efp.prepare() - -efp.set_opts({ - "elec": True, - "elec_damp": "screen", - "xr": True, - "pol": True, - "disp": False, -}) - -efp.compute() -ene = efp.get_energy() -print(ene) diff --git a/lib/pylibefp/wrapper.py b/lib/pylibefp/wrapper.py index 7325a5c9..6921ff67 100644 --- a/lib/pylibefp/wrapper.py +++ b/lib/pylibefp/wrapper.py @@ -918,109 +918,109 @@ def get_frag_count(efpobj): return nfrag -# # -# # def get_multipole_count(efpobj): -# # """Gets the number of multipoles in `efpobj` computation. -# # -# # Returns -# # ------- -# # int -# # Total number of multipoles from electrostatics calculation. -# # -# # """ -# # (res, nmult) = efpobj._efp_get_multipole_count() -# # _result_to_error(res) -# # -# # return nmult -# # -# -# def get_multipole_coordinates(efpobj, verbose=1): -# """Gets the coordinates of `efpobj` electrostatics multipoles. -# -# Parameters -# ---------- -# verbose : int, optional -# Whether to print out the multipole coordinates. 0: no printing. 1: -# print charges and dipoles. 2: additionally print quadrupoles -# and octupoles. -# -# Returns -# ------- -# list -# ``3 n_mult`` (flat) array of multipole locations. -# -# Examples -# -------- -# -# >>> # Use with NumPy -# >>> n_mult = efpobj.get_multipole_count() -# >>> xyz_mult = np.asarray(efpobj.get_multipole_coordinates()).reshape(n_mult, 3) -# -# """ -# nmult = efpobj.get_multipole_count() -# (res, xyz) = efpobj._efp_get_multipole_coordinates(nmult) -# _result_to_error(res) -# -# if verbose >= 1: -# xyz3 = list(map(list, zip(*[iter(xyz)] * 3))) -# -# text = '\n ==> EFP Multipole Coordinates <==\n\n' -# for mu in range(nmult): -# text += '{:6d} {:14.8f} {:14.8f} {:14.8f}\n'.format(mu, *xyz3[mu]) -# print(text) -# -# return xyz -# -# -# def get_multipole_values(efpobj, verbose=1): -# """Gets the computed per-point multipoles of `efpobj`. -# -# Parameters -# ---------- -# verbose : int, optional -# Whether to print out the multipole arrays. 0: no printing. 1: -# print charges and dipoles. ``2``: additionally print quadrupoles -# and octupoles. -# -# Returns -# ------- -# list -# ``20 n_mult`` (flat) array of per-point multipole values including -# charges + dipoles + quadrupoles + octupoles. -# Dipoles stored as x, y, z. -# Quadrupoles stored as xx, yy, zz, xy, xz, yz. -# Octupoles stored as xxx, yyy, zzz, xxy, xxz, xyy, yyz, xzz, yzz, xyz. -# -# Examples -# -------- -# >>> # Use with NumPy -# >>> n_mult = efpobj.get_multipole_count() -# >>> val_mult = np.asarray(efpobj.get_multipole_values()).reshape(n_mult, 20) -# -# """ -# nmult = efpobj.get_multipole_count() -# (res, mult) = efpobj._efp_get_multipole_values(nmult) -# _result_to_error(res) -# -# if verbose >= 1: -# mult20 = list(map(list, zip(*[iter(mult)] * 20))) -# -# text = '\n ==> EFP Multipoles: Charge & Dipole <==\n\n' -# for mu in range(nmult): -# text += '{:6d} {:14.8f} {:14.8f} {:14.8f} {:14.8f}\n'.format(mu, *mult20[mu][:4]) -# -# if verbose >= 2: -# text += '\n ==> EFP Multipoles: Quadrupole <==\n\n' -# for mu in range(nmult): -# text += '{:6d} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f}\n'.format(mu, *mult20[mu][4:10]) -# text += '\n ==> EFP Multipoles: Octupole <==\n\n' -# for mu in range(nmult): -# text += '{:6d} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f}\n'.format( -# mu, *mult20[mu][10:]) -# print(text) -# -# return mult -# + +def get_multipole_count(efpobj): + """Gets the number of multipoles in `efpobj` computation. + + Returns + ------- + int + Total number of multipoles from electrostatics calculation. + + """ + (res, nmult) = efpobj._efp_get_multipole_count() + _result_to_error(res) + + return nmult + + +def get_multipole_coordinates(efpobj, verbose=1): + """Gets the coordinates of `efpobj` electrostatics multipoles. + + Parameters + ---------- + verbose : int, optional + Whether to print out the multipole coordinates. 0: no printing. 1: + print charges and dipoles. 2: additionally print quadrupoles + and octupoles. + + Returns + ------- + list + ``3 n_mult`` (flat) array of multipole locations. + + Examples + -------- + + >>> # Use with NumPy + >>> n_mult = efpobj.get_multipole_count() + >>> xyz_mult = np.asarray(efpobj.get_multipole_coordinates()).reshape(n_mult, 3) + + """ + nmult = efpobj.get_multipole_count() + (res, xyz) = efpobj._efp_get_multipole_coordinates(nmult) + _result_to_error(res) + + if verbose >= 1: + xyz3 = list(map(list, zip(*[iter(xyz)] * 3))) + + text = '\n ==> EFP Multipole Coordinates <==\n\n' + for mu in range(nmult): + text += '{:6d} {:14.8f} {:14.8f} {:14.8f}\n'.format(mu, *xyz3[mu]) + print(text) + + return xyz + + +def get_multipole_values(efpobj, verbose=1): + """Gets the computed per-point multipoles of `efpobj`. + + Parameters + ---------- + verbose : int, optional + Whether to print out the multipole arrays. 0: no printing. 1: + print charges and dipoles. ``2``: additionally print quadrupoles + and octupoles. + + Returns + ------- + list + ``20 n_mult`` (flat) array of per-point multipole values including + charges + dipoles + quadrupoles + octupoles. + Dipoles stored as x, y, z. + Quadrupoles stored as xx, yy, zz, xy, xz, yz. + Octupoles stored as xxx, yyy, zzz, xxy, xxz, xyy, yyz, xzz, yzz, xyz. + + Examples + -------- + >>> # Use with NumPy + >>> n_mult = efpobj.get_multipole_count() + >>> val_mult = np.asarray(efpobj.get_multipole_values()).reshape(n_mult, 20) + + """ + nmult = efpobj.get_multipole_count() + (res, mult) = efpobj._efp_get_multipole_values(nmult) + _result_to_error(res) + + if verbose >= 1: + mult20 = list(map(list, zip(*[iter(mult)] * 20))) + + text = '\n ==> EFP Multipoles: Charge & Dipole <==\n\n' + for mu in range(nmult): + text += '{:6d} {:14.8f} {:14.8f} {:14.8f} {:14.8f}\n'.format(mu, *mult20[mu][:4]) + + if verbose >= 2: + text += '\n ==> EFP Multipoles: Quadrupole <==\n\n' + for mu in range(nmult): + text += '{:6d} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f}\n'.format(mu, *mult20[mu][4:10]) + text += '\n ==> EFP Multipoles: Octupole <==\n\n' + for mu in range(nmult): + text += '{:6d} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f}\n'.format( + mu, *mult20[mu][10:]) + print(text) + + return mult + def get_induced_dipole_count(efpobj): """Gets the number of polarization induced dipoles in `efpobj` computation. @@ -1930,9 +1930,9 @@ def old_to_dict(efpobj): core.efp.get_point_charge_count = get_point_charge_count core.efp.get_point_charge_coordinates = get_point_charge_coordinates core.efp.get_point_charge_values = get_point_charge_values -#core.efp.get_multipole_count = get_multipole_count -#core.efp.get_multipole_coordinates = get_multipole_coordinates -#core.efp.get_multipole_values = get_multipole_values +core.efp.get_multipole_count = get_multipole_count +core.efp.get_multipole_coordinates = get_multipole_coordinates +core.efp.get_multipole_values = get_multipole_values core.efp.get_induced_dipole_count = get_induced_dipole_count core.efp.get_induced_dipole_coordinates = get_induced_dipole_coordinates core.efp.get_induced_dipole_values = get_induced_dipole_values diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 97525d64..4ee3e825 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -166,7 +166,12 @@ install( PATTERN "systems.py" PATTERN "__pycache__" EXCLUDE PATTERN "*pyc" EXCLUDE - PATTERN "test_libefp.py" EXCLUDE # until xr_cutoff resolved - PATTERN "test_opts.py" EXCLUDE # until xr_cutoff resolved + PATTERN "test_dict.py" # EXCLUDE + PATTERN "test_libefp.py" # EXCLUDE # until xr_cutoff resolved + PATTERN "test_opts.py" # EXCLUDE # until xr_cutoff resolved + PATTERN "test_efpefp.py" + PATTERN "test_efpefp2.py" EXCLUDE + PATTERN "libefp2py.py" + PATTERN "test_psi.py" #EXCLUDE + PATTERN "test_scf.py" #EXCLUDE ) - diff --git a/python/core.cc b/python/core.cc index 1cb806b9..a6e37993 100644 --- a/python/core.cc +++ b/python/core.cc @@ -227,7 +227,7 @@ py::tuple _efp_get_periodic_box(efp* efp) { return rets; } -/* + py::tuple _efp_get_multipole_count(efp* efp) { enum efp_result res; size_t n_mult = 0; @@ -269,7 +269,7 @@ py::tuple _efp_get_multipole_values(efp* efp, size_t n_mult) { py::tuple rets = py::make_tuple(res, mult); return rets; } -*/ + py::tuple _efp_get_induced_dipole_count(efp* efp) { enum efp_result res; size_t n_dip = 0; @@ -763,12 +763,12 @@ PYBIND11_MODULE(core, m) { .def("_efp_get_frag_charge", &_efp_get_frag_charge, "Gets total charge on fragment", py::arg("frag_idx")) .def("_efp_get_frag_multiplicity", &_efp_get_frag_multiplicity, "Gets spin multiplicity on fragment") // Multipoles & Induced Dipoles - //.def("_efp_get_multipole_count", &_efp_get_multipole_count, - // "Wrapped gets total number of multipoles from EFP electrostatics") - //.def("_efp_get_multipole_coordinates", &_efp_get_multipole_coordinates, - // "Wrapped gets coordinates of electrostatics multipoles") - //.def("_efp_get_multipole_values", &_efp_get_multipole_values, - // "Wrapped gets electrostatics multipoles from EFP fragments") + .def("_efp_get_multipole_count", &_efp_get_multipole_count, + "Wrapped gets total number of multipoles from EFP electrostatics") + .def("_efp_get_multipole_coordinates", &_efp_get_multipole_coordinates, + "Wrapped gets coordinates of electrostatics multipoles") + .def("_efp_get_multipole_values", &_efp_get_multipole_values, + "Wrapped gets electrostatics multipoles from EFP fragments") .def("_efp_get_induced_dipole_count", &_efp_get_induced_dipole_count, "Wrapped gets the number of polarization induced dipoles") .def("_efp_get_induced_dipole_coordinates", &_efp_get_induced_dipole_coordinates, diff --git a/python/wrapper.py b/python/wrapper.py index b6aad04c..1f2c671d 100644 --- a/python/wrapper.py +++ b/python/wrapper.py @@ -918,109 +918,109 @@ def get_frag_count(efpobj): return nfrag -# # -# # def get_multipole_count(efpobj): -# # """Gets the number of multipoles in `efpobj` computation. -# # -# # Returns -# # ------- -# # int -# # Total number of multipoles from electrostatics calculation. -# # -# # """ -# # (res, nmult) = efpobj._efp_get_multipole_count() -# # _result_to_error(res) -# # -# # return nmult -# # -# -# def get_multipole_coordinates(efpobj, verbose=1): -# """Gets the coordinates of `efpobj` electrostatics multipoles. -# -# Parameters -# ---------- -# verbose : int, optional -# Whether to print out the multipole coordinates. 0: no printing. 1: -# print charges and dipoles. 2: additionally print quadrupoles -# and octupoles. -# -# Returns -# ------- -# list -# ``3 n_mult`` (flat) array of multipole locations. -# -# Examples -# -------- -# -# >>> # Use with NumPy -# >>> n_mult = efpobj.get_multipole_count() -# >>> xyz_mult = np.asarray(efpobj.get_multipole_coordinates()).reshape(n_mult, 3) -# -# """ -# nmult = efpobj.get_multipole_count() -# (res, xyz) = efpobj._efp_get_multipole_coordinates(nmult) -# _result_to_error(res) -# -# if verbose >= 1: -# xyz3 = list(map(list, zip(*[iter(xyz)] * 3))) -# -# text = '\n ==> EFP Multipole Coordinates <==\n\n' -# for mu in range(nmult): -# text += '{:6d} {:14.8f} {:14.8f} {:14.8f}\n'.format(mu, *xyz3[mu]) -# print(text) -# -# return xyz -# -# -# def get_multipole_values(efpobj, verbose=1): -# """Gets the computed per-point multipoles of `efpobj`. -# -# Parameters -# ---------- -# verbose : int, optional -# Whether to print out the multipole arrays. 0: no printing. 1: -# print charges and dipoles. ``2``: additionally print quadrupoles -# and octupoles. -# -# Returns -# ------- -# list -# ``20 n_mult`` (flat) array of per-point multipole values including -# charges + dipoles + quadrupoles + octupoles. -# Dipoles stored as x, y, z. -# Quadrupoles stored as xx, yy, zz, xy, xz, yz. -# Octupoles stored as xxx, yyy, zzz, xxy, xxz, xyy, yyz, xzz, yzz, xyz. -# -# Examples -# -------- -# >>> # Use with NumPy -# >>> n_mult = efpobj.get_multipole_count() -# >>> val_mult = np.asarray(efpobj.get_multipole_values()).reshape(n_mult, 20) -# -# """ -# nmult = efpobj.get_multipole_count() -# (res, mult) = efpobj._efp_get_multipole_values(nmult) -# _result_to_error(res) -# -# if verbose >= 1: -# mult20 = list(map(list, zip(*[iter(mult)] * 20))) -# -# text = '\n ==> EFP Multipoles: Charge & Dipole <==\n\n' -# for mu in range(nmult): -# text += '{:6d} {:14.8f} {:14.8f} {:14.8f} {:14.8f}\n'.format(mu, *mult20[mu][:4]) -# -# if verbose >= 2: -# text += '\n ==> EFP Multipoles: Quadrupole <==\n\n' -# for mu in range(nmult): -# text += '{:6d} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f}\n'.format(mu, *mult20[mu][4:10]) -# text += '\n ==> EFP Multipoles: Octupole <==\n\n' -# for mu in range(nmult): -# text += '{:6d} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f}\n'.format( -# mu, *mult20[mu][10:]) -# print(text) -# -# return mult -# + +def get_multipole_count(efpobj): + """Gets the number of multipoles in `efpobj` computation. + + Returns + ------- + int + Total number of multipoles from electrostatics calculation. + + """ + (res, nmult) = efpobj._efp_get_multipole_count() + _result_to_error(res) + + return nmult + + +def get_multipole_coordinates(efpobj, verbose=1): + """Gets the coordinates of `efpobj` electrostatics multipoles. + + Parameters + ---------- + verbose : int, optional + Whether to print out the multipole coordinates. 0: no printing. 1: + print charges and dipoles. 2: additionally print quadrupoles + and octupoles. + + Returns + ------- + list + ``3 n_mult`` (flat) array of multipole locations. + + Examples + -------- + + >>> # Use with NumPy + >>> n_mult = efpobj.get_multipole_count() + >>> xyz_mult = np.asarray(efpobj.get_multipole_coordinates()).reshape(n_mult, 3) + + """ + nmult = efpobj.get_multipole_count() + (res, xyz) = efpobj._efp_get_multipole_coordinates(nmult) + _result_to_error(res) + + if verbose >= 1: + xyz3 = list(map(list, zip(*[iter(xyz)] * 3))) + + text = '\n ==> EFP Multipole Coordinates <==\n\n' + for mu in range(nmult): + text += '{:6d} {:14.8f} {:14.8f} {:14.8f}\n'.format(mu, *xyz3[mu]) + print(text) + + return xyz + + +def get_multipole_values(efpobj, verbose=1): + """Gets the computed per-point multipoles of `efpobj`. + + Parameters + ---------- + verbose : int, optional + Whether to print out the multipole arrays. 0: no printing. 1: + print charges and dipoles. ``2``: additionally print quadrupoles + and octupoles. + + Returns + ------- + list + ``20 n_mult`` (flat) array of per-point multipole values including + charges + dipoles + quadrupoles + octupoles. + Dipoles stored as x, y, z. + Quadrupoles stored as xx, yy, zz, xy, xz, yz. + Octupoles stored as xxx, yyy, zzz, xxy, xxz, xyy, yyz, xzz, yzz, xyz. + + Examples + -------- + >>> # Use with NumPy + >>> n_mult = efpobj.get_multipole_count() + >>> val_mult = np.asarray(efpobj.get_multipole_values()).reshape(n_mult, 20) + + """ + nmult = efpobj.get_multipole_count() + (res, mult) = efpobj._efp_get_multipole_values(nmult) + _result_to_error(res) + + if verbose >= 1: + mult20 = list(map(list, zip(*[iter(mult)] * 20))) + + text = '\n ==> EFP Multipoles: Charge & Dipole <==\n\n' + for mu in range(nmult): + text += '{:6d} {:14.8f} {:14.8f} {:14.8f} {:14.8f}\n'.format(mu, *mult20[mu][:4]) + + if verbose >= 2: + text += '\n ==> EFP Multipoles: Quadrupole <==\n\n' + for mu in range(nmult): + text += '{:6d} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f}\n'.format(mu, *mult20[mu][4:10]) + text += '\n ==> EFP Multipoles: Octupole <==\n\n' + for mu in range(nmult): + text += '{:6d} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f}\n'.format( + mu, *mult20[mu][10:]) + print(text) + + return mult + def get_induced_dipole_count(efpobj): """Gets the number of polarization induced dipoles in `efpobj` computation. @@ -1930,9 +1930,9 @@ def old_to_dict(efpobj): core.efp.get_point_charge_count = get_point_charge_count core.efp.get_point_charge_coordinates = get_point_charge_coordinates core.efp.get_point_charge_values = get_point_charge_values -#core.efp.get_multipole_count = get_multipole_count -#core.efp.get_multipole_coordinates = get_multipole_coordinates -#core.efp.get_multipole_values = get_multipole_values +core.efp.get_multipole_count = get_multipole_count +core.efp.get_multipole_coordinates = get_multipole_coordinates +core.efp.get_multipole_values = get_multipole_values core.efp.get_induced_dipole_count = get_induced_dipole_count core.efp.get_induced_dipole_coordinates = get_induced_dipole_coordinates core.efp.get_induced_dipole_values = get_induced_dipole_values diff --git a/setup.sh b/setup.sh index 16dca756..c32afcff 100644 --- a/setup.sh +++ b/setup.sh @@ -19,9 +19,9 @@ if [[ "$TORCH_SWITCH" == "ON" ]]; then echo "TORCHANI_DIR=$TORCHANI_DIR" echo "PYTHON_REQS=$PYTHON_REQS" else - unsetenv LIBTORCH_INCLUDE_DIRS - unsetenv TORCH_INSTALLED_DIR - unsetenv TORCHANI_DIR + unset LIBTORCH_INCLUDE_DIRS + unset TORCH_INSTALLED_DIR + unset TORCHANI_DIR echo "Torch integration is disabled. Only basic environment variables are set:" echo "LIBEFP_DIR=$LIBEFP_DIR" diff --git a/src/efp.c b/src/efp.c index a3d95f35..144ff0a9 100644 --- a/src/efp.c +++ b/src/efp.c @@ -341,36 +341,6 @@ copy_frag(struct frag *dest, const struct frag *src) return EFP_RESULT_SUCCESS; } -/* -static enum efp_result -copy_ligand(struct ligand *dest, const struct ligand *src) { - size_t size; - - memcpy(dest, src, sizeof(*dest)); - - if (src->ligand_frag) - copy_frag(dest->ligand_frag, src->ligand_frag); - - if (src->ligand_pts) { - size = src->n_ligand_pts * sizeof(struct ligand_pt); - dest->ligand_pts = (struct ligand_pt *) malloc(size); - if (!dest->ligand_pts) - return EFP_RESULT_NO_MEMORY; - memcpy(dest->ligand_pts, src->ligand_pts, size); - - for (size_t i=0; in_ligand_pts; i++) { - const struct ligand_pt *src_pt = src->ligand_pts + i; - struct ligand_pt *dest_pt = dest->ligand_pts + i; - size = src_pt->n_frag * sizeof(vec_t); - dest_pt->fragment_field = (vec_t *)malloc(size); - if (!dest_pt->fragment_field) - return EFP_RESULT_NO_MEMORY; - memcpy(dest_pt->fragment_field, src_pt->fragment_field, size); - } - } -} -*/ - // updates (shifts) parameters of fragment based on coordinates of fragment atoms static enum efp_result update_params(struct efp_atom *atoms, const struct frag *lib_orig, const struct frag *lib_current) { @@ -546,8 +516,8 @@ static void compute_two_body_range(struct efp *efp, size_t frag_from, size_t frag_to, void *data) { - double e_elec = 0.0, e_disp = 0.0, e_xr = 0.0, e_cp = 0.0, e_elec_tmp = 0.0, e_disp_tmp = 0.0; - double e_lj = 0.0, e_qq = 0.0, e_qq_tmp = 0.0; + double e_elec = 0.0, e_disp = 0.0, e_xr = 0.0, e_cp = 0.0; + double e_lj = 0.0, e_qq = 0.0; (void)data; @@ -565,13 +535,19 @@ compute_two_body_range(struct efp *efp, size_t frag_from, size_t frag_to, for (size_t j = i + 1; j < i + 1 + cnt; j++) { size_t fr_j = j % efp->n_frag; + double e_elec_tmp = 0.0, e_disp_tmp = 0.0, e_qq_tmp = 0.0; + // special fragment - additional check on special terms - bool if_special_fragment = efp->opts.special_fragment == i || efp->opts.special_fragment == fr_j; - bool special_xr = if_special_fragment && (efp->opts.special_terms & EFP_SPEC_TERM_XR); - bool special_elec = if_special_fragment && (efp->opts.special_terms & EFP_SPEC_TERM_ELEC); - bool special_disp = if_special_fragment && (efp->opts.special_terms & EFP_SPEC_TERM_DISP); - bool special_qq = (!if_special_fragment) || (if_special_fragment && (efp->opts.special_terms & EFP_SPEC_TERM_QQ)); - bool special_lj = if_special_fragment && (efp->opts.special_terms & EFP_SPEC_TERM_LJ); + bool if_special_fragment=false, special_xr=false, special_elec=false, + special_disp=false, special_qq=false, special_lj=false; + if (efp->opts.special_fragment >=0) { + if_special_fragment = (size_t)efp->opts.special_fragment == i || (size_t)efp->opts.special_fragment == fr_j; + } + special_xr = if_special_fragment && (efp->opts.special_terms & EFP_SPEC_TERM_XR); + special_elec = if_special_fragment && (efp->opts.special_terms & EFP_SPEC_TERM_ELEC); + special_disp = if_special_fragment && (efp->opts.special_terms & EFP_SPEC_TERM_DISP); + special_qq = (!if_special_fragment) || (if_special_fragment && (efp->opts.special_terms & EFP_SPEC_TERM_QQ)); + special_lj = if_special_fragment && (efp->opts.special_terms & EFP_SPEC_TERM_LJ); if (!efp_skip_frag_pair(efp, i, fr_j)) { double *s; @@ -584,8 +560,8 @@ compute_two_body_range(struct efp *efp, size_t frag_from, size_t frag_to, s = (double *) calloc(n_lmo_ij, sizeof(double)); ds = (six_t *) calloc(n_lmo_ij, sizeof(six_t)); - if (do_xr(&efp->opts) || special_xr) { - double exr, ecp; + if ((do_xr(&efp->opts) && !if_special_fragment) || special_xr) { + double exr = 0.0, ecp = 0.0; efp_frag_frag_xr(efp, i, fr_j, s, ds, &exr, &ecp); @@ -605,16 +581,17 @@ compute_two_body_range(struct efp *efp, size_t frag_from, size_t frag_to, } } } - if ((do_elec(&efp->opts) || special_elec) && efp->frags[i].n_multipole_pts > 0 && - efp->frags[fr_j].n_multipole_pts > 0) { + if (((do_elec(&efp->opts) && !if_special_fragment) || special_elec) + && efp->frags[i].n_multipole_pts > 0 && efp->frags[fr_j].n_multipole_pts > 0) { e_elec_tmp = efp_frag_frag_elec(efp, i, fr_j); if (efp->opts.print > 0 && fabs(e_elec_tmp) > 1.0) printf(" WARNING: elec energy between fragments %zu and %zu is %lf \n", i, fr_j, e_elec_tmp); // zeroing the energy contribution on the special fragment in torch custom models - if (efp->opts.enable_elpot && if_special_fragment) e_elec_tmp = 0.0; - //e_elec_tmp = efp_frag_frag_elec(efp, i, fr_j); + if (efp->opts.enable_elpot && if_special_fragment) + e_elec_tmp = 0.0; + e_elec += e_elec_tmp; /* */ if (if_pairwise) { @@ -624,10 +601,9 @@ compute_two_body_range(struct efp *efp, size_t frag_from, size_t frag_to, efp->pair_energies[i].electrostatic = e_elec_tmp; } } - if ((do_disp(&efp->opts) || special_disp) && efp->frags[i].n_dynamic_polarizable_pts > 0 && - efp->frags[fr_j].n_dynamic_polarizable_pts > 0) { - e_disp_tmp = efp_frag_frag_disp(efp, - i, fr_j, s, ds); + if (((do_disp(&efp->opts) && !if_special_fragment) || special_disp) + && efp->frags[i].n_dynamic_polarizable_pts > 0 && efp->frags[fr_j].n_dynamic_polarizable_pts > 0) { + e_disp_tmp = efp_frag_frag_disp(efp, i, fr_j, s, ds); e_disp += e_disp_tmp; /* */ if (if_pairwise) { @@ -638,11 +614,11 @@ compute_two_body_range(struct efp *efp, size_t frag_from, size_t frag_to, } } // LJ terms - if (do_lj(&efp->opts) || special_lj) { + if ((do_lj(&efp->opts) && !if_special_fragment) || special_lj) { e_lj += efp_frag_frag_lj(efp, i, fr_j); } // MM-like charge-charge interactions - if (do_qq(&efp->opts) && special_qq) { + if ((do_qq(&efp->opts) && !if_special_fragment) || special_qq) { e_qq_tmp = efp_frag_frag_qq(efp, i, fr_j); // zeroing the energy contribution on the special fragment in torch custom models @@ -800,7 +776,7 @@ EFP_EXPORT enum efp_result efp_get_atomic_gradient(struct efp *efp, double *grad) { six_t *efpgrad = NULL; /* Calculated EFP gradient */ - vec_t *pgrad; /* Conversion of grad to vec_t type */ + vec_t *pgrad = NULL; /* Conversion of grad to vec_t type */ size_t i, j, k, l; size_t nr; /* Number of atoms in the current fragment */ size_t maxa; /* Maximum number of size of m, Ia, r arrays */ @@ -833,15 +809,15 @@ efp_get_atomic_gradient(struct efp *efp, double *grad) res = EFP_RESULT_NO_MEMORY; /* Create and initialize some arrays for work */ - if ((r = (vec_t *)malloc(maxa * sizeof(*r))) == NULL) + if ((r = (vec_t *)calloc(maxa, sizeof(*r))) == NULL) goto error; - if ((m = (double *)malloc(maxa * sizeof(*m))) == NULL) + if ((m = (double *)calloc(maxa, sizeof(*m))) == NULL) goto error; - if ((Ia = (double *)malloc(maxa * sizeof(*Ia))) == NULL) + if ((Ia = (double *)calloc(maxa, sizeof(*Ia))) == NULL) goto error; /* Copy computed efp->grad */ - if ((efpgrad = (six_t *)malloc(efp->n_frag * sizeof(*efpgrad))) == NULL) + if ((efpgrad = (six_t *)calloc(efp->n_frag, sizeof(*efpgrad))) == NULL) goto error; memcpy(efpgrad, efp->grad, efp->n_frag * sizeof(*efpgrad)); @@ -963,7 +939,7 @@ EFP_EXPORT enum efp_result efp_get_frag_atomic_gradient(struct efp *efp, size_t frag_id, double *grad) { six_t *efpgrad = NULL; /* Calculated EFP gradient */ - vec_t *pgrad; /* Conversion of grad to vec_t type */ + vec_t *pgrad = NULL; /* Conversion of grad to vec_t type */ size_t i, j, k, l; size_t nr; /* Number of atoms in the current fragment */ vec_t *r = NULL; /* Radius-vector of each atom inside current fragment @@ -992,11 +968,11 @@ efp_get_frag_atomic_gradient(struct efp *efp, size_t frag_id, double *grad) res = EFP_RESULT_NO_MEMORY; /* Create and initialize some arrays for work */ - if ((r = (vec_t *)malloc(nr * sizeof(*r))) == NULL) + if ((r = (vec_t *)calloc(nr, sizeof(*r))) == NULL) goto error; - if ((m = (double *)malloc(nr * sizeof(*m))) == NULL) + if ((m = (double *)calloc(nr, sizeof(*m))) == NULL) goto error; - if ((Ia = (double *)malloc(nr * sizeof(*Ia))) == NULL) + if ((Ia = (double *)calloc(nr, sizeof(*Ia))) == NULL) goto error; //for (size_t i=0; in_frag; i++) { @@ -1005,7 +981,7 @@ efp_get_frag_atomic_gradient(struct efp *efp, size_t frag_id, double *grad) /* Copy computed efp->grad */ - if ((efpgrad = (six_t *)malloc(sizeof(*efpgrad))) == NULL) + if ((efpgrad = (six_t *)calloc(1,sizeof(*efpgrad))) == NULL) goto error; memcpy(efpgrad, efp->grad + frag_id, sizeof(*efpgrad)); @@ -1208,7 +1184,11 @@ update_special_fragment(struct efp *efp, const double *coord) assert(efp); assert(coord); //printf("Inside update_special_fragment\n"); - size_t fr_i = efp->opts.special_fragment; + if (efp->opts.special_fragment < 0) { + efp_log("special fragment not set"); + return EFP_RESULT_FATAL; + } + size_t fr_i = (size_t)efp->opts.special_fragment; struct frag *spec_frag = efp->frags + fr_i; if (set_coord_atoms(spec_frag, coord)) { @@ -1224,11 +1204,11 @@ update_gradient_special_fragment(struct efp *efp) { assert(efp); if (efp->opts.special_fragment < 0) { - efp_log("special fragment not set, continue"); - return EFP_RESULT_SUCCESS; + efp_log("special fragment not set"); + return EFP_RESULT_FATAL; } - size_t fr_i = efp->opts.special_fragment; + size_t fr_i = (size_t)efp->opts.special_fragment; struct frag *spec_frag = efp->frags + fr_i; for (size_t i=0; in_atoms; i++) { @@ -1893,7 +1873,7 @@ efp_get_frag_rank(struct efp *efp, size_t frag_idx, int *rank) return EFP_RESULT_SUCCESS; } -/* + EFP_EXPORT enum efp_result efp_get_multipole_count(struct efp *efp, size_t *n_mult) { @@ -1908,7 +1888,7 @@ efp_get_multipole_count(struct efp *efp, size_t *n_mult) *n_mult = sum; return EFP_RESULT_SUCCESS; } -*/ + EFP_EXPORT enum efp_result efp_get_ho_multipole_count(struct efp *efp, size_t *n_mult) @@ -1942,7 +1922,7 @@ efp_get_mm_multipole_count(struct efp *efp, size_t *n_mult) return EFP_RESULT_SUCCESS; } -/* + EFP_EXPORT enum efp_result efp_get_multipole_coordinates(struct efp *efp, double *xyz) { @@ -1960,7 +1940,7 @@ efp_get_multipole_coordinates(struct efp *efp, double *xyz) } return EFP_RESULT_SUCCESS; } -*/ + EFP_EXPORT enum efp_result efp_get_ho_multipole_coordinates(struct efp *efp, double *xyz) @@ -2000,7 +1980,7 @@ efp_get_mm_multipole_coordinates(struct efp *efp, double *xyz) return EFP_RESULT_SUCCESS; } -/* + EFP_EXPORT enum efp_result efp_get_multipole_values(struct efp *efp, double *mult) { @@ -2028,7 +2008,7 @@ efp_get_multipole_values(struct efp *efp, double *mult) } return EFP_RESULT_SUCCESS; } -*/ + EFP_EXPORT enum efp_result efp_get_ho_multipole_values(struct efp *efp, double *mult) @@ -2484,6 +2464,7 @@ efp_add_fragment(struct efp *efp, const char *name) enum efp_result res; struct frag *frag = efp->frags + efp->n_frag - 1; + memset(frag, 0, sizeof(*frag)); // if update/rotate parameters if (efp->opts.update_params == 1) { @@ -2536,17 +2517,14 @@ efp_add_ligand(struct efp *efp, int ligand_index) { lig->ligand_frag = &efp->frags[ligand_index]; lig->n_ligand_pts = efp->frags[ligand_index].n_polarizable_pts; - size_t size; - size = lig->n_ligand_pts * sizeof(struct ligand_pt); - lig->ligand_pts = (struct ligand_pt *) malloc(size); + lig->ligand_pts = (struct ligand_pt *) calloc(lig->n_ligand_pts, sizeof(struct ligand_pt)); if (lig->ligand_pts == NULL) return EFP_RESULT_NO_MEMORY; for (size_t i = 0; i < lig->n_ligand_pts; i++) { struct ligand_pt *pt = lig->ligand_pts + i; - size = efp->n_frag * sizeof(vec_t); pt->n_frag = efp->n_frag; - pt->fragment_field = (vec_t *) malloc(size); + pt->fragment_field = (vec_t *) calloc(efp->n_frag, sizeof(vec_t)); if (!pt->fragment_field) return EFP_RESULT_NO_MEMORY; } @@ -2745,10 +2723,14 @@ EFP_EXPORT enum efp_result efp_get_atom_mm_info(struct efp *efp, double *charges assert(coords); assert(charges); + if (efp->opts.special_fragment < 0) { + efp_log("special fragment is not defined"); + return EFP_RESULT_FATAL; + } int natom = 0; for (size_t fr_i = 0; fr_i < efp->n_frag; fr_i++) { // assume that this function is always used in QM/MM like models with QM fragment - if (efp->opts.special_fragment == fr_i) continue; + if ((size_t)efp->opts.special_fragment == fr_i) continue; struct frag *frag = efp->frags + fr_i; for (size_t j = 0; j < frag->n_atoms; j++) { struct efp_atom *atom = frag->atoms + j; @@ -2981,7 +2963,7 @@ efp_set_symmlist(struct efp *efp) // this needs to be changed for list settings of symmetry!!! efp->nsymm_frag = efp->n_lib; char name[32]; - char** unique_names=malloc(efp->n_lib * sizeof(name)); + char** unique_names=calloc(efp->n_lib, sizeof(name)); for (int i = 0; i < efp->n_lib; i++){ unique_names[i] = NULL; } diff --git a/src/efp.h b/src/efp.h index 391039bf..f2dd1a5b 100644 --- a/src/efp.h +++ b/src/efp.h @@ -963,7 +963,7 @@ efp_get_frag_rank(struct efp *efp, size_t frag_idx, int *rank); * * \return ::EFP_RESULT_SUCCESS on success or error code otherwise. */ -//enum efp_result efp_get_multipole_count(struct efp *efp, size_t *n_mult); +enum efp_result efp_get_multipole_count(struct efp *efp, size_t *n_mult); /** * Get the total number of high-order (large than monopoles) multipole points from EFP electrostatics. @@ -993,7 +993,7 @@ enum efp_result efp_get_mm_multipole_count(struct efp *efp, size_t *n_mult); * * \return ::EFP_RESULT_SUCCESS on success or error code otherwise. */ -//enum efp_result efp_get_multipole_coordinates(struct efp *efp, double *xyz); +enum efp_result efp_get_multipole_coordinates(struct efp *efp, double *xyz); /** * Get coordinates of high-order (higher than monopoles) electrostatics multipoles. @@ -1038,7 +1038,7 @@ enum efp_result efp_get_mm_multipole_coordinates(struct efp *efp, double *xyz); * * \return ::EFP_RESULT_SUCCESS on success or error code otherwise. */ -//enum efp_result efp_get_multipole_values(struct efp *efp, double *mult); +enum efp_result efp_get_multipole_values(struct efp *efp, double *mult); /** * Get high-order (higher than monopoles) electrostatics multipoles from EFP fragments. diff --git a/src/elec.c b/src/elec.c index 53924746..3382029b 100644 --- a/src/elec.c +++ b/src/elec.c @@ -327,15 +327,31 @@ mult_mult_grad(struct efp *efp, size_t fr_i_idx, size_t fr_j_idx, // gradient is not added to the special fragment in this case // this assumes that we use ml/efp fragment that induces field to other fragments due to its efp nature (multipoles and ind dipoles) // this might need to be changed if ml fragment uses ml-predicted charges instead - - if (!efp->opts.enable_elpot || (efp->opts.special_fragment != fr_i_idx && efp->opts.special_fragment != fr_j_idx)) { + if (!efp->opts.enable_elpot) { efp_add_force(efp->grad + fr_i_idx, CVEC(fr_i->x), CVEC(pt_i->x), &force, &torque_i); efp_sub_force(efp->grad + fr_j_idx, CVEC(fr_j->x), CVEC(pt_j->x), &force, &torque_j); - } - else if (efp->opts.enable_elpot && efp->opts.special_fragment == fr_i_idx) - efp_sub_force(efp->grad + fr_j_idx, CVEC(fr_j->x), CVEC(pt_j->x), &force, &torque_j); - else if (efp->opts.enable_elpot && efp->opts.special_fragment == fr_j_idx) - efp_add_force(efp->grad + fr_i_idx, CVEC(fr_i->x), CVEC(pt_i->x), &force, &torque_i); + } + else if ( efp->opts.special_fragment >= 0) { + if ((size_t)efp->opts.special_fragment != fr_i_idx && (size_t)efp->opts.special_fragment != fr_j_idx) { + efp_add_force(efp->grad + fr_i_idx, CVEC(fr_i->x), CVEC(pt_i->x), &force, &torque_i); + efp_sub_force(efp->grad + fr_j_idx, CVEC(fr_j->x), CVEC(pt_j->x), &force, &torque_j); + } + else if ((size_t)efp->opts.special_fragment == fr_i_idx) + efp_sub_force(efp->grad + fr_j_idx, CVEC(fr_j->x), CVEC(pt_j->x), &force, &torque_j); + else if ((size_t)efp->opts.special_fragment == fr_j_idx) + efp_add_force(efp->grad + fr_i_idx, CVEC(fr_i->x), CVEC(pt_i->x), &force, &torque_i); + } + else + printf("\nERROR: special fragment is not defined in elpot model!\n"); + + // if (!efp->opts.enable_elpot || (efp->opts.special_fragment != fr_i_idx && efp->opts.special_fragment != fr_j_idx)) { + // efp_add_force(efp->grad + fr_i_idx, CVEC(fr_i->x), CVEC(pt_i->x), &force, &torque_i); + // efp_sub_force(efp->grad + fr_j_idx, CVEC(fr_j->x), CVEC(pt_j->x), &force, &torque_j); + // } + // else if (efp->opts.enable_elpot && efp->opts.special_fragment == fr_i_idx) + // efp_sub_force(efp->grad + fr_j_idx, CVEC(fr_j->x), CVEC(pt_j->x), &force, &torque_j); + // else if (efp->opts.enable_elpot && efp->opts.special_fragment == fr_j_idx) + // efp_add_force(efp->grad + fr_i_idx, CVEC(fr_i->x), CVEC(pt_i->x), &force, &torque_i); } double @@ -373,13 +389,23 @@ efp_frag_frag_elec(struct efp *efp, size_t fr_i_idx, size_t fr_j_idx) // a check for torch special model with elpot on special fragment // gradient is not added to the special fragment in this case - if (!efp->opts.enable_elpot || (efp->opts.special_fragment != fr_i_idx && efp->opts.special_fragment != fr_j_idx)) { + if (!efp->opts.enable_elpot) { six_atomic_add_xyz(efp->grad + fr_i_idx, &force); six_atomic_sub_xyz(efp->grad + fr_j_idx, &force); } - else if (efp->opts.enable_elpot && efp->opts.special_fragment == fr_i_idx) six_atomic_sub_xyz(efp->grad + fr_j_idx, &force); - else if (efp->opts.enable_elpot && efp->opts.special_fragment == fr_j_idx) six_atomic_add_xyz(efp->grad + fr_i_idx, &force); - + else if ( efp->opts.special_fragment >= 0) { + if ((size_t)efp->opts.special_fragment != fr_i_idx && (size_t)efp->opts.special_fragment != fr_j_idx) { + six_atomic_add_xyz(efp->grad + fr_i_idx, &force); + six_atomic_sub_xyz(efp->grad + fr_j_idx, &force); + } + else if ((size_t)efp->opts.special_fragment == fr_i_idx) + six_atomic_sub_xyz(efp->grad + fr_j_idx, &force); + else if ((size_t)efp->opts.special_fragment == fr_j_idx) + six_atomic_add_xyz(efp->grad + fr_i_idx, &force); + } + else { + printf("\nERROR: special fragment is not defined in elpot model!\n"); + } return energy * swf.swf; } } @@ -387,7 +413,7 @@ efp_frag_frag_elec(struct efp *efp, size_t fr_i_idx, size_t fr_j_idx) static void rotate_quadrupole(const mat_t *rotmat, const double *in, double *out) { - double full_in[9], full_out[9]; + double full_in[9] = {0.0}, full_out[9] = {0.0}; for (size_t a = 0; a < 3; a++) for (size_t b = 0; b < 3; b++) @@ -403,7 +429,7 @@ rotate_quadrupole(const mat_t *rotmat, const double *in, double *out) static void rotate_octupole(const mat_t *rotmat, const double *in, double *out) { - double full_in[27], full_out[27]; + double full_in[27]={0.0}, full_out[27]={0.0}; for (size_t a = 0; a < 3; a++) for (size_t b = 0; b < 3; b++) @@ -601,8 +627,8 @@ compute_ai_elec_frag(struct efp *efp, size_t frag_idx) } */ for (size_t i = 0; i < fr_i->n_multipole_pts; i++) { + struct multipole_pt *pt_i = fr_i->multipole_pts + i; for (size_t j = 0; j < efp->n_ptc; j++) { - struct multipole_pt *pt_i = fr_i->multipole_pts + i; vec_t dr = vec_sub(CVEC(pt_i->x), efp->ptc_xyz + j); /* charge - monopole */ @@ -704,7 +730,6 @@ static void compute_ai_elec_range(struct efp *efp, size_t from, size_t to, void *data) { double energy = 0.0; - double energy_tmp = 0.0; (void)data; @@ -712,9 +737,11 @@ compute_ai_elec_range(struct efp *efp, size_t from, size_t to, void *data) #pragma omp parallel for schedule(dynamic) reduction(+:energy) #endif for (size_t i = from; i < to; i++) { + double energy_tmp = 0.0; // skip special fragment - if (i == efp->opts.special_fragment) - continue; + if (efp->opts.special_fragment >= -1) + if (i == (size_t)efp->opts.special_fragment) + continue; energy_tmp = compute_ai_elec_frag(efp, i); energy += energy_tmp; @@ -827,7 +854,7 @@ lj_energy(struct efp *efp, size_t fr_i_idx, size_t fr_j_idx, double energy = 0.0; double g = 0.0, g1 = 0.0, g2 = 0.0; - double r2, sr6 = 0.0; + double r2 = 0.0, sr6 = 0.0; if (combination_rule == 1) { // E_lj = (C_12/r^12 - C_6/r^6) @@ -960,12 +987,23 @@ efp_frag_frag_qq(struct efp *efp, size_t fr_i_idx, size_t fr_j_idx) // a check for torch special model with elpot on special fragment // gradient is not added to the special fragment in this case - if (!efp->opts.enable_elpot || (efp->opts.special_fragment != fr_i_idx && efp->opts.special_fragment != fr_j_idx)) { + if (!efp->opts.enable_elpot) { six_atomic_add_xyz(efp->grad + fr_i_idx, &force); six_atomic_sub_xyz(efp->grad + fr_j_idx, &force); } - else if (efp->opts.enable_elpot && efp->opts.special_fragment == fr_i_idx) six_atomic_sub_xyz(efp->grad + fr_j_idx, &force); - else if (efp->opts.enable_elpot && efp->opts.special_fragment == fr_j_idx) six_atomic_add_xyz(efp->grad + fr_i_idx, &force); + else if ( efp->opts.special_fragment >= 0) { + if ((size_t)efp->opts.special_fragment != fr_i_idx && (size_t)efp->opts.special_fragment != fr_j_idx) { + six_atomic_add_xyz(efp->grad + fr_i_idx, &force); + six_atomic_sub_xyz(efp->grad + fr_j_idx, &force); + } + else if ((size_t)efp->opts.special_fragment == fr_i_idx) + six_atomic_sub_xyz(efp->grad + fr_j_idx, &force); + else if ((size_t)efp->opts.special_fragment == fr_j_idx) + six_atomic_add_xyz(efp->grad + fr_i_idx, &force); + } + else { + printf("\nERROR: special fragment is not defined in elpot model!\n"); + } } return energy * swf.swf; @@ -1055,7 +1093,6 @@ static void compute_ai_qq_range(struct efp *efp, size_t from, size_t to, void *data) { double energy = 0.0; - double energy_tmp = 0.0; (void)data; @@ -1063,9 +1100,11 @@ compute_ai_qq_range(struct efp *efp, size_t from, size_t to, void *data) #pragma omp parallel for schedule(dynamic) reduction(+:energy) #endif for (size_t i = from; i < to; i++) { + double energy_tmp = 0.0; // skip special fragment - if (i == efp->opts.special_fragment) - continue; + if (efp->opts.special_fragment >= -1) + if (i == (size_t)efp->opts.special_fragment) + continue; // where are switching functions??? diff --git a/src/mathutil.h b/src/mathutil.h index 3c613d66..1a3a7468 100644 --- a/src/mathutil.h +++ b/src/mathutil.h @@ -369,7 +369,7 @@ euler_to_matrix(double a, double b, double c, mat_t *out) static inline void matrix_to_euler(const mat_t *rotmat, double *ea, double *eb, double *ec) { - double a, b, c, sinb; + double a=0.0, b=0.0, c=0.0, sinb=0.0; if (fabs(rotmat->zz - 1.0) < 1.0e-7) { b = 0.0; diff --git a/src/parse.c b/src/parse.c index 1c1b3a90..b53f03bf 100644 --- a/src/parse.c +++ b/src/parse.c @@ -34,7 +34,7 @@ #include "private.h" static void init_multipole_pt(struct multipole_pt *pt) { - memset(pt, 0, sizeof(*pt)); + memset(pt, 0, sizeof(struct multipole_pt)); pt->screen2 = 10.0; pt->screen0 = 10.0; pt->if_znuc = false; @@ -46,10 +46,6 @@ static void init_multipole_pt(struct multipole_pt *pt) { pt->if_scr0 = false; } -static void init_pol_pt(struct polarizable_pt *pt) { - memset(pt, 0, sizeof(*pt)); -} - static int tok(struct stream *stream, const char *id) { @@ -158,9 +154,7 @@ parse_coordinates(struct frag *frag, struct stream *stream) return EFP_RESULT_SUCCESS; } - struct efp_atom atom; - - memset(&atom, 0, sizeof(atom)); + struct efp_atom atom = {0.0}; if (!tok_label(stream, sizeof(atom.label), atom.label) || !tok_double(stream, &atom.x) || !tok_double(stream, &atom.y) || @@ -488,8 +482,8 @@ parse_polarizable_pts(struct frag *frag, struct stream *stream) struct polarizable_pt *pt = frag->polarizable_pts + frag->n_polarizable_pts - 1; // zero out all entries - init_pol_pt(pt); - + memset(pt, 0, sizeof(struct polarizable_pt)); + if (!efp_stream_advance(stream, 4)){ efp_log("parse_polarizable_pts() failure for fragment %s", frag->name); return EFP_RESULT_SYNTAX_ERROR; @@ -503,7 +497,7 @@ parse_polarizable_pts(struct frag *frag, struct stream *stream) } efp_stream_next_line(stream); - double m[9]; + double m[9]={0.0}; for (size_t i = 0; i < 9; i++) if (!tok_double(stream, m + i)){ @@ -531,7 +525,7 @@ parse_polarizable_pts(struct frag *frag, struct stream *stream) static enum efp_result parse_dynamic_polarizable_pts(struct frag *frag, struct stream *stream) { - double m[9]; + double m[9]={0.0}; efp_stream_next_line(stream); @@ -549,6 +543,8 @@ parse_dynamic_polarizable_pts(struct frag *frag, struct stream *stream) struct dynamic_polarizable_pt *pt = frag->dynamic_polarizable_pts + frag->n_dynamic_polarizable_pts - 1; + memset(pt, 0, sizeof(*pt)); + if (!efp_stream_advance(stream, 5)){ efp_log("parse_dynamic_polarizable_pts() failure for fragment %s", frag->name); @@ -699,6 +695,7 @@ parse_projection_basis(struct frag *frag, struct stream *stream) return EFP_RESULT_NO_MEMORY; struct shell *shell = atom->shells + atom->n_shells - 1; + memset(shell, 0, sizeof(*shell)); shell->type = efp_stream_get_char(stream); //printf("\n shell->type %c", shell->type); @@ -718,7 +715,7 @@ parse_projection_basis(struct frag *frag, struct stream *stream) efp_stream_next_line(stream); size_t cnt = (shell->type == 'L' ? 3 : 2) * shell->n_funcs; - shell->coef = (double *)malloc(cnt * sizeof(double)); + shell->coef = (double *)calloc(cnt, sizeof(double)); if (shell->coef == NULL) return EFP_RESULT_NO_MEMORY; double *ptr = shell->coef; @@ -778,8 +775,8 @@ parse_projection_wf(struct frag *frag, struct stream *stream) return EFP_RESULT_SYNTAX_ERROR; } - frag->xr_wf = (double *)malloc( - frag->n_lmo * frag->xr_wf_size * sizeof(double)); + frag->xr_wf = (double *)calloc( + frag->n_lmo * frag->xr_wf_size, sizeof(double)); if (frag->xr_wf == NULL) return EFP_RESULT_NO_MEMORY; @@ -832,7 +829,7 @@ parse_fock_mat(struct frag *frag, struct stream *stream) efp_stream_next_line(stream); size_t size = frag->n_lmo * (frag->n_lmo + 1) / 2; - frag->xr_fock_mat = (double *)malloc(size * sizeof(double)); + frag->xr_fock_mat = (double *)calloc(size, sizeof(double)); if (frag->xr_fock_mat == NULL) return EFP_RESULT_NO_MEMORY; @@ -865,7 +862,7 @@ parse_lmo_centroids(struct frag *frag, struct stream *stream) efp_log("number of LMO centroids is zero"); return EFP_RESULT_SYNTAX_ERROR; } - frag->lmo_centroids = (vec_t *)malloc(frag->n_lmo * sizeof(vec_t)); + frag->lmo_centroids = (vec_t *)calloc(frag->n_lmo, sizeof(vec_t)); if (frag->lmo_centroids == NULL) return EFP_RESULT_NO_MEMORY; @@ -979,7 +976,7 @@ skip_ctfok(struct frag *frag, struct stream *stream) static enum efp_result parse_dipquad_polarizable_pts(struct frag *frag, struct stream *stream) { - double m[27]; + double m[27]={0.0}; efp_stream_next_line(stream); @@ -997,6 +994,7 @@ parse_dipquad_polarizable_pts(struct frag *frag, struct stream *stream) struct dipquad_polarizable_pt *pt = frag->dipquad_polarizable_pts + frag->n_dynamic_polarizable_pts - 1; + memset(pt, 0, sizeof(*pt)); if (!efp_stream_advance(stream, 5)){ printf("problem with fragment %s", frag->name); @@ -1142,7 +1140,7 @@ parse_screen(struct frag *frag, struct stream *stream) double *scr; char type; - scr = (double *)malloc(frag->n_multipole_pts * sizeof(double)); + scr = (double *)calloc(frag->n_multipole_pts, sizeof(double)); if (scr == NULL) return EFP_RESULT_NO_MEMORY; type = efp_stream_get_char(stream); @@ -1227,7 +1225,7 @@ parse_screen(struct frag *frag, struct stream *stream) struct multipole_pt tmp_pt; memset(&tmp_pt, 0, sizeof(tmp_pt)); - double tmp_screen; + double tmp_screen = 0.0; if (!tok_label(stream, sizeof(tmp_pt.label), tmp_pt.label) || !tok_double(stream, NULL) || !tok_double(stream, &tmp_screen)){ @@ -1278,7 +1276,7 @@ parse_xrfit(struct frag *frag, struct stream *stream) return EFP_RESULT_SYNTAX_ERROR; } - frag->xrfit = (double *)malloc(frag->n_lmo * 4 * sizeof(double)); + frag->xrfit = (double *)calloc(frag->n_lmo * 4, sizeof(double)); if (frag->xrfit == NULL) return EFP_RESULT_NO_MEMORY; efp_stream_next_line(stream); @@ -1445,8 +1443,7 @@ parse_mm_atomtype(struct frag *frag, struct stream *stream) counter, frag->n_atoms, frag->name); return EFP_RESULT_SUCCESS; } - struct efp_atom atom; - memset(&atom, 0, sizeof(atom)); + struct efp_atom atom = {0.0}; if (!tok_label(stream, sizeof(atom.label), atom.label) || !tok_label(stream, sizeof(atom.ff_label), atom.ff_label)){ printf("problem with fragment %s", frag->name); diff --git a/src/pol.c b/src/pol.c index f5b938ce..3309d2bb 100644 --- a/src/pol.c +++ b/src/pol.c @@ -68,14 +68,14 @@ get_multipole_field(const vec_t *xyz, const struct multipole_pt *mult_pt, xyz->z - mult_pt->z - swf->cell.z }; - double t1, t2; + double t1=0.0, t2=0.0; double r = vec_len(&dr); double r3 = r * r * r; double r5 = r3 * r * r; double r7 = r5 * r * r; /* charge */ - if (mult_pt->if_mon || mult_pt->if_dip) { + if (mult_pt->if_mon || mult_pt->if_znuc) { field.x += swf->swf * (mult_pt->monopole + mult_pt->znuc) * dr.x / r3; field.y += swf->swf * (mult_pt->monopole + mult_pt->znuc) * dr.y / r3; field.z += swf->swf * (mult_pt->monopole + mult_pt->znuc) * dr.z / r3; @@ -178,6 +178,8 @@ get_elec_field(const struct efp *efp, size_t frag_idx, size_t pt_idx) const struct polarizable_pt *pt = fr_j->polarizable_pts + pt_idx; vec_t elec_field = vec_zero; + // ATTENTION! field due to MM atoms is not added here. Is it needed in any circumstances? + for (size_t i = 0; i < efp->n_frag; i++) { if (i == frag_idx ) continue; @@ -186,8 +188,9 @@ get_elec_field(const struct efp *efp, size_t frag_idx, size_t pt_idx) continue; // this might need to be changed to a more careful separation of // elec and pol contributions to the field - if (i == efp->opts.special_fragment && !(efp->opts.special_terms & EFP_SPEC_TERM_POL)) - continue; + if (efp->opts.special_fragment >=0) + if (i == (size_t)efp->opts.special_fragment && !(efp->opts.special_terms & EFP_SPEC_TERM_POL)) + continue; const struct frag *fr_i = efp->frags + i; struct swf swf = efp_make_swf(efp, fr_i, fr_j, 0); if (swf.swf == 0.0) @@ -247,8 +250,8 @@ get_ligand_field(const struct efp *efp, size_t frag_idx, size_t pt_idx, int liga /* ligand is not QM */ if (ligand_idx != -1) { - const struct frag *fr_i = efp->frags + ligand_idx; - if (efp_skip_frag_pair(efp, ligand_idx, frag_idx)) + const struct frag *fr_i = efp->frags + (size_t)ligand_idx; + if (efp_skip_frag_pair(efp, (size_t)ligand_idx, frag_idx)) return elec_field; struct swf swf = efp_make_swf(efp, fr_i, fr_j, 0); @@ -306,8 +309,8 @@ add_electron_density_field(struct efp *efp) if (efp->get_electron_density_field == NULL) return EFP_RESULT_SUCCESS; - xyz = (vec_t *)malloc(efp->n_polarizable_pts * sizeof(vec_t)); - field = (vec_t *)malloc(efp->n_polarizable_pts * sizeof(vec_t)); + xyz = (vec_t *)calloc(efp->n_polarizable_pts, sizeof(vec_t)); + field = (vec_t *)calloc(efp->n_polarizable_pts, sizeof(vec_t)); for (size_t i = 0, idx = 0; i < efp->n_frag; i++) { struct frag *frag = efp->frags + i; @@ -349,9 +352,10 @@ compute_elec_field_range(struct efp *efp, size_t from, size_t to, void *data) #pragma omp parallel for schedule(dynamic) #endif for (size_t i = from; i < to; i++) { - if (i == efp->opts.special_fragment && - !(efp->opts.special_terms & EFP_SPEC_TERM_POL)) - continue; + if (efp->opts.special_fragment >=0) + if (i == (size_t)efp->opts.special_fragment && + !(efp->opts.special_terms & EFP_SPEC_TERM_POL)) + continue; // const struct frag *frag = efp->frags + i; struct frag *frag = efp->frags + i; @@ -386,13 +390,12 @@ compute_fragment_field_range(struct efp *efp, size_t from, size_t to, void *data { (void)data; - struct ligand *ligand = efp->ligand; - #ifdef _OPENMP #pragma omp parallel for schedule(dynamic) #endif for (size_t i = from; i < to; i++) { - struct frag *frag = efp->frags + i; + // struct frag *frag = efp->frags + i; + struct ligand *ligand = efp->ligand; for (size_t j = 0; j < ligand->n_ligand_pts; j++) { struct ligand_pt *pt = ligand->ligand_pts + j; @@ -488,9 +491,10 @@ get_induced_dipole_field(struct efp *efp, size_t frag_idx, if (efp->opts.symmetry == 0 && efp_skip_frag_pair(efp, frag_idx, j)) continue; - if (j == efp->opts.special_fragment && - !(efp->opts.special_terms & EFP_SPEC_TERM_POL)) - continue; + if (efp->opts.special_fragment >=0) + if (j == (size_t)efp->opts.special_fragment && + !(efp->opts.special_terms & EFP_SPEC_TERM_POL)) + continue; struct frag *fr_j = efp->frags + j; struct swf swf = efp_make_swf(efp, fr_i, fr_j, 0); @@ -616,9 +620,10 @@ compute_id_range(struct efp *efp, size_t from, size_t to, void *data) #endif for (size_t i = from; i < to; i++) { - if (i == efp->opts.special_fragment && - !(efp->opts.special_terms & EFP_SPEC_TERM_POL)) - continue; + if (efp->opts.special_fragment >=0) + if (i == (size_t)efp->opts.special_fragment && + !(efp->opts.special_terms & EFP_SPEC_TERM_POL)) + continue; struct frag *frag = efp->frags + i; @@ -767,25 +772,27 @@ compute_energy_range(struct efp *efp, size_t from, size_t to, void *data) // if_pairwise tells whether pairwise calculations with EFP (not QM) ligand bool if_pairwise = efp->opts.enable_pairwise && efp->opts.ligand > -1; - struct ligand *ligand; - if (if_pairwise) - ligand = efp->ligand; - #ifdef _OPENMP #pragma omp parallel for schedule(dynamic) reduction(+:energy) #endif for (size_t i = from; i < to; i++) { + // zeroing out polarization pair energies is a must + efp->pair_energies[i].polarization = 0.0; + + struct ligand *ligand; + if (if_pairwise) + ligand = efp->ligand; + // skip energy contribution for a special fragment in case of torch model with elpot // this assumes that we use ml/efp fragment that induces field to other fragments due to its efp nature (multipoles and ind dipoles) // this needs to be changed if ml fragment uses ml-predicted charges instead - if (efp->opts.enable_elpot && efp->opts.special_fragment == i) continue; + if (efp->opts.enable_elpot && efp->opts.special_fragment >=0) + if ((size_t)efp->opts.special_fragment == i) + continue; struct frag *frag = efp->frags + i; - // zeroing out polarization pair energies is a must - efp->pair_energies[i].polarization = 0.0; - for (size_t j = 0; j < frag->n_polarizable_pts; j++) { struct polarizable_pt *pt = frag->polarizable_pts + j; @@ -1079,7 +1086,7 @@ compute_grad_point(struct efp *efp, size_t frag_idx, size_t pt_idx) const struct frag *fr_i = efp->frags + frag_idx; const struct polarizable_pt *pt_i = fr_i->polarizable_pts + pt_idx; vec_t force, add_i, add_j, force_, add_i_, add_j_; - double e; + double ene = 0.0, energy = 0.0; vec_t dipole_i = { 0.5 * (pt_i->indip.x + pt_i->indipconj.x), @@ -1088,6 +1095,10 @@ compute_grad_point(struct efp *efp, size_t frag_idx, size_t pt_idx) }; for (size_t j = 0; j < efp->n_frag; j++) { + + /* energy without switching applied */ + energy = 0.0; + if (j == frag_idx || efp_skip_frag_pair(efp, frag_idx, j)) continue; @@ -1096,24 +1107,32 @@ compute_grad_point(struct efp *efp, size_t frag_idx, size_t pt_idx) // this needs to be changed if ml fragment uses ml-predicted charges instead // this is true for normal cases not related to torch model with elpot - bool not_torch_elpot = !efp->opts.enable_elpot || (efp->opts.special_fragment != frag_idx && efp->opts.special_fragment != j); + bool not_torch_elpot = true; // true when torch with elpot is invoked and fr_idx is the ml fragment - bool torch_elpot_i = efp->opts.enable_elpot && (efp->opts.special_fragment == frag_idx); + bool torch_elpot_i = false; // true when torch with elpot is invoked and j is the ml fragment - bool torch_elpot_j = efp->opts.enable_elpot && (efp->opts.special_fragment == j); + bool torch_elpot_j = false; + + if (efp->opts.special_fragment >= 0) { + not_torch_elpot = !efp->opts.enable_elpot || + ((size_t)efp->opts.special_fragment != frag_idx && + (size_t)efp->opts.special_fragment != j); + torch_elpot_i = efp->opts.enable_elpot && ((size_t)efp->opts.special_fragment == frag_idx); + torch_elpot_j = efp->opts.enable_elpot && ((size_t)efp->opts.special_fragment == j); + } struct frag *fr_j = efp->frags + j; struct swf swf = efp_make_swf(efp, fr_i, fr_j, 0); if (swf.swf == 0.0) continue; - /* energy without switching applied */ - double energy = 0.0; /* induced dipole - multipoles */ for (size_t k = 0; k < fr_j->n_multipole_pts; k++) { struct multipole_pt *pt_j = fr_j->multipole_pts + k; + ene = 0.0; + vec_t dr = { pt_j->x - pt_i->x - swf.cell.x, pt_j->y - pt_i->y - swf.cell.y, @@ -1138,7 +1157,7 @@ compute_grad_point(struct efp *efp, size_t frag_idx, size_t pt_idx) /* induced dipole - charge+monopole */ if (pt_j->if_mon || pt_j->if_znuc) { double qj = pt_j->monopole + pt_j->znuc; - e = -efp_charge_dipole_energy(qj, &dipole_i, &dr); + ene = -efp_charge_dipole_energy(qj, &dipole_i, &dr); efp_charge_dipole_grad(qj, &dipole_i, &dr, &force_, &add_j_, &add_i_); vec_negate(&force_); @@ -1148,7 +1167,7 @@ compute_grad_point(struct efp *efp, size_t frag_idx, size_t pt_idx) /* induced dipole - dipole */ if (pt_j->if_dip) { - e += efp_dipole_dipole_energy(&dipole_i, + ene += efp_dipole_dipole_energy(&dipole_i, &pt_j->dipole, &dr); efp_dipole_dipole_grad(&dipole_i, &pt_j->dipole, &dr, &force_, &add_i_, &add_j_); @@ -1159,7 +1178,7 @@ compute_grad_point(struct efp *efp, size_t frag_idx, size_t pt_idx) /* induced dipole - quadrupole */ if (pt_j->if_quad) { - e += efp_dipole_quadrupole_energy(&dipole_i, + ene += efp_dipole_quadrupole_energy(&dipole_i, pt_j->quadrupole, &dr); efp_dipole_quadrupole_grad(&dipole_i, pt_j->quadrupole, &dr, &force_, &add_i_, &add_j_); @@ -1173,9 +1192,9 @@ compute_grad_point(struct efp *efp, size_t frag_idx, size_t pt_idx) vec_scale(&add_i, p1); vec_scale(&add_j, p1); - force.x += p2 * e * dr.x; - force.y += p2 * e * dr.y; - force.z += p2 * e * dr.z; + force.x += p2 * ene * dr.x; + force.y += p2 * ene * dr.y; + force.z += p2 * ene * dr.z; vec_scale(&force, swf.swf); vec_scale(&add_i, swf.swf); @@ -1187,7 +1206,7 @@ compute_grad_point(struct efp *efp, size_t frag_idx, size_t pt_idx) if (not_torch_elpot) { efp_add_force(efp->grad + frag_idx, CVEC(fr_i->x), CVEC(pt_i->x), &force, &add_i); efp_sub_force(efp->grad + j, CVEC(fr_j->x), CVEC(pt_j->x), &force, &add_j); - energy += p1 * e; + energy += p1 * ene; } // adding gradients to non-ML fragment only in torch elpot model @@ -1200,6 +1219,8 @@ compute_grad_point(struct efp *efp, size_t frag_idx, size_t pt_idx) struct polarizable_pt *pt_j = fr_j->polarizable_pts+jj; // size_t idx_j = fr_j->polarizable_offset+jj; + ene = 0.0; + vec_t dr = { pt_j->x - pt_i->x - swf.cell.x, pt_j->y - pt_i->y - swf.cell.y, @@ -1223,7 +1244,7 @@ compute_grad_point(struct efp *efp, size_t frag_idx, size_t pt_idx) fr_j->pol_damp); } - e = efp_dipole_dipole_energy(&half_dipole_i, &pt_j->indipconj, &dr); + ene = efp_dipole_dipole_energy(&half_dipole_i, &pt_j->indipconj, &dr); efp_dipole_dipole_grad(&half_dipole_i, &pt_j->indipconj, &dr, &force, &add_i, &add_j); vec_negate(&add_j); @@ -1231,9 +1252,9 @@ compute_grad_point(struct efp *efp, size_t frag_idx, size_t pt_idx) vec_scale(&add_i, p1); vec_scale(&add_j, p1); - force.x += p2 * e * dr.x; - force.y += p2 * e * dr.y; - force.z += p2 * e * dr.z; + force.x += p2 * ene * dr.x; + force.y += p2 * ene * dr.y; + force.z += p2 * ene * dr.z; vec_scale(&force, swf.swf); vec_scale(&add_i, swf.swf); @@ -1245,7 +1266,7 @@ compute_grad_point(struct efp *efp, size_t frag_idx, size_t pt_idx) if (not_torch_elpot) { efp_add_force(efp->grad + frag_idx, CVEC(fr_i->x), CVEC(pt_i->x), &force, &add_i); efp_sub_force(efp->grad + j, CVEC(fr_j->x), CVEC(pt_j->x), &force, &add_j); - energy += p1 * e; + energy += p1 * ene; } // adding gradients to non-ML fragment only in torch elpot model @@ -1340,6 +1361,9 @@ efp_get_electric_field(struct efp *efp, size_t frag_idx, const double *xyz, assert(xyz); assert(field); + // ATTENTION! This function does not include the electric field due to MM atoms + // might need to reconsider it! + const struct frag *frag = efp->frags + frag_idx; vec_t elec_field = vec_zero; diff --git a/src/util.c b/src/util.c index c38e0ab8..d7eec956 100644 --- a/src/util.c +++ b/src/util.c @@ -182,12 +182,48 @@ efp_find_lib(struct efp *efp, const char *name) return NULL; } +// void +// efp_add_stress(const vec_t *dr, const vec_t *force, mat_t *stress) +// { +// #ifdef _OPENMP +// #pragma omp critical +// #endif +// { +// stress->xx += dr->x * force->x; +// stress->xy += dr->x * force->y; +// stress->xz += dr->x * force->z; +// stress->yx += dr->y * force->x; +// stress->yy += dr->y * force->y; +// stress->yz += dr->y * force->z; +// stress->zx += dr->z * force->x; +// stress->zy += dr->z * force->y; +// stress->zz += dr->z * force->z; +// } +// } + void efp_add_stress(const vec_t *dr, const vec_t *force, mat_t *stress) { #ifdef _OPENMP -#pragma omp critical -#endif +#pragma omp atomic + stress->xx += dr->x * force->x; +#pragma omp atomic + stress->xy += dr->x * force->y; +#pragma omp atomic + stress->xz += dr->x * force->z; +#pragma omp atomic + stress->yx += dr->y * force->x; +#pragma omp atomic + stress->yy += dr->y * force->y; +#pragma omp atomic + stress->yz += dr->y * force->z; +#pragma omp atomic + stress->zx += dr->z * force->x; +#pragma omp atomic + stress->zy += dr->z * force->y; +#pragma omp atomic + stress->zz += dr->z * force->z; +#else { stress->xx += dr->x * force->x; stress->xy += dr->x * force->y; @@ -199,6 +235,7 @@ efp_add_stress(const vec_t *dr, const vec_t *force, mat_t *stress) stress->zy += dr->z * force->y; stress->zz += dr->z * force->z; } +#endif } void diff --git a/src/xr.c b/src/xr.c index 552622b9..b5c1ebe3 100644 --- a/src/xr.c +++ b/src/xr.c @@ -209,7 +209,7 @@ lmo_lmo_xr_grad(struct efp *efp, size_t fr_i_idx, size_t fr_j_idx, six_t ds_ij = lmo_ds[ij]; six_t dt_ij = lmo_dt[ij]; - double t1, t2; + double t1=0.0, t2=0.0; vec_t force = vec_zero, torque_i = vec_zero; /* first part */ @@ -535,12 +535,12 @@ efp_frag_frag_xr(struct efp *efp, size_t frag_i, size_t frag_j, double *lmo_s, size_t ij_wf_size = fr_i->xr_wf_size * fr_j->xr_wf_size; size_t ij_nlmo = fr_i->n_lmo * fr_j->n_lmo; size_t ij_nlmo_wf_size = fr_i->n_lmo * fr_j->xr_wf_size; - double *s = (double *) malloc(ij_wf_size * sizeof(double)); - double *t = (double *) malloc(ij_wf_size * sizeof(double)); - double *lmo_t = (double *) malloc(ij_nlmo * sizeof(double)); - double *tmp = (double *) malloc(ij_nlmo_wf_size * sizeof(double)); - struct xr_atom *atoms_j = (struct xr_atom *) malloc( - fr_j->n_xr_atoms * sizeof(struct xr_atom)); + double *s = (double *) calloc(ij_wf_size, sizeof(double)); + double *t = (double *) calloc(ij_wf_size, sizeof(double)); + double *lmo_t = (double *) calloc(ij_nlmo, sizeof(double)); + double *tmp = (double *) calloc(ij_nlmo_wf_size, sizeof(double)); + struct xr_atom *atoms_j = (struct xr_atom *) calloc( + fr_j->n_xr_atoms, sizeof(struct xr_atom)); for (size_t j = 0; j < fr_j->n_xr_atoms; j++) { atoms_j[j] = fr_j->xr_atoms[j]; @@ -600,11 +600,11 @@ efp_frag_frag_xr(struct efp *efp, size_t frag_i, size_t frag_j, double *lmo_s, /* compute gradient */ - six_t *ds = (six_t *) malloc(ij_wf_size * sizeof(six_t)); - six_t *dt = (six_t *) malloc(ij_wf_size * sizeof(six_t)); - six_t *lmo_dt = (six_t *) malloc(ij_nlmo * sizeof(six_t)); - six_t *sixtmp = (six_t *) malloc(ij_nlmo_wf_size * sizeof(six_t)); - double *lmo_tmp = (double *) malloc(ij_nlmo * sizeof(double)); + six_t *ds = (six_t *) calloc(ij_wf_size, sizeof(six_t)); + six_t *dt = (six_t *) calloc(ij_wf_size, sizeof(six_t)); + six_t *lmo_dt = (six_t *) calloc(ij_nlmo, sizeof(six_t)); + six_t *sixtmp = (six_t *) calloc(ij_nlmo_wf_size, sizeof(six_t)); + double *lmo_tmp = (double *) calloc(ij_nlmo, sizeof(double)); efp_st_int_deriv(fr_i->n_xr_atoms, fr_i->xr_atoms, fr_j->n_xr_atoms, atoms_j, @@ -731,7 +731,7 @@ rotate_func_f(const mat_t *rotmat, const double *in, double *out) const double norm1 = sqrt(5.0) / 3.0; const double norm2 = sqrt(3.0) / 2.0; - double full_in[27], full_out[27]; + double full_in[27]={0.0}, full_out[27]={0.0}; for (size_t a = 0; a < 3; a++) for (size_t b = 0; b < 3; b++) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 3ea0dd5c..8eded856 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -65,8 +65,11 @@ foreach(test_name elec_3a elec_3b elpot + elpot_frag grad_1 hess_1 + lj_1 + lj_2 md_1 md_2 md_3 @@ -83,6 +86,8 @@ foreach(test_name pol_2b pol_3a pol_3b + polab_0 + polab_1 print qm_1a qm_1b @@ -90,6 +95,10 @@ foreach(test_name qm_2b reduced sp_1 + sp_2 + spec_frag_1 + spec_frag_2 + spec_frag_base symm_1 symm_2full symm_2 diff --git a/tests/atom_coord_2.in b/tests/atom_coord_2.in index 9d45be73..f564abdd 100644 --- a/tests/atom_coord_2.in +++ b/tests/atom_coord_2.in @@ -2,7 +2,7 @@ run_type grad coord atoms elec_damp screen fraglib_path ../fraglib -print 3 +print 1 fragment h2o_l A01O1 -3.394000 -1.900000 -3.700000 diff --git a/tests/pytests/conftest.py b/tests/pytests/conftest.py index c5875acb..9ef25c1e 100644 --- a/tests/pytests/conftest.py +++ b/tests/pytests/conftest.py @@ -1,4 +1,11 @@ import pytest +from libefp2py import read_libefp_input + + +@pytest.fixture +def pyjob_prepper(): + """Converts efpmd input into py format.""" + return read_libefp_input @pytest.fixture(scope="session", autouse=True) diff --git a/tests/pytests/output.dat b/tests/pytests/output.dat new file mode 100644 index 00000000..6d40fd62 --- /dev/null +++ b/tests/pytests/output.dat @@ -0,0 +1,8 @@ + => Loading Basis Set <= + + Name: 6-31G* + Role: ORBITAL + Keyword: BASIS + atoms 1 entry O line 145 file /Users/lyuda/anaconda3/envs/psi4_2026/share/psi4/basis/6-31gs.gbs + atoms 2-3 entry H line 44 file /Users/lyuda/anaconda3/envs/psi4_2026/share/psi4/basis/6-31gs.gbs + diff --git a/tests/pytests/test_coverage.py b/tests/pytests/test_coverage.py index 97c80b63..82c15541 100644 --- a/tests/pytests/test_coverage.py +++ b/tests/pytests/test_coverage.py @@ -6,12 +6,12 @@ from qcelemental.testing import compare -def test_grad_fail(): - asdf = system_1() - asdf.compute(do_gradient=False) - - with pytest.raises(pylibefp.Fatal) as e_info: - grad = asdf.get_gradient() +#def test_grad_fail(): +# asdf = system_1() +# asdf.compute(do_gradient=False) +# +# with pytest.raises(pylibefp.Fatal) as e_info: +# grad = asdf.get_gradient() #def test_frag_file_fail(): diff --git a/tests/pytests/test_efpefp.py b/tests/pytests/test_efpefp.py index ab87085c..fa276a89 100644 --- a/tests/pytests/test_efpefp.py +++ b/tests/pytests/test_efpefp.py @@ -118,7 +118,7 @@ def test_total_1a(): 'pol': True, # 'pol_damp': 'tt', 'disp': True, 'disp_damp': 'tt', - 'print': 2 + 'print': 1 }) asdf.compute() ene = asdf.get_energy() diff --git a/tests/pytests/test_efpefp2.py b/tests/pytests/test_efpefp2.py new file mode 100644 index 00000000..6b125f2e --- /dev/null +++ b/tests/pytests/test_efpefp2.py @@ -0,0 +1,72 @@ +import pylibefp +from qcelemental.testing import compare, compare_values +import pprint +import pytest +import os + +FILES = [ + 'atom_coord.in', 'atom_coord_2.in', 'grad_1.in', 'lj_1.in', 'lj_2.in', + 'pairwise_0.in', 'pairwise_1.in', 'pairwise_2.in', 'pairwise_x.in', 'pbc_1.in', 'pbc_2.in', + 'reduced.in', 'spec_frag_1.in', 'spec_frag_2.in', 'spec_frag_base.in', + 'symm_1.in', 'symm_2.in', 'symm_2full.in', 'symm_2pw.in' +] + +b2a = 0.529177 +a2b = 1.0 / b2a + +def frag_setup(test_name, pyjob_prepper): + # coordinates in Bohr + coord_type, frags, frag_coords, efp_options, if_gradient, ref_energy, periodic_box = pyjob_prepper(test_name) + #print(frag_coords) + + efp = pylibefp.core.efp() + efp.add_potential(frags) + efp.add_fragment(frags) + for i in range(len(frags)): + efp.set_frag_coordinates(i, coord_type, frag_coords[i]) + efp.prepare() + + efp.set_opts(efp_options) + if periodic_box: + #print('box1', periodic_box) + efp.set_periodic_box(periodic_box) + #print('box2', efp.get_periodic_box()) + + #print(frag_coords) + #pprint.pprint(efp_options) + efp.compute(do_gradient = if_gradient) + ene = efp.get_energy() + + # print pairwise components + #if 'enable_pairwise' in efp_options.keys(): + # if efp_options['enable_pairwise'] in [True, 'true', 1]: + # efp.print_pairwise_energies() + + print(efp.energy_summary()) + if if_gradient: + print(efp.gradient_summary()) + if ref_energy != 0.0: + assert compare_values(ref_energy, ene['total'], atol=1.e-5, return_message=True), 'FAILED' + +# ##### +# if __name__ == '__main__': +# files = ['atom_coord.in', 'atom_coord_2.in', 'grad_1.in', 'lj_1.in', 'lj_2.in', +# 'pairwise_0.in', 'pairwise_1.in', 'pairwise_2.in', 'pairwise_x.in', 'pbc_1.in', 'pbc_2.in', +# 'reduced.in', 'spec_frag_0.in', 'spec_frag_1.in', 'spec_frag_2.in', 'spec_frag_base.in', 'spec_frag_ref.in', +# 'symm_1.in', 'symm_2.in', 'symm_2full.in', 'symm_2pw.in'] + +# # running for all tests in files list +# for f in files: +# print(f'\nComputing {f}...') +# frag_setup('../'+f) + +# # single test execution +# frag_setup('../symm_2pw.in') + +BASE_DIR = os.path.dirname(os.path.abspath(__file__)) + +@pytest.mark.parametrize("filename", FILES) +def test_frag_setup(filename, pyjob_prepper): + print(f'\nComputing {filename}...') + full_path = os.path.join(BASE_DIR, '..', filename) + frag_setup(full_path, pyjob_prepper) \ No newline at end of file diff --git a/tests/pytests/test_efpefp_new.py b/tests/pytests/test_efpefp_new.py deleted file mode 100644 index 44284c98..00000000 --- a/tests/pytests/test_efpefp_new.py +++ /dev/null @@ -1,57 +0,0 @@ -import libefp2py -import pylibefp -from qcelemental.testing import compare, compare_values -import pprint - - -b2a = 0.529177 -a2b = 1.0 / b2a - -def frag_setup(test_name): - # coordinates in Bohr - coord_type, frags, frag_coords, efp_options, if_gradient, ref_energy, periodic_box = libefp2py.read_libefp_input(test_name) - #print(frag_coords) - - efp = pylibefp.core.efp() - efp.add_potential(frags) - efp.add_fragment(frags) - for i in range(len(frags)): - efp.set_frag_coordinates(i, coord_type, frag_coords[i]) - efp.prepare() - - efp.set_opts(efp_options) - if periodic_box: - #print('box1', periodic_box) - efp.set_periodic_box(periodic_box) - #print('box2', efp.get_periodic_box()) - - #print(frag_coords) - #pprint.pprint(efp_options) - efp.compute(do_gradient = if_gradient) - ene = efp.get_energy() - - # print pairwise components - #if 'enable_pairwise' in efp_options.keys(): - # if efp_options['enable_pairwise'] in [True, 'true', 1]: - # efp.print_pairwise_energies() - - print(efp.energy_summary()) - if if_gradient: - print(efp.gradient_summary()) - if ref_energy != 0.0: - assert compare_values(ref_energy, ene['total'], atol=1.e-5, return_message=True), 'FAILED' - -##### -if __name__ == '__main__': - files = ['atom_coord.in', 'atom_coord_2.in', 'grad_1.in', 'lj_1.in', 'lj_2.in', - 'pairwise_0.in', 'pairwise_1.in', 'pairwise_2.in', 'pairwise_x.in', 'pbc_1.in', 'pbc_2.in', - 'reduced.in', 'spec_frag_0.in', 'spec_frag_1.in', 'spec_frag_2.in', 'spec_frag_base.in', 'spec_frag_ref.in', - 'symm_1.in', 'symm_2.in', 'symm_2full.in', 'symm_2pw.in'] - - # running for all tests in files list - for f in files: - print(f'\nComputing {f}...') - frag_setup('../'+f) - - # single test execution - frag_setup('../symm_2pw.in') diff --git a/tests/pytests/test_scf.py b/tests/pytests/test_scf.py index 4d1ca339..9cf9c4d6 100644 --- a/tests/pytests/test_scf.py +++ b/tests/pytests/test_scf.py @@ -102,13 +102,14 @@ def modify_Fock_permanent(mol, nbf, efpobj): for pole in range(20): efp_ints[pole] = np.asarray(p4_efp_ints[pole]) - # add frag atom Z into multipole charge (when pos'n of atom matches mp) - for ifr in range(n_fr): - atoms = efpobj.get_frag_atoms(ifr) - for iat in range(natoms[ifr]): - xyz_atom = [atoms[iat]['x'], atoms[iat]['y'], atoms[iat]['z']] - if np.allclose(xyz_atom, origin, atol=1e-10): - val_mp[imp, 0] += atoms[iat]['Z'] + # LVS: atom Z charges are already included in multipoles as of 2025 + # # add frag atom Z into multipole charge (when pos'n of atom matches mp) + # for ifr in range(n_fr): + # atoms = efpobj.get_frag_atoms(ifr) + # for iat in range(natoms[ifr]): + # xyz_atom = [atoms[iat]['x'], atoms[iat]['y'], atoms[iat]['z']] + # if np.allclose(xyz_atom, origin, atol=1e-10): + # val_mp[imp, 0] += atoms[iat]['Z'] # scale multipole integrals by multipole magnitudes. result goes into V for pole in range(20): @@ -380,6 +381,10 @@ def field_fn(xyz): # <-- efp: add in permanent moment contribution and cache Vefp = modify_Fock_permanent(mol, nbf, efpmol) + + print("Vefp") + print(Vefp) + assert compare(1, np.allclose(Vefp, ref_V2), 'EFP permanent Fock contrib') H = H + Vefp Horig = H.copy() diff --git a/tests/pytests/timer.dat b/tests/pytests/timer.dat new file mode 100644 index 00000000..b6fc97b2 --- /dev/null +++ b/tests/pytests/timer.dat @@ -0,0 +1,128 @@ + +Host: MacBook-Pro-5.local + +Timers On : Sun Apr 19 22:37:26 2026 +Timers Off: Sun Apr 19 22:37:35 2026 + +Wall Time: 9.11 seconds + + Time (seconds) +Module User System Wall Calls +Libint2ERI::Libint2ERI : 0.000u 0.000s 0.001w 1 calls + +-------------------------------------------------------------------------------------- +Libint2ERI::Libint2ERI : 0.000u 0.000s 0.001w 1 calls + +************************************************************************************** + +Host: MacBook-Pro-5.local + +Timers On : Sun Apr 19 22:45:19 2026 +Timers Off: Sun Apr 19 22:45:20 2026 + +Wall Time: 0.46 seconds + + Time (seconds) +Module User System Wall Calls +Libint2ERI::Libint2ERI : 0.000u 0.000s 0.007w 1 calls + +-------------------------------------------------------------------------------------- +Libint2ERI::Libint2ERI : 0.000u 0.000s 0.007w 1 calls + +************************************************************************************** + +Host: MacBook-Pro-5.local + +Timers On : Sun Apr 19 22:47:15 2026 +Timers Off: Sun Apr 19 22:47:16 2026 + +Wall Time: 0.44 seconds + + Time (seconds) +Module User System Wall Calls +Libint2ERI::Libint2ERI : 0.000u 0.000s 0.005w 1 calls + +-------------------------------------------------------------------------------------- +Libint2ERI::Libint2ERI : 0.000u 0.000s 0.005w 1 calls + +************************************************************************************** + +Host: MacBook-Pro-5.local + +Timers On : Sun Apr 19 23:26:35 2026 +Timers Off: Sun Apr 19 23:26:36 2026 + +Wall Time: 0.45 seconds + + Time (seconds) +Module User System Wall Calls +Libint2ERI::Libint2ERI : 0.000u 0.000s 0.006w 1 calls + +-------------------------------------------------------------------------------------- +Libint2ERI::Libint2ERI : 0.000u 0.000s 0.006w 1 calls + +************************************************************************************** + +Host: MacBook-Pro-5.local + +Timers On : Sun Apr 19 23:27:15 2026 +Timers Off: Sun Apr 19 23:27:15 2026 + +Wall Time: 0.36 seconds + + Time (seconds) +Module User System Wall Calls +Libint2ERI::Libint2ERI : 0.000u 0.000s 0.006w 1 calls + +-------------------------------------------------------------------------------------- +Libint2ERI::Libint2ERI : 0.000u 0.000s 0.006w 1 calls + +************************************************************************************** + +Host: MacBook-Pro-5.local + +Timers On : Mon Apr 20 00:02:21 2026 +Timers Off: Mon Apr 20 00:02:30 2026 + +Wall Time: 9.12 seconds + + Time (seconds) +Module User System Wall Calls +Libint2ERI::Libint2ERI : 0.000u 0.017s 0.006w 1 calls + +-------------------------------------------------------------------------------------- +Libint2ERI::Libint2ERI : 0.000u 0.017s 0.006w 1 calls + +************************************************************************************** + +Host: MacBook-Pro-5.local + +Timers On : Mon Apr 20 20:30:41 2026 +Timers Off: Mon Apr 20 20:30:50 2026 + +Wall Time: 9.24 seconds + + Time (seconds) +Module User System Wall Calls +Libint2ERI::Libint2ERI : 0.000u 0.000s 0.006w 1 calls + +-------------------------------------------------------------------------------------- +Libint2ERI::Libint2ERI : 0.000u 0.000s 0.006w 1 calls + +************************************************************************************** + +Host: MacBook-Pro-5.local + +Timers On : Mon Apr 20 20:32:58 2026 +Timers Off: Mon Apr 20 20:33:07 2026 + +Wall Time: 8.90 seconds + + Time (seconds) +Module User System Wall Calls +Libint2ERI::Libint2ERI : 0.000u 0.000s 0.005w 1 calls + +-------------------------------------------------------------------------------------- +Libint2ERI::Libint2ERI : 0.000u 0.000s 0.005w 1 calls + +************************************************************************************** diff --git a/tests/spec_frag_0.in b/tests/spec_frag_0.in deleted file mode 100644 index 9ad9244e..00000000 --- a/tests/spec_frag_0.in +++ /dev/null @@ -1,12 +0,0 @@ -run_type sp -coord xyzabc -terms elec pol disp xr -disp_damp off -fraglib_path ../fraglib - -fragment h2o_l - 0.000 0.000 0.000 0.000 0.000 0.000 -fragment h2o_l - 5.000 0.000 0.000 2.000 0.000 2.000 -fragment nh3_l - 0.000 5.000 0.000 5.000 0.000 -3.000 diff --git a/tests/spec_frag_1.in b/tests/spec_frag_1.in index c8917f9f..0ad2d341 100644 --- a/tests/spec_frag_1.in +++ b/tests/spec_frag_1.in @@ -1,4 +1,5 @@ -run_type sp +run_type etest +ref_energy 0.0000841901 coord xyzabc terms elec pol disp xr disp_damp off diff --git a/tests/spec_frag_2.in b/tests/spec_frag_2.in index 2e0dd0fe..1985c0dc 100644 --- a/tests/spec_frag_2.in +++ b/tests/spec_frag_2.in @@ -1,4 +1,5 @@ -run_type sp +run_type etest +ref_energy -0.0000982356 coord xyzabc terms elec pol disp xr disp_damp off diff --git a/tests/spec_frag_base.in b/tests/spec_frag_base.in index 03db67fb..53973cbf 100644 --- a/tests/spec_frag_base.in +++ b/tests/spec_frag_base.in @@ -1,4 +1,5 @@ -run_type sp +run_type etest +ref_energy -0.0000850092 coord xyzabc terms elec pol disp xr disp_damp off