ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'master' into master-small-fixes

Browse Source
This commit is contained in:
Axel Kohlmeyer
2010-06-19 07:34:10 -04:00
parent 18bfaf1b26 b3de0db9db
commit d46562c6ef
37 changed files with 569 additions and 251 deletions
Download Patch File Download Diff File Expand all files Collapse all files

20
doc/Section_commands.html
View File

@ -323,16 +323,16 @@ in the command's documentation.
of each style or click on the style itself for a full description:
</P>
<DIV ALIGN=center><TABLE BORDER=1 >
<TR ALIGN="center"><TD ><A HREF = "fix_addforce.html">addforce</A></TD><TD ><A HREF = "fix_aveforce.html">aveforce</A></TD><TD ><A HREF = "fix_ave_atom.html">ave/atom</A></TD><TD ><A HREF = "fix_ave_histo.html">ave/histo</A></TD><TD ><A HREF = "fix_ave_spatial.html">ave/spatial</A></TD><TD ><A HREF = "fix_ave_time.html">ave/time</A></TD><TD ><A HREF = "fix_bond_break.html">bond/break</A></TD><TD ><A HREF = "fix_bond_create.html">bond/create</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "fix_bond_swap.html">bond/swap</A></TD><TD ><A HREF = "fix_box_relax.html">box/relax</A></TD><TD ><A HREF = "fix_deform.html">deform</A></TD><TD ><A HREF = "fix_deposit.html">deposit</A></TD><TD ><A HREF = "fix_drag.html">drag</A></TD><TD ><A HREF = "fix_dt_reset.html">dt/reset</A></TD><TD ><A HREF = "fix_efield.html">efield</A></TD><TD ><A HREF = "fix_enforce2d.html">enforce2d</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "fix_evaporate.html">evaporate</A></TD><TD ><A HREF = "fix_freeze.html">freeze</A></TD><TD ><A HREF = "fix_gravity.html">gravity</A></TD><TD ><A HREF = "fix_heat.html">heat</A></TD><TD ><A HREF = "fix_indent.html">indent</A></TD><TD ><A HREF = "fix_langevin.html">langevin</A></TD><TD ><A HREF = "fix_lineforce.html">lineforce</A></TD><TD ><A HREF = "fix_momentum.html">momentum</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "fix_move.html">move</A></TD><TD ><A HREF = "fix_msst.html">msst</A></TD><TD ><A HREF = "fix_nh.html">nph</A></TD><TD ><A HREF = "fix_nph_asphere.html">nph/asphere</A></TD><TD ><A HREF = "fix_nph_sphere.html">nph/sphere</A></TD><TD ><A HREF = "fix_nh.html">npt</A></TD><TD ><A HREF = "fix_npt_asphere.html">npt/asphere</A></TD><TD ><A HREF = "fix_npt_sphere.html">npt/sphere</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "fix_nve.html">nve</A></TD><TD ><A HREF = "fix_nve_asphere.html">nve/asphere</A></TD><TD ><A HREF = "fix_nve_limit.html">nve/limit</A></TD><TD ><A HREF = "fix_nve_noforce.html">nve/noforce</A></TD><TD ><A HREF = "fix_nve_sphere.html">nve/sphere</A></TD><TD ><A HREF = "fix_nh.html">nvt</A></TD><TD ><A HREF = "fix_nvt_asphere.html">nvt/asphere</A></TD><TD ><A HREF = "fix_nvt_sllod.html">nvt/sllod</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "fix_nvt_sphere.html">nvt/sphere</A></TD><TD ><A HREF = "fix_orient_fcc.html">orient/fcc</A></TD><TD ><A HREF = "fix_planeforce.html">planeforce</A></TD><TD ><A HREF = "fix_poems.html">poems</A></TD><TD ><A HREF = "fix_pour.html">pour</A></TD><TD ><A HREF = "fix_press_berendsen.html">press/berendsen</A></TD><TD ><A HREF = "fix_print.html">print</A></TD><TD ><A HREF = "fix_reax_bonds.html">reax/bonds</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "fix_recenter.html">recenter</A></TD><TD ><A HREF = "fix_rigid.html">rigid</A></TD><TD ><A HREF = "fix_rigid.html">rigid/nve</A></TD><TD ><A HREF = "fix_rigid.html">rigid/nvt</A></TD><TD ><A HREF = "fix_setforce.html">setforce</A></TD><TD ><A HREF = "fix_shake.html">shake</A></TD><TD ><A HREF = "fix_spring.html">spring</A></TD><TD ><A HREF = "fix_spring_rg.html">spring/rg</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "fix_spring_self.html">spring/self</A></TD><TD ><A HREF = "fix_store_force.html">store/force</A></TD><TD ><A HREF = "fix_store_state.html">store/state</A></TD><TD ><A HREF = "fix_temp_berendsen.html">temp/berendsen</A></TD><TD ><A HREF = "fix_temp_rescale.html">temp/rescale</A></TD><TD ><A HREF = "fix_thermal_conductivity.html">thermal/conductivity</A></TD><TD ><A HREF = "fix_tmd.html">tmd</A></TD><TD ><A HREF = "fix_ttm.html">ttm</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "fix_viscosity.html">viscosity</A></TD><TD ><A HREF = "fix_viscous.html">viscous</A></TD><TD ><A HREF = "fix_wall.html">wall/colloid</A></TD><TD ><A HREF = "fix_wall_gran.html">wall/gran</A></TD><TD ><A HREF = "fix_wall.html">wall/harmonic</A></TD><TD ><A HREF = "fix_wall.html">wall/lj126</A></TD><TD ><A HREF = "fix_wall.html">wall/lj93</A></TD><TD ><A HREF = "fix_wall_reflect.html">wall/reflect</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "fix_wall_region.html">wall/region</A>
<TR ALIGN="center"><TD ><A HREF = "fix_adapt.html">adapt</A></TD><TD ><A HREF = "fix_addforce.html">addforce</A></TD><TD ><A HREF = "fix_aveforce.html">aveforce</A></TD><TD ><A HREF = "fix_ave_atom.html">ave/atom</A></TD><TD ><A HREF = "fix_ave_histo.html">ave/histo</A></TD><TD ><A HREF = "fix_ave_spatial.html">ave/spatial</A></TD><TD ><A HREF = "fix_ave_time.html">ave/time</A></TD><TD ><A HREF = "fix_bond_break.html">bond/break</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "fix_bond_create.html">bond/create</A></TD><TD ><A HREF = "fix_bond_swap.html">bond/swap</A></TD><TD ><A HREF = "fix_box_relax.html">box/relax</A></TD><TD ><A HREF = "fix_deform.html">deform</A></TD><TD ><A HREF = "fix_deposit.html">deposit</A></TD><TD ><A HREF = "fix_drag.html">drag</A></TD><TD ><A HREF = "fix_dt_reset.html">dt/reset</A></TD><TD ><A HREF = "fix_efield.html">efield</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "fix_enforce2d.html">enforce2d</A></TD><TD ><A HREF = "fix_evaporate.html">evaporate</A></TD><TD ><A HREF = "fix_freeze.html">freeze</A></TD><TD ><A HREF = "fix_gravity.html">gravity</A></TD><TD ><A HREF = "fix_heat.html">heat</A></TD><TD ><A HREF = "fix_indent.html">indent</A></TD><TD ><A HREF = "fix_langevin.html">langevin</A></TD><TD ><A HREF = "fix_lineforce.html">lineforce</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "fix_momentum.html">momentum</A></TD><TD ><A HREF = "fix_move.html">move</A></TD><TD ><A HREF = "fix_msst.html">msst</A></TD><TD ><A HREF = "fix_nh.html">nph</A></TD><TD ><A HREF = "fix_nph_asphere.html">nph/asphere</A></TD><TD ><A HREF = "fix_nph_sphere.html">nph/sphere</A></TD><TD ><A HREF = "fix_nh.html">npt</A></TD><TD ><A HREF = "fix_npt_asphere.html">npt/asphere</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "fix_npt_sphere.html">npt/sphere</A></TD><TD ><A HREF = "fix_nve.html">nve</A></TD><TD ><A HREF = "fix_nve_asphere.html">nve/asphere</A></TD><TD ><A HREF = "fix_nve_limit.html">nve/limit</A></TD><TD ><A HREF = "fix_nve_noforce.html">nve/noforce</A></TD><TD ><A HREF = "fix_nve_sphere.html">nve/sphere</A></TD><TD ><A HREF = "fix_nh.html">nvt</A></TD><TD ><A HREF = "fix_nvt_asphere.html">nvt/asphere</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "fix_nvt_sllod.html">nvt/sllod</A></TD><TD ><A HREF = "fix_nvt_sphere.html">nvt/sphere</A></TD><TD ><A HREF = "fix_orient_fcc.html">orient/fcc</A></TD><TD ><A HREF = "fix_planeforce.html">planeforce</A></TD><TD ><A HREF = "fix_poems.html">poems</A></TD><TD ><A HREF = "fix_pour.html">pour</A></TD><TD ><A HREF = "fix_press_berendsen.html">press/berendsen</A></TD><TD ><A HREF = "fix_print.html">print</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "fix_reax_bonds.html">reax/bonds</A></TD><TD ><A HREF = "fix_recenter.html">recenter</A></TD><TD ><A HREF = "fix_rigid.html">rigid</A></TD><TD ><A HREF = "fix_rigid.html">rigid/nve</A></TD><TD ><A HREF = "fix_rigid.html">rigid/nvt</A></TD><TD ><A HREF = "fix_setforce.html">setforce</A></TD><TD ><A HREF = "fix_shake.html">shake</A></TD><TD ><A HREF = "fix_spring.html">spring</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "fix_spring_rg.html">spring/rg</A></TD><TD ><A HREF = "fix_spring_self.html">spring/self</A></TD><TD ><A HREF = "fix_store_force.html">store/force</A></TD><TD ><A HREF = "fix_store_state.html">store/state</A></TD><TD ><A HREF = "fix_temp_berendsen.html">temp/berendsen</A></TD><TD ><A HREF = "fix_temp_rescale.html">temp/rescale</A></TD><TD ><A HREF = "fix_thermal_conductivity.html">thermal/conductivity</A></TD><TD ><A HREF = "fix_tmd.html">tmd</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "fix_ttm.html">ttm</A></TD><TD ><A HREF = "fix_viscosity.html">viscosity</A></TD><TD ><A HREF = "fix_viscous.html">viscous</A></TD><TD ><A HREF = "fix_wall.html">wall/colloid</A></TD><TD ><A HREF = "fix_wall_gran.html">wall/gran</A></TD><TD ><A HREF = "fix_wall.html">wall/harmonic</A></TD><TD ><A HREF = "fix_wall.html">wall/lj126</A></TD><TD ><A HREF = "fix_wall.html">wall/lj93</A></TD></TR>
<TR ALIGN="center"><TD ><A HREF = "fix_wall_reflect.html">wall/reflect</A></TD><TD ><A HREF = "fix_wall_region.html">wall/region</A>
</TD></TR></TABLE></DIV>
<P>These are fix styles contributed by users, which can be used if

