USER-INTEL upgrade from M Brown

Fixes: nve, npt, nvt, nvt/sllod :l

Improper Styles: cvff, harmonic :l

Pair Styles: buck/coul/cut, buck/coul/long, buck, eam, gayberne,

charmm/coul/long, lj/cut, lj/cut/coul/long, sw, tersoff :l

K-Space Styles: pppm :l

charmm/coul/long, lj/cut, lj/cut/coul/long, lj/long/coul/long, sw, tersoff :l

K-Space Styles: pppm, pppm/disp :l

:ule

[Speed-ups to expect:]

LAMMPS {without using other acceleration packages} as these are

under active development (and subject to performance changes). The

measurements were performed using the input files available in

the src/USER-INTEL/TEST directory. These are scalable in size; the

results given are with 512K particles (524K for Liquid Crystal).

Most of the simulations are standard LAMMPS benchmarks (indicated

by the filename extension in parenthesis) with modifications to the

run length and to add a warmup run (for use with offload

benchmarks).

the src/USER-INTEL/TEST directory with the provided run script. 

These are scalable in size; the results given are with 512K 

particles (524K for Liquid Crystal). Most of the simulations are 

standard LAMMPS benchmarks (indicated by the filename extension in

parenthesis) with modifications to the run length and to add a 

warmup run (for use with offload benchmarks).

:c,image(JPG/user_intel.png)

Results are speedups obtained on Intel Xeon E5-2697v4 processors

(code-named Broadwell) and Intel Xeon Phi 7250 processors

(code-named Knights Landing) with "18 Jun 2016" LAMMPS built with

Intel Parallel Studio 2016 update 3. Results are with 1 MPI task

(code-named Knights Landing) with "June 2017" LAMMPS built with

Intel Parallel Studio 2017 update 2. Results are with 1 MPI task

per physical core. See {src/USER-INTEL/TEST/README} for the raw

simulation rates and instructions to reproduce.

:line

[Accuracy and order of operations:]

In most molecular dynamics software, parallelization parameters

