Fix qeq/reax with atom_style charge: Possible ?

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/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

#-------------------------------------------------------------------------------------
  # 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 O

  pair_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

1

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