running soft potential : model from chain tool

LAMMPS LAMMPS Mailing List Mirror
1

Hi all,

Using the chain tool, I am trying to build a polymer for my coarse grained simulation. there are total 100 beads in a chain and total of 100 chains. mass of one bead is 68. To un overlap it using the soft potential, I am using the following code but its not running. There is one warning and one error:

warning: FENE bond too long
Error: Bad fene bond.

I gave few trial to change the pair coefficient to make more stronger bond but then I got another error: Neighbor list overflow, boost neigh_modify one. Can any body please help me to solve it?

units lj
atom_style bond
special_bonds fene

read_data 100chain.file

neighbor 0.3 bin
neigh_modify every 1 delay 1

bond_style fene
bond_coeff 1 30.0 1.5 1.0 1.0

pair_style soft 1.0
pair_coeff * * 0.0 1.0
variable prefactor1 equal ramp(0,60)
fix 1 all adapt 1 pair soft a * * v_prefactor1

velocity all create 1.0 2349852 rot yes dist gaussian

fix 2 all langevin 1.0 1.0 10.0 904297
fix 3 all nve

thermo 100
timestep 0.1

dump 1 all xyz 1000 relax.xyz

run 100000

unfix 1
unfix 2
unfix 3

Thanks
Best Regards

Md Salah Uddin

2

Hi all,

Using the chain tool, I am trying to build a polymer for my coarse grained simulation. there are total 100 beads in a chain and total of 100 chains. mass of one bead is 68. To un overlap it using the soft potential, I am using the following code but its not running. There is one warning and one error:

warning: FENE bond too long
Error: Bad fene bond.

I gave few trial to change the pair coefficient to make more stronger bond but then I got another error: Neighbor list overflow, boost neigh_modify one. Can any body please help me to solve it?

there is the micelle example, that shows how to do what you are trying to do.

axel.

3

Make sure there is no unit mismatch between your data file and your inputs script. Note that your input script is using LJ units.

Carlos

4

Thanks for your suggestions. This same code is running for another model which is built with no change of the chain definition file i.e. 320 chains and 100 monomer per chains.

Okay, I will look at the example micelle.

Thanks
Best Regards

Md Salah Uddin

5

Hi,

Following the micelle example file now, I am trying to use the soft potential push off for the built polymer model from chain tool. The only thing I have changed in the example code are bond coefficents and units. The distance between the atoms in chain definition file was 4.03. This time the program is running but at the end of the total run, it shows the error:

"Bond/angle/dihedral extent> half of periodic box length"

As the run is complete, I am running my main simulation with this better final model (discarding the error). But it is showing "Neighbor list overflow, boost neigh_modify one" . I gave several try to modify the the neighbor list command in my main simulation (initially:neigh_modify delay 0 every 1 check yes ) including the page memory increasing.

But is not solved. Would you please help me to identify the problems? (N.B: atom style is bond for soft running and molecule for the main running)

units real
neighbor 0.3 bin
neigh_modify delay 5
atom_style bond
read_data 100chain.file

special_bonds fene
pair_style soft 1.12246
pair_coeff * * 0.0 1.12246
bond_style harmonic
bond_coeff 1 6.22 4.03

velocity all create 0.45 2349852
variable prefactor equal ramp(1.0,20.0)

fix 1 all nve
fix 2 all temp/rescale 100 0.45 0.45 0.02 1.0
fix 3 all adapt 1 pair soft a * * v_prefactor

dump a all xyz 1000 softpushoff.xyz

thermo 100
run 10000

unfix 3

Thanks
Best Regards

Md Salah Uddin

6

Hi,

Following the micelle example file now, I am trying to use the soft potential push off for the built polymer model from chain tool. The only thing I have changed in the example code are bond coefficents and units. The distance between the atoms in chain definition file was 4.03. This time the program is running but at the end of the total run, it shows the error:

"Bond/angle/dihedral extent> half of periodic box length"

have you visualized the structure based on the topology data and
checked if your bonds are of proper and desired length?

As the run is complete, I am running my main simulation with this better final model (discarding the error).

it is always a bad idea to discard an error without knowing what its origin is.

But it is showing "Neighbor list overflow, boost neigh_modify one" . I gave several try to modify the the neighbor list command in my main simulation (initially:neigh_modify delay 0 every 1 check yes ) including the page memory increasing.