(# of MPI, OpenMP, and vectorization) can change the results due

to changing the order of operations with finite-precision 

calculations. The USER-INTEL package is deterministic. This means

that the results should be reproducible from run to run with the

{same} parallel configurations and when using determinstic 

libraries or library settings (MPI, OpenMP, FFT). However, there

are differences in the USER-INTEL package that can change the

order of operations compared to LAMMPS without acceleration:

Neighbor lists can be created in a different order :ulb,l

Bins used for sorting atoms can be oriented differently :l

The default stencil order for PPPM is 7. By default, LAMMPS will 

calculate other PPPM parameters to fit the desired acuracy with 

this order :l

The {newton} setting applies to all atoms, not just atoms shared

between MPI tasks :l

Vectorization can change the order for adding pairwise forces :l

:ule

The precision mode (described below) used with the USER-INTEL 

package can change the {accuracy} of the calculations. For the 

default {mixed} precision option, calculations between pairs or 

triplets of atoms are performed in single precision, intended to 

be within the inherent error of MD simulations. All accumulation

is performed in double precision to prevent the error from growing 

with the number of atoms in the simulation. {Single} precision

mode should not be used without appropriate validation.

:line

[Quick Start for Experienced Users:]

LAMMPS should be built with the USER-INTEL package installed.

Simulations should be run with 1 MPI task per physical {core},

not {hardware thread}.

For Intel Xeon CPUs:

Edit src/MAKE/OPTIONS/Makefile.intel_cpu_intelmpi as necessary. :ulb,l

If using {kspace_style pppm} in the input script, add "neigh_modify binsize cutoff" and "kspace_modify diff ad" to the input script for better

performance.  Cutoff should be roughly the neighbor list cutoff.  By

default the binsize is half the neighbor list cutoff.  :l

"-pk intel 0 omp 2 -sf intel" added to LAMMPS command-line :l

Set the environment variable KMP_BLOCKTIME=0 :l

"-pk intel 0 omp $t -sf intel" added to LAMMPS command-line :l

$t should be 2 for Intel Xeon CPUs and 2 or 4 for Intel Xeon Phi :l

For some of the simple 2-body potentials without long-range

electrostatics, performance and scalability can be better with

the "newton off" setting added to the input script :l

If using {kspace_style pppm} in the input script, add 

"kspace_modify diff ad" for better performance :l

:ule

For Intel Xeon Phi CPUs for simulations without {kspace_style

pppm} in the input script :

For Intel Xeon Phi CPUs:

Edit src/MAKE/OPTIONS/Makefile.knl as necessary. :ulb,l

Runs should be performed using MCDRAM. :l

"-pk intel 0 omp 2 -sf intel" {or} "-pk intel 0 omp 4 -sf intel"

should be added to the LAMMPS command-line. Choice for best

performance will depend on the simulation. :l

Runs should be performed using MCDRAM. :ulb,l

:ule

For Intel Xeon Phi CPUs for simulations with {kspace_style

pppm} in the input script:

Edit src/MAKE/OPTIONS/Makefile.knl as necessary. :ulb,l

Runs should be performed using MCDRAM. :l

Add "neigh_modify binsize 3" to the input script for better

performance. :l

Add "kspace_modify diff ad" to the input script for better

performance. :l

export KMP_AFFINITY=none :l

"-pk intel 0 omp 3 lrt yes -sf intel" or "-pk intel 0 omp 1 lrt yes

-sf intel" added to LAMMPS command-line. Choice for best performance

will depend on the simulation. :l

For simulations using {kspace_style pppm} on Intel CPUs 

supporting AVX-512:

Add "kspace_modify diff ad" to the input script :ulb,l

The command-line option should be changed to 

"-pk intel 0 omp $r lrt yes -sf intel" where $r is the number of 

threads minus 1. :l

Do not use thread affinity (set KMP_AFFINITY=none) :l

The "newton off" setting may provide better scalability :l

:ule

For Intel Xeon Phi coprocessors (Offload):

[Building LAMMPS with the USER-INTEL package:]

NOTE: See the src/USER-INTEL/README file for additional flags that

might be needed for best performance on Intel server processors

code-named "Skylake".

The USER-INTEL package must be installed into the source directory:

make yes-user-intel :pre

NOTE: The USER-INTEL package will perform better with modifications

to the input script when "PPPM"_kspace_style.html is used:

"kspace_modify diff ad"_kspace_modify.html and "neigh_modify binsize

3"_neigh_modify.html should be added to the input script.

"kspace_modify diff ad"_kspace_modify.html should be added to the 

input script.

Long-Range Thread (LRT) mode is an option to the "package

intel"_package.html command that can improve performance when using

environment variable "KMP_AFFINITY=none". LRT mode is not supported

when using offload.

NOTE: Changing the "newton"_newton.html setting to off can improve

performance and/or scalability for simple 2-body potentials such as

lj/cut or when using LRT mode on processors supporting AVX-512.

Not all styles are supported in the USER-INTEL package. You can mix

the USER-INTEL package with styles from the "OPT"_accelerate_opt.html

package or the "USER-OMP package"_accelerate_omp.html. Of course,

Brown, W.M., Carrillo, J.-M.Y., Mishra, B., Gavhane, N., Thakker, F.M., De Kraker, A.R., Yamada, M., Ang, J.A., Plimpton, S.J., "Optimizing Classical Molecular Dynamics in LAMMPS," in Intel Xeon Phi Processor High Performance Programming: Knights Landing Edition, J. Jeffers, J. Reinders, A. Sodani, Eds. Morgan Kaufmann. :ulb,l

Brown, W. M., Semin, A., Hebenstreit, M., Khvostov, S., Raman, K., Plimpton, S.J. Increasing Molecular Dynamics Simulation Rates with an 8-Fold Increase in Electrical Power Efficiency. 2016 International Conference for High Performance Computing. In press. :l

Brown, W. M., Semin, A., Hebenstreit, M., Khvostov, S., Raman, K., Plimpton, S.J. "Increasing Molecular Dynamics Simulation Rates with an 8-Fold Increase in Electrical Power Efficiency."_http://dl.acm.org/citation.cfm?id=3014915 2016 High Performance Computing, Networking, Storage and Analysis, SC16: International Conference (pp. 82-95). :l

Brown, W.M., Carrillo, J.-M.Y., Gavhane, N., Thakkar, F.M., Plimpton, S.J. Optimizing Legacy Molecular Dynamics Software with Directive-Based Offload. Computer Physics Communications. 2015. 195: p. 95-101. :l

:ule

ID, group-ID are documented in "fix"_fix.html command :ulb,l

neb = style name of this fix command :l

Kspring = parallel spring constant (force/distance units or force units) :l

Kspring = parallel spring constant (force/distance units or force units, see nudge keyword) :l

zero or more keyword/value pairs may be appended :l

keyword = {nudg_style} or {perp} or {freend} or {freend_k_spring} :l

  {nudg_style} value = {neigh} or {idealpos}

    {neigh} = the parallel nudging force is calculated from the distances to neighbouring replicas (in this case, Kspring is in force/distance units)

    {idealpos} = the parallel nudging force is proportional to the distance between the replica and its interpolated ideal position (in this case Kspring is in force units)

  {perp} value {none} or kspring2

    {none} = no perpendicular spring force is applied

    {kspring2} = spring constant for the perpendicular nudging force (in force/distance units)

  {freeend} value = {none} or {ini} or {final} or {finaleini} or {final2eini}

    {none} = no nudging force is applied to the first and last replicas

    {ini} = set the first replica to be a free end 

    {final} = set the last replica to be a free end

    {finaleini} = set the last replica to be a free end and set its target energy as that of the first replica

    {final2eini} = same as {finaleini} plus prevent intermediate replicas to have a lower energy than the first replica

  {freeend_kspring}  value = kspring3

    kspring3 = spring constant of the perpendicular spring force (per distance units)    

   :pre

keyword = {nudge} or {perp} or {ends} :l

  {nudge} value = {neigh} or {ideal}

    {neigh} = parallel nudging force based on distance to neighbor replicas (Kspring = force/distance units)

    {ideal} = parallel nudging force based on interpolated ideal position (Kspring = force units)

  {perp} value = {Kspring2}

    {Kspring2} = spring constant for perpendicular nudging force (force/distance units)

  {end} values = estyle Kspring3

    {estyle} = {first} or {last} or {last/efirst} or {last/efirst/middle}

      {first} = apply force to first replica

      {last} = apply force to last replica

      {last/efirst} = apply force to last replica and set its target energy to that of first replica

      {last/efirst/middle} = same as {last/efirst} plus prevent middle replicas having lower energy than first replica

    {Kspring3} = spring constant for target energy term (1/distance units) :pre

[Examples:]

fix 1 active neb 10.0

fix 2 all neb 1.0 perp 1.0 freeend final

fix 1 all neb 1.0 nudg_style idealpos freeend final2eini freend_kspring 1:pre

fix 2 all neb 1.0 perp 1.0 end last

fix 2 all neb 1.0 perp 1.0 end first end last

fix 1 all neb 1.0 nudge ideal end last/efirst 1 :pre

[Description:]

Add a nudging force to atoms in the group for a multi-replica

Add nudging forces to atoms in the group for a multi-replica

simulation run via the "neb"_neb.html command to perform a nudged

elastic band (NEB) calculation for finding the transition state.

Hi-level explanations of NEB are given with the "neb"_neb.html command

and in "Section_howto 5"_Section_howto.html#howto_5 of the manual.

The fix neb command must be used with the "neb" command and defines

how nudging inter-replica forces are computed.  A NEB calculation is

how inter-replica nudging forces are computed.  A NEB calculation is

divided in two stages. In the first stage n replicas are relaxed

toward a MEP and in a second stage, the climbing image scheme (see

"(Henkelman2)"_#Henkelman2) is turned on so that the replica having

the highest energy relaxes toward the saddle point (i.e. the point of

highest energy along the MEP).

One purpose of the nudging forces is to keep the replicas equally

spaced.  During the NEB, the 3N-length vector of interatomic force Fi

= -Grad(V) of replicas i is altered. For all intermediate replicas

(i.e. for 1<i<n) but the climbing replica the force vector

becomes:

Fi = -Grad(V) + (Grad(V) dot That) That + Fnudgparallel + Fspringperp :pre

That is the unit "tangent" vector for replica i and is a function of

Ri, Ri-1, Ri+1, and the potential energy of the 3 replicas; it points

roughly in the direction of (Ri+i - Ri-1) (see the

"(Henkelman1)"_#Henkelman1 paper for details).  Ri are the atomic

coordinates of replica i; Ri-1 and Ri+1 are the coordinates of its

neighbor replicas.  The term (Grad(V) dot That) is used to remove the

toward a MEP until convergence.  In the second stage, the climbing

image scheme (see "(Henkelman2)"_#Henkelman2) is enabled, so that the

replica having the highest energy relaxes toward the saddle point

(i.e. the point of highest energy along the MEP), and a second

relaxation is performed.

A key purpose of the nudging forces is to keep the replicas equally

spaced.  During the NEB calculation, the 3N-length vector of

interatomic force Fi = -Grad(V) for each replica I is altered.  For

all intermediate replicas (i.e. for 1 < I < N, except the climbing

replica) the force vector becomes:

Fi = -Grad(V) + (Grad(V) dot T') T' + Fnudge_parallel + Fspring_perp :pre

T' is the unit "tangent" vector for replica I and is a function of Ri,

Ri-1, Ri+1, and the potential energy of the 3 replicas; it points

roughly in the direction of (Ri+i - Ri-1); see the

"(Henkelman1)"_#Henkelman1 paper for details.  Ri are the atomic

coordinates of replica I; Ri-1 and Ri+1 are the coordinates of its

neighbor replicas.  The term (Grad(V) dot T') is used to remove the

component of the gradient parallel to the path which would tend to

distribute the replica unevenly along the path.  Fnudgparallel is an

artificial nudging force which is applied only in the tangent direction

and which maintains the replicas equally spaced (see below for more

information).  Fspringperp is an optinal artificial spring which is

applied only perpendicular to the tangent and which prevent the paths

from forming too acute kinks (see below for more information).

distribute the replica unevenly along the path.  Fnudge_parallel is an

artificial nudging force which is applied only in the tangent

direction and which maintains the equal spacing between replicas (see

below for more information).  Fspring_perp is an optional artificial

spring which is applied only perpendicular to the tangent and which

prevent the paths from forming acute kinks (see below for more

information).

The keyword {nudg_style} allow to specify how to parallel

nudging force is computed. With a value of idealpos, the spring 

force is computed as suggested in "(E)"_#E :

In the second stage of the NEB calculation, the interatomic force Fi

for the climbing replica (the replica of highest energy after the

first stage) is changed to:

Fnudgparallel=-{Kspring}* (RD-RDideal)/(2 meanDist) :pre

Fi = -Grad(V) + 2 (Grad(V) dot T') T' :pre

where RD is the "reaction coordinate" see "neb"_neb.html section, and

RDideal is the ideal RD for which all the images are equally spaced

(i.e. RDideal = (i-1)*meanDist when the climbing image is off, where i

is the replica number). The meanDist is the average distance between

replicas.

and the relaxation procedure is continued to a new converged MEP.

:line

When {nudg_style} has a value of neigh (or by default), the parallel 

nudging force is computed as in "(Henkelman1)"_#Henkelman1 by 

connecting each intermediate replica with the previous and the next 

image:

The keyword {nudge} specifies how the parallel nudging force is

computed.  With a value of {neigh}, the parallel nudging force is

computed as in "(Henkelman1)"_#Henkelman1 by connecting each

intermediate replica with the previous and the next image:

Fnudgparallel= {Kspring}* (|Ri+1 - Ri| - |Ri - Ri-1|) :pre

Fnudge_parallel = {Kspring} * (|Ri+1 - Ri| - |Ri - Ri-1|) :pre

The parallel nudging force associated with the key word idealpos should

usually be more efficient at keeping the images equally spaced.

Note that in this case the specified {Kspring) is in force/distance

units.

:line

With a value of {ideal}, the spring force is computed as suggested in

"(WeinenE)"_#WeinenE :

The keyword {perp} allows to add a spring force perpendicular to the

path in order to prevent the path from becoming too kinky. It can

improve significantly the convergence of the NEB when the resolution

is poor (i.e. when too few images are used) (see "(Maras)"_#Maras1).

The perpendicular spring force is given by

Fnudge_parallel = -{Kspring} * (RD-RDideal) / (2 * meanDist) :pre

Fspringperp = {Kspringperp} * f(Ri-1,Ri,Ri+1) (Ri+1 + Ri-1 - 2 Ri) :pre

where RD is the "reaction coordinate" see "neb"_neb.html section, and

RDideal is the ideal RD for which all the images are equally spaced.

I.e. RDideal = (I-1)*meanDist when the climbing replica is off, where

I is the replica number).  The meanDist is the average distance

between replicas.  Note that in this case the specified {Kspring) is

in force units.

f(Ri-1 Ri R+1) is a smooth scalar function of the angle Ri-1 Ri

Ri+1. It is equal to 0 when the path is straight and is equal to 1

when the angle Ri-1 Ri Ri+1 is accute. f(Ri-1 Ri R+1) is defined in

"(Jonsson)"_#Jonsson

Note that the {ideal} form of nudging can often be more effective at

keeping the replicas equally spaced.

:line

By default, the force acting on the first and last replicas is not

altered so that during the NEB relaxation, these ending replicas relax

toward local minima. However it is possible to use the key word

{freeend} to allow either the initial or the final replica to relax

toward a MEP while constraining its energy.  The interatomic force Fi

for the free end image becomes :

The keyword {perp} adds a spring force perpendicular to the path in

order to prevent the path from becoming too kinky, with magnitude It

can significantly improve the convergence of the NEB calculation when

the resolution is poor.  I.e. when too few replicas are used; see

"(Maras)"_#Maras1 for details.

Fi = -Grad(V)+ (Grad(V) dot That + (E-ETarget)*kspring3) That,  {when} Grad(V) dot That < 0

Fi = -Grad(V)+ (Grad(V) dot That + (ETarget- E)*kspring3) That, {when} Grad(V) dot That > 0

:pre

The perpendicular spring force is given by

where E is the energy of the free end replica and ETarget is the

target energy.

When the value {ini} ({final}) is used after the keyword {freeend},

the first (last) replica is considered as a free end. The target

energy is set to the energy of the replica at starting of the NEB

calculation. When the value {finaleini} or {final2eini} is used the

last image is considered as a free end and the target energy is equal

to the energy of the first replica (which can evolve during the NEB

relaxation).  With the value {finaleini}, when the initial path is too

far from the MEP, an intermediate repilica might relax "faster" and

get a lower energy than the last replica. The benefit of the free end

is then lost since this intermediate replica will relax toward a local

minima. This behavior can be prevented by using the value {final2eini}

which remove entirely the contribution of the gradient for all

intermediate replica which have a lower energy than the initial one

thus preventing these replicae to over-relax.  After converging a NEB

with the {final2eini} value it is recommended to check that all

intermediate replica have a larger energy than the initial

replica. Finally note that if the last replica converges toward a

local minimum with a larger energy than the energy of the first

replica, a free end neb calculation with the value {finaleini} or

{final2eini} cannot reach the convergence criteria.

Fspring_perp = {Kspring2} * F(Ri-1,Ri,Ri+1) (Ri+1 + Ri-1 - 2 Ri) :pre

:line

where {Kspring2} is the specified value.  F(Ri-1 Ri R+1) is a smooth

scalar function of the angle Ri-1 Ri Ri+1.  It is equal to 0.0 when

the path is straight and is equal to 1 when the angle Ri-1 Ri Ri+1 is

acute.  F(Ri-1 Ri R+1) is defined in "(Jonsson)"_#Jonsson.

If {Kspring2} is set to 0.0 (the default) then no perpendicular spring

force is added.

:line

In the second stage of the NEB, the interatomic force Fi for the

climbing replica (which is the replica of highest energy) becomes:

By default, no forces act on the first and last replicas during the

NEB relaxation, so these replicas simply relax toward their respective

local minima.  By using the key word {end}, additional forces can be

applied to the first or last replica, to enable them to relax toward a

MEP while constraining their energy.

Fi = -Grad(V) + 2 (Grad(V) dot That) That :pre

The interatomic force Fi for the specified replica becomes:

Fi = -Grad(V) + (Grad(V) dot T' + (E-ETarget)*Kspring3) T',  {when} Grad(V) dot T' < 0

Fi = -Grad(V) + (Grad(V) dot T' + (ETarget- E)*Kspring3) T', {when} Grad(V) dot T' > 0

:pre

where E is the current energy of the replica and ETarget is the target

energy.  The "spring" constant on the difference in energies is the

specified {Kspring3} value.

When {estyle} is specified as {first}, the force is applied to the

first replica.  When {estyle} is specified as {last}, the force is

applied to the last replica.  Note that the {end} keyword can be used

twice to add forces to both the first and last replicas.

For both these {estyle} settings, the target energy {ETarget} is set

to the initial energy of the replica (at the start of the NEB

calculation).

If the {estyle} is specified as {last/efirst} or {last/efirst/middle},

force is applied to the last replica, but the target energy {ETarget}

is continuously set to the energy of the first replica, as it evolves

during the NEB relaxation.

The difference between these two {estyle} options is as follows.  When

{estyle} is specified as {last/efirst}, no change is made to the

inter-replica force applied to the intermediate replicas (neither

first or last).  If the initial path is too far from the MEP, an

intermediate repilica may relax "faster" and reach a lower energy than

the last replica.  In this case the intermediate replica will be

relaxing toward its own local minima.  This behavior can be prevented

by specifying {estyle} as {last/efirst/middle} which will alter the

inter-replica force applied to intermediate replicas by removing the

contribution of the gradient to the inter-replica force.  This will

only be done if a particular intermediate replica has a lower energy

than the first replica.  This should effectively prevent the

intermediate replicas from over-relaxing.

After converging a NEB calculation using an {estyle} of {last/efirst},

you should check that all intermediate replicas have a larger energy

than the first replica.  If not, then repeat the calculation with an

{estyle} of {last/efirst/middle}.

Finally, note that if the last replica converges toward a local

minimum which has a larger energy than the energy of the first

replica, a NEB calculation using an {estyle} of {last/efirst} or

{last/efirst/middle} cannot reach final convergence.

[Restart, fix_modify, output, run start/stop, minimize info:]

[Default:]

The option defaults are nudg_style = neigh, perp = none, freeend = none and freend_kspring = 1.

The option defaults are nudge = neigh, perp = 0.0, ends is not

specified (no inter-replica force on the end replicas).

:line

[(Henkelman2)] Henkelman, Uberuaga, Jonsson, J Chem Phys, 113,

9901-9904 (2000).

:link(E)

[(E)] E, Ren, Vanden-Eijnden, Phys Rev B, 66, 052301 (2002)

:link(WeinenE)

[(WeinenE)] E, Ren, Vanden-Eijnden, Phys Rev B, 66, 052301 (2002).

:link(Jonsson)

[(Jonsson)] Jonsson, Mills and Jacobsen, in Classical and Quantum

Dynamics in Condensed Phase Simulations, edited by Berne, Ciccotti, and Coker

World Scientific, Singapore, 1998, p. 385

Dynamics in Condensed Phase Simulations, edited by Berne, Ciccotti,

and Coker World Scientific, Singapore, 1998, p 385.

:link(Maras1)

[(Maras)] Maras, Trushin, Stukowski, Ala-Nissila, Jonsson,

Comp Phys Comm, 205, 13-21 (2016)

Comp Phys Comm, 205, 13-21 (2016).

gewald = gewald/disp = 0.0, slab = 1.0, compute = yes, cutoff/adjust =

yes (MSM), pressure/scalar = yes (MSM), fftbench = yes (PPPM), diff = ik

(PPPM), mix/disp = pair, force/disp/real = -1.0, force/disp/kspace = -1.0,

split = 0, tol = 1.0e-6, and disp/auto = no.

split = 0, tol = 1.0e-6, and disp/auto = no. For pppm/intel, order =

order/disp = 7.

:line