2 * This file is part of the GROMACS molecular simulation package.
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, 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.
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.
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.
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.
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.
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.
52 #include "gromacs/domdec/domdec_struct.h"
53 #include "gromacs/fileio/confio.h"
54 #include "gromacs/gmxlib/network.h"
55 #include "gromacs/gmxlib/nrnb.h"
56 #include "gromacs/math/functions.h"
57 #include "gromacs/math/units.h"
58 #include "gromacs/math/vec.h"
59 #include "gromacs/mdlib/force.h"
60 #include "gromacs/mdlib/ns.h"
61 #include "gromacs/mdtypes/commrec.h"
62 #include "gromacs/mdtypes/inputrec.h"
63 #include "gromacs/mdtypes/md_enums.h"
64 #include "gromacs/mdtypes/mdatom.h"
65 #include "gromacs/mdtypes/nblist.h"
66 #include "gromacs/pbcutil/ishift.h"
67 #include "gromacs/pbcutil/pbc.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"
73 /* declarations of the interfaces to the QM packages. The _SH indicate
74 * the QM interfaces can be used for Surface Hopping simulations
77 /* GAMESS interface */
80 init_gamess(t_commrec *cr, t_QMrec *qm, t_MMrec *mm);
83 call_gamess(t_forcerec *fr,
84 t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[]);
90 init_mopac(t_QMrec *qm);
93 call_mopac(t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[]);
96 call_mopac_SH(t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[]);
98 #elif GMX_QMMM_GAUSSIAN
99 /* GAUSSIAN interface */
102 init_gaussian(t_QMrec *qm);
105 call_gaussian_SH(t_forcerec *fr, t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[]);
108 call_gaussian(t_forcerec *fr, t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[]);
114 init_orca(t_QMrec *qm);
117 call_orca(t_forcerec *fr, t_QMrec *qm,
118 t_MMrec *mm, rvec f[], rvec fshift[]);
125 /* this struct and these comparison functions are needed for creating
126 * a QMMM input for the QM routines from the QMMM neighbor list.
134 static int struct_comp(const void *a, const void *b)
137 return (int)(((t_j_particle *)a)->j)-(int)(((t_j_particle *)b)->j);
141 real call_QMroutine(t_commrec gmx_unused *cr, t_forcerec gmx_unused *fr, t_QMrec gmx_unused *qm,
142 t_MMrec gmx_unused *mm, rvec gmx_unused f[], rvec gmx_unused fshift[])
144 /* makes a call to the requested QM routine (qm->QMmethod)
145 * Note that f is actually the gradient, i.e. -f
150 /* do a semi-empiprical calculation */
152 if (qm->QMmethod < eQMmethodRHF && !(mm->nrMMatoms))
157 QMener = call_mopac_SH(qm, mm, f, fshift);
161 QMener = call_mopac(qm, mm, f, fshift);
164 gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac.");
169 /* do an ab-initio calculation */
170 if (qm->bSH && qm->QMmethod == eQMmethodCASSCF)
172 #if GMX_QMMM_GAUSSIAN
173 QMener = call_gaussian_SH(fr, qm, mm, f, fshift);
175 gmx_fatal(FARGS, "Ab-initio Surface-hopping only supported with Gaussian.");
181 QMener = call_gamess(fr, qm, mm, f, fshift);
182 #elif GMX_QMMM_GAUSSIAN
183 QMener = call_gaussian(fr, qm, mm, f, fshift);
185 QMener = call_orca(fr, qm, mm, f, fshift);
187 gmx_fatal(FARGS, "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
194 void init_QMroutine(t_commrec gmx_unused *cr, t_QMrec gmx_unused *qm, t_MMrec gmx_unused *mm)
196 /* makes a call to the requested QM routine (qm->QMmethod)
198 if (qm->QMmethod < eQMmethodRHF)
201 /* do a semi-empiprical calculation */
204 gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac.");
209 /* do an ab-initio calculation */
211 init_gamess(cr, qm, mm);
212 #elif GMX_QMMM_GAUSSIAN
217 gmx_fatal(FARGS, "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
220 } /* init_QMroutine */
222 void update_QMMM_coord(rvec x[], t_forcerec *fr, t_QMrec *qm, t_MMrec *mm)
224 /* shifts the QM and MM particles into the central box and stores
225 * these shifted coordinates in the coordinate arrays of the
226 * QMMMrec. These coordinates are passed on the QM subroutines.
231 /* shift the QM atoms into the central box
233 for (i = 0; i < qm->nrQMatoms; i++)
235 rvec_sub(x[qm->indexQM[i]], fr->shift_vec[qm->shiftQM[i]], qm->xQM[i]);
237 /* also shift the MM atoms into the central box, if any
239 for (i = 0; i < mm->nrMMatoms; i++)
241 rvec_sub(x[mm->indexMM[i]], fr->shift_vec[mm->shiftMM[i]], mm->xMM[i]);
243 } /* update_QMMM_coord */
245 static void punch_QMMM_excl(t_QMrec *qm, t_MMrec *mm, t_blocka *excls)
247 /* punch a file containing the bonded interactions of each QM
248 * atom with MM atoms. These need to be excluded in the QM routines
249 * Only needed in case of QM/MM optimizations
254 i, j, k, nrexcl = 0, *excluded = NULL, max_excl = 0;
257 out = fopen("QMMMexcl.dat", "w");
259 /* this can be done more efficiently I think
261 for (i = 0; i < qm->nrQMatoms; i++)
264 for (j = excls->index[qm->indexQM[i]];
265 j < excls->index[qm->indexQM[i]+1];
268 for (k = 0; k < mm->nrMMatoms; k++)
270 if (mm->indexMM[k] == excls->a[j]) /* the excluded MM atom */
272 if (nrexcl >= max_excl)
275 srenew(excluded, max_excl);
277 excluded[nrexcl++] = k;
283 fprintf(out, "%5d %5d\n", i+1, nrexcl);
284 for (j = 0; j < nrexcl; j++)
286 fprintf(out, "%5d ", excluded[j]);
292 } /* punch_QMMM_excl */
295 /* end of QMMM subroutines */
297 /* QMMM core routines */
299 t_QMrec *mk_QMrec(void)
306 t_MMrec *mk_MMrec(void)
313 static void init_QMrec(int grpnr, t_QMrec *qm, int nr, int *atomarray,
314 gmx_mtop_t *mtop, t_inputrec *ir)
316 /* fills the t_QMrec struct of QM group grpnr
319 gmx_mtop_atomlookup_t alook;
325 snew(qm->indexQM, nr);
326 snew(qm->shiftQM, nr); /* the shifts */
327 for (i = 0; i < nr; i++)
329 qm->indexQM[i] = atomarray[i];
332 alook = gmx_mtop_atomlookup_init(mtop);
334 snew(qm->atomicnumberQM, nr);
335 for (i = 0; i < qm->nrQMatoms; i++)
337 gmx_mtop_atomnr_to_atom(alook, qm->indexQM[i], &atom);
338 qm->nelectrons += mtop->atomtypes.atomnumber[atom->type];
339 qm->atomicnumberQM[i] = mtop->atomtypes.atomnumber[atom->type];
342 gmx_mtop_atomlookup_destroy(alook);
344 qm->QMcharge = ir->opts.QMcharge[grpnr];
345 qm->multiplicity = ir->opts.QMmult[grpnr];
346 qm->nelectrons -= ir->opts.QMcharge[grpnr];
348 qm->QMmethod = ir->opts.QMmethod[grpnr];
349 qm->QMbasis = ir->opts.QMbasis[grpnr];
350 /* trajectory surface hopping setup (Gaussian only) */
351 qm->bSH = ir->opts.bSH[grpnr];
352 qm->CASorbitals = ir->opts.CASorbitals[grpnr];
353 qm->CASelectrons = ir->opts.CASelectrons[grpnr];
354 qm->SAsteps = ir->opts.SAsteps[grpnr];
355 qm->SAon = ir->opts.SAon[grpnr];
356 qm->SAoff = ir->opts.SAoff[grpnr];
357 /* hack to prevent gaussian from reinitializing all the time */
358 qm->nQMcpus = 0; /* number of CPU's to be used by g01, is set
359 * upon initializing gaussian
362 /* print the current layer to allow users to check their input */
363 fprintf(stderr, "Layer %d\nnr of QM atoms %d\n", grpnr, nr);
364 fprintf(stderr, "QMlevel: %s/%s\n\n",
365 eQMmethod_names[qm->QMmethod], eQMbasis_names[qm->QMbasis]);
368 snew(qm->frontatoms, nr);
369 /* Lennard-Jones coefficients */
372 /* do we optimize the QM separately using the algorithms of the QM program??
374 qm->bTS = ir->opts.bTS[grpnr];
375 qm->bOPT = ir->opts.bOPT[grpnr];
379 t_QMrec *copy_QMrec(t_QMrec *qm)
381 /* copies the contents of qm into a new t_QMrec struct */
388 qmcopy->nrQMatoms = qm->nrQMatoms;
389 snew(qmcopy->xQM, qmcopy->nrQMatoms);
390 snew(qmcopy->indexQM, qmcopy->nrQMatoms);
391 snew(qmcopy->atomicnumberQM, qm->nrQMatoms);
392 snew(qmcopy->shiftQM, qmcopy->nrQMatoms); /* the shifts */
393 for (i = 0; i < qmcopy->nrQMatoms; i++)
395 qmcopy->shiftQM[i] = qm->shiftQM[i];
396 qmcopy->indexQM[i] = qm->indexQM[i];
397 qmcopy->atomicnumberQM[i] = qm->atomicnumberQM[i];
399 qmcopy->nelectrons = qm->nelectrons;
400 qmcopy->multiplicity = qm->multiplicity;
401 qmcopy->QMcharge = qm->QMcharge;
402 qmcopy->nelectrons = qm->nelectrons;
403 qmcopy->QMmethod = qm->QMmethod;
404 qmcopy->QMbasis = qm->QMbasis;
405 /* trajectory surface hopping setup (Gaussian only) */
406 qmcopy->bSH = qm->bSH;
407 qmcopy->CASorbitals = qm->CASorbitals;
408 qmcopy->CASelectrons = qm->CASelectrons;
409 qmcopy->SAsteps = qm->SAsteps;
410 qmcopy->SAon = qm->SAon;
411 qmcopy->SAoff = qm->SAoff;
412 qmcopy->bOPT = qm->bOPT;
414 /* Gaussian init. variables */
415 qmcopy->nQMcpus = qm->nQMcpus;
416 for (i = 0; i < DIM; i++)
418 qmcopy->SHbasis[i] = qm->SHbasis[i];
420 qmcopy->QMmem = qm->QMmem;
421 qmcopy->accuracy = qm->accuracy;
422 qmcopy->cpmcscf = qm->cpmcscf;
423 qmcopy->SAstep = qm->SAstep;
424 snew(qmcopy->frontatoms, qm->nrQMatoms);
425 snew(qmcopy->c12, qmcopy->nrQMatoms);
426 snew(qmcopy->c6, qmcopy->nrQMatoms);
427 if (qmcopy->bTS || qmcopy->bOPT)
429 for (i = 1; i < qmcopy->nrQMatoms; i++)
431 qmcopy->frontatoms[i] = qm->frontatoms[i];
432 qmcopy->c12[i] = qm->c12[i];
433 qmcopy->c6[i] = qm->c6[i];
441 t_QMMMrec *mk_QMMMrec(void)
452 void init_QMMMrec(t_commrec *cr,
457 /* we put the atomsnumbers of atoms that belong to the QMMM group in
458 * an array that will be copied later to QMMMrec->indexQM[..]. Also
459 * it will be used to create an QMMMrec->bQMMM index array that
460 * simply contains true/false for QM and MM (the other) atoms.
463 gmx_groups_t *groups;
464 int *qm_arr = NULL, vsite, ai, aj;
465 int qm_max = 0, qm_nr = 0, i, j, jmax, k, l, nrvsite2 = 0;
470 gmx_mtop_atomloop_all_t aloop;
472 gmx_mtop_ilistloop_all_t iloop;
475 gmx_mtop_atomlookup_t alook;
477 if (ir->cutoff_scheme != ecutsGROUP)
479 gmx_fatal(FARGS, "QMMM is currently only supported with cutoff-scheme=group");
481 if (!EI_DYNAMICS(ir->eI))
483 gmx_fatal(FARGS, "QMMM is only supported with dynamics");
486 c6au = (HARTREE2KJ*AVOGADRO*gmx::power6(BOHR2NM));
487 c12au = (HARTREE2KJ*AVOGADRO*gmx::power12(BOHR2NM));
488 /* issue a fatal if the user wants to run with more than one node */
491 gmx_fatal(FARGS, "QM/MM does not work in parallel, use a single rank instead\n");
494 /* Make a local copy of the QMMMrec */
497 /* bQMMM[..] is an array containing TRUE/FALSE for atoms that are
498 * QM/not QM. We first set all elemenst at false. Afterwards we use
499 * the qm_arr (=MMrec->indexQM) to changes the elements
500 * corresponding to the QM atoms at TRUE. */
502 qr->QMMMscheme = ir->QMMMscheme;
504 /* we take the possibility into account that a user has
505 * defined more than one QM group:
507 /* an ugly work-around in case there is only one group In this case
508 * the whole system is treated as QM. Otherwise the second group is
509 * always the rest of the total system and is treated as MM.
512 /* small problem if there is only QM.... so no MM */
514 jmax = ir->opts.ngQM;
516 if (qr->QMMMscheme == eQMMMschemeoniom)
518 qr->nrQMlayers = jmax;
525 groups = &mtop->groups;
527 /* there are jmax groups of QM atoms. In case of multiple QM groups
528 * I assume that the users wants to do ONIOM. However, maybe it
529 * should also be possible to define more than one QM subsystem with
530 * independent neighbourlists. I have to think about
534 for (j = 0; j < jmax; j++)
537 aloop = gmx_mtop_atomloop_all_init(mtop);
538 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
543 srenew(qm_arr, qm_max);
545 if (ggrpnr(groups, egcQMMM, i) == j)
551 if (qr->QMMMscheme == eQMMMschemeoniom)
553 /* add the atoms to the bQMMM array
556 /* I assume that users specify the QM groups from small to
557 * big(ger) in the mdp file
559 qr->qm[j] = mk_QMrec();
560 /* we need to throw out link atoms that in the previous layer
561 * existed to separate this QMlayer from the previous
562 * QMlayer. We use the iatoms array in the idef for that
563 * purpose. If all atoms defining the current Link Atom (Dummy2)
564 * are part of the current QM layer it needs to be removed from
567 iloop = gmx_mtop_ilistloop_all_init(mtop);
568 while (gmx_mtop_ilistloop_all_next(iloop, &ilist_mol, &a_offset))
570 nrvsite2 = ilist_mol[F_VSITE2].nr;
571 iatoms = ilist_mol[F_VSITE2].iatoms;
573 for (k = 0; k < nrvsite2; k += 4)
575 vsite = a_offset + iatoms[k+1]; /* the vsite */
576 ai = a_offset + iatoms[k+2]; /* constructing atom */
577 aj = a_offset + iatoms[k+3]; /* constructing atom */
578 if (ggrpnr(groups, egcQMMM, vsite) == ggrpnr(groups, egcQMMM, ai)
580 ggrpnr(groups, egcQMMM, vsite) == ggrpnr(groups, egcQMMM, aj))
582 /* this dummy link atom needs to be removed from the qm_arr
583 * before making the QMrec of this layer!
585 for (i = 0; i < qm_nr; i++)
587 if (qm_arr[i] == vsite)
589 /* drop the element */
590 for (l = i; l < qm_nr; l++)
592 qm_arr[l] = qm_arr[l+1];
601 /* store QM atoms in this layer in the QMrec and initialise layer
603 init_QMrec(j, qr->qm[j], qm_nr, qm_arr, mtop, ir);
605 /* we now store the LJ C6 and C12 parameters in QM rec in case
606 * we need to do an optimization
608 if (qr->qm[j]->bOPT || qr->qm[j]->bTS)
610 for (i = 0; i < qm_nr; i++)
612 /* nbfp now includes the 6.0/12.0 derivative prefactors */
613 qr->qm[j]->c6[i] = C6(fr->nbfp, mtop->ffparams.atnr, atom->type, atom->type)/c6au/6.0;
614 qr->qm[j]->c12[i] = C12(fr->nbfp, mtop->ffparams.atnr, atom->type, atom->type)/c12au/12.0;
617 /* now we check for frontier QM atoms. These occur in pairs that
618 * construct the vsite
620 iloop = gmx_mtop_ilistloop_all_init(mtop);
621 while (gmx_mtop_ilistloop_all_next(iloop, &ilist_mol, &a_offset))
623 nrvsite2 = ilist_mol[F_VSITE2].nr;
624 iatoms = ilist_mol[F_VSITE2].iatoms;
626 for (k = 0; k < nrvsite2; k += 4)
628 vsite = a_offset + iatoms[k+1]; /* the vsite */
629 ai = a_offset + iatoms[k+2]; /* constructing atom */
630 aj = a_offset + iatoms[k+3]; /* constructing atom */
631 if (ggrpnr(groups, egcQMMM, ai) < (groups->grps[egcQMMM].nr-1) &&
632 (ggrpnr(groups, egcQMMM, aj) >= (groups->grps[egcQMMM].nr-1)))
634 /* mark ai as frontier atom */
635 for (i = 0; i < qm_nr; i++)
637 if ( (qm_arr[i] == ai) || (qm_arr[i] == vsite) )
639 qr->qm[j]->frontatoms[i] = TRUE;
643 else if (ggrpnr(groups, egcQMMM, aj) < (groups->grps[egcQMMM].nr-1) &&
644 (ggrpnr(groups, egcQMMM, ai) >= (groups->grps[egcQMMM].nr-1)))
646 /* mark aj as frontier atom */
647 for (i = 0; i < qm_nr; i++)
649 if ( (qm_arr[i] == aj) || (qm_arr[i] == vsite))
651 qr->qm[j]->frontatoms[i] = TRUE;
659 if (qr->QMMMscheme != eQMMMschemeoniom)
662 /* standard QMMM, all layers are merged together so there is one QM
663 * subsystem and one MM subsystem.
664 * Also we set the charges to zero in the md->charge arrays to prevent
665 * the innerloops from doubly counting the electostatic QM MM interaction
668 alook = gmx_mtop_atomlookup_init(mtop);
670 for (k = 0; k < qm_nr; k++)
672 gmx_mtop_atomnr_to_atom(alook, qm_arr[k], &atom);
676 qr->qm[0] = mk_QMrec();
677 /* store QM atoms in the QMrec and initialise
679 init_QMrec(0, qr->qm[0], qm_nr, qm_arr, mtop, ir);
680 if (qr->qm[0]->bOPT || qr->qm[0]->bTS)
682 for (i = 0; i < qm_nr; i++)
684 gmx_mtop_atomnr_to_atom(alook, qm_arr[i], &atom);
685 /* nbfp now includes the 6.0/12.0 derivative prefactors */
686 qr->qm[0]->c6[i] = C6(fr->nbfp, mtop->ffparams.atnr, atom->type, atom->type)/c6au/6.0;
687 qr->qm[0]->c12[i] = C12(fr->nbfp, mtop->ffparams.atnr, atom->type, atom->type)/c12au/12.0;
691 /* find frontier atoms and mark them true in the frontieratoms array.
693 for (i = 0; i < qm_nr; i++)
695 gmx_mtop_atomnr_to_ilist(alook, qm_arr[i], &ilist_mol, &a_offset);
696 nrvsite2 = ilist_mol[F_VSITE2].nr;
697 iatoms = ilist_mol[F_VSITE2].iatoms;
699 for (k = 0; k < nrvsite2; k += 4)
701 vsite = a_offset + iatoms[k+1]; /* the vsite */
702 ai = a_offset + iatoms[k+2]; /* constructing atom */
703 aj = a_offset + iatoms[k+3]; /* constructing atom */
704 if (ggrpnr(groups, egcQMMM, ai) < (groups->grps[egcQMMM].nr-1) &&
705 (ggrpnr(groups, egcQMMM, aj) >= (groups->grps[egcQMMM].nr-1)))
707 /* mark ai as frontier atom */
708 if ( (qm_arr[i] == ai) || (qm_arr[i] == vsite) )
710 qr->qm[0]->frontatoms[i] = TRUE;
713 else if (ggrpnr(groups, egcQMMM, aj) < (groups->grps[egcQMMM].nr-1) &&
714 (ggrpnr(groups, egcQMMM, ai) >= (groups->grps[egcQMMM].nr-1)))
716 /* mark aj as frontier atom */
717 if ( (qm_arr[i] == aj) || (qm_arr[i] == vsite) )
719 qr->qm[0]->frontatoms[i] = TRUE;
725 gmx_mtop_atomlookup_destroy(alook);
727 /* MM rec creation */
729 mm->scalefactor = ir->scalefactor;
730 mm->nrMMatoms = (mtop->natoms)-(qr->qm[0]->nrQMatoms); /* rest of the atoms */
734 { /* MM rec creation */
736 mm->scalefactor = ir->scalefactor;
741 /* these variables get updated in the update QMMMrec */
743 if (qr->nrQMlayers == 1)
745 /* with only one layer there is only one initialisation
746 * needed. Multilayer is a bit more complicated as it requires
747 * re-initialisation at every step of the simulation. This is due
748 * to the use of COMMON blocks in the fortran QM subroutines.
750 if (qr->qm[0]->QMmethod < eQMmethodRHF)
753 /* semi-empiprical 1-layer ONIOM calculation requested (mopac93) */
754 init_mopac(qr->qm[0]);
756 gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac.");
761 /* ab initio calculation requested (gamess/gaussian/ORCA) */
763 init_gamess(cr, qr->qm[0], qr->mm);
764 #elif GMX_QMMM_GAUSSIAN
765 init_gaussian(qr->qm[0]);
767 init_orca(qr->qm[0]);
769 gmx_fatal(FARGS, "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
775 void update_QMMMrec(t_commrec *cr,
782 /* updates the coordinates of both QM atoms and MM atoms and stores
783 * them in the QMMMrec.
785 * NOTE: is NOT yet working if there are no PBC. Also in ns.c, simple
786 * ns needs to be fixed!
789 mm_max = 0, mm_nr = 0, mm_nr_new, i, j, is, k, shift;
791 *mm_j_particles = NULL, *qm_i_particles = NULL;
805 *parallelMMarray = NULL;
809 c6au = (HARTREE2KJ*AVOGADRO*gmx::power6(BOHR2NM));
810 c12au = (HARTREE2KJ*AVOGADRO*gmx::power12(BOHR2NM));
812 /* every cpu has this array. On every processor we fill this array
813 * with 1's and 0's. 1's indicate the atoms is a QM atom on the
814 * current cpu in a later stage these arrays are all summed. indexes
815 * > 0 indicate the atom is a QM atom. Every node therefore knows
816 * whcih atoms are part of the QM subsystem.
818 /* copy some pointers */
821 QMMMlist = fr->QMMMlist;
825 /* init_pbc(box); needs to be called first, see pbc.h */
827 clear_ivec(null_ivec);
828 set_pbc_dd(&pbc, fr->ePBC, DOMAINDECOMP(cr) ? cr->dd->nc : null_ivec,
830 /* only in standard (normal) QMMM we need the neighbouring MM
831 * particles to provide a electric field of point charges for the QM
834 if (qr->QMMMscheme == eQMMMschemenormal) /* also implies 1 QM-layer */
836 /* we NOW create/update a number of QMMMrec entries:
838 * 1) the shiftQM, containing the shifts of the QM atoms
840 * 2) the indexMM array, containing the index of the MM atoms
842 * 3) the shiftMM, containing the shifts of the MM atoms
844 * 4) the shifted coordinates of the MM atoms
846 * the shifts are used for computing virial of the QM/MM particles.
848 qm = qr->qm[0]; /* in case of normal QMMM, there is only one group */
849 snew(qm_i_particles, QMMMlist->nri);
852 qm_i_particles[0].shift = XYZ2IS(0, 0, 0);
853 for (i = 0; i < QMMMlist->nri; i++)
855 qm_i_particles[i].j = QMMMlist->iinr[i];
859 qm_i_particles[i].shift = pbc_dx_aiuc(&pbc, x[QMMMlist->iinr[0]],
860 x[QMMMlist->iinr[i]], dx);
863 /* However, since nri >= nrQMatoms, we do a quicksort, and throw
864 * out double, triple, etc. entries later, as we do for the MM
868 /* compute the shift for the MM j-particles with respect to
869 * the QM i-particle and store them.
872 crd[0] = IS2X(QMMMlist->shift[i]) + IS2X(qm_i_particles[i].shift);
873 crd[1] = IS2Y(QMMMlist->shift[i]) + IS2Y(qm_i_particles[i].shift);
874 crd[2] = IS2Z(QMMMlist->shift[i]) + IS2Z(qm_i_particles[i].shift);
875 is = static_cast<int>(XYZ2IS(crd[0], crd[1], crd[2]));
876 for (j = QMMMlist->jindex[i];
877 j < QMMMlist->jindex[i+1];
883 srenew(mm_j_particles, mm_max);
886 mm_j_particles[mm_nr].j = QMMMlist->jjnr[j];
887 mm_j_particles[mm_nr].shift = is;
892 /* quicksort QM and MM shift arrays and throw away multiple entries */
896 qsort(qm_i_particles, QMMMlist->nri,
897 (size_t)sizeof(qm_i_particles[0]),
899 /* The mm_j_particles argument to qsort is not allowed to be NULL */
902 qsort(mm_j_particles, mm_nr,
903 (size_t)sizeof(mm_j_particles[0]),
906 /* remove multiples in the QM shift array, since in init_QMMM() we
907 * went through the atom numbers from 0 to md.nr, the order sorted
908 * here matches the one of QMindex already.
911 for (i = 0; i < QMMMlist->nri; i++)
913 if (i == 0 || qm_i_particles[i].j != qm_i_particles[i-1].j)
915 qm_i_particles[j++] = qm_i_particles[i];
919 if (qm->bTS || qm->bOPT)
921 /* only remove double entries for the MM array */
922 for (i = 0; i < mm_nr; i++)
924 if ((i == 0 || mm_j_particles[i].j != mm_j_particles[i-1].j)
925 && !md->bQM[mm_j_particles[i].j])
927 mm_j_particles[mm_nr_new++] = mm_j_particles[i];
931 /* we also remove mm atoms that have no charges!
932 * actually this is already done in the ns.c
936 for (i = 0; i < mm_nr; i++)
938 if ((i == 0 || mm_j_particles[i].j != mm_j_particles[i-1].j)
939 && !md->bQM[mm_j_particles[i].j]
940 && (md->chargeA[mm_j_particles[i].j]
941 || (md->chargeB && md->chargeB[mm_j_particles[i].j])))
943 mm_j_particles[mm_nr_new++] = mm_j_particles[i];
948 /* store the data retrieved above into the QMMMrec
951 /* Keep the compiler happy,
952 * shift will always be set in the loop for i=0
955 for (i = 0; i < qm->nrQMatoms; i++)
957 /* not all qm particles might have appeared as i
958 * particles. They might have been part of the same charge
959 * group for instance.
961 if (qm->indexQM[i] == qm_i_particles[k].j)
963 shift = qm_i_particles[k++].shift;
965 /* use previous shift, assuming they belong the same charge
969 qm->shiftQM[i] = shift;
972 /* parallel excecution */
975 snew(parallelMMarray, 2*(md->nr));
976 /* only MM particles have a 1 at their atomnumber. The second part
977 * of the array contains the shifts. Thus:
978 * p[i]=1/0 depending on wether atomnumber i is a MM particle in the QM
979 * step or not. p[i+md->nr] is the shift of atomnumber i.
981 for (i = 0; i < 2*(md->nr); i++)
983 parallelMMarray[i] = 0;
986 for (i = 0; i < mm_nr; i++)
988 parallelMMarray[mm_j_particles[i].j] = 1;
989 parallelMMarray[mm_j_particles[i].j+(md->nr)] = mm_j_particles[i].shift;
991 gmx_sumi(md->nr, parallelMMarray, cr);
995 for (i = 0; i < md->nr; i++)
997 if (parallelMMarray[i])
1002 srenew(mm->indexMM, mm_max);
1003 srenew(mm->shiftMM, mm_max);
1005 mm->indexMM[mm_nr] = i;
1006 mm->shiftMM[mm_nr++] = parallelMMarray[i+md->nr]/parallelMMarray[i];
1009 mm->nrMMatoms = mm_nr;
1010 free(parallelMMarray);
1012 /* serial execution */
1015 mm->nrMMatoms = mm_nr;
1016 srenew(mm->shiftMM, mm_nr);
1017 srenew(mm->indexMM, mm_nr);
1018 for (i = 0; i < mm_nr; i++)
1020 mm->indexMM[i] = mm_j_particles[i].j;
1021 mm->shiftMM[i] = mm_j_particles[i].shift;
1025 /* (re) allocate memory for the MM coordiate array. The QM
1026 * coordinate array was already allocated in init_QMMM, and is
1027 * only (re)filled in the update_QMMM_coordinates routine
1029 srenew(mm->xMM, mm->nrMMatoms);
1030 /* now we (re) fill the array that contains the MM charges with
1031 * the forcefield charges. If requested, these charges will be
1032 * scaled by a factor
1034 srenew(mm->MMcharges, mm->nrMMatoms);
1035 for (i = 0; i < mm->nrMMatoms; i++) /* no free energy yet */
1037 mm->MMcharges[i] = md->chargeA[mm->indexMM[i]]*mm->scalefactor;
1039 if (qm->bTS || qm->bOPT)
1041 /* store (copy) the c6 and c12 parameters into the MMrec struct
1043 srenew(mm->c6, mm->nrMMatoms);
1044 srenew(mm->c12, mm->nrMMatoms);
1045 for (i = 0; i < mm->nrMMatoms; i++)
1047 /* nbfp now includes the 6.0/12.0 derivative prefactors */
1048 mm->c6[i] = C6(fr->nbfp, top->idef.atnr, md->typeA[mm->indexMM[i]], md->typeA[mm->indexMM[i]])/c6au/6.0;
1049 mm->c12[i] = C12(fr->nbfp, top->idef.atnr, md->typeA[mm->indexMM[i]], md->typeA[mm->indexMM[i]])/c12au/12.0;
1051 punch_QMMM_excl(qr->qm[0], mm, &(top->excls));
1053 /* the next routine fills the coordinate fields in the QMMM rec of
1054 * both the qunatum atoms and the MM atoms, using the shifts
1058 update_QMMM_coord(x, fr, qr->qm[0], qr->mm);
1059 free(qm_i_particles);
1060 free(mm_j_particles);
1062 else /* ONIOM */ /* ????? */
1065 /* do for each layer */
1066 for (j = 0; j < qr->nrQMlayers; j++)
1069 qm->shiftQM[0] = XYZ2IS(0, 0, 0);
1070 for (i = 1; i < qm->nrQMatoms; i++)
1072 qm->shiftQM[i] = pbc_dx_aiuc(&pbc, x[qm->indexQM[0]], x[qm->indexQM[i]],
1075 update_QMMM_coord(x, fr, qm, mm);
1078 } /* update_QMMM_rec */
1081 real calculate_QMMM(t_commrec *cr,
1087 /* a selection for the QM package depending on which is requested
1088 * (Gaussian, GAMESS-UK, MOPAC or ORCA) needs to be implemented here. Now
1089 * it works through defines.... Not so nice yet
1098 *forces = NULL, *fshift = NULL,
1099 *forces2 = NULL, *fshift2 = NULL; /* needed for multilayer ONIOM */
1102 /* make a local copy the QMMMrec pointer
1107 /* now different procedures are carried out for one layer ONION and
1108 * normal QMMM on one hand and multilayer oniom on the other
1110 if (qr->QMMMscheme == eQMMMschemenormal || qr->nrQMlayers == 1)
1113 snew(forces, (qm->nrQMatoms+mm->nrMMatoms));
1114 snew(fshift, (qm->nrQMatoms+mm->nrMMatoms));
1115 QMener = call_QMroutine(cr, fr, qm, mm, forces, fshift);
1116 for (i = 0; i < qm->nrQMatoms; i++)
1118 for (j = 0; j < DIM; j++)
1120 f[qm->indexQM[i]][j] -= forces[i][j];
1121 fr->fshift[qm->shiftQM[i]][j] += fshift[i][j];
1124 for (i = 0; i < mm->nrMMatoms; i++)
1126 for (j = 0; j < DIM; j++)
1128 f[mm->indexMM[i]][j] -= forces[qm->nrQMatoms+i][j];
1129 fr->fshift[mm->shiftMM[i]][j] += fshift[qm->nrQMatoms+i][j];
1136 else /* Multi-layer ONIOM */
1138 for (i = 0; i < qr->nrQMlayers-1; i++) /* last layer is special */
1141 qm2 = copy_QMrec(qr->qm[i+1]);
1143 qm2->nrQMatoms = qm->nrQMatoms;
1145 for (j = 0; j < qm2->nrQMatoms; j++)
1147 for (k = 0; k < DIM; k++)
1149 qm2->xQM[j][k] = qm->xQM[j][k];
1151 qm2->indexQM[j] = qm->indexQM[j];
1152 qm2->atomicnumberQM[j] = qm->atomicnumberQM[j];
1153 qm2->shiftQM[j] = qm->shiftQM[j];
1156 qm2->QMcharge = qm->QMcharge;
1157 /* this layer at the higher level of theory */
1158 srenew(forces, qm->nrQMatoms);
1159 srenew(fshift, qm->nrQMatoms);
1160 /* we need to re-initialize the QMroutine every step... */
1161 init_QMroutine(cr, qm, mm);
1162 QMener += call_QMroutine(cr, fr, qm, mm, forces, fshift);
1164 /* this layer at the lower level of theory */
1165 srenew(forces2, qm->nrQMatoms);
1166 srenew(fshift2, qm->nrQMatoms);
1167 init_QMroutine(cr, qm2, mm);
1168 QMener -= call_QMroutine(cr, fr, qm2, mm, forces2, fshift2);
1169 /* E = E1high-E1low The next layer includes the current layer at
1170 * the lower level of theory, which provides + E2low
1171 * this is similar for gradients
1173 for (i = 0; i < qm->nrQMatoms; i++)
1175 for (j = 0; j < DIM; j++)
1177 f[qm->indexQM[i]][j] -= (forces[i][j]-forces2[i][j]);
1178 fr->fshift[qm->shiftQM[i]][j] += (fshift[i][j]-fshift2[i][j]);
1183 /* now the last layer still needs to be done: */
1184 qm = qr->qm[qr->nrQMlayers-1]; /* C counts from 0 */
1185 init_QMroutine(cr, qm, mm);
1186 srenew(forces, qm->nrQMatoms);
1187 srenew(fshift, qm->nrQMatoms);
1188 QMener += call_QMroutine(cr, fr, qm, mm, forces, fshift);
1189 for (i = 0; i < qm->nrQMatoms; i++)
1191 for (j = 0; j < DIM; j++)
1193 f[qm->indexQM[i]][j] -= forces[i][j];
1194 fr->fshift[qm->shiftQM[i]][j] += fshift[i][j];
1202 if (qm->bTS || qm->bOPT)
1204 /* qm[0] still contains the largest ONIOM QM subsystem
1205 * we take the optimized coordiates and put the in x[]
1207 for (i = 0; i < qm->nrQMatoms; i++)
1209 for (j = 0; j < DIM; j++)
1211 x[qm->indexQM[i]][j] = qm->xQM[i][j];
1216 } /* calculate_QMMM */
1218 /* end of QMMM core routines */