again, have you visualized what is going on? this kind of error may
indicate that something with overall system is badly wrong, e.g.
missing or wrong potential parameters.

But is not solved. Would you please help me to identify the problems? (N.B: atom style is bond for soft running and molecule for the main running)

short of having any ESP capabilities, it is difficult to give advice
outside of "you have to systematically narrow down sources of errors
and use a process of elimination to reassure yourself that each choice
you make is a correct one". it looks a lot like you are just trying to
randomly try out stuff and copy and paste without verifying what
certain parameters mean and whether they are meaningful in your
specific case.

units real
neighbor 0.3 bin
neigh_modify delay 5

this for example looks like bad combination of parameters for various reasons.

atom_style bond
read_data 100chain.file

special_bonds fene
pair_style soft 1.12246
pair_coeff * * 0.0 1.12246
bond_style harmonic
bond_coeff 1 6.22 4.03

these look quite strange, too.

axel.

7

Its not the issue of changing parameters without understanding. So far, there are very little options to play with the parameters. While using the chain tool, I changed the distance between monomers as 4.03 which I got from several papers for my polymer and the reduced unit (no distance less than....) as 4.08. Then rest of the work has been done by the chain tool!!!

Then for the micelle example file, the soft running portion is pretty much straight forward. I changed the unit real because in my main simulation I am using all the units in real unit (Can I not use the real unit for soft running?) and the bond harmonic potential I got the parameters from literature. you can say that the potentials are wrong but then i have nothing to say because i didn't build the potentials myself!! But based on these potentials that literature has represented further good analysis.

The portion you are saying strange is obtained from the example file which is not very hard to understand. Is the example file strange then!!!

Thanks
Best Regards

Salah

8

now, now. there is no reason to get upset about being corrected.

Its not the issue of changing parameters without understanding. So far, there are very little options to play with the parameters. While using the chain tool, I changed the distance between monomers as 4.03 which I got from several papers for my polymer and the reduced unit (no distance less than....) as 4.08. Then rest of the work has been done by the chain tool!!!

Then for the micelle example file, the soft running portion is pretty much straight forward. I changed the unit real because in my main simulation I am using all the units in real unit (Can I not use the real unit for soft running?) and the bond harmonic potential I got the parameters from literature. you can say that the potentials are wrong but then i have nothing to say because i didn't build the potentials myself!! But based on these potentials that literature has represented further good analysis.

this is exactly what i am talking about. what is your purpose here? to
take a system that is far away from equilibrium and has overlapping
particles to something that is (much) closer to what you need. does it
matter whether parameters are perfect to <1%?
not at all. you just want to get close, so that you can use the proper
parameters and the proper model without the simulation barfing all
over you. so you have to think what is going to happen. if particles
overlap, they will repel each other, quite strongly if you have a
normal interaction. that is supposed to been taken care of by
"growing" the soft potential. but you also want the chains to be at
the proper bond lengths. so what does that mean for your parameters?
- the repulsion of your atom cores should be rather weaker than stronger
- the bond strength should be rather stronger than weaker.
- lengths and sizes matter, energies can be off.

now with *this* kind of thinking you have to look at your parameters again.
- bond length is fine. anything from 4.0 to 4.1 would probably do
- bond strength. way to wimpy. typical bond strengths in classical
models are easily two orders of magnitude larger and i would probably
go even larger just for safety
- on the other hand, the LJ epsilon is often 1-2 orders of magnitude
smaller in real units than in reduced units

have you actually visualized the trajectory of this step? do your
bonds stretch too far or not?

now considering that you may (still) have fast traveling atoms and
probably longer chains than the micelle example you have to anticipate
that it may be more difficult to unoverlap the chains, so that has to
be taken into account with the neighbor list parameters. also the
differences in units, too. with this in mind, it seems like a good
idea to have at least 1/2 sigma as neighbor list skin and update the
neighbor list every step to be save. how much is your average sigma?
how does it compare to 0.3? why the need to use delay 5 since you are
just running short simulation to unoverlap you initial configuration
and performance is really not a big concern.

BUT, you *did* have warnings about very long bonds and subsequent
problems, yet you chose to ignore the warnings and not investigate if
things are behaving as they should. how can you expect people not to
tell you that you act unreasonably?

