Hello, I found in lammps documents about command fix reax/c that “Use of this pair style requires that a charge be defined for every atom”. But the examples given by lammps all set charge as 0 for every atoms. I was confused, so is it means that the atoms charge don’t need to set in advance and will be distributed automatically by command fix reax/c when doing simulation ? Thank you.
The charge will be updated in every step through the charge equilibration fix (QEq). However, this procedure may get trapped in an unphysical state for some cases, so that a decent initial guess can help substantially to get correct results and speed up the first step.
This should be evident from studying the publications describing the ReaxFF method.
That is not correct, the data files in examples/reaxff/HNS and examples/reaxff/water do set charges != 0.0.
Dear akohlmey, many thanks for your generous and instant help !!
How to reasonably set the initial charge ? Thank your very much
The simplest way would be to take partial charges from non-polarizable classical force fields like CHARMM, Amber, OPLS, and so on. This doesn’t have to be very precise. The key point is to break the symmetry, i.e. make sure that atoms have either a positive or negative charge as expected. For example, for a water molecule, the oxygen atom should be negative and the hydrogen positive, but if you set them to -0.2 and 0.1, respectively or -0.8 and 0.4 initially should not make much of a difference. Important is that that are not 0.0. That helps the QEq algorithm a lot in not getting trapped in the wrong minimum.
1 Like
Got it, thank you again
Hi Axel,
I tried with my system with(non-zero) and without(all-zero) initial charge.
The log file gave exacly the same energy output. If the inital charge was different, I assume the initial potential energy has to be different.
Non-zero inital charge gives,
Step Temp Press KinEng PotEng TotEng
0 304.21579 25286.809 905.90267 -66506.359 -65600.456
10 1588.5172 15895.015 4730.3329 -70390.754 -65660.421
20 858.82628 11305.26 2557.438 -68289.665 -65732.227
30 701.53664 21764.819 2089.0563 -67835.757 -65746.7
40 1803.3259 21173.253 5369.9965 -71218.176 -65848.179
50 555.14885 14061.13 1653.1385 -67549.561 -65896.423
And all-zero initial charge gives,
Step Temp Press KinEng PotEng TotEng
0 304.21579 25286.809 905.90267 -66506.359 -65600.456
10 1588.5172 15895.015 4730.3329 -70390.754 -65660.421
20 858.82628 11305.26 2557.438 -68289.665 -65732.227
30 701.53664 21764.819 2089.0563 -67835.757 -65746.7
40 1803.3259 21173.253 5369.9965 -71218.176 -65848.179
50 555.14885 14061.13 1653.1385 -67549.561 -65896.423
Then I checked the output from fix reaxff/bonds command, it seems that no matter how I set the intial charge, reaxff potential sets the charge for each bond by its own, which is the same for my two tests.
So, based on my test result, does it mean that there is no need to set the inital charge manually when simulating with reaxff?
Thanks.
No. It only means that in this specific case the is no problem. It will be the same for a lot of systems, but not all. Setting initial changes is a protective measure.
I see.
Thank you sir.
If there is no classical force fields to obtain the initial charge of atoms, how to determinate the initial charge ?
You can look at electronegativity and give everything that is less electronegative than hydrogen or carbon a (small) positive charge and everything more electronegative a (small) negative charge while making sure that the total charge is 0. Or you can chance it and leave them all zero. As mentioned before, the motivation is the break any symmetries and thus steer the QEq minimizer into the right direction. It doesn’t have to be very precise.
thank you very much
Another reasonable approach is to use an automated force field parameter service, like the Automated Topology Builder, and start with those partial charges. (Whatever you do, document it, or you will have forgotten by the time you need to write your paper!)