Re: [ESPResSo-devel] [bug #46210] "analyze pressure" affecting simulatio

On Mon, Oct 19, 2015 at 9:34 AM, Markus Deserno <address@hidden> wrote:

Folks,

I strongly support Marcello’s comment.

The discussion on “pressure”, pressure tensor”, and “virial” is sometimes a bit vague and misses out on important distinctions and opportunities. Briefly, the scalar pressure is the change of the free energy under ISOTROPIC volume changes. Hence, any energy that only depends on angles (such as bending or dihedral potentials) makes no contribution to the pressure. It does generally contribute to the pressure TENSOR, but not to the scalar pressure, which is 1/3 of the trace of the pressure tensor. I share Marcello’s suspicion that the contribution of the LJ-angle potential to the scalar pressure is probably quite straightforward: keep the angles fixed and only worry about r-derivatives. Further more careful thoughts would be needed, though.

Generally, if one only wishes to know the pressure, calculating the pressure tensor and subsequently taking its trace is not just inefficient. It will be impossible for cases where we haven’t worked out how an interaction enters the tensorial stress, even though it might be much simpler to contemplate, what the contribution to the scalar virial is.

Best,

Markus

--
Dr. Markus Deserno ++1-412-268-4401 (office)
Associate Professor of Physics ++1-412-681-0648 (fax)
Carnegie Mellon University ++1-412-268-8367 (Amanda Bodnar)
5000 Forbes Avenue www.cmu.edu/biolphys/deserno/
Pittsburgh, PA 15213, USA address@hidden

On Oct 19, 2015, at 9:22 AM, Marcello Sega <address@hidden> wrote:

A small addition:

any angle-dependent potential U(\phi) gives a zero contribution to the
*scalar* virial. My guess is that for the lj-angle potential,
U(r,\phi) only the radial part (which, of course, is modulated by the
angle) will contribute to the scalar virial in the usual way. So,
although this might need to be investigated further, I guess that if
there's need for the lj-angle pressure (not pressure tensor elements!)
it could be probably implemented without too much trouble.

M.

On Mon, Oct 19, 2015 at 11:47 AM, Florian Weik <address@hidden> wrote:
Hi,
Just to be clear: I have moved ljangle out of the inner pair force function
that is used to calc the virials. There are now no side effects by analyze
pressure. Lj angle did not enter before and does not enter now into the
virials.

Regards,
Florian

On Oct 19, 2015 11:14 AM, "Peter Košovan" <address@hidden>
wrote:

Hi all,

one of the reasons for the warning in the user guide that "pressure tensor
won't work for some interactions" is that the virial formula assumes central
forces - see the original paper by Irving and Kirkwood, JCP (1950) or Evans,
J. Stat. Phys (1979). For non-central forces, there additional contribution
from torques which is missing in the implementation. A more recent
literature explains that each non-central force can be decomposed to central
forces acting on individual sites of an asymmetric particle [Thompson et al,
JCP (2009) 131, 154107; Chandra and Tadmor, J Elast (2010) 100: 63–143].
Therefore we agreed at one point to give up searching for an non-central
implementation of pressure tensor and simply put a warning that pressure
tensor (and pressure as well!) won't work for all interactions.

What it means for Nick's system: even if all the implementation issues are
fixed and analyze pressure does not interfere with the integration, still
the pressure (tensor) will not be correct for any of the anisotropic
non-bonded interactions (sec. 5.2 of the user guide). You may need to
replace the anisotropic lj-angle by several interaction sites. Or you may
try to find and implement the full anisotropic formula including the
torques. Some time ago I was unable to find such formula, so I am afraid it
might be necessary to derive it from scratch.

Therefore, let me propose a different solution: simply remove any pressure
calculations for the anisotropic potentials, and make [analyze pressure]
produce an error if called with an asymmetric potential. Unless someone
implements the anisotropic calculation including the torques.

Best regards,

peter

P.S.: one can find lots of articles simulating Gay-Berne particles and
happily calculating the pressure with the virial formula.

On Sun, Oct 18, 2015 at 10:29 AM, Ivan Cimrak <address@hidden> wrote:

Hi Marcello,

You are right about two affinity. One is linked to object in fluid and
the second to Shanchen. They do not overlap in the code, and I doubt they
are used together, however, it is a very good idea to make difference
between them and to rename Shanchen affinity.

Thanks,
Ivan

Odoslané pomocou AquaMail pre Android.
http://www.aqua-mail.com

Dňa 18. október 2015 10:10:46 používateľ Marcello Sega
<address@hidden> napísal:

Hi,

As far as I know, the virial can be calculated without problems for
multibody interactions (angle, torsional & Ewald are commonly used).
The problem comes usually with finding a definition for a local
stress, but that's not our case.

Leaving out for the time being the lj-angle contributions from the
pressure calculation seems to me reasonable, at some point I will
check which interactions are calculating properly their pressure
contribution.

Regarding affinity, there might be a clash with the affinity from
SHANCHEN: I see in

tcl/interaction_data_tcl.cpp

#ifdef AFFINITY

   REGISTER_NONBONDED("affinity", tclcommand_inter_parse_affinity);

#endif

*and*

#ifdef SHANCHEN

   REGISTER_NONBONDED("affinity",tclcommand_inter_parse_affinity);

#endif

I will soon change the shanchen affinity's name to something else,
just to avoid problems.

M.

On 17 Oct 2015 16:12, "Florian Weik" <address@hidden> wrote:

Hi,
Ok, I think a clean way to handle that is via the rotational degrees of
freedom. Non-pair forces like this can never be handled in the same way the
pair forces are handled (who are the extra forces for?).
If I remember correctly the pressure calculation via the virials works
only for pair forces, so the multi-body interactions have to be excluded
from the virial calculation.

For now I've moved the lj-angle out of the (inner) pair force loop and
made the particles const for the pair forces for good measure. The pressure
command should now have no side effects,
albeit still not returning the correct pressure (see pull request). The
fact that the pressure calculation is only correct for some interactions is
reflected in the documentation.
Of course this dug up other stuff. The affinity interaction from object
in fluid was also manipulating particles (bond creation). I am in general
not in favor of the force calculation changing the state.
Changing forces is one thing, but integrate 0 actually can do all sorts
of things, like creating and removing bonds. I think it would be much
cleaner to do a second run through the particle pairs to
do that if it is needed or just flag particles for modification. What
do you think?

Cheers,
Florian

On Sat, Oct 17, 2015 at 9:35 AM, Marcello Sega
<address@hidden> wrote:

Hi,

We're using that in Vienna! We might put some effort to implement it
properly in the future, so maybe for the time being we could leave it in as
is...

M

On 16 Oct 2015 23:33, "Florian Weik" <address@hidden> wrote:

Hi,
what Axel said. This is broken and can't be fixed easily in the
current structure as it expects only pair interactions at this points. This
probably won't fix. I am in favor of doping the interaction,
as it could actually be implemented as simple pair interaction with
rotation and torques, and since this is a cruel hack I'm not motivated to
put any effort into fixing it since apparently the original author was not
motivated to do a clean implementation. Is anybody else using that?

See also
https://github.com/espressomd/espresso/issues/439

Cheers,
Florian

On Fri, Oct 16, 2015 at 10:30 PM, Axel Arnold
<address@hidden> wrote:

Hi all,

the culprit is the LJ-angle interaction, which changes particle
forces directly, rather than returning them. That works for force
calculation, but when you calculate the pressure, your pressure is lacking
the LJ-angle interaction, and your forces are f****ed up. You can get the
simulation working by using the recalc_forces option of integrate, but your
pressure values are still be wrong.

So, just don’t use that interaction or fix the interaction to behave
well (that is, like any other two-particle interaction).

Best,
Axel

Am 15.10.2015 um 22:09 schrieb Nicholas Jin <address@hidden>:

Hi Peter,

I have created a barebones tcl script that replicates the simulation
I am trying to run. I have attached it below (demo.tcl).
As a warning, this script writes 1000 pdb files to whatever
directory it's currently in.

Commenting out the set pressure [analyze pressure] causes
qualitatively massive changes in the simulation.

Espresso 3.3.0 is compiled with the following features:
/* Generic features */
#define PARTIAL_PERIODIC
#define EXTERNAL_FORCES
#define CONSTRAINTS
#define EXCLUSIONS
#define COMFORCE
#define COMFIXED
#define MOLFORCES
#define NPT
#define DPD

/* Interaction features */
#define TABULATED
#define LENNARD_JONES
#define LENNARD_JONES_GENERIC
#define LJ_ANGLE
#define SMOOTH_STEP

#define BOND_ANGLE

Thanks for your help!
Nick

On Thu, Oct 15, 2015 at 4:22 AM, Peter Košovan
<address@hidden> wrote:

Hi Nick,

I believe that analyze pressure is used by many people and if none
has reported an exploding simulation upon the call, I suspect this issue is
related to some specific feature of your simulations.

Providing a minimal example to enable others reproduce the bug
would be very helpful. Otherwise, only a person fully familiar with your
simulations can diagnose the problem.

Best,

peter

On Wed, Oct 14, 2015 at 10:15 PM, Nicholas Jin
<address@hidden> wrote:

URL:
<http://savannah.nongnu.org/bugs/?46210>

                Summary: "analyze pressure" affecting simulation
                Project: ESPResSo
           Submitted by: njin
           Submitted on: Wed 14 Oct 2015 08:15:00 PM GMT
               Category: Simulation core
               Severity: 3 - Normal
                 Status: None
            Assigned to: None
            Open/Closed: Open
        Discussion Lock: Any
                Release: None
          Fixed Release: None

   _______________________________________________________

Details:

Hi all,

The tcl command "analyze pressure" is causing misbehavior in my
simulations.
Absent this analysis command everything runs fine, but once the
analysis runs
the simulation explodes.

This occurs in espresso-3.3.0. Let me know what other information
you might
need to diagnose this issue.

Best,
Nick (address@hidden)

   _______________________________________________________

Reply to this item at:

<http://savannah.nongnu.org/bugs/?46210>

_______________________________________________
Message sent via/by Savannah
http://savannah.nongnu.org/

--
Dr. Peter Košovan

Department of Physical and Macromolecular Chemistry
Faculty of Science, Charles University in Prague, Czech Republic

Katedra fyzikální a makromolekulární chemie
Přírodovědecká fakulta Univerzity Karlovy v Praze

www.natur.cuni.cz/chemistry/fyzchem/
Tel. +420221951290
Fax +420224919752

________________________________
Pokud je tento e-mail součástí obchodního jednání, Přírodovědecká
fakulta Univerzity Karlovy v Praze:
a) si vyhrazuje právo jednání kdykoliv ukončit a to i bez uvedení
důvodu,
b) stanovuje, že smlouva musí mít písemnou formu,
c) vylučuje přijetí nabídky s dodatkem či odchylkou,
d) stanovuje, že smlouva je uzavřena teprve výslovným dosažením
shody na všech náležitostech smlouvy.

