DPD fix_nve - Verlet - electrostatic

[Domenico_Marson]

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,
Domenico.

[sjplimp]

There aren’t enough details for any of your Qs. If you have a paper

where someone says they did a DPD simulation with LAMMPS and

you can’t figure out how they did it, I suggest you contact them

directly. Even if they modified the code, they might be willing

to share their changes with you.

Brief comments:

(a) what modifications? who says it is needed? If they

did it with LAMMPS, did they modify something?

(b) using 2 thermostats on the same particles is simply wrong.

that doesn’t mean things may not “look” right

© again, who says you need spatially distributed charge on

a DPD particle. Did someone do that with LAMMPS? If so,

they must have modified the code and implemented their

own DPD pair style to do it.

Steve

[akohlmey]

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.

i don't see how these two issues are related. if your system "blows
up", then this is almost always due to a problem with your model
and/or the simulation settings. the dpd pair style in LAMMPS does some
velocity dependent operations already during the force computation,
how do you know that this is not equivalent to the modifications you
talk about (...and what are they?).

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.

you just hid one error by making another error. happens in MD
simulations all the time. but a double error doesn't make it right,
just more wrong. that a simulation finishes without stopping is a
necessary requirement for a meaningful MD simulation, but it is not a
sufficient proof.

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?

over the years, i've heard several people voice their desire to
implement a DPD variant for charged particles. so far this has not
happened. there are papers describing such models, yet the model
implemented in LAMMPS does not include any of that kind of extension.
if you want that feature, you have multiple choices:
- find an MD code that has it implemented (i.e. not LAMMPS)
- implement it yourself
- failing that, find a competent collaborator that is willing to
implement it for you or contract a capable developer.

axel.

[Domenico_Marson]

Thank you Axel and Steve for your comments and suggestions.

Just to clarify and to answer Steve’s questions, I’m referring in particular to these papers, where the authors state they’re using LAMMPS:

• Y. Li, X. Li, Z. Li, and H. Gao, “Surface-structure-regulated penetration of nanoparticles across a cell membrane,” Nanoscale, vol. 4, no. 12, pp. 3768–8, 2012.

Here the authors say they use the modified velocity-Verlet algorithm. This is needed to have a controlled temperature, as is stated in [M. B. Liu, G. R. Liu, L. W. Zhou, and J. Z. Chang, “Dissipative Particle Dynamics (DPD): An Overview and Recent Developments,” Archives of Computational Methods in Engineering, pp. 1–28, Apr. 2015.]. A simple velocity-Verlet algorithm is insufficient for DPD simulations because the force between particles depend on both position AND velocities of the particles, so the velocities have to be estimated before the computation of the forces.

• H.-M. Ding and Y.-Q. Ma, “Design maps for cellular uptake of gene nanovectors by computer simulation,” Biomaterials, vol. 34, no. 33, pp. 8401–8407, Nov. 2013.

Here the authors took into accounts the electrostatic interactions. The soft potential in DPD can lead to the formation of artificial ion pairs and cause the divergence of the electrostatic potential. To avoid this they spread out the charges with a distribution reported in the paper.

The relation between the modified Verlet algorithm and the “blowing up” of my system was suggested because I saw a slow and constant increase in the temperature/kinetic energy, and the modified velocity-Verlet algorithm was developed to achieve a better control of the temperature, as far as I know. Also I had already double/triple checked all my simulation parameters and I couldn’t find the culprit unfortunately.

I didn’t know that the DPD pair style in LAMMPS does some velocity dependent operations during the force computation, I’m starting to think that this maybe can answer my question about the modification of the integration algorithm… can you clarify on that or point me in some part of the documentation where this is explained? As far as I know what is needed is that the velocities have to be estimated before the computation of the forces, used for the forces computations, and than corrected, at every time-step…

Unfortunately, I see from Axel answer that the treatment of charged particles is far from being implemented officially in LAMMPS, so I’ll follow his suggestions on this aspect.

If anyone has other experiences on DPD simulations with LAMMPS I’m open to every suggestion!

Regards,

Domenico.

[akohlmey]

Thank you Axel and Steve for your comments and suggestions.

Just to clarify and to answer Steve's questions, I'm referring in particular
to these papers, where the authors state they're using LAMMPS:

if these people state that they are using LAMMPS, then they have to
have modified it and never contributed their changes back to the
upstream version. you have to contact the authors for their
modifications. i have strong personal opinions on this kind of
behavior, but that has no place here.

[...]

The relation between the modified Verlet algorithm and the "blowing up" of
my system was suggested because I saw a slow and constant increase in the
temperature/kinetic energy, and the modified velocity-Verlet algorithm was
developed to achieve a better control of the temperature, as far as I know.
Also I had already double/triple checked all my simulation parameters and I
couldn't find the culprit unfortunately.

i disagree with that assessment. a constant increase in temperature,
is almost always a sign for too large a timestep for the given type
and shape(!) of interactions. it is unlikely to be cause by the DPD
implementation, as DPD interactions are deliberately soft, but rather
by other components of your model, which may not be as soft. since DPD
already includes some thermostatting, an increase of temperature is a
*very* bad sign.

i would suggest to start experimenting with a simpler system and
validating the settings for each component independently or
introducing them incrementally.

I didn't know that the DPD pair style in LAMMPS does some velocity dependent
operations during the force computation, I'm starting to think that this
maybe can answer my question about the modification of the integration
algorithm... can you clarify on that or point me in some part of the
documentation where this is explained? As far as I know what is needed is
that the velocities have to be estimated before the computation of the
forces, used for the forces computations, and than corrected, at every
time-step...

i don't think that what you are referring to is done in LAMMPS.
i also am not certain of the impact of this discrepancy. it most
certainly cannot explain the systematic increase in temperature. what
you are referring to is a far more subtle effect and will mostly
affect computation of statistical dynamical properties.

axel.

[sjplimp]

I’m CCing two people who have done DPD with LAMMPS.

They can comment on timestep issues and integration

issues. I think it is possible to do standard velocity-Verlet

with DPD, possibly with a smaller timestep. You might also

need some thermostatting to reduce energy drift. But I think

its an acceptable approximation. There may be alternative

integrations that are more accurate or allow larger timesteps.

Jim Larentzos has one for hard matter, which he can point you to. Their work

is in LAMMPS and will likely be released in the main code

in the next 6 months. Dave Heine has experiemented with

DPD in LAMMPS for various soft-matter models.

Re: point charge vs smeared charge. Again, I think both

can be acceptable models. LAMMPS can do the point

charge (on top of DPD). If you want smeared charge you

would need to add a pair style that computes the interaction

of 2 spherical Gaussians (or the like) as a function

of distance.