re-build the neighbor list

LAMMPS LAMMPS Mailing List Mirror
1

Dear LAMMPS developers and users,

I am currently modifying the temper.cpp to add a density-swap in addition to its own temperature-swap among replicas.

(reference: Pashek and Garcia, PRL 93, 238105 (2004))

To this end, I need to literally swap the molecules between two replicas. I realize that in the original temperature tempering process, molecules stay where they are in each replica and only the velocities are scaled.

However, for the density swap case, molecules are physically swapped (with volume being expended/shrunk to fit to the new domain), so as the neighbor list.

I wonder if I can manually rebuild the neighbor list (right after the swap) by calling a subroutine, or exchange the neighbor list array? (please help to identify since I can not find it yet…)

Any input/suggestion/comment is gratefully appreciated. Thank you for the time and help.

LC Liu

2

Dear LAMMPS developers and users,

I am currently modifying the temper.cpp to add a density-swap in addition to
its own temperature-swap among replicas.

(reference: Pashek and Garcia, PRL 93, 238105 (2004))

To this end, I need to literally swap the molecules between two replicas. I
realize that in the original temperature tempering process, molecules stay
where they are in each replica and only the velocities are scaled.

However, for the density swap case, molecules are physically swapped (with
volume being expended/shrunk to fit to the new domain), so as the neighbor
list.

I wonder if I can manually rebuild the neighbor list (right after the swap)
by calling a subroutine, or exchange the neighbor list array? (please help
to identify since I can not find it yet....)

Any input/suggestion/comment is gratefully appreciated. Thank you for the
time and help.

have a look at the files verlet.cpp, respa.cpp or min.cpp
and search for neighbor->decide()

this function is used to decide whether a neighbor list
rebuild needs to be done. the following code segment
- for the case of Neighbor::decide() returning true -
are the steps required to do a proper rebuild of the
neighbor lists for the given propagators.

HTH,
     axel.

3

Hi, Axel,

Thank you so much for the information. I tried to add “neighbor->build()” inside the temper.cpp source code, but the list seems not being updated.

Please correct me if I am wrong…neighbor->build() seems not construct the list from scratch but does updates…so, if I want a brand new list, should I invoke neighbor-init() or other subroutine?

Thank you again for the help.

LC Liu

4

Hi, Axel,

Thank you so much for the information. I tried to add "neighbor->build()"
inside the temper.cpp source code, but the list seems not being updated.

did i tell you to only call neighbor->build() ?

Please correct me if I am wrong....neighbor->build() seems not construct the
list from scratch but does updates....so, if I want a brand new list, should
I invoke neighbor-init() or other subroutine?

please re-read my previous e-mail.

axel.

5

Hi, Axel,

Thank you so much for your patience. I did not understand the logic until now. I add the following code segment with necessary variable calls, so that it passes the compiler. Yet, this addition does not seem to be functional…(my system explodes after the molecule swap). Is there a way to take a peek on the neighbor list to make sure it is updated? (i.e., which array stores the list?)

Great appreciations to the guidance.

if (n_pre_exchange) modify->pre_exchange();
if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
if (domain->box_change) {
domain->reset_box();
comm->setup();
if (neighbor->style) neighbor->setup_bins();
}
timer->stamp();
comm->exchange();
if (sortflag && ntimestep >= atom->nextsort) atom->sort();
comm->borders();
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
timer->stamp(TIME_COMM);
if (n_pre_neighbor) modify->pre_neighbor();
neighbor->build();
timer->stamp(TIME_NEIGHBOR);

6

Hi, Axel,

Thank you so much for your patience. I did not understand the logic until
now. I add the following code segment with necessary variable calls, so that
it passes the compiler. Yet, this addition does not seem to be
functional....(my system explodes after the molecule swap). Is there a way

that "explosion" is most likely due to the calculation
working *correctly* after the swap. it is not easy to
insert a molecule or even an atom into a dense system
without having instabilities due to overlaps.
have you tried that just manually on a single calculation?

to take a peek on the neighbor list to make sure it is updated? (i.e., which
array stores the list?)

there is not necessarily a single array.
neighborlists are requested through a
callback mechanism.

by default you have this:

void Pair::init_style()
{
  neighbor->request(this);
}

which causes a NeighList class instance
to be added to the pair style with the desired
settings. see pair.h:

class NeighList *list; // standard neighbor list used by most pairs

the neighborlist code then decides how to build
the individual lists. which may be just sublists of
already computed lists or a reference to an existing
list and so on. the class then has one list per local
atom that points to the beginning of the list of neighbors
of this atom.

axel.

7

In addition to what Axel told you, you might also look
at the read_dump command which sets up a new box
and reads in a completely new set of atom coords at each
iteration.

I presume by density swap for tempering that you mean
that all atoms and box settings are swapped between
the 2 replicas. If this is the case, I'm not clear why
you can't avoid the swap of all that info and swap something
else that is simpler (similar to temperature) in order to
do the replica exchange.

Steve

8

Hi, Axel and Steve,

Thank you so much for giving directions. Please forgive me being wordy, since I would like to make a rather detailed description to the problem, and hopefully someone can also benefit from your advice and experience from this thread.

I would like to reproduce the phase diagram of TIP4P-Ew proposed by Paschek et al, 2008 (DOI: 10.1002/cphc.200800539). It covers T ranging from 150 K to 360K, and density from 0.95 g/cm^3 to 1.35 g/cm^3. I used the temper command to simulate the same temperature grids as Paschek, at a couple of density settings. At high temperatures (roughly above 240K), the results look good, but at low temperatures, calculated pressures are way below the literature data. We presume this is because the states are trapped at local minimums. Hence to speed up the sampling, we need density swap in addition to temperature swap among states.

My goal is: retain the main structure of temper.cpp with an addition to density swap.

Say, I have 6 temperatures and 4 densities, so totally 24 states. Here is my idea. I place the replicas as below:

[0 1 2 3 4 5 ] [6 7 8 9 10 11] [12 13 14 15 16 17] [18 19 20 21 22 23]

The states in the bracket are of the same density, and the temperatures are repeated every 6 replicas.

(i.e., I break down the 2D matrix (T-density) into a 1D array)

T-swap: states in the bracket
density-swap: pairs such as 0<->6<->12<->18, 3<->9<->15<->21, and etc.

I was not thinking carefully yesterday…I intended to swap the physical molecule sets between the density pairs, so I need to update the neighbor list as well. But actually it could be done similarly as how you did in the T-swap. Namely, Only rescale the volumes of the states if a density swap occurs. (molecules stay where they are, only the “location of the state labels” are changed). However, this approach also requires the system domains (and other related settings) to be exchanged.

So, could you please direct me to the code piece to be invoke for the domain swap? I looked at domain.cpp…can I call reset_box() to let it do the thing after my molecule topology is expended/shrunk?

Great appreciations to all the valuable help and advice.

LC Liu

ps., I know rescaling the volume is a bit dangerous…so I will give an eye on the density difference and time step.

9

But actually it could be done similarly as how you did in the T-swap. Namely, Only rescale the volumes of the states if a
density swap occurs

yes, this sounds like a simpler strategy. You can look at the change_box
command which is invoked from the input script to change the box
volume, and (optionally) rescale all the atom positions. It has the logic
you need for calls to the Domain class to reset the various
box parameters.

When you have this worked out, please send us the upgraded temper.cpp.
as this sounds like it might be useful for other folks. I assume you will
have an option to run tempering with just T (like now), or just rho, or
both T and rho?

Steve