1
doc/Section_commands.txt
View File

@ -382,6 +382,7 @@ Fix styles :h4
See the "fix"_fix.html command for one-line descriptions
of each style or click on the style itself for a full description:
"adapt"_fix_adapt.html,
"addforce"_fix_addforce.html,
"aveforce"_fix_aveforce.html,
"ave/atom"_fix_ave_atom.html,

23
doc/compute_property_local.html
View File

@ -21,13 +21,15 @@
<LI>input = one or more attributes
<PRE> possible attributes = patom1 patom2
<PRE> possible attributes = natom1 natom2
patom1 patom2
batom1 batom2 btype
aatom1 aatom2 aatom3 atype
datom1 datom2 datom3 dtype
iatom1 iatom2 iatom3 itype
</PRE>
<PRE> patom1, patom2 = IDs of 2 atoms in each pair
<PRE> natom1, natom2 = IDs of 2 atoms in each pair (within neighbor cutoff)
patom1, patom2 = IDs of 2 atoms in each pair (within force cutoff)
batom1, batom2 = IDs of 2 atoms in each bond
btype = bond type of each bond
aatom1, aatom2, aatom3 = IDs of 3 atoms in each angle
@ -63,10 +65,13 @@ property/local command.
<P>If the inputs are pair attributes, the local data is generated by
looping over the pairwise neighbor list. Info about an individual
pairwise interaction will only be included if both atoms in the pair
are in the specified compute group, and if the current pairwise
distance is less than the force cutoff distance for that interaction,
as defined by the <A HREF = "pair_style.html">pair_style</A> and
<A HREF = "pair_coeff.html">pair_coeff</A> commands.
are in the specified compute group. For <I>natom1</I> and <I>natom2</I>, all
atom pairs in the neighbor list are considered (out to the neighbor
cutoff = force cutoff + <A HREF = "neighbor.html">neighbor skin</A>). For <I>patom1</I>
and <I>patom2</I>, the distance between the atoms must be less than the
force cutoff distance for that pair to be included, as defined by the
<A HREF = "pair_style.html">pair_style</A> and <A HREF = "pair_coeff.html">pair_coeff</A>
commands.
</P>
<P>If the inputs are bond, angle, etc attributes, the local data is
generated by looping over all the atoms owned on a processor and
@ -84,9 +89,9 @@ bond/local</A> command can be combined with bond
atom indices from this command and output by the <A HREF = "dump.html">dump
local</A> command in a consistent way.
</P>
<P>The <I>patom1</I> and <I>patom2</I> attributes refer to the atom IDs of the 2
atoms in each pairwise interaction computed by the
<A HREF = "pair_style.html">pair_style</A> command.
<P>The <I>natom1</I> and <I>natom2</I>, or <I>patom1</I> and <I>patom2</I> attributes refer
to the atom IDs of the 2 atoms in each pairwise interaction computed
by the <A HREF = "pair_style.html">pair_style</A> command.
</P>
<P>IMPORTANT NOTE: For pairs, if two atoms I,J are involved in 1-2, 1-3,
1-4 interactions within the molecular topology, their pairwise

23
doc/compute_property_local.txt
View File

@ -15,13 +15,15 @@ compute ID group-ID property/local input1 input2 ... :pre
ID, group-ID are documented in "compute"_compute.html command :ulb,l
property/local = style name of this compute command :l
input = one or more attributes :l
possible attributes = patom1 patom2
possible attributes = natom1 natom2
patom1 patom2
batom1 batom2 btype
aatom1 aatom2 aatom3 atype
datom1 datom2 datom3 dtype
iatom1 iatom2 iatom3 itype :pre
patom1, patom2 = IDs of 2 atoms in each pair
natom1, natom2 = IDs of 2 atoms in each pair (within neighbor cutoff)
patom1, patom2 = IDs of 2 atoms in each pair (within force cutoff)
batom1, batom2 = IDs of 2 atoms in each bond
btype = bond type of each bond
aatom1, aatom2, aatom3 = IDs of 3 atoms in each angle
@ -56,10 +58,13 @@ property/local command.
If the inputs are pair attributes, the local data is generated by
looping over the pairwise neighbor list. Info about an individual
pairwise interaction will only be included if both atoms in the pair
are in the specified compute group, and if the current pairwise
distance is less than the force cutoff distance for that interaction,
as defined by the "pair_style"_pair_style.html and
"pair_coeff"_pair_coeff.html commands.
are in the specified compute group. For {natom1} and {natom2}, all
atom pairs in the neighbor list are considered (out to the neighbor
cutoff = force cutoff + "neighbor skin"_neighbor.html). For {patom1}
and {patom2}, the distance between the atoms must be less than the
force cutoff distance for that pair to be included, as defined by the
"pair_style"_pair_style.html and "pair_coeff"_pair_coeff.html
commands.
If the inputs are bond, angle, etc attributes, the local data is
generated by looping over all the atoms owned on a processor and
@ -77,9 +82,9 @@ bond/local"_compute_bond_local.html command can be combined with bond
atom indices from this command and output by the "dump
local"_dump.html command in a consistent way.
The {patom1} and {patom2} attributes refer to the atom IDs of the 2
atoms in each pairwise interaction computed by the
"pair_style"_pair_style.html command.
The {natom1} and {natom2}, or {patom1} and {patom2} attributes refer
to the atom IDs of the 2 atoms in each pairwise interaction computed
by the "pair_style"_pair_style.html command.
IMPORTANT NOTE: For pairs, if two atoms I,J are involved in 1-2, 1-3,
1-4 interactions within the molecular topology, their pairwise

3
doc/fix.html
View File