criticism is part of the business and it will get much more harsh and
unfair as you move on. don't see it as a personal attack (what would
be the benefit of that?), but as a challenge (which actually has large
benefits in the long run). you won't be motivated to do a better
effort the next time by giving you a pat on the back, even when things
don't look right.

The portion you are saying strange is obtained from the example file which is not very hard to understand. Is the example file strange then!!!

no. it works for what it does and for the system at hand. but there is
a difference between just blindly copying something and adapting
something and that difference has apparently not yet fully settled in
your mind.

axel.

9

Thanks a lot for your help. I solved that. now soft running is done. I visualized the trajectories, the bonds are stretching and there is no error. I extracted the atomic structure of last time step (100000 fs) andi tried to do my main simulation.But in my main simulation that warning (Bond extent > half of periodic box length) and error (Neighbor list overflow, boost neigh_modify one) is showing up again. Does that mean that my structure from soft running is not at equilibrium yet? I need to run the soft potential for longer time? Or I need to follow the same approach for my main code as I did for the soft running to solve the problem?

Thanks
Best Regards

Salah

10

"Bond/angle/dihedral extent> half of periodic box length"

This probably means you have an angle or dihedral interaction between
atoms which are not actually bonded together. (It could also happen
if you have a bond which has been stretched across the simulation box,
but I'm guessing you would have noticed it when you visualized the
simulation.)

I can't speak generally, but the "neighborlist overflow" problem
happens to me when somewhere in the simulation, there are too many
atoms very close together, or when I use a large distance cutoff in my
pair_style or pair_coeffs commands. (When this happens, my
simulations tend to to be extremely slow as well.) Visualization
should help you find where this is happening. (Use a visualization
program which allows you to cut or clip the viewing volume.)

---- modeling coarse-grained chains ---

Although this might not solve the current problem you are
experiencing, have you tried newer modeling tools? The "moltemplate"
program was written to make it easier to do what it sounds like you
are trying to do (build fancy coarse-grained polymers, and assembling
them into larger objects). If you use the "Angles By Type" and
"Dihedrals By Type" feature, you should probably never run into the
"Bond/angle/dihedral extent> half of periodic box length" error.

Cheers
Andrew

P.S.
There is a (somewhat complicated) example online to build a vesicle
(a spherical lipid bilayer) out of coarse-grained lipids (implemented
as short chains). That example uses PACKMOL to arrange the lipids.

http://www.moltemplate.org
(scroll to the bottom, and click on the big sphere with a hole in it)
http://www.ime.unicamp.br/~martinez/packmol/

Unfortunately, this particular example I mention is complicated
because the published coarse-grained models themselves are
complicated. (The "system.lt" file also contains a lot of extra
"pair_coeff" statements to define interactions -between- the various
different molecule types.) You can probably build a simple micelle
much more easily, especially if you only have a single type of lipid
and no other molecules.

11

Thank you Dr. Andrew.

Let me explain what I am trying to do to solve this problem. At first I used very high harmonic constant for bond and soft pair potential and ran it with very small time step 0.001. I saw that atoms are trying to move but its not getting enough space to be relaxed. Then I thought that the initial model from the chain tool was too compact. So I define larger simulation box size compared to the initial atoms position. and use (bond_coeff 1 100.00 10.00). Then I observed that the bonds are stretching and the whole model is expanding. Then on that last trajectory I used my actual harmonic bond potential which is smaller than initial one(bond_coeff 1 6.22 4.03) and now observed that the whole model is now collapsing to gain this original bond equilibrium length 4.03. I ran this same simulation several times with adding my angle table and dihedral table potential with this to make sure that the structure is equilibrium based on this harmonic bond parameters. Now my model is quite equilibrium based on my bond, angle dihedral and redefine the simulation box size based on this new atomic positions. Now when I attempted to use my pair table it showed the neighbor overflow error again which indicates the chains are too closer based on this pair table. Now I defined a larger simulation box again compared to the new atoms position and use the soft pair cut off distance as 10.00 and observed that the model is now exploding like a bomb into a spherical shape. then I used the soft pair cut off distance 2.00 to limit the explosion and allow the chains to be more apart. after this new run I added my original pair table. Now all of the potential is my original potential and the system ran. there's no error but after few time steps the error came "pair table inner cut off distance" which indicates after some time there some chains still who are coming closer. Now I am thinking to use again a larger simuation box size again and use high potential to expand it and then use my original potential to bring to the original position.

