src/gromacs/mdlib/qmmm.c

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