@ -160,7 +160,8 @@ for further info.
arguments and what it does, as listed below. Here is an alphabetic
list of fix styles available in LAMMPS:
</P>
<UL><LI><A HREF = "fix_addforce.html">addforce</A> - add a force to each atom
<UL><LI><A HREF = "fix_adapt.html">adapt</A> - change a simulation parameter over time
<LI><A HREF = "fix_addforce.html">addforce</A> - add a force to each atom
<LI><A HREF = "fix_aveforce.html">aveforce</A> - add an averaged force to each atom
<LI><A HREF = "fix_ave_atom.html">ave/atom</A> - compute per-atom time-averaged quantities
<LI><A HREF = "fix_ave_atom.html">ave/histo</A> - compute/output time-averaged histograms

1
doc/fix.txt
View File

@ -155,6 +155,7 @@ Each fix style has its own documentation page which describes its
arguments and what it does, as listed below. Here is an alphabetic
list of fix styles available in LAMMPS:
"adapt"_fix_adapt.html - change a simulation parameter over time
"addforce"_fix_addforce.html - add a force to each atom
"aveforce"_fix_aveforce.html - add an averaged force to each atom
"ave/atom"_fix_ave_atom.html - compute per-atom time-averaged quantities

4
doc/pair_lubricate.html
View File

@ -43,6 +43,10 @@ their relative velocities in the presence of a background solvent with
viscosity mu. Note that this is dynamic viscosity which has units of
mass/distance/time, not kinematic viscosity.
</P>
<P>The viscosity mu can be varied in a time-dependent manner over the
course of a simluation. See the <A HREF = "fix_adapt.html">fix adapt</A> command
for details.
</P>
<P>Rc is the outer cutoff specified in the pair_style command, the
translational velocities of the 2 particles are v1 and v2, the angular
velocities are w1 and w2, and n is the unit vector in the direction

4
doc/pair_lubricate.txt
View File

@ -40,6 +40,10 @@ their relative velocities in the presence of a background solvent with
viscosity mu. Note that this is dynamic viscosity which has units of
mass/distance/time, not kinematic viscosity.
The viscosity mu can be varied in a time-dependent manner over the
course of a simluation. See the "fix adapt"_fix_adapt.html command
for details.
Rc is the outer cutoff specified in the pair_style command, the
translational velocities of the 2 particles are v1 and v2, the angular
velocities are w1 and w2, and n is the unit vector in the direction

50
doc/pair_soft.html
View File

@ -20,8 +20,8 @@
<P><B>Examples:</B>
</P>
<PRE>pair_style soft 2.5
pair_coeff * * 0.0 60.0
pair_coeff 1 1 0.0 60.0 3.0
pair_coeff * * 10.0
pair_coeff 1 1 10.0 3.0
</PRE>
<P><B>Description:</B>
</P>
@ -30,11 +30,12 @@ pair_coeff 1 1 0.0 60.0 3.0
<CENTER><IMG SRC = "Eqs/pair_soft.jpg">
</CENTER>
<P>It is useful for pushing apart overlapping atoms, since it does not
blow up as r goes to 0. A is a pre-factor that varies in time from
the start to the end of the run. The <A HREF = "run.html">run</A> command documents
how to make the ramping take place across multiple runs. Rc is the
cutoff. See the <A HREF = "fix_nve_limit.html">fix nve/limit</A> command for
another way to push apart overlapping atoms.
blow up as r goes to 0. A is a pre-factor that can be made to vary in
time from the start to the end of the run (see discussion below),
e.g. to start with a very soft potential and slowly harden the
interactions over time. Rc is the cutoff. See the <A HREF = "fix_nve_limit.html">fix
nve/limit</A> command for another way to push apart
overlapping atoms.
</P>
<P>The following coefficients must be defined for each pair of atom types
via the <A HREF = "pair_coeff.html">pair_coeff</A> command as in the examples above,
@ -42,28 +43,34 @@ or in the data file or restart files read by the
<A HREF = "read_data.html">read_data</A> or <A HREF = "read_restart.html">read_restart</A>
commands, or by mixing as described below:
</P>
<UL><LI>Astart (energy units)
<LI>Astop (energy units)
<UL><LI>A (energy units)
<LI>cutoff (distance units)
</UL>
<P>Astart and Astop are the values of the prefactor at the start and end
of the next run. At intermediate times the value of A will be ramped
between these 2 values. Note that before performing a 2nd run, you
will want to adjust the values of Astart and Astop for all type pairs,
or switch to a new pair style.
</P>
<P>The last coefficient is optional. If not specified, the global soft
cutoff is used.
</P>
<P>The <A HREF = "fix_adapt.html">fix adapt</A> command can be used to vary A over the
course of a simulation. For example these commands will vary the
prefactor A for all pairwise interactions from 0.0 at the beginning to
30.0 at the end of a 10,000 step run:
</P>
<PRE>variable prefactor equal 30.0*elapsed/10000
fix 1 all adapt 1 pair soft a * * prefactor
</PRE>
<P>Note that a formula defined by an <A HREF = "variable.html">equal-style variable</A>
can use the current timestep, elapsed time in the current run, elapsed
time since the beginning of a series of runs, as well as access other
variables.
</P>
<HR>
<P><B>Mixing, shift, table, tail correction, restart, rRESPA info</B>:
</P>
<P>For atom type pairs I,J and I != J, the Astart, Astop coefficients and
cutoff distance for this pair style can be mixed. Astart and Atop are
always mixed via a <I>geometric</I> rule. The cutoff is mixed according to
the pair_modify mix value. The default mix value is <I>geometric</I>. See
the "pair_modify" command for details.
<P>For atom type pairs I,J and I != J, the A coefficient and cutoff
distance for this pair style can be mixed. A is always mixed via a
<I>geometric</I> rule. The cutoff is mixed according to the pair_modify
mix value. The default mix value is <I>geometric</I>. See the
"pair_modify" command for details.
</P>
<P>This pair style does not support the <A HREF = "pair_modify.html">pair_modify</A>
shift option, since the pair interaction goes to 0.0 at the cutoff.
@ -85,7 +92,8 @@ to be specified in an input script that reads a restart file.
</P>
<P><B>Related commands:</B>
</P>
<P><A HREF = "pair_coeff.html">pair_coeff</A>, <A HREF = "fix_nve_limit.html">fix nve/limit</A>
<P><A HREF = "pair_coeff.html">pair_coeff</A>, <A HREF = "fix_nve_limit.html">fix nve/limit</A>, <A HREF = "fix_adapt.html">fix
adapt</A>
</P>
<P><B>Default:</B> none
</P>

50
doc/pair_soft.txt
View File

@ -17,8 +17,8 @@ cutoff = global cutoff for soft interactions (distance units) :ul
[Examples:]
pair_style soft 2.5
pair_coeff * * 0.0 60.0
pair_coeff 1 1 0.0 60.0 3.0 :pre
pair_coeff * * 10.0
pair_coeff 1 1 10.0 3.0 :pre
[Description:]
@ -27,11 +27,12 @@ Style {soft} computes pairwise interactions with the formula
:c,image(Eqs/pair_soft.jpg)
It is useful for pushing apart overlapping atoms, since it does not
blow up as r goes to 0. A is a pre-factor that varies in time from
the start to the end of the run. The "run"_run.html command documents
how to make the ramping take place across multiple runs. Rc is the
cutoff. See the "fix nve/limit"_fix_nve_limit.html command for
another way to push apart overlapping atoms.
blow up as r goes to 0. A is a pre-factor that can be made to vary in
time from the start to the end of the run (see discussion below),
e.g. to start with a very soft potential and slowly harden the
interactions over time. Rc is the cutoff. See the "fix
nve/limit"_fix_nve_limit.html command for another way to push apart
overlapping atoms.
The following coefficients must be defined for each pair of atom types
via the "pair_coeff"_pair_coeff.html command as in the examples above,
@ -39,28 +40,34 @@ or in the data file or restart files read by the
"read_data"_read_data.html or "read_restart"_read_restart.html
commands, or by mixing as described below:
Astart (energy units)
Astop (energy units)
A (energy units)
cutoff (distance units) :ul
Astart and Astop are the values of the prefactor at the start and end
of the next run. At intermediate times the value of A will be ramped
between these 2 values. Note that before performing a 2nd run, you
will want to adjust the values of Astart and Astop for all type pairs,
or switch to a new pair style.
The last coefficient is optional. If not specified, the global soft
cutoff is used.
The "fix adapt"_fix_adapt.html command can be used to vary A over the
course of a simulation. For example these commands will vary the
prefactor A for all pairwise interactions from 0.0 at the beginning to
30.0 at the end of a 10,000 step run:
variable prefactor equal 30.0*elapsed/10000
fix 1 all adapt 1 pair soft a * * prefactor :pre
Note that a formula defined by an "equal-style variable"_variable.html
can use the current timestep, elapsed time in the current run, elapsed
time since the beginning of a series of runs, as well as access other
variables.
:line
[Mixing, shift, table, tail correction, restart, rRESPA info]:
For atom type pairs I,J and I != J, the Astart, Astop coefficients and
cutoff distance for this pair style can be mixed. Astart and Atop are
always mixed via a {geometric} rule. The cutoff is mixed according to
the pair_modify mix value. The default mix value is {geometric}. See
the "pair_modify" command for details.
For atom type pairs I,J and I != J, the A coefficient and cutoff
distance for this pair style can be mixed. A is always mixed via a
{geometric} rule. The cutoff is mixed according to the pair_modify
mix value. The default mix value is {geometric}. See the
"pair_modify" command for details.
This pair style does not support the "pair_modify"_pair_modify.html
shift option, since the pair interaction goes to 0.0 at the cutoff.
@ -82,6 +89,7 @@ This pair style can only be used via the {pair} keyword of the
[Related commands:]
"pair_coeff"_pair_coeff.html, "fix nve/limit"_fix_nve_limit.html
"pair_coeff"_pair_coeff.html, "fix nve/limit"_fix_nve_limit.html, "fix
adapt"_fix_adapt.html
[Default:] none

