From 56e1032c1d0fcb7d47e8c46a5aed9e4feef5af7f Mon Sep 17 00:00:00 2001 From: Jared Date: Thu, 27 Jun 2019 17:50:45 +0200 Subject: [PATCH 1/7] Update gcmc to have a max and min --- src/MC/fix_gcmc.cpp | 24 ++++++++++++++++++++---- src/MC/fix_gcmc.h | 2 ++ 2 files changed, 22 insertions(+), 4 deletions(-) diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp index 7ab0879335..5cf655b1de 100644 --- a/src/MC/fix_gcmc.cpp +++ b/src/MC/fix_gcmc.cpp @@ -391,6 +391,14 @@ void FixGCMC::options(int narg, char **arg) overlap_cutoffsq = rtmp*rtmp; overlap_flag = 1; iarg += 2; + } else if (strcmp(arg[iarg],"min") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix gcmc command"); + min_ngas = force->numeric(FLERR,arg[iarg+1]); + iarg += 2; + } else if (strcmp(arg[iarg],"max") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix gcmc command"); + max_ngas = force->numeric(FLERR,arg[iarg+1]); + iarg += 2; } else error->all(FLERR,"Illegal fix gcmc command"); } } @@ -897,7 +905,7 @@ void FixGCMC::attempt_atomic_deletion() { ndeletion_attempts += 1.0; - if (ngas == 0) return; + if (ngas == 0 || ngas == min_ngas) return; int i = pick_random_gas_atom(); @@ -938,6 +946,8 @@ void FixGCMC::attempt_atomic_insertion() ninsertion_attempts += 1.0; + if (ngas == max_ngas) return; + // pick coordinates for insertion point double coord[3]; @@ -1252,7 +1262,7 @@ void FixGCMC::attempt_molecule_deletion() { ndeletion_attempts += 1.0; - if (ngas == 0) return; + if (ngas == 0 || ngas == min_ngas) return; // work-around to avoid n=0 problem with fix rigid/nvt/small @@ -1291,6 +1301,8 @@ void FixGCMC::attempt_molecule_insertion() double lamda[3]; ninsertion_attempts += 1.0; + if (ngas == max_ngas) return; + double com_coord[3]; if (regionflag) { int region_attempt = 0; @@ -1574,7 +1586,7 @@ void FixGCMC::attempt_atomic_deletion_full() ndeletion_attempts += 1.0; - if (ngas == 0) return; + if (ngas == 0 || ngas == min_ngas) return; double energy_before = energy_stored; @@ -1623,6 +1635,8 @@ void FixGCMC::attempt_atomic_insertion_full() double lamda[3]; ninsertion_attempts += 1.0; + if (ngas == max_ngas) return; + double energy_before = energy_stored; double coord[3]; @@ -1918,7 +1932,7 @@ void FixGCMC::attempt_molecule_deletion_full() { ndeletion_attempts += 1.0; - if (ngas == 0) return; + if (ngas == 0 || ngas == min_ngas) return; // work-around to avoid n=0 problem with fix rigid/nvt/small @@ -2001,6 +2015,8 @@ void FixGCMC::attempt_molecule_insertion_full() double lamda[3]; ninsertion_attempts += 1.0; + if (ngas == max_ngas) return; + double energy_before = energy_stored; tagint maxmol = 0; diff --git a/src/MC/fix_gcmc.h b/src/MC/fix_gcmc.h index 5d0b7aab8f..e19a42ef73 100644 --- a/src/MC/fix_gcmc.h +++ b/src/MC/fix_gcmc.h @@ -120,6 +120,8 @@ class FixGCMC : public Fix { imageint imagezero; double overlap_cutoffsq; // square distance cutoff for overlap int overlap_flag; + int max_ngas; + int min_ngas; double energy_intra; From eb9f5cf1286ece343803217c573a72961d5d1cc3 Mon Sep 17 00:00:00 2001 From: Jared Date: Thu, 27 Jun 2019 17:50:45 +0200 Subject: [PATCH 2/7] Update gcmc to have a max and min --- src/MC/fix_gcmc.cpp | 26 ++++++++++++++++++++++---- src/MC/fix_gcmc.h | 2 ++ 2 files changed, 24 insertions(+), 4 deletions(-) diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp index 7ab0879335..51f781e4d9 100644 --- a/src/MC/fix_gcmc.cpp +++ b/src/MC/fix_gcmc.cpp @@ -272,6 +272,8 @@ void FixGCMC::options(int narg, char **arg) tfac_insert = 1.0; overlap_cutoffsq = 0.0; overlap_flag = 0; + min_ngas = -1; + max_ngas = -1; int iarg = 0; while (iarg < narg) { @@ -391,6 +393,14 @@ void FixGCMC::options(int narg, char **arg) overlap_cutoffsq = rtmp*rtmp; overlap_flag = 1; iarg += 2; + } else if (strcmp(arg[iarg],"min") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix gcmc command"); + min_ngas = force->numeric(FLERR,arg[iarg+1]); + iarg += 2; + } else if (strcmp(arg[iarg],"max") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix gcmc command"); + max_ngas = force->numeric(FLERR,arg[iarg+1]); + iarg += 2; } else error->all(FLERR,"Illegal fix gcmc command"); } } @@ -897,7 +907,7 @@ void FixGCMC::attempt_atomic_deletion() { ndeletion_attempts += 1.0; - if (ngas == 0) return; + if (ngas == 0 || ngas == min_ngas) return; int i = pick_random_gas_atom(); @@ -938,6 +948,8 @@ void FixGCMC::attempt_atomic_insertion() ninsertion_attempts += 1.0; + if (ngas == max_ngas) return; + // pick coordinates for insertion point double coord[3]; @@ -1252,7 +1264,7 @@ void FixGCMC::attempt_molecule_deletion() { ndeletion_attempts += 1.0; - if (ngas == 0) return; + if (ngas == 0 || ngas == min_ngas) return; // work-around to avoid n=0 problem with fix rigid/nvt/small @@ -1291,6 +1303,8 @@ void FixGCMC::attempt_molecule_insertion() double lamda[3]; ninsertion_attempts += 1.0; + if (ngas == max_ngas) return; + double com_coord[3]; if (regionflag) { int region_attempt = 0; @@ -1574,7 +1588,7 @@ void FixGCMC::attempt_atomic_deletion_full() ndeletion_attempts += 1.0; - if (ngas == 0) return; + if (ngas == 0 || ngas == min_ngas) return; double energy_before = energy_stored; @@ -1623,6 +1637,8 @@ void FixGCMC::attempt_atomic_insertion_full() double lamda[3]; ninsertion_attempts += 1.0; + if (ngas == max_ngas) return; + double energy_before = energy_stored; double coord[3]; @@ -1918,7 +1934,7 @@ void FixGCMC::attempt_molecule_deletion_full() { ndeletion_attempts += 1.0; - if (ngas == 0) return; + if (ngas == 0 || ngas == min_ngas) return; // work-around to avoid n=0 problem with fix rigid/nvt/small @@ -2001,6 +2017,8 @@ void FixGCMC::attempt_molecule_insertion_full() double lamda[3]; ninsertion_attempts += 1.0; + if (ngas == max_ngas) return; + double energy_before = energy_stored; tagint maxmol = 0; diff --git a/src/MC/fix_gcmc.h b/src/MC/fix_gcmc.h index 5d0b7aab8f..e19a42ef73 100644 --- a/src/MC/fix_gcmc.h +++ b/src/MC/fix_gcmc.h @@ -120,6 +120,8 @@ class FixGCMC : public Fix { imageint imagezero; double overlap_cutoffsq; // square distance cutoff for overlap int overlap_flag; + int max_ngas; + int min_ngas; double energy_intra; From e5dd154366dbb7649640559eef8193e3cef92601 Mon Sep 17 00:00:00 2001 From: Jared Wood Date: Mon, 14 Oct 2019 11:13:58 +1100 Subject: [PATCH 3/7] Make max/min prevent moves already outside the bounds Previously allowed free movement outside the bounds until they were reached. Now forces movement towards the bounds --- src/MC/fix_gcmc.cpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp index 818ed01fba..516fe2521d 100644 --- a/src/MC/fix_gcmc.cpp +++ b/src/MC/fix_gcmc.cpp @@ -269,7 +269,7 @@ void FixGCMC::options(int narg, char **arg) overlap_cutoffsq = 0.0; overlap_flag = 0; min_ngas = -1; - max_ngas = -1; + max_ngas = INT_MAX; int iarg = 0; while (iarg < narg) { @@ -903,7 +903,7 @@ void FixGCMC::attempt_atomic_deletion() { ndeletion_attempts += 1.0; - if (ngas == 0 || ngas == min_ngas) return; + if (ngas == 0 || ngas <= min_ngas) return; int i = pick_random_gas_atom(); @@ -944,7 +944,7 @@ void FixGCMC::attempt_atomic_insertion() ninsertion_attempts += 1.0; - if (ngas == max_ngas) return; + if (ngas >= max_ngas) return; // pick coordinates for insertion point @@ -1260,7 +1260,7 @@ void FixGCMC::attempt_molecule_deletion() { ndeletion_attempts += 1.0; - if (ngas == 0 || ngas == min_ngas) return; + if (ngas == 0 || ngas <= min_ngas) return; // work-around to avoid n=0 problem with fix rigid/nvt/small @@ -1299,7 +1299,7 @@ void FixGCMC::attempt_molecule_insertion() double lamda[3]; ninsertion_attempts += 1.0; - if (ngas == max_ngas) return; + if (ngas >= max_ngas) return; double com_coord[3]; if (regionflag) { @@ -1584,7 +1584,7 @@ void FixGCMC::attempt_atomic_deletion_full() ndeletion_attempts += 1.0; - if (ngas == 0 || ngas == min_ngas) return; + if (ngas == 0 || ngas <= min_ngas) return; double energy_before = energy_stored; @@ -1633,7 +1633,7 @@ void FixGCMC::attempt_atomic_insertion_full() double lamda[3]; ninsertion_attempts += 1.0; - if (ngas == max_ngas) return; + if (ngas >= max_ngas) return; double energy_before = energy_stored; @@ -1930,7 +1930,7 @@ void FixGCMC::attempt_molecule_deletion_full() { ndeletion_attempts += 1.0; - if (ngas == 0 || ngas == min_ngas) return; + if (ngas == 0 || ngas <= min_ngas) return; // work-around to avoid n=0 problem with fix rigid/nvt/small @@ -2013,7 +2013,7 @@ void FixGCMC::attempt_molecule_insertion_full() double lamda[3]; ninsertion_attempts += 1.0; - if (ngas == max_ngas) return; + if (ngas >= max_ngas) return; double energy_before = energy_stored; From a2eec80f259d615567538f84b86d916eae5bacfc Mon Sep 17 00:00:00 2001 From: Jared Wood Date: Tue, 5 Nov 2019 16:32:35 +1100 Subject: [PATCH 4/7] add max and min to documentation --- doc/src/fix_gcmc.txt | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/doc/src/fix_gcmc.txt b/doc/src/fix_gcmc.txt index 3c0f2c2f17..f28cb5771f 100644 --- a/doc/src/fix_gcmc.txt +++ b/doc/src/fix_gcmc.txt @@ -48,7 +48,9 @@ keyword = {mol}, {region}, {maxangle}, {pressure}, {fugacity_coeff}, {full_energ group-ID = group-ID for inserted atoms (string) {intra_energy} value = intramolecular energy (energy units) {tfac_insert} value = scale up/down temperature of inserted atoms (unitless) - {overlap_cutoff} value = maximum pair distance for overlap rejection (distance units) :pre + {overlap_cutoff} value = maximum pair distance for overlap rejection (distance units) + {max} value = Maximum number of molecules allowed in the system + {min} value = Minimum number of molecules allowed in the system :pre :ule [Examples:] @@ -364,6 +366,12 @@ assigning an infinite positive energy to all new configurations that place any pair of atoms closer than the specified overlap cutoff distance. +The {max} and {min} keywords allow for the restriction of the number +of atoms in the simulation. They automatically reject all insertion +or deletion moves that would take the system beyond the set boundaries. +Should the system already be beyond the boundary, only moves that bring +the system closer to the bounds may be accepted. + The {group} keyword adds all inserted atoms to the "group"_group.html of the group-ID value. The {grouptype} keyword adds all inserted atoms of the specified type to the From d37ee59296576c89994350a61c9b81be658d8ae9 Mon Sep 17 00:00:00 2001 From: Jared Wood Date: Thu, 7 Nov 2019 09:24:01 +1100 Subject: [PATCH 5/7] Add example of fix gcmc max behaviour --- examples/gcmc/in.gcmc.lj.max | 70 ++++++++++++++++++++++++++++++++++++ 1 file changed, 70 insertions(+) create mode 100644 examples/gcmc/in.gcmc.lj.max diff --git a/examples/gcmc/in.gcmc.lj.max b/examples/gcmc/in.gcmc.lj.max new file mode 100644 index 0000000000..519ac162fa --- /dev/null +++ b/examples/gcmc/in.gcmc.lj.max @@ -0,0 +1,70 @@ +# GCMC for LJ simple fluid, no dynamics +# T = 2.0 +# rho ~ 0.5 +# p ~ 1.5 +# mu_ex ~ 0.0 +# comparable to Frenkel and Smit GCMC Case Study, Figure 5.8 + +# variables modifiable using -var command line switch + +variable mu index -1.25 +variable temp index 2.0 +variable disp index 1.0 +variable lbox index 5.0 + +# global model settings + +units lj +atom_style atomic +pair_style lj/cut 3.0 +pair_modify tail no # turn of to avoid triggering full_energy + +# box + +region box block 0 ${lbox} 0 ${lbox} 0 ${lbox} +create_box 1 box + +# lj parameters + +pair_coeff * * 1.0 1.0 +mass * 1.0 + +# we recommend setting up a dedicated group for gcmc + +group gcmcgroup type 1 + +# gcmc + +fix mygcmc gcmcgroup gcmc 1 100 100 1 29494 ${temp} ${mu} ${disp} max 50 + +# atom count + +variable type1 atom "type==1" +group type1 dynamic gcmcgroup var type1 +variable n1 equal count(type1) + +# averaging + +variable rho equal density +variable p equal press +variable nugget equal 1.0e-8 +variable lambda equal 1.0 +variable muex equal ${mu}-${temp}*ln(density*${lambda}+${nugget}) +fix ave all ave/time 10 100 1000 v_rho v_p v_muex v_n1 ave one file rho_vs_p.dat +variable rhoav equal f_ave[1] +variable pav equal f_ave[2] +variable muexav equal f_ave[3] +variable n1av equal f_ave[4] + +# output + +variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+${nugget}) +variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+${nugget}) +variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+${nugget}) +compute_modify thermo_temp dynamic yes +thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc v_rhoav v_pav v_muexav v_n1av +thermo 1000 + +# run + +run 10000 From ce6893e7175f838efa22f8ed77a21880855ac444 Mon Sep 17 00:00:00 2001 From: Jared Wood Date: Fri, 8 Nov 2019 12:59:39 +1100 Subject: [PATCH 6/7] Add max/min changes to documentation again --- doc/src/fix_gcmc.rst | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/doc/src/fix_gcmc.rst b/doc/src/fix_gcmc.rst index 0bb59d3afc..9f64be8699 100644 --- a/doc/src/fix_gcmc.rst +++ b/doc/src/fix_gcmc.rst @@ -51,6 +51,8 @@ Syntax *intra_energy* value = intramolecular energy (energy units) *tfac_insert* value = scale up/down temperature of inserted atoms (unitless) *overlap_cutoff* value = maximum pair distance for overlap rejection (distance units) + *max* value = Maximum number of molecules allowed in the system + *min* value = Minimum number of molecules allowed in the system @@ -385,6 +387,12 @@ assigning an infinite positive energy to all new configurations that place any pair of atoms closer than the specified overlap cutoff distance. +The *max* and *min* keywords allow for the restriction of the number +of atoms in the simulation. They automatically reject all insertion +or deletion moves that would take the system beyond the set boundaries. +Should the system already be beyond the boundary, only moves that bring +the system closer to the bounds may be accepted. + The *group* keyword adds all inserted atoms to the :doc:`group ` of the group-ID value. The *grouptype* keyword adds all inserted atoms of the specified type to the From 6623169d97b731e4d672152b10afd3e2b3ef1fca Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 8 Nov 2019 15:00:26 -0500 Subject: [PATCH 7/7] Delete example file on request of @athomps --- examples/gcmc/in.gcmc.lj.max | 70 ------------------------------------ 1 file changed, 70 deletions(-) delete mode 100644 examples/gcmc/in.gcmc.lj.max diff --git a/examples/gcmc/in.gcmc.lj.max b/examples/gcmc/in.gcmc.lj.max deleted file mode 100644 index 519ac162fa..0000000000 --- a/examples/gcmc/in.gcmc.lj.max +++ /dev/null @@ -1,70 +0,0 @@ -# GCMC for LJ simple fluid, no dynamics -# T = 2.0 -# rho ~ 0.5 -# p ~ 1.5 -# mu_ex ~ 0.0 -# comparable to Frenkel and Smit GCMC Case Study, Figure 5.8 - -# variables modifiable using -var command line switch - -variable mu index -1.25 -variable temp index 2.0 -variable disp index 1.0 -variable lbox index 5.0 - -# global model settings - -units lj -atom_style atomic -pair_style lj/cut 3.0 -pair_modify tail no # turn of to avoid triggering full_energy - -# box - -region box block 0 ${lbox} 0 ${lbox} 0 ${lbox} -create_box 1 box - -# lj parameters - -pair_coeff * * 1.0 1.0 -mass * 1.0 - -# we recommend setting up a dedicated group for gcmc - -group gcmcgroup type 1 - -# gcmc - -fix mygcmc gcmcgroup gcmc 1 100 100 1 29494 ${temp} ${mu} ${disp} max 50 - -# atom count - -variable type1 atom "type==1" -group type1 dynamic gcmcgroup var type1 -variable n1 equal count(type1) - -# averaging - -variable rho equal density -variable p equal press -variable nugget equal 1.0e-8 -variable lambda equal 1.0 -variable muex equal ${mu}-${temp}*ln(density*${lambda}+${nugget}) -fix ave all ave/time 10 100 1000 v_rho v_p v_muex v_n1 ave one file rho_vs_p.dat -variable rhoav equal f_ave[1] -variable pav equal f_ave[2] -variable muexav equal f_ave[3] -variable n1av equal f_ave[4] - -# output - -variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+${nugget}) -variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+${nugget}) -variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+${nugget}) -compute_modify thermo_temp dynamic yes -thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc v_rhoav v_pav v_muexav v_n1av -thermo 1000 - -# run - -run 10000