src/gromacs/mdlib/qmmm.c

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