some small formating change but does not work anymore

# SOME DESCRIPTIVE TITLE.

# Copyright (C) YEAR Free Software Foundation, Inc.

# FIRST AUTHOR <EMAIL@ADDRESS>, YEAR.

#

#, fuzzy

msgid ""

msgstr ""

"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c

:link(lws,http://lammps.sandia.gov)

:link(ld,Manual.html)

:link(lws,http://lammps.sandia.gov) :link(ld,Manual.html)

:link(lc,Section_commands.html#comm)

:line

fix ID group-ID neb Kspring keyword value :pre

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

neb = style name of this fix command

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

keyword = {idealpos} or {nearestneigh} or {PerpSpring} or {freeend}

  {idealpos} value = none = each replica is attached with a spring to its interpolated ideal position (default value)

  {nearestneigh} value = none = each replica is attached with a spring with the previous and next replica.

  {PerpSpring} value = KspringPerpend

  {freeend} value = ini or final or finalWithRespToIni or finalAndInterWithRespToIni

ID, group-ID are documented in "fix"_fix.html command neb = style name of this

fix command Kspring = parallel spring constant (force/distance units) :ul

keyword = {idealpos} or {neigh} or {perp} or {freeend} {idealpos} value = none =

each replica is attached with a spring to its interpolated ideal position

(default value) {neigh} value = none = each replica is attached with a spring

with the previous and next replica.  {perp} value = spring constant for the

perpendicular spring {freeend} value = ini or final or finaleini or final2eini

[Examples:]

fix 1 active neb 10.0 :pre

fix 2 all neb 1.0 PerpSpring 1.0 freeend final :pre

fix 1 all neb 1.0 nearestneigh freeend finalAndInterWithRespToIni :pre

fix 1 active neb 10.0 :pre fix 2 all neb 1.0 perp 1.0 freeend final :pre fix 1

all neb 1.0 neigh freeend final2eini :pre

[Description:]

Add a nudging force 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 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).

Add a nudging force 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

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) except for the climbing replica the force vector becomes:

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) except for

the climbing replica the force vector becomes:

Fi = -Grad(V) + (Grad(V) dot That) That + Fspringparallel + 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 component of the gradient parallel to the path which would tend to distribute the replica unevenly along the path.

Fspringparallel is an artificial spring 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). 

In the second stage of the NEB, the interatomic force Fi for the climbing replica (which is the replica of highest energy) becomes :

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 component of the gradient parallel to the path which would tend to

distribute the replica unevenly along the path.  Fspringparallel is an

artificial spring 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).

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

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

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

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 : 

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 :

Fi = -Grad(V)+ (Grad(V) dot That + E-ETarget) That when Grad(V) dot That < 0

Fi = -Grad(V)+ (Grad(V) dot That + ETarget- E) That when Grad(V) dot That > 0

Fi = -Grad(V)+ (Grad(V) dot That + E-ETarget) That when Grad(V) dot That < 0 Fi

= -Grad(V)+ (Grad(V) dot That + ETarget- E) That when Grad(V) dot That > 0

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 {finalWithRespToIni} or {finalAndInterWithRespToIni} 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 {finalWithRespToIni}, 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 {finalAndInterWithRespToIni} 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 {finalAndInterWithRespToIni} 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 {finalWithRespToIni} or {finalAndInterWithRespToIni} cannot reach the convergence criteria.

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.

:line

The keywords {idealpos} and {nearestneigh} allow to specify how to parallel spring force is computed. 

If the keyword {idealpos} is used or by default, the spring force is computed as suggested in "(E)"_#E :

The keywords {idealpos} and {neigh} allow to specify how to parallel spring

force is computed.  If the keyword {idealpos} is used or by default, the spring

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

Fspringparallel=-{Kspring}* (RD-RDideal)/(2 meanDist)

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.

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.

If the keyword  {nearestneigh} is used, the parallel spring force is computed as in "(Henkelman1)"_#Henkelman1 by connecting each intermediate replica with the previous and the next image:

If the keyword {neigh} is used, the parallel spring force is computed as in

"(Henkelman1)"_#Henkelman1 by connecting each intermediate replica with the

previous and the next image:

Fspringparallel= {Kspring}* (|Ri+1 - Ri| - |Ri - Ri-1|)

The parallel spring force associated with the key word idealpos should usually be more efficient at keeping the images equally spaced.

The parallel spring force associated with the key word idealpos should usually

be more efficient at keeping the images equally spaced.

:line

The keyword {PerpSpring} 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)"_#Maras).

The perpendicular spring force is given by 

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)"_#Maras).  The perpendicular spring force is given by

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

 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 

 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

:line

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

No information about this fix is written to "binary restart

files"_restart.html.  None of the "fix_modify"_fix_modify.html options

are relevant to this fix.  No global or per-atom quantities are stored

by this fix for access by various "output

commands"_Section_howto.html#howto_15.  No parameter of this fix can

be used with the {start/stop} keywords of the "run"_run.html command.

No information about this fix is written to "binary restart files"_restart.html.

None of the "fix_modify"_fix_modify.html options are relevant to this fix.  No

global or per-atom quantities are stored by this fix for access by various

"output commands"_Section_howto.html#howto_15.  No parameter of this fix can be

used with the {start/stop} keywords of the "run"_run.html command.

The forces due to this fix are imposed during an energy minimization,

as invoked by the "minimize"_minimize.html command via the

"neb"_neb.html command.

The forces due to this fix are imposed during an energy minimization, as invoked

by the "minimize"_minimize.html command via the "neb"_neb.html command.

[Restrictions:]

This command can only be used if LAMMPS was built with the REPLICA

package.  See the "Making LAMMPS"_Section_start.html#start_3 section

for more info on packages.

This command can only be used if LAMMPS was built with the REPLICA package.  See

the "Making LAMMPS"_Section_start.html#start_3 section for more info on

packages.

[Related commands:]

[Default:] none

:link(Henkelman)

[(Henkelman1)] Henkelman and Jonsson, J Chem Phys, 113, 9978-9985 (2000).

:link(Henkelman1) [(Henkelman1)] Henkelman and Jonsson, J Chem Phys, 113,

9978-9985 (2000).

:link(Henkelman)

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

