Implementing a new command

Dear Shan,

Thank you for your suggestion and reply.

I have another very rapid question about the implementation of QEq in LAMMPS.

It seems that the fictitious charges on the atoms (s and t) are initialized to zero through their history (s_hist and t_hist).

So the initial charges written in the cell are not considered as the starting point in the CG for the charge.

Is my understanding correct ?

I guess that this is done because of the difficulties in creating fictitious charges s and t from the real ones that satisfy their requirements and equations.

Is that correct ?

Thank you
Bests,
Daniele Scopece

And, in that case, could you please tell me

what is the name in the code of the variable “Step” of the outer loop that is printed out in your examples ?

(the one defined in “thermo_style custom step”)

I have analyzed the files in QEQ folder, but maybe this is defined in some other functions external to this folder …

Is it a global variable ?

I haven’t clear in mind all the flow of lammps yet, just some pieces so far …

I need it since I need to set to zero by hand the atom->q just for step 0, without the need that the user nullifies them in the cell file.
Indeed for the term I need to insert I have to compare the real q with qmin and qmax.

Thank you very much in advance for your assistance.

Bests,
Daniele Scopece

Dear Shan,

Thank you for your suggestion and reply.

I have another very rapid question about the implementation of QEq in
LAMMPS.

It seems that the fictitious charges on the atoms (s and t) are
initialized to zero through their history (s_hist and t_hist).
So the initial charges written in the cell are not considered as the
starting point in the CG for the charge.
Is my understanding correct ?

I guess that this is done because of the difficulties in creating
fictitious charges s and t from the real ones that satisfy their
requirements and equations.
Is that correct ?

Yes.

Ray

And, in that case, could you please tell me
what is the name in the code of the variable "Step" of the outer loop that
is printed out in your examples ?
(the one defined in "thermo_style custom step")
I have analyzed the files in QEQ folder, but maybe this is defined in some
other functions external to this folder ...
Is it a global variable ?
I haven't clear in mind all the flow of lammps yet, just some pieces so
far ...

It is the ntimestep in the update class.

Ray

Dear Shan and Mailing list,

I have other 2 questions concerning the structure of the QEQ and how it is implemented.

(1)

I would have expected, as inferred by an inspection of the source code, that the fix is not altering the energy dictated by the pair_style declaration,
i.e. that the additional energetic terms used in the fix (electronegativity & co) were used just in the fix.

Instead in your examples and others I get different values of the energy printed with and without the fix, even though I impose the same charges in the cell obtained from the charge equilibration.

See attached files: WithFix…* and WithoutFix…* with the input cell data.aC.cut (energy = -18 and -28 eV in the two cases).

So this fix is altering the energy of the atoms ?

Is there a way to print out the energy dictated just by the pair_style declaration and leave these additions just when the fix is acting ?

If yes, how ?

(2)

A more tricky question.

In the term that I have to implement (image U-limiting.png attached).

This is needed to prevent the coulombic-catastrophe when atoms are too close so to make disadvantageous to change the charges out of the limits [qmin:qmax].

You suggested me to develop the square and then to insert in chi and eta.

It is smart, but there is a problem.

This term depends on the real charges q_i :
1-(q_i-qmin)/|qi-qmin| that is either 0 or 2 depending on q_i.

But the two CG of the fictitious charges (s and t) are solved separately one after the other.

So the real charges are not known during the minimization.
And after these CGs the system is converged.

I guess that I need to solve the equations self consistently in the following way:

0 - set q_i = 0 (from histories of s and t, previous mails)

1 - compute the matrices with init_matvec

2 - solve CG for s with the real q from point 0

3 - solve CG for t with the real q from point 0

4 - obtain the real converged q

5 - go back to point 1 with this term updated and repeat until q is converged considering also this “limiting” term.

Can this work, considering the structure of the code ?
I guess that this way of proceeding is slowing down the process considerably.

Is there another way / does anybody have any suggestion on how to do that better ?

Is there a variant of the method of the fictitious charges that is valid also when the matrix is q-dependent.
Please let me know if I wasn’t clear enough.

