Reason of using full_energy option for hybrid pair style but no kspace

[Jibao_Lu1]

Dear lammps community,

My systems have bonds, angles, dihedrals, and pair-wise interactions. There is no Coulombic interaction. The pair_style could be either hybrid of lj/cut and yukawa, or just lj/cut, depending on the systems. I want to run grand canonical simulations using gcmc and nvt. The gcmc only have deletion and insertion of atoms.

I have two questions about the full_energy option in fix gcmc:

1. for the system with only lj/cut pair style (with bonds, angles, and dihedrals), can I turn off the full_energy option?

2. for the system with hybrid pair style of lj/cut and yukawa, I know that lammps automatically turn on the full_energy option. However, the simulation is very slow with this option. Could anyone explain to me why one has to use the full_energy option for hybrid pair_style (no long-range electrostatics)?

Actually, I wrote a modified version of fix gcmc, which calls single() function in the energy() function in fix_gcmc.cpp to compute the energy difference of the system before and after the deletion/insertion, with pair_style hybrid (but no kspace). In my implemented energy() function, the energy difference is calculated by summing over all the pair-wise interactions between the inserted/deleted atom with the rest of the system within the cutoff. Using test runs (serial lammps) in which the deletion/insertion of atoms are performed by the full_energy option and the energy difference before and after the deletion/insertion is calculated by both the full_energy method and the energy() function I implemented. The energy difference computed from my implemented energy() function is exactly the same as the one computed by full_energy. However, when bonded interactions are present in the system, the energy difference calculated by my implemented energy() becomes different from the one calculated by full_energy method: sometimes they are exactly the same, but sometimes they are different (the value could be as large as twice). I think this might be because of the bonded interactions. Could anyone give me some advices on that?

Thanks,
Jibao

[akohlmey]

i can only comment on the very last remark.

Dear lammps community,

My systems have bonds, angles, dihedrals, and pair-wise interactions. There
is no Coulombic interaction. The pair_style could be either hybrid of lj/cut
and yukawa, or just lj/cut, depending on the systems. I want to run grand
canonical simulations using gcmc and nvt. The gcmc only have deletion and
insertion of atoms.

I have two questions about the full_energy option in fix gcmc:
1) for the system with only lj/cut pair style (with bonds, angles, and
dihedrals), can I turn off the full_energy option?

2) for the system with hybrid pair style of lj/cut and yukawa, I know that
lammps automatically turn on the full_energy option. However, the simulation
is very slow with this option. Could anyone explain to me why one has to use
the full_energy option for hybrid pair_style (no long-range electrostatics)?

Actually, I wrote a modified version of fix gcmc, which calls single()
function in the energy() function in fix_gcmc.cpp to compute the energy
difference of the system before and after the deletion/insertion, with
pair_style hybrid (but no kspace). In my implemented energy() function, the
energy difference is calculated by summing over all the pair-wise
interactions between the inserted/deleted atom with the rest of the system
within the cutoff. Using test runs (serial lammps) in which the
deletion/insertion of atoms are performed by the full_energy option and the
energy difference before and after the deletion/insertion is calculated by
both the full_energy method and the energy() function I implemented. The
energy difference computed from my implemented energy() function is exactly
the same as the one computed by full_energy. However, when bonded
interactions are present in the system, the energy difference calculated by
my implemented energy() becomes different from the one calculated by
full_energy method: sometimes they are exactly the same, but sometimes they
are different (the value could be as large as twice). I think this might be
because of the bonded interactions. Could anyone give me some advices on
that?

how exactly do you determine which are the pairs to loop over?
do you use a neighbor list? if yes, how do you deal with "special" neighbors?
also, do you use a half or full neighbor list? with or without
newton_pair/bond active?
...and how does your code deal with those settings?

as you can see from these questions, there are lots of details that
can go wrong here.

axel.

[Jibao_Lu1]

Hi Axel,

Thanks for your comments.

I don’t use a neighbor list; instead, I loop over all the particles in the processor, atom->nlocal + atom->nghost, which may possibly form a pair with the inserted/deleted atom within the cutoff of the respect pair potential. (actually, the lammps I used for the test is a serial version.) The sub-pair-style that a specific pair (i,j) points to is determined by !strcmp(keywords[substyle],“str_name_substyle”), where str_name_substyle is the name of a substyle of the hybrid pair style.

Since only atoms, not molecules, are inserted/deleted, which only have pair interactions with the rest of the system, I think I don’t need worry about the special neighbors, right?

Thanks,
Jibao

[Jibao_Lu1]

Dear LAMMPS community,

I fixed the problem. It was because I forgot to comment out a definition of a struct tyle in a derived class, which has already been defined in the base class. Now the code is working.

