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