:link(Henkelman2) [(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(E) [(E)] 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

: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

:link(Maras)

[(Maras)] Maras, Trushin, Stukowski, Ala-Nissila, Jonsson, Comp Phys Comm, 205, 13-21 (2016)

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

Phys Comm, 205, 13-21 (2016)

#include "memory.h"

#include "error.h"

#include "force.h"

#include "math_const.h"

using namespace LAMMPS_NS;

using namespace FixConst;

using namespace MathConst;

enum{SINGLE_PROC_DIRECT,SINGLE_PROC_MAP,MULTI_PROC};

/* ---------------------------------------------------------------------- */

  int iarg =4;

  while (iarg < narg){

    if (strcmp (arg[iarg],"idealpos")==0) 

      {NEBLongRange = true;

    if (strcmp (arg[iarg],"idealpos")==0){

      NEBLongRange = true;

      iarg+=1;}

    else if (strcmp (arg[iarg],"nearestneigh")==0)

      {NEBLongRange = false;

    else if (strcmp (arg[iarg],"neigh")==0){

      NEBLongRange = false;

      StandardNEB = true;

      iarg+=1;}

    else if (strcmp (arg[iarg],"PerpSpring")==0) 

      {PerpSpring=true;

    else if (strcmp (arg[iarg],"perp")==0){

      PerpSpring=true;

      kspringPerp = force->numeric(FLERR,arg[iarg+1]);

      if (kspringPerp < 0.0) error->all(FLERR,"Illegal fix neb command. The perpendicular spring force was not provided properly");

      iarg+=2;

    }

    else if (strcmp (arg[iarg],"freeend")==0) 

      {

    else if (strcmp (arg[iarg],"freeend")==0){

      if (strcmp (arg[iarg+1],"ini")==0)

        FreeEndIni=true;

      else if (strcmp (arg[iarg+1],"final")==0)

        FreeEndFinal=true;

      else if (strcmp (arg[iarg+1],"final")==0)

        FreeEndFinal=true;

	else if (strcmp (arg[iarg+1],"finalWithRespToIni")==0) 

      else if (strcmp (arg[iarg+1],"finaleini")==0)

	FreeEndFinalWithRespToEIni=true;

	else if (strcmp (arg[iarg+1],"finalAndInterWithRespToIni")==0) 

	  {FinalAndInterWithRespToEIni=true;

      else if (strcmp (arg[iarg+1],"final2eini")==0){

        FinalAndInterWithRespToEIni=true;

        FreeEndFinalWithRespToEIni=true;}

        iarg+=2;}

      else {error->all(FLERR,"Illegal fix neb command. Unknown keyword");}

  int *iroots = new int[nreplica];

  MPI_Group uworldgroup,rootgroup;

  if (NEBLongRange ){

    for (int iIm =0; iIm < nreplica;iIm++)

      {

	iroots[iIm]=universe->root_proc[iIm];

      }

    for (int iIm =0; iIm < nreplica;iIm++){

      iroots[iIm]=universe->root_proc[iIm];}

    MPI_Comm_group(uworld, &uworldgroup);

    MPI_Group_incl(uworldgroup, nreplica, iroots, &rootgroup);

    //    MPI_Comm_create_group(uworld, rootgroup, 0, &rootworld);

    MPI_Comm_create(uworld, rootgroup, &rootworld);

  }

  MPI_Status status;

  MPI_Request request;

  double vIni =0.0;

  // veng = PE of this replica

  // vprev,vnext = PEs of adjacent replicas

  // only proc 0 in each replica communicates

  vprev=vnext=veng = pe->compute_scalar();

    MPI_Bcast(&vnext,1,MPI_DOUBLE,0,world);

  }

  if (FreeEndFinal)

    {

      if (update->ntimestep==0)

	{EFinalIni = veng;}

  if (FreeEndFinal){

    if (update->ntimestep==0){EFinalIni = veng;}

  }

  if (ireplica==0)

    vIni=veng;

  if (FreeEndFinalWithRespToEIni ){

    if ( me ==0){

      int procFirst;

      MPI_Bcast(&vIni,1,MPI_DOUBLE,0,world);

    }

  }

    if (FreeEndIni && ireplica==0 )

    {

    if (FreeEndIni && ireplica==0 ){

      if (me == 0 )

      if (update->ntimestep==0) 

	{

	if (update->ntimestep==0){

          EIniIni = veng;

          if (cmode == MULTI_PROC) 

            MPI_Bcast(&EIniIni,1,MPI_DOUBLE,0,world);

    }

  // communicate atoms to/from adjacent replicas to fill xprev,xnext

  inter_replica_comm();

  // trigger potential energy computation on next timestep

  if (ireplica < nreplica-1)

    MPI_Irecv(xnext[0],3*nlocal,MPI_DOUBLE,procnext,0,uworld,&request);

    MMY_PI_Irecv(xnext[0],3*nlocal,MPI_DOUBLE,procnext,0,uworld,&request);

  if (ireplica > 0)

    MPI_Send(x[0],3*nlocal,MPI_DOUBLE,procprev,0,uworld);

  if (ireplica < nreplica-1) MPI_Wait(&request,&status);

  dottangrad = 0.0;

  if (ireplica ==nreplica-1){

    for (int i = 0; i < nlocal; i++)

      if (mask[i] & groupbit) {

        plen += delxp*delxp + delyp*delyp + delzp*delzp;

        dottangrad += delxp* f[i][0]+ delyp*f[i][1]+delzp*f[i][2];

        gradlen += f[i][0]*f[i][0] + f[i][1]*f[i][1] + f[i][2]*f[i][2];

        if (FreeEndFinal||FreeEndFinalWithRespToEIni){

          tangent[i][0]=delxp;

          tangent[i][1]=delyp;

    }

  }

  else if (ireplica == 0){

    for (int i = 0; i < nlocal; i++)

      if (mask[i] & groupbit) {

          dotgrad += f[i][0]*fnext[i][0] + f[i][1]*fnext[i][1] +

            f[i][2]*fnext[i][2];

          springF[i][0]=kspringPerp*(delxn-delxp);

          springF[i][1]=kspringPerp*(delyn-delyp);

          springF[i][2]=kspringPerp*(delzn-delzp);

        }

  }

  double dotpathall;

  double dottangradall;

  double dotgradall;

  }

  if(ireplica==0)     dottangrad = dottangrad/(nlen*gradlen);

  if(ireplica==nreplica-1)     dottangrad = dottangrad/(plen*gradlen);

  if(ireplica < nreplica-1)

    {  

  if(ireplica < nreplica-1){

    dotgrad = dotgrad /(gradlen*gradnextlen);

  }

  if(FreeEndIni&&ireplica == 0)

    {

  if(FreeEndIni&&ireplica == 0) {

    if (tlen > 0.0) {

      double dotall;

      MPI_Allreduce(&dot,&dotall,1,MPI_DOUBLE,MPI_SUM,world);

      dot=dotall;

      double tleninv = 1.0/tlen;

      dot *= tleninv;

      if (dot<0)

          f[i][1] += prefactor *tangent[i][1]; 

          f[i][2] += prefactor *tangent[i][2];

        }

    }

  }

  if(FreeEndFinal&&ireplica == nreplica -1)

    {if (tlen > 0.0) {

  if(FreeEndFinal&&ireplica == nreplica -1){

    if (tlen > 0.0) {

      double dotall;

      MPI_Allreduce(&dot,&dotall,1,MPI_DOUBLE,MPI_SUM,world);

      dot=dotall;

    }

  }

  if(FreeEndFinalWithRespToEIni&&ireplica == nreplica -1)

    {if (tlen > 0.0) {

  if(FreeEndFinalWithRespToEIni&&ireplica == nreplica -1){

    if (tlen > 0.0) {

      double dotall;

      MPI_Allreduce(&dot,&dotall,1,MPI_DOUBLE,MPI_SUM,world);

      dot=dotall;

    }

  }

  double lentot = 0;

  double meanDist,idealPos,lenuntilIm,lenuntilClimber;

    {

      if (cmode == SINGLE_PROC_DIRECT or cmode == SINGLE_PROC_MAP)

        {MPI_Allgather(&nlen,1,MPI_DOUBLE,&nlenall[0],1,MPI_DOUBLE,uworld);}

      else

	{

	  /*	  int procRootiIm;

	  double nlentmp;

	  for (int iIm = 0; i < nreplica; i++)

	    {

	      procRootiIm=universe->root_proc[iIm];

	      if (ireplica == iIm && me ==0)

		{ nlentmp=nlen;

		  MPI_Bcast(&nlentmp,1,MPI_DOUBLE,procRootiIm,uworld);

		}

	      else

		{

		  MPI_Bcast(&nlentmp,1,MPI_DOUBLE,procRootiIm,uworld);

		}

	      nlenall[iIm]=nlen;

	    }

      */

      else{

	if (me == 0)

	  MPI_Allgather(&nlen,1,MPI_DOUBLE,&nlenall[0],1,MPI_DOUBLE,rootworld);

	MPI_Bcast(nlenall,nreplica,MPI_DOUBLE,0,world);

      }

      for (int i = 0; i < ireplica; i++)

        lenuntilIm += nlenall[i];

      for (int i = 0; i < nreplica; i++)

        lentot += nlenall[i];

      meanDist = lentot/(nreplica -1);

      if (rclimber>0)

	{

      if (rclimber>0) {

	for (int i = 0; i < rclimber; i++)

	  lenuntilClimber += nlenall[i];

	double meanDistBeforeClimber = lenuntilClimber/rclimber;

	double meanDistAfterClimber = (lentot-lenuntilClimber)/(nreplica-rclimber-1);

	if (ireplica<rclimber)

	  idealPos = ireplica * meanDistBeforeClimber;

	else

	  idealPos = lenuntilClimber+ (ireplica-rclimber)*meanDistAfterClimber;

      }

      else{

        idealPos = ireplica * meanDist;

      }

  if (ireplica == 0 || ireplica == nreplica-1)     return ;

  double AngularContr;

  double pi;

  double thetapath;

  pi = 4.0*atan(1.0);

dotpath = dotpath/(plen*nlen);

AngularContr = 0.5 *(1+cos(pi * dotpath));

  double thetapath;

dotpath = dotpath/(plen*nlen);

AngularContr = 0.5 *(1+cos(MY_PI * dotpath));

  dotSpringTangent=0;

  for (int i = 0; i < nlocal; i++) {

    if (mask[i] & groupbit)

      {

    if (mask[i] & groupbit) {

      dot += f[i][0]*tangent[i][0] + f[i][1]*tangent[i][1] + f[i][2]*tangent[i][2];

      dotSpringTangent+=springF[i][0]*tangent[i][0]+springF[i][1]*tangent[i][1]+springF[i][2]*tangent[i][2];}

  }

  dot=dotall;

  //  double prefactor, prefSpring;

  double ToDisp;

  if (ireplica == rclimber) {

  }

  else {

    if(NEBLongRange)

      {prefactor = -dot - kspring*(lenuntilIm-idealPos)/(2*meanDist);}

    else if (StandardNEB)

      {prefactor = -dot + kspring*(nlen-plen);}

    if (FinalAndInterWithRespToEIni&& veng<vIni)

      {

    if (FinalAndInterWithRespToEIni&& veng<vIni){

      for (int i = 0; i < nlocal; i++)

	if (mask[i] & groupbit) {

	  f[i][0] = 0;

    }

  }

  for (int i = 0; i < nlocal; i++)

    if (mask[i] & groupbit) {

      f[i][0] += prefactor*tangent[i][0] +AngularContr*(springF[i][0] -dotSpringTangent*tangent[i][0]);

      f[i][1] += prefactor*tangent[i][1]+ AngularContr*(springF[i][1] - dotSpringTangent*tangent[i][1]);

      f[i][2] += prefactor*tangent[i][2]+ AngularContr*(springF[i][2] - dotSpringTangent*tangent[i][2]);

    if (ireplica < nreplica-1)

      MPI_Send(x[0],3*nlocal,MPI_DOUBLE,procnext,0,uworld);

    if (ireplica > 0) MPI_Wait(&request,MPI_STATUS_IGNORE);

    if (ireplica < nreplica-1)

      MPI_Irecv(xnext[0],3*nlocal,MPI_DOUBLE,procnext,0,uworld,&request);

    if (ireplica > 0)

        xprev[m][2] = xrecv[i][2];

      }

    }

    if (ireplica < nreplica-1) {

      MPI_Irecv(xrecv[0],3*nebatoms,MPI_DOUBLE,procnext,0,uworld,&requests[0]);

      MPI_Irecv(frecv[0],3*nebatoms,MPI_DOUBLE,procnext,0,uworld,&requests[0]);

  memory->destroy(fnext);

  memory->destroy(springF);

  if (NEBLongRange)

        {memory->destroy(nlenall);

  if (NEBLongRange){

    memory->destroy(nlenall);

    memory->create(nlenall,nreplica,"neb:nlenall");

  }

  if (cmode != SINGLE_PROC_DIRECT) {

    memory->destroy(xsend);

    memory->destroy(fsend);