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, 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.
43 #include "types/commrec.h"
45 #include "gromacs/utility/smalloc.h"
46 #include "gromacs/math/units.h"
48 #include "gromacs/math/vec.h"
51 #include "gromacs/fileio/confio.h"
63 #include "gromacs/utility/fatalerror.h"
66 #include "mtop_util.h"
69 /* declarations of the interfaces to the QM packages. The _SH indicate
70 * the QM interfaces can be used for Surface Hopping simulations
72 #ifdef GMX_QMMM_GAMESS
73 /* GAMESS interface */
76 init_gamess(t_commrec
*cr
, t_QMrec
*qm
, t_MMrec
*mm
);
79 call_gamess(t_commrec
*cr
, t_forcerec
*fr
,
80 t_QMrec
*qm
, t_MMrec
*mm
, rvec f
[], rvec fshift
[]);
82 #elif defined GMX_QMMM_MOPAC
86 init_mopac(t_commrec
*cr
, t_QMrec
*qm
, t_MMrec
*mm
);
89 call_mopac(t_commrec
*cr
, t_forcerec
*fr
, t_QMrec
*qm
,
90 t_MMrec
*mm
, rvec f
[], rvec fshift
[]);
93 call_mopac_SH(t_commrec
*cr
, t_forcerec
*fr
, t_QMrec
*qm
,
94 t_MMrec
*mm
, rvec f
[], rvec fshift
[]);
96 #elif defined GMX_QMMM_GAUSSIAN
97 /* GAUSSIAN interface */
100 init_gaussian(t_commrec
*cr
, t_QMrec
*qm
, t_MMrec
*mm
);
103 call_gaussian_SH(t_commrec
*cr
, t_forcerec
*fr
, t_QMrec
*qm
,
104 t_MMrec
*mm
, rvec f
[], rvec fshift
[]);
107 call_gaussian(t_commrec
*cr
, t_forcerec
*fr
, t_QMrec
*qm
,
108 t_MMrec
*mm
, rvec f
[], rvec fshift
[]);
110 #elif defined GMX_QMMM_ORCA
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 static int int_comp(const void *a
, const void *b
)
144 return (*(int *)a
) - (*(int *)b
);
148 static int QMlayer_comp(const void *a
, const void *b
)
151 return (int)(((t_QMrec
*)a
)->nrQMatoms
)-(int)(((t_QMrec
*)b
)->nrQMatoms
);
155 real
call_QMroutine(t_commrec gmx_unused
*cr
, t_forcerec gmx_unused
*fr
, t_QMrec gmx_unused
*qm
,
156 t_MMrec gmx_unused
*mm
, rvec gmx_unused f
[], rvec gmx_unused fshift
[])
158 /* makes a call to the requested QM routine (qm->QMmethod)
159 * Note that f is actually the gradient, i.e. -f
164 /* do a semi-empiprical calculation */
166 if (qm
->QMmethod
< eQMmethodRHF
&& !(mm
->nrMMatoms
))
168 #ifdef GMX_QMMM_MOPAC
171 QMener
= call_mopac_SH(cr
, fr
, qm
, mm
, f
, fshift
);
175 QMener
= call_mopac(cr
, fr
, qm
, mm
, f
, fshift
);
178 gmx_fatal(FARGS
, "Semi-empirical QM only supported with Mopac.");
183 /* do an ab-initio calculation */
184 if (qm
->bSH
&& qm
->QMmethod
== eQMmethodCASSCF
)
186 #ifdef GMX_QMMM_GAUSSIAN
187 QMener
= call_gaussian_SH(cr
, fr
, qm
, mm
, f
, fshift
);
189 gmx_fatal(FARGS
, "Ab-initio Surface-hopping only supported with Gaussian.");
194 #ifdef GMX_QMMM_GAMESS
195 QMener
= call_gamess(cr
, fr
, qm
, mm
, f
, fshift
);
196 #elif defined GMX_QMMM_GAUSSIAN
197 QMener
= call_gaussian(cr
, fr
, qm
, mm
, f
, fshift
);
198 #elif defined GMX_QMMM_ORCA
199 QMener
= call_orca(fr
, qm
, mm
, f
, fshift
);
201 gmx_fatal(FARGS
, "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
208 void init_QMroutine(t_commrec gmx_unused
*cr
, t_QMrec gmx_unused
*qm
, t_MMrec gmx_unused
*mm
)
210 /* makes a call to the requested QM routine (qm->QMmethod)
212 if (qm
->QMmethod
< eQMmethodRHF
)
214 #ifdef GMX_QMMM_MOPAC
215 /* do a semi-empiprical calculation */
216 init_mopac(cr
, qm
, mm
);
218 gmx_fatal(FARGS
, "Semi-empirical QM only supported with Mopac.");
223 /* do an ab-initio calculation */
224 #ifdef GMX_QMMM_GAMESS
225 init_gamess(cr
, qm
, mm
);
226 #elif defined GMX_QMMM_GAUSSIAN
227 init_gaussian(cr
, qm
, mm
);
228 #elif defined GMX_QMMM_ORCA
231 gmx_fatal(FARGS
, "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
234 } /* init_QMroutine */
236 void update_QMMM_coord(rvec x
[], t_forcerec
*fr
, t_QMrec
*qm
, t_MMrec
*mm
)
238 /* shifts the QM and MM particles into the central box and stores
239 * these shifted coordinates in the coordinate arrays of the
240 * QMMMrec. These coordinates are passed on the QM subroutines.
245 /* shift the QM atoms into the central box
247 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
249 rvec_sub(x
[qm
->indexQM
[i
]], fr
->shift_vec
[qm
->shiftQM
[i
]], qm
->xQM
[i
]);
251 /* also shift the MM atoms into the central box, if any
253 for (i
= 0; i
< mm
->nrMMatoms
; i
++)
255 rvec_sub(x
[mm
->indexMM
[i
]], fr
->shift_vec
[mm
->shiftMM
[i
]], mm
->xMM
[i
]);
257 } /* update_QMMM_coord */
259 static void punch_QMMM_excl(t_QMrec
*qm
, t_MMrec
*mm
, t_blocka
*excls
)
261 /* punch a file containing the bonded interactions of each QM
262 * atom with MM atoms. These need to be excluded in the QM routines
263 * Only needed in case of QM/MM optimizations
268 i
, j
, k
, nrexcl
= 0, *excluded
= NULL
, max
= 0;
271 out
= fopen("QMMMexcl.dat", "w");
273 /* this can be done more efficiently I think
275 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
278 for (j
= excls
->index
[qm
->indexQM
[i
]];
279 j
< excls
->index
[qm
->indexQM
[i
]+1];
282 for (k
= 0; k
< mm
->nrMMatoms
; k
++)
284 if (mm
->indexMM
[k
] == excls
->a
[j
]) /* the excluded MM atom */
289 srenew(excluded
, max
);
291 excluded
[nrexcl
++] = k
;
297 fprintf(out
, "%5d %5d\n", i
+1, nrexcl
);
298 for (j
= 0; j
< nrexcl
; j
++)
300 fprintf(out
, "%5d ", excluded
[j
]);
306 } /* punch_QMMM_excl */
309 /* end of QMMM subroutines */
311 /* QMMM core routines */
313 t_QMrec
*mk_QMrec(void)
320 t_MMrec
*mk_MMrec(void)
327 static void init_QMrec(int grpnr
, t_QMrec
*qm
, int nr
, int *atomarray
,
328 gmx_mtop_t
*mtop
, t_inputrec
*ir
)
330 /* fills the t_QMrec struct of QM group grpnr
333 gmx_mtop_atomlookup_t alook
;
339 snew(qm
->indexQM
, nr
);
340 snew(qm
->shiftQM
, nr
); /* the shifts */
341 for (i
= 0; i
< nr
; i
++)
343 qm
->indexQM
[i
] = atomarray
[i
];
346 alook
= gmx_mtop_atomlookup_init(mtop
);
348 snew(qm
->atomicnumberQM
, nr
);
349 for (i
= 0; i
< qm
->nrQMatoms
; i
++)
351 gmx_mtop_atomnr_to_atom(alook
, qm
->indexQM
[i
], &atom
);
352 qm
->nelectrons
+= mtop
->atomtypes
.atomnumber
[atom
->type
];
353 qm
->atomicnumberQM
[i
] = mtop
->atomtypes
.atomnumber
[atom
->type
];
356 gmx_mtop_atomlookup_destroy(alook
);
358 qm
->QMcharge
= ir
->opts
.QMcharge
[grpnr
];
359 qm
->multiplicity
= ir
->opts
.QMmult
[grpnr
];
360 qm
->nelectrons
-= ir
->opts
.QMcharge
[grpnr
];
362 qm
->QMmethod
= ir
->opts
.QMmethod
[grpnr
];
363 qm
->QMbasis
= ir
->opts
.QMbasis
[grpnr
];
364 /* trajectory surface hopping setup (Gaussian only) */
365 qm
->bSH
= ir
->opts
.bSH
[grpnr
];
366 qm
->CASorbitals
= ir
->opts
.CASorbitals
[grpnr
];
367 qm
->CASelectrons
= ir
->opts
.CASelectrons
[grpnr
];
368 qm
->SAsteps
= ir
->opts
.SAsteps
[grpnr
];
369 qm
->SAon
= ir
->opts
.SAon
[grpnr
];
370 qm
->SAoff
= ir
->opts
.SAoff
[grpnr
];
371 /* hack to prevent gaussian from reinitializing all the time */
372 qm
->nQMcpus
= 0; /* number of CPU's to be used by g01, is set
373 * upon initializing gaussian
376 /* print the current layer to allow users to check their input */
377 fprintf(stderr
, "Layer %d\nnr of QM atoms %d\n", grpnr
, nr
);
378 fprintf(stderr
, "QMlevel: %s/%s\n\n",
379 eQMmethod_names
[qm
->QMmethod
], eQMbasis_names
[qm
->QMbasis
]);
382 snew(qm
->frontatoms
, nr
);
383 /* Lennard-Jones coefficients */
386 /* do we optimize the QM separately using the algorithms of the QM program??
388 qm
->bTS
= ir
->opts
.bTS
[grpnr
];
389 qm
->bOPT
= ir
->opts
.bOPT
[grpnr
];
393 t_QMrec
*copy_QMrec(t_QMrec
*qm
)
395 /* copies the contents of qm into a new t_QMrec struct */
402 qmcopy
->nrQMatoms
= qm
->nrQMatoms
;
403 snew(qmcopy
->xQM
, qmcopy
->nrQMatoms
);
404 snew(qmcopy
->indexQM
, qmcopy
->nrQMatoms
);
405 snew(qmcopy
->atomicnumberQM
, qm
->nrQMatoms
);
406 snew(qmcopy
->shiftQM
, qmcopy
->nrQMatoms
); /* the shifts */
407 for (i
= 0; i
< qmcopy
->nrQMatoms
; i
++)
409 qmcopy
->shiftQM
[i
] = qm
->shiftQM
[i
];
410 qmcopy
->indexQM
[i
] = qm
->indexQM
[i
];
411 qmcopy
->atomicnumberQM
[i
] = qm
->atomicnumberQM
[i
];
413 qmcopy
->nelectrons
= qm
->nelectrons
;
414 qmcopy
->multiplicity
= qm
->multiplicity
;
415 qmcopy
->QMcharge
= qm
->QMcharge
;
416 qmcopy
->nelectrons
= qm
->nelectrons
;
417 qmcopy
->QMmethod
= qm
->QMmethod
;
418 qmcopy
->QMbasis
= qm
->QMbasis
;
419 /* trajectory surface hopping setup (Gaussian only) */
420 qmcopy
->bSH
= qm
->bSH
;
421 qmcopy
->CASorbitals
= qm
->CASorbitals
;
422 qmcopy
->CASelectrons
= qm
->CASelectrons
;
423 qmcopy
->SAsteps
= qm
->SAsteps
;
424 qmcopy
->SAon
= qm
->SAon
;
425 qmcopy
->SAoff
= qm
->SAoff
;
426 qmcopy
->bOPT
= qm
->bOPT
;
428 /* Gaussian init. variables */
429 qmcopy
->nQMcpus
= qm
->nQMcpus
;
430 for (i
= 0; i
< DIM
; i
++)
432 qmcopy
->SHbasis
[i
] = qm
->SHbasis
[i
];
434 qmcopy
->QMmem
= qm
->QMmem
;
435 qmcopy
->accuracy
= qm
->accuracy
;
436 qmcopy
->cpmcscf
= qm
->cpmcscf
;
437 qmcopy
->SAstep
= qm
->SAstep
;
438 snew(qmcopy
->frontatoms
, qm
->nrQMatoms
);
439 snew(qmcopy
->c12
, qmcopy
->nrQMatoms
);
440 snew(qmcopy
->c6
, qmcopy
->nrQMatoms
);
441 if (qmcopy
->bTS
|| qmcopy
->bOPT
)
443 for (i
= 1; i
< qmcopy
->nrQMatoms
; i
++)
445 qmcopy
->frontatoms
[i
] = qm
->frontatoms
[i
];
446 qmcopy
->c12
[i
] = qm
->c12
[i
];
447 qmcopy
->c6
[i
] = qm
->c6
[i
];
455 t_QMMMrec
*mk_QMMMrec(void)
466 void init_QMMMrec(t_commrec
*cr
,
471 /* we put the atomsnumbers of atoms that belong to the QMMM group in
472 * an array that will be copied later to QMMMrec->indexQM[..]. Also
473 * it will be used to create an QMMMrec->bQMMM index array that
474 * simply contains true/false for QM and MM (the other) atoms.
477 gmx_groups_t
*groups
;
478 atom_id
*qm_arr
= NULL
, vsite
, ai
, aj
;
479 int qm_max
= 0, qm_nr
= 0, i
, j
, jmax
, k
, l
, nrvsite2
= 0;
484 gmx_mtop_atomloop_all_t aloop
;
486 gmx_mtop_ilistloop_all_t iloop
;
489 gmx_mtop_atomlookup_t alook
;
491 c6au
= (HARTREE2KJ
*AVOGADRO
*pow(BOHR2NM
, 6));
492 c12au
= (HARTREE2KJ
*AVOGADRO
*pow(BOHR2NM
, 12));
493 /* issue a fatal if the user wants to run with more than one node */
496 gmx_fatal(FARGS
, "QM/MM does not work in parallel, use a single node instead\n");
499 /* Make a local copy of the QMMMrec */
502 /* bQMMM[..] is an array containing TRUE/FALSE for atoms that are
503 * QM/not QM. We first set all elemenst at false. Afterwards we use
504 * the qm_arr (=MMrec->indexQM) to changes the elements
505 * corresponding to the QM atoms at TRUE. */
507 qr
->QMMMscheme
= ir
->QMMMscheme
;
509 /* we take the possibility into account that a user has
510 * defined more than one QM group:
512 /* an ugly work-around in case there is only one group In this case
513 * the whole system is treated as QM. Otherwise the second group is
514 * always the rest of the total system and is treated as MM.
517 /* small problem if there is only QM.... so no MM */
519 jmax
= ir
->opts
.ngQM
;
521 if (qr
->QMMMscheme
== eQMMMschemeoniom
)
523 qr
->nrQMlayers
= jmax
;
530 groups
= &mtop
->groups
;
532 /* there are jmax groups of QM atoms. In case of multiple QM groups
533 * I assume that the users wants to do ONIOM. However, maybe it
534 * should also be possible to define more than one QM subsystem with
535 * independent neighbourlists. I have to think about
539 for (j
= 0; j
< jmax
; j
++)
542 aloop
= gmx_mtop_atomloop_all_init(mtop
);
543 while (gmx_mtop_atomloop_all_next(aloop
, &i
, &atom
))
548 srenew(qm_arr
, qm_max
);
550 if (ggrpnr(groups
, egcQMMM
, i
) == j
)
556 if (qr
->QMMMscheme
== eQMMMschemeoniom
)
558 /* add the atoms to the bQMMM array
561 /* I assume that users specify the QM groups from small to
562 * big(ger) in the mdp file
564 qr
->qm
[j
] = mk_QMrec();
565 /* we need to throw out link atoms that in the previous layer
566 * existed to separate this QMlayer from the previous
567 * QMlayer. We use the iatoms array in the idef for that
568 * purpose. If all atoms defining the current Link Atom (Dummy2)
569 * are part of the current QM layer it needs to be removed from
572 iloop
= gmx_mtop_ilistloop_all_init(mtop
);
573 while (gmx_mtop_ilistloop_all_next(iloop
, &ilist_mol
, &a_offset
))
575 nrvsite2
= ilist_mol
[F_VSITE2
].nr
;
576 iatoms
= ilist_mol
[F_VSITE2
].iatoms
;
578 for (k
= 0; k
< nrvsite2
; k
+= 4)
580 vsite
= a_offset
+ iatoms
[k
+1]; /* the vsite */
581 ai
= a_offset
+ iatoms
[k
+2]; /* constructing atom */
582 aj
= a_offset
+ iatoms
[k
+3]; /* constructing atom */
583 if (ggrpnr(groups
, egcQMMM
, vsite
) == ggrpnr(groups
, egcQMMM
, ai
)
585 ggrpnr(groups
, egcQMMM
, vsite
) == ggrpnr(groups
, egcQMMM
, aj
))
587 /* this dummy link atom needs to be removed from the qm_arr
588 * before making the QMrec of this layer!
590 for (i
= 0; i
< qm_nr
; i
++)
592 if (qm_arr
[i
] == vsite
)
594 /* drop the element */
595 for (l
= i
; l
< qm_nr
; l
++)
597 qm_arr
[l
] = qm_arr
[l
+1];
606 /* store QM atoms in this layer in the QMrec and initialise layer
608 init_QMrec(j
, qr
->qm
[j
], qm_nr
, qm_arr
, mtop
, ir
);
610 /* we now store the LJ C6 and C12 parameters in QM rec in case
611 * we need to do an optimization
613 if (qr
->qm
[j
]->bOPT
|| qr
->qm
[j
]->bTS
)
615 for (i
= 0; i
< qm_nr
; i
++)
617 /* nbfp now includes the 6.0/12.0 derivative prefactors */
618 qr
->qm
[j
]->c6
[i
] = C6(fr
->nbfp
, mtop
->ffparams
.atnr
, atom
->type
, atom
->type
)/c6au
/6.0;
619 qr
->qm
[j
]->c12
[i
] = C12(fr
->nbfp
, mtop
->ffparams
.atnr
, atom
->type
, atom
->type
)/c12au
/12.0;
622 /* now we check for frontier QM atoms. These occur in pairs that
623 * construct the vsite
625 iloop
= gmx_mtop_ilistloop_all_init(mtop
);
626 while (gmx_mtop_ilistloop_all_next(iloop
, &ilist_mol
, &a_offset
))
628 nrvsite2
= ilist_mol
[F_VSITE2
].nr
;
629 iatoms
= ilist_mol
[F_VSITE2
].iatoms
;
631 for (k
= 0; k
< nrvsite2
; k
+= 4)
633 vsite
= a_offset
+ iatoms
[k
+1]; /* the vsite */
634 ai
= a_offset
+ iatoms
[k
+2]; /* constructing atom */
635 aj
= a_offset
+ iatoms
[k
+3]; /* constructing atom */
636 if (ggrpnr(groups
, egcQMMM
, ai
) < (groups
->grps
[egcQMMM
].nr
-1) &&
637 (ggrpnr(groups
, egcQMMM
, aj
) >= (groups
->grps
[egcQMMM
].nr
-1)))
639 /* mark ai as frontier atom */
640 for (i
= 0; i
< qm_nr
; i
++)
642 if ( (qm_arr
[i
] == ai
) || (qm_arr
[i
] == vsite
) )
644 qr
->qm
[j
]->frontatoms
[i
] = TRUE
;
648 else if (ggrpnr(groups
, egcQMMM
, aj
) < (groups
->grps
[egcQMMM
].nr
-1) &&
649 (ggrpnr(groups
, egcQMMM
, ai
) >= (groups
->grps
[egcQMMM
].nr
-1)))
651 /* mark aj as frontier atom */
652 for (i
= 0; i
< qm_nr
; i
++)
654 if ( (qm_arr
[i
] == aj
) || (qm_arr
[i
] == vsite
))
656 qr
->qm
[j
]->frontatoms
[i
] = TRUE
;
664 if (qr
->QMMMscheme
!= eQMMMschemeoniom
)
667 /* standard QMMM, all layers are merged together so there is one QM
668 * subsystem and one MM subsystem.
669 * Also we set the charges to zero in the md->charge arrays to prevent
670 * the innerloops from doubly counting the electostatic QM MM interaction
673 alook
= gmx_mtop_atomlookup_init(mtop
);
675 for (k
= 0; k
< qm_nr
; k
++)
677 gmx_mtop_atomnr_to_atom(alook
, qm_arr
[k
], &atom
);
681 qr
->qm
[0] = mk_QMrec();
682 /* store QM atoms in the QMrec and initialise
684 init_QMrec(0, qr
->qm
[0], qm_nr
, qm_arr
, mtop
, ir
);
685 if (qr
->qm
[0]->bOPT
|| qr
->qm
[0]->bTS
)
687 for (i
= 0; i
< qm_nr
; i
++)
689 gmx_mtop_atomnr_to_atom(alook
, qm_arr
[i
], &atom
);
690 /* nbfp now includes the 6.0/12.0 derivative prefactors */
691 qr
->qm
[0]->c6
[i
] = C6(fr
->nbfp
, mtop
->ffparams
.atnr
, atom
->type
, atom
->type
)/c6au
/6.0;
692 qr
->qm
[0]->c12
[i
] = C12(fr
->nbfp
, mtop
->ffparams
.atnr
, atom
->type
, atom
->type
)/c12au
/12.0;
696 /* find frontier atoms and mark them true in the frontieratoms array.
698 for (i
= 0; i
< qm_nr
; i
++)
700 gmx_mtop_atomnr_to_ilist(alook
, qm_arr
[i
], &ilist_mol
, &a_offset
);
701 nrvsite2
= ilist_mol
[F_VSITE2
].nr
;
702 iatoms
= ilist_mol
[F_VSITE2
].iatoms
;
704 for (k
= 0; k
< nrvsite2
; k
+= 4)
706 vsite
= a_offset
+ iatoms
[k
+1]; /* the vsite */
707 ai
= a_offset
+ iatoms
[k
+2]; /* constructing atom */
708 aj
= a_offset
+ iatoms
[k
+3]; /* constructing atom */
709 if (ggrpnr(groups
, egcQMMM
, ai
) < (groups
->grps
[egcQMMM
].nr
-1) &&
710 (ggrpnr(groups
, egcQMMM
, aj
) >= (groups
->grps
[egcQMMM
].nr
-1)))
712 /* mark ai as frontier atom */
713 if ( (qm_arr
[i
] == ai
) || (qm_arr
[i
] == vsite
) )
715 qr
->qm
[0]->frontatoms
[i
] = TRUE
;
718 else if (ggrpnr(groups
, egcQMMM
, aj
) < (groups
->grps
[egcQMMM
].nr
-1) &&
719 (ggrpnr(groups
, egcQMMM
, ai
) >= (groups
->grps
[egcQMMM
].nr
-1)))
721 /* mark aj as frontier atom */
722 if ( (qm_arr
[i
] == aj
) || (qm_arr
[i
] == vsite
) )
724 qr
->qm
[0]->frontatoms
[i
] = TRUE
;
730 gmx_mtop_atomlookup_destroy(alook
);
732 /* MM rec creation */
734 mm
->scalefactor
= ir
->scalefactor
;
735 mm
->nrMMatoms
= (mtop
->natoms
)-(qr
->qm
[0]->nrQMatoms
); /* rest of the atoms */
739 { /* MM rec creation */
741 mm
->scalefactor
= ir
->scalefactor
;
746 /* these variables get updated in the update QMMMrec */
748 if (qr
->nrQMlayers
== 1)
750 /* with only one layer there is only one initialisation
751 * needed. Multilayer is a bit more complicated as it requires
752 * re-initialisation at every step of the simulation. This is due
753 * to the use of COMMON blocks in the fortran QM subroutines.
755 if (qr
->qm
[0]->QMmethod
< eQMmethodRHF
)
757 #ifdef GMX_QMMM_MOPAC
758 /* semi-empiprical 1-layer ONIOM calculation requested (mopac93) */
759 init_mopac(cr
, qr
->qm
[0], qr
->mm
);
761 gmx_fatal(FARGS
, "Semi-empirical QM only supported with Mopac.");
766 /* ab initio calculation requested (gamess/gaussian/ORCA) */
767 #ifdef GMX_QMMM_GAMESS
768 init_gamess(cr
, qr
->qm
[0], qr
->mm
);
769 #elif defined GMX_QMMM_GAUSSIAN
770 init_gaussian(cr
, qr
->qm
[0], qr
->mm
);
771 #elif defined GMX_QMMM_ORCA
772 init_orca(qr
->qm
[0]);
774 gmx_fatal(FARGS
, "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
780 void update_QMMMrec(t_commrec
*cr
,
787 /* updates the coordinates of both QM atoms and MM atoms and stores
788 * them in the QMMMrec.
790 * NOTE: is NOT yet working if there are no PBC. Also in ns.c, simple
791 * ns needs to be fixed!
794 mm_max
= 0, mm_nr
= 0, mm_nr_new
, i
, j
, is
, k
, shift
;
796 *mm_j_particles
= NULL
, *qm_i_particles
= NULL
;
812 *parallelMMarray
= NULL
;
816 c6au
= (HARTREE2KJ
*AVOGADRO
*pow(BOHR2NM
, 6));
817 c12au
= (HARTREE2KJ
*AVOGADRO
*pow(BOHR2NM
, 12));
819 /* every cpu has this array. On every processor we fill this array
820 * with 1's and 0's. 1's indicate the atoms is a QM atom on the
821 * current cpu in a later stage these arrays are all summed. indexes
822 * > 0 indicate the atom is a QM atom. Every node therefore knows
823 * whcih atoms are part of the QM subsystem.
825 /* copy some pointers */
828 QMMMlist
= fr
->QMMMlist
;
832 /* init_pbc(box); needs to be called first, see pbc.h */
833 set_pbc_dd(&pbc
, fr
->ePBC
, DOMAINDECOMP(cr
) ? cr
->dd
: NULL
, FALSE
, box
);
834 /* only in standard (normal) QMMM we need the neighbouring MM
835 * particles to provide a electric field of point charges for the QM
838 if (qr
->QMMMscheme
== eQMMMschemenormal
) /* also implies 1 QM-layer */
840 /* we NOW create/update a number of QMMMrec entries:
842 * 1) the shiftQM, containing the shifts of the QM atoms
844 * 2) the indexMM array, containing the index of the MM atoms
846 * 3) the shiftMM, containing the shifts of the MM atoms
848 * 4) the shifted coordinates of the MM atoms
850 * the shifts are used for computing virial of the QM/MM particles.
852 qm
= qr
->qm
[0]; /* in case of normal QMMM, there is only one group */
853 snew(qm_i_particles
, QMMMlist
.nri
);
856 qm_i_particles
[0].shift
= XYZ2IS(0, 0, 0);
857 for (i
= 0; i
< QMMMlist
.nri
; i
++)
859 qm_i_particles
[i
].j
= QMMMlist
.iinr
[i
];
863 qm_i_particles
[i
].shift
= pbc_dx_aiuc(&pbc
, x
[QMMMlist
.iinr
[0]],
864 x
[QMMMlist
.iinr
[i
]], dx
);
867 /* However, since nri >= nrQMatoms, we do a quicksort, and throw
868 * out double, triple, etc. entries later, as we do for the MM
872 /* compute the shift for the MM j-particles with respect to
873 * the QM i-particle and store them.
876 crd
[0] = IS2X(QMMMlist
.shift
[i
]) + IS2X(qm_i_particles
[i
].shift
);
877 crd
[1] = IS2Y(QMMMlist
.shift
[i
]) + IS2Y(qm_i_particles
[i
].shift
);
878 crd
[2] = IS2Z(QMMMlist
.shift
[i
]) + IS2Z(qm_i_particles
[i
].shift
);
879 is
= XYZ2IS(crd
[0], crd
[1], crd
[2]);
880 for (j
= QMMMlist
.jindex
[i
];
881 j
< QMMMlist
.jindex
[i
+1];
887 srenew(mm_j_particles
, mm_max
);
890 mm_j_particles
[mm_nr
].j
= QMMMlist
.jjnr
[j
];
891 mm_j_particles
[mm_nr
].shift
= is
;
896 /* quicksort QM and MM shift arrays and throw away multiple entries */
900 qsort(qm_i_particles
, QMMMlist
.nri
,
901 (size_t)sizeof(qm_i_particles
[0]),
903 qsort(mm_j_particles
, mm_nr
,
904 (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 */