Hi, I am trying to simulate a DPD polymer chain with FENE bond and bending angle potential (with the setup as following).
In Lammps, there is a repulsive LJ part in “bond_style fene”, since I am using DPD pair_style, I don’t want
the LJ part in fene, or more specifically, my goal is to have only 1st part of fene and DPD applying to 1-2 bonded atoms, and DPD applying to all other pairs, and angle potential applying to 1-2-3 atoms.
How should I achieve this? (a solution will be much appreciated.)
The reason I asked lies in the following confusions:
If “bond_style fene” is used, is the DPD pair interaction excluded for bonded 1-2 atoms? if yes,
how should I include it then? I checked the “special_bonds”, but most choices there are for lj or coul, none of them is addressing dpd.
If there is a way to include dpd for 1-2 bonded atoms, I think I can set the epsilon zero in the bond_coeff for fene to achieve what I want.
============ setup ==========
bond_style fene
bond_coeff 1 30.0 1.5 ? 1.0 # what should be put in the quesion mark, zero?
angle_style cosine
angle_coeff 1 32.0
pair_style hybrid dpd 1.0 1.0 12345 srp 0.8 1 mid exclude yes
pair_coeff 1 1 dpd 60.0 4.5 1.0
pair_coeff 1 2 none
pair_coeff 2 2 srp 100.0 0.8
Allan
See the special_bonds command. It is typically
used with FENE bonds to exclude the 1/2 pair
interactions (DPD or otherwise). When special_bonds
says “lj” it means everything but Coulombic.
Steve
Steve, thanks for replying.
Do you think the following setup is correct if I want to disable the repulsive LJ part of fene and keep the dpd interaction computed for 1-2 atoms?
bond_style fene
bond_coeff 1 30.0 1.5 0 1.0 // note the epsilon coeff is zero.
special_bonds lj 1.0 1.0 1.0
Steve, thanks for replying.
Do you think the following setup is correct if I want to disable the repulsive LJ part of fene and keep the dpd interaction computed for 1-2 atoms?
This is a superfluous question. Set up a small test and verify it yourself. It is easy to do and gives you the ultimate confirmation. Better than depending on anybody’s say so.
Axel
Hi Axel, of course I will verify it by myself. But my point is that in the very beginning I want to make sure the code (or the setup ) is absolutely right (doing as intended). I don’t want to be in a situation where after gathering simulating results and every thing looks right but still I am not sure how the code is actually doing internally.
Hi Axel, of course I will verify it by myself. But my point is that in the
very beginning I want to make sure the code (or the setup ) is absolutely
right (doing as intended). I don't want to be in a situation where after
gathering simulating results and every thing looks right but still I am not
sure how the code is actually doing internally.
i don't understand your logic here. a simple test would provide
*exactly* the answer you want and avoid the scenario that you are
afraid of.
The only test i can think of is to dump all the coordinates, velocities, and forces, see if the numbers match according to force fields used. But due to random numbers, i don’t think they match exactly. Or do you have a better idea?
<sigh>
well, but the impact of special_bonds is independent of it and the
same is true for having a FENE bond with or without the epsilon value
set to zero.
if you run with the same random seed the "random" part is not *that*
random after all. so the logical conclusion is:
1) write a data file that has only two atoms and set up an input that
on computes one frame and prints energy and forces.
2) add the FENE bond and set special_bonds to "fene"
3) change the FENE epsilon to zero
4) change special_bonds to lj/coul 1.0 1.0 1.0
repeat with a second pair of atoms without having bonds between the
two molecules.
repeat with increased complexity until you have covered all the bases.
i would have thought that this is rather obvious. at least that is
what i got taught over 20 years ago when i was first trained doing MD
simulations and checking my whether my modifications to the code were
correct.
axel.
Well, no every professor teaches the same stuff and it was 20 years ago. Thanks for your idea, i will try that.
Well, no every professor teaches the same stuff and it was 20 years ago.
20 years ago, software was *by far* not a convenient to use and making
such tests was not as straightforward and easy.
You can use the compute pe/atom command to
dump components of the energy (pair vs bond)
for individual atoms. And compute pair/local to
dump the pair wise force between pairs of
atoms (in the neighbor list) and compute bond/local
to get force between pairs of bonded atoms.
In other words those commands can be used to
debug your model and verity you are computing
what you expect.
Steve