Uncoupling from forcerec, deleting unused files.
Change-Id: Id72281a26b8b39b6d2a37f5f8a43330f399ce5aa
.. Another useful one-liner to find undocumentedvariables:
.. ( export INPUT_FILE=docs/user-guide/environment-variables.rst; GIT_PAGER="cat "; for ss in `for s in $(git grep getenv | sed 's/.*getenv("\(.*\)".*/\1/' | sort -u | grep '^[A-Z]'); do [ $(grep $s $INPUT_FILE -c) -eq 0 ] && echo $s; done `; do git grep $ss ; done )
-.. todo:: still undocumented GMX_QM_GAUSSIAN_NCPUS
-
Environment Variables
=====================
Analysis and Core Functions
---------------------------
-``GMX_QM_ACCURACY``
- accuracy in Gaussian L510 (MC-SCF) component program.
-
-``GMX_QM_ORCA_BASENAME``
- prefix of :ref:`tpr` files, used in Orca calculations
- for input and output file names.
-
-``GMX_QM_CPMCSCF``
- when set to a nonzero value, Gaussian QM calculations will
- iteratively solve the CP-MCSCF equations.
-
-``GMX_QM_MODIFIED_LINKS_DIR``
- location of modified links in Gaussian.
``DSSP``
used by :ref:`gmx do_dssp` to point to the ``dssp``
executable (not just its path).
-``GMX_QM_GAUSS_DIR``
- directory where Gaussian is installed.
-
-``GMX_QM_GAUSS_EXE``
- name of the Gaussian executable.
-
``GMX_DIPOLE_SPACING``
spacing used by :ref:`gmx dipoles`.
the time unit used in output files, can be
anything in fs, ps, ns, us, ms, s, m or h.
-``GMX_QM_GAUSSIAN_MEMORY``
- memory used for Gaussian QM calculation.
-
``MULTIPROT``
name of the ``multiprot`` executable, used by the
contributed program ``do_multiprot``.
``NCPUS``
number of CPUs to be used for Gaussian QM calculation
-``GMX_ORCA_PATH``
- directory where Orca is installed.
-
-``GMX_QM_SA_STEP``
- simulated annealing step size for Gaussian QM calculation.
-
-``GMX_QM_GROUND_STATE``
- defines state for Gaussian surface hopping calculation.
-
``GMX_TOTAL``
name of the ``total`` executable used by the contributed
``do_shift`` program.
#include "gromacs/mdlib/compute_io.h"
#include "gromacs/mdlib/constr.h"
#include "gromacs/mdlib/perf_est.h"
-#include "gromacs/mdlib/qmmm.h"
#include "gromacs/mdlib/vsite.h"
#include "gromacs/mdrun/mdmodules.h"
#include "gromacs/mdtypes/inputrec.h"
double reppow;
const char* mdparin;
bool bNeedVel, bGenVel;
- gmx_bool have_atomnumber;
gmx_output_env_t* oenv;
gmx_bool bVerbose = FALSE;
warninp* wi;
warning_note(wi, "Temperature coupling is ignored with SD integrators.");
}
- /* If we are doing QM/MM, check that we got the atom numbers */
- have_atomnumber = TRUE;
- for (gmx::index i = 0; i < gmx::ssize(atypes); i++)
- {
- have_atomnumber = have_atomnumber && (atypes.atomNumberFromAtomType(i) >= 0);
- }
- if (!have_atomnumber && ir->bQMMM)
- {
- warning_error(
- wi,
- "\n"
- "It appears as if you are trying to run a QM/MM calculation, but the force\n"
- "field you are using does not contain atom numbers fields. This is an\n"
- "optional field (introduced in GROMACS 3.3) for general runs, but mandatory\n"
- "for QM/MM. The good news is that it is easy to add - put the atom number as\n"
- "an integer just before the mass column in ffXXXnb.itp.\n"
- "NB: United atoms have the same atom numbers as normal ones.\n\n");
- }
-
/* Check for errors in the input now, since they might cause problems
* during processing further down.
*/
#include "gromacs/math/vec.h"
#include "gromacs/math/vecdump.h"
#include "gromacs/mdlib/forcerec_threading.h"
-#include "gromacs/mdlib/qmmm.h"
#include "gromacs/mdlib/rf_util.h"
#include "gromacs/mdlib/wall.h"
#include "gromacs/mdtypes/commrec.h"
auto& forceWithVirial = forceOutputs->forceWithVirial();
- /* do QMMM first if requested */
- if (fr->bQMMM)
- {
- enerd->term[F_EQM] = calculate_QMMM(cr, &forceOutputs->forceWithShiftForces(), fr->qr);
- }
-
/* Call the short range functions all in one go. */
if (ir->nwall)
#include "gromacs/mdlib/forcerec_threading.h"
#include "gromacs/mdlib/gmx_omp_nthreads.h"
#include "gromacs/mdlib/md_support.h"
-#include "gromacs/mdlib/qmmm.h"
#include "gromacs/mdlib/rf_util.h"
#include "gromacs/mdlib/wall.h"
#include "gromacs/mdtypes/commrec.h"
}
// QM/MM initialization if requested
- fr->bQMMM = ir->bQMMM;
- if (fr->bQMMM)
+ if (ir->bQMMM)
{
- // Initialize QM/MM if supported
- if (GMX_QMMM)
- {
- GMX_LOG(mdlog.info)
- .asParagraph()
- .appendText(
- "Large parts of the QM/MM support is deprecated, and may be removed in "
- "a future "
- "version. Please get in touch with the developers if you find the "
- "support useful, "
- "as help is needed if the functionality is to continue to be "
- "available.");
- fr->qr = mk_QMMMrec();
- init_QMMMrec(cr, mtop, ir, fr);
- }
- else
- {
- gmx_incons(
- "QM/MM was requested, but is only available when GROMACS "
- "is configured with QM/MM support");
- }
+
+ gmx_incons("QM/MM was requested, but is no longer available in GROMACS");
}
/* Set all the static charge group info */
#include "gromacs/gpu_utils/hostallocator.h"
#include "gromacs/math/functions.h"
#include "gromacs/mdlib/gmx_omp_nthreads.h"
-#include "gromacs/mdlib/qmmm.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/mdtypes/mdatom.h"
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2017,2018 by the GROMACS development team.
- * Copyright (c) 2019,2020, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-#include "gmxpre.h"
-
-#include "qm_gamess.h"
-
-#include "config.h"
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-
-#include <cmath>
-
-#include "gromacs/fileio/confio.h"
-#include "gromacs/gmxlib/network.h"
-#include "gromacs/gmxlib/nrnb.h"
-#include "gromacs/math/units.h"
-#include "gromacs/math/vec.h"
-#include "gromacs/mdlib/qmmm.h"
-#include "gromacs/mdtypes/commrec.h"
-#include "gromacs/mdtypes/forcerec.h"
-#include "gromacs/mdtypes/md_enums.h"
-#include "gromacs/utility/fatalerror.h"
-#include "gromacs/utility/smalloc.h"
-
-// When not built in a configuration with QMMM support, much of this
-// code is unreachable by design. Tell clang not to warn about it.
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wmissing-noreturn"
-
-/* QMMM sub routines */
-/* mopac interface routines */
-
-
-static void F77_FUNC(inigms, IMIGMS)();
-
-static void F77_FUNC(grads, GRADS)(const int* nrqmat,
- real* qmcrd,
- const int* nrmmat,
- const real* mmchrg,
- real* mmcrd,
- real* qmgrad,
- real* mmgrad,
- real* energy);
-
-#if !GMX_QMMM_GAMESS
-// Stub definitions to make compilation succeed when not configured
-// for GAMESS support. In that case, the module gives a fatal error
-// when the initialization function is called, so there is no need to
-// issue fatal errors here, because that introduces problems with
-// tools suggesting and prohibiting noreturn attributes.
-
-void F77_FUNC(inigms, IMIGMS)(){};
-// NOLINTNEXTLINE(readability-named-parameter)
-void F77_FUNC(grads, GRADS)(const int*, real*, const int*, const real*, real*, real*, real*, real*){};
-#endif
-
-void init_gamess(const t_commrec* cr, t_QMrec* qm, t_MMrec* mm)
-{
- /* it works hopelessly complicated :-)
- * first a file is written. Then the standard gamess input/output
- * routine is called (no system()!) to set up all fortran arrays.
- * this routine writes a punch file, like in a normal gamess run.
- * via this punch file the other games routines, needed for gradient
- * and energy evaluations are called. This setup works fine for
- * dynamics simulations. 7-6-2002 (London)
- */
- int i, j;
- FILE* out;
- char periodic_system[37][3] = { "XX", "H ", "He", "Li", "Be", "B ", "C ", "N ", "O ", "F ",
- "Ne", "Na", "Mg", "Al", "Si", "P ", "S ", "Cl", "Ar", "K ",
- "Ca", "Sc", "Ti", "V ", "Cr", "Mn", "Fe", "Co", "Ni", "Cu",
- "Zn", "Ga", "Ge", "As", "Se", "Br", "Kr" };
- if (!GMX_QMMM_GAMESS)
- {
- gmx_fatal(FARGS,
- "Cannot call GAMESS unless linked against it. Use cmake "
- "-DGMX_QMMM_PROGRAM=GAMESS, and ensure that linking will work correctly.");
- }
-
- if (PAR(cr))
- {
-
- if (MASTER(cr))
- {
- out = fopen("FOR009", "w");
- /* of these options I am not completely sure.... the overall
- * preformance on more than 4 cpu's is rather poor at the moment.
- */
- fprintf(out, "memory 48000000\nPARALLEL IOMODE SCREENED\n");
- fprintf(out, "ELEC %d\nMULT %d\nSUPER ON\nNOSYM\nGEOMETRY ANGSTROM\n", qm->nelectrons,
- qm->multiplicity);
- for (i = 0; i < qm->nrQMatoms; i++)
- {
-#ifdef DOUBLE
- fprintf(out, "%10.7lf %10.7lf %10.7lf %5.3lf %2s\n", i / 2., i / 3., i / 4.,
- qm->atomicnumberQM[i] * 1.0, periodic_system[qm->atomicnumberQM[i]]);
-#else
- fprintf(out, "%10.7f %10.7f %10.7f %5.3f %2s\n", i / 2., i / 3., i / 4.,
- qm->atomicnumberQM[i] * 1.0, periodic_system[qm->atomicnumberQM[i]]);
-#endif
- }
- if (mm->nrMMatoms)
- {
- for (j = i; j < i + 2; j++)
- {
-#ifdef DOUBLE
- fprintf(out, "%10.7lf %10.7lf %10.7lf %5.3lf BQ\n", j / 5., j / 6., j / 7., 1.0);
-#else
- fprintf(out, "%10.7f %10.7f %10.7f %5.3f BQ\n", j / 5., j / 6., j / 7., 2.0);
-#endif
- }
- }
- fprintf(out, "END\nBASIS %s\nRUNTYPE GRADIENT\nSCFTYPE %s\n",
- eQMbasis_names[qm->QMbasis], eQMmethod_names[qm->QMmethod]); /* see enum.h */
- fclose(out);
- }
- gmx_barrier(cr);
- F77_FUNC(inigms, IMIGMS)();
- }
- else /* normal serial run */
-
- {
- out = fopen("FOR009", "w");
- /* of these options I am not completely sure.... the overall
- * preformance on more than 4 cpu's is rather poor at the moment.
- */
- fprintf(out, "ELEC %d\nMULT %d\nSUPER ON\nNOSYM\nGEOMETRY ANGSTROM\n", qm->nelectrons,
- qm->multiplicity);
- for (i = 0; i < qm->nrQMatoms; i++)
- {
-#ifdef DOUBLE
- fprintf(out, "%10.7lf %10.7lf %10.7lf %5.3lf %2s\n", i / 2., i / 3., i / 4.,
- qm->atomicnumberQM[i] * 1.0, periodic_system[qm->atomicnumberQM[i]]);
-#else
- fprintf(out, "%10.7f %10.7f %10.7f %5.3f %2s\n", i / 2., i / 3., i / 4.,
- qm->atomicnumberQM[i] * 1.0, periodic_system[qm->atomicnumberQM[i]]);
-#endif
- }
- if (mm->nrMMatoms)
- {
- for (j = i; j < i + 2; j++)
- {
-#ifdef DOUBLE
- fprintf(out, "%10.7lf %10.7lf %10.7lf %5.3lf BQ\n", j / 5., j / 6., j / 7., 1.0);
-#else
- fprintf(out, "%10.7f %10.7f %10.7f %5.3f BQ\n", j / 5., j / 6., j / 7., 2.0);
-#endif
- }
- }
- fprintf(out, "END\nBASIS %s\nRUNTYPE GRADIENT\nSCFTYPE %s\n", eQMbasis_names[qm->QMbasis],
- eQMmethod_names[qm->QMmethod]); /* see enum.h */
- F77_FUNC(inigms, IMIGMS)();
- }
-}
-
-real call_gamess(const t_QMrec* qm, const t_MMrec* mm, rvec f[], rvec fshift[])
-{
- /* do the actual QMMM calculation using GAMESS-UK. In this
- * implementation (3-2001) a system call is made to the GAMESS-UK
- * binary. Now we are working to get the electron integral, SCF, and
- * gradient routines linked directly
- */
- int i, j;
- real QMener = 0.0, *qmgrad, *mmgrad, *mmcrd, *qmcrd, energy = 0;
-
- snew(qmcrd, 3 * (qm->nrQMatoms));
- snew(mmcrd, 3 * (mm->nrMMatoms));
- snew(qmgrad, 3 * (qm->nrQMatoms));
- snew(mmgrad, 3 * (mm->nrMMatoms));
-
- /* copy the data from qr into the arrays that are going to be used
- * in the fortran routines of gamess
- */
- for (i = 0; i < qm->nrQMatoms; i++)
- {
- for (j = 0; j < DIM; j++)
- {
- qmcrd[DIM * i + j] = 1 / BOHR2NM * qm->xQM[i][j];
- }
- }
- for (i = 0; i < mm->nrMMatoms; i++)
- {
- for (j = 0; j < DIM; j++)
- {
- mmcrd[DIM * i + j] = 1 / BOHR2NM * mm->xMM[i][j];
- }
- }
- for (i = 0; i < 3 * qm->nrQMatoms; i += 3)
- {
- fprintf(stderr, "%8.5f, %8.5f, %8.5f\n", qmcrd[i], qmcrd[i + 1], qmcrd[i + 2]);
- }
-
- F77_FUNC(grads, GRADS)
- (&qm->nrQMatoms, qmcrd, &mm->nrMMatoms, mm->MMcharges, mmcrd, qmgrad, mmgrad, &energy);
-
- for (i = 0; i < qm->nrQMatoms; i++)
- {
- for (j = 0; j < DIM; j++)
- {
- f[i][j] = HARTREE_BOHR2MD * qmgrad[3 * i + j];
- fshift[i][j] = HARTREE_BOHR2MD * qmgrad[3 * i + j];
- }
- }
- for (i = 0; i < mm->nrMMatoms; i++)
- {
- for (j = 0; j < DIM; j++)
- {
- f[i][j] = HARTREE_BOHR2MD * mmgrad[3 * i + j];
- fshift[i][j] = HARTREE_BOHR2MD * mmgrad[3 * i + j];
- }
- }
- /* convert a.u to kJ/mol */
- QMener = energy * HARTREE2KJ * AVOGADRO;
- return (QMener);
-}
-
-#pragma GCC diagnostic pop
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 2018,2019, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-#ifndef GMX_MDLIB_QMGAMESS_H
-#define GMX_MDLIB_QMGAMESS_H
-
-#include "gromacs/mdlib/qmmm.h"
-
-/*! \brief
- * Initialize gamess datastructures.
- *
- * \param[in] cr Commrec datastructure.
- * \param[in] qm QM forcerec.
- * \param[in] mm MM part of forcerec.
- */
-void init_gamess(const t_commrec* cr, t_QMrec* qm, t_MMrec* mm);
-
-/*! \brief
- * Run calculation with Gamess.
- *
- * \param[in] qm QM part of forcerec.
- * \param[in] mm MM part of forcerec.
- * \param[in] f Force vector.
- * \param[in] fshift Force shift vector.
- */
-real call_gamess(const t_QMrec* qm, const t_MMrec* mm, rvec f[], rvec fshift[]);
-
-
-#endif
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-#include "gmxpre.h"
-
-#include "qm_gaussian.h"
-
-#include "config.h"
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-
-#include <cmath>
-
-#include "gromacs/fileio/confio.h"
-#include "gromacs/gmxlib/network.h"
-#include "gromacs/gmxlib/nrnb.h"
-#include "gromacs/math/units.h"
-#include "gromacs/math/vec.h"
-#include "gromacs/mdlib/force.h"
-#include "gromacs/mdlib/qmmm.h"
-#include "gromacs/mdtypes/md_enums.h"
-#include "gromacs/utility/cstringutil.h"
-#include "gromacs/utility/fatalerror.h"
-#include "gromacs/utility/smalloc.h"
-
-// When not built in a configuration with QMMM support, much of this
-// code is unreachable by design. Tell clang not to warn about it.
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wmissing-noreturn"
-
-/* TODO: this should be made thread-safe */
-
-/* Gaussian interface routines */
-
-void init_gaussian(t_QMrec* qm)
-{
- ivec basissets[eQMbasisNR] = { { 0, 3, 0 }, { 0, 3, 0 }, /*added for double sto-3g entry in names.c*/
- { 5, 0, 0 }, { 5, 0, 1 }, { 5, 0, 11 }, { 5, 6, 0 },
- { 1, 6, 0 }, { 1, 6, 1 }, { 1, 6, 11 }, { 4, 6, 0 } };
- char* buf = nullptr;
- int i;
-
- if (!GMX_QMMM_GAUSSIAN)
- {
- gmx_fatal(FARGS,
- "Cannot call GAUSSIAN unless linked against it. Use cmake "
- "-DGMX_QMMM_PROGRAM=GAUSSIAN, and ensure that linking will work correctly.");
- }
-
- /* using the ivec above to convert the basis read form the mdp file
- * in a human readable format into some numbers for the gaussian
- * route. This is necessary as we are using non standard routes to
- * do SH.
- */
-
- /* per layer we make a new subdir for integral file, checkpoint
- * files and such. These dirs are stored in the QMrec for
- * convenience
- */
-
-
- if (!qm->nQMcpus) /* this we do only once per layer
- * as we call g01 externally
- */
-
- {
- for (i = 0; i < DIM; i++)
- {
- qm->SHbasis[i] = basissets[qm->QMbasis][i];
- }
-
- /* init gradually switching on of the SA */
- qm->SAstep = 0;
- /* we read the number of cpus and environment from the environment
- * if set.
- */
- buf = getenv("GMX_QM_GAUSSIAN_NCPUS");
- if (buf)
- {
- sscanf(buf, "%d", &qm->nQMcpus);
- }
- else
- {
- qm->nQMcpus = 1;
- }
- fprintf(stderr, "number of CPUs for gaussian = %d\n", qm->nQMcpus);
- buf = getenv("GMX_QM_GAUSSIAN_MEMORY");
- if (buf)
- {
- sscanf(buf, "%d", &qm->QMmem);
- }
- else
- {
- qm->QMmem = 50000000;
- }
- fprintf(stderr, "memory for gaussian = %d\n", qm->QMmem);
- buf = getenv("GMX_QM_ACCURACY");
- if (buf)
- {
- sscanf(buf, "%d", &qm->accuracy);
- }
- else
- {
- qm->accuracy = 8;
- }
- fprintf(stderr, "accuracy in l510 = %d\n", qm->accuracy);
-
- buf = getenv("GMX_QM_CPMCSCF");
- if (buf)
- {
- sscanf(buf, "%d", &i);
- qm->cpmcscf = (i != 0);
- }
- else
- {
- qm->cpmcscf = FALSE;
- }
- if (qm->cpmcscf)
- {
- fprintf(stderr, "using cp-mcscf in l1003\n");
- }
- else
- {
- fprintf(stderr, "NOT using cp-mcscf in l1003\n");
- }
- buf = getenv("GMX_QM_SA_STEP");
- if (buf)
- {
- sscanf(buf, "%d", &qm->SAstep);
- }
- else
- {
- /* init gradually switching on of the SA */
- qm->SAstep = 0;
- }
- /* we read the number of cpus and environment from the environment
- * if set.
- */
- fprintf(stderr, "Level of SA at start = %d\n", qm->SAstep);
- /* gaussian settings on the system */
- buf = getenv("GMX_QM_GAUSS_DIR");
-
- if (buf)
- {
- qm->gauss_dir = gmx_strdup(buf);
- }
- else
- {
- gmx_fatal(FARGS, "no $GMX_QM_GAUSS_DIR, check gaussian manual\n");
- }
-
- buf = getenv("GMX_QM_GAUSS_EXE");
- if (buf)
- {
- qm->gauss_exe = gmx_strdup(buf);
- }
- else
- {
- gmx_fatal(FARGS, "no $GMX_QM_GAUSS_EXE set, check gaussian manual\n");
- }
- buf = getenv("GMX_QM_MODIFIED_LINKS_DIR");
- if (buf)
- {
- qm->devel_dir = gmx_strdup(buf);
- }
- else
- {
- gmx_fatal(FARGS,
- "no $GMX_QM_MODIFIED_LINKS_DIR, this is were the modified links reside.\n");
- }
-
- /* if(fr->bRF){*/
- /* reactionfield, file is needed using gaussian */
- /* rffile=fopen("rf.dat","w");*/
- /* fprintf(rffile,"%f %f\n",fr->epsilon_r,fr->rcoulomb/BOHR2NM);*/
- /* fclose(rffile);*/
- /* }*/
- }
- fprintf(stderr, "gaussian initialised...\n");
-}
-
-
-static void write_gaussian_SH_input(int step, gmx_bool swap, const t_QMMMrec* QMMMrec, t_QMrec* qm, t_MMrec* mm)
-{
- int i;
- bool bSA = (qm->SAstep > 0);
- FILE* out = fopen("input.com", "w");
- /* write the route */
- fprintf(out, "%s", "%scr=input\n");
- fprintf(out, "%s", "%rwf=input\n");
- fprintf(out, "%s", "%int=input\n");
- fprintf(out, "%s", "%d2e=input\n");
- /* if(step)
- * fprintf(out,"%s","%nosave\n");
- */
- fprintf(out, "%s", "%chk=input\n");
- fprintf(out, "%s%d\n", "%mem=", qm->QMmem);
- fprintf(out, "%s%3d\n", "%nprocshare=", qm->nQMcpus);
-
- /* use the versions of
- * l301 that computes the interaction between MM and QM atoms.
- * l510 that can punch the CI coefficients
- * l701 that can do gradients on MM atoms
- */
-
- /* local version */
- fprintf(out, "%s%s%s", "%subst l510 ", qm->devel_dir, "/l510\n");
- fprintf(out, "%s%s%s", "%subst l301 ", qm->devel_dir, "/l301\n");
- fprintf(out, "%s%s%s", "%subst l701 ", qm->devel_dir, "/l701\n");
-
- fprintf(out, "%s%s%s", "%subst l1003 ", qm->devel_dir, "/l1003\n");
- fprintf(out, "%s%s%s", "%subst l9999 ", qm->devel_dir, "/l9999\n");
- /* print the nonstandard route
- */
- fprintf(out, "%s", "#P nonstd\n 1/18=10,20=1,38=1/1;\n");
- fprintf(out, "%s", " 2/9=110,15=1,17=6,18=5,40=1/2;\n");
- if (mm->nrMMatoms)
- {
- fprintf(out, " 3/5=%d,6=%d,7=%d,25=1,32=1,43=1,94=-2/1,2,3;\n", qm->SHbasis[0],
- qm->SHbasis[1], qm->SHbasis[2]); /*basisset stuff */
- }
- else
- {
- fprintf(out, " 3/5=%d,6=%d,7=%d,25=1,32=1,43=0,94=-2/1,2,3;\n", qm->SHbasis[0],
- qm->SHbasis[1], qm->SHbasis[2]); /*basisset stuff */
- }
- /* development */
- if (step + 1) /* fetch initial guess from check point file */
- { /* hack, to alyays read from chk file!!!!! */
- fprintf(out, "%s%d,%s%d%s", " 4/5=1,7=6,17=", qm->CASelectrons, "18=", qm->CASorbitals,
- "/1,5;\n");
- }
- else /* generate the first checkpoint file */
- {
- fprintf(out, "%s%d,%s%d%s", " 4/5=0,7=6,17=", qm->CASelectrons, "18=", qm->CASorbitals,
- "/1,5;\n");
- }
- /* the rest of the input depends on where the system is on the PES
- */
- if (swap && bSA) /* make a slide to the other surface */
- {
- if (qm->CASorbitals > 6) /* use direct and no full diag */
- {
- fprintf(out, " 5/5=2,16=-2,17=10000000,28=2,32=2,38=6,97=100/10;\n");
- }
- else
- {
- if (qm->cpmcscf)
- {
- fprintf(out, " 5/5=2,6=%d,17=31000200,28=2,32=2,38=6,97=100/10;\n", qm->accuracy);
- if (mm->nrMMatoms > 0)
- {
- fprintf(out, " 7/7=1,16=-2,30=1/1;\n");
- }
- fprintf(out, " 11/31=1,42=1,45=1/1;\n");
- fprintf(out, " 10/6=1,10=700006,28=2,29=1,31=1,97=100/3;\n");
- fprintf(out, " 7/30=1/16;\n 99/10=4/99;\n");
- }
- else
- {
- fprintf(out, " 5/5=2,6=%d,17=11000000,28=2,32=2,38=6,97=100/10;\n", qm->accuracy);
- fprintf(out, " 7/7=1,16=-2,30=1/1,2,3,16;\n 99/10=4/99;\n");
- }
- }
- }
- else if (bSA) /* do a "state-averaged" CAS calculation */
- {
- if (qm->CASorbitals > 6) /* no full diag */
- {
- fprintf(out, " 5/5=2,16=-2,17=10000000,28=2,32=2,38=6/10;\n");
- }
- else
- {
- if (qm->cpmcscf)
- {
- fprintf(out, " 5/5=2,6=%d,17=31000200,28=2,32=2,38=6/10;\n", qm->accuracy);
- if (mm->nrMMatoms > 0)
- {
- fprintf(out, " 7/7=1,16=-2,30=1/1;\n");
- }
- fprintf(out, " 11/31=1,42=1,45=1/1;\n");
- fprintf(out, " 10/6=1,10=700006,28=2,29=1,31=1/3;\n");
- fprintf(out, " 7/30=1/16;\n 99/10=4/99;\n");
- }
- else
- {
- fprintf(out, " 5/5=2,6=%d,17=11000000,28=2,32=2,38=6/10;\n", qm->accuracy);
- fprintf(out, " 7/7=1,16=-2,30=1/1,2,3,16;\n 99/10=4/99;\n");
- }
- }
- }
- else if (swap) /* do a "swapped" CAS calculation */
- {
- if (qm->CASorbitals > 6)
- {
- fprintf(out, " 5/5=2,16=-2,17=0,28=2,32=2,38=6,97=100/10;\n");
- }
- else
- {
- fprintf(out, " 5/5=2,6=%d,17=1000000,28=2,32=2,38=6,97=100/10;\n", qm->accuracy);
- }
- fprintf(out, " 7/7=1,16=-2,30=1/1,2,3,16;\n 99/10=4/99;\n");
- }
- else /* do a "normal" CAS calculation */
- {
- if (qm->CASorbitals > 6)
- {
- fprintf(out, " 5/5=2,16=-2,17=0,28=2,32=2,38=6/10;\n");
- }
- else
- {
- fprintf(out, " 5/5=2,6=%d,17=1000000,28=2,32=2,38=6/10;\n", qm->accuracy);
- }
- fprintf(out, " 7/7=1,16=-2,30=1/1,2,3,16;\n 99/10=4/99;\n");
- }
- fprintf(out, "\ninput-file generated by gromacs\n\n");
- fprintf(out, "%2d%2d\n", qm->QMcharge, qm->multiplicity);
- for (i = 0; i < qm->nrQMatoms; i++)
- {
- fprintf(out, "%3d %10.7f %10.7f %10.7f\n", qm->atomicnumberQM[i],
- qm->xQM[i][XX] / BOHR2NM, qm->xQM[i][YY] / BOHR2NM, qm->xQM[i][ZZ] / BOHR2NM);
- }
- /* MM point charge data */
- if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms)
- {
- fprintf(out, "\n");
- for (i = 0; i < mm->nrMMatoms; i++)
- {
- fprintf(out, "%10.7f %10.7f %10.7f %8.4f\n", mm->xMM[i][XX] / BOHR2NM,
- mm->xMM[i][YY] / BOHR2NM, mm->xMM[i][ZZ] / BOHR2NM, mm->MMcharges[i]);
- }
- }
- if (bSA) /* put the SA coefficients at the end of the file */
- {
- fprintf(out, "\n%10.8f %10.8f\n", qm->SAstep * 0.5 / qm->SAsteps, 1 - qm->SAstep * 0.5 / qm->SAsteps);
- fprintf(stderr, "State Averaging level = %d/%d\n", qm->SAstep, qm->SAsteps);
- }
- fprintf(out, "\n");
- fclose(out);
-} /* write_gaussian_SH_input */
-
-static void write_gaussian_input(int step, const t_QMMMrec* QMMMrec, t_QMrec* qm, t_MMrec* mm)
-{
- int i;
-
- FILE* out = fopen("input.com", "w");
- /* write the route */
-
- if (qm->QMmethod >= eQMmethodRHF)
- {
- fprintf(out, "%s", "%chk=input\n");
- }
- else
- {
- fprintf(out, "%s", "%chk=se\n");
- }
- if (qm->nQMcpus > 1)
- {
- fprintf(out, "%s%3d\n", "%nprocshare=", qm->nQMcpus);
- }
- fprintf(out, "%s%d\n", "%mem=", qm->QMmem);
- fprintf(out, "%s%s%s", "%subst l701 ", qm->devel_dir, "/l701\n");
- fprintf(out, "%s%s%s", "%subst l301 ", qm->devel_dir, "/l301\n");
- fprintf(out, "%s%s%s", "%subst l9999 ", qm->devel_dir, "/l9999\n");
- if (step)
- {
- fprintf(out, "%s", "#T ");
- }
- else
- {
- fprintf(out, "%s", "#P ");
- }
- if (qm->QMmethod == eQMmethodB3LYPLAN)
- {
- fprintf(out, " %s", "B3LYP/GEN Pseudo=Read");
- }
- else
- {
- fprintf(out, " %s", eQMmethod_names[qm->QMmethod]);
-
- if (qm->QMmethod >= eQMmethodRHF)
- {
- if (qm->QMmethod == eQMmethodCASSCF)
- {
- /* in case of cas, how many electrons and orbitals do we need?
- */
- fprintf(out, "(%d,%d)", qm->CASelectrons, qm->CASorbitals);
- }
- fprintf(out, "/%s", eQMbasis_names[qm->QMbasis]);
- }
- }
- if (QMMMrec->QMMMscheme == eQMMMschemenormal && mm->nrMMatoms)
- {
- fprintf(out, " %s", "Charge ");
- }
- if (step || qm->QMmethod == eQMmethodCASSCF)
- {
- /* fetch guess from checkpoint file, always for CASSCF */
- fprintf(out, "%s", " guess=read");
- }
- fprintf(out, "\nNosymm units=bohr\n");
-
- fprintf(out, "FORCE Punch=(Derivatives) ");
- fprintf(out, "iop(3/33=1)\n\n");
- fprintf(out, "input-file generated by gromacs\n\n");
- fprintf(out, "%2d%2d\n", qm->QMcharge, qm->multiplicity);
- for (i = 0; i < qm->nrQMatoms; i++)
- {
- fprintf(out, "%3d %10.7f %10.7f %10.7f\n", qm->atomicnumberQM[i],
- qm->xQM[i][XX] / BOHR2NM, qm->xQM[i][YY] / BOHR2NM, qm->xQM[i][ZZ] / BOHR2NM);
- }
-
- /* Pseudo Potential and ECP are included here if selected (MEthod suffix LAN) */
- if (qm->QMmethod == eQMmethodB3LYPLAN)
- {
- fprintf(out, "\n");
- for (i = 0; i < qm->nrQMatoms; i++)
- {
- if (qm->atomicnumberQM[i] < 21)
- {
- fprintf(out, "%d ", i + 1);
- }
- }
- fprintf(out, "\n%s\n****\n", eQMbasis_names[qm->QMbasis]);
-
- for (i = 0; i < qm->nrQMatoms; i++)
- {
- if (qm->atomicnumberQM[i] > 21)
- {
- fprintf(out, "%d ", i + 1);
- }
- }
- fprintf(out, "\n%s\n****\n\n", "lanl2dz");
-
- for (i = 0; i < qm->nrQMatoms; i++)
- {
- if (qm->atomicnumberQM[i] > 21)
- {
- fprintf(out, "%d ", i + 1);
- }
- }
- fprintf(out, "\n%s\n", "lanl2dz");
- }
-
-
- /* MM point charge data */
- if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms)
- {
- fprintf(stderr, "nr mm atoms in gaussian.c = %d\n", mm->nrMMatoms);
- fprintf(out, "\n");
- for (i = 0; i < mm->nrMMatoms; i++)
- {
- fprintf(out, "%10.7f %10.7f %10.7f %8.4f\n", mm->xMM[i][XX] / BOHR2NM,
- mm->xMM[i][YY] / BOHR2NM, mm->xMM[i][ZZ] / BOHR2NM, mm->MMcharges[i]);
- }
- }
- fprintf(out, "\n");
-
-
- fclose(out);
-
-} /* write_gaussian_input */
-
-static real read_gaussian_output(rvec QMgrad[], rvec MMgrad[], t_QMrec* qm, t_MMrec* mm)
-{
- int i;
- char buf[300];
- real QMener;
- FILE* in;
-
- in = fopen("fort.7", "r");
-
- /* (There was additional content in the file in case
- * of QM optimizations / transition state search,
- * which was removed.
- */
- /* the next line is the energy and in the case of CAS, the energy
- * difference between the two states.
- */
- if (nullptr == fgets(buf, 300, in))
- {
- gmx_fatal(FARGS, "Error reading Gaussian output");
- }
-
-#if GMX_DOUBLE
- sscanf(buf, "%lf\n", &QMener);
-#else
- sscanf(buf, "%f\n", &QMener);
-#endif
- /* next lines contain the gradients of the QM atoms */
- for (i = 0; i < qm->nrQMatoms; i++)
- {
- if (nullptr == fgets(buf, 300, in))
- {
- gmx_fatal(FARGS, "Error reading Gaussian output");
- }
-#if GMX_DOUBLE
- sscanf(buf, "%lf %lf %lf\n", &QMgrad[i][XX], &QMgrad[i][YY], &QMgrad[i][ZZ]);
-#else
- sscanf(buf, "%f %f %f\n", &QMgrad[i][XX], &QMgrad[i][YY], &QMgrad[i][ZZ]);
-#endif
- }
- /* the next lines are the gradients of the MM atoms */
- if (qm->QMmethod >= eQMmethodRHF)
- {
- for (i = 0; i < mm->nrMMatoms; i++)
- {
- if (nullptr == fgets(buf, 300, in))
- {
- gmx_fatal(FARGS, "Error reading Gaussian output");
- }
-#if GMX_DOUBLE
- sscanf(buf, "%lf %lf %lf\n", &MMgrad[i][XX], &MMgrad[i][YY], &MMgrad[i][ZZ]);
-#else
- sscanf(buf, "%f %f %f\n", &MMgrad[i][XX], &MMgrad[i][YY], &MMgrad[i][ZZ]);
-#endif
- }
- }
- fclose(in);
- return (QMener);
-}
-
-static real read_gaussian_SH_output(rvec QMgrad[], rvec MMgrad[], int step, t_QMrec* qm, t_MMrec* mm)
-{
- int i;
- char buf[300];
- real QMener, DeltaE;
- FILE* in;
-
- in = fopen("fort.7", "r");
- /* first line is the energy and in the case of CAS, the energy
- * difference between the two states.
- */
- if (nullptr == fgets(buf, 300, in))
- {
- gmx_fatal(FARGS, "Error reading Gaussian output");
- }
-
-#if GMX_DOUBLE
- sscanf(buf, "%lf %lf\n", &QMener, &DeltaE);
-#else
- sscanf(buf, "%f %f\n", &QMener, &DeltaE);
-#endif
-
- /* switch on/off the State Averaging */
-
- if (DeltaE > qm->SAoff)
- {
- if (qm->SAstep > 0)
- {
- qm->SAstep--;
- }
- }
- else if (DeltaE < qm->SAon || (qm->SAstep > 0))
- {
- if (qm->SAstep < qm->SAsteps)
- {
- qm->SAstep++;
- }
- }
-
- /* for debugging: */
- fprintf(stderr, "Gap = %5f,SA = %3d\n", DeltaE, static_cast<int>(qm->SAstep > 0));
- /* next lines contain the gradients of the QM atoms */
- for (i = 0; i < qm->nrQMatoms; i++)
- {
- if (nullptr == fgets(buf, 300, in))
- {
- gmx_fatal(FARGS, "Error reading Gaussian output");
- }
-
-#if GMX_DOUBLE
- sscanf(buf, "%lf %lf %lf\n", &QMgrad[i][XX], &QMgrad[i][YY], &QMgrad[i][ZZ]);
-#else
- sscanf(buf, "%f %f %f\n", &QMgrad[i][XX], &QMgrad[i][YY], &QMgrad[i][ZZ]);
-#endif
- }
- /* the next lines, are the gradients of the MM atoms */
-
- for (i = 0; i < mm->nrMMatoms; i++)
- {
- if (nullptr == fgets(buf, 300, in))
- {
- gmx_fatal(FARGS, "Error reading Gaussian output");
- }
-#if GMX_DOUBLE
- sscanf(buf, "%lf %lf %lf\n", &MMgrad[i][XX], &MMgrad[i][YY], &MMgrad[i][ZZ]);
-#else
- sscanf(buf, "%f %f %f\n", &MMgrad[i][XX], &MMgrad[i][YY], &MMgrad[i][ZZ]);
-#endif
- }
-
- /* the next line contains the two CI eigenvector elements */
- if (nullptr == fgets(buf, 300, in))
- {
- gmx_fatal(FARGS, "Error reading Gaussian output");
- }
- if (!step)
- {
- sscanf(buf, "%d", &qm->CIdim);
- snew(qm->CIvec1, qm->CIdim);
- snew(qm->CIvec1old, qm->CIdim);
- snew(qm->CIvec2, qm->CIdim);
- snew(qm->CIvec2old, qm->CIdim);
- }
- else
- {
- /* before reading in the new current CI vectors, copy the current
- * CI vector into the old one.
- */
- for (i = 0; i < qm->CIdim; i++)
- {
- qm->CIvec1old[i] = qm->CIvec1[i];
- qm->CIvec2old[i] = qm->CIvec2[i];
- }
- }
- /* first vector */
- for (i = 0; i < qm->CIdim; i++)
- {
- if (nullptr == fgets(buf, 300, in))
- {
- gmx_fatal(FARGS, "Error reading Gaussian output");
- }
-#if GMX_DOUBLE
- sscanf(buf, "%lf\n", &qm->CIvec1[i]);
-#else
- sscanf(buf, "%f\n", &qm->CIvec1[i]);
-#endif
- }
- /* second vector */
- for (i = 0; i < qm->CIdim; i++)
- {
- if (nullptr == fgets(buf, 300, in))
- {
- gmx_fatal(FARGS, "Error reading Gaussian output");
- }
-#if GMX_DOUBLE
- sscanf(buf, "%lf\n", &qm->CIvec2[i]);
-#else
- sscanf(buf, "%f\n", &qm->CIvec2[i]);
-#endif
- }
- fclose(in);
- return (QMener);
-}
-
-static real inproduct(const real* a, const real* b, int n)
-{
- int i;
- real dot = 0.0;
-
- /* computes the inner product between two vectors (a.b), both of
- * which have length n.
- */
- for (i = 0; i < n; i++)
- {
- dot += a[i] * b[i];
- }
- return (dot);
-}
-
-static int hop(int step, t_QMrec* qm)
-{
- int swap = 0;
- real d11 = 0.0, d12 = 0.0, d21 = 0.0, d22 = 0.0;
-
- /* calculates the inproduct between the current Ci vector and the
- * previous CI vector. A diabatic hop will be made if d12 and d21
- * are much bigger than d11 and d22. In that case hop returns true,
- * otherwise it returns false.
- */
- if (step) /* only go on if more than one step has been done */
- {
- d11 = inproduct(qm->CIvec1, qm->CIvec1old, qm->CIdim);
- d12 = inproduct(qm->CIvec1, qm->CIvec2old, qm->CIdim);
- d21 = inproduct(qm->CIvec2, qm->CIvec1old, qm->CIdim);
- d22 = inproduct(qm->CIvec2, qm->CIvec2old, qm->CIdim);
- }
- fprintf(stderr, "-------------------\n");
- fprintf(stderr, "d11 = %13.8f\n", d11);
- fprintf(stderr, "d12 = %13.8f\n", d12);
- fprintf(stderr, "d21 = %13.8f\n", d21);
- fprintf(stderr, "d22 = %13.8f\n", d22);
- fprintf(stderr, "-------------------\n");
-
- if ((fabs(d12) > 0.5) && (fabs(d21) > 0.5))
- {
- swap = 1;
- }
-
- return (swap);
-}
-
-static void do_gaussian(int step, char* exe)
-{
- char buf[STRLEN];
-
- /* make the call to the gaussian binary through system()
- * The location of the binary will be picked up from the
- * environment using getenv().
- */
- if (step) /* hack to prevent long inputfiles */
- {
- sprintf(buf, "%s < %s > %s", exe, "input.com", "input.log");
- }
- else
- {
- sprintf(buf, "%s < %s > %s", exe, "input.com", "input.log");
- }
- fprintf(stderr, "Calling '%s'\n", buf);
- if (system(buf) != 0)
- {
- gmx_fatal(FARGS, "Call to '%s' failed\n", buf);
- }
-}
-
-real call_gaussian(const t_QMMMrec* qmmm, t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[])
-{
- /* normal gaussian jobs */
- static int step = 0;
- int i, j;
- real QMener = 0.0;
- rvec * QMgrad, *MMgrad;
- char* exe;
-
- snew(exe, 30);
- sprintf(exe, "%s/%s", qm->gauss_dir, qm->gauss_exe);
- snew(QMgrad, qm->nrQMatoms);
- snew(MMgrad, mm->nrMMatoms);
-
- write_gaussian_input(step, qmmm, qm, mm);
- do_gaussian(step, exe);
- QMener = read_gaussian_output(QMgrad, MMgrad, qm, mm);
- /* put the QMMM forces in the force array and to the fshift
- */
- for (i = 0; i < qm->nrQMatoms; i++)
- {
- for (j = 0; j < DIM; j++)
- {
- f[i][j] = HARTREE_BOHR2MD * QMgrad[i][j];
- fshift[i][j] = HARTREE_BOHR2MD * QMgrad[i][j];
- }
- }
- for (i = 0; i < mm->nrMMatoms; i++)
- {
- for (j = 0; j < DIM; j++)
- {
- f[i + qm->nrQMatoms][j] = HARTREE_BOHR2MD * MMgrad[i][j];
- fshift[i + qm->nrQMatoms][j] = HARTREE_BOHR2MD * MMgrad[i][j];
- }
- }
- QMener = QMener * HARTREE2KJ * AVOGADRO;
- step++;
- free(exe);
- return (QMener);
-
-} /* call_gaussian */
-
-real call_gaussian_SH(const t_QMMMrec* qmmm, t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[])
-{
- /* a gaussian call routine intended for doing diabatic surface
- * "sliding". See the manual for the theoretical background of this
- * TSH method.
- */
- static int step = 0;
- int state, i, j;
- real QMener = 0.0;
- static gmx_bool swapped = FALSE; /* handle for identifying the current PES */
- gmx_bool swap = FALSE; /* the actual swap */
- rvec * QMgrad, *MMgrad;
- char* buf;
- char* exe;
-
- snew(exe, 30);
- sprintf(exe, "%s/%s", qm->gauss_dir, qm->gauss_exe);
- /* hack to do ground state simulations */
- if (!step)
- {
- snew(buf, 20);
- buf = getenv("GMX_QM_GROUND_STATE");
- if (buf)
- {
- sscanf(buf, "%d", &state);
- }
- else
- {
- state = 2;
- }
- if (state == 1)
- {
- swapped = TRUE;
- }
- }
- /* end of hack */
-
-
- /* copy the QMMMrec pointer */
- snew(QMgrad, qm->nrQMatoms);
- snew(MMgrad, mm->nrMMatoms);
- /* at step 0 there should be no SA */
- /* if(!step)
- * qr->bSA=FALSE;*/
- /* temporray set to step + 1, since there is a chk start */
- write_gaussian_SH_input(step, swapped, qmmm, qm, mm);
-
- do_gaussian(step, exe);
- QMener = read_gaussian_SH_output(QMgrad, MMgrad, step, qm, mm);
-
- /* check for a surface hop. Only possible if we were already state
- * averaging.
- */
- if (qm->SAstep > 0)
- {
- if (!swapped)
- {
- swap = ((step != 0) && (hop(step, qm) != 0));
- swapped = swap;
- }
- else /* already on the other surface, so check if we go back */
- {
- swap = ((step != 0) && (hop(step, qm) != 0));
- swapped = !swap; /* so swapped shoud be false again */
- }
- if (swap) /* change surface, so do another call */
- {
- write_gaussian_SH_input(step, swapped, qmmm, qm, mm);
- do_gaussian(step, exe);
- QMener = read_gaussian_SH_output(QMgrad, MMgrad, step, qm, mm);
- }
- }
- /* add the QMMM forces to the gmx force array and fshift
- */
- for (i = 0; i < qm->nrQMatoms; i++)
- {
- for (j = 0; j < DIM; j++)
- {
- f[i][j] = HARTREE_BOHR2MD * QMgrad[i][j];
- fshift[i][j] = HARTREE_BOHR2MD * QMgrad[i][j];
- }
- }
- for (i = 0; i < mm->nrMMatoms; i++)
- {
- for (j = 0; j < DIM; j++)
- {
- f[i + qm->nrQMatoms][j] = HARTREE_BOHR2MD * MMgrad[i][j];
- fshift[i + qm->nrQMatoms][j] = HARTREE_BOHR2MD * MMgrad[i][j];
- }
- }
- QMener = QMener * HARTREE2KJ * AVOGADRO;
- fprintf(stderr, "step %5d, SA = %5d, swap = %5d\n", step, static_cast<int>(qm->SAstep > 0),
- static_cast<int>(swapped));
- step++;
- free(exe);
- return (QMener);
-
-} /* call_gaussian_SH */
-
-#pragma GCC diagnostic pop
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 2018,2019, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-#ifndef GMX_MDLIB_QMGAUSSIAN_H
-#define GMX_MDLIB_QMGAUSSIAN_H
-
-#include "gromacs/math/vectypes.h"
-#include "gromacs/mdlib/qmmm.h"
-
-struct t_forcerec;
-
-/*! \brief
- * Initialize gaussian datastructures.
- *
- * \param[in] qm QM forcerec.
- */
-void init_gaussian(t_QMrec* qm);
-
-/*! \brief
- * Call gaussian to do qm calculation.
- *
- * \param[in] qmmm QMMM part forcerec.
- * \param[in] qm QM part of forcerec.
- * \param[in] mm mm part of forcerec.
- * \param[in] f force vector.
- * \param[in] fshift shift of force vector.
- */
-real call_gaussian(const t_QMMMrec* qmmm, t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[]);
-
-/*! \brief
- * Call gaussian SH(?) to do qm calculation.
- *
- * \param[in] qmmm QMMM part forcerec.
- * \param[in] qm QM part of forcerec.
- * \param[in] mm mm part of forcerec.
- * \param[in] f force vector.
- * \param[in] fshift shift of force vector.
- */
-real call_gaussian_SH(const t_QMMMrec* qmmm, t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[]);
-
-#endif
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2017,2018 by the GROMACS development team.
- * Copyright (c) 2019,2020, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-#include "gmxpre.h"
-
-#include "qm_mopac.h"
-
-#include "config.h"
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-
-#include <cmath>
-
-#include "gromacs/fileio/confio.h"
-#include "gromacs/gmxlib/network.h"
-#include "gromacs/gmxlib/nrnb.h"
-#include "gromacs/math/units.h"
-#include "gromacs/math/vec.h"
-#include "gromacs/mdlib/qmmm.h"
-#include "gromacs/mdtypes/md_enums.h"
-#include "gromacs/utility/fatalerror.h"
-#include "gromacs/utility/smalloc.h"
-
-
-// When not built in a configuration with QMMM support, much of this
-// code is unreachable by design. Tell clang not to warn about it.
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wmissing-noreturn"
-
-#if GMX_QMMM_MOPAC
-/* mopac interface routines */
-void F77_FUNC(domldt, DOMLDT)(int* nrqmat, int labels[], char keywords[]);
-
-void F77_FUNC(domop, DOMOP)(int* nrqmat,
- double qmcrd[],
- int* nrmmat,
- double mmchrg[],
- double mmcrd[],
- double qmgrad[],
- double mmgrad[],
- double* energy,
- double qmcharges[]);
-
-#else /* GMX_QMMM_MOPAC */
-// Stub definitions to make compilation succeed when not configured
-// for MOPAC support. In that case, the module gives a fatal error
-// when the initialization function is called, so there is no need to
-// issue fatal errors here, because that introduces problems with
-// tools suggesting and prohibiting noreturn attributes.
-
-static void F77_FUNC(domldt, DOMLDT)(int* /*unused*/, int /*unused*/[], char /*unused*/[]) {}
-
-static void F77_FUNC(domop, DOMOP)(int* /*unused*/,
- double /*unused*/[],
- int* /*unused*/,
- double /*unused*/[],
- double /*unused*/[],
- double /*unused*/[],
- double /*unused*/[],
- double* /*unused*/,
- double /*unused*/[])
-{
-}
-
-#endif
-
-
-void init_mopac(t_QMrec* qm)
-{
- /* initializes the mopac routines ans sets up the semiempirical
- * computation by calling moldat(). The inline mopac routines can
- * only perform gradient operations. If one would like to optimize a
- * structure or find a transition state at PM3 level, gaussian is
- * used instead.
- */
- char* keywords;
-
- if (!GMX_QMMM_MOPAC)
- {
- gmx_fatal(FARGS,
- "Cannot call MOPAC unless linked against it. Use cmake -DGMX_QMMM_PROGRAM=MOPAC, "
- "and ensure that linking will work correctly.");
- }
-
- snew(keywords, 240);
-
- if (!qm->bSH) /* if rerun then grad should not be done! */
- {
- sprintf(keywords, "PRECISE GEO-OK CHARGE=%d GRAD MMOK ANALYT %s\n", qm->QMcharge,
- eQMmethod_names[qm->QMmethod]);
- }
- else
- {
- sprintf(keywords, "PRECISE GEO-OK CHARGE=%d SINGLET GRAD %s C.I.=(%d,%d) root=2 MECI \n",
- qm->QMcharge, eQMmethod_names[qm->QMmethod], qm->CASorbitals, qm->CASelectrons / 2);
- }
- F77_FUNC(domldt, DOMLDT)(&qm->nrQMatoms, qm->atomicnumberQM, keywords);
- fprintf(stderr, "keywords are: %s\n", keywords);
- free(keywords);
-
-} /* init_mopac */
-
-real call_mopac(t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[])
-{
- /* do the actual QMMM calculation using directly linked mopac subroutines
- */
- double /* always double as the MOPAC routines are always compiled in
- double precission! */
- *qmcrd = nullptr,
- *qmchrg = nullptr, *mmcrd = nullptr, *mmchrg = nullptr, *qmgrad, *mmgrad = nullptr,
- energy = 0;
- int i, j;
- real QMener = 0.0;
- snew(qmcrd, 3 * (qm->nrQMatoms));
- snew(qmgrad, 3 * (qm->nrQMatoms));
- /* copy the data from qr into the arrays that are going to be used
- * in the fortran routines of MOPAC
- */
- for (i = 0; i < qm->nrQMatoms; i++)
- {
- for (j = 0; j < DIM; j++)
- {
- qmcrd[3 * i + j] = static_cast<double>(qm->xQM[i][j]) * 10;
- }
- }
- if (mm->nrMMatoms)
- {
- /* later we will add the point charges here. There are some
- * conceptual problems with semi-empirical QM in combination with
- * point charges that we need to solve first....
- */
- gmx_fatal(FARGS,
- "At present only ONIOM is allowed in combination"
- " with MOPAC QM subroutines\n");
- }
- else
- {
- /* now compute the energy and the gradients.
- */
-
- snew(qmchrg, qm->nrQMatoms);
- F77_FUNC(domop, DOMOP)
- (&qm->nrQMatoms, qmcrd, &mm->nrMMatoms, mmchrg, mmcrd, qmgrad, mmgrad, &energy, qmchrg);
- /* add the gradients to the f[] array, and also to the fshift[].
- * the mopac gradients are in kCal/angstrom.
- */
- for (i = 0; i < qm->nrQMatoms; i++)
- {
- for (j = 0; j < DIM; j++)
- {
- f[i][j] = static_cast<real>(10) * CAL2JOULE * qmgrad[3 * i + j];
- fshift[i][j] = static_cast<real>(10) * CAL2JOULE * qmgrad[3 * i + j];
- }
- }
- QMener = static_cast<real> CAL2JOULE * energy;
- /* do we do something with the mulliken charges?? */
-
- free(qmchrg);
- }
- free(qmgrad);
- free(qmcrd);
- return (QMener);
-}
-
-real call_mopac_SH(t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[])
-{
- /* do the actual SH QMMM calculation using directly linked mopac
- subroutines */
-
- double /* always double as the MOPAC routines are always compiled in
- double precission! */
- *qmcrd = nullptr,
- *qmchrg = nullptr, *mmcrd = nullptr, *mmchrg = nullptr, *qmgrad, *mmgrad = nullptr,
- energy = 0;
- int i, j;
- real QMener = 0.0;
-
- snew(qmcrd, 3 * (qm->nrQMatoms));
- snew(qmgrad, 3 * (qm->nrQMatoms));
- /* copy the data from qr into the arrays that are going to be used
- * in the fortran routines of MOPAC
- */
- for (i = 0; i < qm->nrQMatoms; i++)
- {
- for (j = 0; j < DIM; j++)
- {
- qmcrd[3 * i + j] = static_cast<double>(qm->xQM[i][j]) * 10;
- }
- }
- if (mm->nrMMatoms)
- {
- /* later we will add the point charges here. There are some
- * conceptual problems with semi-empirical QM in combination with
- * point charges that we need to solve first....
- */
- gmx_fatal(FARGS, "At present only ONIOM is allowed in combination with MOPAC\n");
- }
- else
- {
- /* now compute the energy and the gradients.
- */
- snew(qmchrg, qm->nrQMatoms);
-
- F77_FUNC(domop, DOMOP)
- (&qm->nrQMatoms, qmcrd, &mm->nrMMatoms, mmchrg, mmcrd, qmgrad, mmgrad, &energy, qmchrg);
- /* add the gradients to the f[] array, and also to the fshift[].
- * the mopac gradients are in kCal/angstrom.
- */
- for (i = 0; i < qm->nrQMatoms; i++)
- {
- for (j = 0; j < DIM; j++)
- {
- f[i][j] = static_cast<real>(10) * CAL2JOULE * qmgrad[3 * i + j];
- fshift[i][j] = static_cast<real>(10) * CAL2JOULE * qmgrad[3 * i + j];
- }
- }
- QMener = static_cast<real> CAL2JOULE * energy;
- }
- free(qmgrad);
- free(qmcrd);
- return (QMener);
-} /* call_mopac_SH */
-
-#pragma GCC diagnostic pop
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 2018,2019, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-#ifndef GMX_MDLIB_QMMOPAC_H
-#define GMX_MDLIB_QMMOPAC_H
-
-#include "gromacs/mdlib/qmmm.h"
-
-/*! \brief
- * Initialize mopac datastructures.
- *
- * \param[in] qm QM forcerec.
- */
-void init_mopac(t_QMrec* qm);
-
-/*! \brief
- * Run calculation with MOPAC.
- *
- * \param[in] qm QM part of forcerec.
- * \param[in] mm MM part of forcerec.
- * \param[in] f Force vector.
- * \param[in] fshift Force shift vector.
- */
-real call_mopac(t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[]);
-
-/*! \brief
- * Run surface-hopping calculation with MOPAC.
- *
- * \param[in] qm QM part of forcerec.
- * \param[in] mm MM part of forcerec.
- * \param[in] f Force vector.
- * \param[in] fshift Force shift vector.
- */
-real call_mopac_SH(t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[]);
-
-#endif
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2008, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2017,2018 by the GROMACS development team.
- * Copyright (c) 2019,2020, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-#include "gmxpre.h"
-
-#include "qm_orca.h"
-
-#include "config.h"
-
-#include <cmath>
-#include <cstdio>
-#include <cstdlib>
-#include <cstring>
-
-#include "gromacs/fileio/confio.h"
-#include "gromacs/gmxlib/network.h"
-#include "gromacs/gmxlib/nrnb.h"
-#include "gromacs/math/units.h"
-#include "gromacs/math/vec.h"
-#include "gromacs/mdlib/qmmm.h"
-#include "gromacs/mdtypes/md_enums.h"
-#include "gromacs/utility/fatalerror.h"
-#include "gromacs/utility/smalloc.h"
-
-// When not built in a configuration with QMMM support, much of this
-// code is unreachable by design. Tell clang not to warn about it.
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wmissing-noreturn"
-
-/* ORCA interface routines */
-
-void init_orca(t_QMrec* qm)
-{
- char* buf;
- snew(buf, 200);
-
- if (!GMX_QMMM_ORCA)
- {
- gmx_fatal(FARGS,
- "Cannot call ORCA unless linked against it. Use cmake -DGMX_QMMM_PROGRAM=ORCA, "
- "and ensure that linking will work correctly.");
- }
-
- /* ORCA settings on the system */
- buf = getenv("GMX_QM_ORCA_BASENAME");
- if (buf)
- {
- snew(qm->orca_basename, 200);
- sscanf(buf, "%s", qm->orca_basename);
- }
- else
- {
- gmx_fatal(FARGS, "$GMX_QM_ORCA_BASENAME is not set\n");
- }
-
- /* ORCA directory on the system */
- snew(buf, 200);
- buf = getenv("GMX_ORCA_PATH");
-
- if (buf)
- {
- snew(qm->orca_dir, 200);
- sscanf(buf, "%s", qm->orca_dir);
- }
- else
- {
- gmx_fatal(FARGS, "$GMX_ORCA_PATH not set, check manual\n");
- }
-
- fprintf(stderr, "Setting ORCA path to: %s...\n", qm->orca_dir);
- fprintf(stderr, "ORCA initialised...\n\n");
- /* since we append the output to the BASENAME.out file,
- we should delete an existent old out-file here. */
- sprintf(buf, "%s.out", qm->orca_basename);
- remove(buf);
-}
-
-
-static void write_orca_input(const t_QMMMrec* QMMMrec, t_QMrec* qm, t_MMrec* mm)
-{
- int i;
- FILE *pcFile, *addInputFile;
- char *buf, *orcaInput, *addInputFilename, *pcFilename;
-
- /* write the first part of the input-file */
- snew(orcaInput, 200);
- sprintf(orcaInput, "%s.inp", qm->orca_basename);
- FILE* out = fopen(orcaInput, "w");
-
- snew(addInputFilename, 200);
- sprintf(addInputFilename, "%s.ORCAINFO", qm->orca_basename);
- addInputFile = fopen(addInputFilename, "r");
-
- fprintf(out, "#input-file generated by GROMACS\n");
-
- fprintf(out, "!EnGrad TightSCF\n");
-
- /* here we include the insertion of the additional orca-input */
- snew(buf, 200);
- if (addInputFile != nullptr)
- {
- while (!feof(addInputFile))
- {
- if (fgets(buf, 200, addInputFile) != nullptr)
- {
- fputs(buf, out);
- }
- }
- }
- else
- {
- gmx_fatal(FARGS, "No information on the calculation given in %s\n", addInputFilename);
- }
-
- fclose(addInputFile);
-
- /* write charge and multiplicity */
- fprintf(out, "*xyz %2d%2d\n", qm->QMcharge, qm->multiplicity);
-
- /* write the QM coordinates */
- for (i = 0; i < qm->nrQMatoms; i++)
- {
- int atomNr;
- if (qm->atomicnumberQM[i] == 0)
- {
- atomNr = 1;
- }
- else
- {
- atomNr = qm->atomicnumberQM[i];
- }
- fprintf(out, "%3d %10.7f %10.7f %10.7f\n", atomNr, qm->xQM[i][XX] / 0.1,
- qm->xQM[i][YY] / 0.1, qm->xQM[i][ZZ] / 0.1);
- }
- fprintf(out, "*\n");
-
- /* write the MM point charge data */
- if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms)
- {
- /* name of the point charge file */
- snew(pcFilename, 200);
- sprintf(pcFilename, "%s.pc", qm->orca_basename);
- fprintf(out, "%s%s%s\n", "%pointcharges \"", pcFilename, "\"");
- pcFile = fopen(pcFilename, "w");
- fprintf(pcFile, "%d\n", mm->nrMMatoms);
- for (i = 0; i < mm->nrMMatoms; i++)
- {
- fprintf(pcFile, "%8.4f %10.7f %10.7f %10.7f\n", mm->MMcharges[i],
- mm->xMM[i][XX] / 0.1, mm->xMM[i][YY] / 0.1, mm->xMM[i][ZZ] / 0.1);
- }
- fprintf(pcFile, "\n");
- fclose(pcFile);
- }
- fprintf(out, "\n");
-
- fclose(out);
-} /* write_orca_input */
-
-static real read_orca_output(rvec QMgrad[], rvec MMgrad[], const t_QMMMrec* QMMMrec, t_QMrec* qm, t_MMrec* mm)
-{
- int i, j;
- char buf[300], orca_pcgradFilename[300], orca_engradFilename[300];
- real QMener;
- FILE *pcgrad, *engrad;
- int k;
-
- /* the energy and gradients for the QM part are stored in the engrad file
- * and the gradients for the point charges are stored in the pc file.
- */
- sprintf(orca_engradFilename, "%s.engrad", qm->orca_basename);
- engrad = fopen(orca_engradFilename, "r");
- /* we read the energy and the gradient for the qm-atoms from the engrad file
- */
- /* we can skip the first seven lines
- */
- for (j = 0; j < 7; j++)
- {
- if (fgets(buf, 300, engrad) == nullptr)
- {
- gmx_fatal(FARGS, "Unexpected end of ORCA output");
- }
- }
- /* now comes the energy
- */
- if (fgets(buf, 300, engrad) == nullptr)
- {
- gmx_fatal(FARGS, "Unexpected end of ORCA output");
- }
-#if GMX_DOUBLE
- sscanf(buf, "%lf\n", &QMener);
-#else
- sscanf(buf, "%f\n", &QMener);
-#endif
- /* we can skip the next three lines
- */
- for (j = 0; j < 3; j++)
- {
- if (fgets(buf, 300, engrad) == nullptr)
- {
- gmx_fatal(FARGS, "Unexpected end of ORCA output");
- }
- }
- /* next lines contain the gradients of the QM atoms
- * now comes the gradient, one value per line:
- * (atom1 x \n atom1 y \n atom1 z \n atom2 x ...
- */
-
- for (i = 0; i < 3 * qm->nrQMatoms; i++)
- {
- k = i / 3;
- if (fgets(buf, 300, engrad) == nullptr)
- {
- gmx_fatal(FARGS, "Unexpected end of ORCA output");
- }
-#if GMX_DOUBLE
- if (i % 3 == 0)
- {
- sscanf(buf, "%lf\n", &QMgrad[k][XX]);
- }
- else if (i % 3 == 1)
- {
- sscanf(buf, "%lf\n", &QMgrad[k][YY]);
- }
- else if (i % 3 == 2)
- {
- sscanf(buf, "%lf\n", &QMgrad[k][ZZ]);
- }
-#else
- if (i % 3 == 0)
- {
- sscanf(buf, "%f\n", &QMgrad[k][XX]);
- }
- else if (i % 3 == 1)
- {
- sscanf(buf, "%f\n", &QMgrad[k][YY]);
- }
- else if (i % 3 == 2)
- {
- sscanf(buf, "%f\n", &QMgrad[k][ZZ]);
- }
-#endif
- }
- fclose(engrad);
- /* write the MM point charge data
- */
- if (QMMMrec->QMMMscheme != eQMMMschemeoniom && mm->nrMMatoms)
- {
- sprintf(orca_pcgradFilename, "%s.pcgrad", qm->orca_basename);
- pcgrad = fopen(orca_pcgradFilename, "r");
-
- /* we read the gradient for the mm-atoms from the pcgrad file
- */
- /* we can skip the first line
- */
- if (fgets(buf, 300, pcgrad) == nullptr)
- {
- gmx_fatal(FARGS, "Unexpected end of ORCA output");
- }
- for (i = 0; i < mm->nrMMatoms; i++)
- {
- if (fgets(buf, 300, pcgrad) == nullptr)
- {
- gmx_fatal(FARGS, "Unexpected end of ORCA output");
- }
-#if GMX_DOUBLE
- sscanf(buf, "%lf%lf%lf\n", &MMgrad[i][XX], &MMgrad[i][YY], &MMgrad[i][ZZ]);
-#else
- sscanf(buf, "%f%f%f\n", &MMgrad[i][XX], &MMgrad[i][YY], &MMgrad[i][ZZ]);
-#endif
- }
- fclose(pcgrad);
- }
- return (QMener);
-}
-
-static void do_orca(char* orca_dir, char* basename)
-{
-
- /* make the call to the orca binary through system()
- * The location of the binary is set through the
- * environment.
- */
- char buf[100];
- sprintf(buf, "%s/%s %s.inp >> %s.out", orca_dir, "orca", basename, basename);
- fprintf(stderr, "Calling '%s'\n", buf);
- if (system(buf) != 0)
- {
- gmx_fatal(FARGS, "Call to '%s' failed\n", buf);
- }
-}
-
-real call_orca(const t_QMMMrec* qmmm, t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[])
-{
- /* normal orca jobs */
- static int step = 0;
- int i, j;
- real QMener;
- rvec * QMgrad, *MMgrad;
- char* exe;
-
- snew(exe, 30);
- sprintf(exe, "%s", "orca");
- snew(QMgrad, qm->nrQMatoms);
- snew(MMgrad, mm->nrMMatoms);
-
- write_orca_input(qmmm, qm, mm);
- do_orca(qm->orca_dir, qm->orca_basename);
- QMener = read_orca_output(QMgrad, MMgrad, qmmm, qm, mm);
- /* put the QMMM forces in the force array and to the fshift
- */
- for (i = 0; i < qm->nrQMatoms; i++)
- {
- for (j = 0; j < DIM; j++)
- {
- f[i][j] = HARTREE_BOHR2MD * QMgrad[i][j];
- fshift[i][j] = HARTREE_BOHR2MD * QMgrad[i][j];
- }
- }
- for (i = 0; i < mm->nrMMatoms; i++)
- {
- for (j = 0; j < DIM; j++)
- {
- f[i + qm->nrQMatoms][j] = HARTREE_BOHR2MD * MMgrad[i][j];
- fshift[i + qm->nrQMatoms][j] = HARTREE_BOHR2MD * MMgrad[i][j];
- }
- }
- QMener = QMener * HARTREE2KJ * AVOGADRO;
- step++;
- free(exe);
- return (QMener);
-} /* call_orca */
-
-/* end of orca sub routines */
-
-#pragma GCC diagnostic pop
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 2017,2018,2019, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-#include "gromacs/mdlib/qmmm.h"
-
-#ifndef GMX_MDLIB_QM_ORCA_H
-# define GMX_MDLIB_QM_ORCA_H
-
-void init_orca(t_QMrec* qm);
-
-real call_orca(const t_QMMMrec* qmmm, t_QMrec* qm, t_MMrec* mm, rvec f[], rvec fshift[]);
-
-#endif
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-#include "gmxpre.h"
-
-#include "qmmm.h"
-
-#include "config.h"
-
-#include <cmath>
-#include <cstdio>
-#include <cstdlib>
-#include <cstring>
-
-#include <algorithm>
-
-#include "gromacs/domdec/domdec_struct.h"
-#include "gromacs/fileio/confio.h"
-#include "gromacs/gmxlib/network.h"
-#include "gromacs/gmxlib/nrnb.h"
-#include "gromacs/math/functions.h"
-#include "gromacs/math/units.h"
-#include "gromacs/math/vec.h"
-#include "gromacs/mdlib/qm_gamess.h"
-#include "gromacs/mdlib/qm_gaussian.h"
-#include "gromacs/mdlib/qm_mopac.h"
-#include "gromacs/mdlib/qm_orca.h"
-#include "gromacs/mdtypes/commrec.h"
-#include "gromacs/mdtypes/forceoutput.h"
-#include "gromacs/mdtypes/forcerec.h"
-#include "gromacs/mdtypes/inputrec.h"
-#include "gromacs/mdtypes/md_enums.h"
-#include "gromacs/mdtypes/mdatom.h"
-#include "gromacs/mdtypes/nblist.h"
-#include "gromacs/pbcutil/ishift.h"
-#include "gromacs/pbcutil/pbc.h"
-#include "gromacs/topology/mtop_lookup.h"
-#include "gromacs/topology/mtop_util.h"
-#include "gromacs/topology/topology.h"
-#include "gromacs/utility/fatalerror.h"
-#include "gromacs/utility/smalloc.h"
-
-// When not built in a configuration with QMMM support, much of this
-// code is unreachable by design. Tell clang not to warn about it.
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wunreachable-code"
-#pragma GCC diagnostic ignored "-Wmissing-noreturn"
-
-/* this struct and these comparison functions are needed for creating
- * a QMMM input for the QM routines from the QMMM neighbor list.
- */
-
-typedef struct
-{
- int j;
- int shift;
-} t_j_particle;
-
-static bool struct_comp(const t_j_particle& a, const t_j_particle& b)
-{
- return a.j < b.j;
-}
-
-static real call_QMroutine(const t_commrec gmx_unused* cr,
- const t_QMMMrec gmx_unused* qmmm,
- t_QMrec gmx_unused* qm,
- t_MMrec gmx_unused* mm,
- rvec gmx_unused f[],
- rvec gmx_unused fshift[])
-{
- /* makes a call to the requested QM routine (qm->QMmethod)
- * Note that f is actually the gradient, i.e. -f
- */
- /* do a semi-empiprical calculation */
-
- if (qm->QMmethod < eQMmethodRHF && !(mm->nrMMatoms))
- {
- if (GMX_QMMM_MOPAC)
- {
- if (qm->bSH)
- {
- return call_mopac_SH(qm, mm, f, fshift);
- }
- else
- {
- return call_mopac(qm, mm, f, fshift);
- }
- }
- else
- {
- gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac.");
- }
- }
- else
- {
- /* do an ab-initio calculation */
- if (qm->bSH && qm->QMmethod == eQMmethodCASSCF)
- {
- if (GMX_QMMM_GAUSSIAN)
- {
- return call_gaussian_SH(qmmm, qm, mm, f, fshift);
- }
- else
- {
- gmx_fatal(FARGS, "Ab-initio Surface-hopping only supported with Gaussian.");
- }
- }
- else
- {
- if (GMX_QMMM_GAMESS)
- {
- return call_gamess(qm, mm, f, fshift);
- }
- else if (GMX_QMMM_GAUSSIAN)
- {
- return call_gaussian(qmmm, qm, mm, f, fshift);
- }
- else if (GMX_QMMM_ORCA)
- {
- return call_orca(qmmm, qm, mm, f, fshift);
- }
- else
- {
- gmx_fatal(FARGS,
- "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
- }
- }
- }
-}
-
-static void init_QMroutine(const t_commrec gmx_unused* cr, t_QMrec gmx_unused* qm, t_MMrec gmx_unused* mm)
-{
- /* makes a call to the requested QM routine (qm->QMmethod)
- */
- if (qm->QMmethod < eQMmethodRHF)
- {
- if (GMX_QMMM_MOPAC)
- {
- /* do a semi-empiprical calculation */
- init_mopac(qm);
- }
- else
- {
- gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac.");
- }
- }
- else
- {
- /* do an ab-initio calculation */
- if (GMX_QMMM_GAMESS)
- {
- init_gamess(cr, qm, mm);
- }
- else if (GMX_QMMM_GAUSSIAN)
- {
- init_gaussian(qm);
- }
- else if (GMX_QMMM_ORCA)
- {
- init_orca(qm);
- }
- else
- {
- gmx_fatal(FARGS, "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
- }
- }
-} /* init_QMroutine */
-
-static void update_QMMM_coord(const rvec* x, const t_forcerec* fr, t_QMrec* qm, t_MMrec* mm)
-{
- /* shifts the QM and MM particles into the central box and stores
- * these shifted coordinates in the coordinate arrays of the
- * QMMMrec. These coordinates are passed on the QM subroutines.
- */
- int i;
-
- /* shift the QM atoms into the central box
- */
- for (i = 0; i < qm->nrQMatoms; i++)
- {
- rvec_sub(x[qm->indexQM[i]], fr->shift_vec[qm->shiftQM[i]], qm->xQM[i]);
- }
- /* also shift the MM atoms into the central box, if any
- */
- for (i = 0; i < mm->nrMMatoms; i++)
- {
- rvec_sub(x[mm->indexMM[i]], fr->shift_vec[mm->shiftMM[i]], mm->xMM[i]);
- }
-} /* update_QMMM_coord */
-
-/* end of QMMM subroutines */
-
-/* QMMM core routines */
-
-static t_QMrec* mk_QMrec()
-{
- t_QMrec* qm;
- snew(qm, 1);
- return qm;
-} /* mk_QMrec */
-
-static t_MMrec* mk_MMrec()
-{
- t_MMrec* mm;
- snew(mm, 1);
- return mm;
-} /* mk_MMrec */
-
-static void init_QMrec(int grpnr, t_QMrec* qm, int nr, const int* atomarray, const gmx_mtop_t* mtop, const t_inputrec* ir)
-{
- /* fills the t_QMrec struct of QM group grpnr
- */
-
- qm->nrQMatoms = nr;
- snew(qm->xQM, nr);
- snew(qm->indexQM, nr);
- snew(qm->shiftQM, nr); /* the shifts */
- for (int i = 0; i < nr; i++)
- {
- qm->indexQM[i] = atomarray[i];
- }
-
- snew(qm->atomicnumberQM, nr);
- int molb = 0;
- for (int i = 0; i < qm->nrQMatoms; i++)
- {
- const t_atom& atom = mtopGetAtomParameters(mtop, qm->indexQM[i], &molb);
- qm->nelectrons += mtop->atomtypes.atomnumber[atom.type];
- qm->atomicnumberQM[i] = mtop->atomtypes.atomnumber[atom.type];
- }
-
- qm->QMcharge = ir->opts.QMcharge[grpnr];
- qm->multiplicity = ir->opts.QMmult[grpnr];
- qm->nelectrons -= ir->opts.QMcharge[grpnr];
-
- qm->QMmethod = ir->opts.QMmethod[grpnr];
- qm->QMbasis = ir->opts.QMbasis[grpnr];
- /* trajectory surface hopping setup (Gaussian only) */
- qm->bSH = ir->opts.bSH[grpnr];
- qm->CASorbitals = ir->opts.CASorbitals[grpnr];
- qm->CASelectrons = ir->opts.CASelectrons[grpnr];
- qm->SAsteps = ir->opts.SAsteps[grpnr];
- qm->SAon = ir->opts.SAon[grpnr];
- qm->SAoff = ir->opts.SAoff[grpnr];
- /* hack to prevent gaussian from reinitializing all the time */
- qm->nQMcpus = 0; /* number of CPU's to be used by g01, is set
- * upon initializing gaussian
- * (init_gaussian()
- */
- /* print the current layer to allow users to check their input */
- fprintf(stderr, "Layer %d\nnr of QM atoms %d\n", grpnr, nr);
- fprintf(stderr, "QMlevel: %s/%s\n\n", eQMmethod_names[qm->QMmethod], eQMbasis_names[qm->QMbasis]);
-} /* init_QMrec */
-
-static t_QMrec* copy_QMrec(t_QMrec* qm)
-{
- /* copies the contents of qm into a new t_QMrec struct */
- t_QMrec* qmcopy;
- int i;
-
- qmcopy = mk_QMrec();
- qmcopy->nrQMatoms = qm->nrQMatoms;
- snew(qmcopy->xQM, qmcopy->nrQMatoms);
- snew(qmcopy->indexQM, qmcopy->nrQMatoms);
- snew(qmcopy->atomicnumberQM, qm->nrQMatoms);
- snew(qmcopy->shiftQM, qmcopy->nrQMatoms); /* the shifts */
- for (i = 0; i < qmcopy->nrQMatoms; i++)
- {
- qmcopy->shiftQM[i] = qm->shiftQM[i];
- qmcopy->indexQM[i] = qm->indexQM[i];
- qmcopy->atomicnumberQM[i] = qm->atomicnumberQM[i];
- }
- qmcopy->nelectrons = qm->nelectrons;
- qmcopy->multiplicity = qm->multiplicity;
- qmcopy->QMcharge = qm->QMcharge;
- qmcopy->nelectrons = qm->nelectrons;
- qmcopy->QMmethod = qm->QMmethod;
- qmcopy->QMbasis = qm->QMbasis;
- /* trajectory surface hopping setup (Gaussian only) */
- qmcopy->bSH = qm->bSH;
- qmcopy->CASorbitals = qm->CASorbitals;
- qmcopy->CASelectrons = qm->CASelectrons;
- qmcopy->SAsteps = qm->SAsteps;
- qmcopy->SAon = qm->SAon;
- qmcopy->SAoff = qm->SAoff;
-
- /* Gaussian init. variables */
- qmcopy->nQMcpus = qm->nQMcpus;
- for (i = 0; i < DIM; i++)
- {
- qmcopy->SHbasis[i] = qm->SHbasis[i];
- }
- qmcopy->QMmem = qm->QMmem;
- qmcopy->accuracy = qm->accuracy;
- qmcopy->cpmcscf = qm->cpmcscf;
- qmcopy->SAstep = qm->SAstep;
-
- return (qmcopy);
-
-} /*copy_QMrec */
-
-#if GMX_QMMM
-
-t_QMMMrec* mk_QMMMrec()
-{
- t_QMMMrec* qr;
-
- snew(qr, 1);
-
- return qr;
-
-} /* mk_QMMMrec */
-
-#else /* GMX_QMMM */
-
-t_QMMMrec* mk_QMMMrec()
-{
- gmx_incons("Compiled without QMMM");
-} /* mk_QMMMrec */
-#endif
-
-void init_QMMMrec(const t_commrec* cr, const gmx_mtop_t* mtop, const t_inputrec* ir, const t_forcerec* fr)
-{
- /* we put the atomsnumbers of atoms that belong to the QMMM group in
- * an array that will be copied later to QMMMrec->indexQM[..]. Also
- * it will be used to create an QMMMrec->bQMMM index array that
- * simply contains true/false for QM and MM (the other) atoms.
- */
-
- t_QMMMrec* qr;
- t_MMrec* mm;
-
- if (!GMX_QMMM)
- {
- gmx_incons("Compiled without QMMM");
- }
-
- if (ir->cutoff_scheme != ecutsGROUP)
- {
- gmx_fatal(FARGS, "QMMM is currently only supported with cutoff-scheme=group");
- }
- if (!EI_DYNAMICS(ir->eI))
- {
- gmx_fatal(FARGS, "QMMM is only supported with dynamics");
- }
-
- /* issue a fatal if the user wants to run with more than one node */
- if (PAR(cr))
- {
- gmx_fatal(FARGS, "QM/MM does not work in parallel, use a single rank instead\n");
- }
-
- /* Make a local copy of the QMMMrec */
- qr = fr->qr;
-
- /* bQMMM[..] is an array containing TRUE/FALSE for atoms that are
- * QM/not QM. We first set all elemenst at false. Afterwards we use
- * the qm_arr (=MMrec->indexQM) to changes the elements
- * corresponding to the QM atoms at TRUE. */
-
- qr->QMMMscheme = ir->QMMMscheme;
-
- /* we take the possibility into account that a user has
- * defined more than one QM group:
- */
- /* an ugly work-around in case there is only one group In this case
- * the whole system is treated as QM. Otherwise the second group is
- * always the rest of the total system and is treated as MM.
- */
-
- /* small problem if there is only QM.... so no MM */
-
- int numQmmmGroups = ir->opts.ngQM;
-
- if (qr->QMMMscheme == eQMMMschemeoniom)
- {
- qr->nrQMlayers = numQmmmGroups;
- }
- else
- {
- qr->nrQMlayers = 1;
- }
-
- /* there are numQmmmGroups groups of QM atoms. In case of multiple QM groups
- * I assume that the users wants to do ONIOM. However, maybe it
- * should also be possible to define more than one QM subsystem with
- * independent neighbourlists. I have to think about
- * that.. 11-11-2003
- */
- std::vector<int> qmmmAtoms;
- snew(qr->qm, numQmmmGroups);
- for (int i = 0; i < numQmmmGroups; i++)
- {
- /* new layer */
- if (qr->QMMMscheme == eQMMMschemeoniom)
- {
- /* add the atoms to the bQMMM array
- */
-
- /* I assume that users specify the QM groups from small to
- * big(ger) in the mdp file
- */
- qr->qm[i] = mk_QMrec();
- /* store QM atoms in this layer in the QMrec and initialise layer
- */
- init_QMrec(i, qr->qm[i], qmmmAtoms.size(), qmmmAtoms.data(), mtop, ir);
- }
- }
- if (qr->QMMMscheme != eQMMMschemeoniom)
- {
-
- /* standard QMMM, all layers are merged together so there is one QM
- * subsystem and one MM subsystem.
- * Also we set the charges to zero in mtop to prevent the innerloops
- * from doubly counting the electostatic QM MM interaction
- * TODO: Consider doing this in grompp instead.
- */
-
- qr->qm[0] = mk_QMrec();
- /* store QM atoms in the QMrec and initialise
- */
- init_QMrec(0, qr->qm[0], qmmmAtoms.size(), qmmmAtoms.data(), mtop, ir);
-
- /* MM rec creation */
- mm = mk_MMrec();
- mm->scalefactor = ir->scalefactor;
- mm->nrMMatoms = (mtop->natoms) - (qr->qm[0]->nrQMatoms); /* rest of the atoms */
- qr->mm = mm;
- }
- else /* ONIOM */
- { /* MM rec creation */
- mm = mk_MMrec();
- mm->scalefactor = ir->scalefactor;
- mm->nrMMatoms = 0;
- qr->mm = mm;
- }
-
- /* these variables get updated in the update QMMMrec */
-
- if (qr->nrQMlayers == 1)
- {
- /* with only one layer there is only one initialisation
- * needed. Multilayer is a bit more complicated as it requires
- * re-initialisation at every step of the simulation. This is due
- * to the use of COMMON blocks in the fortran QM subroutines.
- */
- if (qr->qm[0]->QMmethod < eQMmethodRHF)
- {
- if (GMX_QMMM_MOPAC)
- {
- /* semi-empiprical 1-layer ONIOM calculation requested (mopac93) */
- init_mopac(qr->qm[0]);
- }
- else
- {
- gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac.");
- }
- }
- else
- {
- /* ab initio calculation requested (gamess/gaussian/ORCA) */
- if (GMX_QMMM_GAMESS)
- {
- init_gamess(cr, qr->qm[0], qr->mm);
- }
- else if (GMX_QMMM_GAUSSIAN)
- {
- init_gaussian(qr->qm[0]);
- }
- else if (GMX_QMMM_ORCA)
- {
- init_orca(qr->qm[0]);
- }
- else
- {
- gmx_fatal(FARGS,
- "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
- }
- }
- }
-} /* init_QMMMrec */
-
-void update_QMMMrec(const t_commrec* cr, const t_forcerec* fr, const rvec* x, const t_mdatoms* md, const matrix box)
-{
- /* updates the coordinates of both QM atoms and MM atoms and stores
- * them in the QMMMrec.
- *
- * NOTE: is NOT yet working if there are no PBC. Also in ns.c, simple
- * ns needs to be fixed!
- */
- int mm_max = 0, mm_nr = 0, mm_nr_new, i, j, is, k, shift;
- t_j_particle *mm_j_particles = nullptr, *qm_i_particles = nullptr;
- t_QMMMrec* qr;
- t_nblist* QMMMlist;
- rvec dx;
- ivec crd;
- t_QMrec* qm;
- t_MMrec* mm;
- t_pbc pbc;
- int* parallelMMarray = nullptr;
-
- if (!GMX_QMMM)
- {
- gmx_incons("Compiled without QMMM");
- }
-
- /* every cpu has this array. On every processor we fill this array
- * with 1's and 0's. 1's indicate the atoms is a QM atom on the
- * current cpu in a later stage these arrays are all summed. indexes
- * > 0 indicate the atom is a QM atom. Every node therefore knows
- * whcih atoms are part of the QM subsystem.
- */
- /* copy some pointers */
- qr = fr->qr;
- mm = qr->mm;
- QMMMlist = fr->QMMMlist;
-
- /* init_pbc(box); needs to be called first, see pbc.h */
- ivec null_ivec;
- clear_ivec(null_ivec);
- set_pbc_dd(&pbc, fr->pbcType, DOMAINDECOMP(cr) ? cr->dd->numCells : null_ivec, FALSE, box);
- /* only in standard (normal) QMMM we need the neighbouring MM
- * particles to provide a electric field of point charges for the QM
- * atoms.
- */
- if (qr->QMMMscheme == eQMMMschemenormal) /* also implies 1 QM-layer */
- {
- /* we NOW create/update a number of QMMMrec entries:
- *
- * 1) the shiftQM, containing the shifts of the QM atoms
- *
- * 2) the indexMM array, containing the index of the MM atoms
- *
- * 3) the shiftMM, containing the shifts of the MM atoms
- *
- * 4) the shifted coordinates of the MM atoms
- *
- * the shifts are used for computing virial of the QM/MM particles.
- */
- qm = qr->qm[0]; /* in case of normal QMMM, there is only one group */
- snew(qm_i_particles, QMMMlist->nri);
- if (QMMMlist->nri)
- {
- qm_i_particles[0].shift = XYZ2IS(0, 0, 0);
- for (i = 0; i < QMMMlist->nri; i++)
- {
- qm_i_particles[i].j = QMMMlist->iinr[i];
-
- if (i)
- {
- qm_i_particles[i].shift =
- pbc_dx_aiuc(&pbc, x[QMMMlist->iinr[0]], x[QMMMlist->iinr[i]], dx);
- }
- /* However, since nri >= nrQMatoms, we do a quicksort, and throw
- * out double, triple, etc. entries later, as we do for the MM
- * list too.
- */
-
- /* compute the shift for the MM j-particles with respect to
- * the QM i-particle and store them.
- */
-
- crd[0] = IS2X(QMMMlist->shift[i]) + IS2X(qm_i_particles[i].shift);
- crd[1] = IS2Y(QMMMlist->shift[i]) + IS2Y(qm_i_particles[i].shift);
- crd[2] = IS2Z(QMMMlist->shift[i]) + IS2Z(qm_i_particles[i].shift);
- is = XYZ2IS(crd[0], crd[1], crd[2]);
- for (j = QMMMlist->jindex[i]; j < QMMMlist->jindex[i + 1]; j++)
- {
- if (mm_nr >= mm_max)
- {
- mm_max += 1000;
- srenew(mm_j_particles, mm_max);
- }
-
- mm_j_particles[mm_nr].j = QMMMlist->jjnr[j];
- mm_j_particles[mm_nr].shift = is;
- mm_nr++;
- }
- }
-
- /* quicksort QM and MM shift arrays and throw away multiple entries */
-
-
- std::sort(qm_i_particles, qm_i_particles + QMMMlist->nri, struct_comp);
- /* The mm_j_particles argument to qsort is not allowed to be nullptr */
- if (mm_nr > 0)
- {
- std::sort(mm_j_particles, mm_j_particles + mm_nr, struct_comp);
- }
- /* remove multiples in the QM shift array, since in init_QMMM() we
- * went through the atom numbers from 0 to md.nr, the order sorted
- * here matches the one of QMindex already.
- */
- j = 0;
- for (i = 0; i < QMMMlist->nri; i++)
- {
- if (i == 0 || qm_i_particles[i].j != qm_i_particles[i - 1].j)
- {
- qm_i_particles[j++] = qm_i_particles[i];
- }
- }
- mm_nr_new = 0;
- /* Remove double entries for the MM array.
- * Also remove mm atoms that have no charges!
- * actually this is already done in the ns.c
- */
- for (i = 0; i < mm_nr; i++)
- {
- if ((i == 0 || mm_j_particles[i].j != mm_j_particles[i - 1].j)
- && !md->bQM[mm_j_particles[i].j]
- && ((md->chargeA[mm_j_particles[i].j] != 0.0_real)
- || (md->chargeB && (md->chargeB[mm_j_particles[i].j] != 0.0_real))))
- {
- mm_j_particles[mm_nr_new++] = mm_j_particles[i];
- }
- }
- mm_nr = mm_nr_new;
- /* store the data retrieved above into the QMMMrec
- */
- k = 0;
- /* Keep the compiler happy,
- * shift will always be set in the loop for i=0
- */
- shift = 0;
- for (i = 0; i < qm->nrQMatoms; i++)
- {
- /* not all qm particles might have appeared as i
- * particles. They might have been part of the same charge
- * group for instance.
- */
- if (qm->indexQM[i] == qm_i_particles[k].j)
- {
- shift = qm_i_particles[k++].shift;
- }
- /* use previous shift, assuming they belong the same charge
- * group anyway,
- */
-
- qm->shiftQM[i] = shift;
- }
- }
- /* parallel excecution */
- if (PAR(cr))
- {
- snew(parallelMMarray, 2 * (md->nr));
- /* only MM particles have a 1 at their atomnumber. The second part
- * of the array contains the shifts. Thus:
- * p[i]=1/0 depending on wether atomnumber i is a MM particle in the QM
- * step or not. p[i+md->nr] is the shift of atomnumber i.
- */
- for (i = 0; i < 2 * (md->nr); i++)
- {
- parallelMMarray[i] = 0;
- }
-
- for (i = 0; i < mm_nr; i++)
- {
- parallelMMarray[mm_j_particles[i].j] = 1;
- parallelMMarray[mm_j_particles[i].j + (md->nr)] = mm_j_particles[i].shift;
- }
- gmx_sumi(md->nr, parallelMMarray, cr);
- mm_nr = 0;
-
- mm_max = 0;
- for (i = 0; i < md->nr; i++)
- {
- if (parallelMMarray[i])
- {
- if (mm_nr >= mm_max)
- {
- mm_max += 1000;
- srenew(mm->indexMM, mm_max);
- srenew(mm->shiftMM, mm_max);
- }
- mm->indexMM[mm_nr] = i;
- mm->shiftMM[mm_nr++] = parallelMMarray[i + md->nr] / parallelMMarray[i];
- }
- }
- mm->nrMMatoms = mm_nr;
- free(parallelMMarray);
- }
- /* serial execution */
- else
- {
- mm->nrMMatoms = mm_nr;
- srenew(mm->shiftMM, mm_nr);
- srenew(mm->indexMM, mm_nr);
- for (i = 0; i < mm_nr; i++)
- {
- mm->indexMM[i] = mm_j_particles[i].j;
- mm->shiftMM[i] = mm_j_particles[i].shift;
- }
- }
- /* (re) allocate memory for the MM coordiate array. The QM
- * coordinate array was already allocated in init_QMMM, and is
- * only (re)filled in the update_QMMM_coordinates routine
- */
- srenew(mm->xMM, mm->nrMMatoms);
- /* now we (re) fill the array that contains the MM charges with
- * the forcefield charges. If requested, these charges will be
- * scaled by a factor
- */
- srenew(mm->MMcharges, mm->nrMMatoms);
- for (i = 0; i < mm->nrMMatoms; i++) /* no free energy yet */
- {
- mm->MMcharges[i] = md->chargeA[mm->indexMM[i]] * mm->scalefactor;
- }
- /* the next routine fills the coordinate fields in the QMMM rec of
- * both the qunatum atoms and the MM atoms, using the shifts
- * calculated above.
- */
-
- update_QMMM_coord(x, fr, qr->qm[0], qr->mm);
- free(qm_i_particles);
- free(mm_j_particles);
- }
- else /* ONIOM */ /* ????? */
- {
- mm->nrMMatoms = 0;
- /* do for each layer */
- for (j = 0; j < qr->nrQMlayers; j++)
- {
- qm = qr->qm[j];
- qm->shiftQM[0] = XYZ2IS(0, 0, 0);
- for (i = 1; i < qm->nrQMatoms; i++)
- {
- qm->shiftQM[i] = pbc_dx_aiuc(&pbc, x[qm->indexQM[0]], x[qm->indexQM[i]], dx);
- }
- update_QMMM_coord(x, fr, qm, mm);
- }
- }
-} /* update_QMMM_rec */
-
-real calculate_QMMM(const t_commrec* cr, gmx::ForceWithShiftForces* forceWithShiftForces, const t_QMMMrec* qr)
-{
- real QMener = 0.0;
- /* a selection for the QM package depending on which is requested
- * (Gaussian, GAMESS-UK, MOPAC or ORCA) needs to be implemented here. Now
- * it works through defines.... Not so nice yet
- */
- t_QMrec *qm, *qm2;
- t_MMrec* mm = nullptr;
- rvec * forces = nullptr, *fshift = nullptr, *forces2 = nullptr,
- *fshift2 = nullptr; /* needed for multilayer ONIOM */
- int i, j, k;
-
- if (!GMX_QMMM)
- {
- gmx_incons("Compiled without QMMM");
- }
-
- /* make a local copy the QMMMrec pointer
- */
- mm = qr->mm;
-
- /* now different procedures are carried out for one layer ONION and
- * normal QMMM on one hand and multilayer oniom on the other
- */
- gmx::ArrayRef<gmx::RVec> fMM = forceWithShiftForces->force();
- gmx::ArrayRef<gmx::RVec> fshiftMM = forceWithShiftForces->shiftForces();
- if (qr->QMMMscheme == eQMMMschemenormal || qr->nrQMlayers == 1)
- {
- qm = qr->qm[0];
- snew(forces, (qm->nrQMatoms + mm->nrMMatoms));
- snew(fshift, (qm->nrQMatoms + mm->nrMMatoms));
- QMener = call_QMroutine(cr, qr, qm, mm, forces, fshift);
- for (i = 0; i < qm->nrQMatoms; i++)
- {
- for (j = 0; j < DIM; j++)
- {
- fMM[qm->indexQM[i]][j] -= forces[i][j];
- fshiftMM[qm->shiftQM[i]][j] += fshift[i][j];
- }
- }
- for (i = 0; i < mm->nrMMatoms; i++)
- {
- for (j = 0; j < DIM; j++)
- {
- fMM[mm->indexMM[i]][j] -= forces[qm->nrQMatoms + i][j];
- fshiftMM[mm->shiftMM[i]][j] += fshift[qm->nrQMatoms + i][j];
- }
- }
- free(forces);
- free(fshift);
- }
- else /* Multi-layer ONIOM */
- {
- for (i = 0; i < qr->nrQMlayers - 1; i++) /* last layer is special */
- {
- qm = qr->qm[i];
- qm2 = copy_QMrec(qr->qm[i + 1]);
-
- qm2->nrQMatoms = qm->nrQMatoms;
-
- for (j = 0; j < qm2->nrQMatoms; j++)
- {
- for (k = 0; k < DIM; k++)
- {
- qm2->xQM[j][k] = qm->xQM[j][k];
- }
- qm2->indexQM[j] = qm->indexQM[j];
- qm2->atomicnumberQM[j] = qm->atomicnumberQM[j];
- qm2->shiftQM[j] = qm->shiftQM[j];
- }
-
- qm2->QMcharge = qm->QMcharge;
- /* this layer at the higher level of theory */
- srenew(forces, qm->nrQMatoms);
- srenew(fshift, qm->nrQMatoms);
- /* we need to re-initialize the QMroutine every step... */
- init_QMroutine(cr, qm, mm);
- QMener += call_QMroutine(cr, qr, qm, mm, forces, fshift);
-
- /* this layer at the lower level of theory */
- srenew(forces2, qm->nrQMatoms);
- srenew(fshift2, qm->nrQMatoms);
- init_QMroutine(cr, qm2, mm);
- QMener -= call_QMroutine(cr, qr, qm2, mm, forces2, fshift2);
- /* E = E1high-E1low The next layer includes the current layer at
- * the lower level of theory, which provides + E2low
- * this is similar for gradients
- */
- for (i = 0; i < qm->nrQMatoms; i++)
- {
- for (j = 0; j < DIM; j++)
- {
- fMM[qm->indexQM[i]][j] -= (forces[i][j] - forces2[i][j]);
- fshiftMM[qm->shiftQM[i]][j] += (fshift[i][j] - fshift2[i][j]);
- }
- }
- free(qm2);
- }
- /* now the last layer still needs to be done: */
- qm = qr->qm[qr->nrQMlayers - 1]; /* C counts from 0 */
- init_QMroutine(cr, qm, mm);
- srenew(forces, qm->nrQMatoms);
- srenew(fshift, qm->nrQMatoms);
- QMener += call_QMroutine(cr, qr, qm, mm, forces, fshift);
- for (i = 0; i < qm->nrQMatoms; i++)
- {
- for (j = 0; j < DIM; j++)
- {
- fMM[qm->indexQM[i]][j] -= forces[i][j];
- fshiftMM[qm->shiftQM[i]][j] += fshift[i][j];
- }
- }
- free(forces);
- free(fshift);
- free(forces2);
- free(fshift2);
- }
- return (QMener);
-} /* calculate_QMMM */
-
-#pragma GCC diagnostic pop
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2008, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2017,2018 by the GROMACS development team.
- * Copyright (c) 2019,2020, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-#ifndef GMX_MDLIB_QMMM_H
-#define GMX_MDLIB_QMMM_H
-
-#include "config.h"
-
-#include <vector>
-
-#include "gromacs/math/vectypes.h"
-#include "gromacs/mdlib/tgroup.h"
-#include "gromacs/utility/arrayref.h"
-
-#define GMX_QMMM (GMX_QMMM_MOPAC || GMX_QMMM_GAMESS || GMX_QMMM_GAUSSIAN || GMX_QMMM_ORCA)
-
-struct gmx_localtop_t;
-struct gmx_mtop_t;
-struct t_commrec;
-struct t_forcerec;
-struct t_inputrec;
-struct t_mdatoms;
-struct t_QMMMrec;
-
-namespace gmx
-{
-class ForceWithShiftForces;
-}
-
-typedef struct
-{
- int nrQMatoms; /* total nr of QM atoms */
- rvec* xQM; /* shifted to center of box */
- int* indexQM; /* atom i = atom indexQM[i] in mdrun */
- int* atomicnumberQM; /* atomic numbers of QM atoms */
- real* QMcharges; /* atomic charges of QM atoms(ONIOM) */
- int* shiftQM;
- int QMcharge; /* charge of the QM system */
- int multiplicity; /* multipicity (no of unpaired eln) */
- int QMmethod; /* see enums.h for all methods */
- int QMbasis; /* see enums.h for all bases */
- int nelectrons; /* total number of elecs in QM region*/
- /* Gaussian specific stuff */
- int nQMcpus; /* no. of CPUs used for the QM calc. */
- int QMmem; /* memory for the gaussian calc. */
- int accuracy; /* convergence criterium (E(-x)) */
- gmx_bool cpmcscf; /* using cpmcscf(l1003)*/
- char* gauss_dir;
- char* gauss_exe;
- char* devel_dir;
- char* orca_basename; /* basename for I/O with orca */
- char* orca_dir; /* directory for ORCA */
- /* Surface hopping stuff */
- gmx_bool bSH; /* surface hopping (diabatic only) */
- real SAon; /* at which energy gap the SA starts */
- real SAoff; /* at which energy gap the SA stops */
- int SAsteps; /* stepwise switchinng on the SA */
- int SAstep; /* current state of SA */
- int CIdim;
- real* CIvec1;
- real* CIvec2;
- real* CIvec1old;
- real* CIvec2old;
- ivec SHbasis;
- int CASelectrons;
- int CASorbitals;
-} t_QMrec;
-
-typedef struct
-{
- int nrMMatoms; /* nr of MM atoms, updated every step*/
- rvec* xMM; /* shifted to center of box */
- int* indexMM; /* atom i = atom indexMM[I] in mdrun */
- real* MMcharges; /* MM point charges in std QMMM calc.*/
- int* shiftMM;
- int* MMatomtype; /* only important for semi-emp. */
- real scalefactor;
-} t_MMrec;
-
-
-typedef struct t_QMMMrec
-{
- int QMMMscheme; /* ONIOM (multi-layer) or normal */
- int nrQMlayers; /* number of QM layers (total layers +1 (MM)) */
- t_QMrec** qm; /* atoms and run params for each QM group */
- t_MMrec* mm; /* there can only be one MM subsystem ! */
-} t_QMMMrec;
-
-void atomic_number(int nr, char*** atomtype, int* nucnum);
-
-t_QMMMrec* mk_QMMMrec();
-/* allocates memory for QMMMrec */
-
-void init_QMMMrec(const t_commrec* cr, const gmx_mtop_t* mtop, const t_inputrec* ir, const t_forcerec* fr);
-
-/* init_QMMMrec initializes the QMMM record. From
- * topology->atoms.atomname and topology->atoms.atomtype the atom
- * names and types are read; from inputrec->QMcharge
- * resp. inputrec->QMmult the nelecs and multiplicity are determined
- * and md->cQMMM gives numbers of the MM and QM atoms
- */
-void update_QMMMrec(const t_commrec* cr, const t_forcerec* fr, const rvec* x, const t_mdatoms* md, const matrix box);
-
-/* update_QMMMrec fills the MM stuff in QMMMrec. The MM atoms are
- * taken froom the neighbourlists of the QM atoms. In a QMMM run this
- * routine should be called at every step, since it updates the MM
- * elements of the t_QMMMrec struct.
- */
-real calculate_QMMM(const t_commrec* cr, gmx::ForceWithShiftForces* forceWithShiftForces, const t_QMMMrec* qmmm);
-
-/* QMMM computes the QM forces. This routine makes either function
- * calls to gmx QM routines (derived from MOPAC7 (semi-emp.) and MPQC
- * (ab initio)) or generates input files for an external QM package
- * (listed in QMMMrec.QMpackage). The binary of the QM package is
- * called by system().
- */
-
-#endif
#include "gromacs/mdlib/force_flags.h"
#include "gromacs/mdlib/forcerec.h"
#include "gromacs/mdlib/gmx_omp_nthreads.h"
-#include "gromacs/mdlib/qmmm.h"
#include "gromacs/mdlib/update.h"
#include "gromacs/mdlib/vsite.h"
#include "gromacs/mdtypes/commrec.h"
}
}
- /* update QMMMrec, if necessary */
- if (fr->bQMMM)
- {
- update_QMMMrec(cr, fr, as_rvec_array(x.unpaddedArrayRef().data()), mdatoms, box);
- }
-
// TODO Force flags should include haveFreeEnergyWork for this domain
if (ddUsesGpuDirectCommunication && (domainWork.haveCpuBondedWork || domainWork.haveFreeEnergyWork))
{
#include "gromacs/mdlib/md_support.h"
#include "gromacs/mdlib/mdatoms.h"
#include "gromacs/mdlib/membed.h"
-#include "gromacs/mdlib/qmmm.h"
#include "gromacs/mdlib/sighandler.h"
#include "gromacs/mdlib/stophandler.h"
+#include "gromacs/mdlib/tgroup.h"
#include "gromacs/mdlib/updategroups.h"
#include "gromacs/mdlib/vsite.h"
#include "gromacs/mdrun/mdmodules.h"
*/
int n_tpi = 0;
- /* QMMM stuff */
- gmx_bool bQMMM = FALSE;
- struct t_QMMMrec* qr = nullptr;
-
- /* QM-MM neighborlists */
- struct t_nblist* QMMMlist = nullptr;
-
/* Limit for printing large forces, negative is don't print */
real print_force = 0;
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff