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

Merge pull request #510 from akohlmey/collected-small-changes

Browse Source 
Collected small changes
This commit is contained in:
sjplimp
2017-06-15 11:12:14 -06:00
committed by GitHub
parent 1a29048940 b62d526cc9
commit d3c90f3c14
81 changed files with 633 additions and 1318 deletions
Download Patch File Download Diff File Expand all files Collapse all files

BIN
doc/src/Eqs/pair_lj_sf.jpg
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 15 KiB

11
doc/src/Eqs/pair_lj_sf.tex
View File

@ -1,11 +0,0 @@
\documentclass[12pt]{article}
\begin{document}
\begin{eqnarray*}
F & = & F_{\mathrm{LJ}}(r) - F_{\mathrm{LJ}}(r_{\mathrm{c}}) \qquad r < r_{\mathrm{c}} \\
E & = & E_{\mathrm{LJ}}(r) - E_{\mathrm{LJ}}(r_{\mathrm{c}}) + (r - r_{\mathrm{c}}) F_{\mathrm{LJ}}(r_{\mathrm{c}}) \qquad r < r_{\mathrm{c}} \\
\mathrm{with} \qquad E_{\mathrm{LJ}}(r) & = & 4 \epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right] \qquad \mathrm{and} \qquad F_{\mathrm{LJ}}(r) = - E^\prime_{\mathrm{LJ}}(r)
\end{eqnarray*}
\end{document}

1
doc/src/Section_commands.txt
View File

@ -1039,7 +1039,6 @@ package"_Section_start.html#start_3.
"lj/sdk (gko)"_pair_sdk.html, "lj/sdk (gko)"_pair_sdk.html,
"lj/sdk/coul/long (go)"_pair_sdk.html, "lj/sdk/coul/long (go)"_pair_sdk.html,
"lj/sdk/coul/msm (o)"_pair_sdk.html, "lj/sdk/coul/msm (o)"_pair_sdk.html,
"lj/sf (o)"_pair_lj_sf.html,
"meam/spline (o)"_pair_meam_spline.html, "meam/spline (o)"_pair_meam_spline.html,
"meam/sw/spline"_pair_meam_sw_spline.html, "meam/sw/spline"_pair_meam_sw_spline.html,
"mgpt"_pair_mgpt.html, "mgpt"_pair_mgpt.html,

4
doc/src/Section_packages.txt
View File

@ -2027,8 +2027,8 @@ algorithm to formulate single-particle constraint functions
g(xi,yi,zi) = 0 and their derivative (i.e. the normal of the manifold) g(xi,yi,zi) = 0 and their derivative (i.e. the normal of the manifold)
n = grad(g). n = grad(g).
[Author:] Stefan Paquay (Eindhoven University of Technology (TU/e), The [Author:] Stefan Paquay (until 2017: Eindhoven University of Technology (TU/e), The
Netherlands) Netherlands; since 2017: Brandeis University, Waltham, MA, USA)
[Install or un-install:] [Install or un-install:]

2
doc/src/fix_adapt.txt
View File

@ -47,7 +47,7 @@ keyword = {scale} or {reset} :l
fix 1 all adapt 1 pair soft a 1 1 v_prefactor fix 1 all adapt 1 pair soft a 1 1 v_prefactor
fix 1 all adapt 1 pair soft a 2* 3 v_prefactor fix 1 all adapt 1 pair soft a 2* 3 v_prefactor
fix 1 all adapt 1 pair lj/cut epsilon * * v_scale1 coul/cut scale 3 3 v_scale2 scale yes reset yes fix 1 all adapt 1 pair lj/cut epsilon * * v_scale1 coul/cut scale 3 3 v_scale2 scale yes reset yes
fix 1 all adapt 10 atom diameter v_size fix 1 all adapt 10 atom diameter v_size :pre
variable ramp_up equal "ramp(0.01,0.5)" variable ramp_up equal "ramp(0.01,0.5)"
fix stretch all adapt 1 bond harmonic r0 1 v_ramp_up :pre fix stretch all adapt 1 bond harmonic r0 1 v_ramp_up :pre

6
doc/src/fix_deform.txt
View File

@ -565,8 +565,10 @@ more instructions on how to use the accelerated styles effectively.
[Restart, fix_modify, output, run start/stop, minimize info:] [Restart, fix_modify, output, run start/stop, minimize info:]
No information about this fix is written to "binary restart This fix will restore the initial box settings from "binary restart
files"_restart.html. None of the "fix_modify"_fix_modify.html options files"_restart.html, which allows the fix to be properly continue
deformation, when using the start/stop options of the "run"_run.html
command. None of the "fix_modify"_fix_modify.html options
are relevant to this fix. No global or per-atom quantities are stored are relevant to this fix. No global or per-atom quantities are stored
by this fix for access by various "output by this fix for access by various "output
commands"_Section_howto.html#howto_15. commands"_Section_howto.html#howto_15.

1
doc/src/lammps.book
View File

@ -447,7 +447,6 @@ pair_lj96.html
pair_lj_cubic.html pair_lj_cubic.html
pair_lj_expand.html pair_lj_expand.html
pair_lj_long.html pair_lj_long.html
pair_lj_sf.html
pair_lj_smooth.html pair_lj_smooth.html
pair_lj_smooth_linear.html pair_lj_smooth_linear.html
pair_lj_soft.html pair_lj_soft.html

114
doc/src/pair_lj_sf.txt
View File

@ -1,114 +0,0 @@
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
pair_style lj/sf command :h3
pair_style lj/sf/omp command :h3
[Syntax:]
pair_style lj/sf cutoff :pre
cutoff = global cutoff for Lennard-Jones interactions (distance units) :ul
[Examples:]
pair_style lj/sf 2.5
pair_coeff * * 1.0 1.0
pair_coeff 1 1 1.0 1.0 3.0 :pre
[Description:]
Style {lj/sf} computes a truncated and force-shifted LJ interaction
(Shifted Force Lennard-Jones), so that both the potential and the
force go continuously to zero at the cutoff "(Toxvaerd)"_#Toxvaerd:
:c,image(Eqs/pair_lj_sf.jpg)
The following coefficients must be defined for each pair of atoms
types via the "pair_coeff"_pair_coeff.html command as in the examples
above, 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:
epsilon (energy units)
sigma (distance units)
cutoff (distance units) :ul
The last coefficient is optional. If not specified, the global
LJ cutoff specified in the pair_style command is used.
:line
Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are
functionally the same as the corresponding style without the suffix.
They have been optimized to run faster, depending on your available
hardware, as discussed in "Section 5"_Section_accelerate.html
of the manual. The accelerated styles take the same arguments and
should produce the same results, except for round-off and precision
issues.
These accelerated styles are part of the GPU, USER-INTEL, KOKKOS,
USER-OMP and OPT packages, respectively. They are only enabled if
LAMMPS was built with those packages. See the "Making
LAMMPS"_Section_start.html#start_3 section for more info.
You can specify the accelerated styles explicitly in your input script
by including their suffix, or you can use the "-suffix command-line
switch"_Section_start.html#start_7 when you invoke LAMMPS, or you can
use the "suffix"_suffix.html command in your input script.
See "Section 5"_Section_accelerate.html of the manual for
more instructions on how to use the accelerated styles effectively.
:line
[Mixing, shift, table, tail correction, restart, rRESPA info]:
For atom type pairs I,J and I != J, the epsilon and sigma
coefficients and cutoff distance for this pair style can be mixed.
Rin is a cutoff value and is mixed like the cutoff. The
default mix value is {geometric}. See the "pair_modify" command for
details.
The "pair_modify"_pair_modify.html shift option is not relevant for
this pair style, since the pair interaction goes to 0.0 at the cutoff.
The "pair_modify"_pair_modify.html table option is not relevant
for this pair style.
This pair style does not support the "pair_modify"_pair_modify.html
tail option for adding long-range tail corrections to energy and
pressure, since the energy of the pair interaction is smoothed to 0.0
at the cutoff.
This pair style writes its information to "binary restart
files"_restart.html, so pair_style and pair_coeff commands do not need
to be specified in an input script that reads a restart file.
This pair style can only be used via the {pair} keyword of the
"run_style respa"_run_style.html command. It does not support the
{inner}, {middle}, {outer} keywords.
:line
[Restrictions:]
This pair style is part of the USER-MISC package. It is only enabled
if LAMMPS was built with that package. See the "Making
LAMMPS"_Section_start.html#start_3 section for more info.
[Related commands:]
"pair_coeff"_pair_coeff.html
[Default:] none
:line
:link(Toxvaerd)
[(Toxvaerd)] Toxvaerd, Dyre, J Chem Phys, 134, 081102 (2011).

31
doc/src/pair_lj_smooth_linear.txt
View File

@ -11,26 +11,26 @@ pair_style lj/smooth/linear/omp command :h3
[Syntax:] [Syntax:]
pair_style lj/smooth/linear Rc :pre pair_style lj/smooth/linear cutoff :pre
Rc = cutoff for lj/smooth/linear interactions (distance units) :ul cutoff = global cutoff for Lennard-Jones interactions (distance units) :ul
[Examples:] [Examples:]
pair_style lj/smooth/linear 5.456108274435118 pair_style lj/smooth/linear 2.5
pair_coeff * * 0.7242785984051078 2.598146797350056 pair_coeff * * 1.0 1.0
pair_coeff 1 1 20.0 1.3 9.0 :pre pair_coeff 1 1 0.3 3.0 9.0 :pre
[Description:] [Description:]
Style {lj/smooth/linear} computes a LJ interaction that combines the Style {lj/smooth/linear} computes a truncated and force-shifted LJ
standard 12/6 Lennard-Jones function and subtracts a linear term that interaction (aka Shifted Force Lennard-Jones) that combines the
includes the cutoff distance Rc, as in this formula: standard 12/6 Lennard-Jones function and subtracts a linear term based
on the cutoff distance, so that both, the potential and the force, go
continuously to zero at the cutoff Rc "(Toxvaerd)"_#Toxvaerd:
:c,image(Eqs/pair_lj_smooth_linear.jpg) :c,image(Eqs/pair_lj_smooth_linear.jpg)
At the cutoff Rc, the energy and force (its 1st derivative) will be 0.0.
The following coefficients must be defined for each pair of atoms The following coefficients must be defined for each pair of atoms
types via the "pair_coeff"_pair_coeff.html command as in the examples types via the "pair_coeff"_pair_coeff.html command as in the examples
above, or in the data file or restart files read by the above, or in the data file or restart files read by the
@ -41,8 +41,8 @@ epsilon (energy units)
sigma (distance units) sigma (distance units)
cutoff (distance units) :ul cutoff (distance units) :ul
The last coefficient is optional. If not specified, the global value The last coefficient is optional. If not specified, the global
for Rc is used. LJ cutoff specified in the pair_style command is used.
:line :line
@ -76,10 +76,11 @@ and cutoff distance can be mixed. The default mix value is geometric.
See the "pair_modify" command for details. See the "pair_modify" command for details.
This pair style does not support the "pair_modify"_pair_modify.html This pair style does not support the "pair_modify"_pair_modify.html
shift option for the energy of the pair interaction. shift option for the energy of the pair interaction, since it goes
to 0.0 at the cutoff by construction.
The "pair_modify"_pair_modify.html table option is not relevant for The "pair_modify"_pair_modify.html table option is not relevant
this pair style. for this pair style.
This pair style does not support the "pair_modify"_pair_modify.html This pair style does not support the "pair_modify"_pair_modify.html
tail option for adding long-range tail corrections to energy and tail option for adding long-range tail corrections to energy and

1
doc/src/pairs.txt
View File

@ -49,7 +49,6 @@ Pair Styles :h1
pair_lj_cubic pair_lj_cubic
pair_lj_expand pair_lj_expand
pair_lj_long pair_lj_long
pair_lj_sf
pair_lj_smooth pair_lj_smooth
pair_lj_smooth_linear pair_lj_smooth_linear
pair_lj_soft pair_lj_soft

6
doc/src/set.txt
View File

@ -80,6 +80,7 @@ keyword = {type} or {type/fraction} or {mol} or {x} or {y} or {z} or \
value can be an atom-style variable (see below) value can be an atom-style variable (see below)
{image} nx ny nz {image} nx ny nz
nx,ny,nz = which periodic image of the simulation box the atom is in nx,ny,nz = which periodic image of the simulation box the atom is in
any of nx,ny,nz can be an atom-style variable (see below)
{bond} value = bond type for all bonds between selected atoms {bond} value = bond type for all bonds between selected atoms
{angle} value = angle type for all angles between selected atoms {angle} value = angle type for all angles between selected atoms
{dihedral} value = dihedral type for all dihedrals between selected atoms {dihedral} value = dihedral type for all dihedrals between selected atoms
@ -363,9 +364,8 @@ A value of -1 means subtract 1 box length to get the true value.
LAMMPS updates these flags as atoms cross periodic boundaries during LAMMPS updates these flags as atoms cross periodic boundaries during
the simulation. The flags can be output with atom snapshots via the the simulation. The flags can be output with atom snapshots via the
"dump"_dump.html command. If a value of NULL is specified for any of "dump"_dump.html command. If a value of NULL is specified for any of
nx,ny,nz, then the current image value for that dimension is nx,ny,nz, then the current image value for that dimension is unchanged.
unchanged. For non-periodic dimensions only a value of 0 can be For non-periodic dimensions only a value of 0 can be specified.
specified. This keyword does not allow use of atom-style variables.
This command can be useful after a system has been equilibrated and This command can be useful after a system has been equilibrated and
atoms have diffused one or more box lengths in various directions. atoms have diffused one or more box lengths in various directions.
This command can then reset the image values for atoms so that they This command can then reset the image values for atoms so that they

8
doc/src/special_bonds.txt
View File

@ -65,7 +65,13 @@ sense to define permanent bonds between atoms that interact via these
potentials, though such bonds may exist elsewhere in your system, potentials, though such bonds may exist elsewhere in your system,
e.g. when using the "pair_style hybrid"_pair_hybrid.html command. e.g. when using the "pair_style hybrid"_pair_hybrid.html command.
Thus LAMMPS ignores special_bonds settings when manybody potentials Thus LAMMPS ignores special_bonds settings when manybody potentials
are calculated. are calculated. Please note, that the existence of explicit bonds
for atoms that are described by a manybody potential will alter the
neigborlist and thus can render the computation of those interactions
invalid, since those pairs are not only used to determine direct
pairwise interactions but also neighbors of neighbors and more.
The recommended course of action is to remove such bonds, or - if
that is not possible - use a special bonds setting of 1.0 1.0 1.0.
NOTE: Unlike some commands in LAMMPS, you cannot use this command NOTE: Unlike some commands in LAMMPS, you cannot use this command
multiple times in an incremental fashion: e.g. to first set the LJ multiple times in an incremental fashion: e.g. to first set the LJ

12
lib/meam/meam_dens_final.F
View File

@ -98,13 +98,14 @@ c Complete the calculation of density
call G_gam(Gamma(i),ibar_meam(elti), call G_gam(Gamma(i),ibar_meam(elti),
$ gsmooth_factor,G,errorflag) $ gsmooth_factor,G,errorflag)
if (errorflag.ne.0) return if (errorflag.ne.0) return
call get_shpfcn(shp,lattce_meam(elti,elti))
if (ibar_meam(elti).le.0) then if (ibar_meam(elti).le.0) then
Gbar = 1.d0 Gbar = 1.d0
dGbar = 0.d0
else else
call get_shpfcn(shp,lattce_meam(elti,elti))
if (mix_ref_t.eq.1) then if (mix_ref_t.eq.1) then
gam = (t_ave(1,i)*shpi(1)+t_ave(2,i)*shpi(2) gam = (t_ave(1,i)*shp(1)+t_ave(2,i)*shp(2)
$ +t_ave(3,i)*shpi(3))/(Z*Z) $ +t_ave(3,i)*shp(3))/(Z*Z)
else else
gam = (t1_meam(elti)*shp(1)+t2_meam(elti)*shp(2) gam = (t1_meam(elti)*shp(1)+t2_meam(elti)*shp(2)
$ +t3_meam(elti)*shp(3))/(Z*Z) $ +t3_meam(elti)*shp(3))/(Z*Z)
@ -119,9 +120,8 @@ c Complete the calculation of density
Gbar = 1.d0 Gbar = 1.d0
dGbar = 0.d0 dGbar = 0.d0
else else
call get_shpfcn(shpi,lattce_meam(elti,elti)) gam = (t_ave(1,i)*shp(1)+t_ave(2,i)*shp(2)
gam = (t_ave(1,i)*shpi(1)+t_ave(2,i)*shpi(2) $ +t_ave(3,i)*shp(3))/(Z*Z)
$ +t_ave(3,i)*shpi(3))/(Z*Z)
call dG_gam(gam,ibar_meam(elti),gsmooth_factor, call dG_gam(gam,ibar_meam(elti),gsmooth_factor,
$ Gbar,dGbar) $ Gbar,dGbar)
endif endif

2
src/.gitignore vendored
View File

@ -740,8 +740,6 @@
/pair_lj_sdk_coul_long.h /pair_lj_sdk_coul_long.h
/pair_lj_sdk_coul_msm.cpp /pair_lj_sdk_coul_msm.cpp
/pair_lj_sdk_coul_msm.h /pair_lj_sdk_coul_msm.h
/pair_lj_sf.cpp
/pair_lj_sf.h
/pair_lj_sf_dipole_sf.cpp /pair_lj_sf_dipole_sf.cpp
/pair_lj_sf_dipole_sf.h /pair_lj_sf_dipole_sf.h
/pair_lubricateU.cpp /pair_lubricateU.cpp

18
src/MC/fix_gcmc.cpp
View File

@ -501,7 +501,7 @@ void FixGCMC::init()
if (ifix < 0) error->all(FLERR,"Fix gcmc rigid fix does not exist"); if (ifix < 0) error->all(FLERR,"Fix gcmc rigid fix does not exist");
fixrigid = modify->fix[ifix]; fixrigid = modify->fix[ifix];
int tmp; int tmp;
if (onemols != (Molecule **) fixrigid->extract("onemol",tmp)) if (&onemols[imol] != (Molecule **) fixrigid->extract("onemol",tmp))
error->all(FLERR, error->all(FLERR,
"Fix gcmc and fix rigid/small not using " "Fix gcmc and fix rigid/small not using "
"same molecule template ID"); "same molecule template ID");
@ -516,7 +516,7 @@ void FixGCMC::init()
if (ifix < 0) error->all(FLERR,"Fix gcmc shake fix does not exist"); if (ifix < 0) error->all(FLERR,"Fix gcmc shake fix does not exist");
fixshake = modify->fix[ifix]; fixshake = modify->fix[ifix];
int tmp; int tmp;
if (onemols != (Molecule **) fixshake->extract("onemol",tmp)) if (&onemols[imol] != (Molecule **) fixshake->extract("onemol",tmp))
error->all(FLERR,"Fix gcmc and fix shake not using " error->all(FLERR,"Fix gcmc and fix shake not using "
"same molecule template ID"); "same molecule template ID");
} }
@ -1398,11 +1398,12 @@ void FixGCMC::attempt_molecule_insertion()
// FixRigidSmall::set_molecule stores rigid body attributes // FixRigidSmall::set_molecule stores rigid body attributes
// FixShake::set_molecule stores shake info for molecule // FixShake::set_molecule stores shake info for molecule
for (int submol = 0; submol < nmol; ++submol) {
if (rigidflag) if (rigidflag)
fixrigid->set_molecule(nlocalprev,maxtag_all,imol,com_coord,vnew,quat); fixrigid->set_molecule(nlocalprev,maxtag_all,submol,com_coord,vnew,quat);
else if (shakeflag) else if (shakeflag)
fixshake->set_molecule(nlocalprev,maxtag_all,imol,com_coord,vnew,quat); fixshake->set_molecule(nlocalprev,maxtag_all,submol,com_coord,vnew,quat);
}
atom->natoms += natoms_per_molecule; atom->natoms += natoms_per_molecule;
if (atom->natoms < 0) if (atom->natoms < 0)
error->all(FLERR,"Too many total atoms"); error->all(FLERR,"Too many total atoms");
@ -2058,11 +2059,12 @@ void FixGCMC::attempt_molecule_insertion_full()
// FixRigidSmall::set_molecule stores rigid body attributes // FixRigidSmall::set_molecule stores rigid body attributes
// FixShake::set_molecule stores shake info for molecule // FixShake::set_molecule stores shake info for molecule
for (int submol = 0; submol < nmol; ++submol) {
if (rigidflag) if (rigidflag)
fixrigid->set_molecule(nlocalprev,maxtag_all,imol,com_coord,vnew,quat); fixrigid->set_molecule(nlocalprev,maxtag_all,submol,com_coord,vnew,quat);
else if (shakeflag) else if (shakeflag)
fixshake->set_molecule(nlocalprev,maxtag_all,imol,com_coord,vnew,quat); fixshake->set_molecule(nlocalprev,maxtag_all,submol,com_coord,vnew,quat);
}
atom->natoms += natoms_per_molecule; atom->natoms += natoms_per_molecule;
if (atom->natoms < 0) if (atom->natoms < 0)
error->all(FLERR,"Too many total atoms"); error->all(FLERR,"Too many total atoms");

5
src/Purge.list
View File

@ -16,6 +16,11 @@ style_region.h
style_neigh_bin.h style_neigh_bin.h
style_neigh_pair.h style_neigh_pair.h
style_neigh_stencil.h style_neigh_stencil.h
# deleted on 6 June 2017
pair_lj_sf.cpp
pair_lj_sf.h
pair_lj_sf_omp.cpp
pair_lj_sf_omp.h
# deleted on 4 May 2017 # deleted on 4 May 2017
dump_custom_vtk.cpp dump_custom_vtk.cpp
dump_custom_vtk.h dump_custom_vtk.h

2
src/RIGID/fix_rigid.cpp
View File

@ -50,7 +50,7 @@ enum{ISO,ANISO,TRICLINIC};
#define EPSILON 1.0e-7 #define EPSILON 1.0e-7
#define SINERTIA 0.4 // moment of inertia prefactor for sphere #define SINERTIA 0.4 // moment of inertia prefactor for sphere
#define EINERTIA 0.4 // moment of inertia prefactor for ellipsoid #define EINERTIA 0.2 // moment of inertia prefactor for ellipsoid
#define LINERTIA (1.0/12.0) // moment of inertia prefactor for line segment #define LINERTIA (1.0/12.0) // moment of inertia prefactor for line segment
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */

4
src/RIGID/fix_rigid_small.cpp
View File

@ -54,7 +54,7 @@ FixRigidSmall *FixRigidSmall::frsptr;
#define BIG 1.0e20 #define BIG 1.0e20
#define SINERTIA 0.4 // moment of inertia prefactor for sphere #define SINERTIA 0.4 // moment of inertia prefactor for sphere
#define EINERTIA 0.4 // moment of inertia prefactor for ellipsoid #define EINERTIA 0.2 // moment of inertia prefactor for ellipsoid
#define LINERTIA (1.0/12.0) // moment of inertia prefactor for line segment #define LINERTIA (1.0/12.0) // moment of inertia prefactor for line segment
#define DELTA_BODY 10000 #define DELTA_BODY 10000
@ -1047,7 +1047,7 @@ int FixRigidSmall::dof(int tgroup)
j = atom2body[i]; j = atom2body[i];
counts[j][2]++; counts[j][2]++;
if (mask[i] & tgroupbit) { if (mask[i] & tgroupbit) {
if (extended && eflags[i]) counts[j][1]++; if (extended && (eflags[i] & ~(POINT | DIPOLE))) counts[j][1]++;
else counts[j][0]++; else counts[j][0]++;
} }
} }

8
src/USER-MANIFOLD/README
View File

@ -7,10 +7,14 @@ box). It achieves this using the RATTLE constraint algorithm applied
to single-particle constraint functions g(xi,yi,zi) = 0 and their to single-particle constraint functions g(xi,yi,zi) = 0 and their
derivative (i.e. the normal of the manifold) n = grad(g). derivative (i.e. the normal of the manifold) n = grad(g).
Stefan Paquay, s.paquay@tue.nl Stefan Paquay, stefanpaquay@gmail.com
Applied Physics/Theory of Polymers and Soft Matter,
until 2017: Applied Physics/Theory of Polymers and Soft Matter,
Eindhoven University of Technology (TU/e), The Netherlands Eindhoven University of Technology (TU/e), The Netherlands
since 2017: Brandeis University, Waltham, MA, USA.
Thanks to Remy Kusters at TU/e for testing. Thanks to Remy Kusters at TU/e for testing.
This software is distributed under the GNU General Public License. This software is distributed under the GNU General Public License.

1
src/USER-MISC/README
View File

@ -64,7 +64,6 @@ pair_style gauss/cut, Axel Kohlmeyer, akohlmey at gmail.com, 1 Dec 11
pair_style lennard/mdf, Paolo Raiteri, p.raiteri at curtin.edu.au, 2 Dec 15 pair_style lennard/mdf, Paolo Raiteri, p.raiteri at curtin.edu.au, 2 Dec 15
pair_style list, Axel Kohlmeyer (Temple U), akohlmey at gmail.com, 1 Jun 13 pair_style list, Axel Kohlmeyer (Temple U), akohlmey at gmail.com, 1 Jun 13
pair_style lj/mdf, Paolo Raiteri, p.raiteri at curtin.edu.au, 2 Dec 15 pair_style lj/mdf, Paolo Raiteri, p.raiteri at curtin.edu.au, 2 Dec 15
pair_style lj/sf, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11
pair_style kolmogorov/crespi/z, Jaap Kroes (Radboud U), jaapkroes at gmail dot com, 28 Feb 17 pair_style kolmogorov/crespi/z, Jaap Kroes (Radboud U), jaapkroes at gmail dot com, 28 Feb 17
pair_style meam/spline, Alexander Stukowski (LLNL), alex at stukowski.com, 1 Feb 12 pair_style meam/spline, Alexander Stukowski (LLNL), alex at stukowski.com, 1 Feb 12
pair_style meam/sw/spline, Robert Rudd (LLNL), robert.rudd at llnl.gov, 1 Oct 12 pair_style meam/sw/spline, Robert Rudd (LLNL), robert.rudd at llnl.gov, 1 Oct 12

10
src/USER-MISC/pair_cdeam.cpp
View File

@ -456,11 +456,11 @@ void PairCDEAM::read_h_coeff(char *filename)
{ {
if(comm->me == 0) { if(comm->me == 0) {
// Open potential file // Open potential file
FILE *fp; FILE *fptr;
char line[MAXLINE]; char line[MAXLINE];
char nextline[MAXLINE]; char nextline[MAXLINE];
fp = force->open_potential(filename); fptr = force->open_potential(filename);
if (fp == NULL) { if (fptr == NULL) {
char str[128]; char str[128];
sprintf(str,"Cannot open EAM potential file %s", filename); sprintf(str,"Cannot open EAM potential file %s", filename);
error->one(FLERR,str); error->one(FLERR,str);
@ -468,7 +468,7 @@ void PairCDEAM::read_h_coeff(char *filename)
// h coefficients are stored at the end of the file. // h coefficients are stored at the end of the file.
// Skip to last line of file. // Skip to last line of file.
while(fgets(nextline, MAXLINE, fp) != NULL) { while(fgets(nextline, MAXLINE, fptr) != NULL) {
strcpy(line, nextline); strcpy(line, nextline);
} }
char* ptr = strtok(line, " \t\n\r\f"); char* ptr = strtok(line, " \t\n\r\f");
@ -483,7 +483,7 @@ void PairCDEAM::read_h_coeff(char *filename)
error->one(FLERR,"Failed to read h(x) function coefficients from EAM file."); error->one(FLERR,"Failed to read h(x) function coefficients from EAM file.");
// Close the potential file. // Close the potential file.
fclose(fp); fclose(fptr);
} }
MPI_Bcast(&nhcoeff, 1, MPI_INT, 0, world); MPI_Bcast(&nhcoeff, 1, MPI_INT, 0, world);

355
src/USER-MISC/pair_lj_sf.cpp
View File

@ -1,355 +0,0 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Laurent Joly (U Lyon), ljoly.ulyon@gmail.com
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "pair_lj_sf.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairLJShiftedForce::PairLJShiftedForce(LAMMPS *lmp) : Pair(lmp)
{
respa_enable = 0;
}
/* ---------------------------------------------------------------------- */
PairLJShiftedForce::~PairLJShiftedForce()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(cut);
memory->destroy(epsilon);
memory->destroy(sigma);
memory->destroy(lj1);
memory->destroy(lj2);
memory->destroy(lj3);
memory->destroy(lj4);
memory->destroy(foffset);
memory->destroy(offset);
}
}
/* ---------------------------------------------------------------------- */
void PairLJShiftedForce::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double rsq,r2inv,r6inv,forcelj,factor_lj;
double r,t;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
r2inv = 1.0/rsq;
r6inv = r2inv*r2inv*r2inv;
r = sqrt(rsq);
t = r/cut[itype][jtype];
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]) -
t*foffset[itype][jtype];
fpair = factor_lj*forcelj*r2inv;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
}
if (eflag) {
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) +
(t-1.0)*foffset[itype][jtype] - offset[itype][jtype];
evdwl *= factor_lj;
}
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairLJShiftedForce::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(cut,n+1,n+1,"pair:cut");
memory->create(epsilon,n+1,n+1,"pair:epsilon");
memory->create(sigma,n+1,n+1,"pair:sigma");
memory->create(lj1,n+1,n+1,"pair:lj1");
memory->create(lj2,n+1,n+1,"pair:lj2");
memory->create(lj3,n+1,n+1,"pair:lj3");
memory->create(lj4,n+1,n+1,"pair:lj4");
memory->create(foffset,n+1,n+1,"pair:foffset");
memory->create(offset,n+1,n+1,"pair:offset");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairLJShiftedForce::settings(int narg, char **arg)
{
if (narg != 1) error->all(FLERR,"Illegal pair_style command");
cut_global = force->numeric(FLERR,arg[0]);
if (cut_global <= 0.0)
error->all(FLERR,"Illegal pair_style command");
// reset cutoffs that have been explicitly set
if (allocated) {
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++)
if (setflag[i][j]) cut[i][j] = cut_global;
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairLJShiftedForce::coeff(int narg, char **arg)
{
if (narg < 4 || narg > 5)
error->all(FLERR,"Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
double epsilon_one = force->numeric(FLERR,arg[2]);
double sigma_one = force->numeric(FLERR,arg[3]);
double cut_one = cut_global;
if (narg == 5) cut_one = force->numeric(FLERR,arg[4]);
if (cut_one <= 0.0)
error->all(FLERR,"Incorrect args for pair coefficients");
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
epsilon[i][j] = epsilon_one;
sigma[i][j] = sigma_one;
cut[i][j] = cut_one;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairLJShiftedForce::init_one(int i, int j)
{
if (setflag[i][j] == 0) {
epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j],
sigma[i][i],sigma[j][j]);
sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]);
cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
}
lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
double ratio = sigma[i][j] / cut[i][j];
foffset[i][j] = 4.0 * epsilon[i][j] * (12.0 * pow(ratio,12.0) -
6.0 * pow(ratio,6.0));
offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0));
cut[j][i] = cut[i][j];
lj1[j][i] = lj1[i][j];
lj2[j][i] = lj2[i][j];
lj3[j][i] = lj3[i][j];
lj4[j][i] = lj4[i][j];
foffset[j][i] = foffset[i][j];
offset[j][i] = offset[i][j];
return cut[i][j];
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLJShiftedForce::write_restart(FILE *fp)
{
write_restart_settings(fp);
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp);
if (setflag[i][j]) {
fwrite(&epsilon[i][j],sizeof(double),1,fp);
fwrite(&sigma[i][j],sizeof(double),1,fp);
fwrite(&cut[i][j],sizeof(double),1,fp);
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairLJShiftedForce::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
int i,j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
if (setflag[i][j]) {
if (me == 0) {
fread(&epsilon[i][j],sizeof(double),1,fp);
fread(&sigma[i][j],sizeof(double),1,fp);
fread(&cut[i][j],sizeof(double),1,fp);
}
MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLJShiftedForce::write_restart_settings(FILE *fp)
{
fwrite(&cut_global,sizeof(double),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairLJShiftedForce::read_restart_settings(FILE *fp)
{
int me = comm->me;
if (me == 0) {
fread(&cut_global,sizeof(double),1,fp);
fread(&offset_flag,sizeof(int),1,fp);
fread(&mix_flag,sizeof(int),1,fp);
}
MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
}
/* ---------------------------------------------------------------------- */
double PairLJShiftedForce::single(int i, int j, int itype, int jtype, double rsq,
double factor_coul, double factor_lj,
double &fforce)
{
double r2inv,r6inv,forcelj,philj,r,t;
r2inv = 1.0/rsq;
r6inv = r2inv*r2inv*r2inv;
r = sqrt(rsq);
t = r/cut[itype][jtype];
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]) -
t*foffset[itype][jtype];
fforce = factor_lj*forcelj*r2inv;
philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) +
(t-1.0)*foffset[itype][jtype] - offset[itype][jtype];
return factor_lj*philj;
}

53
src/USER-MISC/pair_lj_sf.h
View File

@ -1,53 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
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 PAIR_CLASS
PairStyle(lj/sf,PairLJShiftedForce)
#else
#ifndef LMP_PAIR_LJ_SF_H
#define LMP_PAIR_LJ_SF_H
#include "pair.h"
namespace LAMMPS_NS {
class PairLJShiftedForce : public Pair {
public:
PairLJShiftedForce(class LAMMPS *);
virtual ~PairLJShiftedForce();
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
double single(int, int, int, int, double, double, double, double &);
protected:
double cut_global;
double **cut;
double **epsilon,**sigma;
double **lj1,**lj2,**lj3,**lj4,**foffset,**offset;
void allocate();
};
}
#endif
#endif

2
src/USER-OMP/fix_rigid_omp.cpp
View File

@ -39,7 +39,7 @@ using namespace MathConst;
enum{SINGLE,MOLECULE,GROUP}; // same as in FixRigid enum{SINGLE,MOLECULE,GROUP}; // same as in FixRigid
#define EINERTIA 0.4 // moment of inertia prefactor for ellipsoid #define EINERTIA 0.2 // moment of inertia prefactor for ellipsoid
typedef struct { double x,y,z; } dbl3_t; typedef struct { double x,y,z; } dbl3_t;

2
src/USER-OMP/fix_rigid_small_omp.cpp
View File

@ -37,7 +37,7 @@ using namespace LAMMPS_NS;
using namespace FixConst; using namespace FixConst;
using namespace MathConst; using namespace MathConst;
#define EINERTIA 0.4 // moment of inertia prefactor for ellipsoid #define EINERTIA 0.2 // moment of inertia prefactor for ellipsoid
enum{FULL_BODY,INITIAL,FINAL,FORCE_TORQUE,VCM_ANGMOM,XCM_MASS,ITENSOR,DOF}; enum{FULL_BODY,INITIAL,FINAL,FORCE_TORQUE,VCM_ANGMOM,XCM_MASS,ITENSOR,DOF};

164
src/USER-OMP/pair_lj_sf_omp.cpp
View File

@ -1,164 +0,0 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */
#include <math.h>
#include "pair_lj_sf_omp.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "suffix.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairLJShiftedForceOMP::PairLJShiftedForceOMP(LAMMPS *lmp) :
PairLJShiftedForce(lmp), ThrOMP(lmp, THR_PAIR)
{
suffix_flag |= Suffix::OMP;
respa_enable = 0;
}
/* ---------------------------------------------------------------------- */
void PairLJShiftedForceOMP::compute(int eflag, int vflag)
{
if (eflag || vflag) {
ev_setup(eflag,vflag);
} else evflag = vflag_fdotr = 0;
const int nall = atom->nlocal + atom->nghost;
const int nthreads = comm->nthreads;
const int inum = list->inum;
#if defined(_OPENMP)
#pragma omp parallel default(none) shared(eflag,vflag)
#endif
{
int ifrom, ito, tid;
loop_setup_thr(ifrom, ito, tid, inum, nthreads);
ThrData *thr = fix->get_thr(tid);
thr->timer(Timer::START);
ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr);
if (evflag) {
if (eflag) {
if (force->newton_pair) eval<1,1,1>(ifrom, ito, thr);
else eval<1,1,0>(ifrom, ito, thr);
} else {
if (force->newton_pair) eval<1,0,1>(ifrom, ito, thr);
else eval<1,0,0>(ifrom, ito, thr);
}
} else {
if (force->newton_pair) eval<0,0,1>(ifrom, ito, thr);
else eval<0,0,0>(ifrom, ito, thr);
}
thr->timer(Timer::PAIR);
reduce_thr(this, eflag, vflag, thr);
} // end of omp parallel region
}
template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
void PairLJShiftedForceOMP::eval(int iifrom, int iito, ThrData * const thr)
{
int i,j,ii,jj,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double t,rsq,r2inv,r6inv,forcelj,factor_lj;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = 0.0;
const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0];
const int * _noalias const type = atom->type;
const int nlocal = atom->nlocal;
const double * _noalias const special_lj = force->special_lj;
double fxtmp,fytmp,fztmp;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// loop over neighbors of my atoms
for (ii = iifrom; ii < iito; ++ii) {
i = ilist[ii];
xtmp = x[i].x;
ytmp = x[i].y;
ztmp = x[i].z;
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
fxtmp=fytmp=fztmp=0.0;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
j &= NEIGHMASK;
delx = xtmp - x[j].x;
dely = ytmp - x[j].y;
delz = ztmp - x[j].z;
rsq = delx*delx + dely*dely + delz*delz;
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
r2inv = 1.0/rsq;
r6inv = r2inv*r2inv*r2inv;
t = sqrt(rsq)/cut[itype][jtype];
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]) -
t*foffset[itype][jtype];
fpair = factor_lj*forcelj*r2inv;
fxtmp += delx*fpair;
fytmp += dely*fpair;
fztmp += delz*fpair;
if (NEWTON_PAIR || j < nlocal) {
f[j].x -= delx*fpair;
f[j].y -= dely*fpair;
f[j].z -= delz*fpair;
}
if (EFLAG) {
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) +
(t-1.0)*foffset[itype][jtype] - offset[itype][jtype];
evdwl *= factor_lj;
}
if (EVFLAG) ev_tally_thr(this,i,j,nlocal,NEWTON_PAIR,
evdwl,0.0,fpair,delx,dely,delz,thr);
}
}
f[i].x += fxtmp;
f[i].y += fytmp;
f[i].z += fztmp;
}
}
/* ---------------------------------------------------------------------- */
double PairLJShiftedForceOMP::memory_usage()
{
double bytes = memory_usage_thr();
bytes += PairLJShiftedForce::memory_usage();
return bytes;
}

48
src/USER-OMP/pair_lj_sf_omp.h
View File

@ -1,48 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(lj/sf/omp,PairLJShiftedForceOMP)
#else
#ifndef LMP_PAIR_LJ_SF_OMP_H
#define LMP_PAIR_LJ_SF_OMP_H
#include "pair_lj_sf.h"
#include "thr_omp.h"
namespace LAMMPS_NS {
class PairLJShiftedForceOMP : public PairLJShiftedForce, public ThrOMP {
public:
PairLJShiftedForceOMP(class LAMMPS *);
virtual void compute(int, int);
virtual double memory_usage();
private:
template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
void eval(int ifrom, int ito, ThrData * const thr);
};
}
#endif
#endif

1
src/USER-OMP/pair_lj_smooth_linear_omp.h
View File

@ -18,6 +18,7 @@
#ifdef PAIR_CLASS #ifdef PAIR_CLASS
PairStyle(lj/smooth/linear/omp,PairLJSmoothLinearOMP) PairStyle(lj/smooth/linear/omp,PairLJSmoothLinearOMP)
PairStyle(lj/sf/omp,PairLJSmoothLinearOMP)
#else #else

2
src/compute.h
View File

@ -153,7 +153,7 @@ class Compute : protected Pointers {
double **vbiasall; // stored velocity bias for all atoms double **vbiasall; // stored velocity bias for all atoms
int maxbias; // size of vbiasall array int maxbias; // size of vbiasall array
inline int sbmask(int j) { inline int sbmask(int j) const {
return j >> SBBITS & 3; return j >> SBBITS & 3;
} }

3
src/compute_orientorder_atom.cpp
View File

@ -102,7 +102,7 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg
if (qlist[iw] > qmax) qmax = qlist[iw]; if (qlist[iw] > qmax) qmax = qlist[iw];
} }
iarg += nqlist; iarg += nqlist;
if (strcmp(arg[iarg],"components") == 0) { } else if (strcmp(arg[iarg],"components") == 0) {
qlcompflag = 1; qlcompflag = 1;
if (iarg+2 > narg) if (iarg+2 > narg)
error->all(FLERR,"Illegal compute orientorder/atom command"); error->all(FLERR,"Illegal compute orientorder/atom command");
@ -118,7 +118,6 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg
if (iqlcomp < 0) if (iqlcomp < 0)
error->all(FLERR,"Illegal compute orientorder/atom command"); error->all(FLERR,"Illegal compute orientorder/atom command");
iarg += 2; iarg += 2;
}
} else if (strcmp(arg[iarg],"cutoff") == 0) { } else if (strcmp(arg[iarg],"cutoff") == 0) {
if (iarg+2 > narg) if (iarg+2 > narg)
error->all(FLERR,"Illegal compute orientorder/atom command"); error->all(FLERR,"Illegal compute orientorder/atom command");

2
src/create_bonds.h
View File

@ -30,7 +30,7 @@ class CreateBonds : protected Pointers {
void command(int, char **); void command(int, char **);
private: private:
inline int sbmask(int j) { inline int sbmask(int j) const {
return j >> SBBITS & 3; return j >> SBBITS & 3;
} }
}; };