Thanks,
Jibao

[athomps]

Jibao,

Because of the complexity of the GCMC code, we adopt a cautionary approach and turn on full_energy whenever the calculation is doing something non-standard, including hybrid. It may be that this is unnecessarily restrictive. You can always remove that restriction in init() and see how it works. Let us know. Also, there is no reason why the presence of bonded interactions in the system should require full_energy.

Aidan

[Jibao_Lu1]

Hi Aidan,

Thank you for your comments.
Actually, I have a question for the FixGCMC::pre_exchange() function: in this function there are two lines:
domain->pbc();
comm->exchange();
which are called after the insertion/deletion when full_energy option is enabled. However, they are not called after the insertion/deletion if the full_energy option is turned off. Could you please give me some comments about the reason? so I can make sure that no unknown problems exist in my implementation.

Thanks,
Jibao

I attached here the code of the pre_exchange() function, in which the two lines I mentioned are bolded.

/* ----------------------------------------------------------------------

attempt Monte Carlo translations, rotations, insertions, and deletions

done before exchange, borders, reneighbor

so that ghost atoms and neighbor lists will be correct

------------------------------------------------------------------------- */

void FixGCMC::pre_exchange()

{

// just return if should not be called on this timestep

if (next_reneighbor != update->ntimestep) return;

xlo = domain->boxlo[0];

xhi = domain->boxhi[0];

ylo = domain->boxlo[1];

yhi = domain->boxhi[1];

zlo = domain->boxlo[2];

zhi = domain->boxhi[2];

if (triclinic) {

sublo = domain->sublo_lamda;

subhi = domain->subhi_lamda;

} else {

sublo = domain->sublo;

subhi = domain->subhi;

}

if (regionflag) volume = region_volume;

else volume = domain->xprd * domain->yprd * domain->zprd;

if (triclinic) domain->x2lamda(atom->nlocal);

domain->pbc();

comm->exchange();

atom->nghost = 0;

comm->borders();

if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);

update_gas_atoms_list();

if (full_flag) {

energy_stored = energy_full();

if (mode == MOLECULE) {

for (int i = 0; i < ncycles; i++) {

int random_int_fraction =

static_cast(random_equal->uniform()*ncycles) + 1;

if (random_int_fraction <= nmcmoves) {

if (random_equal->uniform() < 0.5) attempt_molecule_translation_full();

else attempt_molecule_rotation_full();

} else {

if (random_equal->uniform() < 0.5) attempt_molecule_deletion_full();

else attempt_molecule_insertion_full();

}

}

} else {

for (int i = 0; i < ncycles; i++) {

int random_int_fraction =

static_cast(random_equal->uniform()*ncycles) + 1;

if (random_int_fraction <= nmcmoves) {

attempt_atomic_translation_full();

} else {

if (random_equal->uniform() < 0.5) attempt_atomic_deletion_full();

else attempt_atomic_insertion_full();

}

}

}

if (triclinic) domain->x2lamda(atom->nlocal);

domain->pbc();

comm->exchange();

atom->nghost = 0;

comm->borders();

if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);

} else {

if (mode == MOLECULE) {

for (int i = 0; i < ncycles; i++) {

int random_int_fraction =

static_cast(random_equal->uniform()*ncycles) + 1;

if (random_int_fraction <= nmcmoves) {

if (random_equal->uniform() < 0.5) attempt_molecule_translation();

else attempt_molecule_rotation();

} else {

if (random_equal->uniform() < 0.5) attempt_molecule_deletion();

else attempt_molecule_insertion();

}

}

} else {

for (int i = 0; i < ncycles; i++) {

int random_int_fraction =

static_cast(random_equal->uniform()*ncycles) + 1;

if (random_int_fraction <= nmcmoves) {

attempt_atomic_translation();

} else {

if (random_equal->uniform() < 0.5) attempt_atomic_deletion();

else attempt_atomic_insertion();

}

}

}

}

next_reneighbor = update->ntimestep + nevery;

}

[athomps]

Those function calls ensure that added or removed local atoms are communicated to neighboring processors. I don’t know why it is only done for full_energy.

[Jibao_Lu1]

Hi Didan,

A related question: for the system with only pair-wise interactions, if the full_energy option is turned off (which means that the energy difference before and after the insertion/deletion is calculated by looping over the particles in only processor in which the atom is inserted/deleted), how would the lammps ensure that all the interactions are accounted for? For example, there might be some particles, in the processors other than the one where the insertion/deletion occurs, locating within the cutoff of interactions with the inserted/deleted particles. In this case, how does lammps account for these interactions if it only loop over particles in the processor where insertion/deletion occurs?

Thanks,
Jibao

[athomps]

It loops over local and ghost atoms.