Adding new potential to lammps

Hi all,

I am trying to implement in lammps a new kind of a pair interaction potential which is defined by
U® = 0 if r < L,
U® = (r - L)^2 if r>L.
It also has a hard sphere potential if r is less than a particle radius r0.

Is there any suggestion for me to do it in lammps? Any help would be appreciate.

Unless speed is of paramount importance to you, I would use
"pair_style table". It might not be as fast as writing your own code,
but "table" has gpu support, which is nice in the long run.

I realize that I did not really answer the question you asked.

It's possible you could write your own pair style which might be twice
as fast, but I suspect it is not worth the trouble. (If you write
your own pair style, then in addition to writing your own compute()
function which calculates forces and energies, you probably should
make your own coeff() and read/write_restart() functions which parsing
parameters, and restart files. If you want to do this, the best
example to look at is probably:
src/pair_lj_cut.cpp
src/pair_lj_cut.h

---- documentation for pair_style table ----

http://lammps.sandia.gov/doc/pair_table.html

Here is an example:

---- simple example ----

pair_style table linear 1001
pair_coeff 1 1 table_int.dat INT

---- The format of the "table_int.dat" file is shown here: ----

http://www.moltemplate.org/examples/membrane+protein/table_int.dat

---- hybrid example ----

In this example, a "pair_style table" is used for the interaction
between particles of type 1 with each other. All other interactions
(1-2 and 2-2) are handled with simple Lennard-Jones forces.

pair_style hybrid table linear 1001 lj/cut 15.0
pair_coeff 1 1 table table_int.dat INT
pair_coeff 1 2 lj/cut 0.1643 7.5 # <-- use 1 2 not 2 1
pair_coeff 2 2 lj/cut 0.1643 7.5

(In a hurry to reply. I might be wrong, but I think this example works...)
------------ long range, short range -------------

Either way, it sounds like some of your particles will use
long-distance cutoffs (L), and other's won't. I don't know if it it
would make your simulation faster, but you might want to consider
trying the "multi" option for the "neighbor" and "communicate"
commands.

http://lammps.sandia.gov/doc/neighbor.html
http://lammps.sandia.gov/doc/communicate.html

I hope this helps.

Andrew

Hi all,

I am trying to implement in lammps a new kind of a pair interaction
potential which is defined by
U(r) = 0 if r < L,
U(r) = (r - L)^2 if r>L.
It also has a hard sphere potential if r is less than a particle radius r0.

Is there any suggestion for me to do it in lammps? Any help would be
appreciate.

there are two major problems with your model

- i don't see a cutoff. pairwise potentials in LAMMPS are implemented
under the assumption that they are short ranged. and your potential
cannot be of infinite range, i.e. handled with a long-range solver, as
it would be divergent.

- you would need to implement a different integrator for hard sphere
interactions. with a hard sphere potential, you must not have any
overlap of particles, but that is what would happen with LAMMPS'
existing integrators. you would have to write an integrator that can
detect collisions and then adjust the velocities and positions
accordingly. in general, time integration for hard sphere is not very
efficient with a fixed time step.

axel.

For what it's worth, I actually did implement a potential like this
once. I used pair_style table to keep a small polymer trapped in a
spherical container (which could move by diffusion). When multiple
spherical containers were present, we used simple dumb repulsive
lennard-jones forces to keep them from overlapping (an ad-hoc
solution, but no more ridiculous than the other ad-hoc forces we were
using).

(Images attached) That example is located in
lammps/tools/moltemplate/examples/CG_biomolecules/protein_folding_examples/1bead+chaperone/frustrated+chaperonin/
(Warning: the files are in "moltemplate" format. This is is a
somewhat ugly and complex example.)

Cheers
Andrew

unfolded+chaperonin_t=508750tau_LR.jpg

protein2x2x2+minichaperones2x2x2_t=67500tau_LR.jpg

For what it's worth, I actually did implement a potential like this
once. I used pair_style table to keep a small polymer trapped in a
spherical container (which could move by diffusion). When multiple
spherical containers were present, we used simple dumb repulsive
lennard-jones forces to keep them from overlapping (an ad-hoc
solution, but no more ridiculous than the other ad-hoc forces we were
using).

but your model is _different_ in the two _crucial_ points:

- you have a soft repulsion not a hard sphere
- you have a cutoff which is inferred from the singularity of the wall
of your container

axel.