Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F85863347
fix_shake.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
Wed, Oct 2, 15:13
Size
70 KB
Mime Type
text/x-c
Expires
Fri, Oct 4, 15:13 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
21282983
Attached To
rLAMMPS lammps
fix_shake.cpp
View Options
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
www.cs.sandia.gov/~sjplimp/lammps.html
Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "mpi.h"
#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "stdio.h"
#include "fix_shake.h"
#include "atom.h"
#include "update.h"
#include "respa.h"
#include "modify.h"
#include "domain.h"
#include "force.h"
#include "bond.h"
#include "angle.h"
#include "comm.h"
#include "group.h"
#include "fix_respa.h"
#include "memory.h"
#include "error.h"
#define BIG 1.0e20
#define MASSDELTA 0.1
#define MIN(a,b) ((a) < (b) ? (a) : (b))
#define MAX(a,b) ((a) > (b) ? (a) : (b))
/* ---------------------------------------------------------------------- */
FixShake
::
FixShake
(
int
narg
,
char
**
arg
)
:
Fix
(
narg
,
arg
)
{
MPI_Comm_rank
(
world
,
&
me
);
MPI_Comm_size
(
world
,
&
nprocs
);
PI
=
4.0
*
atan
(
1.0
);
// error check
if
(
atom
->
molecular
==
0
)
error
->
all
(
"Cannot use fix shake with non-molecular system"
);
// perform initial allocation of atom-based arrays
// register with atom class
shake_flag
=
NULL
;
shake_atom
=
shake_type
=
NULL
;
xshake
=
NULL
;
grow_arrays
(
atom
->
nmax
);
atom
->
add_callback
(
0
);
// parse SHAKE args
if
(
narg
<
8
)
error
->
all
(
"Illegal fix shake command"
);
tolerance
=
atof
(
arg
[
3
]);
max_iter
=
atoi
(
arg
[
4
]);
output_every
=
atoi
(
arg
[
5
]);
// parse SHAKE args for bond and angle types
// will be used by find_clusters
// store args for "b" "a" "t" as flags in (1:n) list for fast access
// store args for "m" in list of length nmass for looping over
// for "m" verify that atom masses have been set
bond_flag
=
new
int
[
atom
->
nbondtypes
+
1
];
for
(
int
i
=
1
;
i
<=
atom
->
nbondtypes
;
i
++
)
bond_flag
[
i
]
=
0
;
angle_flag
=
new
int
[
atom
->
nangletypes
+
1
];
for
(
int
i
=
1
;
i
<=
atom
->
nangletypes
;
i
++
)
angle_flag
[
i
]
=
0
;
type_flag
=
new
int
[
atom
->
ntypes
+
1
];
for
(
int
i
=
1
;
i
<=
atom
->
ntypes
;
i
++
)
type_flag
[
i
]
=
0
;
mass_list
=
new
double
[
atom
->
ntypes
];
nmass
=
0
;
char
mode
=
'\0'
;
int
next
=
6
;
while
(
next
<
narg
)
{
if
(
strcmp
(
arg
[
next
],
"b"
)
==
0
)
mode
=
'b'
;
else
if
(
strcmp
(
arg
[
next
],
"a"
)
==
0
)
mode
=
'a'
;
else
if
(
strcmp
(
arg
[
next
],
"t"
)
==
0
)
mode
=
't'
;
else
if
(
strcmp
(
arg
[
next
],
"m"
)
==
0
)
{
mode
=
'm'
;
atom
->
check_mass
();
}
else
if
(
mode
==
'b'
)
{
int
i
=
atoi
(
arg
[
next
]);
if
(
i
<
1
||
i
>
atom
->
nbondtypes
)
error
->
all
(
"Invalid bond type index for fix shake"
);
bond_flag
[
i
]
=
1
;
}
else
if
(
mode
==
'a'
)
{
int
i
=
atoi
(
arg
[
next
]);
if
(
i
<
1
||
i
>
atom
->
nangletypes
)
error
->
all
(
"Invalid angle type index for fix shake"
);
angle_flag
[
i
]
=
1
;
}
else
if
(
mode
==
't'
)
{
int
i
=
atoi
(
arg
[
next
]);
if
(
i
<
1
||
i
>
atom
->
ntypes
)
error
->
all
(
"Invalid atom type index for fix shake"
);
type_flag
[
i
]
=
1
;
}
else
if
(
mode
==
'm'
)
{
double
rmass
=
atof
(
arg
[
next
]);
if
(
rmass
==
0.0
)
error
->
all
(
"Invalid atom mass for fix shake"
);
if
(
nmass
==
atom
->
ntypes
)
error
->
all
(
"Too many masses for fix shake"
);
mass_list
[
nmass
++
]
=
rmass
;
}
else
error
->
all
(
"Illegal fix shake command"
);
next
++
;
}
// allocate bond and angle distance arrays, indexed from 1 to n
bond_distance
=
new
double
[
atom
->
nbondtypes
+
1
];
angle_distance
=
new
double
[
atom
->
nangletypes
+
1
];
// allocate statistics arrays
if
(
output_every
)
{
int
nb
=
atom
->
nbondtypes
+
1
;
b_count
=
new
int
[
nb
];
b_count_all
=
new
int
[
nb
];
b_ave
=
new
double
[
nb
];
b_ave_all
=
new
double
[
nb
];
b_max
=
new
double
[
nb
];
b_max_all
=
new
double
[
nb
];
b_min
=
new
double
[
nb
];
b_min_all
=
new
double
[
nb
];
int
na
=
atom
->
nangletypes
+
1
;
a_count
=
new
int
[
na
];
a_count_all
=
new
int
[
na
];
a_ave
=
new
double
[
na
];
a_ave_all
=
new
double
[
na
];
a_max
=
new
double
[
na
];
a_max_all
=
new
double
[
na
];
a_min
=
new
double
[
na
];
a_min_all
=
new
double
[
na
];
}
// identify all SHAKE clusters
find_clusters
();
// initialize list of SHAKE clusters to constrain
maxlist
=
0
;
list
=
NULL
;
}
/* ---------------------------------------------------------------------- */
FixShake
::~
FixShake
()
{
// if atom class still exists:
// unregister this fix so atom class doesn't invoke it any more
// set bond_type and angle_type back to positive for SHAKE clusters
// must set for all SHAKE bonds and angles stored by each atom
if
(
atom
)
{
atom
->
delete_callback
(
id
,
0
);
int
**
bond_type
=
atom
->
bond_type
;
int
**
angle_type
=
atom
->
angle_type
;
int
nlocal
=
atom
->
nlocal
;
int
n
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
shake_flag
[
i
]
==
0
)
continue
;
else
if
(
shake_flag
[
i
]
==
1
)
{
n
=
bondfind
(
i
,
shake_atom
[
i
][
0
],
shake_atom
[
i
][
1
]);
if
(
n
>=
0
)
bond_type
[
i
][
n
]
=
-
bond_type
[
i
][
n
];
n
=
bondfind
(
i
,
shake_atom
[
i
][
0
],
shake_atom
[
i
][
2
]);
if
(
n
>=
0
)
bond_type
[
i
][
n
]
=
-
bond_type
[
i
][
n
];
n
=
anglefind
(
i
,
shake_atom
[
i
][
1
],
shake_atom
[
i
][
2
]);
if
(
n
>=
0
)
angle_type
[
i
][
n
]
=
-
angle_type
[
i
][
n
];
}
else
if
(
shake_flag
[
i
]
==
2
)
{
n
=
bondfind
(
i
,
shake_atom
[
i
][
0
],
shake_atom
[
i
][
1
]);
if
(
n
>=
0
)
bond_type
[
i
][
n
]
=
-
bond_type
[
i
][
n
];
}
else
if
(
shake_flag
[
i
]
==
3
)
{
n
=
bondfind
(
i
,
shake_atom
[
i
][
0
],
shake_atom
[
i
][
1
]);
if
(
n
>=
0
)
bond_type
[
i
][
n
]
=
-
bond_type
[
i
][
n
];
n
=
bondfind
(
i
,
shake_atom
[
i
][
0
],
shake_atom
[
i
][
2
]);
if
(
n
>=
0
)
bond_type
[
i
][
n
]
=
-
bond_type
[
i
][
n
];
}
else
if
(
shake_flag
[
i
]
==
4
)
{
n
=
bondfind
(
i
,
shake_atom
[
i
][
0
],
shake_atom
[
i
][
1
]);
if
(
n
>=
0
)
bond_type
[
i
][
n
]
=
-
bond_type
[
i
][
n
];
n
=
bondfind
(
i
,
shake_atom
[
i
][
0
],
shake_atom
[
i
][
2
]);
if
(
n
>=
0
)
bond_type
[
i
][
n
]
=
-
bond_type
[
i
][
n
];
n
=
bondfind
(
i
,
shake_atom
[
i
][
0
],
shake_atom
[
i
][
3
]);
if
(
n
>=
0
)
bond_type
[
i
][
n
]
=
-
bond_type
[
i
][
n
];
}
}
}
// delete locally stored arrays
memory
->
sfree
(
shake_flag
);
memory
->
destroy_2d_int_array
(
shake_atom
);
memory
->
destroy_2d_int_array
(
shake_type
);
memory
->
destroy_2d_double_array
(
xshake
);
delete
[]
bond_flag
;
delete
[]
angle_flag
;
delete
[]
type_flag
;
delete
[]
mass_list
;
delete
[]
bond_distance
;
delete
[]
angle_distance
;
if
(
output_every
)
{
delete
[]
b_count
;
delete
[]
b_count_all
;
delete
[]
b_ave
;
delete
[]
b_ave_all
;
delete
[]
b_max
;
delete
[]
b_max_all
;
delete
[]
b_min
;
delete
[]
b_min_all
;
delete
[]
a_count
;
delete
[]
a_count_all
;
delete
[]
a_ave
;
delete
[]
a_ave_all
;
delete
[]
a_max
;
delete
[]
a_max_all
;
delete
[]
a_min
;
delete
[]
a_min_all
;
}
memory
->
sfree
(
list
);
}
/* ---------------------------------------------------------------------- */
int
FixShake
::
setmask
()
{
int
mask
=
0
;
mask
|=
PRE_NEIGHBOR
;
mask
|=
POST_FORCE
;
mask
|=
POST_FORCE_RESPA
;
return
mask
;
}
/* ----------------------------------------------------------------------
set bond and angle distances
this init must happen after force->bond and force->angle inits
------------------------------------------------------------------------- */
void
FixShake
::
init
()
{
int
i
,
m
,
flag
,
flag_all
,
type1
,
type2
,
bond1_type
,
bond2_type
;
double
rsq
,
angle
;
// error if more than one shake fix
int
count
=
0
;
for
(
i
=
0
;
i
<
modify
->
nfix
;
i
++
)
if
(
strcmp
(
modify
->
fix
[
i
]
->
style
,
"shake"
)
==
0
)
count
++
;
if
(
count
>
1
)
error
->
all
(
"More than one shake fix"
);
// error if npt,nph fix comes before shake fix
for
(
i
=
0
;
i
<
modify
->
nfix
;
i
++
)
{
if
(
strcmp
(
modify
->
fix
[
i
]
->
style
,
"npt"
)
==
0
)
break
;
if
(
strcmp
(
modify
->
fix
[
i
]
->
style
,
"nph"
)
==
0
)
break
;
}
if
(
i
<
modify
->
nfix
)
{
for
(
int
j
=
i
;
j
<
modify
->
nfix
;
j
++
)
if
(
strcmp
(
modify
->
fix
[
j
]
->
style
,
"shake"
)
==
0
)
error
->
all
(
"Shake fix must come before NPT/NPH fix"
);
}
// if rRESPA, find associated fix that must exist
// could have changed locations in fix list since created
// set ptrs to rRESPA variables
if
(
strcmp
(
update
->
integrate_style
,
"respa"
)
==
0
)
{
for
(
i
=
0
;
i
<
modify
->
nfix
;
i
++
)
if
(
strcmp
(
modify
->
fix
[
i
]
->
style
,
"RESPA"
)
==
0
)
ifix_respa
=
i
;
nlevels_respa
=
((
Respa
*
)
update
->
integrate
)
->
nlevels
;
loop_respa
=
((
Respa
*
)
update
->
integrate
)
->
loop
;
step_respa
=
((
Respa
*
)
update
->
integrate
)
->
step
;
}
// set communication size in comm class
comm
->
maxforward_fix
=
MAX
(
comm
->
maxforward_fix
,
3
);
// set equilibrium bond distances
if
(
force
->
bond
==
NULL
)
error
->
all
(
"Bond potential must be defined for SHAKE"
);
for
(
i
=
1
;
i
<=
atom
->
nbondtypes
;
i
++
)
bond_distance
[
i
]
=
force
->
bond
->
equilibrium_distance
(
i
);
// set equilibrium angle distances
int
nlocal
=
atom
->
nlocal
;
for
(
i
=
1
;
i
<=
atom
->
nangletypes
;
i
++
)
{
if
(
angle_flag
[
i
]
==
0
)
continue
;
// scan all atoms for a SHAKE angle cluster
// extract bond types for the 2 bonds in the cluster
// bond types must be same in all clusters of this angle type,
// else set error flag
flag
=
0
;
bond1_type
=
bond2_type
=
0
;
for
(
m
=
0
;
m
<
nlocal
;
m
++
)
{
if
(
shake_flag
[
m
]
!=
1
)
continue
;
if
(
shake_type
[
m
][
2
]
!=
i
)
continue
;
type1
=
MIN
(
shake_type
[
m
][
0
],
shake_type
[
m
][
1
]);
type2
=
MAX
(
shake_type
[
m
][
0
],
shake_type
[
m
][
1
]);
if
(
bond1_type
>
0
)
{
if
(
type1
!=
bond1_type
||
type2
!=
bond2_type
)
{
flag
=
1
;
break
;
}
}
bond1_type
=
type1
;
bond2_type
=
type2
;
}
// error check for any bond types that are not the same
MPI_Allreduce
(
&
flag
,
&
flag_all
,
1
,
MPI_INT
,
MPI_MAX
,
world
);
if
(
flag_all
)
error
->
all
(
"Shake angles have different bond types"
);
// insure all procs have bond types
MPI_Allreduce
(
&
bond1_type
,
&
flag_all
,
1
,
MPI_INT
,
MPI_MAX
,
world
);
bond1_type
=
flag_all
;
MPI_Allreduce
(
&
bond2_type
,
&
flag_all
,
1
,
MPI_INT
,
MPI_MAX
,
world
);
bond2_type
=
flag_all
;
// if bond types are 0, no SHAKE angles of this type exist
// just skip this angle
if
(
bond1_type
==
0
)
{
angle_distance
[
i
]
=
0.0
;
continue
;
}
// compute the angle distance as a function of 2 bond distances
angle
=
force
->
angle
->
equilibrium_angle
(
i
);
rsq
=
2.0
*
bond_distance
[
bond1_type
]
*
bond_distance
[
bond2_type
]
*
(
1.0
-
cos
(
angle
));
angle_distance
[
i
]
=
sqrt
(
rsq
);
}
}
/* ----------------------------------------------------------------------
SHAKE as pre-integrator constraint
------------------------------------------------------------------------- */
void
FixShake
::
setup
()
{
pre_neighbor
();
if
(
output_every
)
stats
();
// setup SHAKE output
int
ntimestep
=
update
->
ntimestep
;
next_output
=
ntimestep
+
output_every
;
if
(
output_every
==
0
)
next_output
=
update
->
laststep
+
1
;
if
(
output_every
&&
ntimestep
%
output_every
!=
0
)
next_output
=
(
ntimestep
/
output_every
)
*
output_every
+
output_every
;
// half timestep constraint on pre-step, full timestep thereafter
if
(
strcmp
(
update
->
integrate_style
,
"verlet"
)
==
0
)
{
dtv
=
update
->
dt
;
dtfsq
=
0.5
*
update
->
dt
*
update
->
dt
*
force
->
ftm2v
;
post_force
(
1
);
dtfsq
=
update
->
dt
*
update
->
dt
*
force
->
ftm2v
;
}
else
{
dtv
=
step_respa
[
0
];
dtf_innerhalf
=
0.5
*
step_respa
[
0
]
*
force
->
ftm2v
;
dtf_inner
=
dtf_innerhalf
;
((
Respa
*
)
update
->
integrate
)
->
copy_flevel_f
(
nlevels_respa
-
1
);
post_force_respa
(
1
,
nlevels_respa
-
1
,
0
);
((
Respa
*
)
update
->
integrate
)
->
copy_f_flevel
(
nlevels_respa
-
1
);
dtf_inner
=
step_respa
[
0
]
*
force
->
ftm2v
;
}
}
/* ----------------------------------------------------------------------
build list of SHAKE clusters to constrain
if one or more atoms in cluster are on this proc,
this proc lists the cluster exactly once
------------------------------------------------------------------------- */
void
FixShake
::
pre_neighbor
()
{
int
atom1
,
atom2
,
atom3
,
atom4
;
// local copies of atom quantities
// used by SHAKE until next re-neighboring
x
=
atom
->
x
;
v
=
atom
->
v
;
f
=
atom
->
f
;
mass
=
atom
->
mass
;
type
=
atom
->
type
;
nlocal
=
atom
->
nlocal
;
// extend size of SHAKE list if necessary
if
(
nlocal
>
maxlist
)
{
maxlist
=
nlocal
;
memory
->
sfree
(
list
);
list
=
(
int
*
)
memory
->
smalloc
(
maxlist
*
sizeof
(
int
),
"shake:list"
);
}
// build list of SHAKE clusters I compute
nlist
=
0
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
shake_flag
[
i
])
{
if
(
shake_flag
[
i
]
==
2
)
{
atom1
=
atom
->
map
(
shake_atom
[
i
][
0
]);
atom2
=
atom
->
map
(
shake_atom
[
i
][
1
]);
if
(
atom1
==
-
1
||
atom2
==
-
1
)
{
char
str
[
128
];
sprintf
(
str
,
"Shake atoms %d %d missing on proc %d at step %d"
,
shake_atom
[
i
][
0
],
shake_atom
[
i
][
1
],
me
,
update
->
ntimestep
);
error
->
one
(
str
);
}
if
(
i
<=
atom1
&&
i
<=
atom2
)
list
[
nlist
++
]
=
i
;
}
else
if
(
shake_flag
[
i
]
%
2
==
1
)
{
atom1
=
atom
->
map
(
shake_atom
[
i
][
0
]);
atom2
=
atom
->
map
(
shake_atom
[
i
][
1
]);
atom3
=
atom
->
map
(
shake_atom
[
i
][
2
]);
if
(
atom1
==
-
1
||
atom2
==
-
1
||
atom3
==
-
1
)
{
char
str
[
128
];
sprintf
(
str
,
"Shake atoms %d %d %d missing on proc %d at step %d"
,
shake_atom
[
i
][
0
],
shake_atom
[
i
][
1
],
shake_atom
[
i
][
2
],
me
,
update
->
ntimestep
);
error
->
one
(
str
);
}
if
(
i
<=
atom1
&&
i
<=
atom2
&&
i
<=
atom3
)
list
[
nlist
++
]
=
i
;
}
else
{
atom1
=
atom
->
map
(
shake_atom
[
i
][
0
]);
atom2
=
atom
->
map
(
shake_atom
[
i
][
1
]);
atom3
=
atom
->
map
(
shake_atom
[
i
][
2
]);
atom4
=
atom
->
map
(
shake_atom
[
i
][
3
]);
if
(
atom1
==
-
1
||
atom2
==
-
1
||
atom3
==
-
1
||
atom4
==
-
1
)
{
char
str
[
128
];
sprintf
(
str
,
"Shake atoms %d %d %d %d missing on proc %d at step %d"
,
shake_atom
[
i
][
0
],
shake_atom
[
i
][
1
],
shake_atom
[
i
][
2
],
shake_atom
[
i
][
3
],
me
,
update
->
ntimestep
);
error
->
one
(
str
);
}
if
(
i
<=
atom1
&&
i
<=
atom2
&&
i
<=
atom3
&&
i
<=
atom4
)
list
[
nlist
++
]
=
i
;
}
}
}
/* ----------------------------------------------------------------------
compute the force adjustment for SHAKE constraint
------------------------------------------------------------------------- */
void
FixShake
::
post_force
(
int
vflag_in
)
{
if
(
update
->
ntimestep
==
next_output
)
stats
();
// xshake = unconstrained move with current v,f
unconstrained_update
();
// communicate results if necessary
if
(
nprocs
>
1
)
comm
->
comm_fix
(
this
);
// zero out SHAKE contribution to virial
vflag
=
vflag_in
;
if
(
vflag
)
for
(
int
n
=
0
;
n
<
6
;
n
++
)
virial
[
n
]
=
0.0
;
// loop over clusters
int
m
;
for
(
int
i
=
0
;
i
<
nlist
;
i
++
)
{
m
=
list
[
i
];
if
(
shake_flag
[
m
]
==
2
)
shake2
(
m
);
else
if
(
shake_flag
[
m
]
==
3
)
shake3
(
m
);
else
if
(
shake_flag
[
m
]
==
4
)
shake4
(
m
);
else
shake3angle
(
m
);
}
}
/* ----------------------------------------------------------------------
count # of degrees-of-freedom removed by SHAKE for atoms in igroup
------------------------------------------------------------------------- */
int
FixShake
::
dof
(
int
igroup
)
{
int
groupbit
=
group
->
bitmask
[
igroup
];
int
*
mask
=
atom
->
mask
;
int
*
tag
=
atom
->
tag
;
int
nlocal
=
atom
->
nlocal
;
// count dof in a cluster if and only if the central atom is in group
// and this atom i is the central atom
int
n
=
0
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
!
(
mask
[
i
]
&
groupbit
))
continue
;
if
(
shake_flag
[
i
]
==
0
)
continue
;
if
(
shake_atom
[
i
][
0
]
!=
tag
[
i
])
continue
;
if
(
shake_flag
[
i
]
==
1
)
n
+=
3
;
else
if
(
shake_flag
[
i
]
==
2
)
n
+=
1
;
else
if
(
shake_flag
[
i
]
==
3
)
n
+=
2
;
else
if
(
shake_flag
[
i
]
==
4
)
n
+=
3
;
}
int
nall
;
MPI_Allreduce
(
&
n
,
&
nall
,
1
,
MPI_INT
,
MPI_SUM
,
world
);
return
nall
;
}
/* ----------------------------------------------------------------------
identify whether each atom is in a SHAKE cluster
only include atoms in fix group and those bonds/angles specified in input
test whether all clusters are valid
set shake_flag, shake_atom, shake_type values
set bond,angle types negative so will be ignored in neighbor lists
------------------------------------------------------------------------- */
void
FixShake
::
find_clusters
()
{
int
i
,
j
,
m
,
n
;
int
flag
,
flag_all
,
messtag
,
loop
,
nbuf
,
nbufmax
,
size
;
double
imass
,
jmass
,
rmass
;
int
*
buf
,
*
bufcopy
;
MPI_Request
request
;
MPI_Status
status
;
if
(
me
==
0
&&
screen
)
fprintf
(
screen
,
"Finding SHAKE clusters ...
\n
"
);
// local copies of atom ptrs
int
*
tag
=
atom
->
tag
;
int
*
type
=
atom
->
type
;
int
*
mask
=
atom
->
mask
;
double
*
mass
=
atom
->
mass
;
int
**
bond_type
=
atom
->
bond_type
;
int
**
angle_type
=
atom
->
angle_type
;
int
**
nspecial
=
atom
->
nspecial
;
int
**
special
=
atom
->
special
;
int
nlocal
=
atom
->
nlocal
;
// setup ring of procs
int
next
=
me
+
1
;
int
prev
=
me
-
1
;
if
(
next
==
nprocs
)
next
=
0
;
if
(
prev
<
0
)
prev
=
nprocs
-
1
;
// -----------------------------------------------------
// allocate arrays for self (1d) and bond partners (2d)
// max = max # of bond partners for owned atoms = 2nd dim of partner arrays
// npartner[i] = # of bonds attached to atom i
// nshake[i] = # of SHAKE bonds attached to atom i
// partner_tag[i][] = global IDs of each partner
// partner_mask[i][] = mask of each partner
// partner_type[i][] = type of each partner
// partner_bondtype[i][] = type of bond attached to each partner
// partner_shake[i][] = 1 if SHAKE bonded to partner, 0 if not
// partner_nshake[i][] = nshake value for each partner
// -----------------------------------------------------
int
max
=
0
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
max
=
MAX
(
max
,
nspecial
[
i
][
0
]);
int
*
npartner
=
(
int
*
)
memory
->
smalloc
(
nlocal
*
sizeof
(
double
),
"shake:npartner"
);
int
*
nshake
=
(
int
*
)
memory
->
smalloc
(
nlocal
*
sizeof
(
double
),
"shake:nshake"
);
int
**
partner_tag
=
memory
->
create_2d_int_array
(
nlocal
,
max
,
"shake:partner_tag"
);
int
**
partner_mask
=
memory
->
create_2d_int_array
(
nlocal
,
max
,
"shake:partner_mask"
);
int
**
partner_type
=
memory
->
create_2d_int_array
(
nlocal
,
max
,
"shake:partner_type"
);
int
**
partner_bondtype
=
memory
->
create_2d_int_array
(
nlocal
,
max
,
"shake:partner_bondtype"
);
int
**
partner_shake
=
memory
->
create_2d_int_array
(
nlocal
,
max
,
"shake:partner_shake"
);
int
**
partner_nshake
=
memory
->
create_2d_int_array
(
nlocal
,
max
,
"shake:partner_nshake"
);
// -----------------------------------------------------
// set npartner and partner_tag from special arrays
// -----------------------------------------------------
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
npartner
[
i
]
=
nspecial
[
i
][
0
];
for
(
j
=
0
;
j
<
npartner
[
i
];
j
++
)
partner_tag
[
i
][
j
]
=
special
[
i
][
j
];
}
// -----------------------------------------------------
// set partner_mask, partner_type, partner_bondtype for bonded partners
// requires communication for off-proc partners
// -----------------------------------------------------
// fill in mask, type, bondtype if own bond partner
// info to store in buf for each off-proc bond =
// 2 atoms IDs in bond, space for mask, type, bondtype
// nbufmax = largest buffer needed to hold info from any proc
nbuf
=
0
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
for
(
j
=
0
;
j
<
npartner
[
i
];
j
++
)
{
partner_mask
[
i
][
j
]
=
0
;
partner_type
[
i
][
j
]
=
0
;
partner_bondtype
[
i
][
j
]
=
0
;
m
=
atom
->
map
(
partner_tag
[
i
][
j
]);
if
(
m
>=
0
&&
m
<
nlocal
)
{
partner_mask
[
i
][
j
]
=
mask
[
m
];
partner_type
[
i
][
j
]
=
type
[
m
];
n
=
bondfind
(
i
,
tag
[
i
],
partner_tag
[
i
][
j
]);
if
(
n
>=
0
)
partner_bondtype
[
i
][
j
]
=
bond_type
[
i
][
n
];
else
{
n
=
bondfind
(
m
,
tag
[
i
],
partner_tag
[
i
][
j
]);
if
(
n
>=
0
)
partner_bondtype
[
i
][
j
]
=
bond_type
[
m
][
n
];
}
}
else
nbuf
+=
5
;
}
}
MPI_Allreduce
(
&
nbuf
,
&
nbufmax
,
1
,
MPI_INT
,
MPI_MAX
,
world
);
buf
=
new
int
[
nbufmax
];
bufcopy
=
new
int
[
nbufmax
];
// fill buffer with info
size
=
0
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
for
(
j
=
0
;
j
<
npartner
[
i
];
j
++
)
{
m
=
atom
->
map
(
partner_tag
[
i
][
j
]);
if
(
m
<
0
||
m
>=
nlocal
)
{
buf
[
size
]
=
tag
[
i
];
buf
[
size
+
1
]
=
partner_tag
[
i
][
j
];
buf
[
size
+
2
]
=
0
;
buf
[
size
+
3
]
=
0
;
n
=
bondfind
(
i
,
tag
[
i
],
partner_tag
[
i
][
j
]);
if
(
n
>=
0
)
buf
[
size
+
4
]
=
bond_type
[
i
][
n
];
else
buf
[
size
+
4
]
=
0
;
size
+=
5
;
}
}
}
// cycle buffer around ring of procs back to self
// when receive buffer, scan bond partner IDs for atoms I own
// if I own partner, fill in mask and type, search for bond with 1st atom
// and fill in bondtype
messtag
=
1
;
for
(
loop
=
0
;
loop
<
nprocs
;
loop
++
)
{
i
=
0
;
while
(
i
<
size
)
{
m
=
atom
->
map
(
buf
[
i
+
1
]);
if
(
m
>=
0
&&
m
<
nlocal
)
{
buf
[
i
+
2
]
=
mask
[
m
];
buf
[
i
+
3
]
=
type
[
m
];
if
(
buf
[
i
+
4
]
==
0
)
{
n
=
bondfind
(
m
,
buf
[
i
],
buf
[
i
+
1
]);
if
(
n
>=
0
)
buf
[
i
+
4
]
=
bond_type
[
m
][
n
];
}
}
i
+=
5
;
}
if
(
me
!=
next
)
{
MPI_Irecv
(
bufcopy
,
nbufmax
,
MPI_INT
,
prev
,
messtag
,
world
,
&
request
);
MPI_Send
(
buf
,
size
,
MPI_INT
,
next
,
messtag
,
world
);
MPI_Wait
(
&
request
,
&
status
);
MPI_Get_count
(
&
status
,
MPI_INT
,
&
size
);
for
(
j
=
0
;
j
<
size
;
j
++
)
buf
[
j
]
=
bufcopy
[
j
];
}
}
// store partner info returned to me
m
=
0
;
while
(
m
<
size
)
{
i
=
atom
->
map
(
buf
[
m
]);
for
(
j
=
0
;
j
<
npartner
[
i
];
j
++
)
if
(
buf
[
m
+
1
]
==
partner_tag
[
i
][
j
])
break
;
partner_mask
[
i
][
j
]
=
buf
[
m
+
2
];
partner_type
[
i
][
j
]
=
buf
[
m
+
3
];
partner_bondtype
[
i
][
j
]
=
buf
[
m
+
4
];
m
+=
5
;
}
delete
[]
buf
;
delete
[]
bufcopy
;
// error check for unfilled partner info
// if partner_type not set, is an error
// partner_bondtype may not be set if special list is not consistent
// with bondatom (e.g. due to delete_bonds command)
// this is OK if one or both atoms are not in fix group, since
// bond won't be SHAKEn anyway
// else it's an error
flag
=
0
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
for
(
j
=
0
;
j
<
npartner
[
i
];
j
++
)
{
if
(
partner_type
[
i
][
j
]
==
0
)
flag
=
1
;
if
(
!
(
mask
[
i
]
&
groupbit
))
continue
;
if
(
!
(
partner_mask
[
i
][
j
]
&
groupbit
))
continue
;
if
(
partner_bondtype
[
i
][
j
]
==
0
)
flag
=
1
;
}
MPI_Allreduce
(
&
flag
,
&
flag_all
,
1
,
MPI_INT
,
MPI_SUM
,
world
);
if
(
flag_all
)
error
->
all
(
"Did not find fix shake partner info"
);
// -----------------------------------------------------
// identify SHAKEable bonds
// set nshake[i] = # of SHAKE bonds attached to atom i
// set partner_shake[i][] = 1 if SHAKE bonded to partner, 0 if not
// both atoms must be in group, bondtype must be > 0
// check if bondtype is in input bond_flag
// check if type of either atom is in input type_flag
// check if mass of either atom is in input mass_list
// -----------------------------------------------------
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
nshake
[
i
]
=
0
;
for
(
j
=
0
;
j
<
npartner
[
i
];
j
++
)
{
partner_shake
[
i
][
j
]
=
0
;
if
(
!
(
mask
[
i
]
&
groupbit
))
continue
;
if
(
!
(
partner_mask
[
i
][
j
]
&
groupbit
))
continue
;
if
(
partner_bondtype
[
i
][
j
]
<=
0
)
continue
;
if
(
bond_flag
[
partner_bondtype
[
i
][
j
]])
{
partner_shake
[
i
][
j
]
=
1
;
nshake
[
i
]
++
;
continue
;
}
if
(
type_flag
[
type
[
i
]]
||
type_flag
[
partner_type
[
i
][
j
]])
{
partner_shake
[
i
][
j
]
=
1
;
nshake
[
i
]
++
;
continue
;
}
if
(
nmass
)
{
imass
=
mass
[
type
[
i
]];
jmass
=
mass
[
partner_type
[
i
][
j
]];
for
(
m
=
0
;
m
<
nmass
;
m
++
)
{
rmass
=
mass_list
[
m
];
if
(
fabs
(
rmass
-
imass
)
<=
MASSDELTA
||
fabs
(
rmass
-
jmass
)
<=
MASSDELTA
)
{
partner_shake
[
i
][
j
]
=
1
;
nshake
[
i
]
++
;
break
;
}
}
}
}
}
// -----------------------------------------------------
// set partner_nshake for bonded partners
// requires communication for off-proc partners
// -----------------------------------------------------
// fill in partner_nshake if own bond partner
// info to store in buf for each off-proc bond =
// 2 atoms IDs in bond, space for nshake value
// nbufmax = largest buffer needed to hold info from any proc
nbuf
=
0
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
for
(
j
=
0
;
j
<
npartner
[
i
];
j
++
)
{
m
=
atom
->
map
(
partner_tag
[
i
][
j
]);
if
(
m
>=
0
&&
m
<
nlocal
)
partner_nshake
[
i
][
j
]
=
nshake
[
m
];
else
nbuf
+=
3
;
}
}
MPI_Allreduce
(
&
nbuf
,
&
nbufmax
,
1
,
MPI_INT
,
MPI_MAX
,
world
);
buf
=
new
int
[
nbufmax
];
bufcopy
=
new
int
[
nbufmax
];
// fill buffer with info
size
=
0
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
for
(
j
=
0
;
j
<
npartner
[
i
];
j
++
)
{
m
=
atom
->
map
(
partner_tag
[
i
][
j
]);
if
(
m
<
0
||
m
>=
nlocal
)
{
buf
[
size
]
=
tag
[
i
];
buf
[
size
+
1
]
=
partner_tag
[
i
][
j
];
size
+=
3
;
}
}
}
// cycle buffer around ring of procs back to self
// when receive buffer, scan bond partner IDs for atoms I own
// if I own partner, fill in nshake value
messtag
=
2
;
for
(
loop
=
0
;
loop
<
nprocs
;
loop
++
)
{
i
=
0
;
while
(
i
<
size
)
{
m
=
atom
->
map
(
buf
[
i
+
1
]);
if
(
m
>=
0
&&
m
<
nlocal
)
buf
[
i
+
2
]
=
nshake
[
m
];
i
+=
3
;
}
if
(
me
!=
next
)
{
MPI_Irecv
(
bufcopy
,
nbufmax
,
MPI_INT
,
prev
,
messtag
,
world
,
&
request
);
MPI_Send
(
buf
,
size
,
MPI_INT
,
next
,
messtag
,
world
);
MPI_Wait
(
&
request
,
&
status
);
MPI_Get_count
(
&
status
,
MPI_INT
,
&
size
);
for
(
j
=
0
;
j
<
size
;
j
++
)
buf
[
j
]
=
bufcopy
[
j
];
}
}
// store partner info returned to me
m
=
0
;
while
(
m
<
size
)
{
i
=
atom
->
map
(
buf
[
m
]);
for
(
j
=
0
;
j
<
npartner
[
i
];
j
++
)
if
(
buf
[
m
+
1
]
==
partner_tag
[
i
][
j
])
break
;
partner_nshake
[
i
][
j
]
=
buf
[
m
+
2
];
m
+=
3
;
}
delete
[]
buf
;
delete
[]
bufcopy
;
// -----------------------------------------------------
// error checks
// no atom with nshake > 3
// no connected atoms which both have nshake > 1
// -----------------------------------------------------
flag
=
0
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
nshake
[
i
]
>
3
)
flag
=
1
;
MPI_Allreduce
(
&
flag
,
&
flag_all
,
1
,
MPI_INT
,
MPI_SUM
,
world
);
if
(
flag_all
)
error
->
all
(
"Shake cluster of more than 4 atoms"
);
flag
=
0
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
nshake
[
i
]
<=
1
)
continue
;
for
(
j
=
0
;
j
<
npartner
[
i
];
j
++
)
if
(
partner_shake
[
i
][
j
]
&&
partner_nshake
[
i
][
j
]
>
1
)
flag
=
1
;
}
MPI_Allreduce
(
&
flag
,
&
flag_all
,
1
,
MPI_INT
,
MPI_SUM
,
world
);
if
(
flag_all
)
error
->
all
(
"Shake clusters are connected"
);
// -----------------------------------------------------
// set SHAKE arrays that are stored with atoms & add angle constraints
// zero shake arrays for all owned atoms
// if I am central atom set shake_flag & shake_atom & shake_type
// for 2-atom clusters, I am central atom if my atom ID < partner ID
// for 3-atom clusters, test for angle constraint
// angle will be stored by this atom if it exists
// if angle type matches angle_flag, then it is angle-constrained
// shake_flag[] = 0 if atom not in SHAKE cluster
// 2,3,4 = size of bond-only cluster
// 1 = 3-atom angle cluster
// shake_atom[][] = global IDs of 2,3,4 atoms in cluster
// central atom is 1st
// for 2-atom cluster, lowest ID is 1st
// shake_type[][] = bondtype of each bond in cluster
// for 3-atom angle cluster, 3rd value is angletype
// -----------------------------------------------------
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
shake_flag
[
i
]
=
0
;
shake_atom
[
i
][
0
]
=
0
;
shake_atom
[
i
][
1
]
=
0
;
shake_atom
[
i
][
2
]
=
0
;
shake_atom
[
i
][
3
]
=
0
;
shake_type
[
i
][
0
]
=
0
;
shake_type
[
i
][
1
]
=
0
;
shake_type
[
i
][
2
]
=
0
;
if
(
nshake
[
i
]
==
1
)
{
for
(
j
=
0
;
j
<
npartner
[
i
];
j
++
)
if
(
partner_shake
[
i
][
j
])
break
;
if
(
partner_nshake
[
i
][
j
]
==
1
&&
tag
[
i
]
<
partner_tag
[
i
][
j
])
{
shake_flag
[
i
]
=
2
;
shake_atom
[
i
][
0
]
=
tag
[
i
];
shake_atom
[
i
][
1
]
=
partner_tag
[
i
][
j
];
shake_type
[
i
][
0
]
=
partner_bondtype
[
i
][
j
];
}
}
if
(
nshake
[
i
]
>
1
)
{
shake_flag
[
i
]
=
1
;
shake_atom
[
i
][
0
]
=
tag
[
i
];
for
(
j
=
0
;
j
<
npartner
[
i
];
j
++
)
if
(
partner_shake
[
i
][
j
])
{
m
=
shake_flag
[
i
];
shake_atom
[
i
][
m
]
=
partner_tag
[
i
][
j
];
shake_type
[
i
][
m
-
1
]
=
partner_bondtype
[
i
][
j
];
shake_flag
[
i
]
++
;
}
}
if
(
nshake
[
i
]
==
2
)
{
n
=
anglefind
(
i
,
shake_atom
[
i
][
1
],
shake_atom
[
i
][
2
]);
if
(
n
<
0
)
continue
;
if
(
angle_type
[
i
][
n
]
<
0
)
continue
;
if
(
angle_flag
[
angle_type
[
i
][
n
]])
{
shake_flag
[
i
]
=
1
;
shake_type
[
i
][
2
]
=
angle_type
[
i
][
n
];
}
}
}
// -----------------------------------------------------
// set shake_flag,shake_atom,shake_type for non-central atoms
// requires communication for off-proc atoms
// -----------------------------------------------------
// fill in shake arrays for each bond partner I own
// info to store in buf for each off-proc bond =
// all values from shake_flag, shake_atom, shake_type
// nbufmax = largest buffer needed to hold info from any proc
nbuf
=
0
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
shake_flag
[
i
]
==
0
)
continue
;
for
(
j
=
0
;
j
<
npartner
[
i
];
j
++
)
{
if
(
partner_shake
[
i
][
j
]
==
0
)
continue
;
m
=
atom
->
map
(
partner_tag
[
i
][
j
]);
if
(
m
>=
0
&&
m
<
nlocal
)
{
shake_flag
[
m
]
=
shake_flag
[
i
];
shake_atom
[
m
][
0
]
=
shake_atom
[
i
][
0
];
shake_atom
[
m
][
1
]
=
shake_atom
[
i
][
1
];
shake_atom
[
m
][
2
]
=
shake_atom
[
i
][
2
];
shake_atom
[
m
][
3
]
=
shake_atom
[
i
][
3
];
shake_type
[
m
][
0
]
=
shake_type
[
i
][
0
];
shake_type
[
m
][
1
]
=
shake_type
[
i
][
1
];
shake_type
[
m
][
2
]
=
shake_type
[
i
][
2
];
}
else
nbuf
+=
9
;
}
}
MPI_Allreduce
(
&
nbuf
,
&
nbufmax
,
1
,
MPI_INT
,
MPI_MAX
,
world
);
buf
=
new
int
[
nbufmax
];
bufcopy
=
new
int
[
nbufmax
];
// fill buffer with info
size
=
0
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
shake_flag
[
i
]
==
0
)
continue
;
for
(
j
=
0
;
j
<
npartner
[
i
];
j
++
)
{
if
(
partner_shake
[
i
][
j
]
==
0
)
continue
;
m
=
atom
->
map
(
partner_tag
[
i
][
j
]);
if
(
m
<
0
||
m
>=
nlocal
)
{
buf
[
size
]
=
partner_tag
[
i
][
j
];
buf
[
size
+
1
]
=
shake_flag
[
i
];
buf
[
size
+
2
]
=
shake_atom
[
i
][
0
];
buf
[
size
+
3
]
=
shake_atom
[
i
][
1
];
buf
[
size
+
4
]
=
shake_atom
[
i
][
2
];
buf
[
size
+
5
]
=
shake_atom
[
i
][
3
];
buf
[
size
+
6
]
=
shake_type
[
i
][
0
];
buf
[
size
+
7
]
=
shake_type
[
i
][
1
];
buf
[
size
+
8
]
=
shake_type
[
i
][
2
];
size
+=
9
;
}
}
}
// cycle buffer around ring of procs back to self
// when receive buffer, scan for ID that I own
// if I own ID, fill in shake array values
messtag
=
3
;
for
(
loop
=
0
;
loop
<
nprocs
;
loop
++
)
{
i
=
0
;
while
(
i
<
size
)
{
m
=
atom
->
map
(
buf
[
i
]);
if
(
m
>=
0
&&
m
<
nlocal
)
{
shake_flag
[
m
]
=
buf
[
i
+
1
];
shake_atom
[
m
][
0
]
=
buf
[
i
+
2
];
shake_atom
[
m
][
1
]
=
buf
[
i
+
3
];
shake_atom
[
m
][
2
]
=
buf
[
i
+
4
];
shake_atom
[
m
][
3
]
=
buf
[
i
+
5
];
shake_type
[
m
][
0
]
=
buf
[
i
+
6
];
shake_type
[
m
][
1
]
=
buf
[
i
+
7
];
shake_type
[
m
][
2
]
=
buf
[
i
+
8
];
}
i
+=
9
;
}
if
(
me
!=
next
)
{
MPI_Irecv
(
bufcopy
,
nbufmax
,
MPI_INT
,
prev
,
messtag
,
world
,
&
request
);
MPI_Send
(
buf
,
size
,
MPI_INT
,
next
,
messtag
,
world
);
MPI_Wait
(
&
request
,
&
status
);
MPI_Get_count
(
&
status
,
MPI_INT
,
&
size
);
for
(
j
=
0
;
j
<
size
;
j
++
)
buf
[
j
]
=
bufcopy
[
j
];
}
}
delete
[]
buf
;
delete
[]
bufcopy
;
// -----------------------------------------------------
// free local memory
// -----------------------------------------------------
memory
->
sfree
(
npartner
);
memory
->
sfree
(
nshake
);
memory
->
destroy_2d_int_array
(
partner_tag
);
memory
->
destroy_2d_int_array
(
partner_mask
);
memory
->
destroy_2d_int_array
(
partner_type
);
memory
->
destroy_2d_int_array
(
partner_bondtype
);
memory
->
destroy_2d_int_array
(
partner_shake
);
memory
->
destroy_2d_int_array
(
partner_nshake
);
// -----------------------------------------------------
// set bond_type and angle_type negative for SHAKE clusters
// must set for all SHAKE bonds and angles stored by each atom
// -----------------------------------------------------
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
shake_flag
[
i
]
==
0
)
continue
;
else
if
(
shake_flag
[
i
]
==
1
)
{
n
=
bondfind
(
i
,
shake_atom
[
i
][
0
],
shake_atom
[
i
][
1
]);
if
(
n
>=
0
)
bond_type
[
i
][
n
]
=
-
bond_type
[
i
][
n
];
n
=
bondfind
(
i
,
shake_atom
[
i
][
0
],
shake_atom
[
i
][
2
]);
if
(
n
>=
0
)
bond_type
[
i
][
n
]
=
-
bond_type
[
i
][
n
];
n
=
anglefind
(
i
,
shake_atom
[
i
][
1
],
shake_atom
[
i
][
2
]);
if
(
n
>=
0
)
angle_type
[
i
][
n
]
=
-
angle_type
[
i
][
n
];
}
else
if
(
shake_flag
[
i
]
==
2
)
{
n
=
bondfind
(
i
,
shake_atom
[
i
][
0
],
shake_atom
[
i
][
1
]);
if
(
n
>=
0
)
bond_type
[
i
][
n
]
=
-
bond_type
[
i
][
n
];
}
else
if
(
shake_flag
[
i
]
==
3
)
{
n
=
bondfind
(
i
,
shake_atom
[
i
][
0
],
shake_atom
[
i
][
1
]);
if
(
n
>=
0
)
bond_type
[
i
][
n
]
=
-
bond_type
[
i
][
n
];
n
=
bondfind
(
i
,
shake_atom
[
i
][
0
],
shake_atom
[
i
][
2
]);
if
(
n
>=
0
)
bond_type
[
i
][
n
]
=
-
bond_type
[
i
][
n
];
}
else
if
(
shake_flag
[
i
]
==
4
)
{
n
=
bondfind
(
i
,
shake_atom
[
i
][
0
],
shake_atom
[
i
][
1
]);
if
(
n
>=
0
)
bond_type
[
i
][
n
]
=
-
bond_type
[
i
][
n
];
n
=
bondfind
(
i
,
shake_atom
[
i
][
0
],
shake_atom
[
i
][
2
]);
if
(
n
>=
0
)
bond_type
[
i
][
n
]
=
-
bond_type
[
i
][
n
];
n
=
bondfind
(
i
,
shake_atom
[
i
][
0
],
shake_atom
[
i
][
3
]);
if
(
n
>=
0
)
bond_type
[
i
][
n
]
=
-
bond_type
[
i
][
n
];
}
}
// -----------------------------------------------------
// print info on SHAKE clusters
// -----------------------------------------------------
int
count1
,
count2
,
count3
,
count4
;
count1
=
count2
=
count3
=
count4
=
0
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
shake_flag
[
i
]
==
1
)
count1
++
;
else
if
(
shake_flag
[
i
]
==
2
)
count2
++
;
else
if
(
shake_flag
[
i
]
==
3
)
count3
++
;
else
if
(
shake_flag
[
i
]
==
4
)
count4
++
;
}
int
tmp
;
tmp
=
count1
;
MPI_Allreduce
(
&
tmp
,
&
count1
,
1
,
MPI_INT
,
MPI_SUM
,
world
);
tmp
=
count2
;
MPI_Allreduce
(
&
tmp
,
&
count2
,
1
,
MPI_INT
,
MPI_SUM
,
world
);
tmp
=
count3
;
MPI_Allreduce
(
&
tmp
,
&
count3
,
1
,
MPI_INT
,
MPI_SUM
,
world
);
tmp
=
count4
;
MPI_Allreduce
(
&
tmp
,
&
count4
,
1
,
MPI_INT
,
MPI_SUM
,
world
);
if
(
me
==
0
)
{
if
(
screen
)
{
fprintf
(
screen
,
" %d = # of size 2 clusters
\n
"
,
count2
/
2
);
fprintf
(
screen
,
" %d = # of size 3 clusters
\n
"
,
count3
/
3
);
fprintf
(
screen
,
" %d = # of size 4 clusters
\n
"
,
count4
/
4
);
fprintf
(
screen
,
" %d = # of frozen angles
\n
"
,
count1
/
3
);
}
if
(
logfile
)
{
fprintf
(
logfile
,
" %d = # of size 2 clusters
\n
"
,
count2
/
2
);
fprintf
(
logfile
,
" %d = # of size 3 clusters
\n
"
,
count3
/
3
);
fprintf
(
logfile
,
" %d = # of size 4 clusters
\n
"
,
count4
/
4
);
fprintf
(
logfile
,
" %d = # of frozen angles
\n
"
,
count1
/
3
);
}
}
}
/* ----------------------------------------------------------------------
update the unconstrained position of each atom
only for SHAKE clusters, else set to 0.0
assumes NVE update, seems to be accurate enough for NVT,NPT,NPH as well
------------------------------------------------------------------------- */
void
FixShake
::
unconstrained_update
()
{
double
dtfmsq
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
shake_flag
[
i
])
{
dtfmsq
=
dtfsq
/
mass
[
type
[
i
]];
xshake
[
i
][
0
]
=
x
[
i
][
0
]
+
dtv
*
v
[
i
][
0
]
+
dtfmsq
*
f
[
i
][
0
];
xshake
[
i
][
1
]
=
x
[
i
][
1
]
+
dtv
*
v
[
i
][
1
]
+
dtfmsq
*
f
[
i
][
1
];
xshake
[
i
][
2
]
=
x
[
i
][
2
]
+
dtv
*
v
[
i
][
2
]
+
dtfmsq
*
f
[
i
][
2
];
}
else
xshake
[
i
][
2
]
=
xshake
[
i
][
1
]
=
xshake
[
i
][
0
]
=
0.0
;
}
}
/* ---------------------------------------------------------------------- */
void
FixShake
::
shake2
(
int
m
)
{
// local atom IDs and constraint distances
int
i0
=
atom
->
map
(
shake_atom
[
m
][
0
]);
int
i1
=
atom
->
map
(
shake_atom
[
m
][
1
]);
double
bond1
=
bond_distance
[
shake_type
[
m
][
0
]];
// r01 = distance vec between atoms, with PBC
double
r01
[
3
];
r01
[
0
]
=
x
[
i0
][
0
]
-
x
[
i1
][
0
];
r01
[
1
]
=
x
[
i0
][
1
]
-
x
[
i1
][
1
];
r01
[
2
]
=
x
[
i0
][
2
]
-
x
[
i1
][
2
];
domain
->
minimum_image
(
r01
);
// s01 = distance vec after unconstrained update, with PBC
double
s01
[
3
];
s01
[
0
]
=
xshake
[
i0
][
0
]
-
xshake
[
i1
][
0
];
s01
[
1
]
=
xshake
[
i0
][
1
]
-
xshake
[
i1
][
1
];
s01
[
2
]
=
xshake
[
i0
][
2
]
-
xshake
[
i1
][
2
];
domain
->
minimum_image
(
s01
);
// scalar distances between atoms
double
r01sq
=
r01
[
0
]
*
r01
[
0
]
+
r01
[
1
]
*
r01
[
1
]
+
r01
[
2
]
*
r01
[
2
];
double
s01sq
=
s01
[
0
]
*
s01
[
0
]
+
s01
[
1
]
*
s01
[
1
]
+
s01
[
2
]
*
s01
[
2
];
// a,b,c = coeffs in quadratic equation for lamda
double
invmass0
=
1.0
/
mass
[
type
[
i0
]];
double
invmass1
=
1.0
/
mass
[
type
[
i1
]];
double
a
=
(
invmass0
+
invmass1
)
*
(
invmass0
+
invmass1
)
*
r01sq
;
double
b
=
2.0
*
(
invmass0
+
invmass1
)
*
(
s01
[
0
]
*
r01
[
0
]
+
s01
[
1
]
*
r01
[
1
]
+
s01
[
2
]
*
r01
[
2
]);
double
c
=
s01sq
-
bond1
*
bond1
;
// error check
double
determ
=
b
*
b
-
4.0
*
a
*
c
;
if
(
determ
<
0.0
)
{
error
->
warning
(
"Shake determinant < 0.0"
);
determ
=
0.0
;
}
// exact quadratic solution for lamda
double
lamda
,
lamda1
,
lamda2
;
lamda1
=
(
-
b
+
sqrt
(
determ
))
/
(
2.0
*
a
);
lamda2
=
(
-
b
-
sqrt
(
determ
))
/
(
2.0
*
a
);
if
(
fabs
(
lamda1
)
<=
fabs
(
lamda2
))
lamda
=
lamda1
;
else
lamda
=
lamda2
;
// update forces if atom is owned by this processor
lamda
/=
dtfsq
;
if
(
i0
<
nlocal
)
{
f
[
i0
][
0
]
+=
lamda
*
r01
[
0
];
f
[
i0
][
1
]
+=
lamda
*
r01
[
1
];
f
[
i0
][
2
]
+=
lamda
*
r01
[
2
];
}
if
(
i1
<
nlocal
)
{
f
[
i1
][
0
]
-=
lamda
*
r01
[
0
];
f
[
i1
][
1
]
-=
lamda
*
r01
[
1
];
f
[
i1
][
2
]
-=
lamda
*
r01
[
2
];
}
if
(
vflag
)
{
int
factor
=
0
;
if
(
i0
<
nlocal
)
factor
++
;
if
(
i1
<
nlocal
)
factor
++
;
double
rfactor
=
0.5
*
factor
;
virial
[
0
]
+=
rfactor
*
lamda
*
r01
[
0
]
*
r01
[
0
];
virial
[
1
]
+=
rfactor
*
lamda
*
r01
[
1
]
*
r01
[
1
];
virial
[
2
]
+=
rfactor
*
lamda
*
r01
[
2
]
*
r01
[
2
];
virial
[
3
]
+=
rfactor
*
lamda
*
r01
[
0
]
*
r01
[
1
];
virial
[
4
]
+=
rfactor
*
lamda
*
r01
[
0
]
*
r01
[
2
];
virial
[
5
]
+=
rfactor
*
lamda
*
r01
[
1
]
*
r01
[
2
];
}
}
/* ---------------------------------------------------------------------- */
void
FixShake
::
shake3
(
int
m
)
{
// local atom IDs and constraint distances
int
i0
=
atom
->
map
(
shake_atom
[
m
][
0
]);
int
i1
=
atom
->
map
(
shake_atom
[
m
][
1
]);
int
i2
=
atom
->
map
(
shake_atom
[
m
][
2
]);
double
bond1
=
bond_distance
[
shake_type
[
m
][
0
]];
double
bond2
=
bond_distance
[
shake_type
[
m
][
1
]];
// r01,r02 = distance vec between atoms, with PBC
double
r01
[
3
];
r01
[
0
]
=
x
[
i0
][
0
]
-
x
[
i1
][
0
];
r01
[
1
]
=
x
[
i0
][
1
]
-
x
[
i1
][
1
];
r01
[
2
]
=
x
[
i0
][
2
]
-
x
[
i1
][
2
];
domain
->
minimum_image
(
r01
);
double
r02
[
3
];
r02
[
0
]
=
x
[
i0
][
0
]
-
x
[
i2
][
0
];
r02
[
1
]
=
x
[
i0
][
1
]
-
x
[
i2
][
1
];
r02
[
2
]
=
x
[
i0
][
2
]
-
x
[
i2
][
2
];
domain
->
minimum_image
(
r02
);
// s01,s02 = distance vec after unconstrained update, with PBC
double
s01
[
3
];
s01
[
0
]
=
xshake
[
i0
][
0
]
-
xshake
[
i1
][
0
];
s01
[
1
]
=
xshake
[
i0
][
1
]
-
xshake
[
i1
][
1
];
s01
[
2
]
=
xshake
[
i0
][
2
]
-
xshake
[
i1
][
2
];
domain
->
minimum_image
(
s01
);
double
s02
[
3
];
s02
[
0
]
=
xshake
[
i0
][
0
]
-
xshake
[
i2
][
0
];
s02
[
1
]
=
xshake
[
i0
][
1
]
-
xshake
[
i2
][
1
];
s02
[
2
]
=
xshake
[
i0
][
2
]
-
xshake
[
i2
][
2
];
domain
->
minimum_image
(
s02
);
// scalar distances between atoms
double
r01sq
=
r01
[
0
]
*
r01
[
0
]
+
r01
[
1
]
*
r01
[
1
]
+
r01
[
2
]
*
r01
[
2
];
double
r02sq
=
r02
[
0
]
*
r02
[
0
]
+
r02
[
1
]
*
r02
[
1
]
+
r02
[
2
]
*
r02
[
2
];
double
s01sq
=
s01
[
0
]
*
s01
[
0
]
+
s01
[
1
]
*
s01
[
1
]
+
s01
[
2
]
*
s01
[
2
];
double
s02sq
=
s02
[
0
]
*
s02
[
0
]
+
s02
[
1
]
*
s02
[
1
]
+
s02
[
2
]
*
s02
[
2
];
// matrix coeffs and rhs for lamda equations
double
invmass0
=
1.0
/
mass
[
type
[
i0
]];
double
invmass1
=
1.0
/
mass
[
type
[
i1
]];
double
invmass2
=
1.0
/
mass
[
type
[
i2
]];
double
a11
=
2.0
*
(
invmass0
+
invmass1
)
*
(
s01
[
0
]
*
r01
[
0
]
+
s01
[
1
]
*
r01
[
1
]
+
s01
[
2
]
*
r01
[
2
]);
double
a12
=
2.0
*
invmass0
*
(
s01
[
0
]
*
r02
[
0
]
+
s01
[
1
]
*
r02
[
1
]
+
s01
[
2
]
*
r02
[
2
]);
double
a21
=
2.0
*
invmass0
*
(
s02
[
0
]
*
r01
[
0
]
+
s02
[
1
]
*
r01
[
1
]
+
s02
[
2
]
*
r01
[
2
]);
double
a22
=
2.0
*
(
invmass0
+
invmass2
)
*
(
s02
[
0
]
*
r02
[
0
]
+
s02
[
1
]
*
r02
[
1
]
+
s02
[
2
]
*
r02
[
2
]);
// inverse of matrix
double
determ
=
a11
*
a22
-
a12
*
a21
;
if
(
determ
==
0.0
)
error
->
one
(
"Shake determinant = 0.0"
);
double
determinv
=
1.0
/
determ
;
double
a11inv
=
a22
*
determinv
;
double
a12inv
=
-
a12
*
determinv
;
double
a21inv
=
-
a21
*
determinv
;
double
a22inv
=
a11
*
determinv
;
// quadratic correction coeffs
double
r0102
=
(
r01
[
0
]
*
r02
[
0
]
+
r01
[
1
]
*
r02
[
1
]
+
r01
[
2
]
*
r02
[
2
]);
double
quad1_0101
=
(
invmass0
+
invmass1
)
*
(
invmass0
+
invmass1
)
*
r01sq
;
double
quad1_0202
=
invmass0
*
invmass0
*
r02sq
;
double
quad1_0102
=
2.0
*
(
invmass0
+
invmass1
)
*
invmass0
*
r0102
;
double
quad2_0202
=
(
invmass0
+
invmass2
)
*
(
invmass0
+
invmass2
)
*
r02sq
;
double
quad2_0101
=
invmass0
*
invmass0
*
r01sq
;
double
quad2_0102
=
2.0
*
(
invmass0
+
invmass2
)
*
invmass0
*
r0102
;
// iterate until converged
double
lamda01
=
0.0
;
double
lamda02
=
0.0
;
int
niter
=
0
;
int
done
=
0
;
double
quad1
,
quad2
,
b1
,
b2
,
lamda01_new
,
lamda02_new
;
while
(
!
done
&&
niter
<
max_iter
)
{
quad1
=
quad1_0101
*
lamda01
*
lamda01
+
quad1_0202
*
lamda02
*
lamda02
+
quad1_0102
*
lamda01
*
lamda02
;
quad2
=
quad2_0101
*
lamda01
*
lamda01
+
quad2_0202
*
lamda02
*
lamda02
+
quad2_0102
*
lamda01
*
lamda02
;
b1
=
bond1
*
bond1
-
s01sq
-
quad1
;
b2
=
bond2
*
bond2
-
s02sq
-
quad2
;
lamda01_new
=
a11inv
*
b1
+
a12inv
*
b2
;
lamda02_new
=
a21inv
*
b1
+
a22inv
*
b2
;
done
=
1
;
if
(
fabs
(
lamda01_new
-
lamda01
)
>
tolerance
)
done
=
0
;
if
(
fabs
(
lamda02_new
-
lamda02
)
>
tolerance
)
done
=
0
;
lamda01
=
lamda01_new
;
lamda02
=
lamda02_new
;
niter
++
;
}
// update forces if atom is owned by this processor
lamda01
=
lamda01
/
dtfsq
;
lamda02
=
lamda02
/
dtfsq
;
if
(
i0
<
nlocal
)
{
f
[
i0
][
0
]
+=
lamda01
*
r01
[
0
]
+
lamda02
*
r02
[
0
];
f
[
i0
][
1
]
+=
lamda01
*
r01
[
1
]
+
lamda02
*
r02
[
1
];
f
[
i0
][
2
]
+=
lamda01
*
r01
[
2
]
+
lamda02
*
r02
[
2
];
}
if
(
i1
<
nlocal
)
{
f
[
i1
][
0
]
-=
lamda01
*
r01
[
0
];
f
[
i1
][
1
]
-=
lamda01
*
r01
[
1
];
f
[
i1
][
2
]
-=
lamda01
*
r01
[
2
];
}
if
(
i2
<
nlocal
)
{
f
[
i2
][
0
]
-=
lamda02
*
r02
[
0
];
f
[
i2
][
1
]
-=
lamda02
*
r02
[
1
];
f
[
i2
][
2
]
-=
lamda02
*
r02
[
2
];
}
if
(
vflag
)
{
int
factor
=
0
;
if
(
i0
<
nlocal
)
factor
++
;
if
(
i1
<
nlocal
)
factor
++
;
if
(
i2
<
nlocal
)
factor
++
;
double
rfactor
=
factor
/
3.0
;
virial
[
0
]
+=
rfactor
*
lamda01
*
r01
[
0
]
*
r01
[
0
];
virial
[
1
]
+=
rfactor
*
lamda01
*
r01
[
1
]
*
r01
[
1
];
virial
[
2
]
+=
rfactor
*
lamda01
*
r01
[
2
]
*
r01
[
2
];
virial
[
3
]
+=
rfactor
*
lamda01
*
r01
[
0
]
*
r01
[
1
];
virial
[
4
]
+=
rfactor
*
lamda01
*
r01
[
0
]
*
r01
[
2
];
virial
[
5
]
+=
rfactor
*
lamda01
*
r01
[
1
]
*
r01
[
2
];
virial
[
0
]
+=
rfactor
*
lamda02
*
r02
[
0
]
*
r02
[
0
];
virial
[
1
]
+=
rfactor
*
lamda02
*
r02
[
1
]
*
r02
[
1
];
virial
[
2
]
+=
rfactor
*
lamda02
*
r02
[
2
]
*
r02
[
2
];
virial
[
3
]
+=
rfactor
*
lamda02
*
r02
[
0
]
*
r02
[
1
];
virial
[
4
]
+=
rfactor
*
lamda02
*
r02
[
0
]
*
r02
[
2
];
virial
[
5
]
+=
rfactor
*
lamda02
*
r02
[
1
]
*
r02
[
2
];
}
}
/* ---------------------------------------------------------------------- */
void
FixShake
::
shake4
(
int
m
)
{
// local atom IDs and constraint distances
int
i0
=
atom
->
map
(
shake_atom
[
m
][
0
]);
int
i1
=
atom
->
map
(
shake_atom
[
m
][
1
]);
int
i2
=
atom
->
map
(
shake_atom
[
m
][
2
]);
int
i3
=
atom
->
map
(
shake_atom
[
m
][
3
]);
double
bond1
=
bond_distance
[
shake_type
[
m
][
0
]];
double
bond2
=
bond_distance
[
shake_type
[
m
][
1
]];
double
bond3
=
bond_distance
[
shake_type
[
m
][
2
]];
// r01,r02,r03 = distance vec between atoms, with PBC
double
r01
[
3
];
r01
[
0
]
=
x
[
i0
][
0
]
-
x
[
i1
][
0
];
r01
[
1
]
=
x
[
i0
][
1
]
-
x
[
i1
][
1
];
r01
[
2
]
=
x
[
i0
][
2
]
-
x
[
i1
][
2
];
domain
->
minimum_image
(
r01
);
double
r02
[
3
];
r02
[
0
]
=
x
[
i0
][
0
]
-
x
[
i2
][
0
];
r02
[
1
]
=
x
[
i0
][
1
]
-
x
[
i2
][
1
];
r02
[
2
]
=
x
[
i0
][
2
]
-
x
[
i2
][
2
];
domain
->
minimum_image
(
r02
);
double
r03
[
3
];
r03
[
0
]
=
x
[
i0
][
0
]
-
x
[
i3
][
0
];
r03
[
1
]
=
x
[
i0
][
1
]
-
x
[
i3
][
1
];
r03
[
2
]
=
x
[
i0
][
2
]
-
x
[
i3
][
2
];
domain
->
minimum_image
(
r03
);
// s01,s02,s03 = distance vec after unconstrained update, with PBC
double
s01
[
3
];
s01
[
0
]
=
xshake
[
i0
][
0
]
-
xshake
[
i1
][
0
];
s01
[
1
]
=
xshake
[
i0
][
1
]
-
xshake
[
i1
][
1
];
s01
[
2
]
=
xshake
[
i0
][
2
]
-
xshake
[
i1
][
2
];
domain
->
minimum_image
(
s01
);
double
s02
[
3
];
s02
[
0
]
=
xshake
[
i0
][
0
]
-
xshake
[
i2
][
0
];
s02
[
1
]
=
xshake
[
i0
][
1
]
-
xshake
[
i2
][
1
];
s02
[
2
]
=
xshake
[
i0
][
2
]
-
xshake
[
i2
][
2
];
domain
->
minimum_image
(
s02
);
double
s03
[
3
];
s03
[
0
]
=
xshake
[
i0
][
0
]
-
xshake
[
i3
][
0
];
s03
[
1
]
=
xshake
[
i0
][
1
]
-
xshake
[
i3
][
1
];
s03
[
2
]
=
xshake
[
i0
][
2
]
-
xshake
[
i3
][
2
];
domain
->
minimum_image
(
s03
);
// scalar distances between atoms
double
r01sq
=
r01
[
0
]
*
r01
[
0
]
+
r01
[
1
]
*
r01
[
1
]
+
r01
[
2
]
*
r01
[
2
];
double
r02sq
=
r02
[
0
]
*
r02
[
0
]
+
r02
[
1
]
*
r02
[
1
]
+
r02
[
2
]
*
r02
[
2
];
double
r03sq
=
r03
[
0
]
*
r03
[
0
]
+
r03
[
1
]
*
r03
[
1
]
+
r03
[
2
]
*
r03
[
2
];
double
s01sq
=
s01
[
0
]
*
s01
[
0
]
+
s01
[
1
]
*
s01
[
1
]
+
s01
[
2
]
*
s01
[
2
];
double
s02sq
=
s02
[
0
]
*
s02
[
0
]
+
s02
[
1
]
*
s02
[
1
]
+
s02
[
2
]
*
s02
[
2
];
double
s03sq
=
s03
[
0
]
*
s03
[
0
]
+
s03
[
1
]
*
s03
[
1
]
+
s03
[
2
]
*
s03
[
2
];
// matrix coeffs and rhs for lamda equations
double
invmass0
=
1.0
/
mass
[
type
[
i0
]];
double
invmass1
=
1.0
/
mass
[
type
[
i1
]];
double
invmass2
=
1.0
/
mass
[
type
[
i2
]];
double
invmass3
=
1.0
/
mass
[
type
[
i3
]];
double
a11
=
2.0
*
(
invmass0
+
invmass1
)
*
(
s01
[
0
]
*
r01
[
0
]
+
s01
[
1
]
*
r01
[
1
]
+
s01
[
2
]
*
r01
[
2
]);
double
a12
=
2.0
*
invmass0
*
(
s01
[
0
]
*
r02
[
0
]
+
s01
[
1
]
*
r02
[
1
]
+
s01
[
2
]
*
r02
[
2
]);
double
a13
=
2.0
*
invmass0
*
(
s01
[
0
]
*
r03
[
0
]
+
s01
[
1
]
*
r03
[
1
]
+
s01
[
2
]
*
r03
[
2
]);
double
a21
=
2.0
*
invmass0
*
(
s02
[
0
]
*
r01
[
0
]
+
s02
[
1
]
*
r01
[
1
]
+
s02
[
2
]
*
r01
[
2
]);
double
a22
=
2.0
*
(
invmass0
+
invmass2
)
*
(
s02
[
0
]
*
r02
[
0
]
+
s02
[
1
]
*
r02
[
1
]
+
s02
[
2
]
*
r02
[
2
]);
double
a23
=
2.0
*
invmass0
*
(
s02
[
0
]
*
r03
[
0
]
+
s02
[
1
]
*
r03
[
1
]
+
s02
[
2
]
*
r03
[
2
]);
double
a31
=
2.0
*
invmass0
*
(
s03
[
0
]
*
r01
[
0
]
+
s03
[
1
]
*
r01
[
1
]
+
s03
[
2
]
*
r01
[
2
]);
double
a32
=
2.0
*
invmass0
*
(
s03
[
0
]
*
r02
[
0
]
+
s03
[
1
]
*
r02
[
1
]
+
s03
[
2
]
*
r02
[
2
]);
double
a33
=
2.0
*
(
invmass0
+
invmass3
)
*
(
s03
[
0
]
*
r03
[
0
]
+
s03
[
1
]
*
r03
[
1
]
+
s03
[
2
]
*
r03
[
2
]);
// inverse of matrix;
double
determ
=
a11
*
a22
*
a33
+
a12
*
a23
*
a31
+
a13
*
a21
*
a32
-
a11
*
a23
*
a32
-
a12
*
a21
*
a33
-
a13
*
a22
*
a31
;
if
(
determ
==
0.0
)
error
->
one
(
"Shake determinant = 0.0"
);
double
determinv
=
1.0
/
determ
;
double
a11inv
=
determinv
*
(
a22
*
a33
-
a23
*
a32
);
double
a12inv
=
-
determinv
*
(
a12
*
a33
-
a13
*
a32
);
double
a13inv
=
determinv
*
(
a12
*
a23
-
a13
*
a22
);
double
a21inv
=
-
determinv
*
(
a21
*
a33
-
a23
*
a31
);
double
a22inv
=
determinv
*
(
a11
*
a33
-
a13
*
a31
);
double
a23inv
=
-
determinv
*
(
a11
*
a23
-
a13
*
a21
);
double
a31inv
=
determinv
*
(
a21
*
a32
-
a22
*
a31
);
double
a32inv
=
-
determinv
*
(
a11
*
a32
-
a12
*
a31
);
double
a33inv
=
determinv
*
(
a11
*
a22
-
a12
*
a21
);
// quadratic correction coeffs
double
r0102
=
(
r01
[
0
]
*
r02
[
0
]
+
r01
[
1
]
*
r02
[
1
]
+
r01
[
2
]
*
r02
[
2
]);
double
r0103
=
(
r01
[
0
]
*
r03
[
0
]
+
r01
[
1
]
*
r03
[
1
]
+
r01
[
2
]
*
r03
[
2
]);
double
r0203
=
(
r02
[
0
]
*
r03
[
0
]
+
r02
[
1
]
*
r03
[
1
]
+
r02
[
2
]
*
r03
[
2
]);
double
quad1_0101
=
(
invmass0
+
invmass1
)
*
(
invmass0
+
invmass1
)
*
r01sq
;
double
quad1_0202
=
invmass0
*
invmass0
*
r02sq
;
double
quad1_0303
=
invmass0
*
invmass0
*
r03sq
;
double
quad1_0102
=
2.0
*
(
invmass0
+
invmass1
)
*
invmass0
*
r0102
;
double
quad1_0103
=
2.0
*
(
invmass0
+
invmass1
)
*
invmass0
*
r0103
;
double
quad1_0203
=
2.0
*
invmass0
*
invmass0
*
r0203
;
double
quad2_0101
=
invmass0
*
invmass0
*
r01sq
;
double
quad2_0202
=
(
invmass0
+
invmass2
)
*
(
invmass0
+
invmass2
)
*
r02sq
;
double
quad2_0303
=
invmass0
*
invmass0
*
r03sq
;
double
quad2_0102
=
2.0
*
(
invmass0
+
invmass2
)
*
invmass0
*
r0102
;
double
quad2_0103
=
2.0
*
invmass0
*
invmass0
*
r0103
;
double
quad2_0203
=
2.0
*
(
invmass0
+
invmass2
)
*
invmass0
*
r0203
;
double
quad3_0101
=
invmass0
*
invmass0
*
r01sq
;
double
quad3_0202
=
invmass0
*
invmass0
*
r02sq
;
double
quad3_0303
=
(
invmass0
+
invmass3
)
*
(
invmass0
+
invmass3
)
*
r03sq
;
double
quad3_0102
=
2.0
*
invmass0
*
invmass0
*
r0102
;
double
quad3_0103
=
2.0
*
(
invmass0
+
invmass3
)
*
invmass0
*
r0103
;
double
quad3_0203
=
2.0
*
(
invmass0
+
invmass3
)
*
invmass0
*
r0203
;
// iterate until converged
double
lamda01
=
0.0
;
double
lamda02
=
0.0
;
double
lamda03
=
0.0
;
int
niter
=
0
;
int
done
=
0
;
double
quad1
,
quad2
,
quad3
,
b1
,
b2
,
b3
,
lamda01_new
,
lamda02_new
,
lamda03_new
;
while
(
!
done
&&
niter
<
max_iter
)
{
quad1
=
quad1_0101
*
lamda01
*
lamda01
+
quad1_0202
*
lamda02
*
lamda02
+
quad1_0303
*
lamda03
*
lamda03
+
quad1_0102
*
lamda01
*
lamda02
+
quad1_0103
*
lamda01
*
lamda03
+
quad1_0203
*
lamda02
*
lamda03
;
quad2
=
quad2_0101
*
lamda01
*
lamda01
+
quad2_0202
*
lamda02
*
lamda02
+
quad2_0303
*
lamda03
*
lamda03
+
quad2_0102
*
lamda01
*
lamda02
+
quad2_0103
*
lamda01
*
lamda03
+
quad2_0203
*
lamda02
*
lamda03
;
quad3
=
quad3_0101
*
lamda01
*
lamda01
+
quad3_0202
*
lamda02
*
lamda02
+
quad3_0303
*
lamda03
*
lamda03
+
quad3_0102
*
lamda01
*
lamda02
+
quad3_0103
*
lamda01
*
lamda03
+
quad3_0203
*
lamda02
*
lamda03
;
b1
=
bond1
*
bond1
-
s01sq
-
quad1
;
b2
=
bond2
*
bond2
-
s02sq
-
quad2
;
b3
=
bond3
*
bond3
-
s03sq
-
quad3
;
lamda01_new
=
a11inv
*
b1
+
a12inv
*
b2
+
a13inv
*
b3
;
lamda02_new
=
a21inv
*
b1
+
a22inv
*
b2
+
a23inv
*
b3
;
lamda03_new
=
a31inv
*
b1
+
a32inv
*
b2
+
a33inv
*
b3
;
done
=
1
;
if
(
fabs
(
lamda01_new
-
lamda01
)
>
tolerance
)
done
=
0
;
if
(
fabs
(
lamda02_new
-
lamda02
)
>
tolerance
)
done
=
0
;
if
(
fabs
(
lamda03_new
-
lamda03
)
>
tolerance
)
done
=
0
;
lamda01
=
lamda01_new
;
lamda02
=
lamda02_new
;
lamda03
=
lamda03_new
;
niter
++
;
}
// update forces if atom is owned by this processor
lamda01
=
lamda01
/
dtfsq
;
lamda02
=
lamda02
/
dtfsq
;
lamda03
=
lamda03
/
dtfsq
;
if
(
i0
<
nlocal
)
{
f
[
i0
][
0
]
+=
lamda01
*
r01
[
0
]
+
lamda02
*
r02
[
0
]
+
lamda03
*
r03
[
0
];
f
[
i0
][
1
]
+=
lamda01
*
r01
[
1
]
+
lamda02
*
r02
[
1
]
+
lamda03
*
r03
[
1
];
f
[
i0
][
2
]
+=
lamda01
*
r01
[
2
]
+
lamda02
*
r02
[
2
]
+
lamda03
*
r03
[
2
];
}
if
(
i1
<
nlocal
)
{
f
[
i1
][
0
]
-=
lamda01
*
r01
[
0
];
f
[
i1
][
1
]
-=
lamda01
*
r01
[
1
];
f
[
i1
][
2
]
-=
lamda01
*
r01
[
2
];
}
if
(
i2
<
nlocal
)
{
f
[
i2
][
0
]
-=
lamda02
*
r02
[
0
];
f
[
i2
][
1
]
-=
lamda02
*
r02
[
1
];
f
[
i2
][
2
]
-=
lamda02
*
r02
[
2
];
}
if
(
i3
<
nlocal
)
{
f
[
i3
][
0
]
-=
lamda03
*
r03
[
0
];
f
[
i3
][
1
]
-=
lamda03
*
r03
[
1
];
f
[
i3
][
2
]
-=
lamda03
*
r03
[
2
];
}
if
(
vflag
)
{
int
factor
=
0
;
if
(
i0
<
nlocal
)
factor
++
;
if
(
i1
<
nlocal
)
factor
++
;
if
(
i2
<
nlocal
)
factor
++
;
if
(
i3
<
nlocal
)
factor
++
;
double
rfactor
=
0.25
*
factor
;
virial
[
0
]
+=
rfactor
*
lamda01
*
r01
[
0
]
*
r01
[
0
];
virial
[
1
]
+=
rfactor
*
lamda01
*
r01
[
1
]
*
r01
[
1
];
virial
[
2
]
+=
rfactor
*
lamda01
*
r01
[
2
]
*
r01
[
2
];
virial
[
3
]
+=
rfactor
*
lamda01
*
r01
[
0
]
*
r01
[
1
];
virial
[
4
]
+=
rfactor
*
lamda01
*
r01
[
0
]
*
r01
[
2
];
virial
[
5
]
+=
rfactor
*
lamda01
*
r01
[
1
]
*
r01
[
2
];
virial
[
0
]
+=
rfactor
*
lamda02
*
r02
[
0
]
*
r02
[
0
];
virial
[
1
]
+=
rfactor
*
lamda02
*
r02
[
1
]
*
r02
[
1
];
virial
[
2
]
+=
rfactor
*
lamda02
*
r02
[
2
]
*
r02
[
2
];
virial
[
3
]
+=
rfactor
*
lamda02
*
r02
[
0
]
*
r02
[
1
];
virial
[
4
]
+=
rfactor
*
lamda02
*
r02
[
0
]
*
r02
[
2
];
virial
[
5
]
+=
rfactor
*
lamda02
*
r02
[
1
]
*
r02
[
2
];
virial
[
0
]
+=
rfactor
*
lamda03
*
r03
[
0
]
*
r03
[
0
];
virial
[
1
]
+=
rfactor
*
lamda03
*
r03
[
1
]
*
r03
[
1
];
virial
[
2
]
+=
rfactor
*
lamda03
*
r03
[
2
]
*
r03
[
2
];
virial
[
3
]
+=
rfactor
*
lamda03
*
r03
[
0
]
*
r03
[
1
];
virial
[
4
]
+=
rfactor
*
lamda03
*
r03
[
0
]
*
r03
[
2
];
virial
[
5
]
+=
rfactor
*
lamda03
*
r03
[
1
]
*
r03
[
2
];
}
}
/* ---------------------------------------------------------------------- */
void
FixShake
::
shake3angle
(
int
m
)
{
// local atom IDs and constraint distances
int
i0
=
atom
->
map
(
shake_atom
[
m
][
0
]);
int
i1
=
atom
->
map
(
shake_atom
[
m
][
1
]);
int
i2
=
atom
->
map
(
shake_atom
[
m
][
2
]);
double
bond1
=
bond_distance
[
shake_type
[
m
][
0
]];
double
bond2
=
bond_distance
[
shake_type
[
m
][
1
]];
double
bond12
=
angle_distance
[
shake_type
[
m
][
2
]];
// r01,r02,r12 = distance vec between atoms, with PBC
double
r01
[
3
];
r01
[
0
]
=
x
[
i0
][
0
]
-
x
[
i1
][
0
];
r01
[
1
]
=
x
[
i0
][
1
]
-
x
[
i1
][
1
];
r01
[
2
]
=
x
[
i0
][
2
]
-
x
[
i1
][
2
];
domain
->
minimum_image
(
r01
);
double
r02
[
3
];
r02
[
0
]
=
x
[
i0
][
0
]
-
x
[
i2
][
0
];
r02
[
1
]
=
x
[
i0
][
1
]
-
x
[
i2
][
1
];
r02
[
2
]
=
x
[
i0
][
2
]
-
x
[
i2
][
2
];
domain
->
minimum_image
(
r02
);
double
r12
[
3
];
r12
[
0
]
=
x
[
i1
][
0
]
-
x
[
i2
][
0
];
r12
[
1
]
=
x
[
i1
][
1
]
-
x
[
i2
][
1
];
r12
[
2
]
=
x
[
i1
][
2
]
-
x
[
i2
][
2
];
domain
->
minimum_image
(
r12
);
// s01,s02,s12 = distance vec after unconstrained update, with PBC
double
s01
[
3
];
s01
[
0
]
=
xshake
[
i0
][
0
]
-
xshake
[
i1
][
0
];
s01
[
1
]
=
xshake
[
i0
][
1
]
-
xshake
[
i1
][
1
];
s01
[
2
]
=
xshake
[
i0
][
2
]
-
xshake
[
i1
][
2
];
domain
->
minimum_image
(
s01
);
double
s02
[
3
];
s02
[
0
]
=
xshake
[
i0
][
0
]
-
xshake
[
i2
][
0
];
s02
[
1
]
=
xshake
[
i0
][
1
]
-
xshake
[
i2
][
1
];
s02
[
2
]
=
xshake
[
i0
][
2
]
-
xshake
[
i2
][
2
];
domain
->
minimum_image
(
s02
);
double
s12
[
3
];
s12
[
0
]
=
xshake
[
i1
][
0
]
-
xshake
[
i2
][
0
];
s12
[
1
]
=
xshake
[
i1
][
1
]
-
xshake
[
i2
][
1
];
s12
[
2
]
=
xshake
[
i1
][
2
]
-
xshake
[
i2
][
2
];
domain
->
minimum_image
(
s12
);
// scalar distances between atoms
double
r01sq
=
r01
[
0
]
*
r01
[
0
]
+
r01
[
1
]
*
r01
[
1
]
+
r01
[
2
]
*
r01
[
2
];
double
r02sq
=
r02
[
0
]
*
r02
[
0
]
+
r02
[
1
]
*
r02
[
1
]
+
r02
[
2
]
*
r02
[
2
];
double
r12sq
=
r12
[
0
]
*
r12
[
0
]
+
r12
[
1
]
*
r12
[
1
]
+
r12
[
2
]
*
r12
[
2
];
double
s01sq
=
s01
[
0
]
*
s01
[
0
]
+
s01
[
1
]
*
s01
[
1
]
+
s01
[
2
]
*
s01
[
2
];
double
s02sq
=
s02
[
0
]
*
s02
[
0
]
+
s02
[
1
]
*
s02
[
1
]
+
s02
[
2
]
*
s02
[
2
];
double
s12sq
=
s12
[
0
]
*
s12
[
0
]
+
s12
[
1
]
*
s12
[
1
]
+
s12
[
2
]
*
s12
[
2
];
// matrix coeffs and rhs for lamda equations
double
invmass0
=
1.0
/
mass
[
type
[
i0
]];
double
invmass1
=
1.0
/
mass
[
type
[
i1
]];
double
invmass2
=
1.0
/
mass
[
type
[
i2
]];
double
a11
=
2.0
*
(
invmass0
+
invmass1
)
*
(
s01
[
0
]
*
r01
[
0
]
+
s01
[
1
]
*
r01
[
1
]
+
s01
[
2
]
*
r01
[
2
]);
double
a12
=
2.0
*
invmass0
*
(
s01
[
0
]
*
r02
[
0
]
+
s01
[
1
]
*
r02
[
1
]
+
s01
[
2
]
*
r02
[
2
]);
double
a13
=
-
2.0
*
invmass1
*
(
s01
[
0
]
*
r12
[
0
]
+
s01
[
1
]
*
r12
[
1
]
+
s01
[
2
]
*
r12
[
2
]);
double
a21
=
2.0
*
invmass0
*
(
s02
[
0
]
*
r01
[
0
]
+
s02
[
1
]
*
r01
[
1
]
+
s02
[
2
]
*
r01
[
2
]);
double
a22
=
2.0
*
(
invmass0
+
invmass2
)
*
(
s02
[
0
]
*
r02
[
0
]
+
s02
[
1
]
*
r02
[
1
]
+
s02
[
2
]
*
r02
[
2
]);
double
a23
=
2.0
*
invmass2
*
(
s02
[
0
]
*
r12
[
0
]
+
s02
[
1
]
*
r12
[
1
]
+
s02
[
2
]
*
r12
[
2
]);
double
a31
=
-
2.0
*
invmass1
*
(
s12
[
0
]
*
r01
[
0
]
+
s12
[
1
]
*
r01
[
1
]
+
s12
[
2
]
*
r01
[
2
]);
double
a32
=
2.0
*
invmass2
*
(
s12
[
0
]
*
r02
[
0
]
+
s12
[
1
]
*
r02
[
1
]
+
s12
[
2
]
*
r02
[
2
]);
double
a33
=
2.0
*
(
invmass1
+
invmass2
)
*
(
s12
[
0
]
*
r12
[
0
]
+
s12
[
1
]
*
r12
[
1
]
+
s12
[
2
]
*
r12
[
2
]);
// inverse of matrix
double
determ
=
a11
*
a22
*
a33
+
a12
*
a23
*
a31
+
a13
*
a21
*
a32
-
a11
*
a23
*
a32
-
a12
*
a21
*
a33
-
a13
*
a22
*
a31
;
if
(
determ
==
0.0
)
error
->
one
(
"Shake determinant = 0.0"
);
double
determinv
=
1.0
/
determ
;
double
a11inv
=
determinv
*
(
a22
*
a33
-
a23
*
a32
);
double
a12inv
=
-
determinv
*
(
a12
*
a33
-
a13
*
a32
);
double
a13inv
=
determinv
*
(
a12
*
a23
-
a13
*
a22
);
double
a21inv
=
-
determinv
*
(
a21
*
a33
-
a23
*
a31
);
double
a22inv
=
determinv
*
(
a11
*
a33
-
a13
*
a31
);
double
a23inv
=
-
determinv
*
(
a11
*
a23
-
a13
*
a21
);
double
a31inv
=
determinv
*
(
a21
*
a32
-
a22
*
a31
);
double
a32inv
=
-
determinv
*
(
a11
*
a32
-
a12
*
a31
);
double
a33inv
=
determinv
*
(
a11
*
a22
-
a12
*
a21
);
// quadratic correction coeffs
double
r0102
=
(
r01
[
0
]
*
r02
[
0
]
+
r01
[
1
]
*
r02
[
1
]
+
r01
[
2
]
*
r02
[
2
]);
double
r0112
=
(
r01
[
0
]
*
r12
[
0
]
+
r01
[
1
]
*
r12
[
1
]
+
r01
[
2
]
*
r12
[
2
]);
double
r0212
=
(
r02
[
0
]
*
r12
[
0
]
+
r02
[
1
]
*
r12
[
1
]
+
r02
[
2
]
*
r12
[
2
]);
double
quad1_0101
=
(
invmass0
+
invmass1
)
*
(
invmass0
+
invmass1
)
*
r01sq
;
double
quad1_0202
=
invmass0
*
invmass0
*
r02sq
;
double
quad1_1212
=
invmass1
*
invmass1
*
r12sq
;
double
quad1_0102
=
2.0
*
(
invmass0
+
invmass1
)
*
invmass0
*
r0102
;
double
quad1_0112
=
-
2.0
*
(
invmass0
+
invmass1
)
*
invmass1
*
r0112
;
double
quad1_0212
=
-
2.0
*
invmass0
*
invmass1
*
r0212
;
double
quad2_0101
=
invmass0
*
invmass0
*
r01sq
;
double
quad2_0202
=
(
invmass0
+
invmass2
)
*
(
invmass0
+
invmass2
)
*
r02sq
;
double
quad2_1212
=
invmass2
*
invmass2
*
r12sq
;
double
quad2_0102
=
2.0
*
(
invmass0
+
invmass2
)
*
invmass0
*
r0102
;
double
quad2_0112
=
2.0
*
invmass0
*
invmass2
*
r0112
;
double
quad2_0212
=
2.0
*
(
invmass0
+
invmass2
)
*
invmass2
*
r0212
;
double
quad3_0101
=
invmass1
*
invmass1
*
r01sq
;
double
quad3_0202
=
invmass2
*
invmass2
*
r02sq
;
double
quad3_1212
=
(
invmass1
+
invmass2
)
*
(
invmass1
+
invmass2
)
*
r12sq
;
double
quad3_0102
=
-
2.0
*
invmass1
*
invmass2
*
r0102
;
double
quad3_0112
=
-
2.0
*
(
invmass1
+
invmass2
)
*
invmass1
*
r0112
;
double
quad3_0212
=
2.0
*
(
invmass1
+
invmass2
)
*
invmass2
*
r0212
;
// iterate until converged
double
lamda01
=
0.0
;
double
lamda02
=
0.0
;
double
lamda12
=
0.0
;
int
niter
=
0
;
int
done
=
0
;
double
quad1
,
quad2
,
quad3
,
b1
,
b2
,
b3
,
lamda01_new
,
lamda02_new
,
lamda12_new
;
while
(
!
done
&&
niter
<
max_iter
)
{
quad1
=
quad1_0101
*
lamda01
*
lamda01
+
quad1_0202
*
lamda02
*
lamda02
+
quad1_1212
*
lamda12
*
lamda12
+
quad1_0102
*
lamda01
*
lamda02
+
quad1_0112
*
lamda01
*
lamda12
+
quad1_0212
*
lamda02
*
lamda12
;
quad2
=
quad2_0101
*
lamda01
*
lamda01
+
quad2_0202
*
lamda02
*
lamda02
+
quad2_1212
*
lamda12
*
lamda12
+
quad2_0102
*
lamda01
*
lamda02
+
quad2_0112
*
lamda01
*
lamda12
+
quad2_0212
*
lamda02
*
lamda12
;
quad3
=
quad3_0101
*
lamda01
*
lamda01
+
quad3_0202
*
lamda02
*
lamda02
+
quad3_1212
*
lamda12
*
lamda12
+
quad3_0102
*
lamda01
*
lamda02
+
quad3_0112
*
lamda01
*
lamda12
+
quad3_0212
*
lamda02
*
lamda12
;
b1
=
bond1
*
bond1
-
s01sq
-
quad1
;
b2
=
bond2
*
bond2
-
s02sq
-
quad2
;
b3
=
bond12
*
bond12
-
s12sq
-
quad3
;
lamda01_new
=
a11inv
*
b1
+
a12inv
*
b2
+
a13inv
*
b3
;
lamda02_new
=
a21inv
*
b1
+
a22inv
*
b2
+
a23inv
*
b3
;
lamda12_new
=
a31inv
*
b1
+
a32inv
*
b2
+
a33inv
*
b3
;
done
=
1
;
if
(
fabs
(
lamda01_new
-
lamda01
)
>
tolerance
)
done
=
0
;
if
(
fabs
(
lamda02_new
-
lamda02
)
>
tolerance
)
done
=
0
;
if
(
fabs
(
lamda12_new
-
lamda12
)
>
tolerance
)
done
=
0
;
lamda01
=
lamda01_new
;
lamda02
=
lamda02_new
;
lamda12
=
lamda12_new
;
niter
++
;
}
// update forces if atom is owned by this processor
lamda01
=
lamda01
/
dtfsq
;
lamda02
=
lamda02
/
dtfsq
;
lamda12
=
lamda12
/
dtfsq
;
if
(
i0
<
nlocal
)
{
f
[
i0
][
0
]
+=
lamda01
*
r01
[
0
]
+
lamda02
*
r02
[
0
];
f
[
i0
][
1
]
+=
lamda01
*
r01
[
1
]
+
lamda02
*
r02
[
1
];
f
[
i0
][
2
]
+=
lamda01
*
r01
[
2
]
+
lamda02
*
r02
[
2
];
}
if
(
i1
<
nlocal
)
{
f
[
i1
][
0
]
-=
lamda01
*
r01
[
0
]
-
lamda12
*
r12
[
0
];
f
[
i1
][
1
]
-=
lamda01
*
r01
[
1
]
-
lamda12
*
r12
[
1
];
f
[
i1
][
2
]
-=
lamda01
*
r01
[
2
]
-
lamda12
*
r12
[
2
];
}
if
(
i2
<
nlocal
)
{
f
[
i2
][
0
]
-=
lamda02
*
r02
[
0
]
+
lamda12
*
r12
[
0
];
f
[
i2
][
1
]
-=
lamda02
*
r02
[
1
]
+
lamda12
*
r12
[
1
];
f
[
i2
][
2
]
-=
lamda02
*
r02
[
2
]
+
lamda12
*
r12
[
2
];
}
if
(
vflag
)
{
int
factor
=
0
;
if
(
i0
<
nlocal
)
factor
++
;
if
(
i1
<
nlocal
)
factor
++
;
if
(
i2
<
nlocal
)
factor
++
;
double
rfactor
=
factor
/
3.0
;
virial
[
0
]
+=
rfactor
*
lamda01
*
r01
[
0
]
*
r01
[
0
];
virial
[
1
]
+=
rfactor
*
lamda01
*
r01
[
1
]
*
r01
[
1
];
virial
[
2
]
+=
rfactor
*
lamda01
*
r01
[
2
]
*
r01
[
2
];
virial
[
3
]
+=
rfactor
*
lamda01
*
r01
[
0
]
*
r01
[
1
];
virial
[
4
]
+=
rfactor
*
lamda01
*
r01
[
0
]
*
r01
[
2
];
virial
[
5
]
+=
rfactor
*
lamda01
*
r01
[
1
]
*
r01
[
2
];
virial
[
0
]
+=
rfactor
*
lamda02
*
r02
[
0
]
*
r02
[
0
];
virial
[
1
]
+=
rfactor
*
lamda02
*
r02
[
1
]
*
r02
[
1
];
virial
[
2
]
+=
rfactor
*
lamda02
*
r02
[
2
]
*
r02
[
2
];
virial
[
3
]
+=
rfactor
*
lamda02
*
r02
[
0
]
*
r02
[
1
];
virial
[
4
]
+=
rfactor
*
lamda02
*
r02
[
0
]
*
r02
[
2
];
virial
[
5
]
+=
rfactor
*
lamda02
*
r02
[
1
]
*
r02
[
2
];
virial
[
0
]
+=
rfactor
*
lamda12
*
r12
[
0
]
*
r12
[
0
];
virial
[
1
]
+=
rfactor
*
lamda12
*
r12
[
1
]
*
r12
[
1
];
virial
[
2
]
+=
rfactor
*
lamda12
*
r12
[
2
]
*
r12
[
2
];
virial
[
3
]
+=
rfactor
*
lamda12
*
r12
[
0
]
*
r12
[
1
];
virial
[
4
]
+=
rfactor
*
lamda12
*
r12
[
0
]
*
r12
[
2
];
virial
[
5
]
+=
rfactor
*
lamda12
*
r12
[
1
]
*
r12
[
2
];
}
}
/* ----------------------------------------------------------------------
print-out bond & angle statistics
------------------------------------------------------------------------- */
void
FixShake
::
stats
()
{
int
i
,
j
,
m
,
n
,
iatom
,
jatom
,
katom
;
double
delx
,
dely
,
delz
;
double
r
,
r1
,
r2
,
r3
,
angle
;
// zero out accumulators
int
nb
=
atom
->
nbondtypes
+
1
;
int
na
=
atom
->
nangletypes
+
1
;
for
(
i
=
0
;
i
<
nb
;
i
++
)
{
b_count
[
i
]
=
0
;
b_ave
[
i
]
=
b_max
[
i
]
=
0.0
;
b_min
[
i
]
=
BIG
;
}
for
(
i
=
0
;
i
<
na
;
i
++
)
{
a_count
[
i
]
=
0
;
a_ave
[
i
]
=
a_max
[
i
]
=
0.0
;
a_min
[
i
]
=
BIG
;
}
// log stats for each bond & angle
// OK to double count since are just averaging
double
**
x
=
atom
->
x
;
int
nlocal
=
atom
->
nlocal
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
shake_flag
[
i
]
==
0
)
continue
;
// bond stats
n
=
shake_flag
[
i
];
if
(
n
==
1
)
n
=
3
;
iatom
=
atom
->
map
(
shake_atom
[
i
][
0
]);
for
(
j
=
1
;
j
<
n
;
j
++
)
{
jatom
=
atom
->
map
(
shake_atom
[
i
][
j
]);
delx
=
x
[
iatom
][
0
]
-
x
[
jatom
][
0
];
dely
=
x
[
iatom
][
1
]
-
x
[
jatom
][
1
];
delz
=
x
[
iatom
][
2
]
-
x
[
jatom
][
2
];
domain
->
minimum_image
(
&
delx
,
&
dely
,
&
delz
);
r
=
sqrt
(
delx
*
delx
+
dely
*
dely
+
delz
*
delz
);
m
=
shake_type
[
i
][
j
-
1
];
b_count
[
m
]
++
;
b_ave
[
m
]
+=
r
;
b_max
[
m
]
=
MAX
(
b_max
[
m
],
r
);
b_min
[
m
]
=
MIN
(
b_min
[
m
],
r
);
}
// angle stats
if
(
shake_flag
[
i
]
==
1
)
{
iatom
=
atom
->
map
(
shake_atom
[
i
][
0
]);
jatom
=
atom
->
map
(
shake_atom
[
i
][
1
]);
katom
=
atom
->
map
(
shake_atom
[
i
][
2
]);
delx
=
x
[
iatom
][
0
]
-
x
[
jatom
][
0
];
dely
=
x
[
iatom
][
1
]
-
x
[
jatom
][
1
];
delz
=
x
[
iatom
][
2
]
-
x
[
jatom
][
2
];
domain
->
minimum_image
(
&
delx
,
&
dely
,
&
delz
);
r1
=
sqrt
(
delx
*
delx
+
dely
*
dely
+
delz
*
delz
);
delx
=
x
[
iatom
][
0
]
-
x
[
katom
][
0
];
dely
=
x
[
iatom
][
1
]
-
x
[
katom
][
1
];
delz
=
x
[
iatom
][
2
]
-
x
[
katom
][
2
];
domain
->
minimum_image
(
&
delx
,
&
dely
,
&
delz
);
r2
=
sqrt
(
delx
*
delx
+
dely
*
dely
+
delz
*
delz
);
delx
=
x
[
jatom
][
0
]
-
x
[
katom
][
0
];
dely
=
x
[
jatom
][
1
]
-
x
[
katom
][
1
];
delz
=
x
[
jatom
][
2
]
-
x
[
katom
][
2
];
domain
->
minimum_image
(
&
delx
,
&
dely
,
&
delz
);
r3
=
sqrt
(
delx
*
delx
+
dely
*
dely
+
delz
*
delz
);
angle
=
acos
((
r1
*
r1
+
r2
*
r2
-
r3
*
r3
)
/
(
2.0
*
r1
*
r2
));
angle
*=
180.0
/
PI
;
m
=
shake_type
[
i
][
2
];
a_count
[
m
]
++
;
a_ave
[
m
]
+=
angle
;
a_max
[
m
]
=
MAX
(
a_max
[
m
],
angle
);
a_min
[
m
]
=
MIN
(
a_min
[
m
],
angle
);
}
}
// sum across all procs
MPI_Allreduce
(
b_count
,
b_count_all
,
nb
,
MPI_INT
,
MPI_SUM
,
world
);
MPI_Allreduce
(
b_ave
,
b_ave_all
,
nb
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
MPI_Allreduce
(
b_max
,
b_max_all
,
nb
,
MPI_DOUBLE
,
MPI_MAX
,
world
);
MPI_Allreduce
(
b_min
,
b_min_all
,
nb
,
MPI_DOUBLE
,
MPI_MIN
,
world
);
MPI_Allreduce
(
a_count
,
a_count_all
,
na
,
MPI_INT
,
MPI_SUM
,
world
);
MPI_Allreduce
(
a_ave
,
a_ave_all
,
na
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
MPI_Allreduce
(
a_max
,
a_max_all
,
na
,
MPI_DOUBLE
,
MPI_MAX
,
world
);
MPI_Allreduce
(
a_min
,
a_min_all
,
na
,
MPI_DOUBLE
,
MPI_MIN
,
world
);
// print stats only for non-zero counts
if
(
me
==
0
)
{
if
(
screen
)
{
fprintf
(
screen
,
"SHAKE stats (type/ave/delta) on step %d
\n
"
,
update
->
ntimestep
);
for
(
i
=
1
;
i
<
nb
;
i
++
)
if
(
b_count_all
[
i
])
fprintf
(
screen
,
" %d %g %g
\n
"
,
i
,
b_ave_all
[
i
]
/
b_count_all
[
i
],
b_max_all
[
i
]
-
b_min_all
[
i
]);
for
(
i
=
1
;
i
<
na
;
i
++
)
if
(
a_count_all
[
i
])
fprintf
(
screen
,
" %d %g %g
\n
"
,
i
,
a_ave_all
[
i
]
/
a_count_all
[
i
],
a_max_all
[
i
]
-
a_min_all
[
i
]);
}
if
(
logfile
)
{
fprintf
(
logfile
,
"SHAKE stats (type/ave/delta) on step %d
\n
"
,
update
->
ntimestep
);
for
(
i
=
0
;
i
<
nb
;
i
++
)
if
(
b_count_all
[
i
])
fprintf
(
logfile
,
" %d %g %g
\n
"
,
i
,
b_ave_all
[
i
]
/
b_count_all
[
i
],
b_max_all
[
i
]
-
b_min_all
[
i
]);
for
(
i
=
0
;
i
<
na
;
i
++
)
if
(
a_count_all
[
i
])
fprintf
(
logfile
,
" %d %g %g
\n
"
,
i
,
a_ave_all
[
i
]
/
a_count_all
[
i
],
a_max_all
[
i
]
-
a_min_all
[
i
]);
}
}
// next timestep for stats
next_output
+=
output_every
;
}
/* ----------------------------------------------------------------------
find a bond between global tags n1 and n2 stored with local atom i
return -1 if don't find it
return bond index if do find it
------------------------------------------------------------------------- */
int
FixShake
::
bondfind
(
int
i
,
int
n1
,
int
n2
)
{
int
*
tag
=
atom
->
tag
;
int
**
bond_atom
=
atom
->
bond_atom
;
int
nbonds
=
atom
->
num_bond
[
i
];
int
m
;
for
(
m
=
0
;
m
<
nbonds
;
m
++
)
{
if
(
n1
==
tag
[
i
]
&&
n2
==
bond_atom
[
i
][
m
])
break
;
if
(
n1
==
bond_atom
[
i
][
m
]
&&
n2
==
tag
[
i
])
break
;
}
if
(
m
<
nbonds
)
return
m
;
return
-
1
;
}
/* ----------------------------------------------------------------------
find an angle with global end atoms n1 and n2 stored with local atom i
return -1 if don't find it
return angle index if do find it
------------------------------------------------------------------------- */
int
FixShake
::
anglefind
(
int
i
,
int
n1
,
int
n2
)
{
int
**
angle_atom1
=
atom
->
angle_atom1
;
int
**
angle_atom3
=
atom
->
angle_atom3
;
int
nangles
=
atom
->
num_angle
[
i
];
int
m
;
for
(
m
=
0
;
m
<
nangles
;
m
++
)
{
if
(
n1
==
angle_atom1
[
i
][
m
]
&&
n2
==
angle_atom3
[
i
][
m
])
break
;
if
(
n1
==
angle_atom3
[
i
][
m
]
&&
n2
==
angle_atom1
[
i
][
m
])
break
;
}
if
(
m
<
nangles
)
return
m
;
return
-
1
;
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
int
FixShake
::
memory_usage
()
{
int
nmax
=
atom
->
nmax
;
int
bytes
=
nmax
*
sizeof
(
int
);
bytes
+=
nmax
*
4
*
sizeof
(
int
);
bytes
+=
nmax
*
3
*
sizeof
(
int
);
bytes
+=
nmax
*
3
*
sizeof
(
double
);
return
bytes
;
}
/* ----------------------------------------------------------------------
allocate local atom-based arrays
------------------------------------------------------------------------- */
void
FixShake
::
grow_arrays
(
int
nmax
)
{
shake_flag
=
(
int
*
)
memory
->
srealloc
(
shake_flag
,
nmax
*
sizeof
(
int
),
"shake:shake_flag"
);
shake_atom
=
memory
->
grow_2d_int_array
(
shake_atom
,
nmax
,
4
,
"shake:shake_atom"
);
shake_type
=
memory
->
grow_2d_int_array
(
shake_type
,
nmax
,
3
,
"shake:shake_type"
);
memory
->
destroy_2d_double_array
(
xshake
);
xshake
=
memory
->
create_2d_double_array
(
nmax
,
3
,
"shake:xshake"
);
}
/* ----------------------------------------------------------------------
copy values within local atom-based arrays
------------------------------------------------------------------------- */
void
FixShake
::
copy_arrays
(
int
i
,
int
j
)
{
int
flag
=
shake_flag
[
j
]
=
shake_flag
[
i
];
if
(
flag
==
1
)
{
shake_atom
[
j
][
0
]
=
shake_atom
[
i
][
0
];
shake_atom
[
j
][
1
]
=
shake_atom
[
i
][
1
];
shake_atom
[
j
][
2
]
=
shake_atom
[
i
][
2
];
shake_type
[
j
][
0
]
=
shake_type
[
i
][
0
];
shake_type
[
j
][
1
]
=
shake_type
[
i
][
1
];
shake_type
[
j
][
2
]
=
shake_type
[
i
][
2
];
}
else
if
(
flag
==
2
)
{
shake_atom
[
j
][
0
]
=
shake_atom
[
i
][
0
];
shake_atom
[
j
][
1
]
=
shake_atom
[
i
][
1
];
shake_type
[
j
][
0
]
=
shake_type
[
i
][
0
];
}
else
if
(
flag
==
3
)
{
shake_atom
[
j
][
0
]
=
shake_atom
[
i
][
0
];
shake_atom
[
j
][
1
]
=
shake_atom
[
i
][
1
];
shake_atom
[
j
][
2
]
=
shake_atom
[
i
][
2
];
shake_type
[
j
][
0
]
=
shake_type
[
i
][
0
];
shake_type
[
j
][
1
]
=
shake_type
[
i
][
1
];
}
else
if
(
flag
==
4
)
{
shake_atom
[
j
][
0
]
=
shake_atom
[
i
][
0
];
shake_atom
[
j
][
1
]
=
shake_atom
[
i
][
1
];
shake_atom
[
j
][
2
]
=
shake_atom
[
i
][
2
];
shake_atom
[
j
][
3
]
=
shake_atom
[
i
][
3
];
shake_type
[
j
][
0
]
=
shake_type
[
i
][
0
];
shake_type
[
j
][
1
]
=
shake_type
[
i
][
1
];
shake_type
[
j
][
2
]
=
shake_type
[
i
][
2
];
}
}
/* ----------------------------------------------------------------------
pack values in local atom-based arrays for exchange with another proc
------------------------------------------------------------------------- */
int
FixShake
::
pack_exchange
(
int
i
,
double
*
buf
)
{
int
m
=
0
;
buf
[
m
++
]
=
shake_flag
[
i
];
int
flag
=
shake_flag
[
i
];
if
(
flag
==
1
)
{
buf
[
m
++
]
=
shake_atom
[
i
][
0
];
buf
[
m
++
]
=
shake_atom
[
i
][
1
];
buf
[
m
++
]
=
shake_atom
[
i
][
2
];
buf
[
m
++
]
=
shake_type
[
i
][
0
];
buf
[
m
++
]
=
shake_type
[
i
][
1
];
buf
[
m
++
]
=
shake_type
[
i
][
2
];
}
else
if
(
flag
==
2
)
{
buf
[
m
++
]
=
shake_atom
[
i
][
0
];
buf
[
m
++
]
=
shake_atom
[
i
][
1
];
buf
[
m
++
]
=
shake_type
[
i
][
0
];
}
else
if
(
flag
==
3
)
{
buf
[
m
++
]
=
shake_atom
[
i
][
0
];
buf
[
m
++
]
=
shake_atom
[
i
][
1
];
buf
[
m
++
]
=
shake_atom
[
i
][
2
];
buf
[
m
++
]
=
shake_type
[
i
][
0
];
buf
[
m
++
]
=
shake_type
[
i
][
1
];
}
else
if
(
flag
==
4
)
{
buf
[
m
++
]
=
shake_atom
[
i
][
0
];
buf
[
m
++
]
=
shake_atom
[
i
][
1
];
buf
[
m
++
]
=
shake_atom
[
i
][
2
];
buf
[
m
++
]
=
shake_atom
[
i
][
3
];
buf
[
m
++
]
=
shake_type
[
i
][
0
];
buf
[
m
++
]
=
shake_type
[
i
][
1
];
buf
[
m
++
]
=
shake_type
[
i
][
2
];
}
return
m
;
}
/* ----------------------------------------------------------------------
unpack values in local atom-based arrays from exchange with another proc
------------------------------------------------------------------------- */
int
FixShake
::
unpack_exchange
(
int
nlocal
,
double
*
buf
)
{
int
m
=
0
;
int
flag
=
shake_flag
[
nlocal
]
=
static_cast
<
int
>
(
buf
[
m
++
]);
if
(
flag
==
1
)
{
shake_atom
[
nlocal
][
0
]
=
static_cast
<
int
>
(
buf
[
m
++
]);
shake_atom
[
nlocal
][
1
]
=
static_cast
<
int
>
(
buf
[
m
++
]);
shake_atom
[
nlocal
][
2
]
=
static_cast
<
int
>
(
buf
[
m
++
]);
shake_type
[
nlocal
][
0
]
=
static_cast
<
int
>
(
buf
[
m
++
]);
shake_type
[
nlocal
][
1
]
=
static_cast
<
int
>
(
buf
[
m
++
]);
shake_type
[
nlocal
][
2
]
=
static_cast
<
int
>
(
buf
[
m
++
]);
}
else
if
(
flag
==
2
)
{
shake_atom
[
nlocal
][
0
]
=
static_cast
<
int
>
(
buf
[
m
++
]);
shake_atom
[
nlocal
][
1
]
=
static_cast
<
int
>
(
buf
[
m
++
]);
shake_type
[
nlocal
][
0
]
=
static_cast
<
int
>
(
buf
[
m
++
]);
}
else
if
(
flag
==
3
)
{
shake_atom
[
nlocal
][
0
]
=
static_cast
<
int
>
(
buf
[
m
++
]);
shake_atom
[
nlocal
][
1
]
=
static_cast
<
int
>
(
buf
[
m
++
]);
shake_atom
[
nlocal
][
2
]
=
static_cast
<
int
>
(
buf
[
m
++
]);
shake_type
[
nlocal
][
0
]
=
static_cast
<
int
>
(
buf
[
m
++
]);
shake_type
[
nlocal
][
1
]
=
static_cast
<
int
>
(
buf
[
m
++
]);
}
else
if
(
flag
==
4
)
{
shake_atom
[
nlocal
][
0
]
=
static_cast
<
int
>
(
buf
[
m
++
]);
shake_atom
[
nlocal
][
1
]
=
static_cast
<
int
>
(
buf
[
m
++
]);
shake_atom
[
nlocal
][
2
]
=
static_cast
<
int
>
(
buf
[
m
++
]);
shake_atom
[
nlocal
][
3
]
=
static_cast
<
int
>
(
buf
[
m
++
]);
shake_type
[
nlocal
][
0
]
=
static_cast
<
int
>
(
buf
[
m
++
]);
shake_type
[
nlocal
][
1
]
=
static_cast
<
int
>
(
buf
[
m
++
]);
shake_type
[
nlocal
][
2
]
=
static_cast
<
int
>
(
buf
[
m
++
]);
}
return
m
;
}
/* ----------------------------------------------------------------------
enforce SHAKE constraints from rRESPA
prediction portion is different than Verlet
rRESPA updating of atom coords is done with full v, but only portions of f
------------------------------------------------------------------------- */
void
FixShake
::
post_force_respa
(
int
vflag_in
,
int
ilevel
,
int
iloop
)
{
// call stats only on outermost level
if
(
ilevel
==
nlevels_respa
-
1
&&
update
->
ntimestep
==
next_output
)
stats
();
// perform SHAKE on every loop iteration of every rRESPA level
// except last loop iteration of inner levels
if
(
ilevel
<
nlevels_respa
-
1
&&
iloop
==
loop_respa
[
ilevel
]
-
1
)
return
;
// xshake = atom coords after next x update in innermost loop
// depends on rRESPA level
// for levels > 0 this includes more than one velocity update
// xshake = predicted position from call to this routine at level N =
// x + dt0 (v + dtN/m fN + 1/2 dt(N-1)/m f(N-1) + ... + 1/2 dt0/m f0)
double
***
f_level
=
((
FixRespa
*
)
modify
->
fix
[
ifix_respa
])
->
f_level
;
dtfsq
=
dtf_inner
*
step_respa
[
ilevel
];
double
invmass
,
dtfmsq
;
int
jlevel
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
shake_flag
[
i
])
{
invmass
=
1.0
/
mass
[
type
[
i
]];
dtfmsq
=
dtfsq
*
invmass
;
xshake
[
i
][
0
]
=
x
[
i
][
0
]
+
dtv
*
v
[
i
][
0
]
+
dtfmsq
*
f
[
i
][
0
];
xshake
[
i
][
1
]
=
x
[
i
][
1
]
+
dtv
*
v
[
i
][
1
]
+
dtfmsq
*
f
[
i
][
1
];
xshake
[
i
][
2
]
=
x
[
i
][
2
]
+
dtv
*
v
[
i
][
2
]
+
dtfmsq
*
f
[
i
][
2
];
for
(
jlevel
=
0
;
jlevel
<
ilevel
;
jlevel
++
)
{
dtfmsq
=
dtf_innerhalf
*
step_respa
[
jlevel
]
*
invmass
;
xshake
[
i
][
0
]
+=
dtfmsq
*
f_level
[
i
][
jlevel
][
0
];
xshake
[
i
][
1
]
+=
dtfmsq
*
f_level
[
i
][
jlevel
][
1
];
xshake
[
i
][
2
]
+=
dtfmsq
*
f_level
[
i
][
jlevel
][
2
];
}
}
else
xshake
[
i
][
2
]
=
xshake
[
i
][
1
]
=
xshake
[
i
][
0
]
=
0.0
;
}
// communicate results if necessary
if
(
nprocs
>
1
)
comm
->
comm_fix
(
this
);
// zero out SHAKE contribution to virial
vflag
=
vflag_in
;
if
(
vflag
)
for
(
int
n
=
0
;
n
<
6
;
n
++
)
virial
[
n
]
=
0.0
;
// loop over clusters
int
m
;
for
(
int
i
=
0
;
i
<
nlist
;
i
++
)
{
m
=
list
[
i
];
if
(
shake_flag
[
m
]
==
2
)
shake2
(
m
);
else
if
(
shake_flag
[
m
]
==
3
)
shake3
(
m
);
else
if
(
shake_flag
[
m
]
==
4
)
shake4
(
m
);
else
shake3angle
(
m
);
}
}
/* ---------------------------------------------------------------------- */
int
FixShake
::
pack_comm
(
int
n
,
int
*
list
,
double
*
buf
,
int
*
pbc_flags
)
{
int
i
,
j
,
m
;
m
=
0
;
if
(
pbc_flags
[
0
]
==
0
)
{
for
(
i
=
0
;
i
<
n
;
i
++
)
{
j
=
list
[
i
];
buf
[
m
++
]
=
xshake
[
j
][
0
];
buf
[
m
++
]
=
xshake
[
j
][
1
];
buf
[
m
++
]
=
xshake
[
j
][
2
];
}
}
else
{
for
(
i
=
0
;
i
<
n
;
i
++
)
{
j
=
list
[
i
];
buf
[
m
++
]
=
xshake
[
j
][
0
]
+
pbc_flags
[
1
]
*
domain
->
xprd
;
buf
[
m
++
]
=
xshake
[
j
][
1
]
+
pbc_flags
[
2
]
*
domain
->
yprd
;
buf
[
m
++
]
=
xshake
[
j
][
2
]
+
pbc_flags
[
3
]
*
domain
->
zprd
;
}
}
return
3
;
}
/* ---------------------------------------------------------------------- */
void
FixShake
::
unpack_comm
(
int
n
,
int
first
,
double
*
buf
)
{
int
i
,
m
,
last
;
m
=
0
;
last
=
first
+
n
;
for
(
i
=
first
;
i
<
last
;
i
++
)
{
xshake
[
i
][
0
]
=
buf
[
m
++
];
xshake
[
i
][
1
]
=
buf
[
m
++
];
xshake
[
i
][
2
]
=
buf
[
m
++
];
}
}
Event Timeline
Log In to Comment