Problem in splitting dispersion coefficients for LB mixing rule for long range dispersion

Dear all,

I am doing an NVT simulation with 2000 Tip4p/2005 water molecule in presence of 2 stationary atomistic surfaces having 2016 and 572 atoms respectively. The system is non-periodic in the z-direction with two surfaces placed at the top and bottom of the simulation box. I am using hybrid/overlay pair style with lj/long/tip4p/long for water and a surface interaction and lj/sf for water and other surface interaction. I am using LB (arithmetic) mixing rule.

The simulation blows with the error “ERROR on proc 0: Out of range atoms - cannot compute PPPMDisp (…/pppm_disp.cpp:4258)”. I looked into the PPPM dispersion parameters, there are “nan” values for all the accuracies, and dispersion grid is small.

The forcefield parameters are-

pair_style hybrid/overlay lj/long/tip4p/long long long 2 3 1 1 0.1546e-10 16.5e-10 lj/sf 6.0e-10

pair_modify mix arithmetic

pair_coeff 1 3 lj/long/tip4p/long 0.0 0.0

pair_coeff 1 4 lj/long/tip4p/long 0.0 0.0

pair_coeff 3 3 lj/long/tip4p/long 0.0 0.0

pair_coeff 2 3 lj/long/tip4p/long 0.0 0.0

pair_coeff 3 4* lj/long/tip4p/long 0.0 0.0

pair_coeff 2 2 lj/long/tip4p/long 1.28676442064E-21 3.1589e-10

pair_coeff 2 4 lj/sf 1.4e-21 3.5e-10 6.0e-10

pair_coeff 2 4 lj/long/tip4p/long 0.0 0.0

pair_coeff 1 2 lj/long/tip4p/long 6.948111596014770E-22 3.526921654E-10

pair_coeff 2 5 lj/long/tip4p/long 9.118247568989207E-22 3.90855875E-10

pair_coeff 4 4* lj/long/tip4p/long 0.0 0.0

pair_coeff 1 1 lj/long/tip4p/long 3.751755486576454E-22 3.894943308E-10

pair_coeff 1 5 lj/long/tip4p/long 4.923558706877879E-22 4.276580404E-10

pair_coeff 5 5 lj/long/tip4p/long 6.461356670712493E-22 4.6582175E-10

kspace_style pppm/disp/tip4p 5.000E-07

kspace_modify slab 3.0

kspace_modify force/disp/real 6.94786e-15

kspace_modify force/disp/kspace 1.389572e-13

PPPM parameters printed in the log file-

Optimizing splitting of Dispersion coefficients

Using 0 instead of 7 structure factors

Using 0 structure factors

Coulomb G vector (1/distance)= 2.06856e+09

Coulomb grid = 27 30 450

Coulomb stencil order = 5

Coulomb estimated absolute RMS force accuracy = 1.522e-14

Coulomb estimated relative force accuracy = 6.59705e-07

using double precision FFTs

3d grid and FFT values/proc = 569874 364500

Dispersion G vector (1/distance)= 3.03723e-43

Dispersion grid = 2 2 2

Dispersion stencil order = 5

Dispersion estimated absolute RMS force accuracy = -nan

Dispersion estimated absolute real space RMS force accuracy = -nan

Dispersion estimated absolute kspace RMS force accuracy = 0

Dispersion estimated relative force accuracy = -nan

I did troubleshooting and found that I am getting “nan” for eigenvalues and eigenvector for this system from the QR algorithm for eigenvalue decomposition. So, just to check I computed eigenvalues and eigenvectors using MATLAB and hard coded into the “pppm_disp.cpp”. Using these values, I am able to run the simulation without error and get correct results. There might be a bug in the code. Could someone look into this issue?

If required, I can provide the input script and data files.

Regards,