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