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