Thank you in advance for your help and suggestion.

Best regards,
Daniele Scopece

data.aC.cut (752 Bytes)

WithFix…in.qeq.buck (1.08 KB)

WithFix…out.qeq.buck (1.14 KB)

WithoutFix…in.qeq.buck (1.09 KB)

WithoutFix…out.qeq.buck (1.08 KB)

U-limiting.png

Hi Scopece,

Existing QEq pair_styles (ReaxFF and COMB) already include self ionization energy (chiq+0.5eta*q^2) in their pair_styles so that fix qeq does not include this energy. So yes, fix qeq does not change energies, but just calculates new charges.

Hi Scopece,

Existing QEq pair_styles (ReaxFF and COMB) already include self ionization
energy (chi*q+0.5*eta*q^2) in their pair_styles so that fix qeq does not
include this energy. So yes, fix qeq does not change energies, but just
calculates new charges.

Dear Shan and Mailing list,

I have other 2 questions concerning the structure of the QEQ and how it
is implemented.

(1)
I would have expected, as inferred by an inspection of the source code,
that the fix is not altering the energy dictated by the pair_style
declaration,
i.e. that the additional energetic terms used in the fix
(electronegativity & co) were used just in the fix.
Instead in your examples and others I get different values of the energy
printed with and without the fix, even though I impose the same charges in
the cell obtained from the charge equilibration.
See attached files: WithFix...* and WithoutFix...* with the input cell
data.aC.cut (energy = -18 and -28 eV in the two cases).
So this fix is altering the energy of the atoms ?
Is there a way to print out the energy dictated just by the pair_style
declaration and leave these additions just when the fix is acting ?
If yes, how ?

Fix qeq does not change the energy. However, there is a mistake in the
code when it handles the kspace contribution - kspace energy is calculated
using the initial charges even after fix qeq changes the partial charge on
atoms. I will fix this bug.

Meanwhile, if you use a pair_style without a ksapce_style, you will see
identical energies.

Thanks to Axel and Steve, this bug is now fixed in the 8 Oct 14
distribution.

Cheers,
Ray

Dear Shan and mailing list,

Thank you for working on it.

I have installed the version of 9Oct14 from the tar download website and tested it with your example.

I hate to be the one pointing out things, but …
apparently there is still a discrepancy between with and without fix (though much lower than before):

-18.98001 eV vs -19.256448 eV in this example

It seems not to be determined by the cutoff in the fix (see attachments).

Should we consider this discrepancy as “necessary” or did I impose sth wrong ?

Thank you
Bests,
Daniele Scopece

data.aC.cut (800 Bytes)

WithFix…in.qeq.buck (1.08 KB)

With-Fix-141009…01.cutoff-10…out.qeq.buck (1.12 KB)

With-Fix-141009…02.cutoff-100…in.qeq.buck (1.15 KB)

With-Fix-141009…02.cutoff-100…out.qeq.buck (1.14 KB)

With-Fix-141009…03.cutoff-1000…in.qeq.buck (1.15 KB)

With-Fix-141009…03.cutoff-1000…out.qeq.buck (1.14 KB)

WithoutFix…in.qeq.buck (1.09 KB)

Without-Fix-141009…out.qeq.buck (1.08 KB)

Hi Daniele,

Just pointing out that in your tests, there are twice more neighbors when using fix qeq/slater.

Cheers,
Arthur

Scopece,

There is too much files to sort. Can you just post two input files where you think they should generate the same energies but they don’t?

Ray

Scopece,

I hate to be the one telling you that you are wrong again, but … apparently there is something you have to figure out first.

Have you tried this new test with a FULLY PERIODIC system? Say data.aC? Do you get the same energy with and without fix qeq? I bet you have not tested it, and I bet you will get the same energy.

Now, have you visualized the 8 atom data.aC.cut structure of yours? How are the charges distributed? Are all the Si atoms uniformly charged at ? Apparently not. What effect does it have on Coulombic energy? Please do examine it more carefully.

Cheers,
Ray