2
src/delete_atoms.h
View File

@ -45,7 +45,7 @@ class DeleteAtoms : protected Pointers {
void recount_topology(); void recount_topology();
void options(int, char **); void options(int, char **);
inline int sbmask(int j) { inline int sbmask(int j) const {
return j >> SBBITS & 3; return j >> SBBITS & 3;
} }

2
src/domain.cpp
View File

@ -60,6 +60,7 @@ enum{LAYOUT_UNIFORM,LAYOUT_NONUNIFORM,LAYOUT_TILED}; // several files
Domain::Domain(LAMMPS *lmp) : Pointers(lmp) Domain::Domain(LAMMPS *lmp) : Pointers(lmp)
{ {
box_exist = 0; box_exist = 0;
box_change = 0;
dimension = 3; dimension = 3;
nonperiodic = 0; nonperiodic = 0;
@ -1697,6 +1698,7 @@ int Domain::ownatom(int id, double *x, imageint *image, int shrinkexceed)
void Domain::set_lattice(int narg, char **arg) void Domain::set_lattice(int narg, char **arg)
{ {
if (lattice) delete lattice; if (lattice) delete lattice;
lattice = NULL;
lattice = new Lattice(lmp,narg,arg); lattice = new Lattice(lmp,narg,arg);
} }

55
src/fix_deform.cpp
View File

@ -37,7 +37,7 @@ using namespace LAMMPS_NS;
using namespace FixConst; using namespace FixConst;
using namespace MathConst; using namespace MathConst;
enum{NONE,FINAL,DELTA,SCALE,VEL,ERATE,TRATE,VOLUME,WIGGLE,VARIABLE}; enum{NONE=0,FINAL,DELTA,SCALE,VEL,ERATE,TRATE,VOLUME,WIGGLE,VARIABLE};
enum{ONE_FROM_ONE,ONE_FROM_TWO,TWO_FROM_ONE}; enum{ONE_FROM_ONE,ONE_FROM_TWO,TWO_FROM_ONE};
// same as domain.cpp, fix_nvt_sllod.cpp, compute_temp_deform.cpp // same as domain.cpp, fix_nvt_sllod.cpp, compute_temp_deform.cpp
@ -52,6 +52,7 @@ rfix(NULL), irregular(NULL), set(NULL)
if (narg < 4) error->all(FLERR,"Illegal fix deform command"); if (narg < 4) error->all(FLERR,"Illegal fix deform command");
no_change_box = 1; no_change_box = 1;
restart_global = 1;
nevery = force->inumeric(FLERR,arg[3]); nevery = force->inumeric(FLERR,arg[3]);
if (nevery <= 0) error->all(FLERR,"Illegal fix deform command"); if (nevery <= 0) error->all(FLERR,"Illegal fix deform command");
@ -59,12 +60,7 @@ rfix(NULL), irregular(NULL), set(NULL)
// set defaults // set defaults
set = new Set[6]; set = new Set[6];
set[0].style = set[1].style = set[2].style = memset(set,0,6*sizeof(Set));
set[3].style = set[4].style = set[5].style = NONE;
set[0].hstr = set[1].hstr = set[2].hstr =
set[3].hstr = set[4].hstr = set[5].hstr = NULL;
set[0].hratestr = set[1].hratestr = set[2].hratestr =
set[3].hratestr = set[4].hratestr = set[5].hratestr = NULL;
// parse arguments // parse arguments
@ -343,11 +339,9 @@ rfix(NULL), irregular(NULL), set(NULL)
set[i].hi_initial = domain->boxhi[i]; set[i].hi_initial = domain->boxhi[i];
set[i].vol_initial = domain->xprd * domain->yprd * domain->zprd; set[i].vol_initial = domain->xprd * domain->yprd * domain->zprd;
} }
for (int i = 3; i < 6; i++) { set[3].tilt_initial = domain->yz;
if (i == 5) set[i].tilt_initial = domain->xy; set[4].tilt_initial = domain->xz;
else if (i == 4) set[i].tilt_initial = domain->xz; set[5].tilt_initial = domain->xy;
else if (i == 3) set[i].tilt_initial = domain->yz;
}
// reneighboring only forced if flips can occur due to shape changes // reneighboring only forced if flips can occur due to shape changes
@ -955,6 +949,43 @@ void FixDeform::end_of_step()
if (kspace_flag) force->kspace->setup(); if (kspace_flag) force->kspace->setup();
} }
/* ----------------------------------------------------------------------
write Set data to restart file
------------------------------------------------------------------------- */
void FixDeform::write_restart(FILE *fp)
{
if (comm->me == 0) {
int size = 6*sizeof(Set);
fwrite(&size,sizeof(int),1,fp);
fwrite(set,sizeof(Set),6,fp);
}
}
/* ----------------------------------------------------------------------
use selected state info from restart file to restart the Fix
------------------------------------------------------------------------- */
void FixDeform::restart(char *buf)
{
int samestyle = 1;
Set *set_restart = (Set *) buf;
for (int i=0; i<6; ++i) {
// restore data from initial state
set[i].lo_initial = set_restart[i].lo_initial;
set[i].hi_initial = set_restart[i].hi_initial;
set[i].vol_initial = set_restart[i].vol_initial;
set[i].tilt_initial = set_restart[i].tilt_initial;
// check if style settings are consitent (should do the whole set?)
if (set[i].style != set_restart[i].style)
samestyle = 0;
if (set[i].substyle != set_restart[i].substyle)
samestyle = 0;
}
if (!samestyle)
error->all(FLERR,"Fix deform settings not consistent with restart");
}
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixDeform::options(int narg, char **arg) void FixDeform::options(int narg, char **arg)

