"Deviation too large for solvent" in NEB calculation
Clicked A Few Times
9:04:58 PM PST - Sun, Dec 18th 2011
Hi everyone...

I'm trying to run a NEB calculation on a solvated system with 52 atoms in the QM region; it does the DFT calculation and solvent optimization for several of the beads and then crashes with:
53:53: *  53: Deviation too large for solvent  483 *:: 483
(rank:53 hostname:mpc0894 pid:4688):ARMCI DASSERT fail. armci.c:ARMCI_Error():260 cond:0
Last System Error Message from Task 53:: Numerical result out of range

It looks like this happens during a solvent minimization phase that has run longer than the previous ones (347 steps), but not as long as specified by maxiter. My input script is:
start react

permanent_dir /home/Ueda-DNA/QBIC/neb/perm
scratch_dir /home/Ueda-DNA/QBIC/neb/scratch

charge -1

basis 
  C library 6-31G*
  H library 6-31G*
  N library 6-31G*
  O library 6-31G*
end

dft
  xc b3lyp
  direct
  iterations 500
end

md
  system react_neb
  msa 100
  cutoff 1.0
end

qmmm
  bqzone 10.0
  region solvent
  eatoms -1373.8747207761
  maxiter 200
  ncycles 1
  density espfit
  xyz neb_out
end

set qmmm:neb_path_limits react_opt.rst prod_ref.rst
 
set neb:nbeads 20
set neb:stepsize 1.0
set neb:steps 10

task qmmm dft neb

I've seen this problem described in the context of QM/MM optimization (http://www.emsl.pnl.gov/docs/nwchem/nwchem-support/2001/10/0008.Re:_deviation_too_large_fo...), where the recommendation was to try an MM optimization of the starting structure, to arrive at a more reasonable solvent configuration. I'd actually done an extensive QM/MM optimization of my starting coordinates, but (if I understand correctly) the initial bead configurations are a linear interpolation of the qmlink region between the starting and final coordinates, combined with the solvent coordinates from the initial configuration. So, especially for the beads closer to the final coordinates, the initial solvent/solute interactions are inevitably going to be horrible, and the beads near the middle of the trajectory will involve fairly large motions of the atoms in the qm region during optimization.

This seems to be a really common error -- I've probably tried this run about a dozen times now with different initial/final positions, trying out different options, etc. and this is usually what brings it down. Is there a good way to avoid this? My best guess at this point is to try and cook up a set of restart files that won't move as much during minimization rather than using qmmm:neb_path_limits to do the linear interpolation for me.

Any ideas would help a lot.

Thanks,

--craig

Clicked A Few Times
5:33:31 PM PST - Thu, Dec 22nd 2011
I would first check the pdb structures for your NEB beads for any obvious problems. To do this just run prepare command over numbered rst files in your perm directory

...

prepare

  read ./perm/sys_neb002.rst

  write sys_neb002.pdb

end

task prepare

Depending on what you see, you could run a few constant energy dynamics steps on just solvent part for the offending bead(s). Once you fix these problems you could "restart" your NEB calculations using these new rst files.

Regarding the solvent configuration here is how the code works. For new NEB calculation:
1. The solvent configuration is taken from the common rst file specified by the "system" directive in the MD block
(in your case react_neb.rst)
2. The QM region geometries for end points are extracted from the two rst files specified by set qmmm:neb_path_limits directive ( in your case react_opt.rst prod_ref.rst). Note that none of the solvent information is used from these files.
3. The interpolation is done to generate the intermediate QM region structures
4. In the first gradient calculation for the i-th bead the code check for the presence of rst file named sys_neb00i.rst in the perm directory. If it is not there (as typical for new NEB calculation) this rst file is created by taking common rst file (see step 1) and replacing its QM region with that of the i-th bead. Using the resulting rst file, the solvent optimization is performed around the fixed QM region, QM gradient is calculated, and the updated rst file is saved to the perm directory to be used in the next round of NEB calculation. And so on.

So the initial solvent configuration can be affected in two ways - either through common rst file (which will apply the same solvent configuration to all beads) and by depositing numbered rst files into the perm directory prior to calculation to affect individual beads. So for example, you could decide that reactant solvent configuration will work for first 5 beads, but that the product solvent configuration would be more appropriate to use for last 5 beads in you calculation. To accomplish this you copy react_opt.rst into common rst file, react_neb.rst, in the top run directory, and copy prod_ref.rst into sys_neb006.rst, ...., sys_neb010.rst in the perm directory.

Marat

Clicked A Few Times
7:31:22 AM PDT - Mon, Aug 27th 2012
learn more about NEB and its implemenation
Thanks, Marat

Clicked A Few Times
7:50:59 AM PDT - Mon, Aug 27th 2012
Are the product QM struectures replaced by the QM interpolation structure?
Hi Marat,

According your above explaination, for the rst files sys_neb006.rst, ...., sys_neb010.rst, the QM structure in these files are replaced by the interpolation QM structure?

Thanks,

wjb0920

Clicked A Few Times
2:11:39 AM PST - Wed, Mar 6th 2013
neb restart
Hi,

If we set "set qmmm:neb_path_limits react_opt.rst prod_ref.rst", then NWChem will rewrite the qm region coordinates for the depositing rst files sys_neb006.rst, ...., sys_neb010.rst (copy from prod_ref.rst) and reserve the mm region coordinates from product rst file.

If we do not set "set qmmm:neb_path_limits react_opt.rst prod_ref.rst" this line in input file, but specify want to do neb calculation, then the depositing rst file will be found by NWChem as the restart file. Therefore, we can restart or continue to run neb calculation.

Both of all facilitate neb calculation. A good design. Thanks.


Forum >> NWChem's corner >> QMMM