DPD fix_nve - Verlet - electrostatic

Dear LAMMPS users,

I’m running some simulation with the Dissipative Particle Dynamics (DPD) method (using pair_style dpd) and I have three questions about its implementation in LAMMPS.
I think all these questions can be answered with: “Write your own code to do it”, but no one in my team is a developer and knows C/C++… I only have learnt Python to write my own scripts…
I ask these questions because I found many papers in which DPD was successfully performed with LAMMPS recently, and being my questions very basic for DPD simulations I think there must be a simple solution for at least the first two problems…

Here I describe briefly my system for reference:

  • The units are lj;

  • The density is 3 beads/Rc^3;

  • I have a solute composed by 4 different beads/atom types (a gold nanoparticle coated with ligands) in a water solution;

  • The gold nanoparticle (NP) is treated as a rigid body with fix rigid/nve.

  • The ligands can’t enter the NP, and neither can escape it. This is obtained using two spherical “fix indent” centered on the NP center of mass;

  • The beads forming a ligand are kept together with “bond_style harmonic” and “angle_style harmonic”;

  • The timestep is 0.01 tau.

Here are the questions:

A) I’m using “fix nve” to perform the time integration, but the temperature is not well controlled and ultimately my system blows up after 2-3 million steps. In DPD a modified version of the Velocity Verlet algorithm is needed, but I can’t find it anywhere in the documentation! This seems strange to me, considering that LAMMPS has the DPD pair_style well implemented and optimized.

B) Accidentally I ran some of my simulations with “fix nvt” instead of nve, I know this is wrong because I’m thermostating two times my system, but those were the only simulations that were stable and finished without errors. Also, visually the trajectory was fine and the behaviour of the ligands on the nanoparticle was as expected… I don’t know how to interpret that.

C) I want to simulate some charged particles, but I can’t find a way to use a charge distribution on DPD particles instead of point charges, as is requested for DPD. Maybe some kspace_style can help on this regard?

I know these are many questions and I thank you all for spending some time reading and maybe answering to them.
I hope this thread can be helpful for others wanting to use the amazing capabilities of LAMMPS for DPD simulations.

Best regards,