problem with fene and purely repulsive potential

[Deepika Srivastva]

Hello,

I am new to simulating protein using hoomd blue. I a, facing two problems:

1. Fene potential: for fene I am using k=20 and r0=2.0, epsilon=0.0, sigma, and delta values are different for different pairs which I am reading from a file. The problem is my system crashes with an error stating fene bond is out of order, but as soon as I am using r0=2.38 and larger it is running fine. 

2. purely repulsive potential: I am using LJ potential with different sigma values and my epsilon is 1 for all pairs. I have two issues incorporating this potential, first I am getting negative energy values (since it is a purely repulsive case, it should be positive), and secondly, I do not wish to have the repulsive potential for all pairs, but for the ones which are not more than two beads apart and are not bonded.

I am using two bead model where I have 110 backbone beads that are the same and 110 side-chain beads that are different as per the protein sequence. 

Can somebody please help me with this?

I am attaching my sample script.

I really appreciate any help you can provide.

[antoni...@googlemail.com]

Hi,

Because you are reading information from files you have not attached, it is hard to reproduce your issue. You will find that people on the mailing list will reply more likely if you can provide a minimal working example we can execute ourselves easily to find the problem. 

1. Are you really intending to run with an epsilon of 0? That would mean that the minimum of the FENE potential is at 0, regardless of your sigma value. Please check the documentation for the functional form: https://hoomd-blue.readthedocs.io/en/latest/module-md-bond.html#hoomd.md.bond.FENEWCA  Does it say "out of order" or "out of bounds" ?

2. this can be archived with setting the exclusions in your neighbor list. See the documentation for more details: https://hoomd-blue.readthedocs.io/en/latest/module-md-nlist.html 

I recommend that you create a minimal example where you do not read values from a file and initialize a system with 2 particles, and one bond between them to pinpoint the issue. Then you can extend your example to 3 or 4 particles and bonds and set the exclusions. Once that works, you can return to your 110 bead example. 

Regards,

Antonia

[Deepika Srivastva]

Dear Antonia,

Thank you for the reply.

