Hello again :),
I will copy the previous conversation so people can follow the complete discussion:
Jon:
I would like to ask you about another issue which is directly connected with this problem. We had made the same calculations employing the TIP4P/2005 water model. We used the EW3D summations method with fix Efield, substracting by hand the berkowitz correction. The results are surprisingly good if they are compared with other results we have based on DFT. As far as I understood (this was commented in another issue), when fix efield is employed with the TIP4P model, the forces are not correctly applied, instead of applying the force to the massless particle, it is done into the oxygen atom (which has no charge in this particular model), is this correct?
But if this is the case, it is easy to calculate that the torque that would feel each water molecule under the external electric field is 1,35 times higher than the one that would feel the correct TIP4P/2005 dipole with the same external electric field. Thereby, it would be as applying a 1,35 times higher external electrical field (torque = dipole_moment x External field) into the correct system. My big concern with these estimations is that the dynamics of our calculations are not correct, as we are applying the external force into an atom which should not feel it (although it is very close to the massles particle, which is the one that should!). In conclusion, I think that the calculations are not correct, but the results agree surprisingly good, which makes me feel hope :). I would really appreciate to know your opinion about this issue.
Sorry for bothering again, and thank you very much!
Jon
srtee:
To simulate a uniform electric field acting on the virtual charges of TIP4P molecules, you can manually do “field-shifting” by applying one fix efield on the TIP4P oxygens and another on the TIP4P hydrogens, of different field strength , so that the resultant molecular forces and torques are preserved. See pppm_tip4p.cpp:242-252 for some idea of the exact expressions. A third fix efield would have to be applied to non-TIP4P charges.
I do not see how a different external efield would solve the problem. Just to make the problem simple, we are only simulating water, so we have only the TIP4P/2005 molecule. If we apply different external electrical fields to the oxygen and hydrogen atoms, then the molecule will never feel a net external electric force which is zero. This is because all the charges of model sum up to zero. And thereby, I would expect that the dynamics of the system would differ even more with respect to the “real” one, as there would be a net force acting on each molecule, and therefore, directly modifying the motion of the center of mass of each molecule. Maybe there is something I am missing, if this is the case I would really appreciate your clarification or help.
Let me reformulate the particular questions that we have in a more direct way, i think they are going to be easy to answer for you:
1.- If fix efield is applied during a simulation with TIP4P model, the external force that in theory should be applied into the massless particle, in practice is applied into the oxygen atom. This sentence is true, right?
2.- If 1 is true, then the dynamics of the system are not correct, right?
3.- We are confining water into 1 nm scale slab, so the forces are quite strong along the confining direction. The massless particle is close to the oxygen atom, this could suggest that although the dynamics are not correct, the error introduced by applying the external efield force at different (but close) atoms could be small. Do you know if there is something we could do to estimate whether the error on the dynamics is small?
Thank you in advance, you guys do a great job for this community.
Jon
If I understand the source code correctly, the electric field will be applied to the oxygen atom. You can check that by setting up two simulations with fix efield; one with a single TIP4P molecule and a second one with a single four-atom rigid body, where the massless charge site is an atom with a really low mass (like 10-e4) and the oxygen atom itself does not posses any charge.
This way you should be able to estimate the influence of the shifted charge position.
TIP4P is a rigid molecule, so any combination of individual forces gives the same dynamics as long as the total molecular forces and torques are the same. But you don’t have to take my word for it – just check the source code and you will see for yourself that this is how the P3M part of it is handled.
If you scale the oxygen efield by 1-a (so downwards) and the hydrogen efield by 1+a, you can calculate for yourself that applying these at the oxygen and hydrogen positions gives the same molecular force and torque as applying your intended efield at the virtual site and hydrogens. Here 1-a is the ratio of virtual-site molecular dipole to oxygen-site molecular dipole.
Thank you so much, you are completely right.
Just to confirm my hope: The oxygen-site molecular dipole is 1,35 times higher than the virtual-site dipole. In our calculations we employed fix efield with 0.1 V/A. If I understand correctly your argument (I think I do and I agree) and given that the torque in each molecule is the dipole moment times the electric field, the dynamics of the system that we simulated is EXACTLY the same as doing it for the virtual-site dipole (correct TIP4P/2005) with an external field of 0.135 V/A. Is this correct?
Jon
You’re right, and in addition I was wrong about how to set up efield: the PPPM code deals directly with the force corrections, which are imposed in opposite directions on oxygen and hydrogen but therefore represent corrections to electric field in the same direction.
So simply downscaling the (real) electric field by the dipole ratio (and feeding that to efield) should give correct dynamics, equivalently your post is correct. In principle the correction is more involved for inhomogeneous fields – in practice a water molecule is small enough that they shouldn’t matter for external fields.
