3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * GROwing Monsters And Cloning Shrimps
49 #include "gromacs/fileio/confio.h"
61 #include "gmx_fatal.h"
64 #include "mtop_util.h"
67 /* declarations of the interfaces to the QM packages. The _SH indicate
68 * the QM interfaces can be used for Surface Hopping simulations
70 #ifdef GMX_QMMM_GAMESS
71 /* GAMESS interface */
74 init_gamess(t_commrec *cr, t_QMrec *qm, t_MMrec *mm);
77 call_gamess(t_commrec *cr, t_forcerec *fr,
78 t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[]);
80 #elif defined GMX_QMMM_MOPAC
84 init_mopac(t_commrec *cr, t_QMrec *qm, t_MMrec *mm);
87 call_mopac(t_commrec *cr, t_forcerec *fr, t_QMrec *qm,
88 t_MMrec *mm, rvec f[], rvec fshift[]);
91 call_mopac_SH(t_commrec *cr, t_forcerec *fr, t_QMrec *qm,
92 t_MMrec *mm, rvec f[], rvec fshift[]);
94 #elif defined GMX_QMMM_GAUSSIAN
95 /* GAUSSIAN interface */
98 init_gaussian(t_commrec *cr, t_QMrec *qm, t_MMrec *mm);
101 call_gaussian_SH(t_commrec *cr, t_forcerec *fr, t_QMrec *qm,
102 t_MMrec *mm, rvec f[], rvec fshift[]);
105 call_gaussian(t_commrec *cr, t_forcerec *fr, t_QMrec *qm,
106 t_MMrec *mm, rvec f[], rvec fshift[]);
108 #elif defined GMX_QMMM_ORCA
112 init_orca(t_QMrec *qm);
115 call_orca(t_forcerec *fr, t_QMrec *qm,
116 t_MMrec *mm, rvec f[], rvec fshift[]);
123 /* this struct and these comparison functions are needed for creating
124 * a QMMM input for the QM routines from the QMMM neighbor list.
132 static int struct_comp(const void *a, const void *b)
135 return (int)(((t_j_particle *)a)->j)-(int)(((t_j_particle *)b)->j);
139 static int int_comp(const void *a, const void *b)
142 return (*(int *)a) - (*(int *)b);
146 static int QMlayer_comp(const void *a, const void *b)
149 return (int)(((t_QMrec *)a)->nrQMatoms)-(int)(((t_QMrec *)b)->nrQMatoms);
153 real call_QMroutine(t_commrec gmx_unused *cr, t_forcerec gmx_unused *fr, t_QMrec gmx_unused *qm,
154 t_MMrec gmx_unused *mm, rvec gmx_unused f[], rvec gmx_unused fshift[])
156 /* makes a call to the requested QM routine (qm->QMmethod)
157 * Note that f is actually the gradient, i.e. -f
162 /* do a semi-empiprical calculation */
164 if (qm->QMmethod < eQMmethodRHF && !(mm->nrMMatoms))
166 #ifdef GMX_QMMM_MOPAC
169 QMener = call_mopac_SH(cr, fr, qm, mm, f, fshift);
173 QMener = call_mopac(cr, fr, qm, mm, f, fshift);
176 gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac.");
181 /* do an ab-initio calculation */
182 if (qm->bSH && qm->QMmethod == eQMmethodCASSCF)
184 #ifdef GMX_QMMM_GAUSSIAN
185 QMener = call_gaussian_SH(cr, fr, qm, mm, f, fshift);
187 gmx_fatal(FARGS, "Ab-initio Surface-hopping only supported with Gaussian.");
192 #ifdef GMX_QMMM_GAMESS
193 QMener = call_gamess(cr, fr, qm, mm, f, fshift);
194 #elif defined GMX_QMMM_GAUSSIAN
195 QMener = call_gaussian(cr, fr, qm, mm, f, fshift);
196 #elif defined GMX_QMMM_ORCA
197 QMener = call_orca(cr, fr, qm, mm, f, fshift);
199 gmx_fatal(FARGS, "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
206 void init_QMroutine(t_commrec gmx_unused *cr, t_QMrec gmx_unused *qm, t_MMrec gmx_unused *mm)
208 /* makes a call to the requested QM routine (qm->QMmethod)
210 if (qm->QMmethod < eQMmethodRHF)
212 #ifdef GMX_QMMM_MOPAC
213 /* do a semi-empiprical calculation */
214 init_mopac(cr, qm, mm);
216 gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac.");
221 /* do an ab-initio calculation */
222 #ifdef GMX_QMMM_GAMESS
223 init_gamess(cr, qm, mm);
224 #elif defined GMX_QMMM_GAUSSIAN
225 init_gaussian(cr, qm, mm);
226 #elif defined GMX_QMMM_ORCA
229 gmx_fatal(FARGS, "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
232 } /* init_QMroutine */
234 void update_QMMM_coord(rvec x[], t_forcerec *fr, t_QMrec *qm, t_MMrec *mm)
236 /* shifts the QM and MM particles into the central box and stores
237 * these shifted coordinates in the coordinate arrays of the
238 * QMMMrec. These coordinates are passed on the QM subroutines.
243 /* shift the QM atoms into the central box
245 for (i = 0; i < qm->nrQMatoms; i++)
247 rvec_sub(x[qm->indexQM[i]], fr->shift_vec[qm->shiftQM[i]], qm->xQM[i]);
249 /* also shift the MM atoms into the central box, if any
251 for (i = 0; i < mm->nrMMatoms; i++)
253 rvec_sub(x[mm->indexMM[i]], fr->shift_vec[mm->shiftMM[i]], mm->xMM[i]);
255 } /* update_QMMM_coord */
257 static void punch_QMMM_excl(t_QMrec *qm, t_MMrec *mm, t_blocka *excls)
259 /* punch a file containing the bonded interactions of each QM
260 * atom with MM atoms. These need to be excluded in the QM routines
261 * Only needed in case of QM/MM optimizations
266 i, j, k, nrexcl = 0, *excluded = NULL, max = 0;
269 out = fopen("QMMMexcl.dat", "w");
271 /* this can be done more efficiently I think
273 for (i = 0; i < qm->nrQMatoms; i++)
276 for (j = excls->index[qm->indexQM[i]];
277 j < excls->index[qm->indexQM[i]+1];
280 for (k = 0; k < mm->nrMMatoms; k++)
282 if (mm->indexMM[k] == excls->a[j]) /* the excluded MM atom */
287 srenew(excluded, max);
289 excluded[nrexcl++] = k;
295 fprintf(out, "%5d %5d\n", i+1, nrexcl);
296 for (j = 0; j < nrexcl; j++)
298 fprintf(out, "%5d ", excluded[j]);
304 } /* punch_QMMM_excl */
307 /* end of QMMM subroutines */
309 /* QMMM core routines */
311 t_QMrec *mk_QMrec(void)
318 t_MMrec *mk_MMrec(void)
325 static void init_QMrec(int grpnr, t_QMrec *qm, int nr, int *atomarray,
326 gmx_mtop_t *mtop, t_inputrec *ir)
328 /* fills the t_QMrec struct of QM group grpnr
331 gmx_mtop_atomlookup_t alook;
337 snew(qm->indexQM, nr);
338 snew(qm->shiftQM, nr); /* the shifts */
339 for (i = 0; i < nr; i++)
341 qm->indexQM[i] = atomarray[i];
344 alook = gmx_mtop_atomlookup_init(mtop);
346 snew(qm->atomicnumberQM, nr);
347 for (i = 0; i < qm->nrQMatoms; i++)
349 gmx_mtop_atomnr_to_atom(alook, qm->indexQM[i], &atom);
350 qm->nelectrons += mtop->atomtypes.atomnumber[atom->type];
351 qm->atomicnumberQM[i] = mtop->atomtypes.atomnumber[atom->type];
354 gmx_mtop_atomlookup_destroy(alook);
356 qm->QMcharge = ir->opts.QMcharge[grpnr];
357 qm->multiplicity = ir->opts.QMmult[grpnr];
358 qm->nelectrons -= ir->opts.QMcharge[grpnr];
360 qm->QMmethod = ir->opts.QMmethod[grpnr];
361 qm->QMbasis = ir->opts.QMbasis[grpnr];
362 /* trajectory surface hopping setup (Gaussian only) */
363 qm->bSH = ir->opts.bSH[grpnr];
364 qm->CASorbitals = ir->opts.CASorbitals[grpnr];
365 qm->CASelectrons = ir->opts.CASelectrons[grpnr];
366 qm->SAsteps = ir->opts.SAsteps[grpnr];
367 qm->SAon = ir->opts.SAon[grpnr];
368 qm->SAoff = ir->opts.SAoff[grpnr];
369 /* hack to prevent gaussian from reinitializing all the time */
370 qm->nQMcpus = 0; /* number of CPU's to be used by g01, is set
371 * upon initializing gaussian
374 /* print the current layer to allow users to check their input */
375 fprintf(stderr, "Layer %d\nnr of QM atoms %d\n", grpnr, nr);
376 fprintf(stderr, "QMlevel: %s/%s\n\n",
377 eQMmethod_names[qm->QMmethod], eQMbasis_names[qm->QMbasis]);
380 snew(qm->frontatoms, nr);
381 /* Lennard-Jones coefficients */
384 /* do we optimize the QM separately using the algorithms of the QM program??
386 qm->bTS = ir->opts.bTS[grpnr];
387 qm->bOPT = ir->opts.bOPT[grpnr];
391 t_QMrec *copy_QMrec(t_QMrec *qm)
393 /* copies the contents of qm into a new t_QMrec struct */
400 qmcopy->nrQMatoms = qm->nrQMatoms;
401 snew(qmcopy->xQM, qmcopy->nrQMatoms);
402 snew(qmcopy->indexQM, qmcopy->nrQMatoms);
403 snew(qmcopy->atomicnumberQM, qm->nrQMatoms);
404 snew(qmcopy->shiftQM, qmcopy->nrQMatoms); /* the shifts */
405 for (i = 0; i < qmcopy->nrQMatoms; i++)
407 qmcopy->shiftQM[i] = qm->shiftQM[i];
408 qmcopy->indexQM[i] = qm->indexQM[i];
409 qmcopy->atomicnumberQM[i] = qm->atomicnumberQM[i];
411 qmcopy->nelectrons = qm->nelectrons;
412 qmcopy->multiplicity = qm->multiplicity;
413 qmcopy->QMcharge = qm->QMcharge;
414 qmcopy->nelectrons = qm->nelectrons;
415 qmcopy->QMmethod = qm->QMmethod;
416 qmcopy->QMbasis = qm->QMbasis;
417 /* trajectory surface hopping setup (Gaussian only) */
418 qmcopy->bSH = qm->bSH;
419 qmcopy->CASorbitals = qm->CASorbitals;
420 qmcopy->CASelectrons = qm->CASelectrons;
421 qmcopy->SAsteps = qm->SAsteps;
422 qmcopy->SAon = qm->SAon;
423 qmcopy->SAoff = qm->SAoff;
424 qmcopy->bOPT = qm->bOPT;
426 /* Gaussian init. variables */
427 qmcopy->nQMcpus = qm->nQMcpus;
428 for (i = 0; i < DIM; i++)
430 qmcopy->SHbasis[i] = qm->SHbasis[i];
432 qmcopy->QMmem = qm->QMmem;
433 qmcopy->accuracy = qm->accuracy;
434 qmcopy->cpmcscf = qm->cpmcscf;
435 qmcopy->SAstep = qm->SAstep;
436 snew(qmcopy->frontatoms, qm->nrQMatoms);
437 snew(qmcopy->c12, qmcopy->nrQMatoms);
438 snew(qmcopy->c6, qmcopy->nrQMatoms);
439 if (qmcopy->bTS || qmcopy->bOPT)
441 for (i = 1; i < qmcopy->nrQMatoms; i++)
443 qmcopy->frontatoms[i] = qm->frontatoms[i];
444 qmcopy->c12[i] = qm->c12[i];
445 qmcopy->c6[i] = qm->c6[i];
453 t_QMMMrec *mk_QMMMrec(void)
464 void init_QMMMrec(t_commrec *cr,
469 /* we put the atomsnumbers of atoms that belong to the QMMM group in
470 * an array that will be copied later to QMMMrec->indexQM[..]. Also
471 * it will be used to create an QMMMrec->bQMMM index array that
472 * simply contains true/false for QM and MM (the other) atoms.
475 gmx_groups_t *groups;
476 atom_id *qm_arr = NULL, vsite, ai, aj;
477 int qm_max = 0, qm_nr = 0, i, j, jmax, k, l, nrvsite2 = 0;
482 gmx_mtop_atomloop_all_t aloop;
484 gmx_mtop_ilistloop_all_t iloop;
487 gmx_mtop_atomlookup_t alook;
489 c6au = (HARTREE2KJ*AVOGADRO*pow(BOHR2NM, 6));
490 c12au = (HARTREE2KJ*AVOGADRO*pow(BOHR2NM, 12));
491 /* issue a fatal if the user wants to run with more than one node */
494 gmx_fatal(FARGS, "QM/MM does not work in parallel, use a single node instead\n");
497 /* Make a local copy of the QMMMrec */
500 /* bQMMM[..] is an array containing TRUE/FALSE for atoms that are
501 * QM/not QM. We first set all elemenst at false. Afterwards we use
502 * the qm_arr (=MMrec->indexQM) to changes the elements
503 * corresponding to the QM atoms at TRUE. */
505 qr->QMMMscheme = ir->QMMMscheme;
507 /* we take the possibility into account that a user has
508 * defined more than one QM group:
510 /* an ugly work-around in case there is only one group In this case
511 * the whole system is treated as QM. Otherwise the second group is
512 * always the rest of the total system and is treated as MM.
515 /* small problem if there is only QM.... so no MM */
517 jmax = ir->opts.ngQM;
519 if (qr->QMMMscheme == eQMMMschemeoniom)
521 qr->nrQMlayers = jmax;
528 groups = &mtop->groups;
530 /* there are jmax groups of QM atoms. In case of multiple QM groups
531 * I assume that the users wants to do ONIOM. However, maybe it
532 * should also be possible to define more than one QM subsystem with
533 * independent neighbourlists. I have to think about
537 for (j = 0; j < jmax; j++)
540 aloop = gmx_mtop_atomloop_all_init(mtop);
541 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
546 srenew(qm_arr, qm_max);
548 if (ggrpnr(groups, egcQMMM, i) == j)
554 if (qr->QMMMscheme == eQMMMschemeoniom)
556 /* add the atoms to the bQMMM array
559 /* I assume that users specify the QM groups from small to
560 * big(ger) in the mdp file
562 qr->qm[j] = mk_QMrec();
563 /* we need to throw out link atoms that in the previous layer
564 * existed to separate this QMlayer from the previous
565 * QMlayer. We use the iatoms array in the idef for that
566 * purpose. If all atoms defining the current Link Atom (Dummy2)
567 * are part of the current QM layer it needs to be removed from
570 iloop = gmx_mtop_ilistloop_all_init(mtop);
571 while (gmx_mtop_ilistloop_all_next(iloop, &ilist_mol, &a_offset))
573 nrvsite2 = ilist_mol[F_VSITE2].nr;
574 iatoms = ilist_mol[F_VSITE2].iatoms;
576 for (k = 0; k < nrvsite2; k += 4)
578 vsite = a_offset + iatoms[k+1]; /* the vsite */
579 ai = a_offset + iatoms[k+2]; /* constructing atom */
580 aj = a_offset + iatoms[k+3]; /* constructing atom */
581 if (ggrpnr(groups, egcQMMM, vsite) == ggrpnr(groups, egcQMMM, ai)
583 ggrpnr(groups, egcQMMM, vsite) == ggrpnr(groups, egcQMMM, aj))
585 /* this dummy link atom needs to be removed from the qm_arr
586 * before making the QMrec of this layer!
588 for (i = 0; i < qm_nr; i++)
590 if (qm_arr[i] == vsite)
592 /* drop the element */
593 for (l = i; l < qm_nr; l++)
595 qm_arr[l] = qm_arr[l+1];
604 /* store QM atoms in this layer in the QMrec and initialise layer
606 init_QMrec(j, qr->qm[j], qm_nr, qm_arr, mtop, ir);
608 /* we now store the LJ C6 and C12 parameters in QM rec in case
609 * we need to do an optimization
611 if (qr->qm[j]->bOPT || qr->qm[j]->bTS)
613 for (i = 0; i < qm_nr; i++)
615 /* nbfp now includes the 6.0/12.0 derivative prefactors */
616 qr->qm[j]->c6[i] = C6(fr->nbfp, mtop->ffparams.atnr, atom->type, atom->type)/c6au/6.0;
617 qr->qm[j]->c12[i] = C12(fr->nbfp, mtop->ffparams.atnr, atom->type, atom->type)/c12au/12.0;
620 /* now we check for frontier QM atoms. These occur in pairs that
621 * construct the vsite
623 iloop = gmx_mtop_ilistloop_all_init(mtop);
624 while (gmx_mtop_ilistloop_all_next(iloop, &ilist_mol, &a_offset))
626 nrvsite2 = ilist_mol[F_VSITE2].nr;
627 iatoms = ilist_mol[F_VSITE2].iatoms;
629 for (k = 0; k < nrvsite2; k += 4)
631 vsite = a_offset + iatoms[k+1]; /* the vsite */
632 ai = a_offset + iatoms[k+2]; /* constructing atom */
633 aj = a_offset + iatoms[k+3]; /* constructing atom */
634 if (ggrpnr(groups, egcQMMM, ai) < (groups->grps[egcQMMM].nr-1) &&
635 (ggrpnr(groups, egcQMMM, aj) >= (groups->grps[egcQMMM].nr-1)))
637 /* mark ai as frontier atom */
638 for (i = 0; i < qm_nr; i++)
640 if ( (qm_arr[i] == ai) || (qm_arr[i] == vsite) )
642 qr->qm[j]->frontatoms[i] = TRUE;
646 else if (ggrpnr(groups, egcQMMM, aj) < (groups->grps[egcQMMM].nr-1) &&
647 (ggrpnr(groups, egcQMMM, ai) >= (groups->grps[egcQMMM].nr-1)))
649 /* mark aj as frontier atom */
650 for (i = 0; i < qm_nr; i++)
652 if ( (qm_arr[i] == aj) || (qm_arr[i] == vsite))
654 qr->qm[j]->frontatoms[i] = TRUE;
662 if (qr->QMMMscheme != eQMMMschemeoniom)
665 /* standard QMMM, all layers are merged together so there is one QM
666 * subsystem and one MM subsystem.
667 * Also we set the charges to zero in the md->charge arrays to prevent
668 * the innerloops from doubly counting the electostatic QM MM interaction
671 alook = gmx_mtop_atomlookup_init(mtop);
673 for (k = 0; k < qm_nr; k++)
675 gmx_mtop_atomnr_to_atom(alook, qm_arr[k], &atom);
679 qr->qm[0] = mk_QMrec();
680 /* store QM atoms in the QMrec and initialise
682 init_QMrec(0, qr->qm[0], qm_nr, qm_arr, mtop, ir);
683 if (qr->qm[0]->bOPT || qr->qm[0]->bTS)
685 for (i = 0; i < qm_nr; i++)
687 gmx_mtop_atomnr_to_atom(alook, qm_arr[i], &atom);
688 /* nbfp now includes the 6.0/12.0 derivative prefactors */
689 qr->qm[0]->c6[i] = C6(fr->nbfp, mtop->ffparams.atnr, atom->type, atom->type)/c6au/6.0;
690 qr->qm[0]->c12[i] = C12(fr->nbfp, mtop->ffparams.atnr, atom->type, atom->type)/c12au/12.0;
694 /* find frontier atoms and mark them true in the frontieratoms array.
696 for (i = 0; i < qm_nr; i++)
698 gmx_mtop_atomnr_to_ilist(alook, qm_arr[i], &ilist_mol, &a_offset);
699 nrvsite2 = ilist_mol[F_VSITE2].nr;
700 iatoms = ilist_mol[F_VSITE2].iatoms;
702 for (k = 0; k < nrvsite2; k += 4)
704 vsite = a_offset + iatoms[k+1]; /* the vsite */
705 ai = a_offset + iatoms[k+2]; /* constructing atom */
706 aj = a_offset + iatoms[k+3]; /* constructing atom */
707 if (ggrpnr(groups, egcQMMM, ai) < (groups->grps[egcQMMM].nr-1) &&
708 (ggrpnr(groups, egcQMMM, aj) >= (groups->grps[egcQMMM].nr-1)))
710 /* mark ai as frontier atom */
711 if ( (qm_arr[i] == ai) || (qm_arr[i] == vsite) )
713 qr->qm[0]->frontatoms[i] = TRUE;
716 else if (ggrpnr(groups, egcQMMM, aj) < (groups->grps[egcQMMM].nr-1) &&
717 (ggrpnr(groups, egcQMMM, ai) >= (groups->grps[egcQMMM].nr-1)))
719 /* mark aj as frontier atom */
720 if ( (qm_arr[i] == aj) || (qm_arr[i] == vsite) )
722 qr->qm[0]->frontatoms[i] = TRUE;
728 gmx_mtop_atomlookup_destroy(alook);
730 /* MM rec creation */
732 mm->scalefactor = ir->scalefactor;
733 mm->nrMMatoms = (mtop->natoms)-(qr->qm[0]->nrQMatoms); /* rest of the atoms */
737 { /* MM rec creation */
739 mm->scalefactor = ir->scalefactor;
744 /* these variables get updated in the update QMMMrec */
746 if (qr->nrQMlayers == 1)
748 /* with only one layer there is only one initialisation
749 * needed. Multilayer is a bit more complicated as it requires
750 * re-initialisation at every step of the simulation. This is due
751 * to the use of COMMON blocks in the fortran QM subroutines.
753 if (qr->qm[0]->QMmethod < eQMmethodRHF)
755 #ifdef GMX_QMMM_MOPAC
756 /* semi-empiprical 1-layer ONIOM calculation requested (mopac93) */
757 init_mopac(cr, qr->qm[0], qr->mm);
759 gmx_fatal(FARGS, "Semi-empirical QM only supported with Mopac.");
764 /* ab initio calculation requested (gamess/gaussian/ORCA) */
765 #ifdef GMX_QMMM_GAMESS
766 init_gamess(cr, qr->qm[0], qr->mm);
767 #elif defined GMX_QMMM_GAUSSIAN
768 init_gaussian(cr, qr->qm[0], qr->mm);
769 #elif defined GMX_QMMM_ORCA
770 init_orca(cr, qr->qm[0], qr->mm);
772 gmx_fatal(FARGS, "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
778 void update_QMMMrec(t_commrec *cr,
785 /* updates the coordinates of both QM atoms and MM atoms and stores
786 * them in the QMMMrec.
788 * NOTE: is NOT yet working if there are no PBC. Also in ns.c, simple
789 * ns needs to be fixed!
792 mm_max = 0, mm_nr = 0, mm_nr_new, i, j, is, k, shift;
794 *mm_j_particles = NULL, *qm_i_particles = NULL;
810 *parallelMMarray = NULL;
814 c6au = (HARTREE2KJ*AVOGADRO*pow(BOHR2NM, 6));
815 c12au = (HARTREE2KJ*AVOGADRO*pow(BOHR2NM, 12));
817 /* every cpu has this array. On every processor we fill this array
818 * with 1's and 0's. 1's indicate the atoms is a QM atom on the
819 * current cpu in a later stage these arrays are all summed. indexes
820 * > 0 indicate the atom is a QM atom. Every node therefore knows
821 * whcih atoms are part of the QM subsystem.
823 /* copy some pointers */
826 QMMMlist = fr->QMMMlist;
830 /* init_pbc(box); needs to be called first, see pbc.h */
831 set_pbc_dd(&pbc, fr->ePBC, DOMAINDECOMP(cr) ? cr->dd : NULL, FALSE, box);
832 /* only in standard (normal) QMMM we need the neighbouring MM
833 * particles to provide a electric field of point charges for the QM
836 if (qr->QMMMscheme == eQMMMschemenormal) /* also implies 1 QM-layer */
838 /* we NOW create/update a number of QMMMrec entries:
840 * 1) the shiftQM, containing the shifts of the QM atoms
842 * 2) the indexMM array, containing the index of the MM atoms
844 * 3) the shiftMM, containing the shifts of the MM atoms
846 * 4) the shifted coordinates of the MM atoms
848 * the shifts are used for computing virial of the QM/MM particles.
850 qm = qr->qm[0]; /* in case of normal QMMM, there is only one group */
851 snew(qm_i_particles, QMMMlist.nri);
854 qm_i_particles[0].shift = XYZ2IS(0, 0, 0);
855 for (i = 0; i < QMMMlist.nri; i++)
857 qm_i_particles[i].j = QMMMlist.iinr[i];
861 qm_i_particles[i].shift = pbc_dx_aiuc(&pbc, x[QMMMlist.iinr[0]],
862 x[QMMMlist.iinr[i]], dx);
865 /* However, since nri >= nrQMatoms, we do a quicksort, and throw
866 * out double, triple, etc. entries later, as we do for the MM
870 /* compute the shift for the MM j-particles with respect to
871 * the QM i-particle and store them.
874 crd[0] = IS2X(QMMMlist.shift[i]) + IS2X(qm_i_particles[i].shift);
875 crd[1] = IS2Y(QMMMlist.shift[i]) + IS2Y(qm_i_particles[i].shift);
876 crd[2] = IS2Z(QMMMlist.shift[i]) + IS2Z(qm_i_particles[i].shift);
877 is = XYZ2IS(crd[0], crd[1], crd[2]);
878 for (j = QMMMlist.jindex[i];
879 j < QMMMlist.jindex[i+1];
885 srenew(mm_j_particles, mm_max);
888 mm_j_particles[mm_nr].j = QMMMlist.jjnr[j];
889 mm_j_particles[mm_nr].shift = is;
894 /* quicksort QM and MM shift arrays and throw away multiple entries */
898 qsort(qm_i_particles, QMMMlist.nri,
899 (size_t)sizeof(qm_i_particles[0]),
901 qsort(mm_j_particles, mm_nr,
902 (size_t)sizeof(mm_j_particles[0]),
904 /* remove multiples in the QM shift array, since in init_QMMM() we
905 * went through the atom numbers from 0 to md.nr, the order sorted
906 * here matches the one of QMindex already.
909 for (i = 0; i < QMMMlist.nri; i++)
911 if (i == 0 || qm_i_particles[i].j != qm_i_particles[i-1].j)
913 qm_i_particles[j++] = qm_i_particles[i];
917 if (qm->bTS || qm->bOPT)
919 /* only remove double entries for the MM array */
920 for (i = 0; i < mm_nr; i++)
922 if ((i == 0 || mm_j_particles[i].j != mm_j_particles[i-1].j)
923 && !md->bQM[mm_j_particles[i].j])
925 mm_j_particles[mm_nr_new++] = mm_j_particles[i];
929 /* we also remove mm atoms that have no charges!
930 * actually this is already done in the ns.c
934 for (i = 0; i < mm_nr; i++)
936 if ((i == 0 || mm_j_particles[i].j != mm_j_particles[i-1].j)
937 && !md->bQM[mm_j_particles[i].j]
938 && (md->chargeA[mm_j_particles[i].j]
939 || (md->chargeB && md->chargeB[mm_j_particles[i].j])))
941 mm_j_particles[mm_nr_new++] = mm_j_particles[i];
946 /* store the data retrieved above into the QMMMrec
949 /* Keep the compiler happy,
950 * shift will always be set in the loop for i=0
953 for (i = 0; i < qm->nrQMatoms; i++)
955 /* not all qm particles might have appeared as i
956 * particles. They might have been part of the same charge
957 * group for instance.
959 if (qm->indexQM[i] == qm_i_particles[k].j)
961 shift = qm_i_particles[k++].shift;
963 /* use previous shift, assuming they belong the same charge
967 qm->shiftQM[i] = shift;
970 /* parallel excecution */
973 snew(parallelMMarray, 2*(md->nr));
974 /* only MM particles have a 1 at their atomnumber. The second part
975 * of the array contains the shifts. Thus:
976 * p[i]=1/0 depending on wether atomnumber i is a MM particle in the QM
977 * step or not. p[i+md->nr] is the shift of atomnumber i.
979 for (i = 0; i < 2*(md->nr); i++)
981 parallelMMarray[i] = 0;
984 for (i = 0; i < mm_nr; i++)
986 parallelMMarray[mm_j_particles[i].j] = 1;
987 parallelMMarray[mm_j_particles[i].j+(md->nr)] = mm_j_particles[i].shift;
989 gmx_sumi(md->nr, parallelMMarray, cr);
993 for (i = 0; i < md->nr; i++)
995 if (parallelMMarray[i])
1000 srenew(mm->indexMM, mm_max);
1001 srenew(mm->shiftMM, mm_max);
1003 mm->indexMM[mm_nr] = i;
1004 mm->shiftMM[mm_nr++] = parallelMMarray[i+md->nr]/parallelMMarray[i];
1007 mm->nrMMatoms = mm_nr;
1008 free(parallelMMarray);
1010 /* serial execution */
1013 mm->nrMMatoms = mm_nr;
1014 srenew(mm->shiftMM, mm_nr);
1015 srenew(mm->indexMM, mm_nr);
1016 for (i = 0; i < mm_nr; i++)
1018 mm->indexMM[i] = mm_j_particles[i].j;
1019 mm->shiftMM[i] = mm_j_particles[i].shift;
1023 /* (re) allocate memory for the MM coordiate array. The QM
1024 * coordinate array was already allocated in init_QMMM, and is
1025 * only (re)filled in the update_QMMM_coordinates routine
1027 srenew(mm->xMM, mm->nrMMatoms);
1028 /* now we (re) fill the array that contains the MM charges with
1029 * the forcefield charges. If requested, these charges will be
1030 * scaled by a factor
1032 srenew(mm->MMcharges, mm->nrMMatoms);
1033 for (i = 0; i < mm->nrMMatoms; i++) /* no free energy yet */
1035 mm->MMcharges[i] = md->chargeA[mm->indexMM[i]]*mm->scalefactor;
1037 if (qm->bTS || qm->bOPT)
1039 /* store (copy) the c6 and c12 parameters into the MMrec struct
1041 srenew(mm->c6, mm->nrMMatoms);
1042 srenew(mm->c12, mm->nrMMatoms);
1043 for (i = 0; i < mm->nrMMatoms; i++)
1045 /* nbfp now includes the 6.0/12.0 derivative prefactors */
1046 mm->c6[i] = C6(fr->nbfp, top->idef.atnr, md->typeA[mm->indexMM[i]], md->typeA[mm->indexMM[i]])/c6au/6.0;
1047 mm->c12[i] = C12(fr->nbfp, top->idef.atnr, md->typeA[mm->indexMM[i]], md->typeA[mm->indexMM[i]])/c12au/12.0;
1049 punch_QMMM_excl(qr->qm[0], mm, &(top->excls));
1051 /* the next routine fills the coordinate fields in the QMMM rec of
1052 * both the qunatum atoms and the MM atoms, using the shifts
1056 update_QMMM_coord(x, fr, qr->qm[0], qr->mm);
1057 free(qm_i_particles);
1058 free(mm_j_particles);
1060 else /* ONIOM */ /* ????? */
1063 /* do for each layer */
1064 for (j = 0; j < qr->nrQMlayers; j++)
1067 qm->shiftQM[0] = XYZ2IS(0, 0, 0);
1068 for (i = 1; i < qm->nrQMatoms; i++)
1070 qm->shiftQM[i] = pbc_dx_aiuc(&pbc, x[qm->indexQM[0]], x[qm->indexQM[i]],
1073 update_QMMM_coord(x, fr, qm, mm);
1076 } /* update_QMMM_rec */
1079 real calculate_QMMM(t_commrec *cr,
1085 /* a selection for the QM package depending on which is requested
1086 * (Gaussian, GAMESS-UK, MOPAC or ORCA) needs to be implemented here. Now
1087 * it works through defines.... Not so nice yet
1096 *forces = NULL, *fshift = NULL,
1097 *forces2 = NULL, *fshift2 = NULL; /* needed for multilayer ONIOM */
1100 /* make a local copy the QMMMrec pointer
1105 /* now different procedures are carried out for one layer ONION and
1106 * normal QMMM on one hand and multilayer oniom on the other
1108 if (qr->QMMMscheme == eQMMMschemenormal || qr->nrQMlayers == 1)
1111 snew(forces, (qm->nrQMatoms+mm->nrMMatoms));
1112 snew(fshift, (qm->nrQMatoms+mm->nrMMatoms));
1113 QMener = call_QMroutine(cr, fr, qm, mm, forces, fshift);
1114 for (i = 0; i < qm->nrQMatoms; i++)
1116 for (j = 0; j < DIM; j++)
1118 f[qm->indexQM[i]][j] -= forces[i][j];
1119 fr->fshift[qm->shiftQM[i]][j] += fshift[i][j];
1122 for (i = 0; i < mm->nrMMatoms; i++)
1124 for (j = 0; j < DIM; j++)
1126 f[mm->indexMM[i]][j] -= forces[qm->nrQMatoms+i][j];
1127 fr->fshift[mm->shiftMM[i]][j] += fshift[qm->nrQMatoms+i][j];
1134 else /* Multi-layer ONIOM */
1136 for (i = 0; i < qr->nrQMlayers-1; i++) /* last layer is special */
1139 qm2 = copy_QMrec(qr->qm[i+1]);
1141 qm2->nrQMatoms = qm->nrQMatoms;
1143 for (j = 0; j < qm2->nrQMatoms; j++)
1145 for (k = 0; k < DIM; k++)
1147 qm2->xQM[j][k] = qm->xQM[j][k];
1149 qm2->indexQM[j] = qm->indexQM[j];
1150 qm2->atomicnumberQM[j] = qm->atomicnumberQM[j];
1151 qm2->shiftQM[j] = qm->shiftQM[j];
1154 qm2->QMcharge = qm->QMcharge;
1155 /* this layer at the higher level of theory */
1156 srenew(forces, qm->nrQMatoms);
1157 srenew(fshift, qm->nrQMatoms);
1158 /* we need to re-initialize the QMroutine every step... */
1159 init_QMroutine(cr, qm, mm);
1160 QMener += call_QMroutine(cr, fr, qm, mm, forces, fshift);
1162 /* this layer at the lower level of theory */
1163 srenew(forces2, qm->nrQMatoms);
1164 srenew(fshift2, qm->nrQMatoms);
1165 init_QMroutine(cr, qm2, mm);
1166 QMener -= call_QMroutine(cr, fr, qm2, mm, forces2, fshift2);
1167 /* E = E1high-E1low The next layer includes the current layer at
1168 * the lower level of theory, which provides + E2low
1169 * this is similar for gradients
1171 for (i = 0; i < qm->nrQMatoms; i++)
1173 for (j = 0; j < DIM; j++)
1175 f[qm->indexQM[i]][j] -= (forces[i][j]-forces2[i][j]);
1176 fr->fshift[qm->shiftQM[i]][j] += (fshift[i][j]-fshift2[i][j]);
1181 /* now the last layer still needs to be done: */
1182 qm = qr->qm[qr->nrQMlayers-1]; /* C counts from 0 */
1183 init_QMroutine(cr, qm, mm);
1184 srenew(forces, qm->nrQMatoms);
1185 srenew(fshift, qm->nrQMatoms);
1186 QMener += call_QMroutine(cr, fr, qm, mm, forces, fshift);
1187 for (i = 0; i < qm->nrQMatoms; i++)
1189 for (j = 0; j < DIM; j++)
1191 f[qm->indexQM[i]][j] -= forces[i][j];
1192 fr->fshift[qm->shiftQM[i]][j] += fshift[i][j];
1200 if (qm->bTS || qm->bOPT)
1202 /* qm[0] still contains the largest ONIOM QM subsystem
1203 * we take the optimized coordiates and put the in x[]
1205 for (i = 0; i < qm->nrQMatoms; i++)
1207 for (j = 0; j < DIM; j++)
1209 x[qm->indexQM[i]][j] = qm->xQM[i][j];
1214 } /* calculate_QMMM */
1216 /* end of QMMM core routines */