3
doc/run.html
View File

@ -69,8 +69,7 @@ performed and you want a <A HREF = "fix.html">fix</A> command that changes some
value over time (e.g. temperature) to make the change across the
entire set of runs and not just a single run. See the doc page for
individual fixes to see which ones can be used with the <I>start/stop</I>
keywords. The <A HREF = "pair_style.html">pair_style soft</A> potential also
changes its pair potential coefficients in this manner.
keywords.
</P>
<P>For example, consider this fix followed by 10 run commands:
</P>

3
doc/run.txt
View File

@ -62,8 +62,7 @@ performed and you want a "fix"_fix.html command that changes some
value over time (e.g. temperature) to make the change across the
entire set of runs and not just a single run. See the doc page for
individual fixes to see which ones can be used with the {start/stop}
keywords. The "pair_style soft"_pair_style.html potential also
changes its pair potential coefficients in this manner.
keywords.
For example, consider this fix followed by 10 run commands:

9
examples/micelle/in.micelle
View File

@ -13,20 +13,25 @@ read_data data.micelle
special_bonds fene
pair_style soft 1.12246
pair_coeff * * 1.0 20.0 1.12246
pair_coeff * * 0.0 1.12246
bond_style harmonic
bond_coeff 1 50.0 0.75
velocity all create 0.45 2349852
variable prefactor equal 1.0+elapsed*(20.0-1.0)/1000
fix 1 all nve
fix 2 all temp/rescale 100 0.45 0.45 0.02 1.0
fix 3 all enforce2d
fix 3 all adapt 1 pair soft a * * prefactor
fix 4 all enforce2d
thermo 50
run 1000
unfix 3
# Main run
pair_style lj/cut 2.5

2
src/COLLOID/pair_colloid.cpp
View File

@ -471,7 +471,7 @@ double PairColloid::single(int i, int j, int itype, int jtype, double rsq,
c2 = a2[itype][jtype];
K[1] = c2*c2;
K[2] = rsq;
K[0] = K[1] - rsq;
K[0] = K[1] - rsq;
K[4] = rsq*rsq;
K[3] = K[1] - K[2];
K[3] *= K[3]*K[3];

24
src/COLLOID/pair_lubricate.cpp
View File

@ -533,3 +533,27 @@ void PairLubricate::read_restart_settings(FILE *fp)
delete random;
random = new RanMars(lmp,seed + comm->me);
}
/* ----------------------------------------------------------------------
check if name is recognized, return integer index for that name
if name not recognized, return -1
if type pair setting, return -2 if no type pairs are set
------------------------------------------------------------------------- */
int PairLubricate::pre_adapt(char *name, int ilo, int ihi, int jlo, int jhi)
{
if (strcmp(name,"mu") == 0) return 0;
return -1;
}
/* ----------------------------------------------------------------------
adapt parameter indexed by which
change all pair variables affected by the reset parameter
if type pair setting, set I-J and J-I coeffs
------------------------------------------------------------------------- */
void PairLubricate::adapt(int which, int ilo, int ihi, int jlo, int jhi,
double value)
{
mu = value;
}

2
src/COLLOID/pair_lubricate.h
View File

@ -37,6 +37,8 @@ class PairLubricate : public Pair {
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
int pre_adapt(char *, int, int, int, int);
void adapt(int, int, int, int, int, double);
protected:
double cut_inner_global,cut_global;

2
src/MOLECULE/bond_table.cpp
View File

@ -375,7 +375,7 @@ void BondTable::compute_table(Table *tb)
// delta = table spacing for N-1 bins
int tlm1 = tablength-1;
tb->delta = (tb->hi - tb->lo)/ nm1;
tb->delta = (tb->hi - tb->lo)/ tlm1;
tb->invdelta = 1.0/tb->delta;
tb->deltasq6 = tb->delta*tb->delta / 6.0;

2
src/Makefile
View File

@ -14,7 +14,7 @@ OBJ = $(SRC:.cpp=.o)
# Package variables
PACKAGE = asphere class2 colloid dipole dsmc gpu granular \
kspace manybody meam molecule opt peri poems prd reax xtc
kspace manybody meam molecule opt peri poems prd reax shock xtc
PACKUSER = user-ackland user-atc user-cd-eam user-cg-cmm user-ewaldn \
user-imd user-smd

69
src/compute_dihedral_local.cpp
View File

@ -109,9 +109,8 @@ int ComputeDihedralLocal::compute_dihedrals(int flag)
{
int i,m,n,atom1,atom2,atom3,atom4;
double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm;
double sb1,sb2,sb3,rb1,rb3,c0,b1mag2,b1mag,b2mag2;
double b2mag,b3mag2,b3mag,ctmp,r12c1,c1mag,r12c2;
double c2mag,sin2,sc1,sc2,s1,s2,s12,c;
double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,rginv,ra2inv,rb2inv,rabinv;
double s,c;
double *pbuf;
double **x = atom->x;
@ -148,71 +147,53 @@ int ComputeDihedralLocal::compute_dihedrals(int flag)
if (flag) {
// phi calculation from dihedral style OPLS
// phi calculation from dihedral style harmonic
if (pflag >= 0) {
vb1x = x[atom1][0] - x[atom2][0];
vb1y = x[atom1][1] - x[atom2][1];
vb1z = x[atom1][2] - x[atom2][2];
domain->minimum_image(vb1x,vb1y,vb1z);
vb2x = x[atom3][0] - x[atom2][0];
vb2y = x[atom3][1] - x[atom2][1];
vb2z = x[atom3][2] - x[atom2][2];
domain->minimum_image(vb2x,vb2y,vb2z);
vb2xm = -vb2x;
vb2ym = -vb2y;
vb2zm = -vb2z;
domain->minimum_image(vb2xm,vb2ym,vb2zm);
vb3x = x[atom4][0] - x[atom3][0];
vb3y = x[atom4][1] - x[atom3][1];
vb3z = x[atom4][2] - x[atom3][2];
domain->minimum_image(vb3x,vb3y,vb3z);
sb1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z);
sb2 = 1.0 / (vb2x*vb2x + vb2y*vb2y + vb2z*vb2z);
sb3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z);
rb1 = sqrt(sb1);
rb3 = sqrt(sb3);
c0 = (vb1x*vb3x + vb1y*vb3y + vb1z*vb3z) * rb1*rb3;
b1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z;
b1mag = sqrt(b1mag2);
b2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z;
b2mag = sqrt(b2mag2);
b3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z;
b3mag = sqrt(b3mag2);
ctmp = vb1x*vb2x + vb1y*vb2y + vb1z*vb2z;
r12c1 = 1.0 / (b1mag*b2mag);
c1mag = ctmp * r12c1;
ax = vb1y*vb2zm - vb1z*vb2ym;
ay = vb1z*vb2xm - vb1x*vb2zm;
az = vb1x*vb2ym - vb1y*vb2xm;
bx = vb3y*vb2zm - vb3z*vb2ym;
by = vb3z*vb2xm - vb3x*vb2zm;
bz = vb3x*vb2ym - vb3y*vb2xm;
ctmp = vb2xm*vb3x + vb2ym*vb3y + vb2zm*vb3z;
r12c2 = 1.0 / (b2mag*b3mag);
c2mag = ctmp * r12c2;
rasq = ax*ax + ay*ay + az*az;
rbsq = bx*bx + by*by + bz*bz;
rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm;
rg = sqrt(rgsq);
rginv = ra2inv = rb2inv = 0.0;
if (rg > 0) rginv = 1.0/rg;
if (rasq > 0) ra2inv = 1.0/rasq;
if (rbsq > 0) rb2inv = 1.0/rbsq;
rabinv = sqrt(ra2inv*rb2inv);
sin2 = MAX(1.0 - c1mag*c1mag,0.0);
sc1 = sqrt(sin2);
if (sc1 < SMALL) sc1 = SMALL;
sc1 = 1.0/sc1;
sin2 = MAX(1.0 - c2mag*c2mag,0.0);
sc2 = sqrt(sin2);
if (sc2 < SMALL) sc2 = SMALL;
sc2 = 1.0/sc2;
s1 = sc1 * sc1;
s2 = sc2 * sc2;
s12 = sc1 * sc2;
c = (c0 + c1mag*c2mag) * s12;
c = (ax*bx + ay*by + az*bz)*rabinv;
s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z);
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;
pbuf[n] = 180.0*acos(c)/PI;
pbuf[n] = 180.0*atan2(s,c)/PI;
}
n += nvalues;
}

