Pair_style hybrid/overlay + HDNNP

Dear all,
I am trying to use hybrid/overlay pair_style with lj/charmmfsw/coul/long and hdnnp potentials. The system has 10 atom types and the hdnnp is used for the interaction between the last two (for which the model is trained). However, when I try to run the setup I get the following error:

Element mapping string from LAMMPS to n2p2: ",9:H,10:O"
terminate called after throwing an instance of 'std::runtime_error'
  what():  ERROR: Element mapping is inconsistent, NNP elements: 2, emap elements: 3.

Please, can you help me, what am I doing wrong?
Or is it a bug? It looks the same as the bug in this post: Support of multiple pair styles via hybrid/overlay in lammps interface · Issue #29 · CompPhysVienna/n2p2 · GitHub
I used the 2 August 2023 LAMMPS version. The setup file looks like:

variable T         equal 300.0
variable tau_T     equal 100.0
variable p         equal 1.0
variable tau_p     equal 1000.0
variable s_thermo  equal 10000
variable steps     equal 2000000

units              real

neighbor           1.0 bin
neigh_modify       every 1 delay 0 check yes

atom_style         full
bond_style         harmonic
improper_style     harmonic
angle_style        charmm
dihedral_style     charmm
special_bonds      charmm
pair_style         hybrid/overlay lj/charmmfsw/coul/long 8 10 & 
                   hdnnp 6.36 dir "../nnp" cflength 1.8897259 cfenergy 0.0015936

read_data          system.data

pair_coeff         1 1 lj/charmmfsw/coul/long 0.12 2.43572
pair_coeff         2 2 lj/charmmfsw/coul/long 0.15 4.04468
pair_coeff         3 3 lj/charmmfsw/coul/long 0.024 2.38761
pair_coeff         4 4 lj/charmmfsw/coul/long 0.046 0.400014
pair_coeff         5 5 lj/charmmfsw/coul/long 0.11 3.56359
pair_coeff         6 6 lj/charmmfsw/coul/long 0.01 3.38542
pair_coeff         7 7 lj/charmmfsw/coul/long 0.2 2.76179
pair_coeff         8 8 lj/charmmfsw/coul/long 0.12 2.49452

pair_coeff         1 9 lj/charmmfsw/coul/long 0.074 1.417867
pair_coeff         2 9 lj/charmmfsw/coul/long 0.083 2.222347
pair_coeff         3 9 lj/charmmfsw/coul/long 0.083 1.393812
pair_coeff         4 9 lj/charmmfsw/coul/long 0.046 0.400014
pair_coeff         5 9 lj/charmmfsw/coul/long 0.071 1.981802
pair_coeff         6 9 lj/charmmfsw/coul/long 0.021 1.892717
pair_coeff         7 9 lj/charmmfsw/coul/long 0.021 1.580902
pair_coeff         8 9 lj/charmmfsw/coul/long 0.021 1.580902

pair_coeff         1 10 lj/charmmfsw/coul/long 0.135 2.793145
pair_coeff         2 10 lj/charmmfsw/coul/long 0.151 3.597625
pair_coeff         3 10 lj/charmmfsw/coul/long 0.060 2.76909
pair_coeff         4 10 lj/charmmfsw/coul/long 0.0836 1.775292
pair_coeff         5 10 lj/charmmfsw/coul/long 0.129 3.35708
pair_coeff         6 10 lj/charmmfsw/coul/long 0.039 3.267995
pair_coeff         7 10 lj/charmmfsw/coul/long 0.174 2.95618
pair_coeff         8 10 lj/charmmfsw/coul/long 0.135 2.822545

pair_coeff         * * hdnnp NULL NULL NULL NULL NULL NULL NULL NULL H O

pair_modify        mix arithmetic
kspace_style       ewald 1e-6

thermo             ${s_thermo}
thermo_style       multi
timestep           0.5

fix                npt all npt temp ${T} ${T} ${tau_T} tchain 3 iso ${p} ${p} ${tau_p} mtk yes
velocity           all create ${T} 42 dist gaussian

run                ${steps}

Thank you,
Pavol

You should not be using pair_style hybrid/overlay for this but plain pair_style hybrid.
However, that is unrelated to the issue at hand.

Yes, it is a bug.

No, it is not the same bug. It is rather a bug in how the bugfix was implemented. It implicitly assumed that the first atom type was to be used with hdnnp.

Please try the following one-line change, recompile, and let us know if that makes a difference.

  diff --git a/src/ML-HDNNP/pair_hdnnp.cpp b/src/ML-HDNNP/pair_hdnnp.cpp
  index 35b0dd0ae8..cfe1e0bb64 100644
  --- a/src/ML-HDNNP/pair_hdnnp.cpp
  +++ b/src/ML-HDNNP/pair_hdnnp.cpp
  @@ -200,7 +200,7 @@ void PairHDNNP::coeff(int narg, char **arg)
     emap = "";
     for (int i = 2; i < narg; i++) {
       if (strcmp(arg[i], "NULL") != 0) {
  -      if (i != 2) emap += ",";
  +      if (!emap.empty()) emap += ",";
         emap += std::to_string(i - 1) + ":" + arg[i];
         map[i - 1] = 1;
       }

Thank you very much for such a quick response. I corrected and recompiled the code and now it works perfectly.
Also thank you for the other comment, corrected too.

1 Like

thanks for confirming.