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

[Ulf Schiller]

Hi,

I agree with Marcello and Markus as well. I'd like to advocate
transparency though with respect to potential differences between the
pressure calculation and the pressure tensor, i.e., add clear
notes/warnings that the trace of the pressure tensor may be different
from the scalar pressure and why.

Cheers,
Ulf

On 19/10/15 14:34, Markus Deserno 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 Avenuewww.cmu.edu/biolphys/deserno/
> <http://www.cmu.edu/biolphys/deserno/>
> Pittsburgh, PA 15213, address@hidden
> 
> On Oct 19, 2015, at 9:22 AM, Marcello Sega <address@hidden
> <mailto: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
>> <mailto: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 <mailto: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
>>>> <mailto: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/
>>>> <http://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
>>
> 


-- 
Dr Ulf D Schiller
Centre for Computational Science
University College London
20 Gordon Street
London WC1H 0AJ
United Kingdom

Phone: +44 (0)20 7679 5300