LAMMPS and DREIDING hydrogen bonding potential: confusion remains

LAMMPS LAMMPS Mailing List Mirror
1

Dear Steve and Tod

Even with the most recent version of lammps (10th Feb 2012 checkout) with all the patches on this potential, I have the same asnwers as below except that now hydrogen bonding energy is -2.1879901 kcal/mol. So my question still remains:

Why coulomb and vander waals’ energy are changing when DREIDING hydrogen bonding potential is switched ON for the SAME atomic configuration? Shouldn’t they remain same since positions of atoms are exactly same in both the cases ?(hydrogen bonding potential ON and OFF)

Thanks
Chetan

2

Dear Steve and Tod

Even with the most recent version of lammps (10th Feb 2012 checkout) with
all the patches on this potential, I have the same asnwers as below except
that now hydrogen bonding energy is -2.1879901 kcal/mol. So my question
still remains:

Why coulomb and vander waals' energy are changing when DREIDING hydrogen
bonding potential is switched ON for the SAME atomic configuration?
Shouldn't they remain same since positions of atoms are exactly same in both
the cases ?(hydrogen bonding potential ON and OFF)

chetan,

please produce a simple test case.
it is difficult to argue without seeing
exactly what you are doing.

thanks,
    axel.

3

Hi Axel

Thanks. Here are the details that I wrote in my very first email:

I am running some calculations to output energy values for a same system (read as same atomic configurations) with and without DREIDING hydrogen bonding potential (henceforth denoted as hbp). When hbp is ON, lammps outputs total potential energy to be less than that when hbp OFF by 140 kcal/mol. As expected this difference comes from the difference in pair energy in the two cases. In pair energy, elong is same, whereas surprisingly ecoul is less by 111 kcal/mol and evdw less by 29 kcal/mol in case of hbp ON as compared to OFF case. There are two questions here:

  1. Since DREIDING 1990 article clearly says that all electrostatic and VdW interactions are to be included in hbp ON case, I wonder why there is such a difference with respect to these two energy values. Since atomic configurations are same, I would expect vander waals and coulomb contribution to remain same regardless of whether hydrogen bonding potential is on or off.

  2. I initially thought that difference in epair is due to difference in hydrogen bonding energy, however, that value (around 2 kcal/mol) does not seem to be contributing to epair, since epair is directly obtained by summing vander waals, coulomb’s and long contribution.

Could you clarify please? Thank you very much.

regards

Chetan

4

Hi Axel

Thanks. Here are the details that I wrote in my very first email:

I did read your emails, but you didn’t read mine. Please don’t give a description that requires people to first come up with an input by themselves, but provide an actual input deck that somebody can verify to be correct. Best is a very minimal input that contains just the minimum of commands to reproduce your problem.

Axel

5

Hi Axel

I have attached my input files herewith. As said, I am just printing out energy details when DREIDING hydrogen bonding potential is applied and then when that potential is turned off.

Thanks
Chetan

in.ABIm-both-90k (4.36 KB)

restartf2.9000000 (234 KB)

6

Dear All

I just noticed that lammps document says this hydrogen bonding potential does not support mixing, whereas I am using geometric mixing. I wonder if it is due to that. However, I quite vividly recall discussing with Steve that this mixing constraint was removed and also lammps did not give me any error while I used that.

Thanks
Chetan

7

All pair styles tally their energy into what
LAMMPS calls evdwl and ecoul for output with thermo.
So that output is always the total energy of the potential.
Whether it is really Vdwl and Coulombic energy or not.

Some potentials like Dreiding keep track of a further
breakdown, e.g, the HBE component you are asking about.

But you have to manually extract it from the potential via the compute pair
command and output it separately.

E.g. on the pair hbond doc page it says:

These pair styles tally a count of how many hydrogen bonding
interactions they calculate each timestep and the hbond energy. These
quantities can be accessed via the "compute pair"_compute_pair.html
command as a vector of values of length 2.

To print these quantities to the log file (with a descriptive column
heading) the following commands could be included in an input script:

compute hb all pair hbond/dreiding/lj
variable n_hbond equal c_hb\[1\] #number hbonds
variable E_hbond equal c_hb\[2\] #hbond energy
thermo_style custom step temp epair v_E_hbond

Steve

8

Dear All

I just noticed that lammps document says this hydrogen bonding potential
does not support mixing, whereas I am using geometric mixing. I wonder if it
is due to that. However, I quite vividly recall discussing with Steve that

yes. it is.

this mixing constraint was removed and also lammps did not give me any error
while I used that.

as far is i remember from the code, the situation is the following:
for hybrid pair styles mixing is supported, if the mixing is
required between pairs of atom types that all have the same
substyle and that substyle does support mixing. for hybrid/overlay
the rule is basically the same, but now you can have more than
one pair style contribute to a pair, thus making pairs that were
equal before, suddenly no longer equal.

i am attaching a variant of your input file that highlights that,
by trying to do the transition from your first configuration to
your second configuration in smaller steps (it is possible).
also, i don't rely on the system remaining the same, but use
clear and read in the restart again. if you want to compare
things to find the cause of a difference, you have to make
those as similar as possible. so you'll see now where the
change is coming from.

a) the hybrid/overlay setup with the hbond potential added

