Hi LAMMPS Users,
I saw in the manual that fix qeq/reax command "is typically used in conjunction with the ReaxFF force field model [...], but it can be used with any potential in LAMMPS, so long as it defines and uses charges on each atom."
I tried to use for simulating MoS2 (using QEq parameters taken from [Dellavalle et al., Langmuir 2012, 28, 7393-7400]) with atom_style charge and I got "ERROR: Must use pair_style reax/c with fix qeq/reax (../fix_qeq_reax.cpp:322)"
I found a message in the mailing list (http://ehc.ac/p/lammps/mailman/message/31696023/), explaining which lines should be commented for allowing fix qeq/reax with any pair styles. Thus I did this and then I tried to minimize the crystal structure and it worked. But I got atomic charges significantly lower than in the original paper (for Mo 0.25 against 0.74 in the paper). NVT simulations gave me a slightly larger value (0.30) but still half of the reference.
Moreover, I tried several times to perform MD simulations of sliding and/or starting from distorted configurations and I got always fatal errors "ERROR: Lost atoms: original 1152 current 1150 (../thermo.cpp:392)".
Anyone has any experience with fix qeq/reax without ReaxFF and/or with this force field and could give me some any suggestion?
Thanks!
Paolo
One way to troubleshoot your problem is to implement your system in GULP and compare the results to those you get from lammps. QEq is implemented in gulp and you can only feed the final structure to GULP and run QEq on it for a quick charge comparison exercise. Perhaps not the quickest solution but having the numbers from another code will help you attain more robust conclusions.
Carlos
I’m CCing Ray b/c I don’t think that line in the documentation
is correct. I believe It cannot (currently) be used with anything
but pair reax/c. Ray is working on a more general QeQ, but
it has not been released yet.
Steve
Yes, a more generalized fix qeq command will allow variable charge equilibration with any pair style with charges, though it must be used with caution. This work is in progress.
Ray
Hi LAMMPS Users,
I saw in the manual that fix qeq/reax command "is typically used in
conjunction with the ReaxFF force field model [...], but it can be used
with any potential in LAMMPS, so long as it defines and uses charges on
each atom."
I tried to use for simulating MoS2 (using QEq parameters taken from
[Dellavalle et al., Langmuir 2012, 28, 7393-7400]) with atom_style
charge and I got "ERROR: Must use pair_style reax/c with fix qeq/reax
(../fix_qeq_reax.cpp:322)"
I found a message in the mailing list
(http://ehc.ac/p/lammps/mailman/message/31696023/), explaining which
lines should be commented for allowing fix qeq/reax with any pair
styles. Thus I did this and then I tried to minimize the crystal
structure and it worked. But I got atomic charges significantly lower
than in the original paper (for Mo 0.25 against 0.74 in the paper). NVT
simulations gave me a slightly larger value (0.30) but still half of the
reference.
are you sure that the units are set correctly? please note that units for
fix qeq/reax are hardcoded...
Moreover, I tried several times to perform MD simulations of sliding
and/or starting from distorted configurations and I got always fatal
errors "ERROR: Lost atoms: original 1152 current 1150 (../thermo.cpp:392)".
charge equilibration and other kinds of ways to include polarization into
classical MD force fields are a tricky thing, especially for dynamics,
since you can easily induce unphysical overpolarization. thus you often
need to reduce the time step significantly.
Anyone has any experience with fix qeq/reax without ReaxFF and/or with
this force field and could give me some any suggestion?
not anything specific, but i've implemented a polarizable water potential
as an undergrad into a local MD code and i've had a hell of a time to get
the system to behave reasonably.
axel.
Dear Carlos, Steve, Ray and Axel,
thank you a lot for your replies.
@Carlos: Thank you for your suggestion... I will give a try with GULP.
@Steve: That line in documentation is actually misleading... Please let me know when QEq with general potentials will be implemented (even in alpha versions...)
@ Axel: - quite sure about units… I already asked this point (http://ehc.ac/p/lammps/mailman/lammps-users/thread/[email protected].../) to the mailing list… Anyway, units are eV, eV and Ang
- I reduced the timestep to 0.1 fs and I still got problems… Should I go down further to lower values?
Best regards,
Paolo
Dear Carlos, Steve, Ray and Axel,
thank you a lot for your replies.
@Carlos: Thank you for your suggestion... I will give a try with GULP.
@Steve: That line in documentation is actually misleading... Please let me
know when QEq with general potentials will be implemented (even in alpha
versions...)
@ Axel: - quite sure about units... I already asked this point (
http://ehc.ac/p/lammps/mailman/lammps-users/thread/[email protected].../)
to the mailing list... Anyway, units are eV, eV and Ang
- I reduced the timestep to 0.1 fs and I still got problems... Should I go
down further to lower values?
it is pointless to play with the timestep, if you cannot reproduce static
results. if you are dead certain about having entered the correct units,
you may need to check for any differences in how parameters are entered,
since the reference is likely done with a different code. there may be
factors of 2 or similar that are either included or not included in
parameters (cf. harmonic force constants)
axel.
Dear Carlos, Steve, Ray and Axel,
thank you a lot for your replies.
@Carlos: Thank you for your suggestion... I will give a try with GULP.
@Steve: That line in documentation is actually misleading... Please let me
know when QEq with general potentials will be implemented (even in alpha
versions...)
@ Axel: - quite sure about units... I already asked this point (
http://ehc.ac/p/lammps/mailman/lammps-users/thread/[email protected].../)
to the mailing list... Anyway, units are eV, eV and Ang
- I reduced the timestep to 0.1 fs and I still got problems... Should I go
down further to lower values?
Not that just that, the unit style must be in real units too.
Also, when you vary the charges of your original potential that is fitted
with fixed charges (albeit zero, fractional or formal charges) using a QEq
method, you are straying away from the original parameterization into some
unknown territory and it is not guaranteed you get reasonable and physical
answers. Even with small timesteps and reasonable QEq convergence, you may
still suffer instability resulted from charges.
Ray
@Steve: That line in documentation is actually misleading… Please let me know when QEq with general potentials will be >
implemented (even in alpha versions…)
The doc line is correct - I misunderstood it. You can use fix qeq/reax with other pair styles,
but you have to use the filename arg, not “reax/c”.
Steve
Going through reference papers I saw that there are differences in the implementation (I guess). The results I want to reproduce are obtained using QEq method as in Rappe and Goddard paper, that is calculating energy with E(Q) = E0 + chiQ + 1/2JQ^2. In LAMMPS, I guess, reference should be [A. Nakano, Comp. Phys. Comm., 104, 59-69 (1997)] where energy is calculated with the same formula. Thus, chi and J should be exactly the same. I think the point is the charge distributions. In QEq they are modeled with ns Slater-type orbitals ( Nr^(n-1)exp(-zetar) ), while in Nakano any atom is represented by a 1s orbital (N’exp(-2zeta*r)). Thus I think I can’t go ahead… I will try with GULP… Thanks, Paolo
This was OK. I am not sure to understand your point. Anyway, in principle it should work, since the potential I want to use was intended to use with QEq. I mean, someone else had already used this potential with QEq getting physical results. Thus, about this, I should be on the safe side… Paolo
Going through reference papers I saw that there are differences in the
implementation (I guess). The results I want to reproduce are obtained
using QEq method as in Rappe and Goddard paper, that is calculating energy
with E(Q) = E0 + chi*Q + 1/2*J*Q^2. In LAMMPS, I guess, reference should be
[A. Nakano, Comp. Phys. Comm., 104, 59-69 (1997)] where energy is
calculated with the same formula. Thus, chi and J should be exactly the
same.
I think the point is the charge distributions. In QEq they are modeled
with ns Slater-type orbitals ( N*r^(n-1)*exp(-zeta*r) ), while in Nakano
any atom is represented by a 1s orbital (N'*exp(-2*zeta*r)). Thus I think I
can't go ahead...
I am pretty sure they use the same 7th order taper function to approximate
the Slater-type orbitals, hence they are the same.
There may be other reasons as to why you obtained different charges -
structure, QEq parameters, shielding parameters, etc.
Ray
I am getting more confused... If one can use fix qeq/reax with any pair styles, why there is an error message (line 327, fix_qeq_reax.cpp) preventing this possibility in the current distribution of LAMMPS...?
Paolo
This line should be removed. You were correct in commenting this line to enable fix qeq/reax with other pair styles.
Ray
yes, you can remove that one (will be in next patch)
there is a similar error message at line 162 that is correct.
It is for when you specify “reax/c” instead of NULL.
That’s the one I thought you were talking about.
Steve