SRP potential with fix rigid

Axel,

Actually, the problem is not solved. It only runs for a certain number of steps and then stopped due to segmentation fault.
I prepared a small system and it does not even run on one core with both new and old versions of lammps. (files are attached)
I would greatly appreciate it if you can help figure out where the bug may be.

Thanks.

Decai

box.lmps (539 KB)

dpd_box_0.in (2.5 KB)

Axel,

Actually, the problem is not solved. It only runs for a certain number of
steps and then stopped due to segmentation fault.
I prepared a small system and it does not even run on one core with both
new and old versions of lammps. (files are attached)
I would greatly appreciate it if you can help figure out where the bug may
be.

​it is not exactly a bug, but when looking at what pair_style srp does and
what fix rigid does, it becomes clear, that there isn't really a need to
use them both. when you define rigid objects, you can just add your own
repulsive particles at the right locations. it is also possible to assign
them an (almost) zero mass. due to having them included in the rigid body,
the rigid body integrator will properly update those positions, without any
need for the SRP internal position updates, which matter only for bonds
that are flexible.

​axel​

But I do have bonds that are flexible and I can also have another molecule in the system that is flexible.
The rigid body does not prevent bond crossing each other, which SRP prevents.

please explain the reason for using: special_bonds lj 1.0 0.0 0.0

also, why are there no angles and dihedrals in your (partially) rigid molecules?
when i visualize the dynamics, of your system, it looks quite bogus to me.
it looks to me, as if that this is the real problem of your system and that using pair style srp is the wrong way to address it.

apart from that, it looks as if the cause for the segfault is that due to your bogus exclusion parameters, you are confusing the code that tries to assemble the rigid objects. if i correct them, the simulation can run (but still results in bogus conformations).

looks like this is another case of GI-GO.

axel.

also, why are there no angles and dihedrals in your (partially) rigid
molecules?

This is DPD simulation. So there are no angles and dihedrals.

when i visualize the dynamics, of your system, it looks quite bogus to me.
it looks to me, as if that this is the *real* problem of your system and
that using pair style srp is the wrong way to address it.

This is not my real system. I was simplifying it to show that SRP and fix
rigid cannot be used together.

apart from that, it looks as if the cause for the segfault is that due to
your bogus exclusion parameters, you are confusing the code that tries to
assemble the rigid objects. if i correct them, the simulation can run (but
still results in bogus conformations).

What do you mean by bogus exclusion parameters?

also, why are there no angles and dihedrals in your (partially) rigid
molecules?

This is DPD simulation. So there are no angles and dihedrals.

​i don't see how one determines the other. if you want to preserve a
certain structural ​integrity of your molecules, you will have to have
angles and dihedrals. otherwise, why are you using fix rigid? that has the
same effect, but even more drastic.

when i visualize the dynamics, of your system, it looks quite bogus to me.

it looks to me, as if that this is the *real* problem of your system and
that using pair style srp is the wrong way to address it.

This is not my real system. I was simplifying it to show that SRP and fix
rigid cannot be used together.

apart from that, it looks as if the cause for the segfault is that due to
your bogus exclusion parameters, you are confusing the code that tries to
assemble the rigid objects. if i correct them, the simulation can run (but
still results in bogus conformations).

What do you mean by bogus exclusion parameters?

​please see my other question.​

looks like this is another case of GI-GO.

axel.

please explain the reason for using: special_bonds lj 1.0 0.0 0.0

I used this to prevent the bonds from collapsing.

​this makes no sense to me at all. i'd say, you are achieving the exact
opposite. you are turning off non-bonded interactions for 1-3 and 1-4
pairs, which will - due to the lack of angles and dihedrals - *promote* the
system collapsing.

​you are not convincing me.

axel.

also, why are there no angles and dihedrals in your (partially) rigid
molecules?

This is DPD simulation. So there are no angles and dihedrals.

​i don't see how one determines the other. if you want to preserve a
certain structural ​integrity of your molecules, you will have to have
angles and dihedrals. otherwise, why are you using fix rigid? that has the
same effect, but even more drastic.

