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 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.
  11  *
  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.
  16  *
  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.
  21  *
  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.
  26  *
  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.
  34  *
  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.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "qmmm.h"
  41
  42 #include "config.h"
  43
  44 #include <cmath>
  45 #include <cstdio>
  46 #include <cstdlib>
  47 #include <cstring>
  48
  49 #include <algorithm>
  50
  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"
  76
  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"
  82
  83 /* this struct and these comparison functions are needed for creating
  84  * a QMMM input for the QM routines from the QMMM neighbor list.
  85  */
  86
  87 typedef struct
  88 {
  89     int j;
  90     int shift;
  91 } t_j_particle;
  92
  93 static bool struct_comp(const t_j_particle& a, const t_j_particle& b)
  94 {
  95     return a.j < b.j;
  96 }
  97
  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,
 102                            rvec gmx_unused f[],
 103                            rvec gmx_unused fshift[])
 104 {
 105     /* makes a call to the requested QM routine (qm->QMmethod)
 106      * Note that f is actually the gradient, i.e. -f
 107      */
 108     /* do a semi-empiprical calculation */
 109
 110     if (qm->QMmethod < eQMmethodRHF && !(mm->nrMMatoms))
 111     {
 112         if (GMX_QMMM_MOPAC)
 113         {
 114             if (qm->bSH)
 115             {
 116                 return call_mopac_SH(qm, mm, f, fshift);
 117             }
 118             else
 119             {
 120                 return call_mopac(qm, mm, f, fshift);
 121             }
 122         }
 123         else
 124         {
 125             gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac.");
 126         }
 127     }
 128     else
 129     {
 130         /* do an ab-initio calculation */
 131         if (qm->bSH && qm->QMmethod == eQMmethodCASSCF)
 132         {
 133             if (GMX_QMMM_GAUSSIAN)
 134             {
 135                 return call_gaussian_SH(qmmm, qm, mm, f, fshift);
 136             }
 137             else
 138             {
 139                 gmx_fatal(FARGS, "Ab-initio Surface-hopping only supported with Gaussian.");
 140             }
 141         }
 142         else
 143         {
 144             if (GMX_QMMM_GAMESS)
 145             {
 146                 return call_gamess(qm, mm, f, fshift);
 147             }
 148             else if (GMX_QMMM_GAUSSIAN)
 149             {
 150                 return call_gaussian(qmmm, qm, mm, f, fshift);
 151             }
 152             else if (GMX_QMMM_ORCA)
 153             {
 154                 return call_orca(qmmm, qm, mm, f, fshift);
 155             }
 156             else
 157             {
 158                 gmx_fatal(FARGS,
 159                           "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
 160             }
 161         }
 162     }
 163 }
 164
 165 static void init_QMroutine(const t_commrec gmx_unused* cr, t_QMrec gmx_unused* qm, t_MMrec gmx_unused* mm)
 166 {
 167     /* makes a call to the requested QM routine (qm->QMmethod)
 168      */
 169     if (qm->QMmethod < eQMmethodRHF)
 170     {
 171         if (GMX_QMMM_MOPAC)
 172         {
 173             /* do a semi-empiprical calculation */
 174             init_mopac(qm);
 175         }
 176         else
 177         {
 178             gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac.");
 179         }
 180     }
 181     else
 182     {
 183         /* do an ab-initio calculation */
 184         if (GMX_QMMM_GAMESS)
 185         {
 186             init_gamess(cr, qm, mm);
 187         }
 188         else if (GMX_QMMM_GAUSSIAN)
 189         {
 190             init_gaussian(qm);
 191         }
 192         else if (GMX_QMMM_ORCA)
 193         {
 194             init_orca(qm);
 195         }
 196         else
 197         {
 198             gmx_fatal(FARGS, "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
 199         }
 200     }
 201 } /* init_QMroutine */
 202
 203 static void update_QMMM_coord(const rvec* x, const t_forcerec* fr, t_QMrec* qm, t_MMrec* mm)
 204 {
 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.
 208      */
 209     int i;
 210
 211     /* shift the QM atoms into the central box
 212      */
 213     for (i = 0; i < qm->nrQMatoms; i++)
 214     {
 215         rvec_sub(x[qm->indexQM[i]], fr->shift_vec[qm->shiftQM[i]], qm->xQM[i]);
 216     }
 217     /* also shift the MM atoms into the central box, if any
 218      */
 219     for (i = 0; i < mm->nrMMatoms; i++)
 220     {
 221         rvec_sub(x[mm->indexMM[i]], fr->shift_vec[mm->shiftMM[i]], mm->xMM[i]);
 222     }
 223 } /* update_QMMM_coord */
 224
 225 /* end of QMMM subroutines */
 226
 227 /* QMMM core routines */
 228
 229 static t_QMrec* mk_QMrec()
 230 {
 231     t_QMrec* qm;
 232     snew(qm, 1);
 233     return qm;
 234 } /* mk_QMrec */
 235
 236 static t_MMrec* mk_MMrec()
 237 {
 238     t_MMrec* mm;
 239     snew(mm, 1);
 240     return mm;
 241 } /* mk_MMrec */
 242
 243 static void init_QMrec(int grpnr, t_QMrec* qm, int nr, const int* atomarray, const gmx_mtop_t* mtop, const t_inputrec* ir)
 244 {
 245     /* fills the t_QMrec struct of QM group grpnr
 246      */
 247
 248     qm->nrQMatoms = nr;
 249     snew(qm->xQM, nr);
 250     snew(qm->indexQM, nr);
 251     snew(qm->shiftQM, nr); /* the shifts */
 252     for (int i = 0; i < nr; i++)
 253     {
 254         qm->indexQM[i] = atomarray[i];
 255     }
 256
 257     snew(qm->atomicnumberQM, nr);
 258     int molb = 0;
 259     for (int i = 0; i < qm->nrQMatoms; i++)
 260     {
 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];
 264     }
 265
 266     qm->QMcharge     = ir->opts.QMcharge[grpnr];
 267     qm->multiplicity = ir->opts.QMmult[grpnr];
 268     qm->nelectrons -= ir->opts.QMcharge[grpnr];
 269
 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
 282                       * (init_gaussian()
 283                       */
 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]);
 287 } /* init_QMrec */
 288
 289 static t_QMrec* copy_QMrec(t_QMrec* qm)
 290 {
 291     /* copies the contents of qm into a new t_QMrec struct */
 292     t_QMrec* qmcopy;
 293     int      i;
 294
 295     qmcopy            = mk_QMrec();
 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++)
 302     {
 303         qmcopy->shiftQM[i]        = qm->shiftQM[i];
 304         qmcopy->indexQM[i]        = qm->indexQM[i];
 305         qmcopy->atomicnumberQM[i] = qm->atomicnumberQM[i];
 306     }
 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;
 320
 321     /* Gaussian init. variables */
 322     qmcopy->nQMcpus = qm->nQMcpus;
 323     for (i = 0; i < DIM; i++)
 324     {
 325         qmcopy->SHbasis[i] = qm->SHbasis[i];
 326     }
 327     qmcopy->QMmem    = qm->QMmem;
 328     qmcopy->accuracy = qm->accuracy;
 329     qmcopy->cpmcscf  = qm->cpmcscf;
 330     qmcopy->SAstep   = qm->SAstep;
 331
 332     return (qmcopy);
 333
 334 } /*copy_QMrec */
 335
 336 #if GMX_QMMM
 337
 338 t_QMMMrec* mk_QMMMrec()
 339 {
 340     t_QMMMrec* qr;
 341
 342     snew(qr, 1);
 343
 344     return qr;
 345
 346 } /* mk_QMMMrec */
 347
 348 #else /* GMX_QMMM */
 349
 350 t_QMMMrec* mk_QMMMrec()
 351 {
 352     gmx_incons("Compiled without QMMM");
 353 } /* mk_QMMMrec */
 354 #endif
 355
 356 void init_QMMMrec(const t_commrec* cr, const gmx_mtop_t* mtop, const t_inputrec* ir, const t_forcerec* fr)
 357 {
 358     /* we put the atomsnumbers of atoms that belong to the QMMM group in
 359      * an array that will be copied later to QMMMrec->indexQM[..]. Also
 360      * it will be used to create an QMMMrec->bQMMM index array that
 361      * simply contains true/false for QM and MM (the other) atoms.
 362      */
 363
 364     t_QMMMrec* qr;
 365     t_MMrec*   mm;
 366
 367     if (!GMX_QMMM)
 368     {
 369         gmx_incons("Compiled without QMMM");
 370     }
 371
 372     if (ir->cutoff_scheme != ecutsGROUP)
 373     {
 374         gmx_fatal(FARGS, "QMMM is currently only supported with cutoff-scheme=group");
 375     }
 376     if (!EI_DYNAMICS(ir->eI))
 377     {
 378         gmx_fatal(FARGS, "QMMM is only supported with dynamics");
 379     }
 380
 381     /* issue a fatal if the user wants to run with more than one node */
 382     if (PAR(cr))
 383     {
 384         gmx_fatal(FARGS, "QM/MM does not work in parallel, use a single rank instead\n");
 385     }
 386
 387     /* Make a local copy of the QMMMrec */
 388     qr = fr->qr;
 389
 390     /* bQMMM[..] is an array containing TRUE/FALSE for atoms that are
 391      * QM/not QM. We first set all elemenst at false. Afterwards we use
 392      * the qm_arr (=MMrec->indexQM) to changes the elements
 393      * corresponding to the QM atoms at TRUE.  */
 394
 395     qr->QMMMscheme = ir->QMMMscheme;
 396
 397     /* we take the possibility into account that a user has
 398      * defined more than one QM group:
 399      */
 400     /* an ugly work-around in case there is only one group In this case
 401      * the whole system is treated as QM. Otherwise the second group is
 402      * always the rest of the total system and is treated as MM.
 403      */
 404
 405     /* small problem if there is only QM.... so no MM */
 406
 407     int numQmmmGroups = ir->opts.ngQM;
 408
 409     if (qr->QMMMscheme == eQMMMschemeoniom)
 410     {
 411         qr->nrQMlayers = numQmmmGroups;
 412     }
 413     else
 414     {
 415         qr->nrQMlayers = 1;
 416     }
 417
 418     /* there are numQmmmGroups groups of QM atoms. In case of multiple QM groups
 419      * I assume that the users wants to do ONIOM. However, maybe it
 420      * should also be possible to define more than one QM subsystem with
 421      * independent neighbourlists. I have to think about
 422      * that.. 11-11-2003
 423      */
 424     std::vector<int> qmmmAtoms;
 425     snew(qr->qm, numQmmmGroups);
 426     for (int i = 0; i < numQmmmGroups; i++)
 427     {
 428         /* new layer */
 429         if (qr->QMMMscheme == eQMMMschemeoniom)
 430         {
 431             /* add the atoms to the bQMMM array
 432              */
 433
 434             /* I assume that users specify the QM groups from small to
 435              * big(ger) in the mdp file
 436              */
 437             qr->qm[i] = mk_QMrec();
 438             /* store QM atoms in this layer in the QMrec and initialise layer
 439              */
 440             init_QMrec(i, qr->qm[i], qmmmAtoms.size(), qmmmAtoms.data(), mtop, ir);
 441         }
 442     }
 443     if (qr->QMMMscheme != eQMMMschemeoniom)
 444     {
 445
 446         /* standard QMMM, all layers are merged together so there is one QM
 447          * subsystem and one MM subsystem.
 448          * Also we set the charges to zero in mtop to prevent the innerloops
 449          * from doubly counting the electostatic QM MM interaction
 450          * TODO: Consider doing this in grompp instead.
 451          */
 452
 453         qr->qm[0] = mk_QMrec();
 454         /* store QM atoms in the QMrec and initialise
 455          */
 456         init_QMrec(0, qr->qm[0], qmmmAtoms.size(), qmmmAtoms.data(), mtop, ir);
 457
 458         /* MM rec creation */
 459         mm              = mk_MMrec();
 460         mm->scalefactor = ir->scalefactor;
 461         mm->nrMMatoms   = (mtop->natoms) - (qr->qm[0]->nrQMatoms); /* rest of the atoms */
 462         qr->mm          = mm;
 463     }
 464     else /* ONIOM */
 465     {    /* MM rec creation */
 466         mm              = mk_MMrec();
 467         mm->scalefactor = ir->scalefactor;
 468         mm->nrMMatoms   = 0;
 469         qr->mm          = mm;
 470     }
 471
 472     /* these variables get updated in the update QMMMrec */
 473
 474     if (qr->nrQMlayers == 1)
 475     {
 476         /* with only one layer there is only one initialisation
 477          * needed. Multilayer is a bit more complicated as it requires
 478          * re-initialisation at every step of the simulation. This is due
 479          * to the use of COMMON blocks in the fortran QM subroutines.
 480          */
 481         if (qr->qm[0]->QMmethod < eQMmethodRHF)
 482         {
 483             if (GMX_QMMM_MOPAC)
 484             {
 485                 /* semi-empiprical 1-layer ONIOM calculation requested (mopac93) */
 486                 init_mopac(qr->qm[0]);
 487             }
 488             else
 489             {
 490                 gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac.");
 491             }
 492         }
 493         else
 494         {
 495             /* ab initio calculation requested (gamess/gaussian/ORCA) */
 496             if (GMX_QMMM_GAMESS)
 497             {
 498                 init_gamess(cr, qr->qm[0], qr->mm);
 499             }
 500             else if (GMX_QMMM_GAUSSIAN)
 501             {
 502                 init_gaussian(qr->qm[0]);
 503             }
 504             else if (GMX_QMMM_ORCA)
 505             {
 506                 init_orca(qr->qm[0]);
 507             }
 508             else
 509             {
 510                 gmx_fatal(FARGS,
 511                           "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
 512             }
 513         }
 514     }
 515 } /* init_QMMMrec */
 516
 517 void update_QMMMrec(const t_commrec* cr, const t_forcerec* fr, const rvec* x, const t_mdatoms* md, const matrix box)
 518 {
 519     /* updates the coordinates of both QM atoms and MM atoms and stores
 520      * them in the QMMMrec.
 521      *
 522      * NOTE: is NOT yet working if there are no PBC. Also in ns.c, simple
 523      * ns needs to be fixed!
 524      */
 525     int           mm_max = 0, mm_nr = 0, mm_nr_new, i, j, is, k, shift;
 526     t_j_particle *mm_j_particles = nullptr, *qm_i_particles = nullptr;
 527     t_QMMMrec*    qr;
 528     t_nblist*     QMMMlist;
 529     rvec          dx;
 530     ivec          crd;
 531     t_QMrec*      qm;
 532     t_MMrec*      mm;
 533     t_pbc         pbc;
 534     int*          parallelMMarray = nullptr;
 535
 536     if (!GMX_QMMM)
 537     {
 538         gmx_incons("Compiled without QMMM");
 539     }
 540
 541     /* every cpu has this array. On every processor we fill this array
 542      * with 1's and 0's. 1's indicate the atoms is a QM atom on the
 543      * current cpu in a later stage these arrays are all summed. indexes
 544      * > 0 indicate the atom is a QM atom. Every node therefore knows
 545      * whcih atoms are part of the QM subsystem.
 546      */
 547     /* copy some pointers */
 548     qr       = fr->qr;
 549     mm       = qr->mm;
 550     QMMMlist = fr->QMMMlist;
 551
 552     /*  init_pbc(box);  needs to be called first, see pbc.h */
 553     ivec null_ivec;
 554     clear_ivec(null_ivec);
 555     set_pbc_dd(&pbc, fr->pbcType, DOMAINDECOMP(cr) ? cr->dd->numCells : null_ivec, FALSE, box);
 556     /* only in standard (normal) QMMM we need the neighbouring MM
 557      * particles to provide a electric field of point charges for the QM
 558      * atoms.
 559      */
 560     if (qr->QMMMscheme == eQMMMschemenormal) /* also implies 1 QM-layer */
 561     {
 562         /* we NOW create/update a number of QMMMrec entries:
 563          *
 564          * 1) the shiftQM, containing the shifts of the QM atoms
 565          *
 566          * 2) the indexMM array, containing the index of the MM atoms
 567          *
 568          * 3) the shiftMM, containing the shifts of the MM atoms
 569          *
 570          * 4) the shifted coordinates of the MM atoms
 571          *
 572          * the shifts are used for computing virial of the QM/MM particles.
 573          */
 574         qm = qr->qm[0]; /* in case of normal QMMM, there is only one group */
 575         snew(qm_i_particles, QMMMlist->nri);
 576         if (QMMMlist->nri)
 577         {
 578             qm_i_particles[0].shift = XYZ2IS(0, 0, 0);
 579             for (i = 0; i < QMMMlist->nri; i++)
 580             {
 581                 qm_i_particles[i].j = QMMMlist->iinr[i];
 582
 583                 if (i)
 584                 {
 585                     qm_i_particles[i].shift =
 586                             pbc_dx_aiuc(&pbc, x[QMMMlist->iinr[0]], x[QMMMlist->iinr[i]], dx);
 587                 }
 588                 /* However, since nri >= nrQMatoms, we do a quicksort, and throw
 589                  * out double, triple, etc. entries later, as we do for the MM
 590                  * list too.
 591                  */
 592
 593                 /* compute the shift for the MM j-particles with respect to
 594                  * the QM i-particle and store them.
 595                  */
 596
 597                 crd[0] = IS2X(QMMMlist->shift[i]) + IS2X(qm_i_particles[i].shift);
 598                 crd[1] = IS2Y(QMMMlist->shift[i]) + IS2Y(qm_i_particles[i].shift);
 599                 crd[2] = IS2Z(QMMMlist->shift[i]) + IS2Z(qm_i_particles[i].shift);
 600                 is     = XYZ2IS(crd[0], crd[1], crd[2]);
 601                 for (j = QMMMlist->jindex[i]; j < QMMMlist->jindex[i + 1]; j++)
 602                 {
 603                     if (mm_nr >= mm_max)
 604                     {
 605                         mm_max += 1000;
 606                         srenew(mm_j_particles, mm_max);
 607                     }
 608
 609                     mm_j_particles[mm_nr].j     = QMMMlist->jjnr[j];
 610                     mm_j_particles[mm_nr].shift = is;
 611                     mm_nr++;
 612                 }
 613             }
 614
 615             /* quicksort QM and MM shift arrays and throw away multiple entries */
 616
 617
 618             std::sort(qm_i_particles, qm_i_particles + QMMMlist->nri, struct_comp);
 619             /* The mm_j_particles argument to qsort is not allowed to be nullptr */
 620             if (mm_nr > 0)
 621             {
 622                 std::sort(mm_j_particles, mm_j_particles + mm_nr, struct_comp);
 623             }
 624             /* remove multiples in the QM shift array, since in init_QMMM() we
 625              * went through the atom numbers from 0 to md.nr, the order sorted
 626              * here matches the one of QMindex already.
 627              */
 628             j = 0;
 629             for (i = 0; i < QMMMlist->nri; i++)
 630             {
 631                 if (i == 0 || qm_i_particles[i].j != qm_i_particles[i - 1].j)
 632                 {
 633                     qm_i_particles[j++] = qm_i_particles[i];
 634                 }
 635             }
 636             mm_nr_new = 0;
 637             /* Remove double entries for the MM array.
 638              * Also remove mm atoms that have no charges!
 639              * actually this is already done in the ns.c
 640              */
 641             for (i = 0; i < mm_nr; i++)
 642             {
 643                 if ((i == 0 || mm_j_particles[i].j != mm_j_particles[i - 1].j)
 644                     && !md->bQM[mm_j_particles[i].j]
 645                     && ((md->chargeA[mm_j_particles[i].j] != 0.0_real)
 646                         || (md->chargeB && (md->chargeB[mm_j_particles[i].j] != 0.0_real))))
 647                 {
 648                     mm_j_particles[mm_nr_new++] = mm_j_particles[i];
 649                 }
 650             }
 651             mm_nr = mm_nr_new;
 652             /* store the data retrieved above into the QMMMrec
 653              */
 654             k = 0;
 655             /* Keep the compiler happy,
 656              * shift will always be set in the loop for i=0
 657              */
 658             shift = 0;
 659             for (i = 0; i < qm->nrQMatoms; i++)
 660             {
 661                 /* not all qm particles might have appeared as i
 662                  * particles. They might have been part of the same charge
 663                  * group for instance.
 664                  */
 665                 if (qm->indexQM[i] == qm_i_particles[k].j)
 666                 {
 667                     shift = qm_i_particles[k++].shift;
 668                 }
 669                 /* use previous shift, assuming they belong the same charge
 670                  * group anyway,
 671                  */
 672
 673                 qm->shiftQM[i] = shift;
 674             }
 675         }
 676         /* parallel excecution */
 677         if (PAR(cr))
 678         {
 679             snew(parallelMMarray, 2 * (md->nr));
 680             /* only MM particles have a 1 at their atomnumber. The second part
 681              * of the array contains the shifts. Thus:
 682              * p[i]=1/0 depending on wether atomnumber i is a MM particle in the QM
 683              * step or not. p[i+md->nr] is the shift of atomnumber i.
 684              */
 685             for (i = 0; i < 2 * (md->nr); i++)
 686             {
 687                 parallelMMarray[i] = 0;
 688             }
 689
 690             for (i = 0; i < mm_nr; i++)
 691             {
 692                 parallelMMarray[mm_j_particles[i].j]            = 1;
 693                 parallelMMarray[mm_j_particles[i].j + (md->nr)] = mm_j_particles[i].shift;
 694             }
 695             gmx_sumi(md->nr, parallelMMarray, cr);
 696             mm_nr = 0;
 697
 698             mm_max = 0;
 699             for (i = 0; i < md->nr; i++)
 700             {
 701                 if (parallelMMarray[i])
 702                 {
 703                     if (mm_nr >= mm_max)
 704                     {
 705                         mm_max += 1000;
 706                         srenew(mm->indexMM, mm_max);
 707                         srenew(mm->shiftMM, mm_max);
 708                     }
 709                     mm->indexMM[mm_nr]   = i;
 710                     mm->shiftMM[mm_nr++] = parallelMMarray[i + md->nr] / parallelMMarray[i];
 711                 }
 712             }
 713             mm->nrMMatoms = mm_nr;
 714             free(parallelMMarray);
 715         }
 716         /* serial execution */
 717         else
 718         {
 719             mm->nrMMatoms = mm_nr;
 720             srenew(mm->shiftMM, mm_nr);
 721             srenew(mm->indexMM, mm_nr);
 722             for (i = 0; i < mm_nr; i++)
 723             {
 724                 mm->indexMM[i] = mm_j_particles[i].j;
 725                 mm->shiftMM[i] = mm_j_particles[i].shift;
 726             }
 727         }
 728         /* (re) allocate memory for the MM coordiate array. The QM
 729          * coordinate array was already allocated in init_QMMM, and is
 730          * only (re)filled in the update_QMMM_coordinates routine
 731          */
 732         srenew(mm->xMM, mm->nrMMatoms);
 733         /* now we (re) fill the array that contains the MM charges with
 734          * the forcefield charges. If requested, these charges will be
 735          * scaled by a factor
 736          */
 737         srenew(mm->MMcharges, mm->nrMMatoms);
 738         for (i = 0; i < mm->nrMMatoms; i++) /* no free energy yet */
 739         {
 740             mm->MMcharges[i] = md->chargeA[mm->indexMM[i]] * mm->scalefactor;
 741         }
 742         /* the next routine fills the coordinate fields in the QMMM rec of
 743          * both the qunatum atoms and the MM atoms, using the shifts
 744          * calculated above.
 745          */
 746
 747         update_QMMM_coord(x, fr, qr->qm[0], qr->mm);
 748         free(qm_i_particles);
 749         free(mm_j_particles);
 750     }
 751     else /* ONIOM */ /* ????? */
 752     {
 753         mm->nrMMatoms = 0;
 754         /* do for each layer */
 755         for (j = 0; j < qr->nrQMlayers; j++)
 756         {
 757             qm             = qr->qm[j];
 758             qm->shiftQM[0] = XYZ2IS(0, 0, 0);
 759             for (i = 1; i < qm->nrQMatoms; i++)
 760             {
 761                 qm->shiftQM[i] = pbc_dx_aiuc(&pbc, x[qm->indexQM[0]], x[qm->indexQM[i]], dx);
 762             }
 763             update_QMMM_coord(x, fr, qm, mm);
 764         }
 765     }
 766 } /* update_QMMM_rec */
 767
 768 real calculate_QMMM(const t_commrec* cr, gmx::ForceWithShiftForces* forceWithShiftForces, const t_QMMMrec* qr)
 769 {
 770     real QMener = 0.0;
 771     /* a selection for the QM package depending on which is requested
 772      * (Gaussian, GAMESS-UK, MOPAC or ORCA) needs to be implemented here. Now
 773      * it works through defines.... Not so nice yet
 774      */
 775     t_QMrec *qm, *qm2;
 776     t_MMrec* mm     = nullptr;
 777     rvec *   forces = nullptr, *fshift = nullptr, *forces2 = nullptr,
 778          *fshift2 = nullptr; /* needed for multilayer ONIOM */
 779     int i, j, k;
 780
 781     if (!GMX_QMMM)
 782     {
 783         gmx_incons("Compiled without QMMM");
 784     }
 785
 786     /* make a local copy the QMMMrec pointer
 787      */
 788     mm = qr->mm;
 789
 790     /* now different procedures are carried out for one layer ONION and
 791      * normal QMMM on one hand and multilayer oniom on the other
 792      */
 793     gmx::ArrayRef<gmx::RVec> fMM      = forceWithShiftForces->force();
 794     gmx::ArrayRef<gmx::RVec> fshiftMM = forceWithShiftForces->shiftForces();
 795     if (qr->QMMMscheme == eQMMMschemenormal || qr->nrQMlayers == 1)
 796     {
 797         qm = qr->qm[0];
 798         snew(forces, (qm->nrQMatoms + mm->nrMMatoms));
 799         snew(fshift, (qm->nrQMatoms + mm->nrMMatoms));
 800         QMener = call_QMroutine(cr, qr, qm, mm, forces, fshift);
 801         for (i = 0; i < qm->nrQMatoms; i++)
 802         {
 803             for (j = 0; j < DIM; j++)
 804             {
 805                 fMM[qm->indexQM[i]][j] -= forces[i][j];
 806                 fshiftMM[qm->shiftQM[i]][j] += fshift[i][j];
 807             }
 808         }
 809         for (i = 0; i < mm->nrMMatoms; i++)
 810         {
 811             for (j = 0; j < DIM; j++)
 812             {
 813                 fMM[mm->indexMM[i]][j] -= forces[qm->nrQMatoms + i][j];
 814                 fshiftMM[mm->shiftMM[i]][j] += fshift[qm->nrQMatoms + i][j];
 815             }
 816         }
 817         free(forces);
 818         free(fshift);
 819     }
 820     else /* Multi-layer ONIOM */
 821     {
 822         for (i = 0; i < qr->nrQMlayers - 1; i++) /* last layer is special */
 823         {
 824             qm  = qr->qm[i];
 825             qm2 = copy_QMrec(qr->qm[i + 1]);
 826
 827             qm2->nrQMatoms = qm->nrQMatoms;
 828
 829             for (j = 0; j < qm2->nrQMatoms; j++)
 830             {
 831                 for (k = 0; k < DIM; k++)
 832                 {
 833                     qm2->xQM[j][k] = qm->xQM[j][k];
 834                 }
 835                 qm2->indexQM[j]        = qm->indexQM[j];
 836                 qm2->atomicnumberQM[j] = qm->atomicnumberQM[j];
 837                 qm2->shiftQM[j]        = qm->shiftQM[j];
 838             }
 839
 840             qm2->QMcharge = qm->QMcharge;
 841             /* this layer at the higher level of theory */
 842             srenew(forces, qm->nrQMatoms);
 843             srenew(fshift, qm->nrQMatoms);
 844             /* we need to re-initialize the QMroutine every step... */
 845             init_QMroutine(cr, qm, mm);
 846             QMener += call_QMroutine(cr, qr, qm, mm, forces, fshift);
 847
 848             /* this layer at the lower level of theory */
 849             srenew(forces2, qm->nrQMatoms);
 850             srenew(fshift2, qm->nrQMatoms);
 851             init_QMroutine(cr, qm2, mm);
 852             QMener -= call_QMroutine(cr, qr, qm2, mm, forces2, fshift2);
 853             /* E = E1high-E1low The next layer includes the current layer at
 854              * the lower level of theory, which provides + E2low
 855              * this is similar for gradients
 856              */
 857             for (i = 0; i < qm->nrQMatoms; i++)
 858             {
 859                 for (j = 0; j < DIM; j++)
 860                 {
 861                     fMM[qm->indexQM[i]][j] -= (forces[i][j] - forces2[i][j]);
 862                     fshiftMM[qm->shiftQM[i]][j] += (fshift[i][j] - fshift2[i][j]);
 863                 }
 864             }
 865             free(qm2);
 866         }
 867         /* now the last layer still needs to be done: */
 868         qm = qr->qm[qr->nrQMlayers - 1]; /* C counts from 0 */
 869         init_QMroutine(cr, qm, mm);
 870         srenew(forces, qm->nrQMatoms);
 871         srenew(fshift, qm->nrQMatoms);
 872         QMener += call_QMroutine(cr, qr, qm, mm, forces, fshift);
 873         for (i = 0; i < qm->nrQMatoms; i++)
 874         {
 875             for (j = 0; j < DIM; j++)
 876             {
 877                 fMM[qm->indexQM[i]][j] -= forces[i][j];
 878                 fshiftMM[qm->shiftQM[i]][j] += fshift[i][j];
 879             }
 880         }
 881         free(forces);
 882         free(fshift);
 883         free(forces2);
 884         free(fshift2);
 885     }
 886     return (QMener);
 887 } /* calculate_QMMM */
 888
 889 #pragma GCC diagnostic pop