src/gromacs/mdlib/qmmm.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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,2018,2019, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "qmmm.h"
  40
  41 #include "config.h"
  42
  43 #include <cmath>
  44 #include <cstdio>
  45 #include <cstdlib>
  46 #include <cstring>
  47
  48 #include <algorithm>
  49
  50 #include "gromacs/domdec/domdec_struct.h"
  51 #include "gromacs/fileio/confio.h"
  52 #include "gromacs/gmxlib/network.h"
  53 #include "gromacs/gmxlib/nrnb.h"
  54 #include "gromacs/math/functions.h"
  55 #include "gromacs/math/units.h"
  56 #include "gromacs/math/vec.h"
  57 #include "gromacs/mdlib/qm_gamess.h"
  58 #include "gromacs/mdlib/qm_gaussian.h"
  59 #include "gromacs/mdlib/qm_mopac.h"
  60 #include "gromacs/mdlib/qm_orca.h"
  61 #include "gromacs/mdtypes/commrec.h"
  62 #include "gromacs/mdtypes/forceoutput.h"
  63 #include "gromacs/mdtypes/forcerec.h"
  64 #include "gromacs/mdtypes/inputrec.h"
  65 #include "gromacs/mdtypes/md_enums.h"
  66 #include "gromacs/mdtypes/mdatom.h"
  67 #include "gromacs/mdtypes/nblist.h"
  68 #include "gromacs/pbcutil/ishift.h"
  69 #include "gromacs/pbcutil/pbc.h"
  70 #include "gromacs/topology/mtop_lookup.h"
  71 #include "gromacs/topology/mtop_util.h"
  72 #include "gromacs/topology/topology.h"
  73 #include "gromacs/utility/fatalerror.h"
  74 #include "gromacs/utility/smalloc.h"
  75
  76 // When not built in a configuration with QMMM support, much of this
  77 // code is unreachable by design. Tell clang not to warn about it.
  78 #pragma GCC diagnostic push
  79 #pragma GCC diagnostic ignored "-Wunreachable-code"
  80 #pragma GCC diagnostic ignored "-Wmissing-noreturn"
  81
  82 /* this struct and these comparison functions are needed for creating
  83  * a QMMM input for the QM routines from the QMMM neighbor list.
  84  */
  85
  86 typedef struct
  87 {
  88     int j;
  89     int shift;
  90 } t_j_particle;
  91
  92 static bool struct_comp(const t_j_particle& a, const t_j_particle& b)
  93 {
  94     return a.j < b.j;
  95 }
  96
  97 static real call_QMroutine(const t_commrec gmx_unused* cr,
  98                            const t_QMMMrec gmx_unused* qmmm,
  99                            t_QMrec gmx_unused* qm,
 100                            t_MMrec gmx_unused* mm,
 101                            rvec gmx_unused f[],
 102                            rvec gmx_unused fshift[])
 103 {
 104     /* makes a call to the requested QM routine (qm->QMmethod)
 105      * Note that f is actually the gradient, i.e. -f
 106      */
 107     /* do a semi-empiprical calculation */
 108
 109     if (qm->QMmethod < eQMmethodRHF && !(mm->nrMMatoms))
 110     {
 111         if (GMX_QMMM_MOPAC)
 112         {
 113             if (qm->bSH)
 114             {
 115                 return call_mopac_SH(qm, mm, f, fshift);
 116             }
 117             else
 118             {
 119                 return call_mopac(qm, mm, f, fshift);
 120             }
 121         }
 122         else
 123         {
 124             gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac.");
 125         }
 126     }
 127     else
 128     {
 129         /* do an ab-initio calculation */
 130         if (qm->bSH && qm->QMmethod == eQMmethodCASSCF)
 131         {
 132             if (GMX_QMMM_GAUSSIAN)
 133             {
 134                 return call_gaussian_SH(qmmm, qm, mm, f, fshift);
 135             }
 136             else
 137             {
 138                 gmx_fatal(FARGS, "Ab-initio Surface-hopping only supported with Gaussian.");
 139             }
 140         }
 141         else
 142         {
 143             if (GMX_QMMM_GAMESS)
 144             {
 145                 return call_gamess(qm, mm, f, fshift);
 146             }
 147             else if (GMX_QMMM_GAUSSIAN)
 148             {
 149                 return call_gaussian(qmmm, qm, mm, f, fshift);
 150             }
 151             else if (GMX_QMMM_ORCA)
 152             {
 153                 return call_orca(qmmm, qm, mm, f, fshift);
 154             }
 155             else
 156             {
 157                 gmx_fatal(FARGS,
 158                           "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
 159             }
 160         }
 161     }
 162 }
 163
 164 static void init_QMroutine(const t_commrec gmx_unused* cr, t_QMrec gmx_unused* qm, t_MMrec gmx_unused* mm)
 165 {
 166     /* makes a call to the requested QM routine (qm->QMmethod)
 167      */
 168     if (qm->QMmethod < eQMmethodRHF)
 169     {
 170         if (GMX_QMMM_MOPAC)
 171         {
 172             /* do a semi-empiprical calculation */
 173             init_mopac(qm);
 174         }
 175         else
 176         {
 177             gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac.");
 178         }
 179     }
 180     else
 181     {
 182         /* do an ab-initio calculation */
 183         if (GMX_QMMM_GAMESS)
 184         {
 185             init_gamess(cr, qm, mm);
 186         }
 187         else if (GMX_QMMM_GAUSSIAN)
 188         {
 189             init_gaussian(qm);
 190         }
 191         else if (GMX_QMMM_ORCA)
 192         {
 193             init_orca(qm);
 194         }
 195         else
 196         {
 197             gmx_fatal(FARGS, "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
 198         }
 199     }
 200 } /* init_QMroutine */
 201
 202 static void update_QMMM_coord(const rvec* x, const t_forcerec* fr, t_QMrec* qm, t_MMrec* mm)
 203 {
 204     /* shifts the QM and MM particles into the central box and stores
 205      * these shifted coordinates in the coordinate arrays of the
 206      * QMMMrec. These coordinates are passed on the QM subroutines.
 207      */
 208     int i;
 209
 210     /* shift the QM atoms into the central box
 211      */
 212     for (i = 0; i < qm->nrQMatoms; i++)
 213     {
 214         rvec_sub(x[qm->indexQM[i]], fr->shift_vec[qm->shiftQM[i]], qm->xQM[i]);
 215     }
 216     /* also shift the MM atoms into the central box, if any
 217      */
 218     for (i = 0; i < mm->nrMMatoms; i++)
 219     {
 220         rvec_sub(x[mm->indexMM[i]], fr->shift_vec[mm->shiftMM[i]], mm->xMM[i]);
 221     }
 222 } /* update_QMMM_coord */
 223
 224 /* end of QMMM subroutines */
 225
 226 /* QMMM core routines */
 227
 228 static t_QMrec* mk_QMrec()
 229 {
 230     t_QMrec* qm;
 231     snew(qm, 1);
 232     return qm;
 233 } /* mk_QMrec */
 234
 235 static t_MMrec* mk_MMrec()
 236 {
 237     t_MMrec* mm;
 238     snew(mm, 1);
 239     return mm;
 240 } /* mk_MMrec */
 241
 242 static void init_QMrec(int grpnr, t_QMrec* qm, int nr, const int* atomarray, const gmx_mtop_t* mtop, const t_inputrec* ir)
 243 {
 244     /* fills the t_QMrec struct of QM group grpnr
 245      */
 246
 247     qm->nrQMatoms = nr;
 248     snew(qm->xQM, nr);
 249     snew(qm->indexQM, nr);
 250     snew(qm->shiftQM, nr); /* the shifts */
 251     for (int i = 0; i < nr; i++)
 252     {
 253         qm->indexQM[i] = atomarray[i];
 254     }
 255
 256     snew(qm->atomicnumberQM, nr);
 257     int molb = 0;
 258     for (int i = 0; i < qm->nrQMatoms; i++)
 259     {
 260         const t_atom& atom = mtopGetAtomParameters(mtop, qm->indexQM[i], &molb);
 261         qm->nelectrons += mtop->atomtypes.atomnumber[atom.type];
 262         qm->atomicnumberQM[i] = mtop->atomtypes.atomnumber[atom.type];
 263     }
 264
 265     qm->QMcharge     = ir->opts.QMcharge[grpnr];
 266     qm->multiplicity = ir->opts.QMmult[grpnr];
 267     qm->nelectrons -= ir->opts.QMcharge[grpnr];
 268
 269     qm->QMmethod = ir->opts.QMmethod[grpnr];
 270     qm->QMbasis  = ir->opts.QMbasis[grpnr];
 271     /* trajectory surface hopping setup (Gaussian only) */
 272     qm->bSH          = ir->opts.bSH[grpnr];
 273     qm->CASorbitals  = ir->opts.CASorbitals[grpnr];
 274     qm->CASelectrons = ir->opts.CASelectrons[grpnr];
 275     qm->SAsteps      = ir->opts.SAsteps[grpnr];
 276     qm->SAon         = ir->opts.SAon[grpnr];
 277     qm->SAoff        = ir->opts.SAoff[grpnr];
 278     /* hack to prevent gaussian from reinitializing all the time */
 279     qm->nQMcpus = 0; /* number of CPU's to be used by g01, is set
 280                       * upon initializing gaussian
 281                       * (init_gaussian()
 282                       */
 283     /* print the current layer to allow users to check their input */
 284     fprintf(stderr, "Layer %d\nnr of QM atoms %d\n", grpnr, nr);
 285     fprintf(stderr, "QMlevel: %s/%s\n\n", eQMmethod_names[qm->QMmethod], eQMbasis_names[qm->QMbasis]);
 286 } /* init_QMrec */
 287
 288 static t_QMrec* copy_QMrec(t_QMrec* qm)
 289 {
 290     /* copies the contents of qm into a new t_QMrec struct */
 291     t_QMrec* qmcopy;
 292     int      i;
 293
 294     qmcopy            = mk_QMrec();
 295     qmcopy->nrQMatoms = qm->nrQMatoms;
 296     snew(qmcopy->xQM, qmcopy->nrQMatoms);
 297     snew(qmcopy->indexQM, qmcopy->nrQMatoms);
 298     snew(qmcopy->atomicnumberQM, qm->nrQMatoms);
 299     snew(qmcopy->shiftQM, qmcopy->nrQMatoms); /* the shifts */
 300     for (i = 0; i < qmcopy->nrQMatoms; i++)
 301     {
 302         qmcopy->shiftQM[i]        = qm->shiftQM[i];
 303         qmcopy->indexQM[i]        = qm->indexQM[i];
 304         qmcopy->atomicnumberQM[i] = qm->atomicnumberQM[i];
 305     }
 306     qmcopy->nelectrons   = qm->nelectrons;
 307     qmcopy->multiplicity = qm->multiplicity;
 308     qmcopy->QMcharge     = qm->QMcharge;
 309     qmcopy->nelectrons   = qm->nelectrons;
 310     qmcopy->QMmethod     = qm->QMmethod;
 311     qmcopy->QMbasis      = qm->QMbasis;
 312     /* trajectory surface hopping setup (Gaussian only) */
 313     qmcopy->bSH          = qm->bSH;
 314     qmcopy->CASorbitals  = qm->CASorbitals;
 315     qmcopy->CASelectrons = qm->CASelectrons;
 316     qmcopy->SAsteps      = qm->SAsteps;
 317     qmcopy->SAon         = qm->SAon;
 318     qmcopy->SAoff        = qm->SAoff;
 319
 320     /* Gaussian init. variables */
 321     qmcopy->nQMcpus = qm->nQMcpus;
 322     for (i = 0; i < DIM; i++)
 323     {
 324         qmcopy->SHbasis[i] = qm->SHbasis[i];
 325     }
 326     qmcopy->QMmem    = qm->QMmem;
 327     qmcopy->accuracy = qm->accuracy;
 328     qmcopy->cpmcscf  = qm->cpmcscf;
 329     qmcopy->SAstep   = qm->SAstep;
 330
 331     return (qmcopy);
 332
 333 } /*copy_QMrec */
 334
 335 #if GMX_QMMM
 336
 337 t_QMMMrec* mk_QMMMrec()
 338 {
 339     t_QMMMrec* qr;
 340
 341     snew(qr, 1);
 342
 343     return qr;
 344
 345 } /* mk_QMMMrec */
 346
 347 #else /* GMX_QMMM */
 348
 349 t_QMMMrec* mk_QMMMrec()
 350 {
 351     gmx_incons("Compiled without QMMM");
 352 } /* mk_QMMMrec */
 353 #endif
 354
 355 std::vector<int> qmmmAtomIndices(const t_inputrec& ir, const gmx_mtop_t& mtop)
 356 {
 357     const int               numQmmmGroups = ir.opts.ngQM;
 358     const SimulationGroups& groups        = mtop.groups;
 359     std::vector<int>        qmmmAtoms;
 360     for (int i = 0; i < numQmmmGroups; i++)
 361     {
 362         for (const AtomProxy atomP : AtomRange(mtop))
 363         {
 364             int index = atomP.globalAtomNumber();
 365             if (getGroupType(groups, SimulationAtomGroupType::QuantumMechanics, index) == i)
 366             {
 367                 qmmmAtoms.push_back(index);
 368             }
 369         }
 370         if (ir.QMMMscheme == eQMMMschemeoniom)
 371         {
 372             /* I assume that users specify the QM groups from small to
 373              * big(ger) in the mdp file
 374              */
 375             gmx_mtop_ilistloop_all_t iloop      = gmx_mtop_ilistloop_all_init(&mtop);
 376             int                      nral1      = 1 + NRAL(F_VSITE2);
 377             int                      atomOffset = 0;
 378             while (const InteractionLists* ilists = gmx_mtop_ilistloop_all_next(iloop, &atomOffset))
 379             {
 380                 const InteractionList& ilist = (*ilists)[F_VSITE2];
 381                 for (int j = 0; j < ilist.size(); j += nral1)
 382                 {
 383                     const int vsite = atomOffset + ilist.iatoms[j];     /* the vsite         */
 384                     const int ai    = atomOffset + ilist.iatoms[j + 1]; /* constructing atom */
 385                     const int aj    = atomOffset + ilist.iatoms[j + 2]; /* constructing atom */
 386                     if (getGroupType(groups, SimulationAtomGroupType::QuantumMechanics, vsite)
 387                                 == getGroupType(groups, SimulationAtomGroupType::QuantumMechanics, ai)
 388                         && getGroupType(groups, SimulationAtomGroupType::QuantumMechanics, vsite)
 389                                    == getGroupType(groups, SimulationAtomGroupType::QuantumMechanics, aj))
 390                     {
 391                         /* this dummy link atom needs to be removed from qmmmAtoms
 392                          * before making the QMrec of this layer!
 393                          */
 394                         qmmmAtoms.erase(std::remove_if(qmmmAtoms.begin(), qmmmAtoms.end(),
 395                                                        [&vsite](int atom) { return atom == vsite; }),
 396                                         qmmmAtoms.end());
 397                     }
 398                 }
 399             }
 400         }
 401     }
 402     return qmmmAtoms;
 403 }
 404
 405 void removeQmmmAtomCharges(gmx_mtop_t* mtop, gmx::ArrayRef<const int> qmmmAtoms)
 406 {
 407     int molb = 0;
 408     for (gmx::index i = 0; i < qmmmAtoms.ssize(); i++)
 409     {
 410         int indexInMolecule;
 411         mtopGetMolblockIndex(mtop, qmmmAtoms[i], &molb, nullptr, &indexInMolecule);
 412         t_atom* atom = &mtop->moltype[mtop->molblock[molb].type].atoms.atom[indexInMolecule];
 413         atom->q      = 0.0;
 414         atom->qB     = 0.0;
 415     }
 416 }
 417
 418 void init_QMMMrec(const t_commrec* cr, const gmx_mtop_t* mtop, const t_inputrec* ir, const t_forcerec* fr)
 419 {
 420     /* we put the atomsnumbers of atoms that belong to the QMMM group in
 421      * an array that will be copied later to QMMMrec->indexQM[..]. Also
 422      * it will be used to create an QMMMrec->bQMMM index array that
 423      * simply contains true/false for QM and MM (the other) atoms.
 424      */
 425
 426     t_QMMMrec* qr;
 427     t_MMrec*   mm;
 428
 429     if (!GMX_QMMM)
 430     {
 431         gmx_incons("Compiled without QMMM");
 432     }
 433
 434     if (ir->cutoff_scheme != ecutsGROUP)
 435     {
 436         gmx_fatal(FARGS, "QMMM is currently only supported with cutoff-scheme=group");
 437     }
 438     if (!EI_DYNAMICS(ir->eI))
 439     {
 440         gmx_fatal(FARGS, "QMMM is only supported with dynamics");
 441     }
 442
 443     /* issue a fatal if the user wants to run with more than one node */
 444     if (PAR(cr))
 445     {
 446         gmx_fatal(FARGS, "QM/MM does not work in parallel, use a single rank instead\n");
 447     }
 448
 449     /* Make a local copy of the QMMMrec */
 450     qr = fr->qr;
 451
 452     /* bQMMM[..] is an array containing TRUE/FALSE for atoms that are
 453      * QM/not QM. We first set all elemenst at false. Afterwards we use
 454      * the qm_arr (=MMrec->indexQM) to changes the elements
 455      * corresponding to the QM atoms at TRUE.  */
 456
 457     qr->QMMMscheme = ir->QMMMscheme;
 458
 459     /* we take the possibility into account that a user has
 460      * defined more than one QM group:
 461      */
 462     /* an ugly work-around in case there is only one group In this case
 463      * the whole system is treated as QM. Otherwise the second group is
 464      * always the rest of the total system and is treated as MM.
 465      */
 466
 467     /* small problem if there is only QM.... so no MM */
 468
 469     int numQmmmGroups = ir->opts.ngQM;
 470
 471     if (qr->QMMMscheme == eQMMMschemeoniom)
 472     {
 473         qr->nrQMlayers = numQmmmGroups;
 474     }
 475     else
 476     {
 477         qr->nrQMlayers = 1;
 478     }
 479
 480     /* there are numQmmmGroups groups of QM atoms. In case of multiple QM groups
 481      * I assume that the users wants to do ONIOM. However, maybe it
 482      * should also be possible to define more than one QM subsystem with
 483      * independent neighbourlists. I have to think about
 484      * that.. 11-11-2003
 485      */
 486     std::vector<int> qmmmAtoms = qmmmAtomIndices(*ir, *mtop);
 487     snew(qr->qm, numQmmmGroups);
 488     for (int i = 0; i < numQmmmGroups; i++)
 489     {
 490         /* new layer */
 491         if (qr->QMMMscheme == eQMMMschemeoniom)
 492         {
 493             /* add the atoms to the bQMMM array
 494              */
 495
 496             /* I assume that users specify the QM groups from small to
 497              * big(ger) in the mdp file
 498              */
 499             qr->qm[i] = mk_QMrec();
 500             /* store QM atoms in this layer in the QMrec and initialise layer
 501              */
 502             init_QMrec(i, qr->qm[i], qmmmAtoms.size(), qmmmAtoms.data(), mtop, ir);
 503         }
 504     }
 505     if (qr->QMMMscheme != eQMMMschemeoniom)
 506     {
 507
 508         /* standard QMMM, all layers are merged together so there is one QM
 509          * subsystem and one MM subsystem.
 510          * Also we set the charges to zero in mtop to prevent the innerloops
 511          * from doubly counting the electostatic QM MM interaction
 512          * TODO: Consider doing this in grompp instead.
 513          */
 514
 515         qr->qm[0] = mk_QMrec();
 516         /* store QM atoms in the QMrec and initialise
 517          */
 518         init_QMrec(0, qr->qm[0], qmmmAtoms.size(), qmmmAtoms.data(), mtop, ir);
 519
 520         /* MM rec creation */
 521         mm              = mk_MMrec();
 522         mm->scalefactor = ir->scalefactor;
 523         mm->nrMMatoms   = (mtop->natoms) - (qr->qm[0]->nrQMatoms); /* rest of the atoms */
 524         qr->mm          = mm;
 525     }
 526     else /* ONIOM */
 527     {    /* MM rec creation */
 528         mm              = mk_MMrec();
 529         mm->scalefactor = ir->scalefactor;
 530         mm->nrMMatoms   = 0;
 531         qr->mm          = mm;
 532     }
 533
 534     /* these variables get updated in the update QMMMrec */
 535
 536     if (qr->nrQMlayers == 1)
 537     {
 538         /* with only one layer there is only one initialisation
 539          * needed. Multilayer is a bit more complicated as it requires
 540          * re-initialisation at every step of the simulation. This is due
 541          * to the use of COMMON blocks in the fortran QM subroutines.
 542          */
 543         if (qr->qm[0]->QMmethod < eQMmethodRHF)
 544         {
 545             if (GMX_QMMM_MOPAC)
 546             {
 547                 /* semi-empiprical 1-layer ONIOM calculation requested (mopac93) */
 548                 init_mopac(qr->qm[0]);
 549             }
 550             else
 551             {
 552                 gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac.");
 553             }
 554         }
 555         else
 556         {
 557             /* ab initio calculation requested (gamess/gaussian/ORCA) */
 558             if (GMX_QMMM_GAMESS)
 559             {
 560                 init_gamess(cr, qr->qm[0], qr->mm);
 561             }
 562             else if (GMX_QMMM_GAUSSIAN)
 563             {
 564                 init_gaussian(qr->qm[0]);
 565             }
 566             else if (GMX_QMMM_ORCA)
 567             {
 568                 init_orca(qr->qm[0]);
 569             }
 570             else
 571             {
 572                 gmx_fatal(FARGS,
 573                           "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
 574             }
 575         }
 576     }
 577 } /* init_QMMMrec */
 578
 579 void update_QMMMrec(const t_commrec* cr, const t_forcerec* fr, const rvec* x, const t_mdatoms* md, const matrix box)
 580 {
 581     /* updates the coordinates of both QM atoms and MM atoms and stores
 582      * them in the QMMMrec.
 583      *
 584      * NOTE: is NOT yet working if there are no PBC. Also in ns.c, simple
 585      * ns needs to be fixed!
 586      */
 587     int           mm_max = 0, mm_nr = 0, mm_nr_new, i, j, is, k, shift;
 588     t_j_particle *mm_j_particles = nullptr, *qm_i_particles = nullptr;
 589     t_QMMMrec*    qr;
 590     t_nblist*     QMMMlist;
 591     rvec          dx;
 592     ivec          crd;
 593     t_QMrec*      qm;
 594     t_MMrec*      mm;
 595     t_pbc         pbc;
 596     int*          parallelMMarray = nullptr;
 597
 598     if (!GMX_QMMM)
 599     {
 600         gmx_incons("Compiled without QMMM");
 601     }
 602
 603     /* every cpu has this array. On every processor we fill this array
 604      * with 1's and 0's. 1's indicate the atoms is a QM atom on the
 605      * current cpu in a later stage these arrays are all summed. indexes
 606      * > 0 indicate the atom is a QM atom. Every node therefore knows
 607      * whcih atoms are part of the QM subsystem.
 608      */
 609     /* copy some pointers */
 610     qr       = fr->qr;
 611     mm       = qr->mm;
 612     QMMMlist = fr->QMMMlist;
 613
 614     /*  init_pbc(box);  needs to be called first, see pbc.h */
 615     ivec null_ivec;
 616     clear_ivec(null_ivec);
 617     set_pbc_dd(&pbc, fr->ePBC, DOMAINDECOMP(cr) ? cr->dd->numCells : null_ivec, FALSE, box);
 618     /* only in standard (normal) QMMM we need the neighbouring MM
 619      * particles to provide a electric field of point charges for the QM
 620      * atoms.
 621      */
 622     if (qr->QMMMscheme == eQMMMschemenormal) /* also implies 1 QM-layer */
 623     {
 624         /* we NOW create/update a number of QMMMrec entries:
 625          *
 626          * 1) the shiftQM, containing the shifts of the QM atoms
 627          *
 628          * 2) the indexMM array, containing the index of the MM atoms
 629          *
 630          * 3) the shiftMM, containing the shifts of the MM atoms
 631          *
 632          * 4) the shifted coordinates of the MM atoms
 633          *
 634          * the shifts are used for computing virial of the QM/MM particles.
 635          */
 636         qm = qr->qm[0]; /* in case of normal QMMM, there is only one group */
 637         snew(qm_i_particles, QMMMlist->nri);
 638         if (QMMMlist->nri)
 639         {
 640             qm_i_particles[0].shift = XYZ2IS(0, 0, 0);
 641             for (i = 0; i < QMMMlist->nri; i++)
 642             {
 643                 qm_i_particles[i].j = QMMMlist->iinr[i];
 644
 645                 if (i)
 646                 {
 647                     qm_i_particles[i].shift =
 648                             pbc_dx_aiuc(&pbc, x[QMMMlist->iinr[0]], x[QMMMlist->iinr[i]], dx);
 649                 }
 650                 /* However, since nri >= nrQMatoms, we do a quicksort, and throw
 651                  * out double, triple, etc. entries later, as we do for the MM
 652                  * list too.
 653                  */
 654
 655                 /* compute the shift for the MM j-particles with respect to
 656                  * the QM i-particle and store them.
 657                  */
 658
 659                 crd[0] = IS2X(QMMMlist->shift[i]) + IS2X(qm_i_particles[i].shift);
 660                 crd[1] = IS2Y(QMMMlist->shift[i]) + IS2Y(qm_i_particles[i].shift);
 661                 crd[2] = IS2Z(QMMMlist->shift[i]) + IS2Z(qm_i_particles[i].shift);
 662                 is     = XYZ2IS(crd[0], crd[1], crd[2]);
 663                 for (j = QMMMlist->jindex[i]; j < QMMMlist->jindex[i + 1]; j++)
 664                 {
 665                     if (mm_nr >= mm_max)
 666                     {
 667                         mm_max += 1000;
 668                         srenew(mm_j_particles, mm_max);
 669                     }
 670
 671                     mm_j_particles[mm_nr].j     = QMMMlist->jjnr[j];
 672                     mm_j_particles[mm_nr].shift = is;
 673                     mm_nr++;
 674                 }
 675             }
 676
 677             /* quicksort QM and MM shift arrays and throw away multiple entries */
 678
 679
 680             std::sort(qm_i_particles, qm_i_particles + QMMMlist->nri, struct_comp);
 681             /* The mm_j_particles argument to qsort is not allowed to be nullptr */
 682             if (mm_nr > 0)
 683             {
 684                 std::sort(mm_j_particles, mm_j_particles + mm_nr, struct_comp);
 685             }
 686             /* remove multiples in the QM shift array, since in init_QMMM() we
 687              * went through the atom numbers from 0 to md.nr, the order sorted
 688              * here matches the one of QMindex already.
 689              */
 690             j = 0;
 691             for (i = 0; i < QMMMlist->nri; i++)
 692             {
 693                 if (i == 0 || qm_i_particles[i].j != qm_i_particles[i - 1].j)
 694                 {
 695                     qm_i_particles[j++] = qm_i_particles[i];
 696                 }
 697             }
 698             mm_nr_new = 0;
 699             /* Remove double entries for the MM array.
 700              * Also remove mm atoms that have no charges!
 701              * actually this is already done in the ns.c
 702              */
 703             for (i = 0; i < mm_nr; i++)
 704             {
 705                 if ((i == 0 || mm_j_particles[i].j != mm_j_particles[i - 1].j)
 706                     && !md->bQM[mm_j_particles[i].j]
 707                     && ((md->chargeA[mm_j_particles[i].j] != 0.0_real)
 708                         || (md->chargeB && (md->chargeB[mm_j_particles[i].j] != 0.0_real))))
 709                 {
 710                     mm_j_particles[mm_nr_new++] = mm_j_particles[i];
 711                 }
 712             }
 713             mm_nr = mm_nr_new;
 714             /* store the data retrieved above into the QMMMrec
 715              */
 716             k = 0;
 717             /* Keep the compiler happy,
 718              * shift will always be set in the loop for i=0
 719              */
 720             shift = 0;
 721             for (i = 0; i < qm->nrQMatoms; i++)
 722             {
 723                 /* not all qm particles might have appeared as i
 724                  * particles. They might have been part of the same charge
 725                  * group for instance.
 726                  */
 727                 if (qm->indexQM[i] == qm_i_particles[k].j)
 728                 {
 729                     shift = qm_i_particles[k++].shift;
 730                 }
 731                 /* use previous shift, assuming they belong the same charge
 732                  * group anyway,
 733                  */
 734
 735                 qm->shiftQM[i] = shift;
 736             }
 737         }
 738         /* parallel excecution */
 739         if (PAR(cr))
 740         {
 741             snew(parallelMMarray, 2 * (md->nr));
 742             /* only MM particles have a 1 at their atomnumber. The second part
 743              * of the array contains the shifts. Thus:
 744              * p[i]=1/0 depending on wether atomnumber i is a MM particle in the QM
 745              * step or not. p[i+md->nr] is the shift of atomnumber i.
 746              */
 747             for (i = 0; i < 2 * (md->nr); i++)
 748             {
 749                 parallelMMarray[i] = 0;
 750             }
 751
 752             for (i = 0; i < mm_nr; i++)
 753             {
 754                 parallelMMarray[mm_j_particles[i].j]            = 1;
 755                 parallelMMarray[mm_j_particles[i].j + (md->nr)] = mm_j_particles[i].shift;
 756             }
 757             gmx_sumi(md->nr, parallelMMarray, cr);
 758             mm_nr = 0;
 759
 760             mm_max = 0;
 761             for (i = 0; i < md->nr; i++)
 762             {
 763                 if (parallelMMarray[i])
 764                 {
 765                     if (mm_nr >= mm_max)
 766                     {
 767                         mm_max += 1000;
 768                         srenew(mm->indexMM, mm_max);
 769                         srenew(mm->shiftMM, mm_max);
 770                     }
 771                     mm->indexMM[mm_nr]   = i;
 772                     mm->shiftMM[mm_nr++] = parallelMMarray[i + md->nr] / parallelMMarray[i];
 773                 }
 774             }
 775             mm->nrMMatoms = mm_nr;
 776             free(parallelMMarray);
 777         }
 778         /* serial execution */
 779         else
 780         {
 781             mm->nrMMatoms = mm_nr;
 782             srenew(mm->shiftMM, mm_nr);
 783             srenew(mm->indexMM, mm_nr);
 784             for (i = 0; i < mm_nr; i++)
 785             {
 786                 mm->indexMM[i] = mm_j_particles[i].j;
 787                 mm->shiftMM[i] = mm_j_particles[i].shift;
 788             }
 789         }
 790         /* (re) allocate memory for the MM coordiate array. The QM
 791          * coordinate array was already allocated in init_QMMM, and is
 792          * only (re)filled in the update_QMMM_coordinates routine
 793          */
 794         srenew(mm->xMM, mm->nrMMatoms);
 795         /* now we (re) fill the array that contains the MM charges with
 796          * the forcefield charges. If requested, these charges will be
 797          * scaled by a factor
 798          */
 799         srenew(mm->MMcharges, mm->nrMMatoms);
 800         for (i = 0; i < mm->nrMMatoms; i++) /* no free energy yet */
 801         {
 802             mm->MMcharges[i] = md->chargeA[mm->indexMM[i]] * mm->scalefactor;
 803         }
 804         /* the next routine fills the coordinate fields in the QMMM rec of
 805          * both the qunatum atoms and the MM atoms, using the shifts
 806          * calculated above.
 807          */
 808
 809         update_QMMM_coord(x, fr, qr->qm[0], qr->mm);
 810         free(qm_i_particles);
 811         free(mm_j_particles);
 812     }
 813     else /* ONIOM */ /* ????? */
 814     {
 815         mm->nrMMatoms = 0;
 816         /* do for each layer */
 817         for (j = 0; j < qr->nrQMlayers; j++)
 818         {
 819             qm             = qr->qm[j];
 820             qm->shiftQM[0] = XYZ2IS(0, 0, 0);
 821             for (i = 1; i < qm->nrQMatoms; i++)
 822             {
 823                 qm->shiftQM[i] = pbc_dx_aiuc(&pbc, x[qm->indexQM[0]], x[qm->indexQM[i]], dx);
 824             }
 825             update_QMMM_coord(x, fr, qm, mm);
 826         }
 827     }
 828 } /* update_QMMM_rec */
 829
 830 real calculate_QMMM(const t_commrec* cr, gmx::ForceWithShiftForces* forceWithShiftForces, const t_QMMMrec* qr)
 831 {
 832     real QMener = 0.0;
 833     /* a selection for the QM package depending on which is requested
 834      * (Gaussian, GAMESS-UK, MOPAC or ORCA) needs to be implemented here. Now
 835      * it works through defines.... Not so nice yet
 836      */
 837     t_QMrec *qm, *qm2;
 838     t_MMrec* mm     = nullptr;
 839     rvec *   forces = nullptr, *fshift = nullptr, *forces2 = nullptr,
 840          *fshift2 = nullptr; /* needed for multilayer ONIOM */
 841     int i, j, k;
 842
 843     if (!GMX_QMMM)
 844     {
 845         gmx_incons("Compiled without QMMM");
 846     }
 847
 848     /* make a local copy the QMMMrec pointer
 849      */
 850     mm = qr->mm;
 851
 852     /* now different procedures are carried out for one layer ONION and
 853      * normal QMMM on one hand and multilayer oniom on the other
 854      */
 855     gmx::ArrayRef<gmx::RVec> fMM      = forceWithShiftForces->force();
 856     gmx::ArrayRef<gmx::RVec> fshiftMM = forceWithShiftForces->shiftForces();
 857     if (qr->QMMMscheme == eQMMMschemenormal || qr->nrQMlayers == 1)
 858     {
 859         qm = qr->qm[0];
 860         snew(forces, (qm->nrQMatoms + mm->nrMMatoms));
 861         snew(fshift, (qm->nrQMatoms + mm->nrMMatoms));
 862         QMener = call_QMroutine(cr, qr, qm, mm, forces, fshift);
 863         for (i = 0; i < qm->nrQMatoms; i++)
 864         {
 865             for (j = 0; j < DIM; j++)
 866             {
 867                 fMM[qm->indexQM[i]][j] -= forces[i][j];
 868                 fshiftMM[qm->shiftQM[i]][j] += fshift[i][j];
 869             }
 870         }
 871         for (i = 0; i < mm->nrMMatoms; i++)
 872         {
 873             for (j = 0; j < DIM; j++)
 874             {
 875                 fMM[mm->indexMM[i]][j] -= forces[qm->nrQMatoms + i][j];
 876                 fshiftMM[mm->shiftMM[i]][j] += fshift[qm->nrQMatoms + i][j];
 877             }
 878         }
 879         free(forces);
 880         free(fshift);
 881     }
 882     else /* Multi-layer ONIOM */
 883     {
 884         for (i = 0; i < qr->nrQMlayers - 1; i++) /* last layer is special */
 885         {
 886             qm  = qr->qm[i];
 887             qm2 = copy_QMrec(qr->qm[i + 1]);
 888
 889             qm2->nrQMatoms = qm->nrQMatoms;
 890
 891             for (j = 0; j < qm2->nrQMatoms; j++)
 892             {
 893                 for (k = 0; k < DIM; k++)
 894                 {
 895                     qm2->xQM[j][k] = qm->xQM[j][k];
 896                 }
 897                 qm2->indexQM[j]        = qm->indexQM[j];
 898                 qm2->atomicnumberQM[j] = qm->atomicnumberQM[j];
 899                 qm2->shiftQM[j]        = qm->shiftQM[j];
 900             }
 901
 902             qm2->QMcharge = qm->QMcharge;
 903             /* this layer at the higher level of theory */
 904             srenew(forces, qm->nrQMatoms);
 905             srenew(fshift, qm->nrQMatoms);
 906             /* we need to re-initialize the QMroutine every step... */
 907             init_QMroutine(cr, qm, mm);
 908             QMener += call_QMroutine(cr, qr, qm, mm, forces, fshift);
 909
 910             /* this layer at the lower level of theory */
 911             srenew(forces2, qm->nrQMatoms);
 912             srenew(fshift2, qm->nrQMatoms);
 913             init_QMroutine(cr, qm2, mm);
 914             QMener -= call_QMroutine(cr, qr, qm2, mm, forces2, fshift2);
 915             /* E = E1high-E1low The next layer includes the current layer at
 916              * the lower level of theory, which provides + E2low
 917              * this is similar for gradients
 918              */
 919             for (i = 0; i < qm->nrQMatoms; i++)
 920             {
 921                 for (j = 0; j < DIM; j++)
 922                 {
 923                     fMM[qm->indexQM[i]][j] -= (forces[i][j] - forces2[i][j]);
 924                     fshiftMM[qm->shiftQM[i]][j] += (fshift[i][j] - fshift2[i][j]);
 925                 }
 926             }
 927             free(qm2);
 928         }
 929         /* now the last layer still needs to be done: */
 930         qm = qr->qm[qr->nrQMlayers - 1]; /* C counts from 0 */
 931         init_QMroutine(cr, qm, mm);
 932         srenew(forces, qm->nrQMatoms);
 933         srenew(fshift, qm->nrQMatoms);
 934         QMener += call_QMroutine(cr, qr, qm, mm, forces, fshift);
 935         for (i = 0; i < qm->nrQMatoms; i++)
 936         {
 937             for (j = 0; j < DIM; j++)
 938             {
 939                 fMM[qm->indexQM[i]][j] -= forces[i][j];
 940                 fshiftMM[qm->shiftQM[i]][j] += fshift[i][j];
 941             }
 942         }
 943         free(forces);
 944         free(fshift);
 945         free(forces2);
 946         free(fshift2);
 947     }
 948     return (QMener);
 949 } /* calculate_QMMM */
 950
 951 #pragma GCC diagnostic pop