src/gromacs/mdlib/qmmm.cpp

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