mirror of https://github.com/lammps/lammps.git
fix gcmc units change for chemical potential
This commit is contained in:
Steve Plimpton
parent
3dfe4505dd
commit
111d350a22
Binary file not shown.
|
After Width: | Height: | Size: 5.5 KiB |
|
|
@ -0,0 +1,9 @@
|
|||
\documentclass[12pt]{article}
|
||||
|
||||
\begin{document}
|
||||
|
||||
\begin{eqnarray*}
|
||||
\mu &=&\mu^{id} + \mu^{ex}
|
||||
\end{eqnarray*}
|
||||
|
||||
\end{document}
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 10 KiB |
|
|
@ -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}
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 7.3 KiB |
|
|
@ -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}
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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.
|
||||
|
||||
|
|
|
|||
|
|
@ -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.
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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) |
|
||||
|
|
|
|||
|
|
@ -790,7 +790,7 @@ void lammps_gather_atoms(void *ptr, char *name,
|
|||
for (i = 0; i < nlocal; i++) {
|
||||
offset = count*(tag[i]-1);
|
||||
for (j = 0; j < count; j++)
|
||||
copy[offset++] = array[i][0];
|
||||
copy[offset++] = array[i][j];
|
||||
}
|
||||
|
||||
MPI_Allreduce(copy,data,count*natoms,MPI_INT,MPI_SUM,lmp->world);
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
Loading…
Reference in New Issue