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