b) the same thing, but this time a morse potential (doesn't
support mixing!) with a well depth energy of 0.0 (= doesn't
add energy).

c) the hybrid/overlay setup but without a second potential
(yes, it is pointless in real life, but important in this context).

d) the single pair potential without hybrid.

comparing the non-bonded energies, this gives:

E_vdwl E_coul E_pair Press PotEng E_hbond
   372.81497 1831.9196 -1668.1814 -2703.5982 -962.82199 -2.1879901

in.ABIm-both-90k-ak (9.56 KB)

9

Hi Steve

Thanks for clarification. I already have printed that hydrogen bonding component. Then next question is:

That hydrogen bonding component of DREIDING should be equivalent to difference in VdW and Coulomb of two cases (hydrogen bonding potential (hbp) ON and OFF) for the same atomic configuration. When hbp is ON, lammps outputs total potential energy to be less than that when hbp OFF by 140 kcal/mol. As expected this difference comes from the difference in pair energy in the two cases. In pair energy, elong is same, whereas surprisingly ecoul is less by 111 kcal/mol and evdw less by 29 kcal/mol in case of hbp ON as compared to OFF case.

I initially thought that difference in above is due to difference in hydrogen bonding energy, however, that value (around 2 kcal/mol) does not seem to be contributing to difference in VdW and Coulomb above. So question is how to account for this difference in pair energy when hbp is ON as compared to OFF?

I am using geometric mixing and I wonder if it is due to that. I recall talking to you that you had removed that constraint from hydrogen bonding potential. However, lammps document on this pair style still says it does not support mixing. So now I am trying to test it with direct input of pair coefficients instead of using geometric mixing option.

Thanks
Chetan

10

Thanks Axel for shedding more light on hybrid Vs hybrid/overlay. As I said, I have now tested that after I added “I J coefficients” for pair interactions of atoms involved in hydrogen bonding separately with in addition to geometric mixing and that energies now tally perfectly . It proves conclusively that it was me using only geometric mixing while using hybrid/overlay that was giving all problems in energy tallying. I was under impression that geometric mixing was okay while using such hydrogen bonding potential (owing to some past conversations on this mailing lists long time back) and also lammps did not give me an error or warning.

Thanks
Chetan

11

Thanks Axel for shedding more light on hybrid Vs hybrid/overlay. As I said,
I have now tested that after I added "I J coefficients" for pair
interactions of atoms involved in hydrogen bonding separately with in
addition to geometric mixing and that energies now tally perfectly .
It proves conclusively that it was me using only geometric mixing while
using hybrid/overlay that was giving all problems in energy tallying. I was
under impression that geometric mixing was okay while using such hydrogen
bonding potential (owing to some past conversations on this mailing lists
long time back) and also lammps did not give me an error or warning.

given the nature of the issue, it is impossible for LAMMPS to tell
what mixing of potential parameters in hybrid/overlay potentials
is desired or not. the default is to mix whatever it can. this is
probably not a good choice for unsuspecting users, since it leads
to the inconsistency issues that you ran into. it is consistent for
the way how the hyrbid pair styles are implemented and it makes
sense when you know and understand the implementation.
so rather than adding all mixed parameter pairs, it should have
been sufficient to add explicitly mixed pair_coeff statements for
those atom types that are

the hybrid pair styles are both a very powerful and flexible tool,
but also easily a source of confusion and unexpected behavior.

steve,
with this experience in mind, it may be worth considering to add
a diagnostic to the LAMMPS output that notifies a user whenever
potential parameters are not taken directly from the input but
automatically generated using a mixing rule (and which one).
e.g.
pair coeffs for types 1 2 set from geometric mixing
pair coeffs for types 1 3 set from geometric mixing
and so on.

that might be useful feedback for unsuspecting users.

cheers,
     axel.

12

I haven't followed all the details of this discussion.
I don't particularly like the idea of printing this:

pair coeffs for types 1 2 set from geometric mixing

for each mixing, since for a 10 atom type problem,
you would get 100 messages.

My recollection is that mixing under pair hybrid
follows a simple rule, which is explained on the doc
page:

For atom type pairs I,J and I != J, if the sub-style assigned to I,I
and J,J is the same, and if the sub-style allows for mixing, then the
coefficients for I,J can be mixed. This means you do not have to
specify a pair_coeff command for I,J since the I,J type pair will be
assigned automatically to the I,I sub-style and its coefficients
generated by the mixing rule used by that sub-style. For the
{hybrid/overlay} style, there is an additional requirement that both
the I,I and J,J pairs are assigned to a single sub-style. See the
"pair_modify" command for details of mixing rules. See the See the
doc page for the sub-style to see if allows for mixing.

I.e. mixing takes place within a sub-style if it is supported
by that sub-style.

Is this not working correctly, or is there something confusing
about this approach?

Steve

13

Hi Steve

Now everything is clear. I am definitely using geometric mixing (instead of printing them individually) plus printing separately I!=J pair coefficients for atoms involved in DREIDING hydrogen bonding potential. I had not realized earlier to print this extra I!=J pair coefficients, since through some past discussions on the mailing lists, I was under impression that “no geometric mixing constraint” has been removed from DREIDING hydrogen bonding potential. I had ofcourse not noticed the updated lammps page that hydrogen bonding potential does not support mixing now.

Thank you all for your help,
cheers
Chetan