2
src/fix_deform.h
View File

@ -35,6 +35,8 @@ class FixDeform : public Fix {
void init(); void init();
virtual void pre_exchange(); virtual void pre_exchange();
virtual void end_of_step(); virtual void end_of_step();
virtual void write_restart(FILE *);
virtual void restart(char *buf);
double memory_usage(); double memory_usage();
protected: protected:

2
src/fix_lineforce.cpp
View File

@ -29,6 +29,8 @@ using namespace FixConst;
FixLineForce::FixLineForce(LAMMPS *lmp, int narg, char **arg) : FixLineForce::FixLineForce(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg) Fix(lmp, narg, arg)
{ {
dynamic_group_allow = 1;
if (narg != 6) error->all(FLERR,"Illegal fix lineforce command"); if (narg != 6) error->all(FLERR,"Illegal fix lineforce command");
xdir = force->numeric(FLERR,arg[3]); xdir = force->numeric(FLERR,arg[3]);
ydir = force->numeric(FLERR,arg[4]); ydir = force->numeric(FLERR,arg[4]);

2
src/fix_planeforce.cpp
View File

@ -29,6 +29,8 @@ using namespace FixConst;
FixPlaneForce::FixPlaneForce(LAMMPS *lmp, int narg, char **arg) : FixPlaneForce::FixPlaneForce(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg) Fix(lmp, narg, arg)
{ {
dynamic_group_allow = 1;
if (narg != 6) error->all(FLERR,"Illegal fix planeforce command"); if (narg != 6) error->all(FLERR,"Illegal fix planeforce command");
xdir = force->numeric(FLERR,arg[3]); xdir = force->numeric(FLERR,arg[3]);
ydir = force->numeric(FLERR,arg[4]); ydir = force->numeric(FLERR,arg[4]);

2
src/pair.h
View File

@ -245,7 +245,7 @@ class Pair : protected Pointers {
ubuf(int arg) : i(arg) {} ubuf(int arg) : i(arg) {}
}; };
inline int sbmask(int j) { inline int sbmask(int j) const {
return j >> SBBITS & 3; return j >> SBBITS & 3;
} }
}; };

1
src/pair_lj_smooth_linear.h
View File

@ -14,6 +14,7 @@
#ifdef PAIR_CLASS #ifdef PAIR_CLASS
PairStyle(lj/smooth/linear,PairLJSmoothLinear) PairStyle(lj/smooth/linear,PairLJSmoothLinear)
PairStyle(lj/sf,PairLJSmoothLinear)
#else #else

12
src/set.cpp
View File

@ -327,15 +327,18 @@ void Set::command(int narg, char **arg)
ximageflag = yimageflag = zimageflag = 0; ximageflag = yimageflag = zimageflag = 0;
if (strcmp(arg[iarg+1],"NULL") != 0) { if (strcmp(arg[iarg+1],"NULL") != 0) {
ximageflag = 1; ximageflag = 1;
ximage = force->inumeric(FLERR,arg[iarg+1]); if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else ximage = force->inumeric(FLERR,arg[iarg+1]);
} }
if (strcmp(arg[iarg+2],"NULL") != 0) { if (strcmp(arg[iarg+2],"NULL") != 0) {
yimageflag = 1; yimageflag = 1;
yimage = force->inumeric(FLERR,arg[iarg+2]); if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) varparse(arg[iarg+2],2);
else yimage = force->inumeric(FLERR,arg[iarg+2]);
} }
if (strcmp(arg[iarg+3],"NULL") != 0) { if (strcmp(arg[iarg+3],"NULL") != 0) {
zimageflag = 1; zimageflag = 1;
zimage = force->inumeric(FLERR,arg[iarg+3]); if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) varparse(arg[iarg+3],3);
else zimage = force->inumeric(FLERR,arg[iarg+3]);
} }
if (ximageflag && ximage && !domain->xperiodic) if (ximageflag && ximage && !domain->xperiodic)
error->all(FLERR, error->all(FLERR,
@ -811,6 +814,9 @@ void Set::set(int keyword)
int xbox = (atom->image[i] & IMGMASK) - IMGMAX; int xbox = (atom->image[i] & IMGMASK) - IMGMAX;
int ybox = (atom->image[i] >> IMGBITS & IMGMASK) - IMGMAX; int ybox = (atom->image[i] >> IMGBITS & IMGMASK) - IMGMAX;
int zbox = (atom->image[i] >> IMG2BITS) - IMGMAX; int zbox = (atom->image[i] >> IMG2BITS) - IMGMAX;
if (varflag1) ximage = static_cast<int>(xvalue);
if (varflag2) yimage = static_cast<int>(yvalue);
if (varflag3) zimage = static_cast<int>(zvalue);
if (ximageflag) xbox = ximage; if (ximageflag) xbox = ximage;
if (yimageflag) ybox = yimage; if (yimageflag) ybox = yimage;
if (zimageflag) zbox = zimage; if (zimageflag) zbox = zimage;