Bug Summary

File:gromacs/mdlib/qmmm.c
Location:line 903, column 13
Description:Null pointer passed as an argument to a 'nonnull' parameter

Annotated Source Code

1/*
2 * This file is part of the GROMACS molecular simulation package.
3 *
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.
10 *
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.
15 *
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.
20 *
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.
25 *
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.
33 *
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.
36 */
37#ifdef HAVE_CONFIG_H1
38#include <config.h>
39#endif
40
41#include <math.h>
42#include "typedefs.h"
43#include "types/commrec.h"
44#include "macros.h"
45#include "gromacs/utility/smalloc.h"
46#include "physics.h"
47#include "macros.h"
48#include "gromacs/math/vec.h"
49#include "force.h"
50#include "invblock.h"
51#include "gromacs/fileio/confio.h"
52#include "names.h"
53#include "network.h"
54#include "pbc.h"
55#include "ns.h"
56#include "nrnb.h"
57#include "bondf.h"
58#include "mshift.h"
59#include "txtdump.h"
60#include "qmmm.h"
61#include <stdio.h>
62#include <string.h>
63#include "gromacs/utility/fatalerror.h"
64#include "typedefs.h"
65#include <stdlib.h>
66#include "mtop_util.h"
67
68
69/* declarations of the interfaces to the QM packages. The _SH indicate
70 * the QM interfaces can be used for Surface Hopping simulations
71 */
72#ifdef GMX_QMMM_GAMESS
73/* GAMESS interface */
74
75void
76init_gamess(t_commrec *cr, t_QMrec *qm, t_MMrec *mm);
77
78real
79call_gamess(t_commrec *cr, t_forcerec *fr,
80 t_QMrec *qm, t_MMrec *mm, rvec f[], rvec fshift[]);
81
82#elif defined GMX_QMMM_MOPAC
83/* MOPAC interface */
84
85void
86init_mopac(t_commrec *cr, t_QMrec *qm, t_MMrec *mm);
87
88real
89call_mopac(t_commrec *cr, t_forcerec *fr, t_QMrec *qm,
90 t_MMrec *mm, rvec f[], rvec fshift[]);
91
92real
93call_mopac_SH(t_commrec *cr, t_forcerec *fr, t_QMrec *qm,
94 t_MMrec *mm, rvec f[], rvec fshift[]);
95
96#elif defined GMX_QMMM_GAUSSIAN
97/* GAUSSIAN interface */
98
99void
100init_gaussian(t_commrec *cr, t_QMrec *qm, t_MMrec *mm);
101
102real
103call_gaussian_SH(t_commrec *cr, t_forcerec *fr, t_QMrec *qm,
104 t_MMrec *mm, rvec f[], rvec fshift[]);
105
106real
107call_gaussian(t_commrec *cr, t_forcerec *fr, t_QMrec *qm,
108 t_MMrec *mm, rvec f[], rvec fshift[]);
109
110#elif defined GMX_QMMM_ORCA
111/* ORCA interface */
112
113void
114init_orca(t_QMrec *qm);
115
116real
117call_orca(t_forcerec *fr, t_QMrec *qm,
118 t_MMrec *mm, rvec f[], rvec fshift[]);
119
120#endif
121
122
123
124
125/* this struct and these comparison functions are needed for creating
126 * a QMMM input for the QM routines from the QMMM neighbor list.
127 */
128
129typedef struct {
130 int j;
131 int shift;
132} t_j_particle;
133
134static int struct_comp(const void *a, const void *b)
135{
136
137 return (int)(((t_j_particle *)a)->j)-(int)(((t_j_particle *)b)->j);
138
139} /* struct_comp */
140
141static int int_comp(const void *a, const void *b)
142{
143
144 return (*(int *)a) - (*(int *)b);
145
146} /* int_comp */
147
148static int QMlayer_comp(const void *a, const void *b)
149{
150
151 return (int)(((t_QMrec *)a)->nrQMatoms)-(int)(((t_QMrec *)b)->nrQMatoms);
152
153} /* QMlayer_comp */
154
155real call_QMroutine(t_commrec gmx_unused__attribute__ ((unused)) *cr, t_forcerec gmx_unused__attribute__ ((unused)) *fr, t_QMrec gmx_unused__attribute__ ((unused)) *qm,
156 t_MMrec gmx_unused__attribute__ ((unused)) *mm, rvec gmx_unused__attribute__ ((unused)) f[], rvec gmx_unused__attribute__ ((unused)) fshift[])
157{
158 /* makes a call to the requested QM routine (qm->QMmethod)
159 * Note that f is actually the gradient, i.e. -f
160 */
161 real
162 QMener = 0.0;
163
164 /* do a semi-empiprical calculation */
165
166 if (qm->QMmethod < eQMmethodRHF && !(mm->nrMMatoms))
167 {
168#ifdef GMX_QMMM_MOPAC
169 if (qm->bSH)
170 {
171 QMener = call_mopac_SH(cr, fr, qm, mm, f, fshift);
172 }
173 else
174 {
175 QMener = call_mopac(cr, fr, qm, mm, f, fshift);
176 }
177#else
178 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c", 178, "Semi-empirical QM only supported with Mopac.");
179#endif
180 }
181 else
182 {
183 /* do an ab-initio calculation */
184 if (qm->bSH && qm->QMmethod == eQMmethodCASSCF)
185 {
186#ifdef GMX_QMMM_GAUSSIAN
187 QMener = call_gaussian_SH(cr, fr, qm, mm, f, fshift);
188#else
189 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c", 189, "Ab-initio Surface-hopping only supported with Gaussian.");
190#endif
191 }
192 else
193 {
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(cr, fr, qm, mm, f, fshift);
200#else
201 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c", 201, "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
202#endif
203 }
204 }
205 return (QMener);
206}
207
208void init_QMroutine(t_commrec gmx_unused__attribute__ ((unused)) *cr, t_QMrec gmx_unused__attribute__ ((unused)) *qm, t_MMrec gmx_unused__attribute__ ((unused)) *mm)
209{
210 /* makes a call to the requested QM routine (qm->QMmethod)
211 */
212 if (qm->QMmethod < eQMmethodRHF)
213 {
214#ifdef GMX_QMMM_MOPAC
215 /* do a semi-empiprical calculation */
216 init_mopac(cr, qm, mm);
217#else
218 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c", 218, "Semi-empirical QM only supported with Mopac.");
219#endif
220 }
221 else
222 {
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
229 init_orca(qm);
230#else
231 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c", 231, "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
232#endif
233 }
234} /* init_QMroutine */
235
236void update_QMMM_coord(rvec x[], t_forcerec *fr, t_QMrec *qm, t_MMrec *mm)
237{
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.
241 */
242 int
243 i;
244
245 /* shift the QM atoms into the central box
246 */
247 for (i = 0; i < qm->nrQMatoms; i++)
248 {
249 rvec_sub(x[qm->indexQM[i]], fr->shift_vec[qm->shiftQM[i]], qm->xQM[i]);
250 }
251 /* also shift the MM atoms into the central box, if any
252 */
253 for (i = 0; i < mm->nrMMatoms; i++)
254 {
255 rvec_sub(x[mm->indexMM[i]], fr->shift_vec[mm->shiftMM[i]], mm->xMM[i]);
256 }
257} /* update_QMMM_coord */
258
259static void punch_QMMM_excl(t_QMrec *qm, t_MMrec *mm, t_blocka *excls)
260{
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
264 */
265 FILE
266 *out = NULL((void*)0);
267 int
268 i, j, k, nrexcl = 0, *excluded = NULL((void*)0), max = 0;
269
270
271 out = fopen("QMMMexcl.dat", "w");
272
273 /* this can be done more efficiently I think
274 */
275 for (i = 0; i < qm->nrQMatoms; i++)
276 {
277 nrexcl = 0;
278 for (j = excls->index[qm->indexQM[i]];
279 j < excls->index[qm->indexQM[i]+1];
280 j++)
281 {
282 for (k = 0; k < mm->nrMMatoms; k++)
283 {
284 if (mm->indexMM[k] == excls->a[j]) /* the excluded MM atom */
285 {
286 if (nrexcl >= max)
287 {
288 max += 1000;
289 srenew(excluded, max)(excluded) = save_realloc("excluded", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 289, (excluded), (max), sizeof(*(excluded)));
290 }
291 excluded[nrexcl++] = k;
292 continue;
293 }
294 }
295 }
296 /* write to file: */
297 fprintf(out, "%5d %5d\n", i+1, nrexcl);
298 for (j = 0; j < nrexcl; j++)
299 {
300 fprintf(out, "%5d ", excluded[j]);
301 }
302 fprintf(out, "\n");
303 }
304 free(excluded);
305 fclose(out);
306} /* punch_QMMM_excl */
307
308
309/* end of QMMM subroutines */
310
311/* QMMM core routines */
312
313t_QMrec *mk_QMrec(void)
314{
315 t_QMrec *qm;
316 snew(qm, 1)(qm) = save_calloc("qm", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 316, (1), sizeof(*(qm)));
317 return qm;
318} /* mk_QMrec */
319
320t_MMrec *mk_MMrec(void)
321{
322 t_MMrec *mm;
323 snew(mm, 1)(mm) = save_calloc("mm", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 323, (1), sizeof(*(mm)));
324 return mm;
325} /* mk_MMrec */
326
327static void init_QMrec(int grpnr, t_QMrec *qm, int nr, int *atomarray,
328 gmx_mtop_t *mtop, t_inputrec *ir)
329{
330 /* fills the t_QMrec struct of QM group grpnr
331 */
332 int i;
333 gmx_mtop_atomlookup_t alook;
334 t_atom *atom;
335
336
337 qm->nrQMatoms = nr;
338 snew(qm->xQM, nr)(qm->xQM) = save_calloc("qm->xQM", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 338, (nr), sizeof(*(qm->xQM)));
339 snew(qm->indexQM, nr)(qm->indexQM) = save_calloc("qm->indexQM", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 339, (nr), sizeof(*(qm->indexQM)));
340 snew(qm->shiftQM, nr)(qm->shiftQM) = save_calloc("qm->shiftQM", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 340, (nr), sizeof(*(qm->shiftQM))); /* the shifts */
341 for (i = 0; i < nr; i++)
342 {
343 qm->indexQM[i] = atomarray[i];
344 }
345
346 alook = gmx_mtop_atomlookup_init(mtop);
347
348 snew(qm->atomicnumberQM, nr)(qm->atomicnumberQM) = save_calloc("qm->atomicnumberQM"
, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c", 348
, (nr), sizeof(*(qm->atomicnumberQM)));
349 for (i = 0; i < qm->nrQMatoms; i++)
350 {
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];
354 }
355
356 gmx_mtop_atomlookup_destroy(alook);
357
358 qm->QMcharge = ir->opts.QMcharge[grpnr];
359 qm->multiplicity = ir->opts.QMmult[grpnr];
360 qm->nelectrons -= ir->opts.QMcharge[grpnr];
361
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
374 * (init_gaussian()
375 */
376 /* print the current layer to allow users to check their input */
377 fprintf(stderrstderr, "Layer %d\nnr of QM atoms %d\n", grpnr, nr);
378 fprintf(stderrstderr, "QMlevel: %s/%s\n\n",
379 eQMmethod_names[qm->QMmethod], eQMbasis_names[qm->QMbasis]);
380
381 /* frontier atoms */
382 snew(qm->frontatoms, nr)(qm->frontatoms) = save_calloc("qm->frontatoms", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 382, (nr), sizeof(*(qm->frontatoms)));
383 /* Lennard-Jones coefficients */
384 snew(qm->c6, nr)(qm->c6) = save_calloc("qm->c6", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 384, (nr), sizeof(*(qm->c6)));
385 snew(qm->c12, nr)(qm->c12) = save_calloc("qm->c12", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 385, (nr), sizeof(*(qm->c12)));
386 /* do we optimize the QM separately using the algorithms of the QM program??
387 */
388 qm->bTS = ir->opts.bTS[grpnr];
389 qm->bOPT = ir->opts.bOPT[grpnr];
390
391} /* init_QMrec */
392
393t_QMrec *copy_QMrec(t_QMrec *qm)
394{
395 /* copies the contents of qm into a new t_QMrec struct */
396 t_QMrec
397 *qmcopy;
398 int
399 i;
400
401 qmcopy = mk_QMrec();
402 qmcopy->nrQMatoms = qm->nrQMatoms;
403 snew(qmcopy->xQM, qmcopy->nrQMatoms)(qmcopy->xQM) = save_calloc("qmcopy->xQM", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 403, (qmcopy->nrQMatoms), sizeof(*(qmcopy->xQM)));
404 snew(qmcopy->indexQM, qmcopy->nrQMatoms)(qmcopy->indexQM) = save_calloc("qmcopy->indexQM", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 404, (qmcopy->nrQMatoms), sizeof(*(qmcopy->indexQM)));
405 snew(qmcopy->atomicnumberQM, qm->nrQMatoms)(qmcopy->atomicnumberQM) = save_calloc("qmcopy->atomicnumberQM"
, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c", 405
, (qm->nrQMatoms), sizeof(*(qmcopy->atomicnumberQM)));
406 snew(qmcopy->shiftQM, qmcopy->nrQMatoms)(qmcopy->shiftQM) = save_calloc("qmcopy->shiftQM", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 406, (qmcopy->nrQMatoms), sizeof(*(qmcopy->shiftQM))); /* the shifts */
407 for (i = 0; i < qmcopy->nrQMatoms; i++)
408 {
409 qmcopy->shiftQM[i] = qm->shiftQM[i];
410 qmcopy->indexQM[i] = qm->indexQM[i];
411 qmcopy->atomicnumberQM[i] = qm->atomicnumberQM[i];
412 }
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;
427
428 /* Gaussian init. variables */
429 qmcopy->nQMcpus = qm->nQMcpus;
430 for (i = 0; i < DIM3; i++)
431 {
432 qmcopy->SHbasis[i] = qm->SHbasis[i];
433 }
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)(qmcopy->frontatoms) = save_calloc("qmcopy->frontatoms"
, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c", 438
, (qm->nrQMatoms), sizeof(*(qmcopy->frontatoms)));
439 snew(qmcopy->c12, qmcopy->nrQMatoms)(qmcopy->c12) = save_calloc("qmcopy->c12", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 439, (qmcopy->nrQMatoms), sizeof(*(qmcopy->c12)));
440 snew(qmcopy->c6, qmcopy->nrQMatoms)(qmcopy->c6) = save_calloc("qmcopy->c6", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 440, (qmcopy->nrQMatoms), sizeof(*(qmcopy->c6)));
441 if (qmcopy->bTS || qmcopy->bOPT)
442 {
443 for (i = 1; i < qmcopy->nrQMatoms; i++)
444 {
445 qmcopy->frontatoms[i] = qm->frontatoms[i];
446 qmcopy->c12[i] = qm->c12[i];
447 qmcopy->c6[i] = qm->c6[i];
448 }
449 }
450
451 return(qmcopy);
452
453} /*copy_QMrec */
454
455t_QMMMrec *mk_QMMMrec(void)
456{
457
458 t_QMMMrec *qr;
459
460 snew(qr, 1)(qr) = save_calloc("qr", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 460, (1), sizeof(*(qr)));
461
462 return qr;
463
464} /* mk_QMMMrec */
465
466void init_QMMMrec(t_commrec *cr,
467 gmx_mtop_t *mtop,
468 t_inputrec *ir,
469 t_forcerec *fr)
470{
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.
475 */
476
477 gmx_groups_t *groups;
478 atom_id *qm_arr = NULL((void*)0), vsite, ai, aj;
479 int qm_max = 0, qm_nr = 0, i, j, jmax, k, l, nrvsite2 = 0;
480 t_QMMMrec *qr;
481 t_MMrec *mm;
482 t_iatom *iatoms;
483 real c12au, c6au;
484 gmx_mtop_atomloop_all_t aloop;
485 t_atom *atom;
486 gmx_mtop_ilistloop_all_t iloop;
487 int a_offset;
488 t_ilist *ilist_mol;
489 gmx_mtop_atomlookup_t alook;
490
491 c6au = (HARTREE2KJ4.3597482e-21*AVOGADRO(6.0221367e23)*pow(BOHR2NM0.0529177249, 6));
492 c12au = (HARTREE2KJ4.3597482e-21*AVOGADRO(6.0221367e23)*pow(BOHR2NM0.0529177249, 12));
493 /* issue a fatal if the user wants to run with more than one node */
494 if (PAR(cr)((cr)->nnodes > 1))
495 {
496 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c", 496, "QM/MM does not work in parallel, use a single node instead\n");
497 }
498
499 /* Make a local copy of the QMMMrec */
500 qr = fr->qr;
501
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. */
506
507 qr->QMMMscheme = ir->QMMMscheme;
508
509 /* we take the possibility into account that a user has
510 * defined more than one QM group:
511 */
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.
515 */
516
517 /* small problem if there is only QM.... so no MM */
518
519 jmax = ir->opts.ngQM;
520
521 if (qr->QMMMscheme == eQMMMschemeoniom)
522 {
523 qr->nrQMlayers = jmax;
524 }
525 else
526 {
527 qr->nrQMlayers = 1;
528 }
529
530 groups = &mtop->groups;
531
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
536 * that.. 11-11-2003
537 */
538 snew(qr->qm, jmax)(qr->qm) = save_calloc("qr->qm", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 538, (jmax), sizeof(*(qr->qm)));
539 for (j = 0; j < jmax; j++)
540 {
541 /* new layer */
542 aloop = gmx_mtop_atomloop_all_init(mtop);
543 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
544 {
545 if (qm_nr >= qm_max)
546 {
547 qm_max += 1000;
548 srenew(qm_arr, qm_max)(qm_arr) = save_realloc("qm_arr", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 548, (qm_arr), (qm_max), sizeof(*(qm_arr)));
549 }
550 if (ggrpnr(groups, egcQMMM, i)((groups)->grpnr[egcQMMM] ? (groups)->grpnr[egcQMMM][i]
: 0) == j)
551 {
552 /* hack for tip4p */
553 qm_arr[qm_nr++] = i;
554 }
555 }
556 if (qr->QMMMscheme == eQMMMschemeoniom)
557 {
558 /* add the atoms to the bQMMM array
559 */
560
561 /* I assume that users specify the QM groups from small to
562 * big(ger) in the mdp file
563 */
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
570 * qm_arr[]. */
571
572 iloop = gmx_mtop_ilistloop_all_init(mtop);
573 while (gmx_mtop_ilistloop_all_next(iloop, &ilist_mol, &a_offset))
574 {
575 nrvsite2 = ilist_mol[F_VSITE2].nr;
576 iatoms = ilist_mol[F_VSITE2].iatoms;
577
578 for (k = 0; k < nrvsite2; k += 4)
579 {
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)((groups)->grpnr[egcQMMM] ? (groups)->grpnr[egcQMMM][vsite
] : 0) == ggrpnr(groups, egcQMMM, ai)((groups)->grpnr[egcQMMM] ? (groups)->grpnr[egcQMMM][ai
] : 0)
584 &&
585 ggrpnr(groups, egcQMMM, vsite)((groups)->grpnr[egcQMMM] ? (groups)->grpnr[egcQMMM][vsite
] : 0) == ggrpnr(groups, egcQMMM, aj)((groups)->grpnr[egcQMMM] ? (groups)->grpnr[egcQMMM][aj
] : 0))
586 {
587 /* this dummy link atom needs to be removed from the qm_arr
588 * before making the QMrec of this layer!
589 */
590 for (i = 0; i < qm_nr; i++)
591 {
592 if (qm_arr[i] == vsite)
593 {
594 /* drop the element */
595 for (l = i; l < qm_nr; l++)
596 {
597 qm_arr[l] = qm_arr[l+1];
598 }
599 qm_nr--;
600 }
601 }
602 }
603 }
604 }
605
606 /* store QM atoms in this layer in the QMrec and initialise layer
607 */
608 init_QMrec(j, qr->qm[j], qm_nr, qm_arr, mtop, ir);
609
610 /* we now store the LJ C6 and C12 parameters in QM rec in case
611 * we need to do an optimization
612 */
613 if (qr->qm[j]->bOPT || qr->qm[j]->bTS)
614 {
615 for (i = 0; i < qm_nr; i++)
616 {
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)(fr->nbfp)[2*((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)(fr->nbfp)[2*((mtop->ffparams.atnr)*(atom->type)+(atom
->type))+1]/c12au/12.0;
620 }
621 }
622 /* now we check for frontier QM atoms. These occur in pairs that
623 * construct the vsite
624 */
625 iloop = gmx_mtop_ilistloop_all_init(mtop);
626 while (gmx_mtop_ilistloop_all_next(iloop, &ilist_mol, &a_offset))
627 {
628 nrvsite2 = ilist_mol[F_VSITE2].nr;
629 iatoms = ilist_mol[F_VSITE2].iatoms;
630
631 for (k = 0; k < nrvsite2; k += 4)
632 {
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)->grpnr[egcQMMM] ? (groups)->grpnr[egcQMMM][ai
] : 0) < (groups->grps[egcQMMM].nr-1) &&
637 (ggrpnr(groups, egcQMMM, aj)((groups)->grpnr[egcQMMM] ? (groups)->grpnr[egcQMMM][aj
] : 0) >= (groups->grps[egcQMMM].nr-1)))
638 {
639 /* mark ai as frontier atom */
640 for (i = 0; i < qm_nr; i++)
641 {
642 if ( (qm_arr[i] == ai) || (qm_arr[i] == vsite) )
643 {
644 qr->qm[j]->frontatoms[i] = TRUE1;
645 }
646 }
647 }
648 else if (ggrpnr(groups, egcQMMM, aj)((groups)->grpnr[egcQMMM] ? (groups)->grpnr[egcQMMM][aj
] : 0) < (groups->grps[egcQMMM].nr-1) &&
649 (ggrpnr(groups, egcQMMM, ai)((groups)->grpnr[egcQMMM] ? (groups)->grpnr[egcQMMM][ai
] : 0) >= (groups->grps[egcQMMM].nr-1)))
650 {
651 /* mark aj as frontier atom */
652 for (i = 0; i < qm_nr; i++)
653 {
654 if ( (qm_arr[i] == aj) || (qm_arr[i] == vsite))
655 {
656 qr->qm[j]->frontatoms[i] = TRUE1;
657 }
658 }
659 }
660 }
661 }
662 }
663 }
664 if (qr->QMMMscheme != eQMMMschemeoniom)
665 {
666
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
671 */
672
673 alook = gmx_mtop_atomlookup_init(mtop);
674
675 for (k = 0; k < qm_nr; k++)
676 {
677 gmx_mtop_atomnr_to_atom(alook, qm_arr[k], &atom);
678 atom->q = 0.0;
679 atom->qB = 0.0;
680 }
681 qr->qm[0] = mk_QMrec();
682 /* store QM atoms in the QMrec and initialise
683 */
684 init_QMrec(0, qr->qm[0], qm_nr, qm_arr, mtop, ir);
685 if (qr->qm[0]->bOPT || qr->qm[0]->bTS)
686 {
687 for (i = 0; i < qm_nr; i++)
688 {
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)(fr->nbfp)[2*((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)(fr->nbfp)[2*((mtop->ffparams.atnr)*(atom->type)+(atom
->type))+1]/c12au/12.0;
693 }
694 }
695
696 /* find frontier atoms and mark them true in the frontieratoms array.
697 */
698 for (i = 0; i < qm_nr; i++)
699 {
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;
703
704 for (k = 0; k < nrvsite2; k += 4)
705 {
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)->grpnr[egcQMMM] ? (groups)->grpnr[egcQMMM][ai
] : 0) < (groups->grps[egcQMMM].nr-1) &&
710 (ggrpnr(groups, egcQMMM, aj)((groups)->grpnr[egcQMMM] ? (groups)->grpnr[egcQMMM][aj
] : 0) >= (groups->grps[egcQMMM].nr-1)))
711 {
712 /* mark ai as frontier atom */
713 if ( (qm_arr[i] == ai) || (qm_arr[i] == vsite) )
714 {
715 qr->qm[0]->frontatoms[i] = TRUE1;
716 }
717 }
718 else if (ggrpnr(groups, egcQMMM, aj)((groups)->grpnr[egcQMMM] ? (groups)->grpnr[egcQMMM][aj
] : 0) < (groups->grps[egcQMMM].nr-1) &&
719 (ggrpnr(groups, egcQMMM, ai)((groups)->grpnr[egcQMMM] ? (groups)->grpnr[egcQMMM][ai
] : 0) >= (groups->grps[egcQMMM].nr-1)))
720 {
721 /* mark aj as frontier atom */
722 if ( (qm_arr[i] == aj) || (qm_arr[i] == vsite) )
723 {
724 qr->qm[0]->frontatoms[i] = TRUE1;
725 }
726 }
727 }
728 }
729
730 gmx_mtop_atomlookup_destroy(alook);
731
732 /* MM rec creation */
733 mm = mk_MMrec();
734 mm->scalefactor = ir->scalefactor;
735 mm->nrMMatoms = (mtop->natoms)-(qr->qm[0]->nrQMatoms); /* rest of the atoms */
736 qr->mm = mm;
737 }
738 else /* ONIOM */
739 { /* MM rec creation */
740 mm = mk_MMrec();
741 mm->scalefactor = ir->scalefactor;
742 mm->nrMMatoms = 0;
743 qr->mm = mm;
744 }
745
746 /* these variables get updated in the update QMMMrec */
747
748 if (qr->nrQMlayers == 1)
749 {
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.
754 */
755 if (qr->qm[0]->QMmethod < eQMmethodRHF)
756 {
757#ifdef GMX_QMMM_MOPAC
758 /* semi-empiprical 1-layer ONIOM calculation requested (mopac93) */
759 init_mopac(cr, qr->qm[0], qr->mm);
760#else
761 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c", 761, "Semi-empirical QM only supported with Mopac.");
762#endif
763 }
764 else
765 {
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(cr, qr->qm[0], qr->mm);
773#else
774 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c", 774, "Ab-initio calculation only supported with Gamess, Gaussian or ORCA.");
775#endif
776 }
777 }
778} /* init_QMMMrec */
779
780void update_QMMMrec(t_commrec *cr,
781 t_forcerec *fr,
782 rvec x[],
783 t_mdatoms *md,
784 matrix box,
785 gmx_localtop_t *top)
786{
787 /* updates the coordinates of both QM atoms and MM atoms and stores
788 * them in the QMMMrec.
789 *
790 * NOTE: is NOT yet working if there are no PBC. Also in ns.c, simple
791 * ns needs to be fixed!
792 */
793 int
794 mm_max = 0, mm_nr = 0, mm_nr_new, i, j, is, k, shift;
795 t_j_particle
796 *mm_j_particles = NULL((void*)0), *qm_i_particles = NULL((void*)0);
1
'mm_j_particles' initialized to a null pointer value
→
797 t_QMMMrec
798 *qr;
799 t_nblist
800 QMMMlist;
801 rvec
802 dx, crd;
803 int
804 *MMatoms;
805 t_QMrec
806 *qm;
807 t_MMrec
808 *mm;
809 t_pbc
810 pbc;
811 int
812 *parallelMMarray = NULL((void*)0);
813 real
814 c12au, c6au;
815
816 c6au = (HARTREE2KJ4.3597482e-21*AVOGADRO(6.0221367e23)*pow(BOHR2NM0.0529177249, 6));
817 c12au = (HARTREE2KJ4.3597482e-21*AVOGADRO(6.0221367e23)*pow(BOHR2NM0.0529177249, 12));
818
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.
824 */
825 /* copy some pointers */
826 qr = fr->qr;
827 mm = qr->mm;
828 QMMMlist = fr->QMMMlist;
829
830
831
832 /* init_pbc(box); needs to be called first, see pbc.h */
833 set_pbc_dd(&pbc, fr->ePBC, DOMAINDECOMP(cr)(((cr)->dd != ((void*)0)) && ((cr)->nnodes >
1)) ? cr->dd : NULL((void*)0), FALSE0, 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
836 * atoms.
837 */
838 if (qr->QMMMscheme == eQMMMschemenormal) /* also implies 1 QM-layer */
2
←
Taking true branch
→
839 {
840 /* we NOW create/update a number of QMMMrec entries:
841 *
842 * 1) the shiftQM, containing the shifts of the QM atoms
843 *
844 * 2) the indexMM array, containing the index of the MM atoms
845 *
846 * 3) the shiftMM, containing the shifts of the MM atoms
847 *
848 * 4) the shifted coordinates of the MM atoms
849 *
850 * the shifts are used for computing virial of the QM/MM particles.
851 */
852 qm = qr->qm[0]; /* in case of normal QMMM, there is only one group */
853 snew(qm_i_particles, QMMMlist.nri)(qm_i_particles) = save_calloc("qm_i_particles", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 853, (QMMMlist.nri), sizeof(*(qm_i_particles)));
854 if (QMMMlist.nri)
3
←
Taking true branch
→
855 {
856 qm_i_particles[0].shift = XYZ2IS(0, 0, 0)((2*2 +1)*((2*1 +1)*((0)+1)+(0)+1)+(0)+2);
857 for (i = 0; i < QMMMlist.nri; i++)
4
←
Loop condition is false. Execution continues on line 900
→
858 {
859 qm_i_particles[i].j = QMMMlist.iinr[i];
860
861 if (i)
862 {
863 qm_i_particles[i].shift = pbc_dx_aiuc(&pbc, x[QMMMlist.iinr[0]],
864 x[QMMMlist.iinr[i]], dx);
865
866 }
867 /* However, since nri >= nrQMatoms, we do a quicksort, and throw
868 * out double, triple, etc. entries later, as we do for the MM
869 * list too.
870 */
871
872 /* compute the shift for the MM j-particles with respect to
873 * the QM i-particle and store them.
874 */
875
876 crd[0] = IS2X(QMMMlist.shift[i])(((QMMMlist.shift[i]) % (2*2 +1)) - 2) + IS2X(qm_i_particles[i].shift)(((qm_i_particles[i].shift) % (2*2 +1)) - 2);
877 crd[1] = IS2Y(QMMMlist.shift[i])((((QMMMlist.shift[i]) / (2*2 +1)) % (2*1 +1)) - 1) + IS2Y(qm_i_particles[i].shift)((((qm_i_particles[i].shift) / (2*2 +1)) % (2*1 +1)) - 1);
878 crd[2] = IS2Z(QMMMlist.shift[i])((QMMMlist.shift[i]) / ((2*2 +1)*(2*1 +1)) - 1) + IS2Z(qm_i_particles[i].shift)((qm_i_particles[i].shift) / ((2*2 +1)*(2*1 +1)) - 1);
879 is = XYZ2IS(crd[0], crd[1], crd[2])((2*2 +1)*((2*1 +1)*((crd[2])+1)+(crd[1])+1)+(crd[0])+2);
880 for (j = QMMMlist.jindex[i];
881 j < QMMMlist.jindex[i+1];
882 j++)
883 {
884 if (mm_nr >= mm_max)
885 {
886 mm_max += 1000;
887 srenew(mm_j_particles, mm_max)(mm_j_particles) = save_realloc("mm_j_particles", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 887, (mm_j_particles), (mm_max), sizeof(*(mm_j_particles)));
888 }
889
890 mm_j_particles[mm_nr].j = QMMMlist.jjnr[j];
891 mm_j_particles[mm_nr].shift = is;
892 mm_nr++;
893 }
894 }
895
896 /* quicksort QM and MM shift arrays and throw away multiple entries */
897
898
899
900 qsort(qm_i_particles, QMMMlist.nri,
901 (size_t)sizeof(qm_i_particles[0]),
902 struct_comp);
903 qsort(mm_j_particles, mm_nr,
5
←
Null pointer passed as an argument to a 'nonnull' parameter
904 (size_t)sizeof(mm_j_particles[0]),
905 struct_comp);
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.
909 */
910 j = 0;
911 for (i = 0; i < QMMMlist.nri; i++)
912 {
913 if (i == 0 || qm_i_particles[i].j != qm_i_particles[i-1].j)
914 {
915 qm_i_particles[j++] = qm_i_particles[i];
916 }
917 }
918 mm_nr_new = 0;
919 if (qm->bTS || qm->bOPT)
920 {
921 /* only remove double entries for the MM array */
922 for (i = 0; i < mm_nr; i++)
923 {
924 if ((i == 0 || mm_j_particles[i].j != mm_j_particles[i-1].j)
925 && !md->bQM[mm_j_particles[i].j])
926 {
927 mm_j_particles[mm_nr_new++] = mm_j_particles[i];
928 }
929 }
930 }
931 /* we also remove mm atoms that have no charges!
932 * actually this is already done in the ns.c
933 */
934 else
935 {
936 for (i = 0; i < mm_nr; i++)
937 {
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])))
942 {
943 mm_j_particles[mm_nr_new++] = mm_j_particles[i];
944 }
945 }
946 }
947 mm_nr = mm_nr_new;
948 /* store the data retrieved above into the QMMMrec
949 */
950 k = 0;
951 /* Keep the compiler happy,
952 * shift will always be set in the loop for i=0
953 */
954 shift = 0;
955 for (i = 0; i < qm->nrQMatoms; i++)
956 {
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.
960 */
961 if (qm->indexQM[i] == qm_i_particles[k].j)
962 {
963 shift = qm_i_particles[k++].shift;
964 }
965 /* use previous shift, assuming they belong the same charge
966 * group anyway,
967 */
968
969 qm->shiftQM[i] = shift;
970 }
971 }
972 /* parallel excecution */
973 if (PAR(cr)((cr)->nnodes > 1))
974 {
975 snew(parallelMMarray, 2*(md->nr))(parallelMMarray) = save_calloc("parallelMMarray", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 975, (2*(md->nr)), sizeof(*(parallelMMarray)));
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.
980 */
981 for (i = 0; i < 2*(md->nr); i++)
982 {
983 parallelMMarray[i] = 0;
984 }
985
986 for (i = 0; i < mm_nr; i++)
987 {
988 parallelMMarray[mm_j_particles[i].j] = 1;
989 parallelMMarray[mm_j_particles[i].j+(md->nr)] = mm_j_particles[i].shift;
990 }
991 gmx_sumi(md->nr, parallelMMarray, cr);
992 mm_nr = 0;
993
994 mm_max = 0;
995 for (i = 0; i < md->nr; i++)
996 {
997 if (parallelMMarray[i])
998 {
999 if (mm_nr >= mm_max)
1000 {
1001 mm_max += 1000;
1002 srenew(mm->indexMM, mm_max)(mm->indexMM) = save_realloc("mm->indexMM", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 1002, (mm->indexMM), (mm_max), sizeof(*(mm->indexMM))
);
1003 srenew(mm->shiftMM, mm_max)(mm->shiftMM) = save_realloc("mm->shiftMM", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 1003, (mm->shiftMM), (mm_max), sizeof(*(mm->shiftMM))
);
1004 }
1005 mm->indexMM[mm_nr] = i;
1006 mm->shiftMM[mm_nr++] = parallelMMarray[i+md->nr]/parallelMMarray[i];
1007 }
1008 }
1009 mm->nrMMatoms = mm_nr;
1010 free(parallelMMarray);
1011 }
1012 /* serial execution */
1013 else
1014 {
1015 mm->nrMMatoms = mm_nr;
1016 srenew(mm->shiftMM, mm_nr)(mm->shiftMM) = save_realloc("mm->shiftMM", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 1016, (mm->shiftMM), (mm_nr), sizeof(*(mm->shiftMM)));
1017 srenew(mm->indexMM, mm_nr)(mm->indexMM) = save_realloc("mm->indexMM", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 1017, (mm->indexMM), (mm_nr), sizeof(*(mm->indexMM)));
1018 for (i = 0; i < mm_nr; i++)
1019 {
1020 mm->indexMM[i] = mm_j_particles[i].j;
1021 mm->shiftMM[i] = mm_j_particles[i].shift;
1022 }
1023
1024 }
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
1028 */
1029 srenew(mm->xMM, mm->nrMMatoms)(mm->xMM) = save_realloc("mm->xMM", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 1029, (mm->xMM), (mm->nrMMatoms), sizeof(*(mm->xMM
)));
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
1033 */
1034 srenew(mm->MMcharges, mm->nrMMatoms)(mm->MMcharges) = save_realloc("mm->MMcharges", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 1034, (mm->MMcharges), (mm->nrMMatoms), sizeof(*(mm->
MMcharges)));
1035 for (i = 0; i < mm->nrMMatoms; i++) /* no free energy yet */
1036 {
1037 mm->MMcharges[i] = md->chargeA[mm->indexMM[i]]*mm->scalefactor;
1038 }
1039 if (qm->bTS || qm->bOPT)
1040 {
1041 /* store (copy) the c6 and c12 parameters into the MMrec struct
1042 */
1043 srenew(mm->c6, mm->nrMMatoms)(mm->c6) = save_realloc("mm->c6", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 1043, (mm->c6), (mm->nrMMatoms), sizeof(*(mm->c6))
);
1044 srenew(mm->c12, mm->nrMMatoms)(mm->c12) = save_realloc("mm->c12", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 1044, (mm->c12), (mm->nrMMatoms), sizeof(*(mm->c12
)));
1045 for (i = 0; i < mm->nrMMatoms; i++)
1046 {
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]])(fr->nbfp)[2*((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]])(fr->nbfp)[2*((top->idef.atnr)*(md->typeA[mm->indexMM
[i]])+(md->typeA[mm->indexMM[i]]))+1]/c12au/12.0;
1050 }
1051 punch_QMMM_excl(qr->qm[0], mm, &(top->excls));
1052 }
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
1055 * calculated above.
1056 */
1057
1058 update_QMMM_coord(x, fr, qr->qm[0], qr->mm);
1059 free(qm_i_particles);
1060 free(mm_j_particles);
1061 }
1062 else /* ONIOM */ /* ????? */
1063 {
1064 mm->nrMMatoms = 0;
1065 /* do for each layer */
1066 for (j = 0; j < qr->nrQMlayers; j++)
1067 {
1068 qm = qr->qm[j];
1069 qm->shiftQM[0] = XYZ2IS(0, 0, 0)((2*2 +1)*((2*1 +1)*((0)+1)+(0)+1)+(0)+2);
1070 for (i = 1; i < qm->nrQMatoms; i++)
1071 {
1072 qm->shiftQM[i] = pbc_dx_aiuc(&pbc, x[qm->indexQM[0]], x[qm->indexQM[i]],
1073 dx);
1074 }
1075 update_QMMM_coord(x, fr, qm, mm);
1076 }
1077 }
1078} /* update_QMMM_rec */
1079
1080
1081real calculate_QMMM(t_commrec *cr,
1082 rvec x[], rvec f[],
1083 t_forcerec *fr)
1084{
1085 real
1086 QMener = 0.0;
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
1090 */
1091 t_QMMMrec
1092 *qr;
1093 t_QMrec
1094 *qm, *qm2;
1095 t_MMrec
1096 *mm = NULL((void*)0);
1097 rvec
1098 *forces = NULL((void*)0), *fshift = NULL((void*)0),
1099 *forces2 = NULL((void*)0), *fshift2 = NULL((void*)0); /* needed for multilayer ONIOM */
1100 int
1101 i, j, k;
1102 /* make a local copy the QMMMrec pointer
1103 */
1104 qr = fr->qr;
1105 mm = qr->mm;
1106
1107 /* now different procedures are carried out for one layer ONION and
1108 * normal QMMM on one hand and multilayer oniom on the other
1109 */
1110 if (qr->QMMMscheme == eQMMMschemenormal || qr->nrQMlayers == 1)
1111 {
1112 qm = qr->qm[0];
1113 snew(forces, (qm->nrQMatoms+mm->nrMMatoms))(forces) = save_calloc("forces", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 1113, ((qm->nrQMatoms+mm->nrMMatoms)), sizeof(*(forces
)));
1114 snew(fshift, (qm->nrQMatoms+mm->nrMMatoms))(fshift) = save_calloc("fshift", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 1114, ((qm->nrQMatoms+mm->nrMMatoms)), sizeof(*(fshift
)));
1115 QMener = call_QMroutine(cr, fr, qm, mm, forces, fshift);
1116 for (i = 0; i < qm->nrQMatoms; i++)
1117 {
1118 for (j = 0; j < DIM3; j++)
1119 {
1120 f[qm->indexQM[i]][j] -= forces[i][j];
1121 fr->fshift[qm->shiftQM[i]][j] += fshift[i][j];
1122 }
1123 }
1124 for (i = 0; i < mm->nrMMatoms; i++)
1125 {
1126 for (j = 0; j < DIM3; j++)
1127 {
1128 f[mm->indexMM[i]][j] -= forces[qm->nrQMatoms+i][j];
1129 fr->fshift[mm->shiftMM[i]][j] += fshift[qm->nrQMatoms+i][j];
1130 }
1131
1132 }
1133 free(forces);
1134 free(fshift);
1135 }
1136 else /* Multi-layer ONIOM */
1137 {
1138 for (i = 0; i < qr->nrQMlayers-1; i++) /* last layer is special */
1139 {
1140 qm = qr->qm[i];
1141 qm2 = copy_QMrec(qr->qm[i+1]);
1142
1143 qm2->nrQMatoms = qm->nrQMatoms;
1144
1145 for (j = 0; j < qm2->nrQMatoms; j++)
1146 {
1147 for (k = 0; k < DIM3; k++)
1148 {
1149 qm2->xQM[j][k] = qm->xQM[j][k];
1150 }
1151 qm2->indexQM[j] = qm->indexQM[j];
1152 qm2->atomicnumberQM[j] = qm->atomicnumberQM[j];
1153 qm2->shiftQM[j] = qm->shiftQM[j];
1154 }
1155
1156 qm2->QMcharge = qm->QMcharge;
1157 /* this layer at the higher level of theory */
1158 srenew(forces, qm->nrQMatoms)(forces) = save_realloc("forces", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 1158, (forces), (qm->nrQMatoms), sizeof(*(forces)));
1159 srenew(fshift, qm->nrQMatoms)(fshift) = save_realloc("fshift", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 1159, (fshift), (qm->nrQMatoms), sizeof(*(fshift)));
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);
1163
1164 /* this layer at the lower level of theory */
1165 srenew(forces2, qm->nrQMatoms)(forces2) = save_realloc("forces2", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 1165, (forces2), (qm->nrQMatoms), sizeof(*(forces2)));
1166 srenew(fshift2, qm->nrQMatoms)(fshift2) = save_realloc("fshift2", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 1166, (fshift2), (qm->nrQMatoms), sizeof(*(fshift2)));
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
1172 */
1173 for (i = 0; i < qm->nrQMatoms; i++)
1174 {
1175 for (j = 0; j < DIM3; j++)
1176 {
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]);
1179 }
1180 }
1181 free(qm2);
1182 }
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)(forces) = save_realloc("forces", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 1186, (forces), (qm->nrQMatoms), sizeof(*(forces)));
1187 srenew(fshift, qm->nrQMatoms)(fshift) = save_realloc("fshift", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/qmmm.c"
, 1187, (fshift), (qm->nrQMatoms), sizeof(*(fshift)));
1188 QMener += call_QMroutine(cr, fr, qm, mm, forces, fshift);
1189 for (i = 0; i < qm->nrQMatoms; i++)
1190 {
1191 for (j = 0; j < DIM3; j++)
1192 {
1193 f[qm->indexQM[i]][j] -= forces[i][j];
1194 fr->fshift[qm->shiftQM[i]][j] += fshift[i][j];
1195 }
1196 }
1197 free(forces);
1198 free(fshift);
1199 free(forces2);
1200 free(fshift2);
1201 }
1202 if (qm->bTS || qm->bOPT)
1203 {
1204 /* qm[0] still contains the largest ONIOM QM subsystem
1205 * we take the optimized coordiates and put the in x[]
1206 */
1207 for (i = 0; i < qm->nrQMatoms; i++)
1208 {
1209 for (j = 0; j < DIM3; j++)
1210 {
1211 x[qm->indexQM[i]][j] = qm->xQM[i][j];
1212 }
1213 }
1214 }
1215 return(QMener);
1216} /* calculate_QMMM */
1217
1218/* end of QMMM core routines */