diff --git a/doc/src/Eqs/fix_gcmc1.jpg b/doc/src/Eqs/fix_gcmc1.jpg new file mode 100644 index 0000000000..158cf8b61f Binary files /dev/null and b/doc/src/Eqs/fix_gcmc1.jpg differ diff --git a/doc/src/Eqs/fix_gcmc1.tex b/doc/src/Eqs/fix_gcmc1.tex new file mode 100644 index 0000000000..c4b0d62527 --- /dev/null +++ b/doc/src/Eqs/fix_gcmc1.tex @@ -0,0 +1,9 @@ +\documentclass[12pt]{article} + +\begin{document} + +\begin{eqnarray*} +\mu &=&\mu^{id} + \mu^{ex} +\end{eqnarray*} + +\end{document} \ No newline at end of file diff --git a/doc/src/Eqs/fix_gcmc2.jpg b/doc/src/Eqs/fix_gcmc2.jpg new file mode 100644 index 0000000000..d054f58430 Binary files /dev/null and b/doc/src/Eqs/fix_gcmc2.jpg differ diff --git a/doc/src/Eqs/fix_gcmc2.tex b/doc/src/Eqs/fix_gcmc2.tex new file mode 100644 index 0000000000..fc4d90355d --- /dev/null +++ b/doc/src/Eqs/fix_gcmc2.tex @@ -0,0 +1,10 @@ +\documentclass[12pt]{article} + +\begin{document} + +\begin{eqnarray*} +\mu^{id} &=& k T \ln{\rho \Lambda^3} \\ +&=& k T \ln{\frac{\phi P \Lambda^3}{k T}} +\end{eqnarray*} + +\end{document} \ No newline at end of file diff --git a/doc/src/Eqs/fix_gcmc3.jpg b/doc/src/Eqs/fix_gcmc3.jpg new file mode 100644 index 0000000000..e87764afd9 Binary files /dev/null and b/doc/src/Eqs/fix_gcmc3.jpg differ diff --git a/doc/src/Eqs/fix_gcmc3.tex b/doc/src/Eqs/fix_gcmc3.tex new file mode 100644 index 0000000000..6aac2b9e1f --- /dev/null +++ b/doc/src/Eqs/fix_gcmc3.tex @@ -0,0 +1,9 @@ +\documentclass[12pt]{article} + +\begin{document} + +\begin{eqnarray*} +\Lambda &=& \sqrt{ \frac{h^2}{2 \pi m k T}} +\end{eqnarray*} + +\end{document} \ No newline at end of file diff --git a/doc/src/dihedral_charmm.txt b/doc/src/dihedral_charmm.txt index 802cf2c4ad..f4637e776d 100644 --- a/doc/src/dihedral_charmm.txt +++ b/doc/src/dihedral_charmm.txt @@ -40,8 +40,11 @@ field. NOTE: The newer {charmmfsh} style was released in March 2017. We recommend it be used instead of the older {charmm} style when running -a simulation with the CHARMM force field. See the discussion below -and more details on the "pair_style charmm"_pair_charmm.html doc page. +a simulation with the CHARMM force field and Coulomb cutoffs, via the +"pair_style lj/charmmfsw/coul/charmmfsh"_pair_charmm.html command. +Otherwise the older {charmm} style is fine to use. See the discussion +below and more details on the "pair_style charmm"_pair_charmm.html doc +page. The following coefficients must be defined for each dihedral type via the "dihedral_coeff"_dihedral_coeff.html command as in the example above, or in @@ -82,13 +85,18 @@ special_bonds 1-4 scaling factor to 0.0 (which is the default). Otherwise 1-4 non-bonded interactions in dihedrals will be computed twice. -For simulations using the CHARMM force field, the difference between -the {charmm} and {charmmfsh} styles is in the computation of the 1-4 -non-bond interactions, if the distance between the two atoms is within -the switching distance of the pairwise potential defined by the -corresponding CHARMM pair style, i.e. between the inner and outer -cutoffs specified for the pair style. See the discussion on the -"CHARMM pair_style"_pair_charmm.html doc page for details. +For simulations using the CHARMM force field with a Coulomb cutoff, +the difference between the {charmm} and {charmmfsh} styles is in the +computation of the 1-4 non-bond interactions, though only if the +distance between the two atoms is within the switching region of the +pairwise potential defined by the corresponding CHARMM pair style, +i.e. between the inner and outer cutoffs specified for the pair style. +The {charmmfsh} style should only be used when using the "pair_style +lj/charmmfsw/coul/charmmfsh"_pair_charmm.html to make the Coulombic +pairwise calculations consistent. Use the {charmm} style with +long-range Coulombics or the older "pair_style +lj/charmm/coul/charmm"_pair_charmm.html command. See the discussion +on the "CHARMM pair_style"_pair_charmm.html doc page for details. Note that for AMBER force fields, which use pair styles with "lj/cut", the special_bonds 1-4 scaling factor should be set to the AMBER diff --git a/doc/src/fix_gcmc.txt b/doc/src/fix_gcmc.txt index 1a419628e9..723e0ec6d9 100644 --- a/doc/src/fix_gcmc.txt +++ b/doc/src/fix_gcmc.txt @@ -122,11 +122,11 @@ If used with "fix nvt"_fix_nh.html, the temperature of the imaginary reservoir, T, should be set to be equivalent to the target temperature used in fix nvt. Otherwise, the imaginary reservoir will not be in thermal equilibrium with the simulation cell. Also, -it is important that the temperature used by fix nvt be dynamic, +it is important that the temperature used by fix nvt be dynamic/dof, which can be achieved as follows: compute mdtemp mdatoms temp -compute_modify mdtemp dynamic yes +compute_modify mdtemp dynamic/dof yes fix mdnvt mdatoms nvt temp 300.0 300.0 10.0 fix_modify mdnvt temp mdtemp :pre @@ -204,10 +204,43 @@ atoms/molecules are assigned to two groups: the default group "all" and the group specified in the fix gcmc command (which can also be "all"). -The gas reservoir pressure can be specified using the {pressure} +The chemical potential is a user-specified input parameter defined +as: + +:c,image(Eqs/fix_gcmc1.jpg) + +The second term mu_ex is the excess chemical potential due to +energetic interactions and is formally zero for the fictitious gas +reservoir but is non-zero for interacting systems. So, while the chemical +potential of the reservoir and the simulation cell are equal, mu_ex is not, +and as a result, the densities of the two are generally quite different. +The first term mu_id is the ideal gas contribution to the chemical potential. +mu_id can be related to the density or pressure of the fictitious gas reservoir by: + +:c,image(Eqs/fix_gcmc2.jpg) + +where k is Boltzman's constant, +T is the user-specified temperature, rho is the number density, +P is the pressure, and phi is the fugacity coefficient. +The constant Lambda is required for dimensional consistency. +For all unit styles except {lj} it is defined as the thermal +de Broglie wavelength + +:c,image(Eqs/fix_gcmc3.jpg) + +where h is Planck's constant, and m is the mass of the exchanged atom or molecule. +For unit style {lj}, Lambda is simply set to the unity. Note that prior to March 2017 +Lambda for unit style {lj} was calculated using the above formula with h set to +the rather specific value of 0.18292026. Chemical potential +under the old definition can be converted to an equivalent value under the new +definition by subtracting 3kTln(Lambda_old). + +As an alternative to specifying mu directly, the ideal gas reservoir +can be defined by its pressure P using the {pressure} keyword, in which case the user-specified chemical potential is -ignored. For non-ideal gas reservoirs, the user may also specify the -fugacity coefficient using the {fugacity_coeff} keyword. +ignored. The user may also specify the +fugacity coefficient phi using the {fugacity_coeff} keyword, which +defaults to unity. The {full_energy} option means that fix GCMC will compute the total potential energy of the entire simulated system. The total system @@ -224,7 +257,8 @@ potential energy calculations, including the following: many-body pair styles hybrid pair styles eam pair styles - triclinic systems + tail corrections + need to include potential energy contributions from other fixes :ul In these cases, LAMMPS will automatically apply the {full_energy} @@ -276,6 +310,13 @@ therefore, you will want to use the current number of atoms is used as a normalizing factor each time temperature is computed. Here is the necessary command: +NOTE: If the density of the cell is initially very small or zero, +and increases to a much larger density after a period of equilibration, +then certain quantities that are only calculated once at the start +(kspace parameters, tail corrections) may no longer be accurate. +The solution is to start a new simulation after the equilibrium +density has been reached. + With some pair_styles, such as "Buckingham"_pair_buck.html, "Born-Mayer-Huggins"_pair_born.html and "ReaxFF"_pair_reax_c.html, two atoms placed close to each other may have an arbitrary large, @@ -366,7 +407,7 @@ referenced by the user for each subsequent fix gcmc command. [Default:] The option defaults are mol = no, maxangle = 10, overlap_cutoff = 0.0, -and full_energy = no, +fugacity_coeff = 1, and full_energy = no, except for the situations where full_energy is required, as listed above. diff --git a/doc/src/pair_charmm.txt b/doc/src/pair_charmm.txt index e1101b4174..94954a2226 100644 --- a/doc/src/pair_charmm.txt +++ b/doc/src/pair_charmm.txt @@ -84,9 +84,9 @@ CHARMM force field. The styles with {charmm} (not {charmmfsw} or {charmmfsh}) in their name are the older, original LAMMPS implementations. They compute the LJ and Coulombic interactions with an energy switching function (esw, -a cubic polynomial, shown in the formula below), which ramps the -energy smoothly to zero between the inner and outer cutoff. This can -cause irregularities in pair-wise forces (due to the discontinuous 2nd +shown in the formula below as S(r)), which ramps the energy smoothly +to zero between the inner and outer cutoff. This can cause +irregularities in pair-wise forces (due to the discontinuous 2nd derivative of energy at the boundaries of the switching region), which in some cases can result in detectable artifacts in an MD simulation. diff --git a/examples/gcmc/in.gcmc.lj b/examples/gcmc/in.gcmc.lj index 9b1f465f04..fc9afdb7f8 100644 --- a/examples/gcmc/in.gcmc.lj +++ b/examples/gcmc/in.gcmc.lj @@ -1,18 +1,23 @@ # 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 available on command line +# variables modifiable using -var command line switch -variable mu index -21.0 -variable disp index 1.0 +variable mu index -1.25 variable temp index 2.0 -variable lbox index 10.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 yes +pair_style lj/cut 3.0 +pair_modify tail no # turn of to avoid triggering full_energy # box @@ -28,15 +33,27 @@ mass * 1.0 fix mygcmc all gcmc 1 100 100 1 29494 ${temp} ${mu} ${disp} +# 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 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] + # output -variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1) -variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1) -variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1) +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 -thermo 100 +thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc v_rhoav v_pav v_muexav +thermo 1000 # run -run 1000 +run 10000 diff --git a/examples/gcmc/log.24Mar17.gcmc.lj.g++.1 b/examples/gcmc/log.24Mar17.gcmc.lj.g++.1 index 5f967192e4..36a9fe885d 100644 --- a/examples/gcmc/log.24Mar17.gcmc.lj.g++.1 +++ b/examples/gcmc/log.24Mar17.gcmc.lj.g++.1 @@ -1,28 +1,35 @@ LAMMPS (17 Mar 2017) +OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:90) + using 1 OpenMP thread(s) per MPI task # 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 available on command line +# variables modifiable using -var command line switch -variable mu index -21.0 -variable disp index 1.0 +variable mu index -1.25 variable temp index 2.0 -variable lbox index 10.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 yes +pair_modify tail no # turn of to avoid triggering full_energy # box region box block 0 ${lbox} 0 ${lbox} 0 ${lbox} -region box block 0 10.0 0 ${lbox} 0 ${lbox} -region box block 0 10.0 0 10.0 0 ${lbox} -region box block 0 10.0 0 10.0 0 10.0 +region box block 0 5.0 0 ${lbox} 0 ${lbox} +region box block 0 5.0 0 5.0 0 ${lbox} +region box block 0 5.0 0 5.0 0 5.0 create_box 1 box -Created orthogonal box = (0 0 0) to (10 10 10) +Created orthogonal box = (0 0 0) to (5 5 5) 1 by 1 by 1 MPI processor grid # lj parameters @@ -34,70 +41,89 @@ mass * 1.0 fix mygcmc all gcmc 1 100 100 1 29494 ${temp} ${mu} ${disp} fix mygcmc all gcmc 1 100 100 1 29494 2.0 ${mu} ${disp} -fix mygcmc all gcmc 1 100 100 1 29494 2.0 -21.0 ${disp} -fix mygcmc all gcmc 1 100 100 1 29494 2.0 -21.0 1.0 +fix mygcmc all gcmc 1 100 100 1 29494 2.0 -1.25 ${disp} +fix mygcmc all gcmc 1 100 100 1 29494 2.0 -1.25 1.0 + +# 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}) +variable muex equal -1.25-${temp}*ln(density*${lambda}+${nugget}) +variable muex equal -1.25-2.0*ln(density*${lambda}+${nugget}) +variable muex equal -1.25-2.0*ln(density*1+${nugget}) +variable muex equal -1.25-2.0*ln(density*1+1e-08) +fix ave all ave/time 10 100 1000 v_rho v_p v_muex 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] # output -variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1) -variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1) -variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1) +variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+${nugget}) +variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+1e-08) +variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+${nugget}) +variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+1e-08) +variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+${nugget}) +variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+1e-08) compute_modify thermo_temp dynamic yes -thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc -thermo 100 +thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc v_rhoav v_pav v_muexav +thermo 1000 # run -run 1000 +run 10000 Neighbor list info ... update every 1 steps, delay 10 steps, check yes max neighbors/atom: 2000, page size: 100000 master list distance cutoff = 3.3 ghost atom cutoff = 3.3 - binsize = 1.65, bins = 7 7 7 + binsize = 1.65, bins = 4 4 4 1 neighbor lists, perpetual/occasional/extra = 1 0 0 (1) pair lj/cut, perpetual attributes: half, newton on pair build: half/bin/atomonly/newton stencil: half/bin/3d/newton bin: standard -Per MPI rank memory allocation (min/avg/max) = 0.4369 | 0.4369 | 0.4369 Mbytes -Step Temp Press PotEng KinEng Density Atoms v_iacc v_dacc v_tacc - 0 0 0 0 -0 0 0 0 0 0 - 100 1.9042848 0.39026453 -1.7692765 2.8466449 0.292 292 0.3619855 0.30247792 0.40278761 - 200 1.8651924 0.47815517 -1.8494955 2.7886155 0.305 305 0.34021109 0.30357196 0.37759189 - 300 2.0626994 0.52068504 -1.8197295 3.0834166 0.291 291 0.32055605 0.3003043 0.36103862 - 400 2.0394818 0.53751435 -1.7636699 3.0482184 0.278 278 0.31698808 0.29995864 0.35441275 - 500 1.9628066 0.54594742 -1.7145336 2.9339513 0.287 287 0.31211861 0.29724228 0.35161407 - 600 1.9845913 0.40846162 -1.8199325 2.9669308 0.299 299 0.30976643 0.29612711 0.34933559 - 700 1.8582606 0.53445462 -1.7869306 2.777974 0.296 296 0.30642103 0.29446478 0.34633665 - 800 2.0340641 0.66057698 -1.7075279 3.0403148 0.283 283 0.30730979 0.29746793 0.34768045 - 900 2.0830765 0.63731971 -1.894775 3.114911 0.322 322 0.30636338 0.29737705 0.34737644 - 1000 1.9688933 0.50024802 -1.7013944 2.9428299 0.281 281 0.3053174 0.29772245 0.34788254 -Loop time of 3.98286 on 1 procs for 1000 steps with 281 atoms +Per MPI rank memory allocation (min/avg/max) = 0.433 | 0.433 | 0.433 Mbytes +Step Temp Press PotEng KinEng Density Atoms v_iacc v_dacc v_tacc v_rhoav v_pav v_muexav + 0 0 0 0 -0 0 0 0 0 0 0 0 0 + 1000 2.4038954 2.1735585 -2.7041368 3.5476844 0.496 62 0.064790036 0.06313096 0.1081294 0.54304 1.4513524 -0.025479219 + 2000 2.0461168 1.1913842 -2.9880181 3.0212194 0.512 64 0.067416408 0.066335853 0.11306166 0.52736 1.3274665 0.034690004 + 3000 1.7930436 1.3788681 -3.2212667 2.6505861 0.552 69 0.067733191 0.066877836 0.1133516 0.5344 1.3834744 0.0070582537 + 4000 1.981449 1.2541054 -2.8222868 2.9217977 0.472 59 0.068546991 0.067856412 0.11442807 0.52504 1.3815629 0.043309657 + 5000 2.0946818 1.0701629 -3.5213291 3.0977688 0.568 71 0.06813743 0.067567891 0.11342906 0.53824 1.4049567 -0.0054539777 + 6000 1.9793484 0.68224187 -3.410211 2.9247088 0.536 67 0.067797628 0.067420108 0.11295333 0.5384 1.401683 -0.0066894359 + 7000 2.1885798 1.6745012 -3.185499 3.2345922 0.544 68 0.068630201 0.068261832 0.11403705 0.5244 1.449239 0.045987399 + 8000 2.2175324 1.5897263 -3.078898 3.2759002 0.528 66 0.068180395 0.067899629 0.11332691 0.53928 1.5488388 -0.01075766 + 9000 1.8610779 1.0396231 -2.923262 2.7465908 0.496 62 0.068346453 0.068028117 0.1134132 0.52912 1.4352871 0.027082544 + 10000 2.1079271 1.1746643 -2.9112062 3.1091925 0.48 60 0.068352878 0.068054948 0.11335434 0.5316 1.4462327 0.018503094 +Loop time of 13.05 on 1 procs for 10000 steps with 60 atoms -Performance: 108464.750 tau/day, 251.076 timesteps/s -99.9% CPU use with 1 MPI tasks x no OpenMP threads +Performance: 331035.016 tau/day, 766.285 timesteps/s +100.0% CPU use with 1 MPI tasks x 1 OpenMP threads MPI task timing breakdown: Section | min time | avg time | max time |%varavg| %total --------------------------------------------------------------- -Pair | 0.10563 | 0.10563 | 0.10563 | 0.0 | 2.65 -Neigh | 0.33428 | 0.33428 | 0.33428 | 0.0 | 8.39 -Comm | 0.027969 | 0.027969 | 0.027969 | 0.0 | 0.70 -Output | 0.00017285 | 0.00017285 | 0.00017285 | 0.0 | 0.00 -Modify | 3.5096 | 3.5096 | 3.5096 | 0.0 | 88.12 -Other | | 0.005197 | | | 0.13 +Pair | 0.37239 | 0.37239 | 0.37239 | 0.0 | 2.85 +Neigh | 0.94764 | 0.94764 | 0.94764 | 0.0 | 7.26 +Comm | 0.092473 | 0.092473 | 0.092473 | 0.0 | 0.71 +Output | 0.00023365 | 0.00023365 | 0.00023365 | 0.0 | 0.00 +Modify | 11.627 | 11.627 | 11.627 | 0.0 | 89.09 +Other | | 0.01054 | | | 0.08 -Nlocal: 281 ave 281 max 281 min +Nlocal: 60 ave 60 max 60 min Histogram: 1 0 0 0 0 0 0 0 0 0 -Nghost: 977 ave 977 max 977 min +Nghost: 663 ave 663 max 663 min Histogram: 1 0 0 0 0 0 0 0 0 0 -Neighs: 5902 ave 5902 max 5902 min +Neighs: 2133 ave 2133 max 2133 min Histogram: 1 0 0 0 0 0 0 0 0 0 -Total # of neighbors = 5902 -Ave neighs/atom = 21.0036 -Neighbor list builds = 1000 +Total # of neighbors = 2133 +Ave neighs/atom = 35.55 +Neighbor list builds = 10000 Dangerous builds = 0 -Total wall time: 0:00:03 +Total wall time: 0:00:13 diff --git a/examples/gcmc/log.24Mar17.gcmc.lj.g++.4 b/examples/gcmc/log.24Mar17.gcmc.lj.g++.4 index db4cc7fd9a..8694d8b95e 100644 --- a/examples/gcmc/log.24Mar17.gcmc.lj.g++.4 +++ b/examples/gcmc/log.24Mar17.gcmc.lj.g++.4 @@ -1,28 +1,35 @@ LAMMPS (17 Mar 2017) +OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:90) + using 1 OpenMP thread(s) per MPI task # 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 available on command line +# variables modifiable using -var command line switch -variable mu index -21.0 -variable disp index 1.0 +variable mu index -1.25 variable temp index 2.0 -variable lbox index 10.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 yes +pair_modify tail no # turn of to avoid triggering full_energy # box region box block 0 ${lbox} 0 ${lbox} 0 ${lbox} -region box block 0 10.0 0 ${lbox} 0 ${lbox} -region box block 0 10.0 0 10.0 0 ${lbox} -region box block 0 10.0 0 10.0 0 10.0 +region box block 0 5.0 0 ${lbox} 0 ${lbox} +region box block 0 5.0 0 5.0 0 ${lbox} +region box block 0 5.0 0 5.0 0 5.0 create_box 1 box -Created orthogonal box = (0 0 0) to (10 10 10) +Created orthogonal box = (0 0 0) to (5 5 5) 1 by 2 by 2 MPI processor grid # lj parameters @@ -34,70 +41,89 @@ mass * 1.0 fix mygcmc all gcmc 1 100 100 1 29494 ${temp} ${mu} ${disp} fix mygcmc all gcmc 1 100 100 1 29494 2.0 ${mu} ${disp} -fix mygcmc all gcmc 1 100 100 1 29494 2.0 -21.0 ${disp} -fix mygcmc all gcmc 1 100 100 1 29494 2.0 -21.0 1.0 +fix mygcmc all gcmc 1 100 100 1 29494 2.0 -1.25 ${disp} +fix mygcmc all gcmc 1 100 100 1 29494 2.0 -1.25 1.0 + +# 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}) +variable muex equal -1.25-${temp}*ln(density*${lambda}+${nugget}) +variable muex equal -1.25-2.0*ln(density*${lambda}+${nugget}) +variable muex equal -1.25-2.0*ln(density*1+${nugget}) +variable muex equal -1.25-2.0*ln(density*1+1e-08) +fix ave all ave/time 10 100 1000 v_rho v_p v_muex 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] # output -variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1) -variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1) -variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1) +variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+${nugget}) +variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+1e-08) +variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+${nugget}) +variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+1e-08) +variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+${nugget}) +variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+1e-08) compute_modify thermo_temp dynamic yes -thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc -thermo 100 +thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc v_rhoav v_pav v_muexav +thermo 1000 # run -run 1000 +run 10000 Neighbor list info ... update every 1 steps, delay 10 steps, check yes max neighbors/atom: 2000, page size: 100000 master list distance cutoff = 3.3 ghost atom cutoff = 3.3 - binsize = 1.65, bins = 7 7 7 + binsize = 1.65, bins = 4 4 4 1 neighbor lists, perpetual/occasional/extra = 1 0 0 (1) pair lj/cut, perpetual attributes: half, newton on pair build: half/bin/atomonly/newton stencil: half/bin/3d/newton bin: standard -Per MPI rank memory allocation (min/avg/max) = 0.434 | 0.434 | 0.434 Mbytes -Step Temp Press PotEng KinEng Density Atoms v_iacc v_dacc v_tacc - 0 0 0 0 -0 0 0 0 0 0 - 100 2.0328045 0.58661762 -1.6812724 3.0385824 0.287 287 0.35917318 0.30067507 0.38663622 - 200 1.9594279 0.50682399 -1.7308396 2.9287927 0.284 284 0.33788365 0.30337335 0.37300293 - 300 2.0602937 0.7028247 -1.9278541 3.0806296 0.315 315 0.31882007 0.29697498 0.36167185 - 400 1.995183 0.4328246 -1.8715454 2.983026 0.307 307 0.31527654 0.29681901 0.35673374 - 500 2.1390101 0.48232215 -1.554138 3.1960306 0.257 257 0.31372975 0.30003067 0.35558858 - 600 2.0584244 0.4929049 -1.6995569 3.0767263 0.283 283 0.31114213 0.29801665 0.35160109 - 700 1.9155066 0.49654243 -1.5770611 2.8624174 0.265 265 0.31056419 0.29944173 0.35157337 - 800 2.0883562 0.52731947 -1.8261112 3.1220925 0.3 300 0.30730979 0.29704354 0.34898892 - 900 2.0470677 0.5605993 -2.0130053 3.0610656 0.322 322 0.30484441 0.29586719 0.34678883 - 1000 2.004135 0.50642204 -1.6956257 2.9955798 0.283 283 0.30396929 0.29634309 0.34770304 -Loop time of 3.688 on 4 procs for 1000 steps with 283 atoms +Per MPI rank memory allocation (min/avg/max) = 0.4477 | 0.4477 | 0.4477 Mbytes +Step Temp Press PotEng KinEng Density Atoms v_iacc v_dacc v_tacc v_rhoav v_pav v_muexav + 0 0 0 0 -0 0 0 0 0 0 0 0 0 + 1000 1.956397 1.7699101 -2.7889468 2.8864874 0.488 61 0.068894746 0.067229075 0.1141726 0.53288 1.3832798 0.013392866 + 2000 2.040943 0.56060899 -2.8001647 3.0077055 0.456 57 0.069858594 0.068831934 0.11629114 0.5232 1.3587174 0.049995794 + 3000 2.0004866 1.5736515 -3.3098044 2.9572411 0.552 69 0.069594029 0.068727791 0.11592543 0.53096 1.4129434 0.020022578 + 4000 2.1127942 2.642809 -2.8865084 3.1211733 0.528 66 0.070268697 0.069533235 0.11693806 0.52424 1.3444615 0.046884078 + 5000 2.3663648 1.354269 -3.1917346 3.4957662 0.528 66 0.070519633 0.069960064 0.11710321 0.52688 1.3595814 0.036270867 + 6000 1.9224136 0.82756699 -3.1965 2.839257 0.52 65 0.06985018 0.069474645 0.11628632 0.536 1.47062 0.00141549 + 7000 2.0266192 1.5593811 -2.9972341 2.9931606 0.52 65 0.070244693 0.069880791 0.11666541 0.52528 1.3246332 0.040754793 + 8000 1.7790467 1.8680568 -2.8028819 2.6275151 0.52 65 0.070454494 0.070172368 0.11736806 0.524 1.4213649 0.047985191 + 9000 1.7968847 1.3195587 -3.261001 2.6550983 0.536 67 0.069952011 0.069618327 0.11650087 0.53904 1.4624201 -0.01069837 + 10000 2.1566109 1.1015729 -3.4999837 3.1880335 0.552 69 0.069603309 0.069284134 0.11625548 0.53128 1.3587249 0.02075238 +Loop time of 13.0611 on 4 procs for 10000 steps with 69 atoms -Performance: 117136.751 tau/day, 271.150 timesteps/s -99.2% CPU use with 4 MPI tasks x no OpenMP threads +Performance: 330753.007 tau/day, 765.632 timesteps/s +99.7% CPU use with 4 MPI tasks x 1 OpenMP threads MPI task timing breakdown: Section | min time | avg time | max time |%varavg| %total --------------------------------------------------------------- -Pair | 0.024644 | 0.026027 | 0.027483 | 0.6 | 0.71 -Neigh | 0.085449 | 0.088998 | 0.092893 | 0.9 | 2.41 -Comm | 0.045756 | 0.051296 | 0.056578 | 1.7 | 1.39 -Output | 0.00028491 | 0.00030857 | 0.00035262 | 0.0 | 0.01 -Modify | 3.5189 | 3.5191 | 3.5194 | 0.0 | 95.42 -Other | | 0.002221 | | | 0.06 +Pair | 0.08888 | 0.09443 | 0.099889 | 1.4 | 0.72 +Neigh | 0.27721 | 0.28532 | 0.29177 | 1.1 | 2.18 +Comm | 0.27648 | 0.28875 | 0.30268 | 1.9 | 2.21 +Output | 0.00029635 | 0.00043058 | 0.00048113 | 0.0 | 0.00 +Modify | 12.384 | 12.384 | 12.384 | 0.0 | 94.82 +Other | | 0.008055 | | | 0.06 -Nlocal: 70.75 ave 77 max 68 min -Histogram: 1 2 0 0 0 0 0 0 0 1 -Nghost: 514.25 ave 520 max 507 min -Histogram: 1 0 0 0 1 0 0 1 0 1 -Neighs: 1483.5 ave 1715 max 1359 min -Histogram: 2 0 0 1 0 0 0 0 0 1 +Nlocal: 17.25 ave 23 max 10 min +Histogram: 1 0 0 0 0 0 2 0 0 1 +Nghost: 506.5 ave 519 max 490 min +Histogram: 1 0 1 0 0 0 0 0 0 2 +Neighs: 705.75 ave 998 max 369 min +Histogram: 1 0 0 0 0 1 1 0 0 1 -Total # of neighbors = 5934 -Ave neighs/atom = 20.9682 -Neighbor list builds = 1000 +Total # of neighbors = 2823 +Ave neighs/atom = 40.913 +Neighbor list builds = 10000 Dangerous builds = 0 -Total wall time: 0:00:03 +Total wall time: 0:00:13 diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp index a7483a1d00..407c980729 100644 --- a/src/MC/fix_gcmc.cpp +++ b/src/MC/fix_gcmc.cpp @@ -432,7 +432,8 @@ void FixGCMC::init() (force->pair == NULL) || (force->pair->single_enable == 0) || (force->pair_match("hybrid",0)) || - (force->pair_match("eam",0)) + (force->pair_match("eam",0)) || + (force->pair->tail_flag) ) { full_flag = true; if (comm->me == 0) @@ -618,13 +619,18 @@ void FixGCMC::init() } // compute beta, lambda, sigma, and the zz factor - + // For LJ units, lambda=1 beta = 1.0/(force->boltz*reservoir_temperature); - double lambda = sqrt(force->hplanck*force->hplanck/ - (2.0*MY_PI*gas_mass*force->mvv2e* + if (strcmp(update->unit_style,"lj") == 0) + zz = exp(beta*chemical_potential); + else { + double lambda = sqrt(force->hplanck*force->hplanck/ + (2.0*MY_PI*gas_mass*force->mvv2e* force->boltz*reservoir_temperature)); + zz = exp(beta*chemical_potential)/(pow(lambda,3.0)); + } + sigma = sqrt(force->boltz*reservoir_temperature*tfac_insert/gas_mass/force->mvv2e); - zz = exp(beta*chemical_potential)/(pow(lambda,3.0)); if (pressure_flag) zz = pressure*fugacity_coeff*beta/force->nktv2p; imagezero = ((imageint) IMGMAX << IMG2BITS) | diff --git a/src/USER-MANIFOLD/fix_manifoldforce.cpp b/src/USER-MANIFOLD/fix_manifoldforce.cpp index 2fac8e40df..e66b7c9fc2 100644 --- a/src/USER-MANIFOLD/fix_manifoldforce.cpp +++ b/src/USER-MANIFOLD/fix_manifoldforce.cpp @@ -84,8 +84,8 @@ FixManifoldForce::FixManifoldForce(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,msg); } - *(ptr_m->get_params()) = new double[nvars]; - if( ptr_m->get_params() == NULL ){ + ptr_m->params = new double[nvars]; + if( ptr_m->params == NULL ){ error->all(FLERR,"Parameter pointer was NULL!"); } @@ -94,7 +94,7 @@ FixManifoldForce::FixManifoldForce(LAMMPS *lmp, int narg, char **arg) : // and sets the values of those arguments that were _not_ // equal style vars (so that they are not overwritten each time step). - double *params = *(ptr_m->get_params()); + double *params = ptr_m->params; for( int i = 0; i < nvars; ++i ){ if( was_var( arg[i+4] ) ) error->all(FLERR,"Equal-style variables not allowed with fix manifoldforce"); diff --git a/src/USER-MANIFOLD/fix_nve_manifold_rattle.cpp b/src/USER-MANIFOLD/fix_nve_manifold_rattle.cpp index e27762a7d5..0d48e145e0 100644 --- a/src/USER-MANIFOLD/fix_nve_manifold_rattle.cpp +++ b/src/USER-MANIFOLD/fix_nve_manifold_rattle.cpp @@ -131,11 +131,11 @@ FixNVEManifoldRattle::FixNVEManifoldRattle( LAMMPS *lmp, int &narg, char **arg, strcpy( tstrs[i], arg[i+6] + offset ); } - *ptr_m->get_params() = new double[nvars]; - if( !(*ptr_m->get_params()) ) error->all(FLERR,"Failed to allocate params!"); + ptr_m->params = new double[nvars]; + if( !ptr_m->params ) error->all(FLERR,"Failed to allocate params!"); for( int i = 0; i < nvars; ++i ){ // If param i was variable type, it will be set later... - (*ptr_m->get_params())[i] = is_var[i] ? 0.0 : force->numeric( FLERR, arg[i+6] ); + ptr_m->params[i] = is_var[i] ? 0.0 : force->numeric( FLERR, arg[i+6] ); } ptr_m->post_param_init(); @@ -262,6 +262,7 @@ void FixNVEManifoldRattle::init() void FixNVEManifoldRattle::update_var_params() { + if( nevery > 0 ){ stats.x_iters = 0; stats.v_iters = 0; @@ -270,7 +271,8 @@ void FixNVEManifoldRattle::update_var_params() stats.v_iters_per_atom = 0.0; } - double **ptr_params = ptr_m->get_params(); + double *ptr_params = ptr_m->params; + for( int i = 0; i < nvars; ++i ){ if( is_var[i] ){ tvars[i] = input->variable->find(tstrs[i]); @@ -281,8 +283,8 @@ void FixNVEManifoldRattle::update_var_params() if( input->variable->equalstyle(tvars[i]) ){ tstyle[i] = EQUAL; double new_val = input->variable->compute_equal(tvars[i]); - // fprintf( stdout, "New value of var %d is now %f\n", i+1, new_val ); - *(ptr_params[i]) = new_val; + + ptr_params[i] = new_val; }else{ error->all(FLERR, "Variable for fix nve/manifold/rattle is invalid style"); diff --git a/src/USER-MANIFOLD/manifold.h b/src/USER-MANIFOLD/manifold.h index 85692a8b6b..d0ffa214ac 100644 --- a/src/USER-MANIFOLD/manifold.h +++ b/src/USER-MANIFOLD/manifold.h @@ -63,12 +63,12 @@ namespace user_manifold { virtual void set_atom_id( tagint a_id ){} virtual int nparams() = 0; - double **get_params(){ return ¶ms; }; + // double *get_params(){ return params; }; // Overload if any initialization depends on params: virtual void post_param_init(){} virtual void checkup(){} // Some diagnostics... - protected: + double *params; }; diff --git a/src/update.cpp b/src/update.cpp index 35492006cc..5599dc6c88 100644 --- a/src/update.cpp +++ b/src/update.cpp @@ -128,7 +128,7 @@ void Update::set_units(const char *style) if (strcmp(style,"lj") == 0) { force->boltz = 1.0; - force->hplanck = 0.18292026; // using LJ parameters for argon + force->hplanck = 1.0; force->mvv2e = 1.0; force->ftm2v = 1.0; force->mv2d = 1.0;