Dear mailing list,
I have written a potential within the atom_style charge.
I need to find the charges on each atom that minimizes the electrostatic energy.
I am trying to use the
fix qeq/reax
But I get the error:
ERROR: Must use pair_style reax/c with fix qeq/reax (…/fix_qeq_reax.cpp:322)
My input file is this (my questions below).
%%%%%%%%%%%%% start of input file %%%%%%%%%%%%%
units metal
boundary p p p
atom_style charge # dove si definisce qst ???
read_data Cell.now.data
Dear mailing list,
I have written a potential within the atom_style charge.
I need to find the charges on each atom that minimizes the electrostatic
energy.
I am trying to use the
fix qeq/reaxBut I get the error:
ERROR: Must use pair_style reax/c with fix qeq/reax
(../fix_qeq_reax.cpp:322)My input file is this (my questions below).
%%%%%%%%%%%%% start of input file %%%%%%%%%%%%%
units metal
boundary p p p
atom_style charge # dove si definisce qst ???
read_data Cell.now.data#-------------------------------------------------------------------------------------
# newton on # proviamo a metterlo
mass 1 26.98 # Al
mass 2 16.00 # O#---------------------------------------------------------------------------------
kspace_style pppm 1e-5
kspace_modify gewald 0.319182 # to modify g_ewald = 1 / sigma # units
of g_ewald = 1/distance#-------------------------------------------------------------------------------------
pair_style hybrid/overlay rfmeamlegfalse ctip 9.4
pair_coeff * * rfmeamlegfalse
../library.rfmeam.AlO.LazicThijsse.SET-05 Al O
pair_coeff * * ctip ../library.ctip.AlO.LazicThijsse.SET-04 Al Opair_modify table 0 # in order never to enter in tabinner
#-------------------------------------------------------------------------------------
# TRY here charge equilibration
# prova con valori a caso per taper radius : giusto per vedere se gira
fix 1 all qeq/reax 1 0.0 10.0 0.1 params_qeq_AlO_DS#----------------------------------------------------------------------------------------------------------------
# PRINT FORCES ON ATOMS
dump 1 all custom 1 Forces-data id fx fy fz type q #type # x y z xs ys zs
thermo_style custom etotal step fmax evdwl ecoul elong enthalpy#-------------------------------------------------------------------------------------------------------------------
# MINIMIZATION
min_style cg
thermo 1 # == print at every step
minimize 0.2 0.2 0 0#---------------------------------------------------------------------------------------------------------------
# Print the final cell
write_data data.*%%%%%%%%%%%%% end of input file %%%%%%%%%%%%%
I am using version Feb14 without any patch.
My questions are:
1) Is the fix qeq/reax usable for the reaxFF-like potentials only ?
Once you comment "if (!force->pair_match("reax/c",1)) error->all(FLERR,"Must
use pair_style reax/c with fix qeq/reax");" out in the init() function, you
can use it with your own pair_style provided that you supply the qeq
parameters.
Does this case include also the "atom_style full" as in:
http://lammps.sandia.gov/threads/msg40099.html ?
Or am I missing sth and it can be used also with atom_style charge ?
It has to be used with atom_style charge, but since atom_style full is a
superset of charge, fix qeq/reax can be used with full also.
2) Is there anything already implemented (or soon to come) in Lammps that
computes the equilibrium charges in a generic system ?
Yes, soon to come. But a user still needs to supply qeq parameters.
3) In case there is not, could you advice me some publications on how to
implement it ?
Rappe and Goddard; Rick, Stuart and Berne; Streitz and Mintmire; Zhou,
Wadley, Filhol and Neurock, etc.
What would you advice me to change (or to start from) in the file
fix_qeq_reax.cpp to exclude the dihedral and similar parts that are
excluded in my potential so to be able to use it ?
Fix qeq/reax is pure pair wise for all non-bonded pairs.
Ray
Dear Shan,
Thank you for your reply.
But I still have some problems when using it.
I have commented out the line you told me.
I am testing a simple bulk cell with 2 atoms that should have a charge of minimum energy around 0.8e.
Since it is small I am replicating it.
But regardless of the parameters that I use, I always get NaN for the charge with the fix.
Any advice of what I might be the problem ?
Below are my input files.
The code works correctly without the fix.
Should I include in my files also the derivative of the energy wrt the charge ?
If so what is the syntax and the general variable names to be used ?
A correlated question:
Does it make a difference if in my potential the electronegativity and the coulomb hardness is already included ?
Is it enough in this case to insert zero in the parameters file for the fix qeq ?
Thank you in advance for your help.
Bests,
Daniele Scopece
units metal
boundary p p p
atom_style charge
read_data Cell.now.data
replicate 3 3 3
Dear Shan,
Thank you for your reply.
But I still have some problems when using it.I have commented out the line you told me.
I am testing a simple bulk cell with 2 atoms that should have a charge of
minimum energy around 0.8e.
Since it is small I am replicating it.
If I understand you correctly, you have a small unitcell with 2 atoms, and
the equilibrium charges should be +0.8 and -0.8 for the two atoms.
But regardless of the parameters that I use, I always get NaN for the
charge with the fix.Any advice of what I might be the problem ?
There could be a lot of reasons. Your QEq parameters are not from 1994
Streitz-Mintmire nor 2004 Zhou for Al/O, so that could be a big problem.
Your pair_style ctip may be incorrect, etc.
Below are my input files.
The code works correctly without the fix.
Should I include in my files also the derivative of the energy wrt the
charge ?
No, why?
If so what is the syntax and the general variable names to be used ?
A correlated question:
Does it make a difference if in my potential the electronegativity and the
coulomb hardness is already included ?
You have to include those parameters in a file for fix qeq/reax to read.
But you don't need the third parameter (screening) for Streitz-Mintmire.
Your NANs are most likely resulted from bad QEq parameters.
Ray
Dear Shan,
Thank you for your reply.
My comments below
[...]
What do you mean for "incorrect" ? I checked it with another serial code
for fixed charges and it gives the same results.Why do you say that the parameters different from the ones you mention
could be a serious problem ?
Is the method tuned to work just with some set of parameters ?
Is it unstable with other parameters ?
I said "maybe" incorrect - without seeing the code nobody can be sure.
Fix qeq/reax needs the Coulomb screening parameter to actually work. A bad
third parameter can give you really bad results.
Mind you, you cannot mix and match two models and expect them to work right
off the shelf. Streitz-Mintmire's CTIP potential is meant to be used with
its QEq model, not Rappe&Goddard's original QEq model. If you are using
fix qeq/reax, you are using Rappe&Goddards' model, not CTIP's model. Then
you have to experiment with the third parameter to make it work.
That is all I can say. Besides, this is not a LAMMPS question anymore.
Ray