Now my question is :
1. Is this approach right?

2. From the original model of chain tool, I have already increased the box size to a certain level and still planning to expand it more. But as soon as I get an equilibrium condition based on that potential, I redefined the simulation box size based on the maximum atomic position in x, y, and z direction on the last trajectory. By this way, am I adding too much void in my model? then how to select a perfect simulation box size or how to control the density?

3. By this method after getting an equilibrium structure, is there any function in lammps to squeeze the model to make a better solid model having less void based on the potential file. and that function will make sure that this is the maximum possible solid structure on that condition?

Thanks for the information about the moltemplate. right now I have no idea how to use it. I will go through the website and try to learn. then if I face any difficulties, I will email you. Thanks again for your help.

Best Regards
Salah

12

Thank you Dr. Andrew.

this all sounds overly and needlessly complex and hints at some
principal problems with your pair potential. how have you tested that
potential?

Let me explain what I am trying to do to solve this problem. At first I used very high harmonic constant for bond and soft pair potential and ran it with very small time step 0.001. I saw that atoms are trying to move but its not getting enough space to be relaxed. Then I thought that the initial model from the chain tool was too compact. So I define larger simulation box size compared to the initial atoms position. and use (bond_coeff 1 100.00 10.00). Then I observed that the bonds are stretching and the whole model is expanding. Then on that last trajectory I used my actual harmonic bond potential which is smaller than initial one(bond_coeff 1 6.22 4.03) and now observed that the whole model is now collapsing to gain this original bond equilibrium length 4.03. I ran this same simulation several times with adding my angle table and dihedral table potential with this to make sure that the structure is equilibrium based on this harmonic bond parameters. Now my model is quite equilibrium based on my bond, angle dihedral and redefine the simulation box size based on this new atomic positions. Now when I attempted to use my pair table it showed the neighbor overflow error again which indicates the chains are too closer based on this pair table. Now I defined a larger simulation box again compared to the new atoms position and use the soft pair cut off distance as 10.00 and observed that the model is now exploding like a bomb into a spherical shape. then I used the soft pair cut off distance 2.00 to limit the explosion and allow the chains to be more apart. after this new run I added my original pair table. Now all of the potential is my original potential and the system ran. there's no error but after few time steps the error came "pair table inner cut off distance" which indicates after some time there some chains still who are coming closer. Now I am thinking to use again a larger simuation box size again and use high potential to expand it and then use my original potential to bring to the original position.

Now my question is :
1. Is this approach right?

no.

2. From the original model of chain tool, I have already increased the box size to a certain level and still planning to expand it more. But as soon as I get an equilibrium condition based on that potential, I redefined the simulation box size based on the maximum atomic position in x, y, and z direction on the last trajectory. By this way, am I adding too much void in my model? then how to select a perfect simulation box size or how to control the density?

you have to approximately know what density your system is supposed to
have. take that to compute your box volume and make it a little bit
larger to leave some room for equilibration, say 10-20%. then
unoverlap your atoms. you can easily determine from the shape of the
non-bonded potential, what parameters you need to use.

3. By this method after getting an equilibrium structure, is there any function in lammps to squeeze the model to make a better solid model having less void based on the potential file. and that function will make sure that this is the maximum possible solid structure on that condition?

you can use fix deform or fix npt with a suitably chose pressure.
there are advantages and distadvantages for both. i don't understand
what you mean by "maximum possible solid structure". in simulations,
it is always about statistics not absolute truths.

Thanks for the information about the moltemplate. right now I have no idea how to use it. I will go through the website and try to learn. then if I face any difficulties, I will email you. Thanks again for your help.

before trying anything as complex as you are trying to do, practice on
something simpler. build your skills slowly. things almost never just
magically work (and if they do, there is usually something wrong).

axel.

13

Would you please explain more about testing my potential? I don't know how to test the potential. I got these potentials for my material from a literature.

Thanks
Best Regards

Md Salah Uddin

14

Would you please explain more about testing my potential? I don't know how to test the potential. I got these potentials for my material from a literature.

sorry, but i am not your adviser and don't have the time to teach you
computer simulation basics.

(and this is the - by far - most polite answer i can come up with, so
be very grateful),
    axel.