espressomd-users
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [ESPResSo-users] dihedral angle, constraint and pressure calculation


From: Arash Azari
Subject: Re: [ESPResSo-users] dihedral angle, constraint and pressure calculation
Date: Thu, 6 Feb 2014 10:17:29 -0800 (PST)

Hi Axel,

Thank you very much. You are absolutely right. I should change my script.


Best regards,
Arash

 
Arash Azari


On Thursday, February 6, 2014 8:02 PM, Axel Arnold <address@hidden> wrote:
Hi!

As a follow-up: if you use angle bonds to construct a straight polymer chain, then using dihedrals cannot work. The angle potential tries to flatten the triangles between which the dihedral interaction is calculated (see the picture in the user's guide). Therefore, the dihedral angle is actually degenerate and pretty much takes random values, which will probably lead to large forces that tear apart your polymer.

Axel

On 02/06/2014 06:49 PM, Arash Azari wrote:
Dear Axel,

Thank you very much indeed for detailed explanation.
Sorry, I forgot to mentioned that of course I am using the bond potential (FENE) and bond angle potential together with dihedral potential. I will try to adjust the parameters more carefully and physically.
The pressure calculation and spherical constraint points are really useful, thank you.


Best regards,
Arash
 
Arash Azari


On Thursday, February 6, 2014 3:39 PM, Axel Arnold <address@hidden> wrote:
Hi!

On 02/06/2014 01:05 PM, Arash Azari wrote:
> 1- I would like to simulate a rod-like (stiff) polymer and I have
> tried using dihedral angle potential.


What other potentials do you use? Without pair bonds, dihedrals are impossible, since nothing is keeping the particle s together. Also, you need an angular potential which makes sure that the dihedrals do not degenerate.

> “Note that usage of dihedrals increases the interaction range of
> bonded interactions to 2 times the maximal bond length!”
> I cannot fully understand the reason and the meaning of the above
> statement.


It just means that the largest distance at which two particles interact is twice the bond length, which is simply
the largest possible distance from one of the inner particles of a dihedral to the tips.

> To prevent bond broken message, we can increase the
> maximum bond length a little and on the other hand the dihedral
> potential multiplies it by 2.


No interactions are modified. Just the maximal range of all short ranged bonded interactions doubles, as the message states. Also, you cannot just modify the parameters at will, they are determined by the force field you want to implement.

> What is the correct way of using
> dihedral potential for rod-like system? I mean time step, warm up,
> settings, and simulation.


There is no "correct" way, just different models. For example, if you use harmonic bonds with a constant of about 30 to form the pair connections and angular bonds with a constant around 200 and angle pi, that would already stabilize a rod fairly well at the regular timestep of 0.01. I don't understand what you need dihedrals for, but of course you can add them as well. That is just not a standard model, and I cannot tell you how you would "correctly" set it up. Just what I can tell you that it will not work without pair bonds and angles in any case, and then the dihedrals are pretty superfluous.

> I had a look at the manual and I found
> these features there, GAUSSRANDOM and GAUSSRANDOMCUT , only in the
> manual and nowhere else. Did they already implemented in the codes?

First of all, GAUSSRANDOM/CUT will not solve your problemm since that is just a wrong setup, no matter what RNG you use. Both are implemented, but might be only in the git repository, so just check that out.

>  I think one way to simulate rod-like polymer is using the virtual
> sites, but honestly, in my opinion, the user manual for virtual sites
> is poor. It will be great if anyone could add some examples (which
> work!) to this section. I have tried to use virtual sites a number of
> times, but I could not (sorry, I cannot remember the details).

No, using virtual sites is not a good idea, since you potentially create long objects, which experience high torques. In fact, you don't want fully stiff rods in MD, they should be a tiny but flexible. Just like you avoid hard spheres and rather use Weeks-Chandler-Anderson. As for some examples, you can take a look at the testsuite tests for virtual sites. Of course it would be great to have some more examples, you are welcome to contribute!

> 2- When we do not use the periodic boundary condition, why we still
> need the box length? Does the code need the box length for domain
> decomposition for parallel run?

Not just in parallel, but even on a single core, a domain decomposition linked cell scheme is used. All particles outside the given box will be in a boundary box, making interactions N^2. So to be fast, better make sure your particles are in the box.

> 3- I would like to impose spherical confinement constraint on my
> system and I start with a very large sphere, then decreasing the
> radius of the sphere gradually. I thought that the sphere should be
> surrounded by the simulation box, say the diameter of the sphere
> should be smaller than the box length. In this way, the simulation is
> too slow because of the large box. If the box length in this case is
> only used to decompose the domains for parallel run, may I start with
> a smaller box length?


Again, that will be very expensive, since then the particle interactions will be calculated in a N^2 fashion. You can shrink the box with the sphere, and maybe try to limit the number of cells by "setmd max_num_cells 10000" or so. On the other hand, I also don't see a point in setting up a box that is more than twice as large as your target box. We could even create polymer melts like that.

> 4- Another question related to the pressure
> calculation in the above system. Does the pressure calculation
> procedure work in constrained and confined systems? Especially non
> rectangular system where we cannot use the change_volume command.

The pressure is just the virial pressure, using the box volume for normalization. If you have constraints, you need to correct for the accessible volume manually by multiplying the pressure with V_box/V_acc, because Espresso cannot compute the accessible volume. However, if you have soft constraints, it is sometimes not clear what is the accessible volume, that is up to you to decide.


> 5- When we get the force on the constraint, constraint force command,
> is it the force per area or just net force exerted on the constraint,
> e.g. boundary walls.


That is the total force, since force per area would only make sense for a planar wall. For all others, you actually need surface perpendicular force per area, which currently cannot be computed in Espresso.

Regards,
Axel

--
JP Dr. Axel Arnold
ICP, Universität Stuttgart
Allmandring 3
70569 Stuttgart, Germany
Email: address@hidden
Tel: +49 711 685 67609






-- 
JP Dr. Axel Arnold
ICP, Universität Stuttgart
Allmandring 3
70569 Stuttgart, Germany
Email: address@hidden
Tel: +49 711 685 67609



reply via email to

[Prev in Thread] Current Thread [Next in Thread]