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, 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/fileio/confio.h"
53 #include "gromacs/legacyheaders/names.h"
54 #include "gromacs/legacyheaders/network.h"
55 #include "gromacs/legacyheaders/nrnb.h"
56 #include "gromacs/legacyheaders/txtdump.h"
57 #include "gromacs/legacyheaders/typedefs.h"
58 #include "gromacs/legacyheaders/types/commrec.h"
59 #include "gromacs/legacyheaders/types/mdatom.h"
60 #include "gromacs/math/units.h"
61 #include "gromacs/math/vec.h"
62 #include "gromacs/mdlib/force.h"
63 #include "gromacs/mdlib/ns.h"
64 #include "gromacs/pbcutil/ishift.h"
65 #include "gromacs/pbcutil/pbc.h"
66 #include "gromacs/topology/mtop_util.h"
67 #include "gromacs/utility/fatalerror.h"
68 #include "gromacs/utility/smalloc.h"
70 /* declarations of the interfaces to the QM packages. The _SH indicate
71 * the QM interfaces can be used for Surface Hopping simulations
74 /* GAMESS interface */
77 init_gamess(t_commrec
*cr
, t_QMrec
*qm
, t_MMrec
*mm
);
80 call_gamess(t_forcerec
*fr
,
81 t_QMrec
*qm
, t_MMrec
*mm
, rvec f
[], rvec fshift
[]);
87 init_mopac(t_QMrec
*qm
);
90 call_mopac(t_QMrec
*qm
, t_MMrec
*mm
, rvec f
[], rvec fshift
[]);
93 call_mopac_SH(t_QMrec
*qm
, t_MMrec
*mm
, rvec f
[], rvec fshift
[]);
95 #elif GMX_QMMM_GAUSSIAN
96 /* GAUSSIAN interface */
99 init_gaussian(t_QMrec
*qm
);
102 call_gaussian_SH(t_forcerec
*fr
, t_QMrec
*qm
, t_MMrec
*mm
, rvec f
[], rvec fshift
[]);
105 call_gaussian(t_forcerec
*fr
, t_QMrec
*qm
, t_MMrec
*mm
, rvec f
[], rvec fshift
[]);
111 init_orca(t_QMrec
*qm
);
114 call_orca(t_forcerec
*fr
, t_QMrec
*qm
,
115 t_MMrec
*mm
, rvec f
[], rvec fshift
[]);
122 /* this struct and these comparison functions are needed for creating
123 * a QMMM input for the QM routines from the QMMM neighbor list.
131 static int struct_comp(const void *a
, const void *b
)
134 return (int)(((t_j_particle
*)a
)->j
)-(int)(((t_j_particle
*)b
)->j
);
138 static int int_comp(const void *a
, const void *b
)
141 return (*(int *)a
) - (*(int *)b
);
145 static int QMlayer_comp(const void *a
, const void *b
)
148 return (int)(((t_QMrec
*)a
)->nrQMatoms
)-(int)(((t_QMrec
*)b
)->nrQMatoms
);
152 real
call_QMroutine(t_commrec gmx_unused
*cr
, t_forcerec gmx_unused
*fr
, t_QMrec gmx_unused
*qm
,
153 t_MMrec gmx_unused
*mm
, rvec gmx_unused f
[], rvec gmx_unused fshift
[])
155 /* makes a call to the requested QM routine (qm->QMmethod)
156 * Note that f is actually the gradient, i.e. -f
161 /* do a semi-empiprical calculation */
163 if (qm
->QMmethod
< eQMmethodRHF
&& !(mm
->nrMMatoms
))
168 QMener
= call_mopac_SH(qm
, mm
, f
, fshift
);
172 QMener
= call_mopac(qm
, mm
, f
, fshift
);
175 gmx_fatal(FARGS
, "Semi-empirical QM only supported with Mopac.");
180 /* do an ab-initio calculation */
181 if (qm
->bSH
&& qm
->QMmethod
== eQMmethodCASSCF
)
183 #if GMX_QMMM_GAUSSIAN
184 QMener
= call_gaussian_SH(fr
, qm
, mm
, f
, fshift
);
186 gmx_fatal(FARGS
, "Ab-initio Surface-hopping only supported with Gaussian.");
192 QMener
= call_gamess(fr
, qm
, mm
, f
, fshift
);
193 #elif GMX_QMMM_GAUSSIAN
194 QMener
= call_gaussian(fr
, qm
, mm
, f
, fshift
);
196 QMener
= call_orca(fr
, qm
, mm
, f
, fshift
);
198 gmx_fatal(FARGS
, "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
205 void init_QMroutine(t_commrec gmx_unused
*cr
, t_QMrec gmx_unused
*qm
, t_MMrec gmx_unused
*mm
)
207 /* makes a call to the requested QM routine (qm->QMmethod)
209 if (qm
->QMmethod
< eQMmethodRHF
)
212 /* do a semi-empiprical calculation */
215 gmx_fatal(FARGS
, "Semi-empirical QM only supported with Mopac.");
220 /* do an ab-initio calculation */
222 init_gamess(cr
, qm
, mm
);
223 #elif GMX_QMMM_GAUSSIAN
228 gmx_fatal(FARGS
, "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
231 } /* init_QMroutine */
233 void update_QMMM_coord(rvec x
[], t_forcerec
*fr
, t_QMrec
*qm
, t_MMrec
*mm
)
235 /* shifts the QM and MM particles into the central box and stores
236 * these shifted coordinates in the coordinate arrays of the
237 * QMMMrec. These coordinates are passed on the QM subroutines.
242 /* shift the QM atoms into the central box
244 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
246 rvec_sub(x
[qm
->indexQM
[i
]], fr
->shift_vec
[qm
->shiftQM
[i
]], qm
->xQM
[i
]);
248 /* also shift the MM atoms into the central box, if any
250 for (i
= 0; i
< mm
->nrMMatoms
; i
++)
252 rvec_sub(x
[mm
->indexMM
[i
]], fr
->shift_vec
[mm
->shiftMM
[i
]], mm
->xMM
[i
]);
254 } /* update_QMMM_coord */
256 static void punch_QMMM_excl(t_QMrec
*qm
, t_MMrec
*mm
, t_blocka
*excls
)
258 /* punch a file containing the bonded interactions of each QM
259 * atom with MM atoms. These need to be excluded in the QM routines
260 * Only needed in case of QM/MM optimizations
265 i
, j
, k
, nrexcl
= 0, *excluded
= NULL
, max_excl
= 0;
268 out
= fopen("QMMMexcl.dat", "w");
270 /* this can be done more efficiently I think
272 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
275 for (j
= excls
->index
[qm
->indexQM
[i
]];
276 j
< excls
->index
[qm
->indexQM
[i
]+1];
279 for (k
= 0; k
< mm
->nrMMatoms
; k
++)
281 if (mm
->indexMM
[k
] == excls
->a
[j
]) /* the excluded MM atom */
283 if (nrexcl
>= max_excl
)
286 srenew(excluded
, max_excl
);
288 excluded
[nrexcl
++] = k
;
294 fprintf(out
, "%5d %5d\n", i
+1, nrexcl
);
295 for (j
= 0; j
< nrexcl
; j
++)
297 fprintf(out
, "%5d ", excluded
[j
]);
303 } /* punch_QMMM_excl */
306 /* end of QMMM subroutines */
308 /* QMMM core routines */
310 t_QMrec
*mk_QMrec(void)
317 t_MMrec
*mk_MMrec(void)
324 static void init_QMrec(int grpnr
, t_QMrec
*qm
, int nr
, int *atomarray
,
325 gmx_mtop_t
*mtop
, t_inputrec
*ir
)
327 /* fills the t_QMrec struct of QM group grpnr
330 gmx_mtop_atomlookup_t alook
;
336 snew(qm
->indexQM
, nr
);
337 snew(qm
->shiftQM
, nr
); /* the shifts */
338 for (i
= 0; i
< nr
; i
++)
340 qm
->indexQM
[i
] = atomarray
[i
];
343 alook
= gmx_mtop_atomlookup_init(mtop
);
345 snew(qm
->atomicnumberQM
, nr
);
346 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
348 gmx_mtop_atomnr_to_atom(alook
, qm
->indexQM
[i
], &atom
);
349 qm
->nelectrons
+= mtop
->atomtypes
.atomnumber
[atom
->type
];
350 qm
->atomicnumberQM
[i
] = mtop
->atomtypes
.atomnumber
[atom
->type
];
353 gmx_mtop_atomlookup_destroy(alook
);
355 qm
->QMcharge
= ir
->opts
.QMcharge
[grpnr
];
356 qm
->multiplicity
= ir
->opts
.QMmult
[grpnr
];
357 qm
->nelectrons
-= ir
->opts
.QMcharge
[grpnr
];
359 qm
->QMmethod
= ir
->opts
.QMmethod
[grpnr
];
360 qm
->QMbasis
= ir
->opts
.QMbasis
[grpnr
];
361 /* trajectory surface hopping setup (Gaussian only) */
362 qm
->bSH
= ir
->opts
.bSH
[grpnr
];
363 qm
->CASorbitals
= ir
->opts
.CASorbitals
[grpnr
];
364 qm
->CASelectrons
= ir
->opts
.CASelectrons
[grpnr
];
365 qm
->SAsteps
= ir
->opts
.SAsteps
[grpnr
];
366 qm
->SAon
= ir
->opts
.SAon
[grpnr
];
367 qm
->SAoff
= ir
->opts
.SAoff
[grpnr
];
368 /* hack to prevent gaussian from reinitializing all the time */
369 qm
->nQMcpus
= 0; /* number of CPU's to be used by g01, is set
370 * upon initializing gaussian
373 /* print the current layer to allow users to check their input */
374 fprintf(stderr
, "Layer %d\nnr of QM atoms %d\n", grpnr
, nr
);
375 fprintf(stderr
, "QMlevel: %s/%s\n\n",
376 eQMmethod_names
[qm
->QMmethod
], eQMbasis_names
[qm
->QMbasis
]);
379 snew(qm
->frontatoms
, nr
);
380 /* Lennard-Jones coefficients */
383 /* do we optimize the QM separately using the algorithms of the QM program??
385 qm
->bTS
= ir
->opts
.bTS
[grpnr
];
386 qm
->bOPT
= ir
->opts
.bOPT
[grpnr
];
390 t_QMrec
*copy_QMrec(t_QMrec
*qm
)
392 /* copies the contents of qm into a new t_QMrec struct */
399 qmcopy
->nrQMatoms
= qm
->nrQMatoms
;
400 snew(qmcopy
->xQM
, qmcopy
->nrQMatoms
);
401 snew(qmcopy
->indexQM
, qmcopy
->nrQMatoms
);
402 snew(qmcopy
->atomicnumberQM
, qm
->nrQMatoms
);
403 snew(qmcopy
->shiftQM
, qmcopy
->nrQMatoms
); /* the shifts */
404 for (i
= 0; i
< qmcopy
->nrQMatoms
; i
++)
406 qmcopy
->shiftQM
[i
] = qm
->shiftQM
[i
];
407 qmcopy
->indexQM
[i
] = qm
->indexQM
[i
];
408 qmcopy
->atomicnumberQM
[i
] = qm
->atomicnumberQM
[i
];
410 qmcopy
->nelectrons
= qm
->nelectrons
;
411 qmcopy
->multiplicity
= qm
->multiplicity
;
412 qmcopy
->QMcharge
= qm
->QMcharge
;
413 qmcopy
->nelectrons
= qm
->nelectrons
;
414 qmcopy
->QMmethod
= qm
->QMmethod
;
415 qmcopy
->QMbasis
= qm
->QMbasis
;
416 /* trajectory surface hopping setup (Gaussian only) */
417 qmcopy
->bSH
= qm
->bSH
;
418 qmcopy
->CASorbitals
= qm
->CASorbitals
;
419 qmcopy
->CASelectrons
= qm
->CASelectrons
;
420 qmcopy
->SAsteps
= qm
->SAsteps
;
421 qmcopy
->SAon
= qm
->SAon
;
422 qmcopy
->SAoff
= qm
->SAoff
;
423 qmcopy
->bOPT
= qm
->bOPT
;
425 /* Gaussian init. variables */
426 qmcopy
->nQMcpus
= qm
->nQMcpus
;
427 for (i
= 0; i
< DIM
; i
++)
429 qmcopy
->SHbasis
[i
] = qm
->SHbasis
[i
];
431 qmcopy
->QMmem
= qm
->QMmem
;
432 qmcopy
->accuracy
= qm
->accuracy
;
433 qmcopy
->cpmcscf
= qm
->cpmcscf
;
434 qmcopy
->SAstep
= qm
->SAstep
;
435 snew(qmcopy
->frontatoms
, qm
->nrQMatoms
);
436 snew(qmcopy
->c12
, qmcopy
->nrQMatoms
);
437 snew(qmcopy
->c6
, qmcopy
->nrQMatoms
);
438 if (qmcopy
->bTS
|| qmcopy
->bOPT
)
440 for (i
= 1; i
< qmcopy
->nrQMatoms
; i
++)
442 qmcopy
->frontatoms
[i
] = qm
->frontatoms
[i
];
443 qmcopy
->c12
[i
] = qm
->c12
[i
];
444 qmcopy
->c6
[i
] = qm
->c6
[i
];
452 t_QMMMrec
*mk_QMMMrec(void)
463 void init_QMMMrec(t_commrec
*cr
,
468 /* we put the atomsnumbers of atoms that belong to the QMMM group in
469 * an array that will be copied later to QMMMrec->indexQM[..]. Also
470 * it will be used to create an QMMMrec->bQMMM index array that
471 * simply contains true/false for QM and MM (the other) atoms.
474 gmx_groups_t
*groups
;
475 atom_id
*qm_arr
= NULL
, vsite
, ai
, aj
;
476 int qm_max
= 0, qm_nr
= 0, i
, j
, jmax
, k
, l
, nrvsite2
= 0;
481 gmx_mtop_atomloop_all_t aloop
;
483 gmx_mtop_ilistloop_all_t iloop
;
486 gmx_mtop_atomlookup_t alook
;
488 c6au
= (HARTREE2KJ
*AVOGADRO
*std::pow(BOHR2NM
, 6));
489 c12au
= (HARTREE2KJ
*AVOGADRO
*std::pow(BOHR2NM
, 12));
490 /* issue a fatal if the user wants to run with more than one node */
493 gmx_fatal(FARGS
, "QM/MM does not work in parallel, use a single rank instead\n");
496 /* Make a local copy of the QMMMrec */
499 /* bQMMM[..] is an array containing TRUE/FALSE for atoms that are
500 * QM/not QM. We first set all elemenst at false. Afterwards we use
501 * the qm_arr (=MMrec->indexQM) to changes the elements
502 * corresponding to the QM atoms at TRUE. */
504 qr
->QMMMscheme
= ir
->QMMMscheme
;
506 /* we take the possibility into account that a user has
507 * defined more than one QM group:
509 /* an ugly work-around in case there is only one group In this case
510 * the whole system is treated as QM. Otherwise the second group is
511 * always the rest of the total system and is treated as MM.
514 /* small problem if there is only QM.... so no MM */
516 jmax
= ir
->opts
.ngQM
;
518 if (qr
->QMMMscheme
== eQMMMschemeoniom
)
520 qr
->nrQMlayers
= jmax
;
527 groups
= &mtop
->groups
;
529 /* there are jmax groups of QM atoms. In case of multiple QM groups
530 * I assume that the users wants to do ONIOM. However, maybe it
531 * should also be possible to define more than one QM subsystem with
532 * independent neighbourlists. I have to think about
536 for (j
= 0; j
< jmax
; j
++)
539 aloop
= gmx_mtop_atomloop_all_init(mtop
);
540 while (gmx_mtop_atomloop_all_next(aloop
, &i
, &atom
))
545 srenew(qm_arr
, qm_max
);
547 if (ggrpnr(groups
, egcQMMM
, i
) == j
)
553 if (qr
->QMMMscheme
== eQMMMschemeoniom
)
555 /* add the atoms to the bQMMM array
558 /* I assume that users specify the QM groups from small to
559 * big(ger) in the mdp file
561 qr
->qm
[j
] = mk_QMrec();
562 /* we need to throw out link atoms that in the previous layer
563 * existed to separate this QMlayer from the previous
564 * QMlayer. We use the iatoms array in the idef for that
565 * purpose. If all atoms defining the current Link Atom (Dummy2)
566 * are part of the current QM layer it needs to be removed from
569 iloop
= gmx_mtop_ilistloop_all_init(mtop
);
570 while (gmx_mtop_ilistloop_all_next(iloop
, &ilist_mol
, &a_offset
))
572 nrvsite2
= ilist_mol
[F_VSITE2
].nr
;
573 iatoms
= ilist_mol
[F_VSITE2
].iatoms
;
575 for (k
= 0; k
< nrvsite2
; k
+= 4)
577 vsite
= a_offset
+ iatoms
[k
+1]; /* the vsite */
578 ai
= a_offset
+ iatoms
[k
+2]; /* constructing atom */
579 aj
= a_offset
+ iatoms
[k
+3]; /* constructing atom */
580 if (ggrpnr(groups
, egcQMMM
, vsite
) == ggrpnr(groups
, egcQMMM
, ai
)
582 ggrpnr(groups
, egcQMMM
, vsite
) == ggrpnr(groups
, egcQMMM
, aj
))
584 /* this dummy link atom needs to be removed from the qm_arr
585 * before making the QMrec of this layer!
587 for (i
= 0; i
< qm_nr
; i
++)
589 if (qm_arr
[i
] == vsite
)
591 /* drop the element */
592 for (l
= i
; l
< qm_nr
; l
++)
594 qm_arr
[l
] = qm_arr
[l
+1];
603 /* store QM atoms in this layer in the QMrec and initialise layer
605 init_QMrec(j
, qr
->qm
[j
], qm_nr
, qm_arr
, mtop
, ir
);
607 /* we now store the LJ C6 and C12 parameters in QM rec in case
608 * we need to do an optimization
610 if (qr
->qm
[j
]->bOPT
|| qr
->qm
[j
]->bTS
)
612 for (i
= 0; i
< qm_nr
; i
++)
614 /* nbfp now includes the 6.0/12.0 derivative prefactors */
615 qr
->qm
[j
]->c6
[i
] = C6(fr
->nbfp
, mtop
->ffparams
.atnr
, atom
->type
, atom
->type
)/c6au
/6.0;
616 qr
->qm
[j
]->c12
[i
] = C12(fr
->nbfp
, mtop
->ffparams
.atnr
, atom
->type
, atom
->type
)/c12au
/12.0;
619 /* now we check for frontier QM atoms. These occur in pairs that
620 * construct the vsite
622 iloop
= gmx_mtop_ilistloop_all_init(mtop
);
623 while (gmx_mtop_ilistloop_all_next(iloop
, &ilist_mol
, &a_offset
))
625 nrvsite2
= ilist_mol
[F_VSITE2
].nr
;
626 iatoms
= ilist_mol
[F_VSITE2
].iatoms
;
628 for (k
= 0; k
< nrvsite2
; k
+= 4)
630 vsite
= a_offset
+ iatoms
[k
+1]; /* the vsite */
631 ai
= a_offset
+ iatoms
[k
+2]; /* constructing atom */
632 aj
= a_offset
+ iatoms
[k
+3]; /* constructing atom */
633 if (ggrpnr(groups
, egcQMMM
, ai
) < (groups
->grps
[egcQMMM
].nr
-1) &&
634 (ggrpnr(groups
, egcQMMM
, aj
) >= (groups
->grps
[egcQMMM
].nr
-1)))
636 /* mark ai as frontier atom */
637 for (i
= 0; i
< qm_nr
; i
++)
639 if ( (qm_arr
[i
] == ai
) || (qm_arr
[i
] == vsite
) )
641 qr
->qm
[j
]->frontatoms
[i
] = TRUE
;
645 else if (ggrpnr(groups
, egcQMMM
, aj
) < (groups
->grps
[egcQMMM
].nr
-1) &&
646 (ggrpnr(groups
, egcQMMM
, ai
) >= (groups
->grps
[egcQMMM
].nr
-1)))
648 /* mark aj as frontier atom */
649 for (i
= 0; i
< qm_nr
; i
++)
651 if ( (qm_arr
[i
] == aj
) || (qm_arr
[i
] == vsite
))
653 qr
->qm
[j
]->frontatoms
[i
] = TRUE
;
661 if (qr
->QMMMscheme
!= eQMMMschemeoniom
)
664 /* standard QMMM, all layers are merged together so there is one QM
665 * subsystem and one MM subsystem.
666 * Also we set the charges to zero in the md->charge arrays to prevent
667 * the innerloops from doubly counting the electostatic QM MM interaction
670 alook
= gmx_mtop_atomlookup_init(mtop
);
672 for (k
= 0; k
< qm_nr
; k
++)
674 gmx_mtop_atomnr_to_atom(alook
, qm_arr
[k
], &atom
);
678 qr
->qm
[0] = mk_QMrec();
679 /* store QM atoms in the QMrec and initialise
681 init_QMrec(0, qr
->qm
[0], qm_nr
, qm_arr
, mtop
, ir
);
682 if (qr
->qm
[0]->bOPT
|| qr
->qm
[0]->bTS
)
684 for (i
= 0; i
< qm_nr
; i
++)
686 gmx_mtop_atomnr_to_atom(alook
, qm_arr
[i
], &atom
);
687 /* nbfp now includes the 6.0/12.0 derivative prefactors */
688 qr
->qm
[0]->c6
[i
] = C6(fr
->nbfp
, mtop
->ffparams
.atnr
, atom
->type
, atom
->type
)/c6au
/6.0;
689 qr
->qm
[0]->c12
[i
] = C12(fr
->nbfp
, mtop
->ffparams
.atnr
, atom
->type
, atom
->type
)/c12au
/12.0;
693 /* find frontier atoms and mark them true in the frontieratoms array.
695 for (i
= 0; i
< qm_nr
; i
++)
697 gmx_mtop_atomnr_to_ilist(alook
, qm_arr
[i
], &ilist_mol
, &a_offset
);
698 nrvsite2
= ilist_mol
[F_VSITE2
].nr
;
699 iatoms
= ilist_mol
[F_VSITE2
].iatoms
;
701 for (k
= 0; k
< nrvsite2
; k
+= 4)
703 vsite
= a_offset
+ iatoms
[k
+1]; /* the vsite */
704 ai
= a_offset
+ iatoms
[k
+2]; /* constructing atom */
705 aj
= a_offset
+ iatoms
[k
+3]; /* constructing atom */
706 if (ggrpnr(groups
, egcQMMM
, ai
) < (groups
->grps
[egcQMMM
].nr
-1) &&
707 (ggrpnr(groups
, egcQMMM
, aj
) >= (groups
->grps
[egcQMMM
].nr
-1)))
709 /* mark ai as frontier atom */
710 if ( (qm_arr
[i
] == ai
) || (qm_arr
[i
] == vsite
) )
712 qr
->qm
[0]->frontatoms
[i
] = TRUE
;
715 else if (ggrpnr(groups
, egcQMMM
, aj
) < (groups
->grps
[egcQMMM
].nr
-1) &&
716 (ggrpnr(groups
, egcQMMM
, ai
) >= (groups
->grps
[egcQMMM
].nr
-1)))
718 /* mark aj as frontier atom */
719 if ( (qm_arr
[i
] == aj
) || (qm_arr
[i
] == vsite
) )
721 qr
->qm
[0]->frontatoms
[i
] = TRUE
;
727 gmx_mtop_atomlookup_destroy(alook
);
729 /* MM rec creation */
731 mm
->scalefactor
= ir
->scalefactor
;
732 mm
->nrMMatoms
= (mtop
->natoms
)-(qr
->qm
[0]->nrQMatoms
); /* rest of the atoms */
736 { /* MM rec creation */
738 mm
->scalefactor
= ir
->scalefactor
;
743 /* these variables get updated in the update QMMMrec */
745 if (qr
->nrQMlayers
== 1)
747 /* with only one layer there is only one initialisation
748 * needed. Multilayer is a bit more complicated as it requires
749 * re-initialisation at every step of the simulation. This is due
750 * to the use of COMMON blocks in the fortran QM subroutines.
752 if (qr
->qm
[0]->QMmethod
< eQMmethodRHF
)
755 /* semi-empiprical 1-layer ONIOM calculation requested (mopac93) */
756 init_mopac(qr
->qm
[0]);
758 gmx_fatal(FARGS
, "Semi-empirical QM only supported with Mopac.");
763 /* ab initio calculation requested (gamess/gaussian/ORCA) */
765 init_gamess(cr
, qr
->qm
[0], qr
->mm
);
766 #elif GMX_QMMM_GAUSSIAN
767 init_gaussian(qr
->qm
[0]);
769 init_orca(qr
->qm
[0]);
771 gmx_fatal(FARGS
, "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
777 void update_QMMMrec(t_commrec
*cr
,
784 /* updates the coordinates of both QM atoms and MM atoms and stores
785 * them in the QMMMrec.
787 * NOTE: is NOT yet working if there are no PBC. Also in ns.c, simple
788 * ns needs to be fixed!
791 mm_max
= 0, mm_nr
= 0, mm_nr_new
, i
, j
, is
, k
, shift
;
793 *mm_j_particles
= NULL
, *qm_i_particles
= NULL
;
807 *parallelMMarray
= NULL
;
811 c6au
= (HARTREE2KJ
*AVOGADRO
*std::pow(BOHR2NM
, 6));
812 c12au
= (HARTREE2KJ
*AVOGADRO
*std::pow(BOHR2NM
, 12));
814 /* every cpu has this array. On every processor we fill this array
815 * with 1's and 0's. 1's indicate the atoms is a QM atom on the
816 * current cpu in a later stage these arrays are all summed. indexes
817 * > 0 indicate the atom is a QM atom. Every node therefore knows
818 * whcih atoms are part of the QM subsystem.
820 /* copy some pointers */
823 QMMMlist
= fr
->QMMMlist
;
827 /* init_pbc(box); needs to be called first, see pbc.h */
828 set_pbc_dd(&pbc
, fr
->ePBC
, DOMAINDECOMP(cr
) ? cr
->dd
: NULL
, FALSE
, box
);
829 /* only in standard (normal) QMMM we need the neighbouring MM
830 * particles to provide a electric field of point charges for the QM
833 if (qr
->QMMMscheme
== eQMMMschemenormal
) /* also implies 1 QM-layer */
835 /* we NOW create/update a number of QMMMrec entries:
837 * 1) the shiftQM, containing the shifts of the QM atoms
839 * 2) the indexMM array, containing the index of the MM atoms
841 * 3) the shiftMM, containing the shifts of the MM atoms
843 * 4) the shifted coordinates of the MM atoms
845 * the shifts are used for computing virial of the QM/MM particles.
847 qm
= qr
->qm
[0]; /* in case of normal QMMM, there is only one group */
848 snew(qm_i_particles
, QMMMlist
.nri
);
851 qm_i_particles
[0].shift
= XYZ2IS(0, 0, 0);
852 for (i
= 0; i
< QMMMlist
.nri
; i
++)
854 qm_i_particles
[i
].j
= QMMMlist
.iinr
[i
];
858 qm_i_particles
[i
].shift
= pbc_dx_aiuc(&pbc
, x
[QMMMlist
.iinr
[0]],
859 x
[QMMMlist
.iinr
[i
]], dx
);
862 /* However, since nri >= nrQMatoms, we do a quicksort, and throw
863 * out double, triple, etc. entries later, as we do for the MM
867 /* compute the shift for the MM j-particles with respect to
868 * the QM i-particle and store them.
871 crd
[0] = IS2X(QMMMlist
.shift
[i
]) + IS2X(qm_i_particles
[i
].shift
);
872 crd
[1] = IS2Y(QMMMlist
.shift
[i
]) + IS2Y(qm_i_particles
[i
].shift
);
873 crd
[2] = IS2Z(QMMMlist
.shift
[i
]) + IS2Z(qm_i_particles
[i
].shift
);
874 is
= static_cast<int>(XYZ2IS(crd
[0], crd
[1], crd
[2]));
875 for (j
= QMMMlist
.jindex
[i
];
876 j
< QMMMlist
.jindex
[i
+1];
882 srenew(mm_j_particles
, mm_max
);
885 mm_j_particles
[mm_nr
].j
= QMMMlist
.jjnr
[j
];
886 mm_j_particles
[mm_nr
].shift
= is
;
891 /* quicksort QM and MM shift arrays and throw away multiple entries */
895 qsort(qm_i_particles
, QMMMlist
.nri
,
896 (size_t)sizeof(qm_i_particles
[0]),
898 /* The mm_j_particles argument to qsort is not allowed to be NULL */
901 qsort(mm_j_particles
, mm_nr
,
902 (size_t)sizeof(mm_j_particles
[0]),
905 /* remove multiples in the QM shift array, since in init_QMMM() we
906 * went through the atom numbers from 0 to md.nr, the order sorted
907 * here matches the one of QMindex already.
910 for (i
= 0; i
< QMMMlist
.nri
; i
++)
912 if (i
== 0 || qm_i_particles
[i
].j
!= qm_i_particles
[i
-1].j
)
914 qm_i_particles
[j
++] = qm_i_particles
[i
];
918 if (qm
->bTS
|| qm
->bOPT
)
920 /* only remove double entries for the MM array */
921 for (i
= 0; i
< mm_nr
; i
++)
923 if ((i
== 0 || mm_j_particles
[i
].j
!= mm_j_particles
[i
-1].j
)
924 && !md
->bQM
[mm_j_particles
[i
].j
])
926 mm_j_particles
[mm_nr_new
++] = mm_j_particles
[i
];
930 /* we also remove mm atoms that have no charges!
931 * actually this is already done in the ns.c
935 for (i
= 0; i
< mm_nr
; i
++)
937 if ((i
== 0 || mm_j_particles
[i
].j
!= mm_j_particles
[i
-1].j
)
938 && !md
->bQM
[mm_j_particles
[i
].j
]
939 && (md
->chargeA
[mm_j_particles
[i
].j
]
940 || (md
->chargeB
&& md
->chargeB
[mm_j_particles
[i
].j
])))
942 mm_j_particles
[mm_nr_new
++] = mm_j_particles
[i
];
947 /* store the data retrieved above into the QMMMrec
950 /* Keep the compiler happy,
951 * shift will always be set in the loop for i=0
954 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
956 /* not all qm particles might have appeared as i
957 * particles. They might have been part of the same charge
958 * group for instance.
960 if (qm
->indexQM
[i
] == qm_i_particles
[k
].j
)
962 shift
= qm_i_particles
[k
++].shift
;
964 /* use previous shift, assuming they belong the same charge
968 qm
->shiftQM
[i
] = shift
;
971 /* parallel excecution */
974 snew(parallelMMarray
, 2*(md
->nr
));
975 /* only MM particles have a 1 at their atomnumber. The second part
976 * of the array contains the shifts. Thus:
977 * p[i]=1/0 depending on wether atomnumber i is a MM particle in the QM
978 * step or not. p[i+md->nr] is the shift of atomnumber i.
980 for (i
= 0; i
< 2*(md
->nr
); i
++)
982 parallelMMarray
[i
] = 0;
985 for (i
= 0; i
< mm_nr
; i
++)
987 parallelMMarray
[mm_j_particles
[i
].j
] = 1;
988 parallelMMarray
[mm_j_particles
[i
].j
+(md
->nr
)] = mm_j_particles
[i
].shift
;
990 gmx_sumi(md
->nr
, parallelMMarray
, cr
);
994 for (i
= 0; i
< md
->nr
; i
++)
996 if (parallelMMarray
[i
])
1001 srenew(mm
->indexMM
, mm_max
);
1002 srenew(mm
->shiftMM
, mm_max
);
1004 mm
->indexMM
[mm_nr
] = i
;
1005 mm
->shiftMM
[mm_nr
++] = parallelMMarray
[i
+md
->nr
]/parallelMMarray
[i
];
1008 mm
->nrMMatoms
= mm_nr
;
1009 free(parallelMMarray
);
1011 /* serial execution */
1014 mm
->nrMMatoms
= mm_nr
;
1015 srenew(mm
->shiftMM
, mm_nr
);
1016 srenew(mm
->indexMM
, mm_nr
);
1017 for (i
= 0; i
< mm_nr
; i
++)
1019 mm
->indexMM
[i
] = mm_j_particles
[i
].j
;
1020 mm
->shiftMM
[i
] = mm_j_particles
[i
].shift
;
1024 /* (re) allocate memory for the MM coordiate array. The QM
1025 * coordinate array was already allocated in init_QMMM, and is
1026 * only (re)filled in the update_QMMM_coordinates routine
1028 srenew(mm
->xMM
, mm
->nrMMatoms
);
1029 /* now we (re) fill the array that contains the MM charges with
1030 * the forcefield charges. If requested, these charges will be
1031 * scaled by a factor
1033 srenew(mm
->MMcharges
, mm
->nrMMatoms
);
1034 for (i
= 0; i
< mm
->nrMMatoms
; i
++) /* no free energy yet */
1036 mm
->MMcharges
[i
] = md
->chargeA
[mm
->indexMM
[i
]]*mm
->scalefactor
;
1038 if (qm
->bTS
|| qm
->bOPT
)
1040 /* store (copy) the c6 and c12 parameters into the MMrec struct
1042 srenew(mm
->c6
, mm
->nrMMatoms
);
1043 srenew(mm
->c12
, mm
->nrMMatoms
);
1044 for (i
= 0; i
< mm
->nrMMatoms
; i
++)
1046 /* nbfp now includes the 6.0/12.0 derivative prefactors */
1047 mm
->c6
[i
] = C6(fr
->nbfp
, top
->idef
.atnr
, md
->typeA
[mm
->indexMM
[i
]], md
->typeA
[mm
->indexMM
[i
]])/c6au
/6.0;
1048 mm
->c12
[i
] = C12(fr
->nbfp
, top
->idef
.atnr
, md
->typeA
[mm
->indexMM
[i
]], md
->typeA
[mm
->indexMM
[i
]])/c12au
/12.0;
1050 punch_QMMM_excl(qr
->qm
[0], mm
, &(top
->excls
));
1052 /* the next routine fills the coordinate fields in the QMMM rec of
1053 * both the qunatum atoms and the MM atoms, using the shifts
1057 update_QMMM_coord(x
, fr
, qr
->qm
[0], qr
->mm
);
1058 free(qm_i_particles
);
1059 free(mm_j_particles
);
1061 else /* ONIOM */ /* ????? */
1064 /* do for each layer */
1065 for (j
= 0; j
< qr
->nrQMlayers
; j
++)
1068 qm
->shiftQM
[0] = XYZ2IS(0, 0, 0);
1069 for (i
= 1; i
< qm
->nrQMatoms
; i
++)
1071 qm
->shiftQM
[i
] = pbc_dx_aiuc(&pbc
, x
[qm
->indexQM
[0]], x
[qm
->indexQM
[i
]],
1074 update_QMMM_coord(x
, fr
, qm
, mm
);
1077 } /* update_QMMM_rec */
1080 real
calculate_QMMM(t_commrec
*cr
,
1086 /* a selection for the QM package depending on which is requested
1087 * (Gaussian, GAMESS-UK, MOPAC or ORCA) needs to be implemented here. Now
1088 * it works through defines.... Not so nice yet
1097 *forces
= NULL
, *fshift
= NULL
,
1098 *forces2
= NULL
, *fshift2
= NULL
; /* needed for multilayer ONIOM */
1101 /* make a local copy the QMMMrec pointer
1106 /* now different procedures are carried out for one layer ONION and
1107 * normal QMMM on one hand and multilayer oniom on the other
1109 if (qr
->QMMMscheme
== eQMMMschemenormal
|| qr
->nrQMlayers
== 1)
1112 snew(forces
, (qm
->nrQMatoms
+mm
->nrMMatoms
));
1113 snew(fshift
, (qm
->nrQMatoms
+mm
->nrMMatoms
));
1114 QMener
= call_QMroutine(cr
, fr
, qm
, mm
, forces
, fshift
);
1115 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
1117 for (j
= 0; j
< DIM
; j
++)
1119 f
[qm
->indexQM
[i
]][j
] -= forces
[i
][j
];
1120 fr
->fshift
[qm
->shiftQM
[i
]][j
] += fshift
[i
][j
];
1123 for (i
= 0; i
< mm
->nrMMatoms
; i
++)
1125 for (j
= 0; j
< DIM
; j
++)
1127 f
[mm
->indexMM
[i
]][j
] -= forces
[qm
->nrQMatoms
+i
][j
];
1128 fr
->fshift
[mm
->shiftMM
[i
]][j
] += fshift
[qm
->nrQMatoms
+i
][j
];
1135 else /* Multi-layer ONIOM */
1137 for (i
= 0; i
< qr
->nrQMlayers
-1; i
++) /* last layer is special */
1140 qm2
= copy_QMrec(qr
->qm
[i
+1]);
1142 qm2
->nrQMatoms
= qm
->nrQMatoms
;
1144 for (j
= 0; j
< qm2
->nrQMatoms
; j
++)
1146 for (k
= 0; k
< DIM
; k
++)
1148 qm2
->xQM
[j
][k
] = qm
->xQM
[j
][k
];
1150 qm2
->indexQM
[j
] = qm
->indexQM
[j
];
1151 qm2
->atomicnumberQM
[j
] = qm
->atomicnumberQM
[j
];
1152 qm2
->shiftQM
[j
] = qm
->shiftQM
[j
];
1155 qm2
->QMcharge
= qm
->QMcharge
;
1156 /* this layer at the higher level of theory */
1157 srenew(forces
, qm
->nrQMatoms
);
1158 srenew(fshift
, qm
->nrQMatoms
);
1159 /* we need to re-initialize the QMroutine every step... */
1160 init_QMroutine(cr
, qm
, mm
);
1161 QMener
+= call_QMroutine(cr
, fr
, qm
, mm
, forces
, fshift
);
1163 /* this layer at the lower level of theory */
1164 srenew(forces2
, qm
->nrQMatoms
);
1165 srenew(fshift2
, qm
->nrQMatoms
);
1166 init_QMroutine(cr
, qm2
, mm
);
1167 QMener
-= call_QMroutine(cr
, fr
, qm2
, mm
, forces2
, fshift2
);
1168 /* E = E1high-E1low The next layer includes the current layer at
1169 * the lower level of theory, which provides + E2low
1170 * this is similar for gradients
1172 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
1174 for (j
= 0; j
< DIM
; j
++)
1176 f
[qm
->indexQM
[i
]][j
] -= (forces
[i
][j
]-forces2
[i
][j
]);
1177 fr
->fshift
[qm
->shiftQM
[i
]][j
] += (fshift
[i
][j
]-fshift2
[i
][j
]);
1182 /* now the last layer still needs to be done: */
1183 qm
= qr
->qm
[qr
->nrQMlayers
-1]; /* C counts from 0 */
1184 init_QMroutine(cr
, qm
, mm
);
1185 srenew(forces
, qm
->nrQMatoms
);
1186 srenew(fshift
, qm
->nrQMatoms
);
1187 QMener
+= call_QMroutine(cr
, fr
, qm
, mm
, forces
, fshift
);
1188 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
1190 for (j
= 0; j
< DIM
; j
++)
1192 f
[qm
->indexQM
[i
]][j
] -= forces
[i
][j
];
1193 fr
->fshift
[qm
->shiftQM
[i
]][j
] += fshift
[i
][j
];
1201 if (qm
->bTS
|| qm
->bOPT
)
1203 /* qm[0] still contains the largest ONIOM QM subsystem
1204 * we take the optimized coordiates and put the in x[]
1206 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
1208 for (j
= 0; j
< DIM
; j
++)
1210 x
[qm
->indexQM
[i
]][j
] = qm
->xQM
[i
][j
];
1215 } /* calculate_QMMM */
1217 /* end of QMMM core routines */