diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp index bc8e5b4316..c678bee4a7 100644 --- a/src/MC/fix_gcmc.cpp +++ b/src/MC/fix_gcmc.cpp @@ -1288,6 +1288,7 @@ void FixGCMC::attempt_atomic_deletion_full() atom->q[i] = 0.0; } } + if (force->kspace) force->kspace->qsum_qsq(); double energy_after = energy_full(); if (random_equal->uniform() < @@ -1305,6 +1306,7 @@ void FixGCMC::attempt_atomic_deletion_full() atom->mask[i] = tmpmask; if (atom->q_flag) atom->q[i] = q_tmp; } + if (force->kspace) force->kspace->qsum_qsq(); energy_stored = energy_before; } update_gas_atoms_list(); @@ -1360,7 +1362,7 @@ void FixGCMC::attempt_atomic_insertion_full() } atom->nghost = 0; comm->borders(); - + if (force->kspace) force->kspace->qsum_qsq(); double energy_after = energy_full(); if (random_equal->uniform() < @@ -1371,6 +1373,7 @@ void FixGCMC::attempt_atomic_insertion_full() } else { atom->natoms--; if (proc_flag) atom->nlocal--; + if (force->kspace) force->kspace->qsum_qsq(); energy_stored = energy_before; } update_gas_atoms_list(); @@ -1575,10 +1578,10 @@ void FixGCMC::attempt_molecule_deletion_full() m++; atom->q[i] = 0.0; } - + toggle_intramolecular(i); } } - + if (force->kspace) force->kspace->qsum_qsq(); double energy_after = energy_full(); if (random_equal->uniform() < @@ -1600,12 +1603,14 @@ void FixGCMC::attempt_molecule_deletion_full() for (int i = 0; i < atom->nlocal; i++) { if (atom->molecule[i] == deletion_molecule) { atom->mask[i] = tmpmask[i]; + toggle_intramolecular(i); if (atom->q_flag) { atom->q[i] = q_tmp[m]; m++; } } } + if (force->kspace) force->kspace->qsum_qsq(); } update_gas_atoms_list(); } @@ -1720,6 +1725,7 @@ void FixGCMC::attempt_molecule_insertion_full() if (atom->map_style) atom->map_init(); atom->nghost = 0; comm->borders(); + if (force->kspace) force->kspace->qsum_qsq(); double energy_after = energy_full(); if (random_equal->uniform() < zz*volume*natoms_per_molecule* @@ -1744,6 +1750,7 @@ void FixGCMC::attempt_molecule_insertion_full() atom->nlocal--; } else i++; } + if (force->kspace) force->kspace->qsum_qsq(); } update_gas_atoms_list(); } @@ -1823,6 +1830,28 @@ tagint FixGCMC::pick_random_gas_molecule() return gas_molecule_id_all; } +/* ---------------------------------------------------------------------- +------------------------------------------------------------------------- */ + +void FixGCMC::toggle_intramolecular(int i) +{ + if (atom->avec->bonds_allow) + for (int m = 0; m < atom->num_bond[i]; m++) + atom->bond_type[i][m] = -atom->bond_type[i][m]; + + if (atom->avec->angles_allow) + for (int m = 0; m < atom->num_angle[i]; m++) + atom->angle_type[i][m] = -atom->angle_type[i][m]; + + if (atom->avec->dihedrals_allow) + for (int m = 0; m < atom->num_dihedral[i]; m++) + atom->dihedral_type[i][m] = -atom->dihedral_type[i][m]; + + if (atom->avec->impropers_allow) + for (int m = 0; m < atom->num_improper[i]; m++) + atom->improper_type[i][m] = -atom->improper_type[i][m]; +} + /* ---------------------------------------------------------------------- update the list of gas atoms ------------------------------------------------------------------------- */ diff --git a/src/MC/fix_gcmc.h b/src/MC/fix_gcmc.h index 46db9ee41c..7e26d0709d 100644 --- a/src/MC/fix_gcmc.h +++ b/src/MC/fix_gcmc.h @@ -50,6 +50,7 @@ class FixGCMC : public Fix { double energy_full(); int pick_random_gas_atom(); tagint pick_random_gas_molecule(); + void toggle_intramolecular(int); double molecule_energy(tagint); void update_gas_atoms_list(); double compute_vector(int); @@ -245,7 +246,7 @@ E: Illegal fix gcmc gas mass <= 0 The computed mass of the designated gas molecule or atom type was less than or equal to zero. -E: Cannot do GCMC on atoms in atom_modify first group +E: Cannot do gcmc on atoms in atom_modify first group This is a restriction due to the way atoms are organized in a list to enable the atom_modify first command.