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