Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F71778400
fix_neb.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Jul 13, 02:02
Size
27 KB
Mime Type
text/x-c
Expires
Mon, Jul 15, 02:02 (2 d)
Engine
blob
Format
Raw Data
Handle
18994313
Attached To
rLAMMPS lammps
fix_neb.cpp
View Options
/* ----------------------------------------------------------------------
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 for: Emile Maras (CEA, France)
new options for inter-replica forces, first/last replica treatment
------------------------------------------------------------------------- */
#include <mpi.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "fix_neb.h"
#include "universe.h"
#include "update.h"
#include "atom.h"
#include "domain.h"
#include "comm.h"
#include "modify.h"
#include "compute.h"
#include "group.h"
#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
};
#define BUFSIZE 8
/* ---------------------------------------------------------------------- */
FixNEB
::
FixNEB
(
LAMMPS
*
lmp
,
int
narg
,
char
**
arg
)
:
Fix
(
lmp
,
narg
,
arg
),
id_pe
(
NULL
),
pe
(
NULL
),
nlenall
(
NULL
),
xprev
(
NULL
),
xnext
(
NULL
),
fnext
(
NULL
),
springF
(
NULL
),
tangent
(
NULL
),
xsend
(
NULL
),
xrecv
(
NULL
),
fsend
(
NULL
),
frecv
(
NULL
),
tagsend
(
NULL
),
tagrecv
(
NULL
),
xsendall
(
NULL
),
xrecvall
(
NULL
),
fsendall
(
NULL
),
frecvall
(
NULL
),
tagsendall
(
NULL
),
tagrecvall
(
NULL
),
counts
(
NULL
),
displacements
(
NULL
)
{
if
(
narg
<
4
)
error
->
all
(
FLERR
,
"Illegal fix neb command"
);
kspring
=
force
->
numeric
(
FLERR
,
arg
[
3
]);
if
(
kspring
<=
0.0
)
error
->
all
(
FLERR
,
"Illegal fix neb command"
);
// optional params
NEBLongRange
=
false
;
StandardNEB
=
true
;
PerpSpring
=
FreeEndIni
=
FreeEndFinal
=
false
;
FreeEndFinalWithRespToEIni
=
FinalAndInterWithRespToEIni
=
false
;
kspringPerp
=
0.0
;
kspringIni
=
1.0
;
kspringFinal
=
1.0
;
int
iarg
=
4
;
while
(
iarg
<
narg
)
{
if
(
strcmp
(
arg
[
iarg
],
"parallel"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix neb command"
);
if
(
strcmp
(
arg
[
iarg
+
1
],
"ideal"
)
==
0
)
{
NEBLongRange
=
true
;
StandardNEB
=
false
;
}
else
if
(
strcmp
(
arg
[
iarg
+
1
],
"neigh"
)
==
0
)
{
NEBLongRange
=
false
;
StandardNEB
=
true
;
}
else
error
->
all
(
FLERR
,
"Illegal fix neb command"
);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"perp"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix neb command"
);
PerpSpring
=
true
;
kspringPerp
=
force
->
numeric
(
FLERR
,
arg
[
iarg
+
1
]);
if
(
kspringPerp
==
0.0
)
PerpSpring
=
false
;
if
(
kspringPerp
<
0.0
)
error
->
all
(
FLERR
,
"Illegal fix neb command"
);
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"end"
)
==
0
)
{
if
(
iarg
+
3
>
narg
)
error
->
all
(
FLERR
,
"Illegal fix neb command"
);
if
(
strcmp
(
arg
[
iarg
+
1
],
"first"
)
==
0
)
{
FreeEndIni
=
true
;
kspringIni
=
force
->
numeric
(
FLERR
,
arg
[
iarg
+
2
]);
}
else
if
(
strcmp
(
arg
[
iarg
+
1
],
"last"
)
==
0
)
{
FreeEndFinal
=
true
;
FinalAndInterWithRespToEIni
=
false
;
FreeEndFinalWithRespToEIni
=
false
;
kspringFinal
=
force
->
numeric
(
FLERR
,
arg
[
iarg
+
2
]);
}
else
if
(
strcmp
(
arg
[
iarg
+
1
],
"last/efirst"
)
==
0
)
{
FreeEndFinal
=
false
;
FinalAndInterWithRespToEIni
=
false
;
FreeEndFinalWithRespToEIni
=
true
;
kspringFinal
=
force
->
numeric
(
FLERR
,
arg
[
iarg
+
2
]);
}
else
if
(
strcmp
(
arg
[
iarg
+
1
],
"last/efirst/middle"
)
==
0
)
{
FreeEndFinal
=
false
;
FinalAndInterWithRespToEIni
=
true
;
FreeEndFinalWithRespToEIni
=
true
;
kspringFinal
=
force
->
numeric
(
FLERR
,
arg
[
iarg
+
2
]);
}
else
error
->
all
(
FLERR
,
"Illegal fix neb command"
);
iarg
+=
3
;
}
else
error
->
all
(
FLERR
,
"Illegal fix neb command"
);
}
// nreplica = number of partitions
// ireplica = which world I am in universe
// nprocs_universe = # of procs in all replicase
// procprev,procnext = root proc in adjacent replicas
me
=
comm
->
me
;
nprocs
=
comm
->
nprocs
;
nprocs_universe
=
universe
->
nprocs
;
nreplica
=
universe
->
nworlds
;
ireplica
=
universe
->
iworld
;
if
(
ireplica
>
0
)
procprev
=
universe
->
root_proc
[
ireplica
-
1
];
else
procprev
=
-
1
;
if
(
ireplica
<
nreplica
-
1
)
procnext
=
universe
->
root_proc
[
ireplica
+
1
];
else
procnext
=
-
1
;
uworld
=
universe
->
uworld
;
int
*
iroots
=
new
int
[
nreplica
];
MPI_Group
uworldgroup
,
rootgroup
;
if
(
NEBLongRange
)
{
for
(
int
i
=
0
;
i
<
nreplica
;
i
++
)
iroots
[
i
]
=
universe
->
root_proc
[
i
];
MPI_Comm_group
(
uworld
,
&
uworldgroup
);
MPI_Group_incl
(
uworldgroup
,
nreplica
,
iroots
,
&
rootgroup
);
MPI_Comm_create
(
uworld
,
rootgroup
,
&
rootworld
);
}
delete
[]
iroots
;
// create a new compute pe style
// id = fix-ID + pe, compute group = all
int
n
=
strlen
(
id
)
+
4
;
id_pe
=
new
char
[
n
];
strcpy
(
id_pe
,
id
);
strcat
(
id_pe
,
"_pe"
);
char
**
newarg
=
new
char
*
[
3
];
newarg
[
0
]
=
id_pe
;
newarg
[
1
]
=
(
char
*
)
"all"
;
newarg
[
2
]
=
(
char
*
)
"pe"
;
modify
->
add_compute
(
3
,
newarg
);
delete
[]
newarg
;
// initialize local storage
maxlocal
=
-
1
;
ntotal
=
-
1
;
}
/* ---------------------------------------------------------------------- */
FixNEB
::~
FixNEB
()
{
modify
->
delete_compute
(
id_pe
);
delete
[]
id_pe
;
memory
->
destroy
(
xprev
);
memory
->
destroy
(
xnext
);
memory
->
destroy
(
tangent
);
memory
->
destroy
(
fnext
);
memory
->
destroy
(
springF
);
memory
->
destroy
(
xsend
);
memory
->
destroy
(
xrecv
);
memory
->
destroy
(
fsend
);
memory
->
destroy
(
frecv
);
memory
->
destroy
(
tagsend
);
memory
->
destroy
(
tagrecv
);
memory
->
destroy
(
xsendall
);
memory
->
destroy
(
xrecvall
);
memory
->
destroy
(
fsendall
);
memory
->
destroy
(
frecvall
);
memory
->
destroy
(
tagsendall
);
memory
->
destroy
(
tagrecvall
);
memory
->
destroy
(
counts
);
memory
->
destroy
(
displacements
);
if
(
NEBLongRange
)
{
if
(
rootworld
!=
MPI_COMM_NULL
)
MPI_Comm_free
(
&
rootworld
);
memory
->
destroy
(
nlenall
);
}
}
/* ---------------------------------------------------------------------- */
int
FixNEB
::
setmask
()
{
int
mask
=
0
;
mask
|=
MIN_POST_FORCE
;
return
mask
;
}
/* ---------------------------------------------------------------------- */
void
FixNEB
::
init
()
{
int
icompute
=
modify
->
find_compute
(
id_pe
);
if
(
icompute
<
0
)
error
->
all
(
FLERR
,
"Potential energy ID for fix neb does not exist"
);
pe
=
modify
->
compute
[
icompute
];
// turn off climbing mode, NEB command turns it on after init()
rclimber
=
-
1
;
// nebatoms = # of atoms in fix group = atoms with inter-replica forces
bigint
count
=
group
->
count
(
igroup
);
if
(
count
>
MAXSMALLINT
)
error
->
all
(
FLERR
,
"Too many active NEB atoms"
);
nebatoms
=
count
;
// comm mode for inter-replica exchange of coords
if
(
nreplica
==
nprocs_universe
&&
nebatoms
==
atom
->
natoms
&&
atom
->
sortfreq
==
0
)
cmode
=
SINGLE_PROC_DIRECT
;
else
if
(
nreplica
==
nprocs_universe
)
cmode
=
SINGLE_PROC_MAP
;
else
cmode
=
MULTI_PROC
;
// ntotal = total # of atoms in system, NEB atoms or not
if
(
atom
->
natoms
>
MAXSMALLINT
)
error
->
all
(
FLERR
,
"Too many atoms for NEB"
);
ntotal
=
atom
->
natoms
;
if
(
atom
->
nmax
>
maxlocal
)
reallocate
();
if
(
MULTI_PROC
&&
counts
==
NULL
)
{
memory
->
create
(
xsendall
,
ntotal
,
3
,
"neb:xsendall"
);
memory
->
create
(
xrecvall
,
ntotal
,
3
,
"neb:xrecvall"
);
memory
->
create
(
fsendall
,
ntotal
,
3
,
"neb:fsendall"
);
memory
->
create
(
frecvall
,
ntotal
,
3
,
"neb:frecvall"
);
memory
->
create
(
tagsendall
,
ntotal
,
"neb:tagsendall"
);
memory
->
create
(
tagrecvall
,
ntotal
,
"neb:tagrecvall"
);
memory
->
create
(
counts
,
nprocs
,
"neb:counts"
);
memory
->
create
(
displacements
,
nprocs
,
"neb:displacements"
);
}
}
/* ---------------------------------------------------------------------- */
void
FixNEB
::
min_setup
(
int
vflag
)
{
min_post_force
(
vflag
);
// trigger potential energy computation on next timestep
pe
->
addstep
(
update
->
ntimestep
+
1
);
}
/* ---------------------------------------------------------------------- */
void
FixNEB
::
min_post_force
(
int
vflag
)
{
double
vprev
,
vnext
;
double
delxp
,
delyp
,
delzp
,
delxn
,
delyn
,
delzn
;
double
vIni
=
0.0
;
vprev
=
vnext
=
veng
=
pe
->
compute_scalar
();
if
(
ireplica
<
nreplica
-
1
&&
me
==
0
)
MPI_Send
(
&
veng
,
1
,
MPI_DOUBLE
,
procnext
,
0
,
uworld
);
if
(
ireplica
>
0
&&
me
==
0
)
MPI_Recv
(
&
vprev
,
1
,
MPI_DOUBLE
,
procprev
,
0
,
uworld
,
MPI_STATUS_IGNORE
);
if
(
ireplica
>
0
&&
me
==
0
)
MPI_Send
(
&
veng
,
1
,
MPI_DOUBLE
,
procprev
,
0
,
uworld
);
if
(
ireplica
<
nreplica
-
1
&&
me
==
0
)
MPI_Recv
(
&
vnext
,
1
,
MPI_DOUBLE
,
procnext
,
0
,
uworld
,
MPI_STATUS_IGNORE
);
if
(
cmode
==
MULTI_PROC
)
{
MPI_Bcast
(
&
vprev
,
1
,
MPI_DOUBLE
,
0
,
world
);
MPI_Bcast
(
&
vnext
,
1
,
MPI_DOUBLE
,
0
,
world
);
}
if
(
FreeEndFinal
&&
ireplica
==
nreplica
-
1
&&
(
update
->
ntimestep
==
0
))
EFinalIni
=
veng
;
if
(
ireplica
==
0
)
vIni
=
veng
;
if
(
FreeEndFinalWithRespToEIni
)
{
if
(
me
==
0
)
{
int
procFirst
;
procFirst
=
universe
->
root_proc
[
0
];
MPI_Bcast
(
&
vIni
,
1
,
MPI_DOUBLE
,
procFirst
,
uworld
);
}
if
(
cmode
==
MULTI_PROC
)
{
MPI_Bcast
(
&
vIni
,
1
,
MPI_DOUBLE
,
0
,
world
);
}
}
if
(
FreeEndIni
&&
ireplica
==
0
&&
(
update
->
ntimestep
==
0
))
EIniIni
=
veng
;
/* if (FreeEndIni && ireplica == 0) {
// if (me == 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
pe
->
addstep
(
update
->
ntimestep
+
1
);
double
**
x
=
atom
->
x
;
int
*
mask
=
atom
->
mask
;
double
dot
=
0.0
;
double
prefactor
=
0.0
;
double
**
f
=
atom
->
f
;
int
nlocal
=
atom
->
nlocal
;
//calculating separation between images
plen
=
0.0
;
nlen
=
0.0
;
double
tlen
=
0.0
;
double
gradnextlen
=
0.0
;
dotgrad
=
gradlen
=
dotpath
=
dottangrad
=
0.0
;
if
(
ireplica
==
nreplica
-
1
)
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
{
delxp
=
x
[
i
][
0
]
-
xprev
[
i
][
0
];
delyp
=
x
[
i
][
1
]
-
xprev
[
i
][
1
];
delzp
=
x
[
i
][
2
]
-
xprev
[
i
][
2
];
domain
->
minimum_image
(
delxp
,
delyp
,
delzp
);
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
;
tangent
[
i
][
2
]
=
delzp
;
tlen
+=
tangent
[
i
][
0
]
*
tangent
[
i
][
0
]
+
tangent
[
i
][
1
]
*
tangent
[
i
][
1
]
+
tangent
[
i
][
2
]
*
tangent
[
i
][
2
];
dot
+=
f
[
i
][
0
]
*
tangent
[
i
][
0
]
+
f
[
i
][
1
]
*
tangent
[
i
][
1
]
+
f
[
i
][
2
]
*
tangent
[
i
][
2
];
}
}
}
else
if
(
ireplica
==
0
)
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
{
delxn
=
xnext
[
i
][
0
]
-
x
[
i
][
0
];
delyn
=
xnext
[
i
][
1
]
-
x
[
i
][
1
];
delzn
=
xnext
[
i
][
2
]
-
x
[
i
][
2
];
domain
->
minimum_image
(
delxn
,
delyn
,
delzn
);
nlen
+=
delxn
*
delxn
+
delyn
*
delyn
+
delzn
*
delzn
;
gradnextlen
+=
fnext
[
i
][
0
]
*
fnext
[
i
][
0
]
+
fnext
[
i
][
1
]
*
fnext
[
i
][
1
]
+
fnext
[
i
][
2
]
*
fnext
[
i
][
2
];
dotgrad
+=
f
[
i
][
0
]
*
fnext
[
i
][
0
]
+
f
[
i
][
1
]
*
fnext
[
i
][
1
]
+
f
[
i
][
2
]
*
fnext
[
i
][
2
];
dottangrad
+=
delxn
*
f
[
i
][
0
]
+
delyn
*
f
[
i
][
1
]
+
delzn
*
f
[
i
][
2
];
gradlen
+=
f
[
i
][
0
]
*
f
[
i
][
0
]
+
f
[
i
][
1
]
*
f
[
i
][
1
]
+
f
[
i
][
2
]
*
f
[
i
][
2
];
if
(
FreeEndIni
)
{
tangent
[
i
][
0
]
=
delxn
;
tangent
[
i
][
1
]
=
delyn
;
tangent
[
i
][
2
]
=
delzn
;
tlen
+=
tangent
[
i
][
0
]
*
tangent
[
i
][
0
]
+
tangent
[
i
][
1
]
*
tangent
[
i
][
1
]
+
tangent
[
i
][
2
]
*
tangent
[
i
][
2
];
dot
+=
f
[
i
][
0
]
*
tangent
[
i
][
0
]
+
f
[
i
][
1
]
*
tangent
[
i
][
1
]
+
f
[
i
][
2
]
*
tangent
[
i
][
2
];
}
}
}
else
{
// not the first or last replica
double
vmax
=
MAX
(
fabs
(
vnext
-
veng
),
fabs
(
vprev
-
veng
));
double
vmin
=
MIN
(
fabs
(
vnext
-
veng
),
fabs
(
vprev
-
veng
));
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
{
delxp
=
x
[
i
][
0
]
-
xprev
[
i
][
0
];
delyp
=
x
[
i
][
1
]
-
xprev
[
i
][
1
];
delzp
=
x
[
i
][
2
]
-
xprev
[
i
][
2
];
domain
->
minimum_image
(
delxp
,
delyp
,
delzp
);
plen
+=
delxp
*
delxp
+
delyp
*
delyp
+
delzp
*
delzp
;
delxn
=
xnext
[
i
][
0
]
-
x
[
i
][
0
];
delyn
=
xnext
[
i
][
1
]
-
x
[
i
][
1
];
delzn
=
xnext
[
i
][
2
]
-
x
[
i
][
2
];
domain
->
minimum_image
(
delxn
,
delyn
,
delzn
);
if
(
vnext
>
veng
&&
veng
>
vprev
)
{
tangent
[
i
][
0
]
=
delxn
;
tangent
[
i
][
1
]
=
delyn
;
tangent
[
i
][
2
]
=
delzn
;
}
else
if
(
vnext
<
veng
&&
veng
<
vprev
)
{
tangent
[
i
][
0
]
=
delxp
;
tangent
[
i
][
1
]
=
delyp
;
tangent
[
i
][
2
]
=
delzp
;
}
else
{
if
(
vnext
>
vprev
)
{
tangent
[
i
][
0
]
=
vmax
*
delxn
+
vmin
*
delxp
;
tangent
[
i
][
1
]
=
vmax
*
delyn
+
vmin
*
delyp
;
tangent
[
i
][
2
]
=
vmax
*
delzn
+
vmin
*
delzp
;
}
else
{
tangent
[
i
][
0
]
=
vmin
*
delxn
+
vmax
*
delxp
;
tangent
[
i
][
1
]
=
vmin
*
delyn
+
vmax
*
delyp
;
tangent
[
i
][
2
]
=
vmin
*
delzn
+
vmax
*
delzp
;
}
}
nlen
+=
delxn
*
delxn
+
delyn
*
delyn
+
delzn
*
delzn
;
tlen
+=
tangent
[
i
][
0
]
*
tangent
[
i
][
0
]
+
tangent
[
i
][
1
]
*
tangent
[
i
][
1
]
+
tangent
[
i
][
2
]
*
tangent
[
i
][
2
];
gradlen
+=
f
[
i
][
0
]
*
f
[
i
][
0
]
+
f
[
i
][
1
]
*
f
[
i
][
1
]
+
f
[
i
][
2
]
*
f
[
i
][
2
];
dotpath
+=
delxp
*
delxn
+
delyp
*
delyn
+
delzp
*
delzn
;
dottangrad
+=
tangent
[
i
][
0
]
*
f
[
i
][
0
]
+
tangent
[
i
][
1
]
*
f
[
i
][
1
]
+
tangent
[
i
][
2
]
*
f
[
i
][
2
];
gradnextlen
+=
fnext
[
i
][
0
]
*
fnext
[
i
][
0
]
+
fnext
[
i
][
1
]
*
fnext
[
i
][
1
]
+
fnext
[
i
][
2
]
*
fnext
[
i
][
2
];
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
bufin
[
BUFSIZE
],
bufout
[
BUFSIZE
];
bufin
[
0
]
=
nlen
;
bufin
[
1
]
=
plen
;
bufin
[
2
]
=
tlen
;
bufin
[
3
]
=
gradlen
;
bufin
[
4
]
=
gradnextlen
;
bufin
[
5
]
=
dotpath
;
bufin
[
6
]
=
dottangrad
;
bufin
[
7
]
=
dotgrad
;
MPI_Allreduce
(
bufin
,
bufout
,
BUFSIZE
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
nlen
=
sqrt
(
bufout
[
0
]);
plen
=
sqrt
(
bufout
[
1
]);
tlen
=
sqrt
(
bufout
[
2
]);
gradlen
=
sqrt
(
bufout
[
3
]);
gradnextlen
=
sqrt
(
bufout
[
4
]);
dotpath
=
bufout
[
5
];
dottangrad
=
bufout
[
6
];
dotgrad
=
bufout
[
7
];
// normalize tangent vector
if
(
tlen
>
0.0
)
{
double
tleninv
=
1.0
/
tlen
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
{
tangent
[
i
][
0
]
*=
tleninv
;
tangent
[
i
][
1
]
*=
tleninv
;
tangent
[
i
][
2
]
*=
tleninv
;
}
}
// first or last replica has no change to forces, just return
if
(
ireplica
>
0
&&
ireplica
<
nreplica
-
1
)
dottangrad
=
dottangrad
/
(
tlen
*
gradlen
);
if
(
ireplica
==
0
)
dottangrad
=
dottangrad
/
(
nlen
*
gradlen
);
if
(
ireplica
==
nreplica
-
1
)
dottangrad
=
dottangrad
/
(
plen
*
gradlen
);
if
(
ireplica
<
nreplica
-
1
)
dotgrad
=
dotgrad
/
(
gradlen
*
gradnextlen
);
if
(
FreeEndIni
&&
ireplica
==
0
)
{
if
(
tlen
>
0.0
)
{
double
dotall
;
MPI_Allreduce
(
&
dot
,
&
dotall
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
dot
=
dotall
/
tlen
;
if
(
dot
<
0
)
prefactor
=
-
dot
-
kspringIni
*
(
veng
-
EIniIni
);
else
prefactor
=
-
dot
+
kspringIni
*
(
veng
-
EIniIni
);
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
{
f
[
i
][
0
]
+=
prefactor
*
tangent
[
i
][
0
];
f
[
i
][
1
]
+=
prefactor
*
tangent
[
i
][
1
];
f
[
i
][
2
]
+=
prefactor
*
tangent
[
i
][
2
];
}
}
}
if
(
FreeEndFinal
&&
ireplica
==
nreplica
-
1
)
{
if
(
tlen
>
0.0
)
{
double
dotall
;
MPI_Allreduce
(
&
dot
,
&
dotall
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
dot
=
dotall
/
tlen
;
if
(
dot
<
0
)
prefactor
=
-
dot
-
kspringFinal
*
(
veng
-
EFinalIni
);
else
prefactor
=
-
dot
+
kspringFinal
*
(
veng
-
EFinalIni
);
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
{
f
[
i
][
0
]
+=
prefactor
*
tangent
[
i
][
0
];
f
[
i
][
1
]
+=
prefactor
*
tangent
[
i
][
1
];
f
[
i
][
2
]
+=
prefactor
*
tangent
[
i
][
2
];
}
}
}
if
(
FreeEndFinalWithRespToEIni
&&
ireplica
==
nreplica
-
1
)
{
if
(
tlen
>
0.0
)
{
double
dotall
;
MPI_Allreduce
(
&
dot
,
&
dotall
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
dot
=
dotall
/
tlen
;
if
(
dot
<
0
)
prefactor
=
-
dot
-
kspringFinal
*
(
veng
-
vIni
);
else
prefactor
=
-
dot
+
kspringFinal
*
(
veng
-
vIni
);
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
{
f
[
i
][
0
]
+=
prefactor
*
tangent
[
i
][
0
];
f
[
i
][
1
]
+=
prefactor
*
tangent
[
i
][
1
];
f
[
i
][
2
]
+=
prefactor
*
tangent
[
i
][
2
];
}
}
}
double
lentot
=
0
;
double
meanDist
,
idealPos
,
lenuntilIm
,
lenuntilClimber
;
lenuntilClimber
=
0
;
if
(
NEBLongRange
)
{
if
(
cmode
==
SINGLE_PROC_DIRECT
||
cmode
==
SINGLE_PROC_MAP
)
{
MPI_Allgather
(
&
nlen
,
1
,
MPI_DOUBLE
,
&
nlenall
[
0
],
1
,
MPI_DOUBLE
,
uworld
);
}
else
{
if
(
me
==
0
)
MPI_Allgather
(
&
nlen
,
1
,
MPI_DOUBLE
,
&
nlenall
[
0
],
1
,
MPI_DOUBLE
,
rootworld
);
MPI_Bcast
(
nlenall
,
nreplica
,
MPI_DOUBLE
,
0
,
world
);
}
lenuntilIm
=
0
;
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
)
{
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
;
dotpath
=
dotpath
/
(
plen
*
nlen
);
AngularContr
=
0.5
*
(
1
+
cos
(
MY_PI
*
dotpath
));
double
dotSpringTangent
;
dotSpringTangent
=
0
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
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
];}
}
double
dotSpringTangentall
;
MPI_Allreduce
(
&
dotSpringTangent
,
&
dotSpringTangentall
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
dotSpringTangent
=
dotSpringTangentall
;
double
dotall
;
MPI_Allreduce
(
&
dot
,
&
dotall
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
dot
=
dotall
;
if
(
ireplica
==
rclimber
)
prefactor
=
-
2.0
*
dot
;
else
{
if
(
NEBLongRange
)
{
prefactor
=
-
dot
-
kspring
*
(
lenuntilIm
-
idealPos
)
/
(
2
*
meanDist
);
}
else
if
(
StandardNEB
)
{
prefactor
=
-
dot
+
kspring
*
(
nlen
-
plen
);
}
if
(
FinalAndInterWithRespToEIni
&&
veng
<
vIni
)
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
{
f
[
i
][
0
]
=
0
;
f
[
i
][
1
]
=
0
;
f
[
i
][
2
]
=
0
;
}
prefactor
=
kspring
*
(
nlen
-
plen
);
AngularContr
=
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
]);
}
}
/* ----------------------------------------------------------------------
send/recv NEB atoms to/from adjacent replicas
received atoms matching my local atoms are stored in xprev,xnext
replicas 0 and N-1 send but do not receive any atoms
------------------------------------------------------------------------- */
void
FixNEB
::
inter_replica_comm
()
{
int
i
,
m
;
MPI_Request
request
;
MPI_Request
requests
[
2
];
MPI_Status
statuses
[
2
];
// reallocate memory if necessary
if
(
atom
->
nmax
>
maxlocal
)
reallocate
();
double
**
x
=
atom
->
x
;
double
**
f
=
atom
->
f
;
tagint
*
tag
=
atom
->
tag
;
int
*
mask
=
atom
->
mask
;
int
nlocal
=
atom
->
nlocal
;
// -----------------------------------------------------
// 3 cases: two for single proc per replica
// one for multiple procs per replica
// -----------------------------------------------------
// single proc per replica
// all atoms are NEB atoms and no atom sorting
// direct comm of x -> xprev and x -> xnext
if
(
cmode
==
SINGLE_PROC_DIRECT
)
{
if
(
ireplica
>
0
)
MPI_Irecv
(
xprev
[
0
],
3
*
nlocal
,
MPI_DOUBLE
,
procprev
,
0
,
uworld
,
&
request
);
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
)
MPI_Send
(
x
[
0
],
3
*
nlocal
,
MPI_DOUBLE
,
procprev
,
0
,
uworld
);
if
(
ireplica
<
nreplica
-
1
)
MPI_Wait
(
&
request
,
MPI_STATUS_IGNORE
);
if
(
ireplica
<
nreplica
-
1
)
MPI_Irecv
(
fnext
[
0
],
3
*
nlocal
,
MPI_DOUBLE
,
procnext
,
0
,
uworld
,
&
request
);
if
(
ireplica
>
0
)
MPI_Send
(
f
[
0
],
3
*
nlocal
,
MPI_DOUBLE
,
procprev
,
0
,
uworld
);
if
(
ireplica
<
nreplica
-
1
)
MPI_Wait
(
&
request
,
MPI_STATUS_IGNORE
);
return
;
}
// single proc per replica
// but only some atoms are NEB atoms or atom sorting is enabled
// send atom IDs and coords of only NEB atoms to prev/next proc
// recv procs use atom->map() to match received coords to owned atoms
if
(
cmode
==
SINGLE_PROC_MAP
)
{
m
=
0
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
{
tagsend
[
m
]
=
tag
[
i
];
xsend
[
m
][
0
]
=
x
[
i
][
0
];
xsend
[
m
][
1
]
=
x
[
i
][
1
];
xsend
[
m
][
2
]
=
x
[
i
][
2
];
fsend
[
m
][
0
]
=
f
[
i
][
0
];
fsend
[
m
][
1
]
=
f
[
i
][
1
];
fsend
[
m
][
2
]
=
f
[
i
][
2
];
m
++
;
}
if
(
ireplica
>
0
)
{
MPI_Irecv
(
xrecv
[
0
],
3
*
nebatoms
,
MPI_DOUBLE
,
procprev
,
0
,
uworld
,
&
requests
[
0
]);
MPI_Irecv
(
tagrecv
,
nebatoms
,
MPI_LMP_TAGINT
,
procprev
,
0
,
uworld
,
&
requests
[
1
]);
}
if
(
ireplica
<
nreplica
-
1
)
{
MPI_Send
(
xsend
[
0
],
3
*
nebatoms
,
MPI_DOUBLE
,
procnext
,
0
,
uworld
);
MPI_Send
(
tagsend
,
nebatoms
,
MPI_LMP_TAGINT
,
procnext
,
0
,
uworld
);
}
if
(
ireplica
>
0
)
{
MPI_Waitall
(
2
,
requests
,
statuses
);
for
(
i
=
0
;
i
<
nebatoms
;
i
++
)
{
m
=
atom
->
map
(
tagrecv
[
i
]);
xprev
[
m
][
0
]
=
xrecv
[
i
][
0
];
xprev
[
m
][
1
]
=
xrecv
[
i
][
1
];
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
]);
MPI_Irecv
(
tagrecv
,
nebatoms
,
MPI_LMP_TAGINT
,
procnext
,
0
,
uworld
,
&
requests
[
1
]);
}
if
(
ireplica
>
0
)
{
MPI_Send
(
xsend
[
0
],
3
*
nebatoms
,
MPI_DOUBLE
,
procprev
,
0
,
uworld
);
MPI_Send
(
fsend
[
0
],
3
*
nebatoms
,
MPI_DOUBLE
,
procprev
,
0
,
uworld
);
MPI_Send
(
tagsend
,
nebatoms
,
MPI_LMP_TAGINT
,
procprev
,
0
,
uworld
);
}
if
(
ireplica
<
nreplica
-
1
)
{
MPI_Waitall
(
2
,
requests
,
statuses
);
for
(
i
=
0
;
i
<
nebatoms
;
i
++
)
{
m
=
atom
->
map
(
tagrecv
[
i
]);
xnext
[
m
][
0
]
=
xrecv
[
i
][
0
];
xnext
[
m
][
1
]
=
xrecv
[
i
][
1
];
xnext
[
m
][
2
]
=
xrecv
[
i
][
2
];
fnext
[
m
][
0
]
=
frecv
[
i
][
0
];
fnext
[
m
][
1
]
=
frecv
[
i
][
1
];
fnext
[
m
][
2
]
=
frecv
[
i
][
2
];
}
}
return
;
}
// multiple procs per replica
// MPI_Gather all coords and atom IDs to root proc of each replica
// send to root of adjacent replicas
// bcast within each replica
// each proc extracts info for atoms it owns via atom->map()
m
=
0
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
)
{
tagsend
[
m
]
=
tag
[
i
];
xsend
[
m
][
0
]
=
x
[
i
][
0
];
xsend
[
m
][
1
]
=
x
[
i
][
1
];
xsend
[
m
][
2
]
=
x
[
i
][
2
];
fsend
[
m
][
0
]
=
f
[
i
][
0
];
fsend
[
m
][
1
]
=
f
[
i
][
1
];
fsend
[
m
][
2
]
=
f
[
i
][
2
];
m
++
;
}
MPI_Gather
(
&
m
,
1
,
MPI_INT
,
counts
,
1
,
MPI_INT
,
0
,
world
);
displacements
[
0
]
=
0
;
for
(
i
=
0
;
i
<
nprocs
-
1
;
i
++
)
displacements
[
i
+
1
]
=
displacements
[
i
]
+
counts
[
i
];
MPI_Gatherv
(
tagsend
,
m
,
MPI_LMP_TAGINT
,
tagsendall
,
counts
,
displacements
,
MPI_LMP_TAGINT
,
0
,
world
);
for
(
i
=
0
;
i
<
nprocs
;
i
++
)
counts
[
i
]
*=
3
;
for
(
i
=
0
;
i
<
nprocs
-
1
;
i
++
)
displacements
[
i
+
1
]
=
displacements
[
i
]
+
counts
[
i
];
if
(
xsend
)
{
MPI_Gatherv
(
xsend
[
0
],
3
*
m
,
MPI_DOUBLE
,
xsendall
[
0
],
counts
,
displacements
,
MPI_DOUBLE
,
0
,
world
);
MPI_Gatherv
(
fsend
[
0
],
3
*
m
,
MPI_DOUBLE
,
fsendall
[
0
],
counts
,
displacements
,
MPI_DOUBLE
,
0
,
world
);
}
else
{
MPI_Gatherv
(
NULL
,
3
*
m
,
MPI_DOUBLE
,
xsendall
[
0
],
counts
,
displacements
,
MPI_DOUBLE
,
0
,
world
);
MPI_Gatherv
(
NULL
,
3
*
m
,
MPI_DOUBLE
,
fsendall
[
0
],
counts
,
displacements
,
MPI_DOUBLE
,
0
,
world
);
}
if
(
ireplica
>
0
&&
me
==
0
)
{
MPI_Irecv
(
xrecvall
[
0
],
3
*
nebatoms
,
MPI_DOUBLE
,
procprev
,
0
,
uworld
,
&
requests
[
0
]);
MPI_Irecv
(
tagrecvall
,
nebatoms
,
MPI_LMP_TAGINT
,
procprev
,
0
,
uworld
,
&
requests
[
1
]);
}
if
(
ireplica
<
nreplica
-
1
&&
me
==
0
)
{
MPI_Send
(
xsendall
[
0
],
3
*
nebatoms
,
MPI_DOUBLE
,
procnext
,
0
,
uworld
);
MPI_Send
(
tagsendall
,
nebatoms
,
MPI_LMP_TAGINT
,
procnext
,
0
,
uworld
);
}
if
(
ireplica
>
0
)
{
if
(
me
==
0
)
MPI_Waitall
(
2
,
requests
,
statuses
);
MPI_Bcast
(
tagrecvall
,
nebatoms
,
MPI_INT
,
0
,
world
);
MPI_Bcast
(
xrecvall
[
0
],
3
*
nebatoms
,
MPI_DOUBLE
,
0
,
world
);
for
(
i
=
0
;
i
<
nebatoms
;
i
++
)
{
m
=
atom
->
map
(
tagrecvall
[
i
]);
if
(
m
<
0
||
m
>=
nlocal
)
continue
;
xprev
[
m
][
0
]
=
xrecvall
[
i
][
0
];
xprev
[
m
][
1
]
=
xrecvall
[
i
][
1
];
xprev
[
m
][
2
]
=
xrecvall
[
i
][
2
];
}
}
if
(
ireplica
<
nreplica
-
1
&&
me
==
0
)
{
MPI_Irecv
(
xrecvall
[
0
],
3
*
nebatoms
,
MPI_DOUBLE
,
procnext
,
0
,
uworld
,
&
requests
[
0
]);
MPI_Irecv
(
frecvall
[
0
],
3
*
nebatoms
,
MPI_DOUBLE
,
procnext
,
0
,
uworld
,
&
requests
[
0
]);
MPI_Irecv
(
tagrecvall
,
nebatoms
,
MPI_LMP_TAGINT
,
procnext
,
0
,
uworld
,
&
requests
[
1
]);
}
if
(
ireplica
>
0
&&
me
==
0
)
{
MPI_Send
(
xsendall
[
0
],
3
*
nebatoms
,
MPI_DOUBLE
,
procprev
,
0
,
uworld
);
MPI_Send
(
fsendall
[
0
],
3
*
nebatoms
,
MPI_DOUBLE
,
procprev
,
0
,
uworld
);
MPI_Send
(
tagsendall
,
nebatoms
,
MPI_LMP_TAGINT
,
procprev
,
0
,
uworld
);
}
if
(
ireplica
<
nreplica
-
1
)
{
if
(
me
==
0
)
MPI_Waitall
(
2
,
requests
,
statuses
);
MPI_Bcast
(
tagrecvall
,
nebatoms
,
MPI_INT
,
0
,
world
);
MPI_Bcast
(
xrecvall
[
0
],
3
*
nebatoms
,
MPI_DOUBLE
,
0
,
world
);
MPI_Bcast
(
frecvall
[
0
],
3
*
nebatoms
,
MPI_DOUBLE
,
0
,
world
);
for
(
i
=
0
;
i
<
nebatoms
;
i
++
)
{
m
=
atom
->
map
(
tagrecvall
[
i
]);
if
(
m
<
0
||
m
>=
nlocal
)
continue
;
xnext
[
m
][
0
]
=
xrecvall
[
i
][
0
];
xnext
[
m
][
1
]
=
xrecvall
[
i
][
1
];
xnext
[
m
][
2
]
=
xrecvall
[
i
][
2
];
fnext
[
m
][
0
]
=
frecvall
[
i
][
0
];
fnext
[
m
][
1
]
=
frecvall
[
i
][
1
];
fnext
[
m
][
2
]
=
frecvall
[
i
][
2
];
}
}
}
/* ----------------------------------------------------------------------
reallocate xprev,xnext,tangent arrays if necessary
reallocate communication arrays if necessary
------------------------------------------------------------------------- */
void
FixNEB
::
reallocate
()
{
maxlocal
=
atom
->
nmax
;
memory
->
destroy
(
xprev
);
memory
->
destroy
(
xnext
);
memory
->
destroy
(
tangent
);
memory
->
destroy
(
fnext
);
memory
->
destroy
(
springF
);
memory
->
create
(
xprev
,
maxlocal
,
3
,
"neb:xprev"
);
memory
->
create
(
xnext
,
maxlocal
,
3
,
"neb:xnext"
);
memory
->
create
(
tangent
,
maxlocal
,
3
,
"neb:tangent"
);
memory
->
create
(
fnext
,
maxlocal
,
3
,
"neb:fnext"
);
memory
->
create
(
springF
,
maxlocal
,
3
,
"neb:springF"
);
if
(
cmode
!=
SINGLE_PROC_DIRECT
)
{
memory
->
destroy
(
xsend
);
memory
->
destroy
(
fsend
);
memory
->
destroy
(
xrecv
);
memory
->
destroy
(
frecv
);
memory
->
destroy
(
tagsend
);
memory
->
destroy
(
tagrecv
);
memory
->
create
(
xsend
,
maxlocal
,
3
,
"neb:xsend"
);
memory
->
create
(
fsend
,
maxlocal
,
3
,
"neb:fsend"
);
memory
->
create
(
xrecv
,
maxlocal
,
3
,
"neb:xrecv"
);
memory
->
create
(
frecv
,
maxlocal
,
3
,
"neb:frecv"
);
memory
->
create
(
tagsend
,
maxlocal
,
"neb:tagsend"
);
memory
->
create
(
tagrecv
,
maxlocal
,
"neb:tagrecv"
);
}
if
(
NEBLongRange
)
{
memory
->
destroy
(
nlenall
);
memory
->
create
(
nlenall
,
nreplica
,
"neb:nlenall"
);
}
}
Event Timeline
Log In to Comment