src/gromacs/mdlib/qmmm.c

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