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

[Marcello Sega]

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