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,2017, 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_lookup.h"
69 #include "gromacs/topology/mtop_util.h"
70 #include "gromacs/topology/topology.h"
71 #include "gromacs/utility/fatalerror.h"
72 #include "gromacs/utility/smalloc.h"
74 /* declarations of the interfaces to the QM packages. The _SH indicate
75 * the QM interfaces can be used for Surface Hopping simulations
78 /* GAMESS interface */
81 init_gamess(t_commrec
*cr
, t_QMrec
*qm
, t_MMrec
*mm
);
84 call_gamess(t_forcerec
*fr
,
85 t_QMrec
*qm
, t_MMrec
*mm
, rvec f
[], rvec fshift
[]);
91 init_mopac(t_QMrec
*qm
);
94 call_mopac(t_QMrec
*qm
, t_MMrec
*mm
, rvec f
[], rvec fshift
[]);
97 call_mopac_SH(t_QMrec
*qm
, t_MMrec
*mm
, rvec f
[], rvec fshift
[]);
99 #elif GMX_QMMM_GAUSSIAN
100 /* GAUSSIAN interface */
103 init_gaussian(t_QMrec
*qm
);
106 call_gaussian_SH(t_forcerec
*fr
, t_QMrec
*qm
, t_MMrec
*mm
, rvec f
[], rvec fshift
[]);
109 call_gaussian(t_forcerec
*fr
, t_QMrec
*qm
, t_MMrec
*mm
, rvec f
[], rvec fshift
[]);
115 init_orca(t_QMrec
*qm
);
118 call_orca(t_forcerec
*fr
, t_QMrec
*qm
,
119 t_MMrec
*mm
, rvec f
[], rvec fshift
[]);
126 /* this struct and these comparison functions are needed for creating
127 * a QMMM input for the QM routines from the QMMM neighbor list.
135 static int struct_comp(const void *a
, const void *b
)
138 return (int)(((t_j_particle
*)a
)->j
)-(int)(((t_j_particle
*)b
)->j
);
142 real
call_QMroutine(t_commrec gmx_unused
*cr
, t_forcerec gmx_unused
*fr
, t_QMrec gmx_unused
*qm
,
143 t_MMrec gmx_unused
*mm
, rvec gmx_unused f
[], rvec gmx_unused fshift
[])
145 /* makes a call to the requested QM routine (qm->QMmethod)
146 * Note that f is actually the gradient, i.e. -f
151 /* do a semi-empiprical calculation */
153 if (qm
->QMmethod
< eQMmethodRHF
&& !(mm
->nrMMatoms
))
158 QMener
= call_mopac_SH(qm
, mm
, f
, fshift
);
162 QMener
= call_mopac(qm
, mm
, f
, fshift
);
165 gmx_fatal(FARGS
, "Semi-empirical QM only supported with Mopac.");
170 /* do an ab-initio calculation */
171 if (qm
->bSH
&& qm
->QMmethod
== eQMmethodCASSCF
)
173 #if GMX_QMMM_GAUSSIAN
174 QMener
= call_gaussian_SH(fr
, qm
, mm
, f
, fshift
);
176 gmx_fatal(FARGS
, "Ab-initio Surface-hopping only supported with Gaussian.");
182 QMener
= call_gamess(fr
, qm
, mm
, f
, fshift
);
183 #elif GMX_QMMM_GAUSSIAN
184 QMener
= call_gaussian(fr
, qm
, mm
, f
, fshift
);
186 QMener
= call_orca(fr
, qm
, mm
, f
, fshift
);
188 gmx_fatal(FARGS
, "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
195 void init_QMroutine(t_commrec gmx_unused
*cr
, t_QMrec gmx_unused
*qm
, t_MMrec gmx_unused
*mm
)
197 /* makes a call to the requested QM routine (qm->QMmethod)
199 if (qm
->QMmethod
< eQMmethodRHF
)
202 /* do a semi-empiprical calculation */
205 gmx_fatal(FARGS
, "Semi-empirical QM only supported with Mopac.");
210 /* do an ab-initio calculation */
212 init_gamess(cr
, qm
, mm
);
213 #elif GMX_QMMM_GAUSSIAN
218 gmx_fatal(FARGS
, "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
221 } /* init_QMroutine */
223 void update_QMMM_coord(rvec x
[], t_forcerec
*fr
, t_QMrec
*qm
, t_MMrec
*mm
)
225 /* shifts the QM and MM particles into the central box and stores
226 * these shifted coordinates in the coordinate arrays of the
227 * QMMMrec. These coordinates are passed on the QM subroutines.
232 /* shift the QM atoms into the central box
234 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
236 rvec_sub(x
[qm
->indexQM
[i
]], fr
->shift_vec
[qm
->shiftQM
[i
]], qm
->xQM
[i
]);
238 /* also shift the MM atoms into the central box, if any
240 for (i
= 0; i
< mm
->nrMMatoms
; i
++)
242 rvec_sub(x
[mm
->indexMM
[i
]], fr
->shift_vec
[mm
->shiftMM
[i
]], mm
->xMM
[i
]);
244 } /* update_QMMM_coord */
246 static void punch_QMMM_excl(t_QMrec
*qm
, t_MMrec
*mm
, t_blocka
*excls
)
248 /* punch a file containing the bonded interactions of each QM
249 * atom with MM atoms. These need to be excluded in the QM routines
250 * Only needed in case of QM/MM optimizations
255 i
, j
, k
, nrexcl
= 0, *excluded
= nullptr, max_excl
= 0;
258 out
= fopen("QMMMexcl.dat", "w");
260 /* this can be done more efficiently I think
262 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
265 for (j
= excls
->index
[qm
->indexQM
[i
]];
266 j
< excls
->index
[qm
->indexQM
[i
]+1];
269 for (k
= 0; k
< mm
->nrMMatoms
; k
++)
271 if (mm
->indexMM
[k
] == excls
->a
[j
]) /* the excluded MM atom */
273 if (nrexcl
>= max_excl
)
276 srenew(excluded
, max_excl
);
278 excluded
[nrexcl
++] = k
;
284 fprintf(out
, "%5d %5d\n", i
+1, nrexcl
);
285 for (j
= 0; j
< nrexcl
; j
++)
287 fprintf(out
, "%5d ", excluded
[j
]);
293 } /* punch_QMMM_excl */
296 /* end of QMMM subroutines */
298 /* QMMM core routines */
300 t_QMrec
*mk_QMrec(void)
307 t_MMrec
*mk_MMrec(void)
314 static void init_QMrec(int grpnr
, t_QMrec
*qm
, int nr
, int *atomarray
,
315 gmx_mtop_t
*mtop
, t_inputrec
*ir
)
317 /* fills the t_QMrec struct of QM group grpnr
322 snew(qm
->indexQM
, nr
);
323 snew(qm
->shiftQM
, nr
); /* the shifts */
324 for (int i
= 0; i
< nr
; i
++)
326 qm
->indexQM
[i
] = atomarray
[i
];
329 snew(qm
->atomicnumberQM
, nr
);
331 for (int i
= 0; i
< qm
->nrQMatoms
; i
++)
333 const t_atom
&atom
= mtopGetAtomParameters(mtop
, qm
->indexQM
[i
], &molb
);
334 qm
->nelectrons
+= mtop
->atomtypes
.atomnumber
[atom
.type
];
335 qm
->atomicnumberQM
[i
] = mtop
->atomtypes
.atomnumber
[atom
.type
];
338 qm
->QMcharge
= ir
->opts
.QMcharge
[grpnr
];
339 qm
->multiplicity
= ir
->opts
.QMmult
[grpnr
];
340 qm
->nelectrons
-= ir
->opts
.QMcharge
[grpnr
];
342 qm
->QMmethod
= ir
->opts
.QMmethod
[grpnr
];
343 qm
->QMbasis
= ir
->opts
.QMbasis
[grpnr
];
344 /* trajectory surface hopping setup (Gaussian only) */
345 qm
->bSH
= ir
->opts
.bSH
[grpnr
];
346 qm
->CASorbitals
= ir
->opts
.CASorbitals
[grpnr
];
347 qm
->CASelectrons
= ir
->opts
.CASelectrons
[grpnr
];
348 qm
->SAsteps
= ir
->opts
.SAsteps
[grpnr
];
349 qm
->SAon
= ir
->opts
.SAon
[grpnr
];
350 qm
->SAoff
= ir
->opts
.SAoff
[grpnr
];
351 /* hack to prevent gaussian from reinitializing all the time */
352 qm
->nQMcpus
= 0; /* number of CPU's to be used by g01, is set
353 * upon initializing gaussian
356 /* print the current layer to allow users to check their input */
357 fprintf(stderr
, "Layer %d\nnr of QM atoms %d\n", grpnr
, nr
);
358 fprintf(stderr
, "QMlevel: %s/%s\n\n",
359 eQMmethod_names
[qm
->QMmethod
], eQMbasis_names
[qm
->QMbasis
]);
362 snew(qm
->frontatoms
, nr
);
363 /* Lennard-Jones coefficients */
366 /* do we optimize the QM separately using the algorithms of the QM program??
368 qm
->bTS
= ir
->opts
.bTS
[grpnr
];
369 qm
->bOPT
= ir
->opts
.bOPT
[grpnr
];
373 t_QMrec
*copy_QMrec(t_QMrec
*qm
)
375 /* copies the contents of qm into a new t_QMrec struct */
382 qmcopy
->nrQMatoms
= qm
->nrQMatoms
;
383 snew(qmcopy
->xQM
, qmcopy
->nrQMatoms
);
384 snew(qmcopy
->indexQM
, qmcopy
->nrQMatoms
);
385 snew(qmcopy
->atomicnumberQM
, qm
->nrQMatoms
);
386 snew(qmcopy
->shiftQM
, qmcopy
->nrQMatoms
); /* the shifts */
387 for (i
= 0; i
< qmcopy
->nrQMatoms
; i
++)
389 qmcopy
->shiftQM
[i
] = qm
->shiftQM
[i
];
390 qmcopy
->indexQM
[i
] = qm
->indexQM
[i
];
391 qmcopy
->atomicnumberQM
[i
] = qm
->atomicnumberQM
[i
];
393 qmcopy
->nelectrons
= qm
->nelectrons
;
394 qmcopy
->multiplicity
= qm
->multiplicity
;
395 qmcopy
->QMcharge
= qm
->QMcharge
;
396 qmcopy
->nelectrons
= qm
->nelectrons
;
397 qmcopy
->QMmethod
= qm
->QMmethod
;
398 qmcopy
->QMbasis
= qm
->QMbasis
;
399 /* trajectory surface hopping setup (Gaussian only) */
400 qmcopy
->bSH
= qm
->bSH
;
401 qmcopy
->CASorbitals
= qm
->CASorbitals
;
402 qmcopy
->CASelectrons
= qm
->CASelectrons
;
403 qmcopy
->SAsteps
= qm
->SAsteps
;
404 qmcopy
->SAon
= qm
->SAon
;
405 qmcopy
->SAoff
= qm
->SAoff
;
406 qmcopy
->bOPT
= qm
->bOPT
;
408 /* Gaussian init. variables */
409 qmcopy
->nQMcpus
= qm
->nQMcpus
;
410 for (i
= 0; i
< DIM
; i
++)
412 qmcopy
->SHbasis
[i
] = qm
->SHbasis
[i
];
414 qmcopy
->QMmem
= qm
->QMmem
;
415 qmcopy
->accuracy
= qm
->accuracy
;
416 qmcopy
->cpmcscf
= qm
->cpmcscf
;
417 qmcopy
->SAstep
= qm
->SAstep
;
418 snew(qmcopy
->frontatoms
, qm
->nrQMatoms
);
419 snew(qmcopy
->c12
, qmcopy
->nrQMatoms
);
420 snew(qmcopy
->c6
, qmcopy
->nrQMatoms
);
421 if (qmcopy
->bTS
|| qmcopy
->bOPT
)
423 for (i
= 1; i
< qmcopy
->nrQMatoms
; i
++)
425 qmcopy
->frontatoms
[i
] = qm
->frontatoms
[i
];
426 qmcopy
->c12
[i
] = qm
->c12
[i
];
427 qmcopy
->c6
[i
] = qm
->c6
[i
];
435 t_QMMMrec
*mk_QMMMrec(void)
446 void init_QMMMrec(t_commrec
*cr
,
451 /* we put the atomsnumbers of atoms that belong to the QMMM group in
452 * an array that will be copied later to QMMMrec->indexQM[..]. Also
453 * it will be used to create an QMMMrec->bQMMM index array that
454 * simply contains true/false for QM and MM (the other) atoms.
457 gmx_groups_t
*groups
;
458 int *qm_arr
= nullptr, vsite
, ai
, aj
;
459 int qm_max
= 0, qm_nr
= 0, i
, j
, jmax
, k
, l
, nrvsite2
= 0;
464 gmx_mtop_atomloop_all_t aloop
;
465 gmx_mtop_ilistloop_all_t iloop
;
469 if (ir
->cutoff_scheme
!= ecutsGROUP
)
471 gmx_fatal(FARGS
, "QMMM is currently only supported with cutoff-scheme=group");
473 if (!EI_DYNAMICS(ir
->eI
))
475 gmx_fatal(FARGS
, "QMMM is only supported with dynamics");
478 c6au
= (HARTREE2KJ
*AVOGADRO
*gmx::power6(BOHR2NM
));
479 c12au
= (HARTREE2KJ
*AVOGADRO
*gmx::power12(BOHR2NM
));
480 /* issue a fatal if the user wants to run with more than one node */
483 gmx_fatal(FARGS
, "QM/MM does not work in parallel, use a single rank instead\n");
486 /* Make a local copy of the QMMMrec */
489 /* bQMMM[..] is an array containing TRUE/FALSE for atoms that are
490 * QM/not QM. We first set all elemenst at false. Afterwards we use
491 * the qm_arr (=MMrec->indexQM) to changes the elements
492 * corresponding to the QM atoms at TRUE. */
494 qr
->QMMMscheme
= ir
->QMMMscheme
;
496 /* we take the possibility into account that a user has
497 * defined more than one QM group:
499 /* an ugly work-around in case there is only one group In this case
500 * the whole system is treated as QM. Otherwise the second group is
501 * always the rest of the total system and is treated as MM.
504 /* small problem if there is only QM.... so no MM */
506 jmax
= ir
->opts
.ngQM
;
508 if (qr
->QMMMscheme
== eQMMMschemeoniom
)
510 qr
->nrQMlayers
= jmax
;
517 groups
= &mtop
->groups
;
519 /* there are jmax groups of QM atoms. In case of multiple QM groups
520 * I assume that the users wants to do ONIOM. However, maybe it
521 * should also be possible to define more than one QM subsystem with
522 * independent neighbourlists. I have to think about
526 for (j
= 0; j
< jmax
; j
++)
529 aloop
= gmx_mtop_atomloop_all_init(mtop
);
531 while (gmx_mtop_atomloop_all_next(aloop
, &i
, &atom
))
536 srenew(qm_arr
, qm_max
);
538 if (ggrpnr(groups
, egcQMMM
, i
) == j
)
544 if (qr
->QMMMscheme
== eQMMMschemeoniom
)
546 /* add the atoms to the bQMMM array
549 /* I assume that users specify the QM groups from small to
550 * big(ger) in the mdp file
552 qr
->qm
[j
] = mk_QMrec();
553 /* we need to throw out link atoms that in the previous layer
554 * existed to separate this QMlayer from the previous
555 * QMlayer. We use the iatoms array in the idef for that
556 * purpose. If all atoms defining the current Link Atom (Dummy2)
557 * are part of the current QM layer it needs to be removed from
560 iloop
= gmx_mtop_ilistloop_all_init(mtop
);
561 while (gmx_mtop_ilistloop_all_next(iloop
, &ilist_mol
, &a_offset
))
563 nrvsite2
= ilist_mol
[F_VSITE2
].nr
;
564 iatoms
= ilist_mol
[F_VSITE2
].iatoms
;
566 for (k
= 0; k
< nrvsite2
; k
+= 4)
568 vsite
= a_offset
+ iatoms
[k
+1]; /* the vsite */
569 ai
= a_offset
+ iatoms
[k
+2]; /* constructing atom */
570 aj
= a_offset
+ iatoms
[k
+3]; /* constructing atom */
571 if (ggrpnr(groups
, egcQMMM
, vsite
) == ggrpnr(groups
, egcQMMM
, ai
)
573 ggrpnr(groups
, egcQMMM
, vsite
) == ggrpnr(groups
, egcQMMM
, aj
))
575 /* this dummy link atom needs to be removed from the qm_arr
576 * before making the QMrec of this layer!
578 for (i
= 0; i
< qm_nr
; i
++)
580 if (qm_arr
[i
] == vsite
)
582 /* drop the element */
583 for (l
= i
; l
< qm_nr
; l
++)
585 qm_arr
[l
] = qm_arr
[l
+1];
594 /* store QM atoms in this layer in the QMrec and initialise layer
596 init_QMrec(j
, qr
->qm
[j
], qm_nr
, qm_arr
, mtop
, ir
);
598 /* we now store the LJ C6 and C12 parameters in QM rec in case
599 * we need to do an optimization
601 if (qr
->qm
[j
]->bOPT
|| qr
->qm
[j
]->bTS
)
603 for (i
= 0; i
< qm_nr
; i
++)
605 /* nbfp now includes the 6.0/12.0 derivative prefactors */
606 qr
->qm
[j
]->c6
[i
] = C6(fr
->nbfp
, mtop
->ffparams
.atnr
, atom
->type
, atom
->type
)/c6au
/6.0;
607 qr
->qm
[j
]->c12
[i
] = C12(fr
->nbfp
, mtop
->ffparams
.atnr
, atom
->type
, atom
->type
)/c12au
/12.0;
610 /* now we check for frontier QM atoms. These occur in pairs that
611 * construct the vsite
613 iloop
= gmx_mtop_ilistloop_all_init(mtop
);
614 while (gmx_mtop_ilistloop_all_next(iloop
, &ilist_mol
, &a_offset
))
616 nrvsite2
= ilist_mol
[F_VSITE2
].nr
;
617 iatoms
= ilist_mol
[F_VSITE2
].iatoms
;
619 for (k
= 0; k
< nrvsite2
; k
+= 4)
621 vsite
= a_offset
+ iatoms
[k
+1]; /* the vsite */
622 ai
= a_offset
+ iatoms
[k
+2]; /* constructing atom */
623 aj
= a_offset
+ iatoms
[k
+3]; /* constructing atom */
624 if (ggrpnr(groups
, egcQMMM
, ai
) < (groups
->grps
[egcQMMM
].nr
-1) &&
625 (ggrpnr(groups
, egcQMMM
, aj
) >= (groups
->grps
[egcQMMM
].nr
-1)))
627 /* mark ai as frontier atom */
628 for (i
= 0; i
< qm_nr
; i
++)
630 if ( (qm_arr
[i
] == ai
) || (qm_arr
[i
] == vsite
) )
632 qr
->qm
[j
]->frontatoms
[i
] = TRUE
;
636 else if (ggrpnr(groups
, egcQMMM
, aj
) < (groups
->grps
[egcQMMM
].nr
-1) &&
637 (ggrpnr(groups
, egcQMMM
, ai
) >= (groups
->grps
[egcQMMM
].nr
-1)))
639 /* mark aj as frontier atom */
640 for (i
= 0; i
< qm_nr
; i
++)
642 if ( (qm_arr
[i
] == aj
) || (qm_arr
[i
] == vsite
))
644 qr
->qm
[j
]->frontatoms
[i
] = TRUE
;
652 if (qr
->QMMMscheme
!= eQMMMschemeoniom
)
655 /* standard QMMM, all layers are merged together so there is one QM
656 * subsystem and one MM subsystem.
657 * Also we set the charges to zero in mtop to prevent the innerloops
658 * from doubly counting the electostatic QM MM interaction
659 * TODO: Consider doing this in grompp instead.
663 for (k
= 0; k
< qm_nr
; k
++)
666 mtopGetMolblockIndex(mtop
, qm_arr
[k
], &molb
, nullptr, &indexInMolecule
);
667 t_atom
*atom
= &mtop
->moltype
[mtop
->molblock
[molb
].type
].atoms
.atom
[indexInMolecule
];
671 qr
->qm
[0] = mk_QMrec();
672 /* store QM atoms in the QMrec and initialise
674 init_QMrec(0, qr
->qm
[0], qm_nr
, qm_arr
, mtop
, ir
);
675 if (qr
->qm
[0]->bOPT
|| qr
->qm
[0]->bTS
)
677 for (i
= 0; i
< qm_nr
; i
++)
679 const t_atom
&atom
= mtopGetAtomParameters(mtop
, qm_arr
[i
], &molb
);
680 /* nbfp now includes the 6.0/12.0 derivative prefactors */
681 qr
->qm
[0]->c6
[i
] = C6(fr
->nbfp
, mtop
->ffparams
.atnr
, atom
.type
, atom
.type
)/c6au
/6.0;
682 qr
->qm
[0]->c12
[i
] = C12(fr
->nbfp
, mtop
->ffparams
.atnr
, atom
.type
, atom
.type
)/c12au
/12.0;
686 /* find frontier atoms and mark them true in the frontieratoms array.
688 for (i
= 0; i
< qm_nr
; i
++)
690 mtopGetMolblockIndex(mtop
, qm_arr
[i
], &molb
, nullptr, &a_offset
);
691 ilist_mol
= mtop
->moltype
[mtop
->molblock
[molb
].type
].ilist
;
692 nrvsite2
= ilist_mol
[F_VSITE2
].nr
;
693 iatoms
= ilist_mol
[F_VSITE2
].iatoms
;
695 for (k
= 0; k
< nrvsite2
; k
+= 4)
697 vsite
= a_offset
+ iatoms
[k
+1]; /* the vsite */
698 ai
= a_offset
+ iatoms
[k
+2]; /* constructing atom */
699 aj
= a_offset
+ iatoms
[k
+3]; /* constructing atom */
700 if (ggrpnr(groups
, egcQMMM
, ai
) < (groups
->grps
[egcQMMM
].nr
-1) &&
701 (ggrpnr(groups
, egcQMMM
, aj
) >= (groups
->grps
[egcQMMM
].nr
-1)))
703 /* mark ai as frontier atom */
704 if ( (qm_arr
[i
] == ai
) || (qm_arr
[i
] == vsite
) )
706 qr
->qm
[0]->frontatoms
[i
] = TRUE
;
709 else if (ggrpnr(groups
, egcQMMM
, aj
) < (groups
->grps
[egcQMMM
].nr
-1) &&
710 (ggrpnr(groups
, egcQMMM
, ai
) >= (groups
->grps
[egcQMMM
].nr
-1)))
712 /* mark aj as frontier atom */
713 if ( (qm_arr
[i
] == aj
) || (qm_arr
[i
] == vsite
) )
715 qr
->qm
[0]->frontatoms
[i
] = TRUE
;
721 /* MM rec creation */
723 mm
->scalefactor
= ir
->scalefactor
;
724 mm
->nrMMatoms
= (mtop
->natoms
)-(qr
->qm
[0]->nrQMatoms
); /* rest of the atoms */
728 { /* MM rec creation */
730 mm
->scalefactor
= ir
->scalefactor
;
735 /* these variables get updated in the update QMMMrec */
737 if (qr
->nrQMlayers
== 1)
739 /* with only one layer there is only one initialisation
740 * needed. Multilayer is a bit more complicated as it requires
741 * re-initialisation at every step of the simulation. This is due
742 * to the use of COMMON blocks in the fortran QM subroutines.
744 if (qr
->qm
[0]->QMmethod
< eQMmethodRHF
)
747 /* semi-empiprical 1-layer ONIOM calculation requested (mopac93) */
748 init_mopac(qr
->qm
[0]);
750 gmx_fatal(FARGS
, "Semi-empirical QM only supported with Mopac.");
755 /* ab initio calculation requested (gamess/gaussian/ORCA) */
757 init_gamess(cr
, qr
->qm
[0], qr
->mm
);
758 #elif GMX_QMMM_GAUSSIAN
759 init_gaussian(qr
->qm
[0]);
761 init_orca(qr
->qm
[0]);
763 gmx_fatal(FARGS
, "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
769 void update_QMMMrec(t_commrec
*cr
,
776 /* updates the coordinates of both QM atoms and MM atoms and stores
777 * them in the QMMMrec.
779 * NOTE: is NOT yet working if there are no PBC. Also in ns.c, simple
780 * ns needs to be fixed!
783 mm_max
= 0, mm_nr
= 0, mm_nr_new
, i
, j
, is
, k
, shift
;
785 *mm_j_particles
= nullptr, *qm_i_particles
= nullptr;
799 *parallelMMarray
= nullptr;
803 c6au
= (HARTREE2KJ
*AVOGADRO
*gmx::power6(BOHR2NM
));
804 c12au
= (HARTREE2KJ
*AVOGADRO
*gmx::power12(BOHR2NM
));
806 /* every cpu has this array. On every processor we fill this array
807 * with 1's and 0's. 1's indicate the atoms is a QM atom on the
808 * current cpu in a later stage these arrays are all summed. indexes
809 * > 0 indicate the atom is a QM atom. Every node therefore knows
810 * whcih atoms are part of the QM subsystem.
812 /* copy some pointers */
815 QMMMlist
= fr
->QMMMlist
;
819 /* init_pbc(box); needs to be called first, see pbc.h */
821 clear_ivec(null_ivec
);
822 set_pbc_dd(&pbc
, fr
->ePBC
, DOMAINDECOMP(cr
) ? cr
->dd
->nc
: null_ivec
,
824 /* only in standard (normal) QMMM we need the neighbouring MM
825 * particles to provide a electric field of point charges for the QM
828 if (qr
->QMMMscheme
== eQMMMschemenormal
) /* also implies 1 QM-layer */
830 /* we NOW create/update a number of QMMMrec entries:
832 * 1) the shiftQM, containing the shifts of the QM atoms
834 * 2) the indexMM array, containing the index of the MM atoms
836 * 3) the shiftMM, containing the shifts of the MM atoms
838 * 4) the shifted coordinates of the MM atoms
840 * the shifts are used for computing virial of the QM/MM particles.
842 qm
= qr
->qm
[0]; /* in case of normal QMMM, there is only one group */
843 snew(qm_i_particles
, QMMMlist
->nri
);
846 qm_i_particles
[0].shift
= XYZ2IS(0, 0, 0);
847 for (i
= 0; i
< QMMMlist
->nri
; i
++)
849 qm_i_particles
[i
].j
= QMMMlist
->iinr
[i
];
853 qm_i_particles
[i
].shift
= pbc_dx_aiuc(&pbc
, x
[QMMMlist
->iinr
[0]],
854 x
[QMMMlist
->iinr
[i
]], dx
);
857 /* However, since nri >= nrQMatoms, we do a quicksort, and throw
858 * out double, triple, etc. entries later, as we do for the MM
862 /* compute the shift for the MM j-particles with respect to
863 * the QM i-particle and store them.
866 crd
[0] = IS2X(QMMMlist
->shift
[i
]) + IS2X(qm_i_particles
[i
].shift
);
867 crd
[1] = IS2Y(QMMMlist
->shift
[i
]) + IS2Y(qm_i_particles
[i
].shift
);
868 crd
[2] = IS2Z(QMMMlist
->shift
[i
]) + IS2Z(qm_i_particles
[i
].shift
);
869 is
= static_cast<int>(XYZ2IS(crd
[0], crd
[1], crd
[2]));
870 for (j
= QMMMlist
->jindex
[i
];
871 j
< QMMMlist
->jindex
[i
+1];
877 srenew(mm_j_particles
, mm_max
);
880 mm_j_particles
[mm_nr
].j
= QMMMlist
->jjnr
[j
];
881 mm_j_particles
[mm_nr
].shift
= is
;
886 /* quicksort QM and MM shift arrays and throw away multiple entries */
890 qsort(qm_i_particles
, QMMMlist
->nri
,
891 (size_t)sizeof(qm_i_particles
[0]),
893 /* The mm_j_particles argument to qsort is not allowed to be NULL */
896 qsort(mm_j_particles
, mm_nr
,
897 (size_t)sizeof(mm_j_particles
[0]),
900 /* remove multiples in the QM shift array, since in init_QMMM() we
901 * went through the atom numbers from 0 to md.nr, the order sorted
902 * here matches the one of QMindex already.
905 for (i
= 0; i
< QMMMlist
->nri
; i
++)
907 if (i
== 0 || qm_i_particles
[i
].j
!= qm_i_particles
[i
-1].j
)
909 qm_i_particles
[j
++] = qm_i_particles
[i
];
913 if (qm
->bTS
|| qm
->bOPT
)
915 /* only remove double entries for the MM array */
916 for (i
= 0; i
< mm_nr
; i
++)
918 if ((i
== 0 || mm_j_particles
[i
].j
!= mm_j_particles
[i
-1].j
)
919 && !md
->bQM
[mm_j_particles
[i
].j
])
921 mm_j_particles
[mm_nr_new
++] = mm_j_particles
[i
];
925 /* we also remove mm atoms that have no charges!
926 * actually this is already done in the ns.c
930 for (i
= 0; i
< mm_nr
; i
++)
932 if ((i
== 0 || mm_j_particles
[i
].j
!= mm_j_particles
[i
-1].j
)
933 && !md
->bQM
[mm_j_particles
[i
].j
]
934 && (md
->chargeA
[mm_j_particles
[i
].j
]
935 || (md
->chargeB
&& md
->chargeB
[mm_j_particles
[i
].j
])))
937 mm_j_particles
[mm_nr_new
++] = mm_j_particles
[i
];
942 /* store the data retrieved above into the QMMMrec
945 /* Keep the compiler happy,
946 * shift will always be set in the loop for i=0
949 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
951 /* not all qm particles might have appeared as i
952 * particles. They might have been part of the same charge
953 * group for instance.
955 if (qm
->indexQM
[i
] == qm_i_particles
[k
].j
)
957 shift
= qm_i_particles
[k
++].shift
;
959 /* use previous shift, assuming they belong the same charge
963 qm
->shiftQM
[i
] = shift
;
966 /* parallel excecution */
969 snew(parallelMMarray
, 2*(md
->nr
));
970 /* only MM particles have a 1 at their atomnumber. The second part
971 * of the array contains the shifts. Thus:
972 * p[i]=1/0 depending on wether atomnumber i is a MM particle in the QM
973 * step or not. p[i+md->nr] is the shift of atomnumber i.
975 for (i
= 0; i
< 2*(md
->nr
); i
++)
977 parallelMMarray
[i
] = 0;
980 for (i
= 0; i
< mm_nr
; i
++)
982 parallelMMarray
[mm_j_particles
[i
].j
] = 1;
983 parallelMMarray
[mm_j_particles
[i
].j
+(md
->nr
)] = mm_j_particles
[i
].shift
;
985 gmx_sumi(md
->nr
, parallelMMarray
, cr
);
989 for (i
= 0; i
< md
->nr
; i
++)
991 if (parallelMMarray
[i
])
996 srenew(mm
->indexMM
, mm_max
);
997 srenew(mm
->shiftMM
, mm_max
);
999 mm
->indexMM
[mm_nr
] = i
;
1000 mm
->shiftMM
[mm_nr
++] = parallelMMarray
[i
+md
->nr
]/parallelMMarray
[i
];
1003 mm
->nrMMatoms
= mm_nr
;
1004 free(parallelMMarray
);
1006 /* serial execution */
1009 mm
->nrMMatoms
= mm_nr
;
1010 srenew(mm
->shiftMM
, mm_nr
);
1011 srenew(mm
->indexMM
, mm_nr
);
1012 for (i
= 0; i
< mm_nr
; i
++)
1014 mm
->indexMM
[i
] = mm_j_particles
[i
].j
;
1015 mm
->shiftMM
[i
] = mm_j_particles
[i
].shift
;
1019 /* (re) allocate memory for the MM coordiate array. The QM
1020 * coordinate array was already allocated in init_QMMM, and is
1021 * only (re)filled in the update_QMMM_coordinates routine
1023 srenew(mm
->xMM
, mm
->nrMMatoms
);
1024 /* now we (re) fill the array that contains the MM charges with
1025 * the forcefield charges. If requested, these charges will be
1026 * scaled by a factor
1028 srenew(mm
->MMcharges
, mm
->nrMMatoms
);
1029 for (i
= 0; i
< mm
->nrMMatoms
; i
++) /* no free energy yet */
1031 mm
->MMcharges
[i
] = md
->chargeA
[mm
->indexMM
[i
]]*mm
->scalefactor
;
1033 if (qm
->bTS
|| qm
->bOPT
)
1035 /* store (copy) the c6 and c12 parameters into the MMrec struct
1037 srenew(mm
->c6
, mm
->nrMMatoms
);
1038 srenew(mm
->c12
, mm
->nrMMatoms
);
1039 for (i
= 0; i
< mm
->nrMMatoms
; i
++)
1041 /* nbfp now includes the 6.0/12.0 derivative prefactors */
1042 mm
->c6
[i
] = C6(fr
->nbfp
, top
->idef
.atnr
, md
->typeA
[mm
->indexMM
[i
]], md
->typeA
[mm
->indexMM
[i
]])/c6au
/6.0;
1043 mm
->c12
[i
] = C12(fr
->nbfp
, top
->idef
.atnr
, md
->typeA
[mm
->indexMM
[i
]], md
->typeA
[mm
->indexMM
[i
]])/c12au
/12.0;
1045 punch_QMMM_excl(qr
->qm
[0], mm
, &(top
->excls
));
1047 /* the next routine fills the coordinate fields in the QMMM rec of
1048 * both the qunatum atoms and the MM atoms, using the shifts
1052 update_QMMM_coord(x
, fr
, qr
->qm
[0], qr
->mm
);
1053 free(qm_i_particles
);
1054 free(mm_j_particles
);
1056 else /* ONIOM */ /* ????? */
1059 /* do for each layer */
1060 for (j
= 0; j
< qr
->nrQMlayers
; j
++)
1063 qm
->shiftQM
[0] = XYZ2IS(0, 0, 0);
1064 for (i
= 1; i
< qm
->nrQMatoms
; i
++)
1066 qm
->shiftQM
[i
] = pbc_dx_aiuc(&pbc
, x
[qm
->indexQM
[0]], x
[qm
->indexQM
[i
]],
1069 update_QMMM_coord(x
, fr
, qm
, mm
);
1072 } /* update_QMMM_rec */
1075 real
calculate_QMMM(t_commrec
*cr
,
1081 /* a selection for the QM package depending on which is requested
1082 * (Gaussian, GAMESS-UK, MOPAC or ORCA) needs to be implemented here. Now
1083 * it works through defines.... Not so nice yet
1092 *forces
= nullptr, *fshift
= nullptr,
1093 *forces2
= nullptr, *fshift2
= nullptr; /* needed for multilayer ONIOM */
1096 /* make a local copy the QMMMrec pointer
1101 /* now different procedures are carried out for one layer ONION and
1102 * normal QMMM on one hand and multilayer oniom on the other
1104 if (qr
->QMMMscheme
== eQMMMschemenormal
|| qr
->nrQMlayers
== 1)
1107 snew(forces
, (qm
->nrQMatoms
+mm
->nrMMatoms
));
1108 snew(fshift
, (qm
->nrQMatoms
+mm
->nrMMatoms
));
1109 QMener
= call_QMroutine(cr
, fr
, qm
, mm
, forces
, fshift
);
1110 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
1112 for (j
= 0; j
< DIM
; j
++)
1114 f
[qm
->indexQM
[i
]][j
] -= forces
[i
][j
];
1115 fr
->fshift
[qm
->shiftQM
[i
]][j
] += fshift
[i
][j
];
1118 for (i
= 0; i
< mm
->nrMMatoms
; i
++)
1120 for (j
= 0; j
< DIM
; j
++)
1122 f
[mm
->indexMM
[i
]][j
] -= forces
[qm
->nrQMatoms
+i
][j
];
1123 fr
->fshift
[mm
->shiftMM
[i
]][j
] += fshift
[qm
->nrQMatoms
+i
][j
];
1130 else /* Multi-layer ONIOM */
1132 for (i
= 0; i
< qr
->nrQMlayers
-1; i
++) /* last layer is special */
1135 qm2
= copy_QMrec(qr
->qm
[i
+1]);
1137 qm2
->nrQMatoms
= qm
->nrQMatoms
;
1139 for (j
= 0; j
< qm2
->nrQMatoms
; j
++)
1141 for (k
= 0; k
< DIM
; k
++)
1143 qm2
->xQM
[j
][k
] = qm
->xQM
[j
][k
];
1145 qm2
->indexQM
[j
] = qm
->indexQM
[j
];
1146 qm2
->atomicnumberQM
[j
] = qm
->atomicnumberQM
[j
];
1147 qm2
->shiftQM
[j
] = qm
->shiftQM
[j
];
1150 qm2
->QMcharge
= qm
->QMcharge
;
1151 /* this layer at the higher level of theory */
1152 srenew(forces
, qm
->nrQMatoms
);
1153 srenew(fshift
, qm
->nrQMatoms
);
1154 /* we need to re-initialize the QMroutine every step... */
1155 init_QMroutine(cr
, qm
, mm
);
1156 QMener
+= call_QMroutine(cr
, fr
, qm
, mm
, forces
, fshift
);
1158 /* this layer at the lower level of theory */
1159 srenew(forces2
, qm
->nrQMatoms
);
1160 srenew(fshift2
, qm
->nrQMatoms
);
1161 init_QMroutine(cr
, qm2
, mm
);
1162 QMener
-= call_QMroutine(cr
, fr
, qm2
, mm
, forces2
, fshift2
);
1163 /* E = E1high-E1low The next layer includes the current layer at
1164 * the lower level of theory, which provides + E2low
1165 * this is similar for gradients
1167 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
1169 for (j
= 0; j
< DIM
; j
++)
1171 f
[qm
->indexQM
[i
]][j
] -= (forces
[i
][j
]-forces2
[i
][j
]);
1172 fr
->fshift
[qm
->shiftQM
[i
]][j
] += (fshift
[i
][j
]-fshift2
[i
][j
]);
1177 /* now the last layer still needs to be done: */
1178 qm
= qr
->qm
[qr
->nrQMlayers
-1]; /* C counts from 0 */
1179 init_QMroutine(cr
, qm
, mm
);
1180 srenew(forces
, qm
->nrQMatoms
);
1181 srenew(fshift
, qm
->nrQMatoms
);
1182 QMener
+= call_QMroutine(cr
, fr
, qm
, mm
, forces
, fshift
);
1183 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
1185 for (j
= 0; j
< DIM
; j
++)
1187 f
[qm
->indexQM
[i
]][j
] -= forces
[i
][j
];
1188 fr
->fshift
[qm
->shiftQM
[i
]][j
] += fshift
[i
][j
];
1196 if (qm
->bTS
|| qm
->bOPT
)
1198 /* qm[0] still contains the largest ONIOM QM subsystem
1199 * we take the optimized coordiates and put the in x[]
1201 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
1203 for (j
= 0; j
< DIM
; j
++)
1205 x
[qm
->indexQM
[i
]][j
] = qm
->xQM
[i
][j
];
1210 } /* calculate_QMMM */
1212 /* end of QMMM core routines */