I thought DPD simulations can speed up without using of angle and
dihedrals (which MD simulation has). That is one advantage of running DPD
simulation. I am using rigid for part of the system because it would be
really unrealistic to see an polyaromatic ring structure distort and
collapse.

when i visualize the dynamics, of your system, it looks quite bogus to me.

it looks to me, as if that this is the *real* problem of your system and
that using pair style srp is the wrong way to address it.

This is not my real system. I was simplifying it to show that SRP and fix
rigid cannot be used together.

apart from that, it looks as if the cause for the segfault is that due
to your bogus exclusion parameters, you are confusing the code that tries
to assemble the rigid objects. if i correct them, the simulation can run
(but still results in bogus conformations).

What do you mean by bogus exclusion parameters?

​please see my other question.

Can you explain how you correct these "exclusion parameters"? Are you
talking about special bonds?

looks like this is another case of GI-GO.

axel.

please explain the reason for using: special_bonds lj 1.0 0.0 0.0

I used this to prevent the bonds from collapsing.

​this makes no sense to me at all. i'd say, you are achieving the exact
opposite. you are turning off non-bonded interactions for 1-3 and 1-4
pairs, which will - due to the lack of angles and dihedrals - *promote* the
system collapsing.

​you are not convincing me.

Do you think turning on 1-3 and 1-4 pairs would have similar effect as SRP?
Are you suggesting turning on 1-3 and 1-4 pairs would eliminate the need of
using SRP?

also, why are there no angles and dihedrals in your (partially) rigid
molecules?

This is DPD simulation. So there are no angles and dihedrals.

​i don't see how one determines the other. if you want to preserve a
certain structural ​integrity of your molecules, you will have to have
angles and dihedrals. otherwise, why are you using fix rigid? that has the
same effect, but even more drastic.

I thought DPD simulations can speed up without using of angle and
dihedrals (which MD simulation has). That is one advantage of running DPD
simulation. I am using rigid for part of the system because it would be
really unrealistic to see an polyaromatic ring structure distort and
collapse.

​the computational cost of angles/dihedrals is negligible. the reason for
DPD being faster, is that you are representing whole chunks of your system
with soft potentials and ​that allows you to use a much larger time step.
in an extremely simplified fashion, you can rationalize DPD by taking an
atomistic simulation and then do windowed time averaging of the resulting
trajectory. you'll see how a lot of detail vanishes and molecules shrink to
simpler entities and how atoms can move "through" each other. in this
simplified rationalization, DPD essentially allows you to model this
averaged trajectory directly. a single DPD bead would then represent 10s if
not 100s of atoms or more.

​please see my other question.

Can you explain how you correct these "exclusion parameters"? Are you
talking about special bonds?

​yes. how to do that should be obvious based on my explanations. if not,
you still need to learn about MD, and it is going to be a good exercise to
figure it out on your own.

looks like this is another case of GI-GO.

axel.

please explain the reason for using: special_bonds lj 1.0 0.0 0.0

I used this to prevent the bonds from collapsing.

​this makes no sense to me at all. i'd say, you are achieving the exact
opposite. you are turning off non-bonded interactions for 1-3 and 1-4
pairs, which will - due to the lack of angles and dihedrals - *promote* the
system collapsing.

​you are not convincing me.

Do you think turning on 1-3 and 1-4 pairs would have similar effect as SRP?
Are you suggesting turning on 1-3 and 1-4 pairs would eliminate the need
of using SRP?

i​n my personal opinion, your entire model is bogus. you're throwing stuff
together ​that looks and feels inconsistent to me and then try to "fix" it
by adding something meant for different purposes on top of everything.

​i've done some testing with valgrind and other debugging tools, and have
to conclude, that one can (temporarily) trick your test system into running
with both srp and fix rigid, but it works only by chance and is not likely
to work for longer periods, simply because the additional srp ​particles
confuse the rigid fixes, as those must be forcibly removed from all groups.
when switching from fix rigid/small to plain fix rigid. the segmentation
fault is unavoidable.

i am not an expert on DPD, just know some of the basic principles. but my
impressions is, that you'd need to come up with a suitable DPD level
representation (size and time scale) for your molecules that currently seem
to be handled as atomic scale entities.

axel.