[gmx-developers] GROMACS 4.0.7 QM/MM exclusions

[Lee-Ping]

Hi there,

My coworker Tim and I am doing a careful implementation of the Q-Chem
interface in GROMACS 4.0.7.  So far, we have added driver routines to
call Q-Chem and perform the I/O.  We've also modified the "simple"
neighbor search so that it can work with QM/MM.  We are currently in the
middle of testing the code with incrementally larger and more
complicated systems, and I've hit a few roadblocks along the way when
performing simulations on water clusters (either all QM, or part QM /
part MM).

The function "generate_qmexcl" in src/kernel/topio.c seems to have the
following problem.

	  if (molb->nmol > 1) {
	    /* Split the molblock after this molecule */
	    sys->nmolblock++;
988	    srenew(sys->molblock,sys->nmolblock);
	    for(i=mb; i<sys->nmolblock-1; i++) {
	      sys->molblock[i+1] = sys->molblock[i];
	    }
	    sys->molblock[mb  ].nmol  = 1;
	    sys->molblock[mb+1].nmol -= 1;
993         //molb = &sys->molblock[mb];
	  }

	  /* Add a moltype for the QMMM molecule */
	  sys->nmoltype++;
	  srenew(sys->moltype,sys->nmoltype);
	  /* Copy the moltype struct */
1000	  sys->moltype[sys->nmoltype-1] = sys->moltype[molb->type];
	  /* Copy the exclusions to a new array, since this is the only
	   * thing that needs to be modified for QMMM.
	   */
1004	  copy_blocka(&sys->moltype[molb->type     ].excls,
		      &sys->moltype[sys->nmoltype-1].excls);
	  /* Set the molecule type for the QMMM molblock */
	  molb->type = sys->nmoltype - 1;
	}
	generate_qmexcl_moltype(&sys->moltype[molb->type],grpnr,ir);

When I try to run grompp on a system of two water molecules (both in the
QM region, or one QM + one MM), I get a segmentation fault at line 1000.
When I compile with gdb and backtrace, molb->type is gibberish.  This
happens even when I compile the vanilla 4.0.7 distribution with
Gaussian, so I'm sure it's not a bug resulting from our implementation.

The problem must be coming from the "srenew" command at line 998.  I
think I've been able to fix the problem with the insertion at line 993;
this is copied from line 982 and restores the correct molb->type.  

However, there is still another problem.  The entire moltype struct is
not actually copied at line 1000; line 1004 copies the exclusions over
and the comment claims that "this is the only thing that needs to be
modified for QMMM."  But this is untrue; generate_qmexcl_moltype not
only modifies the exclusions, but also deletes entries from the bonded
interactions (ilist) within moltype!  

This will modify the bonded interactions for the MM molecules if they
were originally the same type as the QM molecules.  I have verified
this, and it's a serious problem that will break QM/MM calculations
where QM and MM molecules start from the same type (e.g. if they are
both water).

I think this is a bug that can be avoided by copying over the ilist
struct (or the entire moltype struct) when the QM/MM moltype is created,
but no such functions seem to be built-in to GROMACS 4.0.7.  I could
write these functions myself if people would find it useful.  Please
advise.

- Lee-Ping Wang