2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
51 #include "gromacs/domdec/domdec_struct.h"
52 #include "gromacs/fileio/confio.h"
53 #include "gromacs/gmxlib/network.h"
54 #include "gromacs/gmxlib/nrnb.h"
55 #include "gromacs/math/functions.h"
56 #include "gromacs/math/units.h"
57 #include "gromacs/math/vec.h"
58 #include "gromacs/mdlib/qm_gamess.h"
59 #include "gromacs/mdlib/qm_gaussian.h"
60 #include "gromacs/mdlib/qm_mopac.h"
61 #include "gromacs/mdlib/qm_orca.h"
62 #include "gromacs/mdtypes/commrec.h"
63 #include "gromacs/mdtypes/forceoutput.h"
64 #include "gromacs/mdtypes/forcerec.h"
65 #include "gromacs/mdtypes/inputrec.h"
66 #include "gromacs/mdtypes/md_enums.h"
67 #include "gromacs/mdtypes/mdatom.h"
68 #include "gromacs/mdtypes/nblist.h"
69 #include "gromacs/pbcutil/ishift.h"
70 #include "gromacs/pbcutil/pbc.h"
71 #include "gromacs/topology/mtop_lookup.h"
72 #include "gromacs/topology/mtop_util.h"
73 #include "gromacs/topology/topology.h"
74 #include "gromacs/utility/fatalerror.h"
75 #include "gromacs/utility/smalloc.h"
77 // When not built in a configuration with QMMM support, much of this
78 // code is unreachable by design. Tell clang not to warn about it.
79 #pragma GCC diagnostic push
80 #pragma GCC diagnostic ignored "-Wunreachable-code"
81 #pragma GCC diagnostic ignored "-Wmissing-noreturn"
83 /* this struct and these comparison functions are needed for creating
84 * a QMMM input for the QM routines from the QMMM neighbor list.
93 static bool struct_comp(const t_j_particle& a, const t_j_particle& b)
98 static real call_QMroutine(const t_commrec gmx_unused* cr,
99 const t_QMMMrec gmx_unused* qmmm,
100 t_QMrec gmx_unused* qm,
101 t_MMrec gmx_unused* mm,
103 rvec gmx_unused fshift[])
105 /* makes a call to the requested QM routine (qm->QMmethod)
106 * Note that f is actually the gradient, i.e. -f
108 /* do a semi-empiprical calculation */
110 if (qm->QMmethod < eQMmethodRHF && !(mm->nrMMatoms))
116 return call_mopac_SH(qm, mm, f, fshift);
120 return call_mopac(qm, mm, f, fshift);
125 gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac.");
130 /* do an ab-initio calculation */
131 if (qm->bSH && qm->QMmethod == eQMmethodCASSCF)
133 if (GMX_QMMM_GAUSSIAN)
135 return call_gaussian_SH(qmmm, qm, mm, f, fshift);
139 gmx_fatal(FARGS, "Ab-initio Surface-hopping only supported with Gaussian.");
146 return call_gamess(qm, mm, f, fshift);
148 else if (GMX_QMMM_GAUSSIAN)
150 return call_gaussian(qmmm, qm, mm, f, fshift);
152 else if (GMX_QMMM_ORCA)
154 return call_orca(qmmm, qm, mm, f, fshift);
159 "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
165 static void init_QMroutine(const t_commrec gmx_unused* cr, t_QMrec gmx_unused* qm, t_MMrec gmx_unused* mm)
167 /* makes a call to the requested QM routine (qm->QMmethod)
169 if (qm->QMmethod < eQMmethodRHF)
173 /* do a semi-empiprical calculation */
178 gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac.");
183 /* do an ab-initio calculation */
186 init_gamess(cr, qm, mm);
188 else if (GMX_QMMM_GAUSSIAN)
192 else if (GMX_QMMM_ORCA)
198 gmx_fatal(FARGS, "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
201 } /* init_QMroutine */
203 static void update_QMMM_coord(const rvec* x, const t_forcerec* fr, t_QMrec* qm, t_MMrec* mm)
205 /* shifts the QM and MM particles into the central box and stores
206 * these shifted coordinates in the coordinate arrays of the
207 * QMMMrec. These coordinates are passed on the QM subroutines.
211 /* shift the QM atoms into the central box
213 for (i = 0; i < qm->nrQMatoms; i++)
215 rvec_sub(x[qm->indexQM[i]], fr->shift_vec[qm->shiftQM[i]], qm->xQM[i]);
217 /* also shift the MM atoms into the central box, if any
219 for (i = 0; i < mm->nrMMatoms; i++)
221 rvec_sub(x[mm->indexMM[i]], fr->shift_vec[mm->shiftMM[i]], mm->xMM[i]);
223 } /* update_QMMM_coord */
225 /* end of QMMM subroutines */
227 /* QMMM core routines */
229 static t_QMrec* mk_QMrec()
236 static t_MMrec* mk_MMrec()
243 static void init_QMrec(int grpnr, t_QMrec* qm, int nr, const int* atomarray, const gmx_mtop_t* mtop, const t_inputrec* ir)
245 /* fills the t_QMrec struct of QM group grpnr
250 snew(qm->indexQM, nr);
251 snew(qm->shiftQM, nr); /* the shifts */
252 for (int i = 0; i < nr; i++)
254 qm->indexQM[i] = atomarray[i];
257 snew(qm->atomicnumberQM, nr);
259 for (int i = 0; i < qm->nrQMatoms; i++)
261 const t_atom& atom = mtopGetAtomParameters(mtop, qm->indexQM[i], &molb);
262 qm->nelectrons += mtop->atomtypes.atomnumber[atom.type];
263 qm->atomicnumberQM[i] = mtop->atomtypes.atomnumber[atom.type];
266 qm->QMcharge = ir->opts.QMcharge[grpnr];
267 qm->multiplicity = ir->opts.QMmult[grpnr];
268 qm->nelectrons -= ir->opts.QMcharge[grpnr];
270 qm->QMmethod = ir->opts.QMmethod[grpnr];
271 qm->QMbasis = ir->opts.QMbasis[grpnr];
272 /* trajectory surface hopping setup (Gaussian only) */
273 qm->bSH = ir->opts.bSH[grpnr];
274 qm->CASorbitals = ir->opts.CASorbitals[grpnr];
275 qm->CASelectrons = ir->opts.CASelectrons[grpnr];
276 qm->SAsteps = ir->opts.SAsteps[grpnr];
277 qm->SAon = ir->opts.SAon[grpnr];
278 qm->SAoff = ir->opts.SAoff[grpnr];
279 /* hack to prevent gaussian from reinitializing all the time */
280 qm->nQMcpus = 0; /* number of CPU's to be used by g01, is set
281 * upon initializing gaussian
284 /* print the current layer to allow users to check their input */
285 fprintf(stderr, "Layer %d\nnr of QM atoms %d\n", grpnr, nr);
286 fprintf(stderr, "QMlevel: %s/%s\n\n", eQMmethod_names[qm->QMmethod], eQMbasis_names[qm->QMbasis]);
289 static t_QMrec* copy_QMrec(t_QMrec* qm)
291 /* copies the contents of qm into a new t_QMrec struct */
296 qmcopy->nrQMatoms = qm->nrQMatoms;
297 snew(qmcopy->xQM, qmcopy->nrQMatoms);
298 snew(qmcopy->indexQM, qmcopy->nrQMatoms);
299 snew(qmcopy->atomicnumberQM, qm->nrQMatoms);
300 snew(qmcopy->shiftQM, qmcopy->nrQMatoms); /* the shifts */
301 for (i = 0; i < qmcopy->nrQMatoms; i++)
303 qmcopy->shiftQM[i] = qm->shiftQM[i];
304 qmcopy->indexQM[i] = qm->indexQM[i];
305 qmcopy->atomicnumberQM[i] = qm->atomicnumberQM[i];
307 qmcopy->nelectrons = qm->nelectrons;
308 qmcopy->multiplicity = qm->multiplicity;
309 qmcopy->QMcharge = qm->QMcharge;
310 qmcopy->nelectrons = qm->nelectrons;
311 qmcopy->QMmethod = qm->QMmethod;
312 qmcopy->QMbasis = qm->QMbasis;
313 /* trajectory surface hopping setup (Gaussian only) */
314 qmcopy->bSH = qm->bSH;
315 qmcopy->CASorbitals = qm->CASorbitals;
316 qmcopy->CASelectrons = qm->CASelectrons;
317 qmcopy->SAsteps = qm->SAsteps;
318 qmcopy->SAon = qm->SAon;
319 qmcopy->SAoff = qm->SAoff;
321 /* Gaussian init. variables */
322 qmcopy->nQMcpus = qm->nQMcpus;
323 for (i = 0; i < DIM; i++)
325 qmcopy->SHbasis[i] = qm->SHbasis[i];
327 qmcopy->QMmem = qm->QMmem;
328 qmcopy->accuracy = qm->accuracy;
329 qmcopy->cpmcscf = qm->cpmcscf;
330 qmcopy->SAstep = qm->SAstep;
338 t_QMMMrec* mk_QMMMrec()
350 t_QMMMrec* mk_QMMMrec()
352 gmx_incons("Compiled without QMMM");
356 std::vector<int> qmmmAtomIndices(const t_inputrec& ir, const gmx_mtop_t& mtop)
358 const int numQmmmGroups = ir.opts.ngQM;
359 const SimulationGroups& groups = mtop.groups;
360 std::vector<int> qmmmAtoms;
361 for (int i = 0; i < numQmmmGroups; i++)
363 for (const AtomProxy atomP : AtomRange(mtop))
365 int index = atomP.globalAtomNumber();
366 if (getGroupType(groups, SimulationAtomGroupType::QuantumMechanics, index) == i)
368 qmmmAtoms.push_back(index);
371 if (ir.QMMMscheme == eQMMMschemeoniom)
373 /* I assume that users specify the QM groups from small to
374 * big(ger) in the mdp file
376 gmx_mtop_ilistloop_all_t iloop = gmx_mtop_ilistloop_all_init(&mtop);
377 int nral1 = 1 + NRAL(F_VSITE2);
379 while (const InteractionLists* ilists = gmx_mtop_ilistloop_all_next(iloop, &atomOffset))
381 const InteractionList& ilist = (*ilists)[F_VSITE2];
382 for (int j = 0; j < ilist.size(); j += nral1)
384 const int vsite = atomOffset + ilist.iatoms[j]; /* the vsite */
385 const int ai = atomOffset + ilist.iatoms[j + 1]; /* constructing atom */
386 const int aj = atomOffset + ilist.iatoms[j + 2]; /* constructing atom */
387 if (getGroupType(groups, SimulationAtomGroupType::QuantumMechanics, vsite)
388 == getGroupType(groups, SimulationAtomGroupType::QuantumMechanics, ai)
389 && getGroupType(groups, SimulationAtomGroupType::QuantumMechanics, vsite)
390 == getGroupType(groups, SimulationAtomGroupType::QuantumMechanics, aj))
392 /* this dummy link atom needs to be removed from qmmmAtoms
393 * before making the QMrec of this layer!
395 qmmmAtoms.erase(std::remove_if(qmmmAtoms.begin(), qmmmAtoms.end(),
396 [&vsite](int atom) { return atom == vsite; }),
406 void removeQmmmAtomCharges(gmx_mtop_t* mtop, gmx::ArrayRef<const int> qmmmAtoms)
409 for (gmx::index i = 0; i < qmmmAtoms.ssize(); i++)
412 mtopGetMolblockIndex(mtop, qmmmAtoms[i], &molb, nullptr, &indexInMolecule);
413 t_atom* atom = &mtop->moltype[mtop->molblock[molb].type].atoms.atom[indexInMolecule];
419 void init_QMMMrec(const t_commrec* cr, const gmx_mtop_t* mtop, const t_inputrec* ir, const t_forcerec* fr)
421 /* we put the atomsnumbers of atoms that belong to the QMMM group in
422 * an array that will be copied later to QMMMrec->indexQM[..]. Also
423 * it will be used to create an QMMMrec->bQMMM index array that
424 * simply contains true/false for QM and MM (the other) atoms.
432 gmx_incons("Compiled without QMMM");
435 if (ir->cutoff_scheme != ecutsGROUP)
437 gmx_fatal(FARGS, "QMMM is currently only supported with cutoff-scheme=group");
439 if (!EI_DYNAMICS(ir->eI))
441 gmx_fatal(FARGS, "QMMM is only supported with dynamics");
444 /* issue a fatal if the user wants to run with more than one node */
447 gmx_fatal(FARGS, "QM/MM does not work in parallel, use a single rank instead\n");
450 /* Make a local copy of the QMMMrec */
453 /* bQMMM[..] is an array containing TRUE/FALSE for atoms that are
454 * QM/not QM. We first set all elemenst at false. Afterwards we use
455 * the qm_arr (=MMrec->indexQM) to changes the elements
456 * corresponding to the QM atoms at TRUE. */
458 qr->QMMMscheme = ir->QMMMscheme;
460 /* we take the possibility into account that a user has
461 * defined more than one QM group:
463 /* an ugly work-around in case there is only one group In this case
464 * the whole system is treated as QM. Otherwise the second group is
465 * always the rest of the total system and is treated as MM.
468 /* small problem if there is only QM.... so no MM */
470 int numQmmmGroups = ir->opts.ngQM;
472 if (qr->QMMMscheme == eQMMMschemeoniom)
474 qr->nrQMlayers = numQmmmGroups;
481 /* there are numQmmmGroups groups of QM atoms. In case of multiple QM groups
482 * I assume that the users wants to do ONIOM. However, maybe it
483 * should also be possible to define more than one QM subsystem with
484 * independent neighbourlists. I have to think about
487 std::vector<int> qmmmAtoms = qmmmAtomIndices(*ir, *mtop);
488 snew(qr->qm, numQmmmGroups);
489 for (int i = 0; i < numQmmmGroups; i++)
492 if (qr->QMMMscheme == eQMMMschemeoniom)
494 /* add the atoms to the bQMMM array
497 /* I assume that users specify the QM groups from small to
498 * big(ger) in the mdp file
500 qr->qm[i] = mk_QMrec();
501 /* store QM atoms in this layer in the QMrec and initialise layer
503 init_QMrec(i, qr->qm[i], qmmmAtoms.size(), qmmmAtoms.data(), mtop, ir);
506 if (qr->QMMMscheme != eQMMMschemeoniom)
509 /* standard QMMM, all layers are merged together so there is one QM
510 * subsystem and one MM subsystem.
511 * Also we set the charges to zero in mtop to prevent the innerloops
512 * from doubly counting the electostatic QM MM interaction
513 * TODO: Consider doing this in grompp instead.
516 qr->qm[0] = mk_QMrec();
517 /* store QM atoms in the QMrec and initialise
519 init_QMrec(0, qr->qm[0], qmmmAtoms.size(), qmmmAtoms.data(), mtop, ir);
521 /* MM rec creation */
523 mm->scalefactor = ir->scalefactor;
524 mm->nrMMatoms = (mtop->natoms) - (qr->qm[0]->nrQMatoms); /* rest of the atoms */
528 { /* MM rec creation */
530 mm->scalefactor = ir->scalefactor;
535 /* these variables get updated in the update QMMMrec */
537 if (qr->nrQMlayers == 1)
539 /* with only one layer there is only one initialisation
540 * needed. Multilayer is a bit more complicated as it requires
541 * re-initialisation at every step of the simulation. This is due
542 * to the use of COMMON blocks in the fortran QM subroutines.
544 if (qr->qm[0]->QMmethod < eQMmethodRHF)
548 /* semi-empiprical 1-layer ONIOM calculation requested (mopac93) */
549 init_mopac(qr->qm[0]);
553 gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac.");
558 /* ab initio calculation requested (gamess/gaussian/ORCA) */
561 init_gamess(cr, qr->qm[0], qr->mm);
563 else if (GMX_QMMM_GAUSSIAN)
565 init_gaussian(qr->qm[0]);
567 else if (GMX_QMMM_ORCA)
569 init_orca(qr->qm[0]);
574 "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
580 void update_QMMMrec(const t_commrec* cr, const t_forcerec* fr, const rvec* x, const t_mdatoms* md, const matrix box)
582 /* updates the coordinates of both QM atoms and MM atoms and stores
583 * them in the QMMMrec.
585 * NOTE: is NOT yet working if there are no PBC. Also in ns.c, simple
586 * ns needs to be fixed!
588 int mm_max = 0, mm_nr = 0, mm_nr_new, i, j, is, k, shift;
589 t_j_particle *mm_j_particles = nullptr, *qm_i_particles = nullptr;
597 int* parallelMMarray = nullptr;
601 gmx_incons("Compiled without QMMM");
604 /* every cpu has this array. On every processor we fill this array
605 * with 1's and 0's. 1's indicate the atoms is a QM atom on the
606 * current cpu in a later stage these arrays are all summed. indexes
607 * > 0 indicate the atom is a QM atom. Every node therefore knows
608 * whcih atoms are part of the QM subsystem.
610 /* copy some pointers */
613 QMMMlist = fr->QMMMlist;
615 /* init_pbc(box); needs to be called first, see pbc.h */
617 clear_ivec(null_ivec);
618 set_pbc_dd(&pbc, fr->ePBC, DOMAINDECOMP(cr) ? cr->dd->numCells : null_ivec, FALSE, box);
619 /* only in standard (normal) QMMM we need the neighbouring MM
620 * particles to provide a electric field of point charges for the QM
623 if (qr->QMMMscheme == eQMMMschemenormal) /* also implies 1 QM-layer */
625 /* we NOW create/update a number of QMMMrec entries:
627 * 1) the shiftQM, containing the shifts of the QM atoms
629 * 2) the indexMM array, containing the index of the MM atoms
631 * 3) the shiftMM, containing the shifts of the MM atoms
633 * 4) the shifted coordinates of the MM atoms
635 * the shifts are used for computing virial of the QM/MM particles.
637 qm = qr->qm[0]; /* in case of normal QMMM, there is only one group */
638 snew(qm_i_particles, QMMMlist->nri);
641 qm_i_particles[0].shift = XYZ2IS(0, 0, 0);
642 for (i = 0; i < QMMMlist->nri; i++)
644 qm_i_particles[i].j = QMMMlist->iinr[i];
648 qm_i_particles[i].shift =
649 pbc_dx_aiuc(&pbc, x[QMMMlist->iinr[0]], x[QMMMlist->iinr[i]], dx);
651 /* However, since nri >= nrQMatoms, we do a quicksort, and throw
652 * out double, triple, etc. entries later, as we do for the MM
656 /* compute the shift for the MM j-particles with respect to
657 * the QM i-particle and store them.
660 crd[0] = IS2X(QMMMlist->shift[i]) + IS2X(qm_i_particles[i].shift);
661 crd[1] = IS2Y(QMMMlist->shift[i]) + IS2Y(qm_i_particles[i].shift);
662 crd[2] = IS2Z(QMMMlist->shift[i]) + IS2Z(qm_i_particles[i].shift);
663 is = XYZ2IS(crd[0], crd[1], crd[2]);
664 for (j = QMMMlist->jindex[i]; j < QMMMlist->jindex[i + 1]; j++)
669 srenew(mm_j_particles, mm_max);
672 mm_j_particles[mm_nr].j = QMMMlist->jjnr[j];
673 mm_j_particles[mm_nr].shift = is;
678 /* quicksort QM and MM shift arrays and throw away multiple entries */
681 std::sort(qm_i_particles, qm_i_particles + QMMMlist->nri, struct_comp);
682 /* The mm_j_particles argument to qsort is not allowed to be nullptr */
685 std::sort(mm_j_particles, mm_j_particles + mm_nr, struct_comp);
687 /* remove multiples in the QM shift array, since in init_QMMM() we
688 * went through the atom numbers from 0 to md.nr, the order sorted
689 * here matches the one of QMindex already.
692 for (i = 0; i < QMMMlist->nri; i++)
694 if (i == 0 || qm_i_particles[i].j != qm_i_particles[i - 1].j)
696 qm_i_particles[j++] = qm_i_particles[i];
700 /* Remove double entries for the MM array.
701 * Also remove mm atoms that have no charges!
702 * actually this is already done in the ns.c
704 for (i = 0; i < mm_nr; i++)
706 if ((i == 0 || mm_j_particles[i].j != mm_j_particles[i - 1].j)
707 && !md->bQM[mm_j_particles[i].j]
708 && ((md->chargeA[mm_j_particles[i].j] != 0.0_real)
709 || (md->chargeB && (md->chargeB[mm_j_particles[i].j] != 0.0_real))))
711 mm_j_particles[mm_nr_new++] = mm_j_particles[i];
715 /* store the data retrieved above into the QMMMrec
718 /* Keep the compiler happy,
719 * shift will always be set in the loop for i=0
722 for (i = 0; i < qm->nrQMatoms; i++)
724 /* not all qm particles might have appeared as i
725 * particles. They might have been part of the same charge
726 * group for instance.
728 if (qm->indexQM[i] == qm_i_particles[k].j)
730 shift = qm_i_particles[k++].shift;
732 /* use previous shift, assuming they belong the same charge
736 qm->shiftQM[i] = shift;
739 /* parallel excecution */
742 snew(parallelMMarray, 2 * (md->nr));
743 /* only MM particles have a 1 at their atomnumber. The second part
744 * of the array contains the shifts. Thus:
745 * p[i]=1/0 depending on wether atomnumber i is a MM particle in the QM
746 * step or not. p[i+md->nr] is the shift of atomnumber i.
748 for (i = 0; i < 2 * (md->nr); i++)
750 parallelMMarray[i] = 0;
753 for (i = 0; i < mm_nr; i++)
755 parallelMMarray[mm_j_particles[i].j] = 1;
756 parallelMMarray[mm_j_particles[i].j + (md->nr)] = mm_j_particles[i].shift;
758 gmx_sumi(md->nr, parallelMMarray, cr);
762 for (i = 0; i < md->nr; i++)
764 if (parallelMMarray[i])
769 srenew(mm->indexMM, mm_max);
770 srenew(mm->shiftMM, mm_max);
772 mm->indexMM[mm_nr] = i;
773 mm->shiftMM[mm_nr++] = parallelMMarray[i + md->nr] / parallelMMarray[i];
776 mm->nrMMatoms = mm_nr;
777 free(parallelMMarray);
779 /* serial execution */
782 mm->nrMMatoms = mm_nr;
783 srenew(mm->shiftMM, mm_nr);
784 srenew(mm->indexMM, mm_nr);
785 for (i = 0; i < mm_nr; i++)
787 mm->indexMM[i] = mm_j_particles[i].j;
788 mm->shiftMM[i] = mm_j_particles[i].shift;
791 /* (re) allocate memory for the MM coordiate array. The QM
792 * coordinate array was already allocated in init_QMMM, and is
793 * only (re)filled in the update_QMMM_coordinates routine
795 srenew(mm->xMM, mm->nrMMatoms);
796 /* now we (re) fill the array that contains the MM charges with
797 * the forcefield charges. If requested, these charges will be
800 srenew(mm->MMcharges, mm->nrMMatoms);
801 for (i = 0; i < mm->nrMMatoms; i++) /* no free energy yet */
803 mm->MMcharges[i] = md->chargeA[mm->indexMM[i]] * mm->scalefactor;
805 /* the next routine fills the coordinate fields in the QMMM rec of
806 * both the qunatum atoms and the MM atoms, using the shifts
810 update_QMMM_coord(x, fr, qr->qm[0], qr->mm);
811 free(qm_i_particles);
812 free(mm_j_particles);
814 else /* ONIOM */ /* ????? */
817 /* do for each layer */
818 for (j = 0; j < qr->nrQMlayers; j++)
821 qm->shiftQM[0] = XYZ2IS(0, 0, 0);
822 for (i = 1; i < qm->nrQMatoms; i++)
824 qm->shiftQM[i] = pbc_dx_aiuc(&pbc, x[qm->indexQM[0]], x[qm->indexQM[i]], dx);
826 update_QMMM_coord(x, fr, qm, mm);
829 } /* update_QMMM_rec */
831 real calculate_QMMM(const t_commrec* cr, gmx::ForceWithShiftForces* forceWithShiftForces, const t_QMMMrec* qr)
834 /* a selection for the QM package depending on which is requested
835 * (Gaussian, GAMESS-UK, MOPAC or ORCA) needs to be implemented here. Now
836 * it works through defines.... Not so nice yet
839 t_MMrec* mm = nullptr;
840 rvec * forces = nullptr, *fshift = nullptr, *forces2 = nullptr,
841 *fshift2 = nullptr; /* needed for multilayer ONIOM */
846 gmx_incons("Compiled without QMMM");
849 /* make a local copy the QMMMrec pointer
853 /* now different procedures are carried out for one layer ONION and
854 * normal QMMM on one hand and multilayer oniom on the other
856 gmx::ArrayRef<gmx::RVec> fMM = forceWithShiftForces->force();
857 gmx::ArrayRef<gmx::RVec> fshiftMM = forceWithShiftForces->shiftForces();
858 if (qr->QMMMscheme == eQMMMschemenormal || qr->nrQMlayers == 1)
861 snew(forces, (qm->nrQMatoms + mm->nrMMatoms));
862 snew(fshift, (qm->nrQMatoms + mm->nrMMatoms));
863 QMener = call_QMroutine(cr, qr, qm, mm, forces, fshift);
864 for (i = 0; i < qm->nrQMatoms; i++)
866 for (j = 0; j < DIM; j++)
868 fMM[qm->indexQM[i]][j] -= forces[i][j];
869 fshiftMM[qm->shiftQM[i]][j] += fshift[i][j];
872 for (i = 0; i < mm->nrMMatoms; i++)
874 for (j = 0; j < DIM; j++)
876 fMM[mm->indexMM[i]][j] -= forces[qm->nrQMatoms + i][j];
877 fshiftMM[mm->shiftMM[i]][j] += fshift[qm->nrQMatoms + i][j];
883 else /* Multi-layer ONIOM */
885 for (i = 0; i < qr->nrQMlayers - 1; i++) /* last layer is special */
888 qm2 = copy_QMrec(qr->qm[i + 1]);
890 qm2->nrQMatoms = qm->nrQMatoms;
892 for (j = 0; j < qm2->nrQMatoms; j++)
894 for (k = 0; k < DIM; k++)
896 qm2->xQM[j][k] = qm->xQM[j][k];
898 qm2->indexQM[j] = qm->indexQM[j];
899 qm2->atomicnumberQM[j] = qm->atomicnumberQM[j];
900 qm2->shiftQM[j] = qm->shiftQM[j];
903 qm2->QMcharge = qm->QMcharge;
904 /* this layer at the higher level of theory */
905 srenew(forces, qm->nrQMatoms);
906 srenew(fshift, qm->nrQMatoms);
907 /* we need to re-initialize the QMroutine every step... */
908 init_QMroutine(cr, qm, mm);
909 QMener += call_QMroutine(cr, qr, qm, mm, forces, fshift);
911 /* this layer at the lower level of theory */
912 srenew(forces2, qm->nrQMatoms);
913 srenew(fshift2, qm->nrQMatoms);
914 init_QMroutine(cr, qm2, mm);
915 QMener -= call_QMroutine(cr, qr, qm2, mm, forces2, fshift2);
916 /* E = E1high-E1low The next layer includes the current layer at
917 * the lower level of theory, which provides + E2low
918 * this is similar for gradients
920 for (i = 0; i < qm->nrQMatoms; i++)
922 for (j = 0; j < DIM; j++)
924 fMM[qm->indexQM[i]][j] -= (forces[i][j] - forces2[i][j]);
925 fshiftMM[qm->shiftQM[i]][j] += (fshift[i][j] - fshift2[i][j]);
930 /* now the last layer still needs to be done: */
931 qm = qr->qm[qr->nrQMlayers - 1]; /* C counts from 0 */
932 init_QMroutine(cr, qm, mm);
933 srenew(forces, qm->nrQMatoms);
934 srenew(fshift, qm->nrQMatoms);
935 QMener += call_QMroutine(cr, qr, qm, mm, forces, fshift);
936 for (i = 0; i < qm->nrQMatoms; i++)
938 for (j = 0; j < DIM; j++)
940 fMM[qm->indexQM[i]][j] -= forces[i][j];
941 fshiftMM[qm->shiftQM[i]][j] += fshift[i][j];
950 } /* calculate_QMMM */
952 #pragma GCC diagnostic pop