andysim / helpme Goto Github PK

Hi Andy,

I'm running into energy dependence on the kappa parameter when dipoles are present. I've noticed that helPME does not have kappasweep tests for this case so I'm wondering if you've verified it. When testing, I did add the E_self(mu) but it didn't help much. I also tried with a wide range of k-space grids and with a very large real-space cutoff, and both with zero and non-zero total dipole moment to no avail. Below is a table of one of the test runs for a water dimer from your unit tests (including coordinates, charges and dipole moments) with different kappas (called alpha). All quantities are in atomic units (Hartree and Bohr).

I would appreciate if you'd let me know if I'm missing something and maybe add a kappasweep test with dipoles.

Kind regards,
Alexei

Edit: the GridSp is approximate. It's an input parameter used to compute the grid dimensions.
Edit2: the splineOrder was 7 in all tests but I've also run with higher orders with the same results.

Hi Andy,
I have a question regarding reciprocal-space dipole contributions to the electric field (computePRec method). It's the last term in the Eq (49) https://doi.org/10.1063/1.3511713 for electric field due to dipole moments (I guess it's the field due to interaction the dipole's images in other cells?). I'm not an expert in the Ewald summation to judge if the articles cited in the helPME docs have or not the term but I did run some comparison between helPME and another (non-PME) Ewald implementation where this term can be switched off. The results without the term were identical to helPME, which makes me think that this term might be missing. Am I missing something or is there a mistake in the paper by Sala et. al?

Kind regards,
Alexei

Does this library provide the periodic effective potential matrix (i.e. the 1/r_{ij} matrix for molecular system) at each atomic/charge site in the system?

In other libraries e.g. OpenMM, sampling the potential (i.e. get q_j/r_{ij} at site i) at a charge site just returns infinity and it would be awesome to have this feature.

Many thanks!

Dear developers,

Shouldn't the sign of force be opposite to the default one?
The following code describes how I calculated numerical differentiation and compared it to the analytical one.

helpme::Matrix<double> coordsD(
    {{2.0, 2.0, 2.0}, {2.5, 2.0, 3.0}, {1.5, 2.0, 3.0}, {0.0, 0.0, 0.0}, {0.5, 0.0, 1.0}, {-0.5, 0.0, 1.0}});
helpme::Matrix<double> chargesD({-0.834, 0.417, 0.417, -0.834, 0.417, 0.417});
double scaleFactorD = 332.0716;

double energyD = 0;
helpme::Matrix<double> forcesD(6, 3);
helpme::Matrix<double> virialD(1, 6);
helpme::Matrix<double> potentialAndFieldD(6, 4);

auto pmeD = std::unique_ptr<PMEInstanceD>(new PMEInstanceD);
pmeD->setup(1, 0.3, 5, 32, 32, 32, scaleFactorD, 1);
pmeD->setLatticeVectors(20, 20, 20, 90, 90, 90, PMEInstanceD::LatticeType::XAligned);
// Compute the energy and forces
energyD = pmeD->computeEFRec(0, chargesD, coordsD, forcesD);

std::cout << "force_analytical: " << forcesD[0][0] << '\n';

double alpha = 1e-5;

coordsD[0][0] += alpha;

double ene_p = pmeD->computeEFRec(0, chargesD, coordsD, forcesD);

coordsD[0][0] -= 2. * alpha;

double ene_m = pmeD->computeEFRec(0, chargesD, coordsD, forcesD);

double numerical = -(ene_p - ene_m) * 0.5 / alpha;

std::cout << "force_numerical: " << numerical << '\n';`

The output is::

force_analytical: 0.603153 force_numerical: -0.603153

I'm sorry if I have overlooked something.

Best,
Koichi