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