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