Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F78533929
fix_rigid_nh.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, Aug 21, 11:30
Size
40 KB
Mime Type
text/x-c
Expires
Fri, Aug 23, 11:30 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20060908
Attached To
rLAMMPS lammps
fix_rigid_nh.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: Tony Sheh (U Michigan), Trung Dac Nguyen (U Michigan)
references: Kamberaj et al., J. Chem. Phys. 122, 224114 (2005)
Miller et al., J Chem Phys. 116, 8649-8659 (2002)
------------------------------------------------------------------------- */
#include "math.h"
#include "stdio.h"
#include "string.h"
#include "fix_rigid_nh.h"
#include "math_extra.h"
#include "atom.h"
#include "compute.h"
#include "domain.h"
#include "update.h"
#include "modify.h"
#include "fix_deform.h"
#include "group.h"
#include "comm.h"
#include "force.h"
#include "kspace.h"
#include "output.h"
#include "memory.h"
#include "error.h"
using
namespace
LAMMPS_NS
;
using
namespace
FixConst
;
enum
{
NONE
,
XYZ
,
XY
,
YZ
,
XZ
};
// same as in FixRigid
enum
{
ISO
,
ANISO
,
TRICLINIC
};
// same as in FixRigid
#define EPSILON 1.0e-7
/* ---------------------------------------------------------------------- */
FixRigidNH
::
FixRigidNH
(
LAMMPS
*
lmp
,
int
narg
,
char
**
arg
)
:
FixRigid
(
lmp
,
narg
,
arg
)
{
// error checks: could be moved up to FixRigid
if
((
p_flag
[
0
]
==
1
&&
p_period
[
0
]
<=
0.0
)
||
(
p_flag
[
1
]
==
1
&&
p_period
[
1
]
<=
0.0
)
||
(
p_flag
[
2
]
==
1
&&
p_period
[
2
]
<=
0.0
))
error
->
all
(
FLERR
,
"Fix rigid npt/nph period must be > 0.0"
);
if
(
dimension
==
2
&&
p_flag
[
2
])
error
->
all
(
FLERR
,
"Invalid fix rigid npt/nph command for a 2d simulation"
);
if
(
dimension
==
2
&&
(
pcouple
==
YZ
||
pcouple
==
XZ
))
error
->
all
(
FLERR
,
"Invalid fix rigid npt/nph command for a 2d simulation"
);
if
(
pcouple
==
XYZ
&&
(
p_flag
[
0
]
==
0
||
p_flag
[
1
]
==
0
))
error
->
all
(
FLERR
,
"Invalid fix rigid npt/nph command pressure settings"
);
if
(
pcouple
==
XYZ
&&
dimension
==
3
&&
p_flag
[
2
]
==
0
)
error
->
all
(
FLERR
,
"Invalid fix rigid npt/nph command pressure settings"
);
if
(
pcouple
==
XY
&&
(
p_flag
[
0
]
==
0
||
p_flag
[
1
]
==
0
))
error
->
all
(
FLERR
,
"Invalid fix rigid npt/nph command pressure settings"
);
if
(
pcouple
==
YZ
&&
(
p_flag
[
1
]
==
0
||
p_flag
[
2
]
==
0
))
error
->
all
(
FLERR
,
"Invalid fix rigid npt/nph command pressure settings"
);
if
(
pcouple
==
XZ
&&
(
p_flag
[
0
]
==
0
||
p_flag
[
2
]
==
0
))
error
->
all
(
FLERR
,
"Invalid fix rigid npt/nph command pressure settings"
);
// require periodicity in tensile dimension
if
(
p_flag
[
0
]
&&
domain
->
xperiodic
==
0
)
error
->
all
(
FLERR
,
"Cannot use fix rigid npt/nph on a non-periodic dimension"
);
if
(
p_flag
[
1
]
&&
domain
->
yperiodic
==
0
)
error
->
all
(
FLERR
,
"Cannot use fix rigid npt/nph on a non-periodic dimension"
);
if
(
p_flag
[
2
]
&&
domain
->
zperiodic
==
0
)
error
->
all
(
FLERR
,
"Cannot use fix rigid npt/nph on a non-periodic dimension"
);
if
(
pcouple
==
XYZ
&&
dimension
==
3
&&
(
p_start
[
0
]
!=
p_start
[
1
]
||
p_start
[
0
]
!=
p_start
[
2
]
||
p_stop
[
0
]
!=
p_stop
[
1
]
||
p_stop
[
0
]
!=
p_stop
[
2
]
||
p_period
[
0
]
!=
p_period
[
1
]
||
p_period
[
0
]
!=
p_period
[
2
]))
error
->
all
(
FLERR
,
"Invalid fix rigid npt/nph pressure settings"
);
if
(
pcouple
==
XYZ
&&
dimension
==
2
&&
(
p_start
[
0
]
!=
p_start
[
1
]
||
p_stop
[
0
]
!=
p_stop
[
1
]
||
p_period
[
0
]
!=
p_period
[
1
]))
error
->
all
(
FLERR
,
"Invalid fix rigid npt/nph pressure settings"
);
if
(
pcouple
==
XY
&&
(
p_start
[
0
]
!=
p_start
[
1
]
||
p_stop
[
0
]
!=
p_stop
[
1
]
||
p_period
[
0
]
!=
p_period
[
1
]))
error
->
all
(
FLERR
,
"Invalid fix rigid npt/nph pressure settings"
);
if
(
pcouple
==
YZ
&&
(
p_start
[
1
]
!=
p_start
[
2
]
||
p_stop
[
1
]
!=
p_stop
[
2
]
||
p_period
[
1
]
!=
p_period
[
2
]))
error
->
all
(
FLERR
,
"Invalid fix rigid npt/nph pressure settings"
);
if
(
pcouple
==
XZ
&&
(
p_start
[
0
]
!=
p_start
[
2
]
||
p_stop
[
0
]
!=
p_stop
[
2
]
||
p_period
[
0
]
!=
p_period
[
2
]))
error
->
all
(
FLERR
,
"Invalid fix rigid npt/nph pressure settings"
);
if
((
tstat_flag
&&
t_period
<=
0.0
)
||
(
p_flag
[
0
]
&&
p_period
[
0
]
<=
0.0
)
||
(
p_flag
[
1
]
&&
p_period
[
1
]
<=
0.0
)
||
(
p_flag
[
2
]
&&
p_period
[
2
]
<=
0.0
))
error
->
all
(
FLERR
,
"Fix rigid nvt/npt/nph damping parameters must be > 0.0"
);
// memory allocation and initialization
memory
->
create
(
conjqm
,
nbody
,
4
,
"rigid_nh:conjqm"
);
if
(
tstat_flag
||
pstat_flag
)
{
allocate_chain
();
allocate_order
();
}
if
(
tstat_flag
)
{
eta_t
[
0
]
=
eta_r
[
0
]
=
0.0
;
eta_dot_t
[
0
]
=
eta_dot_r
[
0
]
=
0.0
;
f_eta_t
[
0
]
=
f_eta_r
[
0
]
=
0.0
;
for
(
int
i
=
1
;
i
<
t_chain
;
i
++
)
{
eta_t
[
i
]
=
eta_r
[
i
]
=
0.0
;
eta_dot_t
[
i
]
=
eta_dot_r
[
i
]
=
0.0
;
}
}
if
(
pstat_flag
)
{
epsilon_dot
[
0
]
=
epsilon_dot
[
1
]
=
epsilon_dot
[
2
]
=
0.0
;
eta_b
[
0
]
=
eta_dot_b
[
0
]
=
f_eta_b
[
0
]
=
0.0
;
for
(
int
i
=
1
;
i
<
p_chain
;
i
++
)
eta_b
[
i
]
=
eta_dot_b
[
i
]
=
0.0
;
}
// rigid body pointers
nrigidfix
=
0
;
rfix
=
NULL
;
vol0
=
0.0
;
t0
=
1.0
;
tcomputeflag
=
0
;
pcomputeflag
=
0
;
}
/* ---------------------------------------------------------------------- */
FixRigidNH
::~
FixRigidNH
()
{
memory
->
destroy
(
conjqm
);
if
(
tstat_flag
||
pstat_flag
)
{
deallocate_chain
();
deallocate_order
();
}
if
(
rfix
)
delete
[]
rfix
;
if
(
tcomputeflag
)
{
modify
->
delete_compute
(
id_temp
);
delete
[]
id_temp
;
}
// delete pressure if fix created it
if
(
pstat_flag
)
{
if
(
pcomputeflag
)
modify
->
delete_compute
(
id_press
);
delete
[]
id_press
;
}
}
/* ---------------------------------------------------------------------- */
int
FixRigidNH
::
setmask
()
{
int
mask
=
0
;
mask
=
FixRigid
::
setmask
();
if
(
tstat_flag
||
pstat_flag
)
mask
|=
THERMO_ENERGY
;
return
mask
;
}
/* ---------------------------------------------------------------------- */
void
FixRigidNH
::
init
()
{
FixRigid
::
init
();
// recheck that dilate group has not been deleted
if
(
allremap
==
0
)
{
int
idilate
=
group
->
find
(
id_dilate
);
if
(
idilate
==
-
1
)
error
->
all
(
FLERR
,
"Fix rigid npt/nph dilate group ID does not exist"
);
dilate_group_bit
=
group
->
bitmask
[
idilate
];
}
// initialize thermostats
// set timesteps, constants
// store Yoshida-Suzuki integrator parameters
dtv
=
update
->
dt
;
dtf
=
0.5
*
update
->
dt
*
force
->
ftm2v
;
dtq
=
0.5
*
update
->
dt
;
boltz
=
force
->
boltz
;
nktv2p
=
force
->
nktv2p
;
mvv2e
=
force
->
mvv2e
;
if
(
force
->
kspace
)
kspace_flag
=
1
;
else
kspace_flag
=
0
;
nf_t
=
nf_r
=
dimension
*
nbody
;
for
(
int
ibody
=
0
;
ibody
<
nbody
;
ibody
++
)
for
(
int
k
=
0
;
k
<
domain
->
dimension
;
k
++
)
if
(
fabs
(
inertia
[
ibody
][
k
])
<
EPSILON
)
nf_r
--
;
// see Table 1 in Kamberaj et al
if
(
tstat_flag
||
pstat_flag
)
{
if
(
t_order
==
3
)
{
w
[
0
]
=
1.0
/
(
2.0
-
pow
(
2.0
,
1.0
/
3.0
));
w
[
1
]
=
1.0
-
2.0
*
w
[
0
];
w
[
2
]
=
w
[
0
];
}
else
if
(
t_order
==
5
)
{
w
[
0
]
=
1.0
/
(
4.0
-
pow
(
4.0
,
1.0
/
3.0
));
w
[
1
]
=
w
[
0
];
w
[
2
]
=
1.0
-
4.0
*
w
[
0
];
w
[
3
]
=
w
[
0
];
w
[
4
]
=
w
[
0
];
}
}
g_f
=
nf_t
+
nf_r
;
onednft
=
1.0
+
(
double
)(
dimension
)
/
(
double
)
g_f
;
onednfr
=
(
double
)
(
dimension
)
/
(
double
)
g_f
;
int
icompute
;
if
(
tcomputeflag
)
{
icompute
=
modify
->
find_compute
(
id_temp
);
if
(
icompute
<
0
)
error
->
all
(
FLERR
,
"Temp ID for fix rigid npt/nph does not exist"
);
temperature
=
modify
->
compute
[
icompute
];
}
if
(
pstat_flag
)
{
if
(
domain
->
triclinic
)
error
->
all
(
FLERR
,
"fix rigid npt/nph does not yet allow triclinic box"
);
// ensure no conflict with fix deform
for
(
int
i
=
0
;
i
<
modify
->
nfix
;
i
++
)
if
(
strcmp
(
modify
->
fix
[
i
]
->
style
,
"deform"
)
==
0
)
{
int
*
dimflag
=
((
FixDeform
*
)
modify
->
fix
[
i
])
->
dimflag
;
if
((
p_flag
[
0
]
&&
dimflag
[
0
])
||
(
p_flag
[
1
]
&&
dimflag
[
1
])
||
(
p_flag
[
2
]
&&
dimflag
[
2
]))
error
->
all
(
FLERR
,
"Cannot use fix rigid npt/nph and fix deform on "
"same component of stress tensor"
);
}
// set frequency
p_freq_max
=
0.0
;
p_freq_max
=
MAX
(
p_freq
[
0
],
p_freq
[
1
]);
p_freq_max
=
MAX
(
p_freq_max
,
p_freq
[
2
]);
// tally the number of dimensions that are barostatted
// set initial volume and reference cell, if not already done
pdim
=
p_flag
[
0
]
+
p_flag
[
1
]
+
p_flag
[
2
];
if
(
vol0
==
0.0
)
{
if
(
dimension
==
2
)
vol0
=
domain
->
xprd
*
domain
->
yprd
;
else
vol0
=
domain
->
xprd
*
domain
->
yprd
*
domain
->
zprd
;
}
// set pressure compute ptr
icompute
=
modify
->
find_compute
(
id_press
);
if
(
icompute
<
0
)
error
->
all
(
FLERR
,
"Press ID for fix rigid npt/nph does not exist"
);
pressure
=
modify
->
compute
[
icompute
];
// detect if any rigid fixes exist so rigid bodies move on remap
// rfix[] = indices to each fix rigid
// this will include self
if
(
rfix
)
delete
[]
rfix
;
nrigidfix
=
0
;
rfix
=
NULL
;
for
(
int
i
=
0
;
i
<
modify
->
nfix
;
i
++
)
if
(
modify
->
fix
[
i
]
->
rigid_flag
)
nrigidfix
++
;
if
(
nrigidfix
)
{
rfix
=
new
int
[
nrigidfix
];
nrigidfix
=
0
;
for
(
int
i
=
0
;
i
<
modify
->
nfix
;
i
++
)
if
(
modify
->
fix
[
i
]
->
rigid_flag
)
rfix
[
nrigidfix
++
]
=
i
;
}
}
}
/* ---------------------------------------------------------------------- */
void
FixRigidNH
::
setup
(
int
vflag
)
{
FixRigid
::
setup
(
vflag
);
double
mbody
[
3
];
akin_t
=
akin_r
=
0.0
;
for
(
int
ibody
=
0
;
ibody
<
nbody
;
ibody
++
)
{
MathExtra
::
transpose_matvec
(
ex_space
[
ibody
],
ey_space
[
ibody
],
ez_space
[
ibody
],
angmom
[
ibody
],
mbody
);
MathExtra
::
quatvec
(
quat
[
ibody
],
mbody
,
conjqm
[
ibody
]);
conjqm
[
ibody
][
0
]
*=
2.0
;
conjqm
[
ibody
][
1
]
*=
2.0
;
conjqm
[
ibody
][
2
]
*=
2.0
;
conjqm
[
ibody
][
3
]
*=
2.0
;
if
(
tstat_flag
||
pstat_flag
)
{
akin_t
+=
masstotal
[
ibody
]
*
(
vcm
[
ibody
][
0
]
*
vcm
[
ibody
][
0
]
+
vcm
[
ibody
][
1
]
*
vcm
[
ibody
][
1
]
+
vcm
[
ibody
][
2
]
*
vcm
[
ibody
][
2
]);
akin_r
+=
angmom
[
ibody
][
0
]
*
omega
[
ibody
][
0
]
+
angmom
[
ibody
][
1
]
*
omega
[
ibody
][
1
]
+
angmom
[
ibody
][
2
]
*
omega
[
ibody
][
2
];
}
}
// compute target temperature
if
(
tstat_flag
)
compute_temp_target
();
else
if
(
pstat_flag
)
{
t0
=
temperature
->
compute_scalar
();
if
(
t0
==
0.0
)
{
if
(
strcmp
(
update
->
unit_style
,
"lj"
)
==
0
)
t0
=
1.0
;
else
t0
=
300.0
;
}
t_target
=
t0
;
}
// compute target pressure
// compute current pressure
// trigger virial computation on next timestep
if
(
pstat_flag
)
{
compute_press_target
();
temperature
->
compute_scalar
();
if
(
pstyle
==
ISO
)
pressure
->
compute_scalar
();
else
pressure
->
compute_vector
();
couple
();
pressure
->
addstep
(
update
->
ntimestep
+
1
);
}
// initialize thermostat/barostat settings
double
kt
,
t_mass
,
tb_mass
;
kt
=
boltz
*
t_target
;
if
(
tstat_flag
)
{
t_mass
=
kt
/
(
t_freq
*
t_freq
);
q_t
[
0
]
=
nf_t
*
t_mass
;
q_r
[
0
]
=
nf_r
*
t_mass
;
for
(
int
i
=
1
;
i
<
t_chain
;
i
++
)
q_t
[
i
]
=
q_r
[
i
]
=
t_mass
;
for
(
int
i
=
1
;
i
<
t_chain
;
i
++
)
{
f_eta_t
[
i
]
=
(
q_t
[
i
-
1
]
*
eta_dot_t
[
i
-
1
]
*
eta_dot_t
[
i
-
1
]
-
kt
)
/
q_t
[
i
];
f_eta_r
[
i
]
=
(
q_r
[
i
-
1
]
*
eta_dot_r
[
i
-
1
]
*
eta_dot_r
[
i
-
1
]
-
kt
)
/
q_r
[
i
];
}
}
// initial forces on barostat thermostat variables
if
(
pstat_flag
)
{
for
(
int
i
=
0
;
i
<
3
;
i
++
)
if
(
p_flag
[
i
])
{
epsilon_mass
[
i
]
=
(
g_f
+
dimension
)
*
kt
/
(
p_freq
[
i
]
*
p_freq
[
i
]);
epsilon
[
i
]
=
log
(
vol0
)
/
dimension
;
}
tb_mass
=
kt
/
(
p_freq_max
*
p_freq_max
);
q_b
[
0
]
=
dimension
*
dimension
*
tb_mass
;
for
(
int
i
=
1
;
i
<
p_chain
;
i
++
)
{
q_b
[
i
]
=
tb_mass
;
f_eta_b
[
i
]
=
(
q_b
[
i
]
*
eta_dot_b
[
i
-
1
]
*
eta_dot_b
[
i
-
1
]
-
kt
)
/
q_b
[
i
];
}
}
// update order/timestep dependent coefficients
if
(
tstat_flag
||
pstat_flag
)
{
for
(
int
i
=
0
;
i
<
t_order
;
i
++
)
{
wdti1
[
i
]
=
w
[
i
]
*
dtv
/
t_iter
;
wdti2
[
i
]
=
wdti1
[
i
]
/
2.0
;
wdti4
[
i
]
=
wdti1
[
i
]
/
4.0
;
}
}
}
/* ----------------------------------------------------------------------
perform preforce velocity Verlet integration
see Kamberaj paper for step references
------------------------------------------------------------------------- */
void
FixRigidNH
::
initial_integrate
(
int
vflag
)
{
double
tmp
,
scale_r
,
scale_t
[
3
],
scale_v
[
3
];
double
dtfm
,
mbody
[
3
],
tbody
[
3
],
fquat
[
4
];
double
dtf2
=
dtf
*
2.0
;
// compute target temperature
// update thermostat chains coupled to particles
if
(
tstat_flag
)
{
compute_temp_target
();
nhc_temp_integrate
();
}
// compute target pressure
// update epsilon dot
// update thermostat coupled to barostat
if
(
pstat_flag
)
{
nhc_press_integrate
();
if
(
pstyle
==
ISO
)
{
temperature
->
compute_scalar
();
pressure
->
compute_scalar
();
}
else
{
temperature
->
compute_vector
();
pressure
->
compute_vector
();
}
couple
();
pressure
->
addstep
(
update
->
ntimestep
+
1
);
compute_press_target
();
nh_epsilon_dot
();
}
// compute scale variables
scale_t
[
0
]
=
scale_t
[
1
]
=
scale_t
[
2
]
=
1.0
;
scale_v
[
0
]
=
scale_v
[
1
]
=
scale_v
[
2
]
=
1.0
;
scale_r
=
1.0
;
if
(
tstat_flag
)
{
akin_t
=
akin_r
=
0.0
;
tmp
=
exp
(
-
dtq
*
eta_dot_t
[
0
]);
scale_t
[
0
]
=
scale_t
[
1
]
=
scale_t
[
2
]
=
tmp
;
tmp
=
exp
(
-
dtq
*
eta_dot_r
[
0
]);
scale_r
=
tmp
;
}
if
(
pstat_flag
)
{
akin_t
=
akin_r
=
0.0
;
scale_t
[
0
]
*=
exp
(
-
dtq
*
(
epsilon_dot
[
0
]
+
mtk_term2
));
scale_t
[
1
]
*=
exp
(
-
dtq
*
(
epsilon_dot
[
1
]
+
mtk_term2
));
scale_t
[
2
]
*=
exp
(
-
dtq
*
(
epsilon_dot
[
2
]
+
mtk_term2
));
scale_r
*=
exp
(
-
dtq
*
(
pdim
*
mtk_term2
));
tmp
=
dtq
*
epsilon_dot
[
0
];
scale_v
[
0
]
=
dtv
*
exp
(
tmp
)
*
maclaurin_series
(
tmp
);
tmp
=
dtq
*
epsilon_dot
[
1
];
scale_v
[
1
]
=
dtv
*
exp
(
tmp
)
*
maclaurin_series
(
tmp
);
tmp
=
dtq
*
epsilon_dot
[
2
];
scale_v
[
2
]
=
dtv
*
exp
(
tmp
)
*
maclaurin_series
(
tmp
);
}
// update xcm, vcm, quat, conjqm and angmom
for
(
int
ibody
=
0
;
ibody
<
nbody
;
ibody
++
)
{
// step 1.1 - update vcm by 1/2 step
dtfm
=
dtf
/
masstotal
[
ibody
];
vcm
[
ibody
][
0
]
+=
dtfm
*
fcm
[
ibody
][
0
]
*
fflag
[
ibody
][
0
];
vcm
[
ibody
][
1
]
+=
dtfm
*
fcm
[
ibody
][
1
]
*
fflag
[
ibody
][
1
];
vcm
[
ibody
][
2
]
+=
dtfm
*
fcm
[
ibody
][
2
]
*
fflag
[
ibody
][
2
];
if
(
tstat_flag
||
pstat_flag
)
{
vcm
[
ibody
][
0
]
*=
scale_t
[
0
];
vcm
[
ibody
][
1
]
*=
scale_t
[
1
];
vcm
[
ibody
][
2
]
*=
scale_t
[
2
];
tmp
=
vcm
[
ibody
][
0
]
*
vcm
[
ibody
][
0
]
+
vcm
[
ibody
][
1
]
*
vcm
[
ibody
][
1
]
+
vcm
[
ibody
][
2
]
*
vcm
[
ibody
][
2
];
akin_t
+=
masstotal
[
ibody
]
*
tmp
;
}
// step 1.2 - update xcm by full step
if
(
!
pstat_flag
)
{
xcm
[
ibody
][
0
]
+=
dtv
*
vcm
[
ibody
][
0
];
xcm
[
ibody
][
1
]
+=
dtv
*
vcm
[
ibody
][
1
];
xcm
[
ibody
][
2
]
+=
dtv
*
vcm
[
ibody
][
2
];
}
else
{
xcm
[
ibody
][
0
]
+=
scale_v
[
0
]
*
vcm
[
ibody
][
0
];
xcm
[
ibody
][
1
]
+=
scale_v
[
1
]
*
vcm
[
ibody
][
1
];
xcm
[
ibody
][
2
]
+=
scale_v
[
2
]
*
vcm
[
ibody
][
2
];
}
// step 1.3 - apply torque (body coords) to quaternion momentum
torque
[
ibody
][
0
]
*=
tflag
[
ibody
][
0
];
torque
[
ibody
][
1
]
*=
tflag
[
ibody
][
1
];
torque
[
ibody
][
2
]
*=
tflag
[
ibody
][
2
];
MathExtra
::
transpose_matvec
(
ex_space
[
ibody
],
ey_space
[
ibody
],
ez_space
[
ibody
],
torque
[
ibody
],
tbody
);
MathExtra
::
quatvec
(
quat
[
ibody
],
tbody
,
fquat
);
conjqm
[
ibody
][
0
]
+=
dtf2
*
fquat
[
0
];
conjqm
[
ibody
][
1
]
+=
dtf2
*
fquat
[
1
];
conjqm
[
ibody
][
2
]
+=
dtf2
*
fquat
[
2
];
conjqm
[
ibody
][
3
]
+=
dtf2
*
fquat
[
3
];
if
(
tstat_flag
||
pstat_flag
)
{
conjqm
[
ibody
][
0
]
*=
scale_r
;
conjqm
[
ibody
][
1
]
*=
scale_r
;
conjqm
[
ibody
][
2
]
*=
scale_r
;
conjqm
[
ibody
][
3
]
*=
scale_r
;
}
// step 1.4 to 1.13 - use no_squish rotate to update p and q
no_squish_rotate
(
3
,
conjqm
[
ibody
],
quat
[
ibody
],
inertia
[
ibody
],
dtq
);
no_squish_rotate
(
2
,
conjqm
[
ibody
],
quat
[
ibody
],
inertia
[
ibody
],
dtq
);
no_squish_rotate
(
1
,
conjqm
[
ibody
],
quat
[
ibody
],
inertia
[
ibody
],
dtv
);
no_squish_rotate
(
2
,
conjqm
[
ibody
],
quat
[
ibody
],
inertia
[
ibody
],
dtq
);
no_squish_rotate
(
3
,
conjqm
[
ibody
],
quat
[
ibody
],
inertia
[
ibody
],
dtq
);
// update exyz_space
// transform p back to angmom
// update angular velocity
MathExtra
::
q_to_exyz
(
quat
[
ibody
],
ex_space
[
ibody
],
ey_space
[
ibody
],
ez_space
[
ibody
]);
MathExtra
::
invquatvec
(
quat
[
ibody
],
conjqm
[
ibody
],
mbody
);
MathExtra
::
matvec
(
ex_space
[
ibody
],
ey_space
[
ibody
],
ez_space
[
ibody
],
mbody
,
angmom
[
ibody
]);
angmom
[
ibody
][
0
]
*=
0.5
;
angmom
[
ibody
][
1
]
*=
0.5
;
angmom
[
ibody
][
2
]
*=
0.5
;
MathExtra
::
angmom_to_omega
(
angmom
[
ibody
],
ex_space
[
ibody
],
ey_space
[
ibody
],
ez_space
[
ibody
],
inertia
[
ibody
],
omega
[
ibody
]);
if
(
tstat_flag
||
pstat_flag
)
{
akin_r
+=
angmom
[
ibody
][
0
]
*
omega
[
ibody
][
0
]
+
angmom
[
ibody
][
1
]
*
omega
[
ibody
][
1
]
+
angmom
[
ibody
][
2
]
*
omega
[
ibody
][
2
];
}
}
// virial setup before call to set_xv
if
(
vflag
)
v_setup
(
vflag
);
else
evflag
=
0
;
// remap simulation box by 1/2 step
if
(
pstat_flag
)
remap
();
// set coords/orient and velocity/rotation of atoms in rigid bodies
// from quarternion and omega
set_xv
();
// remap simulation box by full step
// redo KSpace coeffs since volume has changed
if
(
pstat_flag
)
{
remap
();
if
(
kspace_flag
)
force
->
kspace
->
setup
();
}
}
/* ---------------------------------------------------------------------- */
void
FixRigidNH
::
final_integrate
()
{
int
i
,
ibody
;
double
tmp
,
scale_t
[
3
],
scale_r
;
double
dtfm
,
xy
,
xz
,
yz
;
double
mbody
[
3
],
tbody
[
3
],
fquat
[
4
];
double
dtf2
=
dtf
*
2.0
;
// compute scale variables
scale_t
[
0
]
=
scale_t
[
1
]
=
scale_t
[
2
]
=
1.0
;
scale_r
=
1.0
;
if
(
tstat_flag
)
{
tmp
=
exp
(
-
1.0
*
dtq
*
eta_dot_t
[
0
]);
scale_t
[
0
]
=
scale_t
[
1
]
=
scale_t
[
2
]
=
tmp
;
scale_r
=
exp
(
-
1.0
*
dtq
*
eta_dot_r
[
0
]);
}
if
(
pstat_flag
)
{
scale_t
[
0
]
*=
exp
(
-
dtq
*
(
epsilon_dot
[
0
]
+
mtk_term2
));
scale_t
[
1
]
*=
exp
(
-
dtq
*
(
epsilon_dot
[
1
]
+
mtk_term2
));
scale_t
[
2
]
*=
exp
(
-
dtq
*
(
epsilon_dot
[
2
]
+
mtk_term2
));
scale_r
*=
exp
(
-
dtq
*
(
pdim
*
mtk_term2
));
akin_t
=
akin_r
=
0.0
;
}
// sum over atoms to get force and torque on rigid body
int
*
image
=
atom
->
image
;
double
**
x
=
atom
->
x
;
double
**
f
=
atom
->
f
;
int
nlocal
=
atom
->
nlocal
;
double
xprd
=
domain
->
xprd
;
double
yprd
=
domain
->
yprd
;
double
zprd
=
domain
->
zprd
;
if
(
triclinic
)
{
xy
=
domain
->
xy
;
xz
=
domain
->
xz
;
yz
=
domain
->
yz
;
}
int
xbox
,
ybox
,
zbox
;
double
xunwrap
,
yunwrap
,
zunwrap
,
dx
,
dy
,
dz
;
for
(
ibody
=
0
;
ibody
<
nbody
;
ibody
++
)
for
(
i
=
0
;
i
<
6
;
i
++
)
sum
[
ibody
][
i
]
=
0.0
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
body
[
i
]
<
0
)
continue
;
ibody
=
body
[
i
];
sum
[
ibody
][
0
]
+=
f
[
i
][
0
];
sum
[
ibody
][
1
]
+=
f
[
i
][
1
];
sum
[
ibody
][
2
]
+=
f
[
i
][
2
];
xbox
=
(
image
[
i
]
&
IMGMASK
)
-
IMGMAX
;
ybox
=
(
image
[
i
]
>>
IMGBITS
&
IMGMASK
)
-
IMGMAX
;
zbox
=
(
image
[
i
]
>>
IMG2BITS
)
-
IMGMAX
;
if
(
triclinic
==
0
)
{
xunwrap
=
x
[
i
][
0
]
+
xbox
*
xprd
;
yunwrap
=
x
[
i
][
1
]
+
ybox
*
yprd
;
zunwrap
=
x
[
i
][
2
]
+
zbox
*
zprd
;
}
else
{
xunwrap
=
x
[
i
][
0
]
+
xbox
*
xprd
+
ybox
*
xy
+
zbox
*
xz
;
yunwrap
=
x
[
i
][
1
]
+
ybox
*
yprd
+
zbox
*
yz
;
zunwrap
=
x
[
i
][
2
]
+
zbox
*
zprd
;
}
dx
=
xunwrap
-
xcm
[
ibody
][
0
];
dy
=
yunwrap
-
xcm
[
ibody
][
1
];
dz
=
zunwrap
-
xcm
[
ibody
][
2
];
sum
[
ibody
][
3
]
+=
dy
*
f
[
i
][
2
]
-
dz
*
f
[
i
][
1
];
sum
[
ibody
][
4
]
+=
dz
*
f
[
i
][
0
]
-
dx
*
f
[
i
][
2
];
sum
[
ibody
][
5
]
+=
dx
*
f
[
i
][
1
]
-
dy
*
f
[
i
][
0
];
}
// extended particles add their torque to torque of body
if
(
extended
)
{
double
**
torque_one
=
atom
->
torque
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
body
[
i
]
<
0
)
continue
;
ibody
=
body
[
i
];
if
(
eflags
[
i
]
&
TORQUE
)
{
sum
[
ibody
][
3
]
+=
torque_one
[
i
][
0
];
sum
[
ibody
][
4
]
+=
torque_one
[
i
][
1
];
sum
[
ibody
][
5
]
+=
torque_one
[
i
][
2
];
}
}
}
MPI_Allreduce
(
sum
[
0
],
all
[
0
],
6
*
nbody
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
// update vcm and angmom
// include Langevin thermostat forces
// fflag,tflag = 0 for some dimensions in 2d
for
(
ibody
=
0
;
ibody
<
nbody
;
ibody
++
)
{
fcm
[
ibody
][
0
]
=
all
[
ibody
][
0
]
+
langextra
[
ibody
][
0
];
fcm
[
ibody
][
1
]
=
all
[
ibody
][
1
]
+
langextra
[
ibody
][
1
];
fcm
[
ibody
][
2
]
=
all
[
ibody
][
2
]
+
langextra
[
ibody
][
2
];
torque
[
ibody
][
0
]
=
all
[
ibody
][
3
]
+
langextra
[
ibody
][
3
];
torque
[
ibody
][
1
]
=
all
[
ibody
][
4
]
+
langextra
[
ibody
][
4
];
torque
[
ibody
][
2
]
=
all
[
ibody
][
5
]
+
langextra
[
ibody
][
5
];
// update vcm by 1/2 step
dtfm
=
dtf
/
masstotal
[
ibody
];
if
(
tstat_flag
||
pstat_flag
)
{
vcm
[
ibody
][
0
]
*=
scale_t
[
0
];
vcm
[
ibody
][
1
]
*=
scale_t
[
1
];
vcm
[
ibody
][
2
]
*=
scale_t
[
2
];
}
vcm
[
ibody
][
0
]
+=
dtfm
*
fcm
[
ibody
][
0
]
*
fflag
[
ibody
][
0
];
vcm
[
ibody
][
1
]
+=
dtfm
*
fcm
[
ibody
][
1
]
*
fflag
[
ibody
][
1
];
vcm
[
ibody
][
2
]
+=
dtfm
*
fcm
[
ibody
][
2
]
*
fflag
[
ibody
][
2
];
if
(
pstat_flag
)
{
tmp
=
vcm
[
ibody
][
0
]
*
vcm
[
ibody
][
0
]
+
vcm
[
ibody
][
1
]
*
vcm
[
ibody
][
1
]
+
vcm
[
ibody
][
2
]
*
vcm
[
ibody
][
2
];
akin_t
+=
masstotal
[
ibody
]
*
tmp
;
}
// update conjqm, then transform to angmom, set velocity again
// virial is already setup from initial_integrate
torque
[
ibody
][
0
]
*=
tflag
[
ibody
][
0
];
torque
[
ibody
][
1
]
*=
tflag
[
ibody
][
1
];
torque
[
ibody
][
2
]
*=
tflag
[
ibody
][
2
];
MathExtra
::
transpose_matvec
(
ex_space
[
ibody
],
ey_space
[
ibody
],
ez_space
[
ibody
],
torque
[
ibody
],
tbody
);
MathExtra
::
quatvec
(
quat
[
ibody
],
tbody
,
fquat
);
if
(
tstat_flag
||
pstat_flag
)
{
conjqm
[
ibody
][
0
]
=
scale_r
*
conjqm
[
ibody
][
0
]
+
dtf2
*
fquat
[
0
];
conjqm
[
ibody
][
1
]
=
scale_r
*
conjqm
[
ibody
][
1
]
+
dtf2
*
fquat
[
1
];
conjqm
[
ibody
][
2
]
=
scale_r
*
conjqm
[
ibody
][
2
]
+
dtf2
*
fquat
[
2
];
conjqm
[
ibody
][
3
]
=
scale_r
*
conjqm
[
ibody
][
3
]
+
dtf2
*
fquat
[
3
];
}
else
{
conjqm
[
ibody
][
0
]
+=
dtf2
*
fquat
[
0
];
conjqm
[
ibody
][
1
]
+=
dtf2
*
fquat
[
1
];
conjqm
[
ibody
][
2
]
+=
dtf2
*
fquat
[
2
];
conjqm
[
ibody
][
3
]
+=
dtf2
*
fquat
[
3
];
}
MathExtra
::
invquatvec
(
quat
[
ibody
],
conjqm
[
ibody
],
mbody
);
MathExtra
::
matvec
(
ex_space
[
ibody
],
ey_space
[
ibody
],
ez_space
[
ibody
],
mbody
,
angmom
[
ibody
]);
angmom
[
ibody
][
0
]
*=
0.5
;
angmom
[
ibody
][
1
]
*=
0.5
;
angmom
[
ibody
][
2
]
*=
0.5
;
MathExtra
::
angmom_to_omega
(
angmom
[
ibody
],
ex_space
[
ibody
],
ey_space
[
ibody
],
ez_space
[
ibody
],
inertia
[
ibody
],
omega
[
ibody
]);
if
(
pstat_flag
)
{
akin_r
+=
angmom
[
ibody
][
0
]
*
omega
[
ibody
][
0
]
+
angmom
[
ibody
][
1
]
*
omega
[
ibody
][
1
]
+
angmom
[
ibody
][
2
]
*
omega
[
ibody
][
2
];
}
}
// set velocity/rotation of atoms in rigid bodies
// virial is already setup from initial_integrate
set_v
();
// compute temperature and pressure tensor
// couple to compute current pressure components
// trigger virial computation on next timestep
if
(
tcomputeflag
)
t_current
=
temperature
->
compute_scalar
();
if
(
pstat_flag
)
{
if
(
pstyle
==
ISO
)
pressure
->
compute_scalar
();
else
pressure
->
compute_vector
();
couple
();
pressure
->
addstep
(
update
->
ntimestep
+
1
);
}
if
(
pstat_flag
)
nh_epsilon_dot
();
// update eta_dot_t and eta_dot_r
// update eta_dot_b
if
(
tstat_flag
)
nhc_temp_integrate
();
if
(
pstat_flag
)
nhc_press_integrate
();
}
/* ---------------------------------------------------------------------- */
void
FixRigidNH
::
nhc_temp_integrate
()
{
int
i
,
j
,
k
;
double
kt
,
gfkt_t
,
gfkt_r
,
tmp
,
ms
,
s
,
s2
;
kt
=
boltz
*
t_target
;
gfkt_t
=
nf_t
*
kt
;
gfkt_r
=
nf_r
*
kt
;
// update thermostat masses
double
t_mass
=
boltz
*
t_target
/
(
t_freq
*
t_freq
);
q_t
[
0
]
=
nf_t
*
t_mass
;
q_r
[
0
]
=
nf_r
*
t_mass
;
for
(
i
=
1
;
i
<
t_chain
;
i
++
)
q_t
[
i
]
=
q_r
[
i
]
=
t_mass
;
// update force of thermostats coupled to particles
f_eta_t
[
0
]
=
(
akin_t
*
mvv2e
-
gfkt_t
)
/
q_t
[
0
];
f_eta_r
[
0
]
=
(
akin_r
*
mvv2e
-
gfkt_r
)
/
q_r
[
0
];
// multiple timestep iteration
for
(
i
=
0
;
i
<
t_iter
;
i
++
)
{
for
(
j
=
0
;
j
<
t_order
;
j
++
)
{
// update thermostat velocities half step
eta_dot_t
[
t_chain
-
1
]
+=
wdti2
[
j
]
*
f_eta_t
[
t_chain
-
1
];
eta_dot_r
[
t_chain
-
1
]
+=
wdti2
[
j
]
*
f_eta_r
[
t_chain
-
1
];
for
(
k
=
1
;
k
<
t_chain
;
k
++
)
{
tmp
=
wdti4
[
j
]
*
eta_dot_t
[
t_chain
-
k
];
ms
=
maclaurin_series
(
tmp
);
s
=
exp
(
-
1.0
*
tmp
);
s2
=
s
*
s
;
eta_dot_t
[
t_chain
-
k
-
1
]
=
eta_dot_t
[
t_chain
-
k
-
1
]
*
s2
+
wdti2
[
j
]
*
f_eta_t
[
t_chain
-
k
-
1
]
*
s
*
ms
;
tmp
=
wdti4
[
j
]
*
eta_dot_r
[
t_chain
-
k
];
ms
=
maclaurin_series
(
tmp
);
s
=
exp
(
-
1.0
*
tmp
);
s2
=
s
*
s
;
eta_dot_r
[
t_chain
-
k
-
1
]
=
eta_dot_r
[
t_chain
-
k
-
1
]
*
s2
+
wdti2
[
j
]
*
f_eta_r
[
t_chain
-
k
-
1
]
*
s
*
ms
;
}
// update thermostat positions a full step
for
(
k
=
0
;
k
<
t_chain
;
k
++
)
{
eta_t
[
k
]
+=
wdti1
[
j
]
*
eta_dot_t
[
k
];
eta_r
[
k
]
+=
wdti1
[
j
]
*
eta_dot_r
[
k
];
}
// update thermostat forces
for
(
k
=
1
;
k
<
t_chain
;
k
++
)
{
f_eta_t
[
k
]
=
q_t
[
k
-
1
]
*
eta_dot_t
[
k
-
1
]
*
eta_dot_t
[
k
-
1
]
-
kt
;
f_eta_t
[
k
]
/=
q_t
[
k
];
f_eta_r
[
k
]
=
q_r
[
k
-
1
]
*
eta_dot_r
[
k
-
1
]
*
eta_dot_r
[
k
-
1
]
-
kt
;
f_eta_r
[
k
]
/=
q_r
[
k
];
}
// update thermostat velocities a full step
for
(
k
=
0
;
k
<
t_chain
-
1
;
k
++
)
{
tmp
=
wdti4
[
j
]
*
eta_dot_t
[
k
+
1
];
ms
=
maclaurin_series
(
tmp
);
s
=
exp
(
-
1.0
*
tmp
);
s2
=
s
*
s
;
eta_dot_t
[
k
]
=
eta_dot_t
[
k
]
*
s2
+
wdti2
[
j
]
*
f_eta_t
[
k
]
*
s
*
ms
;
tmp
=
q_t
[
k
]
*
eta_dot_t
[
k
]
*
eta_dot_t
[
k
]
-
kt
;
f_eta_t
[
k
+
1
]
=
tmp
/
q_t
[
k
+
1
];
tmp
=
wdti4
[
j
]
*
eta_dot_r
[
k
+
1
];
ms
=
maclaurin_series
(
tmp
);
s
=
exp
(
-
1.0
*
tmp
);
s2
=
s
*
s
;
eta_dot_r
[
k
]
=
eta_dot_r
[
k
]
*
s2
+
wdti2
[
j
]
*
f_eta_r
[
k
]
*
s
*
ms
;
tmp
=
q_r
[
k
]
*
eta_dot_r
[
k
]
*
eta_dot_r
[
k
]
-
kt
;
f_eta_r
[
k
+
1
]
=
tmp
/
q_r
[
k
+
1
];
}
eta_dot_t
[
t_chain
-
1
]
+=
wdti2
[
j
]
*
f_eta_t
[
t_chain
-
1
];
eta_dot_r
[
t_chain
-
1
]
+=
wdti2
[
j
]
*
f_eta_r
[
t_chain
-
1
];
}
}
}
/* ---------------------------------------------------------------------- */
void
FixRigidNH
::
nhc_press_integrate
()
{
int
i
,
k
;
double
tmp
,
s
,
s2
,
ms
,
kecurrent
;
double
kt
=
boltz
*
t_target
;
double
lkt_press
=
kt
;
// update thermostat masses
double
tb_mass
=
kt
/
(
p_freq_max
*
p_freq_max
);
q_b
[
0
]
=
tb_mass
;
for
(
int
i
=
1
;
i
<
p_chain
;
i
++
)
{
q_b
[
i
]
=
tb_mass
;
f_eta_b
[
i
]
=
q_b
[
i
-
1
]
*
eta_dot_b
[
i
-
1
]
*
eta_dot_b
[
i
-
1
]
-
kt
;
f_eta_b
[
i
]
/=
q_b
[
i
];
}
// update forces acting on thermostat
kecurrent
=
0.0
;
for
(
i
=
0
;
i
<
3
;
i
++
)
if
(
p_flag
[
i
])
{
epsilon_mass
[
i
]
=
(
g_f
+
dimension
)
*
kt
/
(
p_freq
[
i
]
*
p_freq
[
i
]);
kecurrent
+=
epsilon_mass
[
i
]
*
epsilon_dot
[
i
]
*
epsilon_dot
[
i
];
}
f_eta_b
[
0
]
=
(
kecurrent
-
lkt_press
)
/
q_b
[
0
];
// update thermostat velocities a half step
eta_dot_b
[
p_chain
-
1
]
+=
0.5
*
dtq
*
f_eta_b
[
p_chain
-
1
];
for
(
k
=
0
;
k
<
p_chain
-
1
;
k
++
)
{
tmp
=
0.5
*
dtq
*
eta_dot_b
[
p_chain
-
k
-
1
];
ms
=
maclaurin_series
(
tmp
);
s
=
exp
(
-
0.5
*
tmp
);
s2
=
s
*
s
;
eta_dot_b
[
p_chain
-
k
-
2
]
=
eta_dot_b
[
p_chain
-
k
-
2
]
*
s2
+
dtq
*
f_eta_b
[
p_chain
-
k
-
2
]
*
s
*
ms
;
}
// update thermostat positions
for
(
k
=
0
;
k
<
p_chain
;
k
++
)
eta_b
[
k
]
+=
dtv
*
eta_dot_b
[
k
];
// update epsilon dot
s
=
exp
(
-
1.0
*
dtq
*
eta_dot_b
[
0
]);
for
(
i
=
0
;
i
<
3
;
i
++
)
if
(
p_flag
[
i
])
epsilon_dot
[
i
]
*=
s
;
kecurrent
=
0.0
;
for
(
i
=
0
;
i
<
3
;
i
++
)
if
(
p_flag
[
i
])
kecurrent
+=
epsilon_mass
[
i
]
*
epsilon_dot
[
i
]
*
epsilon_dot
[
i
];
f_eta_b
[
0
]
=
(
kecurrent
-
lkt_press
)
/
q_b
[
0
];
// update thermostat velocites a full step
for
(
k
=
0
;
k
<
p_chain
-
1
;
k
++
)
{
tmp
=
0.5
*
dtq
*
eta_dot_b
[
k
+
1
];
ms
=
maclaurin_series
(
tmp
);
s
=
exp
(
-
0.5
*
tmp
);
s2
=
s
*
s
;
eta_dot_b
[
k
]
=
eta_dot_b
[
k
]
*
s2
+
dtq
*
f_eta_b
[
k
]
*
s
*
ms
;
tmp
=
q_b
[
k
]
*
eta_dot_b
[
k
]
*
eta_dot_b
[
k
]
-
kt
;
f_eta_b
[
k
+
1
]
=
tmp
/
q_b
[
k
+
1
];
}
eta_dot_b
[
p_chain
-
1
]
+=
0.5
*
dtq
*
f_eta_b
[
p_chain
-
1
];
}
/* ----------------------------------------------------------------------
compute kinetic energy in the extended Hamiltonian
conserved quantity = sum of returned energy and potential energy
-----------------------------------------------------------------------*/
double
FixRigidNH
::
compute_scalar
()
{
int
i
,
k
,
ibody
;
double
kt
=
boltz
*
t_target
;
double
energy
,
ke_t
,
ke_q
,
tmp
,
Pkq
[
4
];
// compute the kinetic parts of H_NVE in Kameraj et al (JCP 2005, pp 224114)
// translational kinetic energy
ke_t
=
0.0
;
for
(
ibody
=
0
;
ibody
<
nbody
;
ibody
++
)
ke_t
+=
0.5
*
masstotal
[
ibody
]
*
(
vcm
[
ibody
][
0
]
*
vcm
[
ibody
][
0
]
+
vcm
[
ibody
][
1
]
*
vcm
[
ibody
][
1
]
+
vcm
[
ibody
][
2
]
*
vcm
[
ibody
][
2
]);
// rotational kinetic energy
ke_q
=
0.0
;
for
(
ibody
=
0
;
ibody
<
nbody
;
ibody
++
)
{
for
(
k
=
1
;
k
<
4
;
k
++
)
{
if
(
k
==
1
)
{
Pkq
[
0
]
=
-
quat
[
ibody
][
1
];
Pkq
[
1
]
=
quat
[
ibody
][
0
];
Pkq
[
2
]
=
quat
[
ibody
][
3
];
Pkq
[
3
]
=
-
quat
[
ibody
][
2
];
}
else
if
(
k
==
2
)
{
Pkq
[
0
]
=
-
quat
[
ibody
][
2
];
Pkq
[
1
]
=
-
quat
[
ibody
][
3
];
Pkq
[
2
]
=
quat
[
ibody
][
0
];
Pkq
[
3
]
=
quat
[
ibody
][
1
];
}
else
if
(
k
==
3
)
{
Pkq
[
0
]
=
-
quat
[
ibody
][
3
];
Pkq
[
1
]
=
quat
[
ibody
][
2
];
Pkq
[
2
]
=
-
quat
[
ibody
][
1
];
Pkq
[
3
]
=
quat
[
ibody
][
0
];
}
tmp
=
conjqm
[
ibody
][
0
]
*
Pkq
[
0
]
+
conjqm
[
ibody
][
1
]
*
Pkq
[
1
]
+
conjqm
[
ibody
][
2
]
*
Pkq
[
2
]
+
conjqm
[
ibody
][
3
]
*
Pkq
[
3
];
tmp
*=
tmp
;
if
(
fabs
(
inertia
[
ibody
][
k
-
1
])
<
1e-6
)
tmp
=
0.0
;
else
tmp
/=
(
8.0
*
inertia
[
ibody
][
k
-
1
]);
ke_q
+=
tmp
;
}
}
energy
=
(
ke_t
+
ke_q
)
*
mvv2e
;
if
(
tstat_flag
)
{
// thermostat chain energy: from equation 12 in Kameraj et al (JCP 2005)
energy
+=
kt
*
(
nf_t
*
eta_t
[
0
]
+
nf_r
*
eta_r
[
0
]);
for
(
i
=
1
;
i
<
t_chain
;
i
++
)
energy
+=
kt
*
(
eta_t
[
i
]
+
eta_r
[
i
]);
for
(
i
=
0
;
i
<
t_chain
;
i
++
)
{
energy
+=
0.5
*
q_t
[
i
]
*
(
eta_dot_t
[
i
]
*
eta_dot_t
[
i
]);
energy
+=
0.5
*
q_r
[
i
]
*
(
eta_dot_r
[
i
]
*
eta_dot_r
[
i
]);
}
}
if
(
pstat_flag
)
{
// using equation 22 in Kameraj et al for H_NPT
for
(
i
=
0
;
i
<
3
;
i
++
)
energy
+=
0.5
*
epsilon_mass
[
i
]
*
epsilon_dot
[
i
]
*
epsilon_dot
[
i
];
double
vol
;
if
(
dimension
==
2
)
vol
=
domain
->
xprd
*
domain
->
yprd
;
else
vol
=
domain
->
xprd
*
domain
->
yprd
*
domain
->
zprd
;
double
p0
=
(
p_target
[
0
]
+
p_target
[
1
]
+
p_target
[
2
])
/
3.0
;
energy
+=
p0
*
vol
/
nktv2p
;
for
(
i
=
0
;
i
<
p_chain
;
i
++
)
{
energy
+=
kt
*
eta_b
[
i
];
energy
+=
0.5
*
q_b
[
i
]
*
(
eta_dot_b
[
i
]
*
eta_dot_b
[
i
]);
}
}
return
energy
;
}
/* ---------------------------------------------------------------------- */
void
FixRigidNH
::
couple
()
{
double
*
tensor
=
pressure
->
vector
;
if
(
pstyle
==
ISO
)
{
p_current
[
0
]
=
p_current
[
1
]
=
p_current
[
2
]
=
pressure
->
scalar
;
}
else
if
(
pcouple
==
XYZ
)
{
double
ave
=
1.0
/
3.0
*
(
tensor
[
0
]
+
tensor
[
1
]
+
tensor
[
2
]);
p_current
[
0
]
=
p_current
[
1
]
=
p_current
[
2
]
=
ave
;
}
else
if
(
pcouple
==
XY
)
{
double
ave
=
0.5
*
(
tensor
[
0
]
+
tensor
[
1
]);
p_current
[
0
]
=
p_current
[
1
]
=
ave
;
p_current
[
2
]
=
tensor
[
2
];
}
else
if
(
pcouple
==
YZ
)
{
double
ave
=
0.5
*
(
tensor
[
1
]
+
tensor
[
2
]);
p_current
[
1
]
=
p_current
[
2
]
=
ave
;
p_current
[
0
]
=
tensor
[
0
];
}
else
if
(
pcouple
==
XZ
)
{
double
ave
=
0.5
*
(
tensor
[
0
]
+
tensor
[
2
]);
p_current
[
0
]
=
p_current
[
2
]
=
ave
;
p_current
[
1
]
=
tensor
[
1
];
}
else
{
p_current
[
0
]
=
tensor
[
0
];
p_current
[
1
]
=
tensor
[
1
];
p_current
[
2
]
=
tensor
[
2
];
}
}
/* ---------------------------------------------------------------------- */
void
FixRigidNH
::
remap
()
{
int
i
;
double
oldlo
,
oldhi
,
ctr
,
expfac
;
double
**
x
=
atom
->
x
;
int
*
mask
=
atom
->
mask
;
int
nlocal
=
atom
->
nlocal
;
// epsilon is not used, except for book-keeping
for
(
i
=
0
;
i
<
3
;
i
++
)
epsilon
[
i
]
+=
dtq
*
epsilon_dot
[
i
];
// convert pertinent atoms and rigid bodies to lamda coords
if
(
allremap
)
domain
->
x2lamda
(
nlocal
);
else
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
dilate_group_bit
)
domain
->
x2lamda
(
x
[
i
],
x
[
i
]);
}
if
(
nrigid
)
for
(
i
=
0
;
i
<
nrigidfix
;
i
++
)
modify
->
fix
[
rfix
[
i
]]
->
deform
(
0
);
// reset global and local box to new size/shape
for
(
i
=
0
;
i
<
3
;
i
++
)
{
if
(
p_flag
[
i
])
{
oldlo
=
domain
->
boxlo
[
i
];
oldhi
=
domain
->
boxhi
[
i
];
ctr
=
0.5
*
(
oldlo
+
oldhi
);
expfac
=
exp
(
dtq
*
epsilon_dot
[
i
]);
domain
->
boxlo
[
i
]
=
(
oldlo
-
ctr
)
*
expfac
+
ctr
;
domain
->
boxhi
[
i
]
=
(
oldhi
-
ctr
)
*
expfac
+
ctr
;
}
}
domain
->
set_global_box
();
domain
->
set_local_box
();
// convert pertinent atoms and rigid bodies back to box coords
if
(
allremap
)
domain
->
lamda2x
(
nlocal
);
else
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
dilate_group_bit
)
domain
->
lamda2x
(
x
[
i
],
x
[
i
]);
}
if
(
nrigid
)
for
(
i
=
0
;
i
<
nrigidfix
;
i
++
)
modify
->
fix
[
rfix
[
i
]]
->
deform
(
1
);
}
/* ----------------------------------------------------------------------
compute target temperature and kinetic energy
-----------------------------------------------------------------------*/
void
FixRigidNH
::
compute_temp_target
()
{
double
delta
=
update
->
ntimestep
-
update
->
beginstep
;
if
(
delta
!=
0.0
)
delta
/=
update
->
endstep
-
update
->
beginstep
;
t_target
=
t_start
+
delta
*
(
t_stop
-
t_start
);
}
/* ----------------------------------------------------------------------
compute hydrostatic target pressure
-----------------------------------------------------------------------*/
void
FixRigidNH
::
compute_press_target
()
{
double
delta
=
update
->
ntimestep
-
update
->
beginstep
;
if
(
delta
!=
0.0
)
delta
/=
update
->
endstep
-
update
->
beginstep
;
p_hydro
=
0.0
;
for
(
int
i
=
0
;
i
<
3
;
i
++
)
if
(
p_flag
[
i
])
{
p_target
[
i
]
=
p_start
[
i
]
+
delta
*
(
p_stop
[
i
]
-
p_start
[
i
]);
p_hydro
+=
p_target
[
i
];
}
p_hydro
/=
pdim
;
}
/* ----------------------------------------------------------------------
update epsilon_dot
-----------------------------------------------------------------------*/
void
FixRigidNH
::
nh_epsilon_dot
()
{
int
i
;
double
volume
,
scale
,
f_epsilon
;
if
(
dimension
==
2
)
volume
=
domain
->
xprd
*
domain
->
yprd
;
else
volume
=
domain
->
xprd
*
domain
->
yprd
*
domain
->
zprd
;
// MTK terms
mtk_term1
=
(
akin_t
+
akin_r
)
*
mvv2e
/
g_f
;
scale
=
exp
(
-
1.0
*
dtq
*
eta_dot_b
[
0
]);
for
(
i
=
0
;
i
<
3
;
i
++
)
if
(
p_flag
[
i
])
{
f_epsilon
=
(
p_current
[
i
]
-
p_hydro
)
*
volume
/
nktv2p
+
mtk_term1
;
f_epsilon
/=
epsilon_mass
[
i
];
epsilon_dot
[
i
]
+=
dtq
*
f_epsilon
;
epsilon_dot
[
i
]
*=
scale
;
}
mtk_term2
=
0.0
;
for
(
i
=
0
;
i
<
3
;
i
++
)
if
(
p_flag
[
i
])
mtk_term2
+=
epsilon_dot
[
i
];
mtk_term2
/=
g_f
;
}
/* ----------------------------------------------------------------------
pack entire state of Fix into one write
------------------------------------------------------------------------- */
void
FixRigidNH
::
write_restart
(
FILE
*
fp
)
{
if
(
tstat_flag
==
0
&&
pstat_flag
==
0
)
return
;
int
nsize
=
2
;
// tstat_flag and pstat_flag
if
(
tstat_flag
)
{
nsize
+=
1
;
// t_chain
nsize
+=
4
*
t_chain
;
// eta_t, eta_r, eta_dot_t, eta_dot_r
}
if
(
pstat_flag
)
{
nsize
+=
7
;
// p_chain, epsilon(3) and epsilon_dot(3)
nsize
+=
2
*
p_chain
;
}
double
*
list
;
memory
->
create
(
list
,
nsize
,
"rigid_nh:list"
);
int
n
=
0
;
list
[
n
++
]
=
tstat_flag
;
if
(
tstat_flag
)
{
list
[
n
++
]
=
t_chain
;
for
(
int
i
=
0
;
i
<
t_chain
;
i
++
)
{
list
[
n
++
]
=
eta_t
[
i
];
list
[
n
++
]
=
eta_r
[
i
];
list
[
n
++
]
=
eta_dot_t
[
i
];
list
[
n
++
]
=
eta_dot_r
[
i
];
}
}
list
[
n
++
]
=
pstat_flag
;
if
(
pstat_flag
)
{
list
[
n
++
]
=
epsilon
[
0
];
list
[
n
++
]
=
epsilon
[
1
];
list
[
n
++
]
=
epsilon
[
2
];
list
[
n
++
]
=
epsilon_dot
[
0
];
list
[
n
++
]
=
epsilon_dot
[
1
];
list
[
n
++
]
=
epsilon_dot
[
2
];
list
[
n
++
]
=
p_chain
;
for
(
int
i
=
0
;
i
<
p_chain
;
i
++
)
{
list
[
n
++
]
=
eta_b
[
i
];
list
[
n
++
]
=
eta_dot_b
[
i
];
}
}
if
(
comm
->
me
==
0
)
{
int
size
=
(
nsize
)
*
sizeof
(
double
);
fwrite
(
&
size
,
sizeof
(
int
),
1
,
fp
);
fwrite
(
list
,
sizeof
(
double
),
nsize
,
fp
);
}
memory
->
destroy
(
list
);
}
/* ----------------------------------------------------------------------
use state info from restart file to restart the Fix
------------------------------------------------------------------------- */
void
FixRigidNH
::
restart
(
char
*
buf
)
{
int
n
=
0
;
double
*
list
=
(
double
*
)
buf
;
int
flag
=
static_cast
<
int
>
(
list
[
n
++
]);
if
(
flag
)
{
int
m
=
static_cast
<
int
>
(
list
[
n
++
]);
if
(
tstat_flag
&&
m
==
t_chain
)
{
for
(
int
i
=
0
;
i
<
t_chain
;
i
++
)
{
eta_t
[
i
]
=
list
[
n
++
];
eta_r
[
i
]
=
list
[
n
++
];
eta_dot_t
[
i
]
=
list
[
n
++
];
eta_dot_r
[
i
]
=
list
[
n
++
];
}
}
else
n
+=
4
*
m
;
}
flag
=
static_cast
<
int
>
(
list
[
n
++
]);
if
(
flag
)
{
epsilon
[
0
]
=
list
[
n
++
];
epsilon
[
1
]
=
list
[
n
++
];
epsilon
[
2
]
=
list
[
n
++
];
epsilon_dot
[
0
]
=
list
[
n
++
];
epsilon_dot
[
1
]
=
list
[
n
++
];
epsilon_dot
[
2
]
=
list
[
n
++
];
int
m
=
static_cast
<
int
>
(
list
[
n
++
]);
if
(
pstat_flag
&&
m
==
p_chain
)
{
for
(
int
i
=
0
;
i
<
p_chain
;
i
++
)
{
eta_b
[
i
]
=
list
[
n
++
];
eta_dot_b
[
i
]
=
list
[
n
++
];
}
}
else
n
+=
2
*
m
;
}
}
/* ---------------------------------------------------------------------- */
int
FixRigidNH
::
modify_param
(
int
narg
,
char
**
arg
)
{
if
(
strcmp
(
arg
[
0
],
"temp"
)
==
0
)
{
if
(
narg
<
2
)
error
->
all
(
FLERR
,
"Illegal fix_modify command"
);
if
(
!
pstat_flag
)
error
->
all
(
FLERR
,
"Illegal fix_modify command"
);
if
(
tcomputeflag
)
{
modify
->
delete_compute
(
id_temp
);
tcomputeflag
=
0
;
}
delete
[]
id_temp
;
int
n
=
strlen
(
arg
[
1
])
+
1
;
id_temp
=
new
char
[
n
];
strcpy
(
id_temp
,
arg
[
1
]);
int
icompute
=
modify
->
find_compute
(
arg
[
1
]);
if
(
icompute
<
0
)
error
->
all
(
FLERR
,
"Could not find fix_modify temperature ID"
);
temperature
=
modify
->
compute
[
icompute
];
if
(
temperature
->
tempflag
==
0
)
error
->
all
(
FLERR
,
"Fix_modify temperature ID does not compute temperature"
);
if
(
temperature
->
igroup
!=
0
&&
comm
->
me
==
0
)
error
->
warning
(
FLERR
,
"Temperature for fix modify is not for group all"
);
// reset id_temp of pressure to new temperature ID
if
(
pstat_flag
)
{
icompute
=
modify
->
find_compute
(
id_press
);
if
(
icompute
<
0
)
error
->
all
(
FLERR
,
"Pressure ID for fix modify does not exist"
);
modify
->
compute
[
icompute
]
->
reset_extra_compute_fix
(
id_temp
);
}
return
2
;
}
else
if
(
strcmp
(
arg
[
0
],
"press"
)
==
0
)
{
if
(
narg
<
2
)
error
->
all
(
FLERR
,
"Illegal fix_modify command"
);
if
(
!
pstat_flag
)
error
->
all
(
FLERR
,
"Illegal fix_modify command"
);
if
(
pcomputeflag
)
{
modify
->
delete_compute
(
id_press
);
pcomputeflag
=
0
;
}
delete
[]
id_press
;
int
n
=
strlen
(
arg
[
1
])
+
1
;
id_press
=
new
char
[
n
];
strcpy
(
id_press
,
arg
[
1
]);
int
icompute
=
modify
->
find_compute
(
arg
[
1
]);
if
(
icompute
<
0
)
error
->
all
(
FLERR
,
"Could not find fix_modify pressure ID"
);
pressure
=
modify
->
compute
[
icompute
];
if
(
pressure
->
pressflag
==
0
)
error
->
all
(
FLERR
,
"Fix_modify pressure ID does not compute pressure"
);
return
2
;
}
return
0
;
}
/* ---------------------------------------------------------------------- */
void
FixRigidNH
::
allocate_chain
()
{
if
(
tstat_flag
)
{
q_t
=
new
double
[
t_chain
];
q_r
=
new
double
[
t_chain
];
eta_t
=
new
double
[
t_chain
];
eta_r
=
new
double
[
t_chain
];
eta_dot_t
=
new
double
[
t_chain
];
eta_dot_r
=
new
double
[
t_chain
];
f_eta_t
=
new
double
[
t_chain
];
f_eta_r
=
new
double
[
t_chain
];
}
if
(
pstat_flag
)
{
q_b
=
new
double
[
p_chain
];
eta_b
=
new
double
[
p_chain
];
eta_dot_b
=
new
double
[
p_chain
];
f_eta_b
=
new
double
[
p_chain
];
}
}
/* ---------------------------------------------------------------------- */
void
FixRigidNH
::
reset_target
(
double
t_new
)
{
t_start
=
t_stop
=
t_new
;
}
/* ---------------------------------------------------------------------- */
void
FixRigidNH
::
allocate_order
()
{
w
=
new
double
[
t_order
];
wdti1
=
new
double
[
t_order
];
wdti2
=
new
double
[
t_order
];
wdti4
=
new
double
[
t_order
];
}
/* ---------------------------------------------------------------------- */
void
FixRigidNH
::
deallocate_chain
()
{
if
(
tstat_flag
)
{
delete
[]
q_t
;
delete
[]
q_r
;
delete
[]
eta_t
;
delete
[]
eta_r
;
delete
[]
eta_dot_t
;
delete
[]
eta_dot_r
;
delete
[]
f_eta_t
;
delete
[]
f_eta_r
;
}
if
(
pstat_flag
)
{
delete
[]
q_b
;
delete
[]
eta_b
;
delete
[]
eta_dot_b
;
delete
[]
f_eta_b
;
}
}
/* ---------------------------------------------------------------------- */
void
FixRigidNH
::
deallocate_order
()
{
delete
[]
w
;
delete
[]
wdti1
;
delete
[]
wdti2
;
delete
[]
wdti4
;
}
Event Timeline
Log In to Comment