1
src/compute_erotate_sphere.cpp
View File

@ -109,6 +109,7 @@ double ComputeERotateSphere::compute_scalar()
radius[i]*radius[i]*mass[itype];
}
}
} else {
if (rmass) {
for (i = 0; i < nlocal; i++)

43
src/compute_property_local.cpp
View File

@ -26,7 +26,7 @@
using namespace LAMMPS_NS;
enum{NONE,PAIR,BOND,ANGLE,DIHEDRAL,IMPROPER};
enum{NONE,NEIGH,PAIR,BOND,ANGLE,DIHEDRAL,IMPROPER};
#define DELTA 10000
@ -50,7 +50,18 @@ ComputePropertyLocal::ComputePropertyLocal(LAMMPS *lmp, int narg, char **arg) :
for (int iarg = 3; iarg < narg; iarg++) {
i = iarg-3;
if (strcmp(arg[iarg],"patom1") == 0) {
if (strcmp(arg[iarg],"natom1") == 0) {
pack_choice[i] = &ComputePropertyLocal::pack_patom1;
if (kindflag != NONE && kindflag != NEIGH)
error->all("Compute property/local cannot use these inputs together");
kindflag = NEIGH;
} else if (strcmp(arg[iarg],"natom2") == 0) {
pack_choice[i] = &ComputePropertyLocal::pack_patom2;
if (kindflag != NONE && kindflag != NEIGH)
error->all("Compute property/local cannot use these inputs together");
kindflag = NEIGH;
} else if (strcmp(arg[iarg],"patom1") == 0) {
pack_choice[i] = &ComputePropertyLocal::pack_patom1;
if (kindflag != NONE && kindflag != PAIR)
error->all("Compute property/local cannot use these inputs together");
@ -184,16 +195,16 @@ ComputePropertyLocal::~ComputePropertyLocal()
void ComputePropertyLocal::init()
{
if (kindflag == PAIR) {
if (kindflag == NEIGH || kindflag == PAIR) {
if (force->pair == NULL)
error->all("No pair style is defined for compute property/local");
if (force->pair->single_enable == 0)
error->all("Pair style does not support compute property/local");
}
// for PAIR< need an occasional half neighbor list
// for NEIGH/PAIR need an occasional half neighbor list
if (kindflag == PAIR) {
if (kindflag == NEIGH || kindflag == PAIR) {
int irequest = neighbor->request((void *) this);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->compute = 1;
@ -201,9 +212,10 @@ void ComputePropertyLocal::init()
}
// do initial memory allocation so that memory_usage() is correct
// cannot be done yet for PAIR, since neigh list does not exist
// cannot be done yet for NEIGH/PAIR, since neigh list does not exist
if (kindflag == PAIR) ncount = 0;
if (kindflag == NEIGH) ncount = 0;
else if (kindflag == PAIR) ncount = 0;
else if (kindflag == BOND) ncount = count_bonds(0);
else if (kindflag == ANGLE) ncount = count_angles(0);
else if (kindflag == DIHEDRAL) ncount = count_dihedrals(0);
@ -228,7 +240,8 @@ void ComputePropertyLocal::compute_local()
// count local entries and generate list of indices
if (kindflag == PAIR) ncount = count_pairs(0);
if (kindflag == NEIGH) ncount = count_pairs(0,0);
else if (kindflag == PAIR) ncount = count_pairs(0,1);
else if (kindflag == BOND) ncount = count_bonds(0);
else if (kindflag == ANGLE) ncount = count_angles(0);
else if (kindflag == DIHEDRAL) ncount = count_dihedrals(0);
@ -237,7 +250,8 @@ void ComputePropertyLocal::compute_local()
if (ncount > nmax) reallocate(ncount);
size_local_rows = ncount;
if (kindflag == PAIR) ncount = count_pairs(1);
if (kindflag == NEIGH) ncount = count_pairs(1,0);
else if (kindflag == PAIR) ncount = count_pairs(1,1);
else if (kindflag == BOND) ncount = count_bonds(1);
else if (kindflag == ANGLE) ncount = count_angles(1);
else if (kindflag == DIHEDRAL) ncount = count_dihedrals(1);
@ -259,10 +273,11 @@ void ComputePropertyLocal::compute_local()
count pairs and compute pair info on this proc
only count pair once if newton_pair is off
both atom I,J must be in group
if flag is set, compute requested info about pair
if allflag is set, compute requested info about pair
if forceflag = 1, pair must be within force cutoff, else neighbor cutoff
------------------------------------------------------------------------- */
int ComputePropertyLocal::count_pairs(int flag)
int ComputePropertyLocal::count_pairs(int allflag, int forceflag)
{
int i,j,m,n,ii,jj,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz;
@ -281,7 +296,7 @@ int ComputePropertyLocal::count_pairs(int flag)
// invoke half neighbor list (will copy or build if necessary)
if (flag == 0) neighbor->build_one(list->index);
if (allflag == 0) neighbor->build_one(list->index);
inum = list->inum;
ilist = list->ilist;
@ -323,9 +338,9 @@ int ComputePropertyLocal::count_pairs(int flag)
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
jtype = type[j];
if (rsq >= cutsq[itype][jtype]) continue;
if (forceflag && rsq >= cutsq[itype][jtype]) continue;
if (flag) {
if (allflag) {
indices[m][0] = tag[i];
indices[m][1] = tag[j];
}

2
src/compute_property_local.h
View File

@ -46,7 +46,7 @@ class ComputePropertyLocal : public Compute {
int ncount;
int **indices;
int count_pairs(int);
int count_pairs(int, int);
int count_bonds(int);
int count_angles(int);
int count_dihedrals(int);

90
src/compute_reduce.cpp
View File

@ -488,8 +488,8 @@ double ComputeReduce::compute_one(int m, int flag)
// only include atoms in group for atom properties and per-atom quantities
index = -1;
int v_idx = value2index[m];
int a_idx = argindex[m];
int vidx = value2index[m];
int aidx = argindex[m];
int *mask = atom->mask;
int nlocal = atom->nlocal;
@ -504,28 +504,28 @@ double ComputeReduce::compute_one(int m, int flag)
double **x = atom->x;
if (flag < 0) {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) combine(one,x[i][a_idx],i);
} else one = x[flag][a_idx];
if (mask[i] & groupbit) combine(one,x[i][aidx],i);
} else one = x[flag][aidx];
} else if (which[m] == V) {
double **v = atom->v;
if (flag < 0) {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) combine(one,v[i][a_idx],i);
} else one = v[flag][a_idx];
if (mask[i] & groupbit) combine(one,v[i][aidx],i);
} else one = v[flag][aidx];
} else if (which[m] == F) {
double **f = atom->f;
if (flag < 0) {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) combine(one,f[i][a_idx],i);
} else one = f[flag][a_idx];
if (mask[i] & groupbit) combine(one,f[i][aidx],i);
} else one = f[flag][aidx];
// invoke compute if not previously invoked
} else if (which[m] == COMPUTE) {
Compute *compute = modify->compute[v_idx];
Compute *compute = modify->compute[vidx];
if (flavor[m] == GLOBAL) {
if (a_idx == 0) {
if (aidx == 0) {
if (!(compute->invoked_flag & INVOKED_VECTOR)) {
compute->compute_vector();
compute->invoked_flag |= INVOKED_VECTOR;
@ -541,13 +541,13 @@ double ComputeReduce::compute_one(int m, int flag)
compute->compute_array();
compute->invoked_flag |= INVOKED_ARRAY;
}
double **comp_arr = compute->array;
double **carray = compute->array;
int n = compute->size_array_rows;
int a_idx_1 = a_idx - 1;
int aidx_1 = aidx - 1;
if (flag < 0)
for (i = 0; i < n; i++)
combine(one,comp_arr[i][a_idx_1],i);
else one = comp_arr[flag][a_idx_1];
combine(one,carray[i][aidx_1],i);
else one = carray[flag][aidx_1];
}
} else if (flavor[m] == PERATOM) {
@ -556,7 +556,7 @@ double ComputeReduce::compute_one(int m, int flag)
compute->invoked_flag |= INVOKED_PERATOM;
}
if (a_idx == 0) {
if (aidx == 0) {
double *comp_vec = compute->vector_atom;
int n = nlocal;
if (flag < 0) {
@ -564,13 +564,13 @@ double ComputeReduce::compute_one(int m, int flag)
if (mask[i] & groupbit) combine(one,comp_vec[i],i);
} else one = comp_vec[flag];
} else {
double **comp_arr = compute->array_atom;
double **carray_atom = compute->array_atom;
int n = nlocal;
int a_idxm1 = a_idx - 1;
int aidxm1 = aidx - 1;
if (flag < 0) {
for (i = 0; i < n; i++)
if (mask[i] & groupbit) combine(one,comp_arr[i][a_idxm1],i);
} else one = comp_arr[flag][a_idxm1];
if (mask[i] & groupbit) combine(one,carray_atom[i][aidxm1],i);
} else one = carray_atom[flag][aidxm1];
}
} else if (flavor[m] == LOCAL) {
@ -579,7 +579,7 @@ double ComputeReduce::compute_one(int m, int flag)
compute->invoked_flag |= INVOKED_LOCAL;
}
if (a_idx == 0) {
if (aidx == 0) {
double *comp_vec = compute->vector_local;
int n = compute->size_local_rows;
if (flag < 0)
@ -587,40 +587,40 @@ double ComputeReduce::compute_one(int m, int flag)
combine(one,comp_vec[i],i);
else one = comp_vec[flag];
} else {
double **comp_arr = compute->array_local;
double **carray_local = compute->array_local;
int n = compute->size_local_rows;
int a_idxm1 = a_idx - 1;
int aidxm1 = aidx - 1;
if (flag < 0)
for (i = 0; i < n; i++)
combine(one,comp_arr[i][a_idxm1],i);
else one = comp_arr[flag][a_idxm1];
combine(one,carray_local[i][aidxm1],i);
else one = carray_local[flag][aidxm1];
}
}
// access fix fields, check if fix frequency is a match
} else if (which[m] == FIX) {
if (update->ntimestep % modify->fix[v_idx]->peratom_freq)
if (update->ntimestep % modify->fix[vidx]->peratom_freq)
error->all("Fix used in compute reduce not computed at compatible time");
Fix *fix = modify->fix[v_idx];
Fix *fix = modify->fix[vidx];
if (flavor[m] == GLOBAL) {
if (a_idx == 0) {
if (aidx == 0) {
int n = fix->size_vector;
if (flag < 0)
for (i = 0; i < n; i++)
combine(one,fix->compute_vector(i),i);
else one = fix->compute_vector(flag);
} else {
int a_idxm1 = a_idx - 1;
int aidxm1 = aidx - 1;
if (flag < 0)
for (i = 0; i < nlocal; i++)
combine(one,fix->compute_array(i,a_idxm1),i);
else one = fix->compute_array(flag,a_idxm1);
combine(one,fix->compute_array(i,aidxm1),i);
else one = fix->compute_array(flag,aidxm1);
}
} else if (flavor[m] == PERATOM) {
if (a_idx == 0) {
if (aidx == 0) {
double *fix_vector = fix->vector_atom;
int n = nlocal;
if (flag < 0) {
@ -629,15 +629,15 @@ double ComputeReduce::compute_one(int m, int flag)
} else one = fix_vector[flag];
} else {
double **fix_array = fix->array_atom;
int a_idxm1 = a_idx - 1;
int aidxm1 = aidx - 1;
if (flag < 0) {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) combine(one,fix_array[i][a_idxm1],i);
} else one = fix_array[flag][a_idxm1];
if (mask[i] & groupbit) combine(one,fix_array[i][aidxm1],i);
} else one = fix_array[flag][aidxm1];
}
} else if (flavor[m] == LOCAL) {
if (a_idx == 0) {
if (aidx == 0) {
double *fix_vector = fix->vector_local;
int n = fix->size_local_rows;
if (flag < 0)
@ -647,11 +647,11 @@ double ComputeReduce::compute_one(int m, int flag)
} else {
double **fix_array = fix->array_local;
int n = fix->size_local_rows;
int a_idxm1 = a_idx - 1;
int aidxm1 = aidx - 1;
if (flag < 0)
for (i = 0; i < n; i++)
combine(one,fix_array[i][a_idxm1],i);
else one = fix_array[flag][a_idxm1];
combine(one,fix_array[i][aidxm1],i);
else one = fix_array[flag][aidxm1];
}
}
@ -665,7 +665,7 @@ double ComputeReduce::compute_one(int m, int flag)
memory->smalloc(maxatom*sizeof(double),"reduce:varatom");
}
input->variable->compute_atom(v_idx,igroup,varatom,1,0);
input->variable->compute_atom(vidx,igroup,varatom,1,0);
if (flag < 0) {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) combine(one,varatom[i],i);
@ -679,15 +679,15 @@ double ComputeReduce::compute_one(int m, int flag)
double ComputeReduce::count(int m)
{
int v_idx = value2index[m];
int a_idx = argindex[m];
int vidx = value2index[m];
int aidx = argindex[m];
if (which[m] == X || which[m] == V || which[m] == F)
return group->count(igroup);
else if (which[m] == COMPUTE) {
Compute *compute = modify->compute[v_idx];
Compute *compute = modify->compute[vidx];
if (flavor[m] == GLOBAL) {
if (a_idx == 0) return(1.0*compute->size_vector);
if (aidx == 0) return(1.0*compute->size_vector);
else return(1.0*compute->size_array_rows);
} else if (flavor[m] == PERATOM) {
return group->count(igroup);
@ -698,9 +698,9 @@ double ComputeReduce::count(int m)
return ncountall;
}
} else if (which[m] == FIX) {
Fix *fix = modify->fix[v_idx];
Fix *fix = modify->fix[vidx];
if (flavor[m] == GLOBAL) {
if (a_idx == 0) return(1.0*fix->size_vector);
if (aidx == 0) return(1.0*fix->size_vector);
else return(1.0*fix->size_array_rows);
} else if (flavor[m] == PERATOM) {
return group->count(igroup);

212
src/fix_adapt.cpp Normal file
View File

@ -0,0 +1,212 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "math.h"
#include "string.h"
#include "stdlib.h"
#include "fix_adapt.h"
#include "atom.h"
#include "force.h"
#include "pair.h"
#include "input.h"
#include "variable.h"
#include "error.h"
using namespace LAMMPS_NS;
enum{PAIR,ATOM};
enum{DIAMETER};
/* ---------------------------------------------------------------------- */
FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
{
if (narg < 4) error->all("Illegal fix adapt command");
nevery = atoi(arg[3]);
if (nevery <= 0) error->all("Illegal fix adapt command");
// count # of adaptations
nadapt = 0;
int iarg = 4;
while (iarg < narg) {
if (strcmp(arg[iarg],"pair") == 0) {
if (iarg+6 > narg) error->all("Illegal fix adapt command");
nadapt++;
iarg += 6;
} else if (strcmp(arg[iarg],"atom") == 0) {
if (iarg+3 > narg) error->all("Illegal fix adapt command");
nadapt++;
iarg += 3;
} else error->all("Illegal fix adapt command");
}
// allocate per-adapt vectors
which = new int[nadapt];
pair = new char*[nadapt];
param = new char*[nadapt];
ilo = new int[nadapt];
ihi = new int[nadapt];
jlo = new int[nadapt];
jhi = new int[nadapt];
var = new char*[nadapt];
ivar = new int[nadapt];
pairptr = new Pair*[nadapt];
pairindex = new int[nadapt];
awhich = new int[nadapt];
// parse keywords
nadapt = 0;
iarg = 4;
while (iarg < narg) {
if (strcmp(arg[iarg],"pair") == 0) {
if (iarg+6 > narg) error->all("Illegal fix adapt command");
which[nadapt] = PAIR;
int n = strlen(arg[iarg+1]) + 1;
pair[nadapt] = new char[n];
strcpy(pair[nadapt],arg[iarg+1]);
n = strlen(arg[iarg+2]) + 1;
param[nadapt] = new char[n];
strcpy(param[nadapt],arg[iarg+2]);
force->bounds(arg[iarg+3],atom->ntypes,ilo[nadapt],ihi[nadapt]);
force->bounds(arg[iarg+4],atom->ntypes,jlo[nadapt],jhi[nadapt]);
n = strlen(arg[iarg+5]) + 1;
var[nadapt] = new char[n];
strcpy(var[nadapt],arg[iarg+5]);
nadapt++;
iarg += 6;
} else if (strcmp(arg[iarg],"atom") == 0) {
if (iarg+3 > narg) error->all("Illegal fix adapt command");
which[nadapt] = ATOM;
int n = strlen(arg[iarg+1]) + 1;
param[nadapt] = new char[n];
strcpy(param[nadapt],arg[iarg+1]);
n = strlen(arg[iarg+2]) + 1;
var[nadapt] = new char[n];
strcpy(var[nadapt],arg[iarg+2]);
nadapt++;
iarg += 3;
} else error->all("Illegal fix adapt command");
}
}
/* ---------------------------------------------------------------------- */
FixAdapt::~FixAdapt()
{
for (int i = 0; i < nadapt; i++) {
if (which[i] == PAIR) delete [] pair[i];
delete [] param[i];
delete [] var[i];
}
delete [] which;
delete [] pair;
delete [] param;
delete [] ilo;
delete [] ihi;
delete [] jlo;
delete [] jhi;
delete [] var;
delete [] ivar;
delete [] pairptr;
delete [] pairindex;
delete [] awhich;
}
/* ---------------------------------------------------------------------- */
int FixAdapt::setmask()
{
int mask = 0;
mask |= PRE_FORCE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixAdapt::init()
{
// error checks
for (int m = 0; m < nadapt; m++) {
if (which[m] == PAIR) {
pairptr[m] = force->pair_match(pair[m],1);
if (pairptr[m] == NULL)
error->all("Fix adapt pair style does not exist");
pairindex[m] =
pairptr[m]->pre_adapt(param[m],ilo[m],ihi[m],jlo[m],jhi[m]);
if (pairindex[m] == -1)
error->all("Fix adapt pair parameter is not recognized");
if (pairindex[m] == -2)
error->all("Fix adapt pair types are not valid");
} else if (which[m] == ATOM) {
if (strcmp(param[m],"diameter") == 0) {
awhich[m] = DIAMETER;
if (!atom->radius_flag)
error->all("Fix adapt requires atom attribute radius");
} else error->all("Fix adapt atom attribute is not recognized");
}
ivar[m] = input->variable->find(var[m]);
if (ivar[m] < 0) error->all("Variable name for fix adapt does not exist");
if (!input->variable->equalstyle(ivar[m]))
error->all("Variable for fix adapt is not equal style");
}
// set params to values for initial force calculation
// needs to happen here in init() instead of setup()
// because modify->setup() is called after pre-Verlet forces are computed
pre_force(0);
}
/* ---------------------------------------------------------------------- */
void FixAdapt::pre_force(int vflag)
{
for (int m = 0; m < nadapt; m++) {
double value = input->variable->compute_equal(ivar[m]);
if (which[m] == PAIR)
pairptr[m]->adapt(pairindex[m],ilo[m],ihi[m],jlo[m],jhi[m],value);
else if (which[m] == ATOM) {
// set radius from diameter
// set rmass if both rmass and density are defined
if (awhich[m] == DIAMETER) {
int mflag = 0;
if (atom->rmass_flag && atom->density_flag) mflag = 1;
double PI = 4.0*atan(1.0);
double *radius = atom->radius;
double *rmass = atom->rmass;
double *density = atom->density;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
radius[i] = 0.5*value;
rmass[i] = 4.0*PI/3.0 * radius[i]*radius[i]*radius[i] * density[i];
}
}
}
}
}

50
src/fix_adapt.h Normal file
View File

@ -0,0 +1,50 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(adapt,FixAdapt)
#else
#ifndef LMP_FIX_ADAPT_H
#define LMP_FIX_ADAPT_H
#include "fix.h"
namespace LAMMPS_NS {
class FixAdapt : public Fix {
public:
FixAdapt(class LAMMPS *, int, char **);
~FixAdapt();
int setmask();
void init();
void pre_force(int);
private:
int nadapt;
int *which;
char **pair,**param,**var;
int *ilo,*ihi,*jlo,*jhi;
int *ivar;
class Pair **pairptr;
int *pairindex;
int *awhich;
};
}
#endif
#endif

2
src/fix_aveforce.cpp
View File

@ -40,6 +40,8 @@ FixAveForce::FixAveForce(LAMMPS *lmp, int narg, char **arg) :
global_freq = 1;
extvector = 1;
xstr = ystr = zstr = NULL;
if (strstr(arg[3],"v_") == arg[3]) {
int n = strlen(&arg[3][2]) + 1;
xstr = new char[n];

8
src/fix_wall_reflect.cpp
View File

@ -74,7 +74,7 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) :
if (iarg+2 > narg) error->all("Illegal fix wall/reflect command");
if (strcmp(arg[iarg+1],"EDGE") == 0) {
yloflag = EDGE;
ylo = domain->boxlo[0];
ylo = domain->boxlo[1];
} else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
yloflag = VARIABLE;
int n = strlen(&arg[iarg+1][2]) + 1;
@ -89,7 +89,7 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) :
if (iarg+2 > narg) error->all("Illegal fix wall/reflect command");
if (strcmp(arg[iarg+1],"EDGE") == 0) {
yhiflag = EDGE;
yhi = domain->boxhi[0];
yhi = domain->boxhi[1];
} else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
yhiflag = VARIABLE;
int n = strlen(&arg[iarg+1][2]) + 1;
@ -104,7 +104,7 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) :
if (iarg+2 > narg) error->all("Illegal fix wall/reflect command");
if (strcmp(arg[iarg+1],"EDGE") == 0) {
zloflag = EDGE;
zlo = domain->boxlo[0];
zlo = domain->boxlo[2];
} else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
zloflag = VARIABLE;
int n = strlen(&arg[iarg+1][2]) + 1;
@ -119,7 +119,7 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) :
if (iarg+2 > narg) error->all("Illegal fix wall/reflect command");
if (strcmp(arg[iarg+1],"EDGE") == 0) {
zhiflag = EDGE;
zhi = domain->boxhi[0];
zhi = domain->boxhi[2];
} else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
zhiflag = VARIABLE;
int n = strlen(&arg[iarg+1][2]) + 1;

8
src/pair.cpp
View File

@ -23,7 +23,6 @@
#include "stdlib.h"
#include "string.h"
#include "pair.h"
#include "pair_soft.h"
#include "atom.h"
#include "neighbor.h"
#include "neigh_list.h"
@ -895,13 +894,6 @@ void Pair::write_file(int narg, char **arg)
force->init();
// if pair style = soft, set prefactor to final value
Pair *spair = force->pair_match("soft",1);
if (spair)
((PairSoft *) spair)->prefactor[itype][jtype] =
((PairSoft *) spair)->prestop[itype][jtype];
// if pair style = any of EAM, swap in dummy fp vector
double eamfp[2];

2
src/pair.h
View File

@ -97,6 +97,8 @@ class Pair : protected Pointers {
virtual void swap_eam(double *, double **) {}
virtual void reset_dt() {}
virtual void min_pointers(double **, double **) {}
virtual int pre_adapt(char *, int, int, int, int) {return -1;}
virtual void adapt(int, int, int, int, int, double) {}
protected:
int allocated; // 0/1 = whether arrays are allocated

2
src/pair_hybrid.h
View File

@ -32,7 +32,7 @@ class PairHybrid : public Pair {
char **keywords; // style name of each Pair style
PairHybrid(class LAMMPS *);
~PairHybrid();
virtual ~PairHybrid();
void compute(int, int);
void settings(int, char **);
virtual void coeff(int, char **);

1
src/pair_hybrid_overlay.h
View File

@ -27,6 +27,7 @@ namespace LAMMPS_NS {
class PairHybridOverlay : public PairHybrid {
public:
PairHybridOverlay(class LAMMPS *);
~PairHybridOverlay() {}
void coeff(int, char **);
private:

87
src/pair_soft.cpp
View File

@ -14,6 +14,7 @@
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "pair_soft.h"
#include "atom.h"
#include "comm.h"
@ -43,8 +44,6 @@ PairSoft::~PairSoft()
memory->destroy_2d_int_array(setflag);
memory->destroy_2d_double_array(cutsq);
memory->destroy_2d_double_array(prestart);
memory->destroy_2d_double_array(prestop);
memory->destroy_2d_double_array(prefactor);
memory->destroy_2d_double_array(cut);
}
@ -63,21 +62,6 @@ void PairSoft::compute(int eflag, int vflag)
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
// set current prefactor
// for minimization, set to prestop
// for dynamics, ramp from prestart to prestop
// for 0-step dynamics, set to prestart
double delta = update->ntimestep - update->beginstep;
if (update->whichflag == 2) delta = 1.0;
else if (update->nsteps) delta /= update->endstep - update->beginstep;
else delta = 0.0;
int ntypes = atom->ntypes;
for (i = 1; i <= ntypes; i++)
for (j = 1; j <= ntypes; j++)
prefactor[i][j] = prestart[i][j] +
delta * (prestop[i][j] - prestart[i][j]);
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
@ -161,20 +145,8 @@ void PairSoft::allocate()
cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq");
prestart = memory->create_2d_double_array(n+1,n+1,"pair:prestart");
prestop = memory->create_2d_double_array(n+1,n+1,"pair:prestop");
prefactor = memory->create_2d_double_array(n+1,n+1,"pair:prefactor");
cut = memory->create_2d_double_array(n+1,n+1,"pair:cut");
// init prestart and prestop to 0.0
// since pair_hybrid can use all types even if pair_soft sub-class
// never sets them
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++) {
prestart[i][j] = 0.0;
prestop[i][j] = 0.0;
}
}
/* ----------------------------------------------------------------------
@ -203,24 +175,22 @@ void PairSoft::settings(int narg, char **arg)
void PairSoft::coeff(int narg, char **arg)
{
if (narg < 4 || narg > 5) error->all("Incorrect args for pair coefficients");
if (narg < 3 || narg > 4) error->all("Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
force->bounds(arg[0],atom->ntypes,ilo,ihi);
force->bounds(arg[1],atom->ntypes,jlo,jhi);
double prestart_one = force->numeric(arg[2]);
double prestop_one = force->numeric(arg[3]);
double prefactor_one = force->numeric(arg[2]);
double cut_one = cut_global;
if (narg == 5) cut_one = force->numeric(arg[4]);
if (narg == 4) cut_one = force->numeric(arg[3]);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
prestart[i][j] = prestart_one;
prestop[i][j] = prestop_one;
prefactor[i][j] = prefactor_one;
cut[i][j] = cut_one;
setflag[i][j] = 1;
count++;
@ -239,13 +209,11 @@ double PairSoft::init_one(int i, int j)
// always mix prefactors geometrically
if (setflag[i][j] == 0) {
prestart[i][j] = sqrt(prestart[i][i]*prestart[j][j]);
prestop[i][j] = sqrt(prestop[i][i]*prestop[j][j]);
prefactor[i][j] = sqrt(prefactor[i][i]*prefactor[j][j]);
cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
}
prestart[j][i] = prestart[i][j];
prestop[j][i] = prestop[i][j];
prefactor[j][i] = prefactor[i][j];
cut[j][i] = cut[i][j];
return cut[i][j];
@ -264,8 +232,7 @@ void PairSoft::write_restart(FILE *fp)
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp);
if (setflag[i][j]) {
fwrite(&prestart[i][j],sizeof(double),1,fp);
fwrite(&prestop[i][j],sizeof(double),1,fp);
fwrite(&prefactor[i][j],sizeof(double),1,fp);
fwrite(&cut[i][j],sizeof(double),1,fp);
}
}
@ -289,12 +256,10 @@ void PairSoft::read_restart(FILE *fp)
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
if (setflag[i][j]) {
if (me == 0) {
fread(&prestart[i][j],sizeof(double),1,fp);
fread(&prestop[i][j],sizeof(double),1,fp);
fread(&prefactor[i][j],sizeof(double),1,fp);
fread(&cut[i][j],sizeof(double),1,fp);
}
MPI_Bcast(&prestart[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&prestop[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&prefactor[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
}
}
@ -324,6 +289,38 @@ void PairSoft::read_restart_settings(FILE *fp)
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
}
/* ----------------------------------------------------------------------
check if name is recognized, return integer index for that name
if name not recognized, return -1
if type pair setting, return -2 if no type pairs are set
------------------------------------------------------------------------- */
int PairSoft::pre_adapt(char *name, int ilo, int ihi, int jlo, int jhi)
{
int count = 0;
for (int i = ilo; i <= ihi; i++)
for (int j = MAX(jlo,i); j <= jhi; j++)
count++;
if (count == 0) return -2;
if (strcmp(name,"a") == 0) return 0;
return -1;
}
/* ----------------------------------------------------------------------
adapt parameter indexed by which
change all pair variables affected by the reset parameter
if type pair setting, set I-J and J-I coeffs
------------------------------------------------------------------------- */
void PairSoft::adapt(int which, int ilo, int ihi, int jlo, int jhi,
double value)
{
for (int i = ilo; i <= ihi; i++)
for (int j = MAX(jlo,i); j <= jhi; j++)
prefactor[i][j] = prefactor[j][i] = value;
}
/* ---------------------------------------------------------------------- */
double PairSoft::single(int i, int j, int itype, int jtype, double rsq,

3
src/pair_soft.h
View File

@ -38,12 +38,13 @@ class PairSoft : public Pair {
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
int pre_adapt(char *, int, int, int, int);
void adapt(int, int, int, int, int, double);
double single(int, int, int, int, double, double, double, double &);
private:
double PI;
double cut_global;
double **prestart,**prestop;
double **prefactor;
double **cut;

4
src/pair_table.cpp
View File

@ -287,15 +287,13 @@ void PairTable::coeff(int narg, char **arg)
// match = 1 if don't need to spline read-in tables
// this is only the case if r values needed by final tables
// exactly match r values read from file
// for tabstyle SPLINE, always need to build spline tables
tb->match = 0;
if (tabstyle == LINEAR && tb->ninput == tablength &&
tb->rflag == RSQ && tb->rhi == tb->cut) tb->match = 1;
if (tabstyle == SPLINE && tb->ninput == tablength &&
tb->rflag == RSQ && tb->rhi == tb->cut) tb->match = 1;
if (tabstyle == BITMAP && tb->ninput == 1 << tablength &&
tb->rflag == BMP && tb->rhi == tb->cut) tb->match = 1;
if (tb->rflag == BMP && tb->match == 0)
error->all("Bitmapped table in file does not match requested table");

4
src/region_cone.cpp
View File

@ -436,7 +436,7 @@ int RegCone::surface_exterior(double *x, double cutoff)
// y is exterior to cone or on its surface
// corner1234 = 4 corner pts of half trapezoid = cone surf in plane of y
// project y to 3 line segments in half trapezoid (4th is axis of cone)
// project x to 3 line segments in half trapezoid (4th is axis of cone)
// nearest = point on surface of cone that y is closest to
// could be edge of cone
// do not add contact point if r >= cutoff
@ -481,7 +481,7 @@ int RegCone::surface_exterior(double *x, double cutoff)
// z is exterior to cone or on its surface
// corner1234 = 4 corner pts of half trapezoid = cone surf in plane of z
// project z to 3 line segments in half trapezoid (4th is axis of cone)
// project x to 3 line segments in half trapezoid (4th is axis of cone)
// nearest = point on surface of cone that z is closest to
// could be edge of cone
// do not add contact point if r >= cutoff

4
src/region_cylinder.cpp
View File

@ -352,7 +352,7 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
if (r < radius && x[1] > lo && x[1] < hi) return 0;
// y is exterior to cylinder or on its surface
// xp,yp,zp = point on surface of cylinder that y is closest to
// xp,yp,zp = point on surface of cylinder that x is closest to
// could be edge of cylinder
// do not add contact point if r >= cutoff
@ -383,7 +383,7 @@ int RegCylinder::surface_exterior(double *x, double cutoff)
if (r < radius && x[2] > lo && x[2] < hi) return 0;
// z is exterior to cylinder or on its surface
// xp,yp,zp = point on surface of cylinder that z is closest to
// xp,yp,zp = point on surface of cylinder that x is closest to
// could be edge of cylinder
// do not add contact point if r >= cutoff

2
src/version.h
View File

@ -1 +1 @@
#define LAMMPS_VERSION "19 Jun 2010"
#define LAMMPS_VERSION "24 Jun 2010"