jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Made insertion work correctly for triclinic 

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14043 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
athomps 2015-09-12 16:41:00 +00:00
parent 67771a01fb
commit ed7559f329
2 changed files with 153 additions and 65 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

217
src/MC/fix_gcmc.cpp
View File

@ -377,6 +377,9 @@ int FixGCMC::setmask()
void FixGCMC::init() void FixGCMC::init()
{ {
triclinic = domain->triclinic;
// decide whether to switch to the full_energy option // decide whether to switch to the full_energy option
if (!full_flag) { if (!full_flag) {
@ -385,7 +388,7 @@ void FixGCMC::init()
(force->pair->single_enable == 0) || (force->pair->single_enable == 0) ||
(force->pair_match("hybrid",0)) || (force->pair_match("hybrid",0)) ||
(force->pair_match("eam",0)) || (force->pair_match("eam",0)) ||
(domain->triclinic == 1) (triclinic == 1)
) { ) {
full_flag = true; full_flag = true;
if (comm->me == 0) if (comm->me == 0)
@ -614,7 +617,7 @@ void FixGCMC::pre_exchange()
yhi = domain->boxhi[1]; yhi = domain->boxhi[1];
zlo = domain->boxlo[2]; zlo = domain->boxlo[2];
zhi = domain->boxhi[2]; zhi = domain->boxhi[2];
if (domain->triclinic) { if (triclinic) {
sublo = domain->sublo_lamda; sublo = domain->sublo_lamda;
subhi = domain->subhi_lamda; subhi = domain->subhi_lamda;
} else { } else {
@ -625,12 +628,12 @@ void FixGCMC::pre_exchange()
if (regionflag) volume = region_volume; if (regionflag) volume = region_volume;
else volume = domain->xprd * domain->yprd * domain->zprd; else volume = domain->xprd * domain->yprd * domain->zprd;
if (domain->triclinic) domain->x2lamda(atom->nlocal); if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc(); domain->pbc();
comm->exchange(); comm->exchange();
atom->nghost = 0; atom->nghost = 0;
comm->borders(); comm->borders();
if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
update_gas_atoms_list(); update_gas_atoms_list();
if (full_flag) { if (full_flag) {
@ -660,12 +663,12 @@ void FixGCMC::pre_exchange()
} }
} }
} }
if (domain->triclinic) domain->x2lamda(atom->nlocal); if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc(); domain->pbc();
comm->exchange(); comm->exchange();
atom->nghost = 0; atom->nghost = 0;
comm->borders(); comm->borders();
if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
} else { } else {
@ -756,12 +759,12 @@ void FixGCMC::attempt_atomic_translation()
MPI_Allreduce(&success,&success_all,1,MPI_INT,MPI_MAX,world); MPI_Allreduce(&success,&success_all,1,MPI_INT,MPI_MAX,world);
if (success_all) { if (success_all) {
if (domain->triclinic) domain->x2lamda(atom->nlocal); if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc(); domain->pbc();
comm->exchange(); comm->exchange();
atom->nghost = 0; atom->nghost = 0;
comm->borders(); comm->borders();
if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
update_gas_atoms_list(); update_gas_atoms_list();
ntranslation_successes += 1.0; ntranslation_successes += 1.0;
} }
@ -809,8 +812,12 @@ void FixGCMC::attempt_atomic_deletion()
void FixGCMC::attempt_atomic_insertion() void FixGCMC::attempt_atomic_insertion()
{ {
double lamda[3];
ninsertion_attempts += 1.0; ninsertion_attempts += 1.0;
// pick coordinates for insertion point
double coord[3]; double coord[3];
if (regionflag) { if (regionflag) {
int region_attempt = 0; int region_attempt = 0;
@ -824,31 +831,49 @@ void FixGCMC::attempt_atomic_insertion()
region_attempt++; region_attempt++;
if (region_attempt >= max_region_attempts) return; if (region_attempt >= max_region_attempts) return;
} }
if (triclinic) domain->x2lamda(coord,lamda);
} else { } else {
coord[0] = xlo + random_equal->uniform() * (xhi-xlo); if (triclinic == 0) {
coord[1] = ylo + random_equal->uniform() * (yhi-ylo); coord[0] = xlo + random_equal->uniform() * (xhi-xlo);
coord[2] = zlo + random_equal->uniform() * (zhi-zlo); coord[1] = ylo + random_equal->uniform() * (yhi-ylo);
coord[2] = zlo + random_equal->uniform() * (zhi-zlo);
} else {
lamda[0] = random_equal->uniform();
lamda[1] = random_equal->uniform();
lamda[2] = random_equal->uniform();
// wasteful, but necessary
if (lamda[0] == 1.0) lamda[0] = 0.0;
if (lamda[1] == 1.0) lamda[1] = 0.0;
if (lamda[2] == 1.0) lamda[2] = 0.0;
domain->lamda2x(lamda,coord);
}
} }
// apply remap to protect against rounding errors
domain->remap(coord);
if (!domain->inside(coord))
error->one(FLERR,"Fix gcmc put atom outside box");
int proc_flag = 0; int proc_flag = 0;
if (coord[0] >= sublo[0] && coord[0] < subhi[0] && if (triclinic == 0) {
coord[1] >= sublo[1] && coord[1] < subhi[1] && domain->remap(coord);
coord[2] >= sublo[2] && coord[2] < subhi[2]) proc_flag = 1; if (!domain->inside(coord))
error->one(FLERR,"Fix gcmc put atom outside box");
if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
coord[1] >= sublo[1] && coord[1] < subhi[1] &&
coord[2] >= sublo[2] && coord[2] < subhi[2]) proc_flag = 1;
} else {
if (lamda[0] >= sublo[0] && lamda[0] < subhi[0] &&
lamda[1] >= sublo[1] && lamda[1] < subhi[1] &&
lamda[2] >= sublo[2] && lamda[2] < subhi[2]) proc_flag = 1;
}
int success = 0; int success = 0;
if (proc_flag) { if (proc_flag) {
int ii = -1; int ii = -1;
if (charge_flag) { if (charge_flag) {
ii = atom->nlocal + atom->nghost; ii = atom->nlocal + atom->nghost;
if (ii >= atom->nmax) atom->avec->grow(0); if (ii >= atom->nmax) atom->avec->grow(0);
atom->q[ii] = charge; atom->q[ii] = charge;
} }
double insertion_energy = energy(ii,ngcmc_type,-1,coord); double insertion_energy = energy(ii,ngcmc_type,-1,coord);
if (random_unequal->uniform() < if (random_unequal->uniform() <
zz*volume*exp(-beta*insertion_energy)/(ngas+1)) { zz*volume*exp(-beta*insertion_energy)/(ngas+1)) {
@ -974,12 +999,12 @@ void FixGCMC::attempt_molecule_translation()
x[i][2] += com_displace[2]; x[i][2] += com_displace[2];
} }
} }
if (domain->triclinic) domain->x2lamda(atom->nlocal); if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc(); domain->pbc();
comm->exchange(); comm->exchange();
atom->nghost = 0; atom->nghost = 0;
comm->borders(); comm->borders();
if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
update_gas_atoms_list(); update_gas_atoms_list();
ntranslation_successes += 1.0; ntranslation_successes += 1.0;
} }
@ -1072,12 +1097,12 @@ void FixGCMC::attempt_molecule_rotation()
n++; n++;
} }
} }
if (domain->triclinic) domain->x2lamda(atom->nlocal); if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc(); domain->pbc();
comm->exchange(); comm->exchange();
atom->nghost = 0; atom->nghost = 0;
comm->borders(); comm->borders();
if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
update_gas_atoms_list(); update_gas_atoms_list();
nrotation_successes += 1.0; nrotation_successes += 1.0;
} }
@ -1120,6 +1145,7 @@ void FixGCMC::attempt_molecule_deletion()
void FixGCMC::attempt_molecule_insertion() void FixGCMC::attempt_molecule_insertion()
{ {
double lamda[3];
ninsertion_attempts += 1.0; ninsertion_attempts += 1.0;
double com_coord[3]; double com_coord[3];
@ -1142,10 +1168,25 @@ void FixGCMC::attempt_molecule_insertion()
region_attempt++; region_attempt++;
if (region_attempt >= max_region_attempts) return; if (region_attempt >= max_region_attempts) return;
} }
if (triclinic) domain->x2lamda(com_coord,lamda);
} else { } else {
com_coord[0] = xlo + random_equal->uniform() * (xhi-xlo); if (triclinic == 0) {
com_coord[1] = ylo + random_equal->uniform() * (yhi-ylo); com_coord[0] = xlo + random_equal->uniform() * (xhi-xlo);
com_coord[2] = zlo + random_equal->uniform() * (zhi-zlo); com_coord[1] = ylo + random_equal->uniform() * (yhi-ylo);
com_coord[2] = zlo + random_equal->uniform() * (zhi-zlo);
} else {
lamda[0] = random_equal->uniform();
lamda[1] = random_equal->uniform();
lamda[2] = random_equal->uniform();
// wasteful, but necessary
if (lamda[0] == 1.0) lamda[0] = 0.0;
if (lamda[1] == 1.0) lamda[1] = 0.0;
if (lamda[2] == 1.0) lamda[2] = 0.0;
domain->lamda2x(lamda,com_coord);
}
} }
// generate point in unit cube // generate point in unit cube
@ -1168,7 +1209,6 @@ void FixGCMC::attempt_molecule_insertion()
double insertion_energy = 0.0; double insertion_energy = 0.0;
bool procflag[natoms_per_molecule]; bool procflag[natoms_per_molecule];
for (int i = 0; i < natoms_per_molecule; i++) { for (int i = 0; i < natoms_per_molecule; i++) {
MathExtra::matvec(rotmat,onemols[imol]->x[i],atom_coord[i]); MathExtra::matvec(rotmat,onemols[imol]->x[i],atom_coord[i]);
atom_coord[i][0] += com_coord[0]; atom_coord[i][0] += com_coord[0];
@ -1187,10 +1227,18 @@ void FixGCMC::attempt_molecule_insertion()
error->one(FLERR,"Fix gcmc put atom outside box"); error->one(FLERR,"Fix gcmc put atom outside box");
procflag[i] = false; procflag[i] = false;
if (xtmp[0] >= sublo[0] && xtmp[0] < subhi[0] && if (triclinic == 0) {
xtmp[1] >= sublo[1] && xtmp[1] < subhi[1] && if (xtmp[0] >= sublo[0] && xtmp[0] < subhi[0] &&
xtmp[2] >= sublo[2] && xtmp[2] < subhi[2]) { xtmp[1] >= sublo[1] && xtmp[1] < subhi[1] &&
procflag[i] = true; xtmp[2] >= sublo[2] && xtmp[2] < subhi[2]) procflag[i] = true;
} else {
domain->x2lamda(xtmp,lamda);
if (lamda[0] >= sublo[0] && lamda[0] < subhi[0] &&
lamda[1] >= sublo[1] && lamda[1] < subhi[1] &&
lamda[2] >= sublo[2] && lamda[2] < subhi[2]) procflag[i] = true;
}
if (procflag[i]) {
int ii = -1; int ii = -1;
if (onemols[imol]->qflag == 1) { if (onemols[imol]->qflag == 1) {
ii = atom->nlocal + atom->nghost; ii = atom->nlocal + atom->nghost;
@ -1416,6 +1464,7 @@ void FixGCMC::attempt_atomic_deletion_full()
void FixGCMC::attempt_atomic_insertion_full() void FixGCMC::attempt_atomic_insertion_full()
{ {
double lamda[3];
ninsertion_attempts += 1.0; ninsertion_attempts += 1.0;
double energy_before = energy_stored; double energy_before = energy_stored;
@ -1433,30 +1482,42 @@ void FixGCMC::attempt_atomic_insertion_full()
region_attempt++; region_attempt++;
if (region_attempt >= max_region_attempts) return; if (region_attempt >= max_region_attempts) return;
} }
if (triclinic) domain->x2lamda(coord,lamda);
} else { } else {
coord[0] = xlo + random_equal->uniform() * (xhi-xlo); if (triclinic == 0) {
coord[1] = ylo + random_equal->uniform() * (yhi-ylo); coord[0] = xlo + random_equal->uniform() * (xhi-xlo);
coord[2] = zlo + random_equal->uniform() * (zhi-zlo); coord[1] = ylo + random_equal->uniform() * (yhi-ylo);
coord[2] = zlo + random_equal->uniform() * (zhi-zlo);
} else {
lamda[0] = random_equal->uniform();
lamda[1] = random_equal->uniform();
lamda[2] = random_equal->uniform();
// wasteful, but necessary
if (lamda[0] == 1.0) lamda[0] = 0.0;
if (lamda[1] == 1.0) lamda[1] = 0.0;
if (lamda[2] == 1.0) lamda[2] = 0.0;
domain->lamda2x(lamda,coord);
}
} }
// apply remap to protect against rounding errors int proc_flag = 0;
if (triclinic == 0) {
domain->remap(coord); domain->remap(coord);
if (!domain->inside(coord)) if (!domain->inside(coord))
error->one(FLERR,"Fix gcmc put atom outside box"); error->one(FLERR,"Fix gcmc put atom outside box");
if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
// if triclinic, convert to lamda coords (0-1) coord[1] >= sublo[1] && coord[1] < subhi[1] &&
coord[2] >= sublo[2] && coord[2] < subhi[2]) proc_flag = 1;
double lamda[3]; } else {
if (domain->triclinic) { if (lamda[0] >= sublo[0] && lamda[0] < subhi[0] &&
domain->x2lamda(coord,coord); lamda[1] >= sublo[1] && lamda[1] < subhi[1] &&
lamda[2] >= sublo[2] && lamda[2] < subhi[2]) proc_flag = 1;
} }
int proc_flag = 0; if (proc_flag) {
if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
coord[1] >= sublo[1] && coord[1] < subhi[1] &&
coord[2] >= sublo[2] && coord[2] < subhi[2]) {
proc_flag = 1;
atom->avec->create_atom(ngcmc_type,coord); atom->avec->create_atom(ngcmc_type,coord);
int m = atom->nlocal - 1; int m = atom->nlocal - 1;
@ -1755,6 +1816,7 @@ void FixGCMC::attempt_molecule_deletion_full()
void FixGCMC::attempt_molecule_insertion_full() void FixGCMC::attempt_molecule_insertion_full()
{ {
double lamda[3];
ninsertion_attempts += 1.0; ninsertion_attempts += 1.0;
double energy_before = energy_stored; double energy_before = energy_stored;
@ -1795,10 +1857,26 @@ void FixGCMC::attempt_molecule_insertion_full()
region_attempt++; region_attempt++;
if (region_attempt >= max_region_attempts) return; if (region_attempt >= max_region_attempts) return;
} }
if (triclinic) domain->x2lamda(com_coord,lamda);
} else { } else {
com_coord[0] = xlo + random_equal->uniform() * (xhi-xlo); if (triclinic == 0) {
com_coord[1] = ylo + random_equal->uniform() * (yhi-ylo); com_coord[0] = xlo + random_equal->uniform() * (xhi-xlo);
com_coord[2] = zlo + random_equal->uniform() * (zhi-zlo); com_coord[1] = ylo + random_equal->uniform() * (yhi-ylo);
com_coord[2] = zlo + random_equal->uniform() * (zhi-zlo);
} else {
lamda[0] = random_equal->uniform();
lamda[1] = random_equal->uniform();
lamda[2] = random_equal->uniform();
// wasteful, but necessary
if (lamda[0] == 1.0) lamda[0] = 0.0;
if (lamda[1] == 1.0) lamda[1] = 0.0;
if (lamda[2] == 1.0) lamda[2] = 0.0;
domain->lamda2x(lamda,com_coord);
}
} }
// generate point in unit cube // generate point in unit cube
@ -1837,10 +1915,19 @@ void FixGCMC::attempt_molecule_insertion_full()
if (!domain->inside(xtmp)) if (!domain->inside(xtmp))
error->one(FLERR,"Fix gcmc put atom outside box"); error->one(FLERR,"Fix gcmc put atom outside box");
if (xtmp[0] >= sublo[0] && xtmp[0] < subhi[0] && int proc_flag = 0;
xtmp[1] >= sublo[1] && xtmp[1] < subhi[1] && if (triclinic == 0) {
xtmp[2] >= sublo[2] && xtmp[2] < subhi[2]) { if (xtmp[0] >= sublo[0] && xtmp[0] < subhi[0] &&
xtmp[1] >= sublo[1] && xtmp[1] < subhi[1] &&
xtmp[2] >= sublo[2] && xtmp[2] < subhi[2]) proc_flag = 1;
} else {
domain->x2lamda(xtmp,lamda);
if (lamda[0] >= sublo[0] && lamda[0] < subhi[0] &&
lamda[1] >= sublo[1] && lamda[1] < subhi[1] &&
lamda[2] >= sublo[2] && lamda[2] < subhi[2]) proc_flag = 1;
}
if (proc_flag) {
atom->avec->create_atom(onemols[imol]->type[i],xtmp); atom->avec->create_atom(onemols[imol]->type[i],xtmp);
int m = atom->nlocal - 1; int m = atom->nlocal - 1;
@ -1980,12 +2067,12 @@ double FixGCMC::molecule_energy(tagint gas_molecule_id)
double FixGCMC::energy_full() double FixGCMC::energy_full()
{ {
if (domain->triclinic) domain->x2lamda(atom->nlocal); if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc(); domain->pbc();
comm->exchange(); comm->exchange();
atom->nghost = 0; atom->nghost = 0;
comm->borders(); comm->borders();
if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
if (modify->n_pre_neighbor) modify->pre_neighbor(); if (modify->n_pre_neighbor) modify->pre_neighbor();
neighbor->build(); neighbor->build();
int eflag = 1; int eflag = 1;

1
src/MC/fix_gcmc.h
View File

@ -128,6 +128,7 @@ class FixGCMC : public Fix {
class Fix *fixshake; class Fix *fixshake;
int shakeflag; int shakeflag;
char *idshake; char *idshake;
int triclinic; // 0 = orthog box, 1 = triclinic
class Compute *c_pe; class Compute *c_pe;