I’m running a molecular dynamics (MD) simulation of a bulk NaCl aqueous solution under an external electric field using LAMMPS, and I’ve encountered some unexpected behavior that I’d like to get feedback on.
Simulation Setup
System: Bulk NaCl aqueous solution
Boundary Conditions:Periodic in all directions
Units: metal
Potential:TIP4P/2005 for water + OPLS_AA for ions
Equilibration and Production Steps:
NPT ensemble: 2 ns (for density equilibration)
NVT + efield: 1 ns (apply external electric field)
NVE + efield: 3 ns (to observe ion dynamics without thermostat)
Electric Field: Applied along the z-direction using (fix 3 all efield 0.0 0.0 0.01).
What I Expected
Under the applied field, I expected:
Na⁺ ions to drift in one direction (along the electric field).
Cl⁻ ions to drift in the opposite direction.
The water molecules show some polarization or alignment along the field.
What I Observed
Ion motion appears random — I don’t observe any clear directional drift of cations and anions.
When I switched from NVT+efield to NVE+efield, the system temperature kept increasing with time, even though total energy should be conserved in principle.
Questions
Is the temperature rise under NVE+efield expected behavior due to Joule heating?
What’s the best practice to apply an external electric field while keeping the temperature stable without destroying the ionic drift? Should I thermostat only the solvent molecules?
Are there recommended thermostat–integration combinations (e.g., fix nve + fix langevin on water only) for such nonequilibrium MD setups?
Any advice on how to properly quantify the net ionic drift under the applied field?
Any insight or experience would be highly appreciated — especially if anyone has simulated ion transport or electrophoresis under electric fields in bulk systems using LAMMPS.
No it should not. If you have an external force this is not an NVE ensemble anymore, and you are actively injecting energy in the system, so a temperature increase is expected.
That expectation is not justifiable by the physics. The forces from the applied electric field are quite small compared to the forces from the Coulomb between the water molecules and the water and the ions, yet the time of your observation is extremely short. Even experimentally, the effect of ions diffusing preferentially to differently charged electrodes is quite slow. The magnitude of an electric field needed to observe the behavior you expect would be so large, in a real experiment there would be a flash all all water would be evaporated. You would need average over a much larger system and a very much longer simulation time to see any effect of the kind you describe and realistic choices of electric fields. Most simulation studies use field strengths one or more orders of magnitude larger than what is experimentally viable. Check the corresponding published literature. There is a lot of that. I did that kind of simulation 30 years ago as an assignment when I was an undergrad.
Why? You add forces to a system, that means adding energy. When you add friction to your hands (i.e. rub them) they get warmer, too, don’t they?
Also, you have to first validate that your system that conserve energy without fix efield to make certain, you don’t have an energy drift due to bad simulation settings.
By how much? A small drift is to be expected for simulations employing floating point math and interactions with a cutoff.
You don’t say how large the system is and how high the NaCl concentration.
This is a bit of a bad idea because for fix efield the potential is inconsistent with the force when crossing a periodic boundary. Also you may induce a drift of the center of mass and eventually may become victim of the flying icecube syndrom in case you are also applying a thermostat.
As mentioned above this is not a question you should ask here, but a topic for a search of the published literature. That is the best way to get access to this kind of information. That is why people publish their findings (and to graduate and/or stay gainfully employed, of course).
Thank you so much for your response.
System size is (35 36 76) in Angstrom unit, and the concentration is 0.5 m.
I would appreciate it if you wanted to add more.
As @akohlmey says, this is a well-established field with important guidelines to observe. I can offer some pointers but you have to read the relevant literature and understand for yourself. (If a reviewer asks you to justify your simulation choices, you will not get very far saying “According to a random poster on MatSci.org …”)
In general it is possible to perform “electrophoresis” in a 3D bulk simulation. See [1] for a fairly recent example. But you must know how to measure your results. When you say you “expect Na ions to move along with the field” it is not clear whether you just performed a visual observation (loaded up the trajectory into VMD and watched for the Na ions to all move one way) or whether you actually statistically averaged the Na ionic velocities. You can only detect “electrophoresis” through the latter numerical approach, not visually (at least not for a short simulation), since any simulation in which the electrophoretic force was so dominant would be one in which most basic near-equilibrium assumptions would no longer be useful.
To me, the most sensible thermostatting setup is to maintain a simple Nose-Hoover (or Langevin) thermostat on only the solvent (water) molecules, and to use pure NVE integration for the ions. The solvent thermostat will automatically thermostat the ions through dissipation, and reasonably well, as long as the imposed electric field is not too large. (You should not apply some kind of directional thermostat, i.e. compute temp/partial, to a system with long-range orientational ordered rigid molecules [2]: it will result in systematic – although small – errors in thermostatted temperature.)
All the best. Do read the papers carefully before proceeding.