<demo.tcl>

--------------------------------------------
Axel Arnold
Richard-Wagner-Ring 16
76437 Rastatt
Email: address@hidden
Telefon: +49 173 870 6659

--
Florian Weik

address@hidden
++49 157 85939252

--
Florian Weik

address@hidden
++49 157 85939252

--
Dr. Peter Košovan

Department of Physical and Macromolecular Chemistry
Faculty of Science, Charles University in Prague, Czech Republic

Katedra fyzikální a makromolekulární chemie
Přírodovědecká fakulta Univerzity Karlovy v Praze

www.natur.cuni.cz/chemistry/fyzchem/
Tel. +420221951290
Fax +420224919752

________________________________
Pokud je tento e-mail součástí obchodního jednání, Přírodovědecká fakulta
Univerzity Karlovy v Praze:
a) si vyhrazuje právo jednání kdykoliv ukončit a to i bez uvedení důvodu,
b) stanovuje, že smlouva musí mít písemnou formu,
c) vylučuje přijetí nabídky s dodatkem či odchylkou,
d) stanovuje, že smlouva je uzavřena teprve výslovným dosažením shody na
všech náležitostech smlouvy.

--
University of Vienna, Institute of Computational Physics

From:	Nicholas Jin
Subject:	Re: [ESPResSo-devel] [bug #46210] "analyze pressure" affecting simulation
Date:	Fri, 13 Nov 2015 16:22:56 -0500