As per your suggestion, I tried running the simulation with a small system of 3 particles and 2 bonds (script attached). For now, I only have bonded potential (FENE). Only with FENE when I run the script my FENE energy values are too high and do not match with my calculated values. The reason why I used epsilon as 0.0 is that I do not want to have WCA along with FENE which is present inFEENEWCA (https://hoomd-blue.readthedocs.io/en/latest/module-md-bond.html#hoomd.md.bond.FENEWCA). Even I tried with epsilon 1.0, still, the energy values are the same.

Can you please help me with this?

Thanks again.

Regards,

Deepika

[Joshua Anderson]

Deepika,

If 3 particles and 2 bonds is still to complex to match the results you expect, then try an even simpler system - 2 particles, 1 bond, one particle at [0,0,0] and the second at [r,0,0]. Can you also share the output of your script and show how you compute what you expect (e.g. a Python script) based on the equations in the documentation: https://hoomd-blue.readthedocs.io/en/v3.0.1/module-md-bond.html#hoomd.md.bond.FENEWCA
------
Joshua A. Anderson, Ph.D.
Research Area Specialist, Chemical Engineering, University of Michigan

> --
> You received this message because you are subscribed to the Google Groups "hoomd-users" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to hoomd-users...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/1e753c91-22ed-4d42-bbb4-a24dc2d5a1d4n%40googlegroups.com.
> <script_hoomd.py>

[Joshua Anderson]

Deepika,

HOOMD does not implement 4*epsilon(sigma/r)^6 as a potential. To make LJ purely repulsive, most research applications use Weeks-Chandler-Anderson (WCA) potentials which is LJ with r_cut=2**(1/6)*sigma.

You don't define exactly how the special pair is "not working", so the best I can suggest is that you follow instructions in the documentation.

------
Joshua A. Anderson, Ph.D.
Research Area Specialist, Chemical Engineering, University of Michigan

> On Apr 27, 2022, at 5:40 AM, Deepika Srivastva <dsri...@gmail.com> wrote:
>
> Hi Joshua,
>
> Thank you so much for your response.
>
> I have fixed my code for FENE potential, the simulation is running fine, and the FENE energy is matching my calculated energy values.
>
> Now, I have another problem.
> Firstly, I want to implement purely repulsive LJ potential of the form (4*epsilon(sigma/r)^6), which I am unable to implement in hoomd. I tried combining Mie potential with Morse potential to obtain the desired potential LJ form, but that does not seem to work.
>
> Secondly, I have a question regarding the use of special pair in the exclusions in the neighbor list (hoomd.md.nlist). I want to have two types of special pairs which I want to use for purely repulsive LJ (the one I mentioned above) and the other one for normal LJ potential. When I am providing different names to the special pair, it is not working.
> Can you please suggest what can I do?
>
> Thanks,
> Deepika

> To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/6819e84c-c6a7-4a5d-8e19-c05cc686b4c6n%40googlegroups.com.

[Deepika Srivastva]

Hi Joshua,

Thank you again.

Here, I would like to elaborate on my problem regarding the special pair.

Suppose I have 10 bead polymer in my system. For the first 5 beads, I would like to have WCA potential and for the next 5, I would like normal expanded LJ.

I have defined two types of special pair.

snapshot.pairs.N=5

snapshot.pairs.types=''special_pair1"

snapshot.pairs.typeid=[0]*5

snapshot.pairs.group=[[0,1],[1,2], [2,3], [3,4], [4,5]] 

snapshot.pairs.N=5

snapshot.pairs.types=''special_pair2"

snapshot.pairs.typeid=[0]*5

snapshot.pairs.group=[[5,6],[6,7], [7,8], [8,9], [9,10]] 

cell=hoomd.md.nlist.Cell(buffer=0.05, exclusions=('special_pair1'))

cell1=hoomd.md.nlist.Cell(buffer=0.05, exclusions=('special_pair2'))

My question is when I am stating special_pair1 in the cell list exclusion it is throwing an error. I don't know any other way by which I can define desired pairs for different potentials.

Please suggest?

Thank you in advance.

Regrads,

Deepika

[Michael Howard]

Hi Deepika,

There are two main errors with your code snippet:

1. You are reassigning all of the pairs attributes in your second code block (overwriting the ones in the first). If you want to have different types of special pairs between different sets of particles, then you need to include all types in snapshot.pairs.types and assign each pair the right typeid. I don't know whether you are actually doing anything with this snapshot later, though.

2. The name of a special pair is not a valid exclusion type. You can exclude *all* special pairs with the 'special_pair' exclusion, but not *some*. See https://hoomd-blue.readthedocs.io/en/v3.1.0/module-md-nlist.html.

Regards,

Mike

[Joshua Anderson]

Deepika,

Have you tried using just pair potentials (no special pair) ExpandedLJ can represent WCA as well expanded LJ interactions. Give your first 5 particles type A and the next 5 type B. Then set the ExpandedLJ A-A, A-B, and B-B interaction terms as needed.

------
Joshua A. Anderson, Ph.D.
Research Area Specialist, Chemical Engineering, University of Michigan

> To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/5792ca7f-5092-492d-94c0-facb93424d4cn%40googlegroups.com.

[Deepika Srivastva]

Hello Joshua and mike,

Thank you for the reply.

For my system Joshua's suggestion is not going to work as I have 21 different types of beads in my. 

As per Mike's suggestion I can define different types of special_pair by specifying different typeid, but in exclusion it  will exclude all the defined special pairs.

Is there any other way by which I can exclude different pairs for different potentials.

I have one more question regarding Yukawa potential. I am attaching a script with 3 bead example where I have two types of beads 'A' and 'B'. What I found is, giving different epsilon values for different pair does not seems to work. The Yukawa energy obtained for my system is not correct.

Can you please help me to understand.

Thank you,

Deepika

[Joshua Anderson]

Deepika,

The number of types is not relevant. Atomistic simulations have many more than 21 types. What is relevant is whether you can express the interaction potential you want using particle types and the pair interaction parameter matrix. Based on your earlier description of the problem where you listed 10 particles, you can.

Special pair forces are designed for use with atomistic force fields, like CHARMM, where some of the pair interactions are special. Specifically, the 1-4 particles in dihedrals scale their pair interaction by 0.5. These are a relatively small fraction of the pair interactions in the entire system.

If you have pairs are not special and should follow the normal pair potential parameters, then do not include them in the special pairs list.

I cannot answer your question on the Yukawa energy because you do not provide any output from your script nor do you sufficiently describe what you expect to be correct or why. I suggest that you refer to HOOMD's documentation on the functional form of the potential and how HOOMD sums pairwise energies per particle and for the total system.

In the future, please start a new thread when you have unrelated questions. In addition to a minimal script that demonstrates the problem, describe the behavior you get (including pasted output), and the behavior you expect (including your reasoning why you expected it based on the details provided in the HOOMD documentation).

------
Joshua A. Anderson, Ph.D.
Research Area Specialist, Chemical Engineering, University of Michigan

> To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/5d51d47a-ef49-49ed-ba16-39176407420bn%40googlegroups.com.
> <hoomd_fene.py>

[Deepika Srivastva]

Joshua,

Thanks again,

I will start a new thread for new topic.

My question regarding Yukawa potential is that, in my script I have 3 beads, 2 A types and 1 B type. When I am giving epsilon value to all pairs A-A, A-B, B-B, than I have certain Yukawa energy value. Since, I do not want to have yukawa for all pair types I set epsilon for B-B and A-B as 0. In that case my Yukawa energy value is 0.0. Here, I do no t understand how Yukawa Potential works.

I am attaching my script file and output log file.

Thank you in advance.

Regards,

Deepika

[Joshua Anderson]

Deepika,

You bond your A-A particles together:
> snapshot.particles.types = ['A','B','A']
> snapshot.particles.typeid = [0,1,0]#[0] * 3
> ...
> snapshot.bonds.group = [[0, 1], [0, 2]]

Pair potentials do not evaluate pairs excluded from the neighbor list and nlist.Cell excludes bonded particles by default, so what you describe is expected and documented behavior.

https://hoomd-blue.readthedocs.io/en/v3.1.0/module-md-pair.html
https://hoomd-blue.readthedocs.io/en/v3.1.0/module-md-nlist.html#hoomd.md.nlist.Cell

To improve performance, set r_cut=0 for the pair interactions that you disable with epsilon=0.

------
Joshua A. Anderson, Ph.D.
Research Area Specialist, Chemical Engineering, University of Michigan

> To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/a546b639-0117-40b0-9557-05e4d8ca9cden%40googlegroups.com.
> <log.dat><hoomd_fene.py>

[Deepika Srivastva]

Hello Joshua,

Thank you again for your reply. 
I am still facing problems with my WCA potential. Here, I have attached my script file and my log file.
I have 2 types of beads A and B with different sigma (total: 3 polymer beads). I have applied FENE and WCA to my system. WCA energy value in the output file is showing 0.0.
Can you please suggest if I am not implementing it properly or logging it incorrectly?

Thanks,
Deepika

[Joshua Anderson]

Deepika,

I cannot suggest that you are implementing something improperly or not logging correctly.

Your HOOMD script behaves exactly as I expect it to based on the documentation. As I stated in the previous message, the neighbor list `bond` exclusion causes Pair.LJ to ignore bonded particle pairs. In your simulation, that leaves the only non-excluded pair interaction between particles 1 and 2. These particles are placed further apart that your cutoff radius, so the pair interaction energy is 0.

>>> v = numpy.array([2.84, -19.69, 32.43]) - numpy.array([-1.10, -14.39, 34.04])
>>> numpy.sqrt(numpy.dot(v,v))
6.797477473298459

I doubt that they could move within the r_cut in only 10 timesteps.

If you don't understand how something works, then simplify it down and experiment until you gain an understanding. I recommended before that you try a two particle system with particles placed at [0,0,0] and [r,0,0] as a minimal test case. In this trivial system, you can experiment with bonds/no bonds, exclusions/no exclusions, and pair interaction parameters.

------
Joshua A. Anderson, Ph.D.
Research Area Specialist, Chemical Engineering, University of Michigan

> --
> You received this message because you are subscribed to the Google Groups "hoomd-users" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to hoomd-users...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/4fe2252f-2036-415d-ac4c-7e8b01f8076dn%40googlegroups.com.
> <hoomd_fene.py><log.dat>

[Deepika Srivastva]

Joshua, thanks a lot. It was helpful.