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