DPD with harmonic bonds

Dear all

I am doing DPD on a polymeric system connected by harmonic bonds.

a) Using special bonds lj 1.0 1.0 1.0 and

b) Without any special bonds at all.

From what I understood both the above scenarios should be identical as in both cases there will not be any scaling of forces or energy (because prefactor is 1.0). In one case (with special bonds), it was terminating with "Bond atoms 39 40 missing on proc 10 at step 173 " and without special bonds it was running.

Is my understanding correct? If not what am I missing here?

I am using the latest 27-Nov-2013 version. Running on 12-processors.

Attached is the script(comment out special bonds command to get part b) and their corresponding output.

Data file is generated by chain.py with num of chains = 100 and mon per chain = 64.

Thanks in advance

Vaidyanathan M S

UT Austin

Dept of Chemical Engineering

dpdequil_1 (1.17 KB)

log.dpd (5.12 KB)

withspl_bonds.dpd (3.89 KB)

Dear all

I am doing DPD on a polymeric system connected by harmonic bonds.

a) Using special bonds lj 1.0 1.0 1.0 and

b) Without any special bonds at all.

From what I understood both the above scenarios should be identical as in

no, they should not.

both cases there will not be any scaling of forces or energy (because
prefactor is 1.0). In one case (with special bonds), it was terminating with
"Bond atoms 39 40 missing on proc 10 at step 173 " and without special bonds
it was running.

Is my understanding correct? If not what am I missing here?

you are wrong. the default setting for special_bonds is lj/coul 0.0 0.0 1.0

why would you want to have something else?

axel.

Dear Axel,

Thanks for the reply.

I used special_bonds lj 1.0 1.0 1.0 because I wanted DPD interactions to be
present between every particle irrespective of whether they are bonded or
not.
I did simulations with and without special_bonds as a check.

Side note: special_bonds lj 0.0 0.0 1.0 will nullify the 1-2 and 1-3 pair
interactions, wont they, which is not I wanted.

Thanks

Vaidyanathan M S

your email is correct, which is why Axel said

you will get different behavior with 2 different
settings. If you are losing atoms it typically
means they are interacting too strongly and one
is getting blown out of the simulation box.

I.e. you have a bad model.

Steve

Dear Axel/Steve

Sorry for the confusion in my lines. Thanks for the replies. Couple of more questions

Q1: In the first mail when Axel comment “it should not”, are you referring to it should not be “equal” or it should not “matter”? I am sorry that I did not understand this. I dont see any reason why they should not be the same if i have the same input configs (apart from the differences that come from stochastisity)

Q2: Even if i override the default values of special_bonds lj/coul, it should work, shouldn’t it? So if i use 0.0 1.0 1.0, it means there will not any DPD interactions between the paired molecules or will it still stay in default?

Q3: Any ideas as to what should be the general communicate cutoff for DPD simulations (i know that this depends on the system size and parameters chosen). But a value of communicate cutoff 2.0 is blowing up where as a value of 4.0 is not. I am not sure whether it is a large value.

Thanks and Regards

Vaidyanathan M S

Dear Axel/Steve

Sorry for the confusion in my lines. Thanks for the replies. Couple of more
questions

Q1: In the first mail when Axel comment "it should not", are you referring
to it should not be "equal" or it should not "matter"? I am sorry that I did

it should not be equal.

not understand this. I dont see any reason why they should not be the same

i do see a good reason.

if i have the same input configs (apart from the differences that come from
stochastisity)

i already explained. the default settings for special_bonds are
lj/coul 0.0 0.0 1.0 which will be applied if you omit the
special_bonds command. *many* commands in LAMMPS that modify *global*
settings have defaults that will be used, if you omit that command.

Q2: Even if i override the default values of special_bonds lj/coul, it
should work, shouldn't it? So if i use 0.0 1.0 1.0, it means there will not

"it should work" is a horribly imprecise statement. nobody knows what
you are referring to exactly. only if you ask precise questions, you
can expect specific answers.

any DPD interactions between the paired molecules or will it still stay in
default?

there is a global setting for exclusions (well, one for lj and one for
coulomb) and the special_bonds command adjusts it.
thus i don't understand what you mean by your statement.

Q3: Any ideas as to what should be the general communicate cutoff for DPD
simulations (i know that this depends on the system size and parameters
chosen). But a value of communicate cutoff 2.0 is blowing up where as a
value of 4.0 is not. I am not sure whether it is a large value.

there is no need to set the communicate cutoff unless your bonds may
be stretched longer than the largest cutoff plus skin.

please re-read the documentation for the communicate command. it
explains what the cutoff does and how it is applied.

axel.

Dear Axel

Thanks. I was/am a bit confused. I shall re-read the documentation really careful before asking any other doubts.

DPD has a cutoff for just 1 (lj units) which is way too less for the bonded interactions as bonds do extend more than that.

If i have to set a communication value of above 2, it means bonds are stretching to more than 2 units. I shall check where is it going bad.

Thanks alot for the patient replies :slight_smile:

Regards

Vaidyanathan M S

Dear all

Just as a continuation, to the above one (i am trying to figure out where am i going wrong)

The output is showing ave special neighs/atom = 5.8125.

I am using special_bonds lj 0.0 1.0 1.0 for the harmonic bonds and DPD pair potential.

Does it signify how many atoms are bonded to each other on an average using the special_bonds command?

This could help me to check whether i am going absurdly wrong using the command and hence i asked.

Thanks in advance

Vaidyanathan M S

special/neighs per atom means the average count
of 1-2 + 1-3 + 1-4 neighbors that each atom has
in the molecule topology of bonds you have defined.

E.g. if you had long linear chains, each atom would
have 2 1-2 neighs, 2 1-3 neighs, and 2 1-4 neighs,

for an “ave” of 6, ignoring end effects. So that’s
probably what you have, since 5.8 includes the end effects.

Steve