Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F106921119
thr_omp.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, Apr 2, 13:29
Size
36 KB
Mime Type
text/x-c
Expires
Fri, Apr 4, 13:29 (2 d)
Engine
blob
Format
Raw Data
Handle
25306209
Attached To
rLAMMPS lammps
thr_omp.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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
OpenMP based threading support for LAMMPS
Contributing author: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "modify.h"
#include "neighbor.h"
#include "thr_omp.h"
#include "pair.h"
#include "bond.h"
#include "angle.h"
#include "dihedral.h"
#include "improper.h"
#include "kspace.h"
#include "math_const.h"
#include <string.h>
using
namespace
LAMMPS_NS
;
using
namespace
MathConst
;
/* ---------------------------------------------------------------------- */
ThrOMP
::
ThrOMP
(
LAMMPS
*
ptr
,
int
style
)
:
lmp
(
ptr
),
fix
(
NULL
),
thr_style
(
style
),
thr_error
(
0
)
{
// register fix omp with this class
int
ifix
=
lmp
->
modify
->
find_fix
(
"package_omp"
);
if
(
ifix
<
0
)
lmp
->
error
->
all
(
FLERR
,
"The 'package omp' command is required for /omp styles"
);
fix
=
static_cast
<
FixOMP
*>
(
lmp
->
modify
->
fix
[
ifix
]);
}
/* ---------------------------------------------------------------------- */
ThrOMP
::~
ThrOMP
()
{
// nothing to do?
}
/* ----------------------------------------------------------------------
Hook up per thread per atom arrays into the tally infrastructure
---------------------------------------------------------------------- */
void
ThrOMP
::
ev_setup_thr
(
int
eflag
,
int
vflag
,
int
nall
,
double
*
eatom
,
double
**
vatom
,
ThrData
*
thr
)
{
const
int
tid
=
thr
->
get_tid
();
if
(
tid
==
0
)
thr_error
=
0
;
if
(
thr_style
&
THR_PAIR
)
{
if
(
eflag
&
2
)
{
thr
->
eatom_pair
=
eatom
+
tid
*
nall
;
if
(
nall
>
0
)
memset
(
&
(
thr
->
eatom_pair
[
0
]),
0
,
nall
*
sizeof
(
double
));
}
if
(
vflag
&
4
)
{
thr
->
vatom_pair
=
vatom
+
tid
*
nall
;
if
(
nall
>
0
)
memset
(
&
(
thr
->
vatom_pair
[
0
][
0
]),
0
,
nall
*
6
*
sizeof
(
double
));
}
}
if
(
thr_style
&
THR_BOND
)
{
if
(
eflag
&
2
)
{
thr
->
eatom_bond
=
eatom
+
tid
*
nall
;
if
(
nall
>
0
)
memset
(
&
(
thr
->
eatom_bond
[
0
]),
0
,
nall
*
sizeof
(
double
));
}
if
(
vflag
&
4
)
{
thr
->
vatom_bond
=
vatom
+
tid
*
nall
;
if
(
nall
>
0
)
memset
(
&
(
thr
->
vatom_bond
[
0
][
0
]),
0
,
nall
*
6
*
sizeof
(
double
));
}
}
if
(
thr_style
&
THR_ANGLE
)
{
if
(
eflag
&
2
)
{
thr
->
eatom_angle
=
eatom
+
tid
*
nall
;
if
(
nall
>
0
)
memset
(
&
(
thr
->
eatom_angle
[
0
]),
0
,
nall
*
sizeof
(
double
));
}
if
(
vflag
&
4
)
{
thr
->
vatom_angle
=
vatom
+
tid
*
nall
;
if
(
nall
>
0
)
memset
(
&
(
thr
->
vatom_angle
[
0
][
0
]),
0
,
nall
*
6
*
sizeof
(
double
));
}
}
if
(
thr_style
&
THR_DIHEDRAL
)
{
if
(
eflag
&
2
)
{
thr
->
eatom_dihed
=
eatom
+
tid
*
nall
;
if
(
nall
>
0
)
memset
(
&
(
thr
->
eatom_dihed
[
0
]),
0
,
nall
*
sizeof
(
double
));
}
if
(
vflag
&
4
)
{
thr
->
vatom_dihed
=
vatom
+
tid
*
nall
;
if
(
nall
>
0
)
memset
(
&
(
thr
->
vatom_dihed
[
0
][
0
]),
0
,
nall
*
6
*
sizeof
(
double
));
}
}
if
(
thr_style
&
THR_IMPROPER
)
{
if
(
eflag
&
2
)
{
thr
->
eatom_imprp
=
eatom
+
tid
*
nall
;
if
(
nall
>
0
)
memset
(
&
(
thr
->
eatom_imprp
[
0
]),
0
,
nall
*
sizeof
(
double
));
}
if
(
vflag
&
4
)
{
thr
->
vatom_imprp
=
vatom
+
tid
*
nall
;
if
(
nall
>
0
)
memset
(
&
(
thr
->
vatom_imprp
[
0
][
0
]),
0
,
nall
*
6
*
sizeof
(
double
));
}
}
// nothing to do for THR_KSPACE
}
/* ----------------------------------------------------------------------
Reduce per thread data into the regular structures
Reduction of global properties is serialized with a "critical"
directive, so that only one thread at a time will access the
global variables. Since we are not synchronized, this should
come with little overhead. The reduction of per-atom properties
in contrast is parallelized over threads in the same way as forces.
---------------------------------------------------------------------- */
void
ThrOMP
::
reduce_thr
(
void
*
style
,
const
int
eflag
,
const
int
vflag
,
ThrData
*
const
thr
)
{
const
int
nlocal
=
lmp
->
atom
->
nlocal
;
const
int
nghost
=
lmp
->
atom
->
nghost
;
const
int
nall
=
nlocal
+
nghost
;
const
int
nfirst
=
lmp
->
atom
->
nfirst
;
const
int
nthreads
=
lmp
->
comm
->
nthreads
;
const
int
evflag
=
eflag
|
vflag
;
const
int
tid
=
thr
->
get_tid
();
double
**
f
=
lmp
->
atom
->
f
;
double
**
x
=
lmp
->
atom
->
x
;
int
need_force_reduce
=
1
;
if
(
evflag
)
sync_threads
();
switch
(
thr_style
)
{
case
THR_PAIR:
{
Pair
*
const
pair
=
lmp
->
force
->
pair
;
if
(
pair
->
vflag_fdotr
)
{
// this is a non-hybrid pair style. compute per thread fdotr
if
(
fix
->
last_pair_hybrid
==
NULL
)
{
if
(
lmp
->
neighbor
->
includegroup
==
0
)
thr
->
virial_fdotr_compute
(
x
,
nlocal
,
nghost
,
-
1
);
else
thr
->
virial_fdotr_compute
(
x
,
nlocal
,
nghost
,
nfirst
);
}
else
{
if
(
style
==
fix
->
last_pair_hybrid
)
{
// pair_style hybrid will compute fdotr for us
// but we first need to reduce the forces
data_reduce_thr
(
&
(
f
[
0
][
0
]),
nall
,
nthreads
,
3
,
tid
);
fix
->
did_reduce
();
need_force_reduce
=
0
;
}
}
}
if
(
evflag
)
{
#if defined(_OPENMP)
#pragma omp critical
#endif
{
if
(
eflag
&
1
)
{
pair
->
eng_vdwl
+=
thr
->
eng_vdwl
;
pair
->
eng_coul
+=
thr
->
eng_coul
;
thr
->
eng_vdwl
=
0.0
;
thr
->
eng_coul
=
0.0
;
}
if
(
vflag
&
3
)
for
(
int
i
=
0
;
i
<
6
;
++
i
)
{
pair
->
virial
[
i
]
+=
thr
->
virial_pair
[
i
];
thr
->
virial_pair
[
i
]
=
0.0
;
}
}
if
(
eflag
&
2
)
{
data_reduce_thr
(
&
(
pair
->
eatom
[
0
]),
nall
,
nthreads
,
1
,
tid
);
}
if
(
vflag
&
4
)
{
data_reduce_thr
(
&
(
pair
->
vatom
[
0
][
0
]),
nall
,
nthreads
,
6
,
tid
);
}
}
}
break
;
case
THR_BOND:
if
(
evflag
)
{
Bond
*
const
bond
=
lmp
->
force
->
bond
;
#if defined(_OPENMP)
#pragma omp critical
#endif
{
if
(
eflag
&
1
)
{
bond
->
energy
+=
thr
->
eng_bond
;
thr
->
eng_bond
=
0.0
;
}
if
(
vflag
&
3
)
{
for
(
int
i
=
0
;
i
<
6
;
++
i
)
{
bond
->
virial
[
i
]
+=
thr
->
virial_bond
[
i
];
thr
->
virial_bond
[
i
]
=
0.0
;
}
}
}
if
(
eflag
&
2
)
{
data_reduce_thr
(
&
(
bond
->
eatom
[
0
]),
nall
,
nthreads
,
1
,
tid
);
}
if
(
vflag
&
4
)
{
data_reduce_thr
(
&
(
bond
->
vatom
[
0
][
0
]),
nall
,
nthreads
,
6
,
tid
);
}
}
break
;
case
THR_ANGLE:
if
(
evflag
)
{
Angle
*
const
angle
=
lmp
->
force
->
angle
;
#if defined(_OPENMP)
#pragma omp critical
#endif
{
if
(
eflag
&
1
)
{
angle
->
energy
+=
thr
->
eng_angle
;
thr
->
eng_angle
=
0.0
;
}
if
(
vflag
&
3
)
{
for
(
int
i
=
0
;
i
<
6
;
++
i
)
{
angle
->
virial
[
i
]
+=
thr
->
virial_angle
[
i
];
thr
->
virial_angle
[
i
]
=
0.0
;
}
}
}
if
(
eflag
&
2
)
{
data_reduce_thr
(
&
(
angle
->
eatom
[
0
]),
nall
,
nthreads
,
1
,
tid
);
}
if
(
vflag
&
4
)
{
data_reduce_thr
(
&
(
angle
->
vatom
[
0
][
0
]),
nall
,
nthreads
,
6
,
tid
);
}
}
break
;
case
THR_DIHEDRAL:
if
(
evflag
)
{
Dihedral
*
const
dihedral
=
lmp
->
force
->
dihedral
;
#if defined(_OPENMP)
#pragma omp critical
#endif
{
if
(
eflag
&
1
)
{
dihedral
->
energy
+=
thr
->
eng_dihed
;
thr
->
eng_dihed
=
0.0
;
}
if
(
vflag
&
3
)
{
for
(
int
i
=
0
;
i
<
6
;
++
i
)
{
dihedral
->
virial
[
i
]
+=
thr
->
virial_dihed
[
i
];
thr
->
virial_dihed
[
i
]
=
0.0
;
}
}
}
if
(
eflag
&
2
)
{
data_reduce_thr
(
&
(
dihedral
->
eatom
[
0
]),
nall
,
nthreads
,
1
,
tid
);
}
if
(
vflag
&
4
)
{
data_reduce_thr
(
&
(
dihedral
->
vatom
[
0
][
0
]),
nall
,
nthreads
,
6
,
tid
);
}
}
break
;
case
THR_DIHEDRAL
|
THR_CHARMM:
// special case for CHARMM dihedrals
if
(
evflag
)
{
Dihedral
*
const
dihedral
=
lmp
->
force
->
dihedral
;
Pair
*
const
pair
=
lmp
->
force
->
pair
;
#if defined(_OPENMP)
#pragma omp critical
#endif
{
if
(
eflag
&
1
)
{
dihedral
->
energy
+=
thr
->
eng_dihed
;
pair
->
eng_vdwl
+=
thr
->
eng_vdwl
;
pair
->
eng_coul
+=
thr
->
eng_coul
;
thr
->
eng_dihed
=
0.0
;
thr
->
eng_vdwl
=
0.0
;
thr
->
eng_coul
=
0.0
;
}
if
(
vflag
&
3
)
{
for
(
int
i
=
0
;
i
<
6
;
++
i
)
{
dihedral
->
virial
[
i
]
+=
thr
->
virial_dihed
[
i
];
pair
->
virial
[
i
]
+=
thr
->
virial_pair
[
i
];
thr
->
virial_dihed
[
i
]
=
0.0
;
thr
->
virial_pair
[
i
]
=
0.0
;
}
}
}
if
(
eflag
&
2
)
{
data_reduce_thr
(
&
(
dihedral
->
eatom
[
0
]),
nall
,
nthreads
,
1
,
tid
);
data_reduce_thr
(
&
(
pair
->
eatom
[
0
]),
nall
,
nthreads
,
1
,
tid
);
}
if
(
vflag
&
4
)
{
data_reduce_thr
(
&
(
dihedral
->
vatom
[
0
][
0
]),
nall
,
nthreads
,
6
,
tid
);
data_reduce_thr
(
&
(
pair
->
vatom
[
0
][
0
]),
nall
,
nthreads
,
6
,
tid
);
}
}
break
;
case
THR_IMPROPER:
if
(
evflag
)
{
Improper
*
improper
=
lmp
->
force
->
improper
;
#if defined(_OPENMP)
#pragma omp critical
#endif
{
if
(
eflag
&
1
)
{
improper
->
energy
+=
thr
->
eng_imprp
;
thr
->
eng_imprp
=
0.0
;
}
if
(
vflag
&
3
)
{
for
(
int
i
=
0
;
i
<
6
;
++
i
)
{
improper
->
virial
[
i
]
+=
thr
->
virial_imprp
[
i
];
thr
->
virial_imprp
[
i
]
=
0.0
;
}
}
}
if
(
eflag
&
2
)
{
data_reduce_thr
(
&
(
improper
->
eatom
[
0
]),
nall
,
nthreads
,
1
,
tid
);
}
if
(
vflag
&
4
)
{
data_reduce_thr
(
&
(
improper
->
vatom
[
0
][
0
]),
nall
,
nthreads
,
6
,
tid
);
}
}
break
;
case
THR_KSPACE:
// nothing to do. XXX may need to add support for per-atom info
break
;
case
THR_INTGR:
// nothing to do
break
;
default
:
printf
(
"tid:%d unhandled thr_style case %d
\n
"
,
tid
,
thr_style
);
break
;
}
if
(
style
==
fix
->
last_omp_style
)
{
if
(
need_force_reduce
)
{
data_reduce_thr
(
&
(
f
[
0
][
0
]),
nall
,
nthreads
,
3
,
tid
);
fix
->
did_reduce
();
}
if
(
lmp
->
atom
->
torque
)
data_reduce_thr
(
&
(
lmp
->
atom
->
torque
[
0
][
0
]),
nall
,
nthreads
,
3
,
tid
);
}
}
/* ----------------------------------------------------------------------
tally eng_vdwl and eng_coul into per thread global and per-atom accumulators
------------------------------------------------------------------------- */
void
ThrOMP
::
e_tally_thr
(
Pair
*
const
pair
,
const
int
i
,
const
int
j
,
const
int
nlocal
,
const
int
newton_pair
,
const
double
evdwl
,
const
double
ecoul
,
ThrData
*
const
thr
)
{
if
(
pair
->
eflag_global
)
{
if
(
newton_pair
)
{
thr
->
eng_vdwl
+=
evdwl
;
thr
->
eng_coul
+=
ecoul
;
}
else
{
const
double
evdwlhalf
=
0.5
*
evdwl
;
const
double
ecoulhalf
=
0.5
*
ecoul
;
if
(
i
<
nlocal
)
{
thr
->
eng_vdwl
+=
evdwlhalf
;
thr
->
eng_coul
+=
ecoulhalf
;
}
if
(
j
<
nlocal
)
{
thr
->
eng_vdwl
+=
evdwlhalf
;
thr
->
eng_coul
+=
ecoulhalf
;
}
}
}
if
(
pair
->
eflag_atom
)
{
const
double
epairhalf
=
0.5
*
(
evdwl
+
ecoul
);
if
(
newton_pair
||
i
<
nlocal
)
thr
->
eatom_pair
[
i
]
+=
epairhalf
;
if
(
newton_pair
||
j
<
nlocal
)
thr
->
eatom_pair
[
j
]
+=
epairhalf
;
}
}
/* helper functions */
static
void
v_tally
(
double
*
const
vout
,
const
double
*
const
vin
)
{
vout
[
0
]
+=
vin
[
0
];
vout
[
1
]
+=
vin
[
1
];
vout
[
2
]
+=
vin
[
2
];
vout
[
3
]
+=
vin
[
3
];
vout
[
4
]
+=
vin
[
4
];
vout
[
5
]
+=
vin
[
5
];
}
static
void
v_tally
(
double
*
const
vout
,
const
double
scale
,
const
double
*
const
vin
)
{
vout
[
0
]
+=
scale
*
vin
[
0
];
vout
[
1
]
+=
scale
*
vin
[
1
];
vout
[
2
]
+=
scale
*
vin
[
2
];
vout
[
3
]
+=
scale
*
vin
[
3
];
vout
[
4
]
+=
scale
*
vin
[
4
];
vout
[
5
]
+=
scale
*
vin
[
5
];
}
/* ----------------------------------------------------------------------
tally virial into per thread global and per-atom accumulators
------------------------------------------------------------------------- */
void
ThrOMP
::
v_tally_thr
(
Pair
*
const
pair
,
const
int
i
,
const
int
j
,
const
int
nlocal
,
const
int
newton_pair
,
const
double
*
const
v
,
ThrData
*
const
thr
)
{
if
(
pair
->
vflag_global
)
{
double
*
const
va
=
thr
->
virial_pair
;
if
(
newton_pair
)
{
v_tally
(
va
,
v
);
}
else
{
if
(
i
<
nlocal
)
v_tally
(
va
,
0.5
,
v
);
if
(
j
<
nlocal
)
v_tally
(
va
,
0.5
,
v
);
}
}
if
(
pair
->
vflag_atom
)
{
if
(
newton_pair
||
i
<
nlocal
)
{
double
*
const
va
=
thr
->
vatom_pair
[
i
];
v_tally
(
va
,
0.5
,
v
);
}
if
(
newton_pair
||
j
<
nlocal
)
{
double
*
const
va
=
thr
->
vatom_pair
[
j
];
v_tally
(
va
,
0.5
,
v
);
}
}
}
/* ----------------------------------------------------------------------
tally eng_vdwl and virial into per thread global and per-atom accumulators
need i < nlocal test since called by bond_quartic and dihedral_charmm
------------------------------------------------------------------------- */
void
ThrOMP
::
ev_tally_thr
(
Pair
*
const
pair
,
const
int
i
,
const
int
j
,
const
int
nlocal
,
const
int
newton_pair
,
const
double
evdwl
,
const
double
ecoul
,
const
double
fpair
,
const
double
delx
,
const
double
dely
,
const
double
delz
,
ThrData
*
const
thr
)
{
if
(
pair
->
eflag_either
)
e_tally_thr
(
pair
,
i
,
j
,
nlocal
,
newton_pair
,
evdwl
,
ecoul
,
thr
);
if
(
pair
->
vflag_either
)
{
double
v
[
6
];
v
[
0
]
=
delx
*
delx
*
fpair
;
v
[
1
]
=
dely
*
dely
*
fpair
;
v
[
2
]
=
delz
*
delz
*
fpair
;
v
[
3
]
=
delx
*
dely
*
fpair
;
v
[
4
]
=
delx
*
delz
*
fpair
;
v
[
5
]
=
dely
*
delz
*
fpair
;
v_tally_thr
(
pair
,
i
,
j
,
nlocal
,
newton_pair
,
v
,
thr
);
}
}
/* ----------------------------------------------------------------------
tally eng_vdwl and virial into global and per-atom accumulators
for virial, have delx,dely,delz and fx,fy,fz
------------------------------------------------------------------------- */
void
ThrOMP
::
ev_tally_xyz_thr
(
Pair
*
const
pair
,
const
int
i
,
const
int
j
,
const
int
nlocal
,
const
int
newton_pair
,
const
double
evdwl
,
const
double
ecoul
,
const
double
fx
,
const
double
fy
,
const
double
fz
,
const
double
delx
,
const
double
dely
,
const
double
delz
,
ThrData
*
const
thr
)
{
if
(
pair
->
eflag_either
)
e_tally_thr
(
pair
,
i
,
j
,
nlocal
,
newton_pair
,
evdwl
,
ecoul
,
thr
);
if
(
pair
->
vflag_either
)
{
double
v
[
6
];
v
[
0
]
=
delx
*
fx
;
v
[
1
]
=
dely
*
fy
;
v
[
2
]
=
delz
*
fz
;
v
[
3
]
=
delx
*
fy
;
v
[
4
]
=
delx
*
fz
;
v
[
5
]
=
dely
*
fz
;
v_tally_thr
(
pair
,
i
,
j
,
nlocal
,
newton_pair
,
v
,
thr
);
}
}
/* ----------------------------------------------------------------------
tally eng_vdwl and virial into global and per-atom accumulators
called by SW and hbond potentials, newton_pair is always on
virial = riFi + rjFj + rkFk = (rj-ri) Fj + (rk-ri) Fk = drji*fj + drki*fk
------------------------------------------------------------------------- */
void
ThrOMP
::
ev_tally3_thr
(
Pair
*
const
pair
,
const
int
i
,
const
int
j
,
const
int
k
,
const
double
evdwl
,
const
double
ecoul
,
const
double
*
const
fj
,
const
double
*
const
fk
,
const
double
*
const
drji
,
const
double
*
const
drki
,
ThrData
*
const
thr
)
{
if
(
pair
->
eflag_either
)
{
if
(
pair
->
eflag_global
)
{
thr
->
eng_vdwl
+=
evdwl
;
thr
->
eng_coul
+=
ecoul
;
}
if
(
pair
->
eflag_atom
)
{
const
double
epairthird
=
THIRD
*
(
evdwl
+
ecoul
);
thr
->
eatom_pair
[
i
]
+=
epairthird
;
thr
->
eatom_pair
[
j
]
+=
epairthird
;
thr
->
eatom_pair
[
k
]
+=
epairthird
;
}
}
if
(
pair
->
vflag_either
)
{
double
v
[
6
];
v
[
0
]
=
drji
[
0
]
*
fj
[
0
]
+
drki
[
0
]
*
fk
[
0
];
v
[
1
]
=
drji
[
1
]
*
fj
[
1
]
+
drki
[
1
]
*
fk
[
1
];
v
[
2
]
=
drji
[
2
]
*
fj
[
2
]
+
drki
[
2
]
*
fk
[
2
];
v
[
3
]
=
drji
[
0
]
*
fj
[
1
]
+
drki
[
0
]
*
fk
[
1
];
v
[
4
]
=
drji
[
0
]
*
fj
[
2
]
+
drki
[
0
]
*
fk
[
2
];
v
[
5
]
=
drji
[
1
]
*
fj
[
2
]
+
drki
[
1
]
*
fk
[
2
];
if
(
pair
->
vflag_global
)
v_tally
(
thr
->
virial_pair
,
v
);
if
(
pair
->
vflag_atom
)
{
v_tally
(
thr
->
vatom_pair
[
i
],
THIRD
,
v
);
v_tally
(
thr
->
vatom_pair
[
j
],
THIRD
,
v
);
v_tally
(
thr
->
vatom_pair
[
k
],
THIRD
,
v
);
}
}
}
/* ----------------------------------------------------------------------
tally eng_vdwl and virial into global and per-atom accumulators
called by AIREBO potential, newton_pair is always on
------------------------------------------------------------------------- */
void
ThrOMP
::
ev_tally4_thr
(
Pair
*
const
pair
,
const
int
i
,
const
int
j
,
const
int
k
,
const
int
m
,
const
double
evdwl
,
const
double
*
const
fi
,
const
double
*
const
fj
,
const
double
*
const
fk
,
const
double
*
const
drim
,
const
double
*
const
drjm
,
const
double
*
const
drkm
,
ThrData
*
const
thr
)
{
double
v
[
6
];
if
(
pair
->
eflag_either
)
{
if
(
pair
->
eflag_global
)
thr
->
eng_vdwl
+=
evdwl
;
if
(
pair
->
eflag_atom
)
{
const
double
epairfourth
=
0.25
*
evdwl
;
thr
->
eatom_pair
[
i
]
+=
epairfourth
;
thr
->
eatom_pair
[
j
]
+=
epairfourth
;
thr
->
eatom_pair
[
k
]
+=
epairfourth
;
thr
->
eatom_pair
[
m
]
+=
epairfourth
;
}
}
if
(
pair
->
vflag_atom
)
{
v
[
0
]
=
0.25
*
(
drim
[
0
]
*
fi
[
0
]
+
drjm
[
0
]
*
fj
[
0
]
+
drkm
[
0
]
*
fk
[
0
]);
v
[
1
]
=
0.25
*
(
drim
[
1
]
*
fi
[
1
]
+
drjm
[
1
]
*
fj
[
1
]
+
drkm
[
1
]
*
fk
[
1
]);
v
[
2
]
=
0.25
*
(
drim
[
2
]
*
fi
[
2
]
+
drjm
[
2
]
*
fj
[
2
]
+
drkm
[
2
]
*
fk
[
2
]);
v
[
3
]
=
0.25
*
(
drim
[
0
]
*
fi
[
1
]
+
drjm
[
0
]
*
fj
[
1
]
+
drkm
[
0
]
*
fk
[
1
]);
v
[
4
]
=
0.25
*
(
drim
[
0
]
*
fi
[
2
]
+
drjm
[
0
]
*
fj
[
2
]
+
drkm
[
0
]
*
fk
[
2
]);
v
[
5
]
=
0.25
*
(
drim
[
1
]
*
fi
[
2
]
+
drjm
[
1
]
*
fj
[
2
]
+
drkm
[
1
]
*
fk
[
2
]);
v_tally
(
thr
->
vatom_pair
[
i
],
v
);
v_tally
(
thr
->
vatom_pair
[
j
],
v
);
v_tally
(
thr
->
vatom_pair
[
k
],
v
);
v_tally
(
thr
->
vatom_pair
[
m
],
v
);
}
}
/* ----------------------------------------------------------------------
tally ecoul and virial into each of n atoms in list
called by TIP4P potential, newton_pair is always on
changes v values by dividing by n
------------------------------------------------------------------------- */
void
ThrOMP
::
ev_tally_list_thr
(
Pair
*
const
pair
,
const
int
key
,
const
int
*
const
list
,
const
double
*
const
v
,
const
double
ecoul
,
const
double
alpha
,
ThrData
*
const
thr
)
{
int
i
;
if
(
pair
->
eflag_either
)
{
if
(
pair
->
eflag_global
)
thr
->
eng_coul
+=
ecoul
;
if
(
pair
->
eflag_atom
)
{
if
(
key
==
0
)
{
thr
->
eatom_pair
[
list
[
0
]]
+=
0.5
*
ecoul
;
thr
->
eatom_pair
[
list
[
1
]]
+=
0.5
*
ecoul
;
}
else
if
(
key
==
1
)
{
thr
->
eatom_pair
[
list
[
0
]]
+=
0.5
*
ecoul
*
(
1
-
alpha
);
thr
->
eatom_pair
[
list
[
1
]]
+=
0.25
*
ecoul
*
alpha
;
thr
->
eatom_pair
[
list
[
2
]]
+=
0.25
*
ecoul
*
alpha
;
thr
->
eatom_pair
[
list
[
3
]]
+=
0.5
*
ecoul
;
}
else
if
(
key
==
2
)
{
thr
->
eatom_pair
[
list
[
0
]]
+=
0.5
*
ecoul
;
thr
->
eatom_pair
[
list
[
1
]]
+=
0.5
*
ecoul
*
(
1
-
alpha
);
thr
->
eatom_pair
[
list
[
2
]]
+=
0.25
*
ecoul
*
alpha
;
thr
->
eatom_pair
[
list
[
3
]]
+=
0.25
*
ecoul
*
alpha
;
}
else
{
thr
->
eatom_pair
[
list
[
0
]]
+=
0.5
*
ecoul
*
(
1
-
alpha
);
thr
->
eatom_pair
[
list
[
1
]]
+=
0.25
*
ecoul
*
alpha
;
thr
->
eatom_pair
[
list
[
2
]]
+=
0.25
*
ecoul
*
alpha
;
thr
->
eatom_pair
[
list
[
3
]]
+=
0.5
*
ecoul
*
(
1
-
alpha
);
thr
->
eatom_pair
[
list
[
4
]]
+=
0.25
*
ecoul
*
alpha
;
thr
->
eatom_pair
[
list
[
5
]]
+=
0.25
*
ecoul
*
alpha
;
}
}
}
if
(
pair
->
vflag_either
)
{
if
(
pair
->
vflag_global
)
v_tally
(
thr
->
virial_pair
,
v
);
if
(
pair
->
vflag_atom
)
{
if
(
key
==
0
)
{
for
(
i
=
0
;
i
<=
5
;
i
++
)
{
thr
->
vatom_pair
[
list
[
0
]][
i
]
+=
0.5
*
v
[
i
];
thr
->
vatom_pair
[
list
[
1
]][
i
]
+=
0.5
*
v
[
i
];
}
}
else
if
(
key
==
1
)
{
for
(
i
=
0
;
i
<=
5
;
i
++
)
{
thr
->
vatom_pair
[
list
[
0
]][
i
]
+=
0.5
*
v
[
i
]
*
(
1
-
alpha
);
thr
->
vatom_pair
[
list
[
1
]][
i
]
+=
0.25
*
v
[
i
]
*
alpha
;
thr
->
vatom_pair
[
list
[
2
]][
i
]
+=
0.25
*
v
[
i
]
*
alpha
;
thr
->
vatom_pair
[
list
[
3
]][
i
]
+=
0.5
*
v
[
i
];
}
}
else
if
(
key
==
2
)
{
for
(
i
=
0
;
i
<=
5
;
i
++
)
{
thr
->
vatom_pair
[
list
[
0
]][
i
]
+=
0.5
*
v
[
i
];
thr
->
vatom_pair
[
list
[
1
]][
i
]
+=
0.5
*
v
[
i
]
*
(
1
-
alpha
);
thr
->
vatom_pair
[
list
[
2
]][
i
]
+=
0.25
*
v
[
i
]
*
alpha
;
thr
->
vatom_pair
[
list
[
3
]][
i
]
+=
0.25
*
v
[
i
]
*
alpha
;
}
}
else
{
for
(
i
=
0
;
i
<=
5
;
i
++
)
{
thr
->
vatom_pair
[
list
[
0
]][
i
]
+=
0.5
*
v
[
i
]
*
(
1
-
alpha
);
thr
->
vatom_pair
[
list
[
1
]][
i
]
+=
0.25
*
v
[
i
]
*
alpha
;
thr
->
vatom_pair
[
list
[
2
]][
i
]
+=
0.25
*
v
[
i
]
*
alpha
;
thr
->
vatom_pair
[
list
[
3
]][
i
]
+=
0.5
*
v
[
i
]
*
(
1
-
alpha
);
thr
->
vatom_pair
[
list
[
4
]][
i
]
+=
0.25
*
v
[
i
]
*
alpha
;
thr
->
vatom_pair
[
list
[
5
]][
i
]
+=
0.25
*
v
[
i
]
*
alpha
;
}
}
}
}
}
/* ----------------------------------------------------------------------
tally energy and virial into global and per-atom accumulators
------------------------------------------------------------------------- */
void
ThrOMP
::
ev_tally_thr
(
Bond
*
const
bond
,
const
int
i
,
const
int
j
,
const
int
nlocal
,
const
int
newton_bond
,
const
double
ebond
,
const
double
fbond
,
const
double
delx
,
const
double
dely
,
const
double
delz
,
ThrData
*
const
thr
)
{
if
(
bond
->
eflag_either
)
{
const
double
ebondhalf
=
0.5
*
ebond
;
if
(
newton_bond
)
{
if
(
bond
->
eflag_global
)
thr
->
eng_bond
+=
ebond
;
if
(
bond
->
eflag_atom
)
{
thr
->
eatom_bond
[
i
]
+=
ebondhalf
;
thr
->
eatom_bond
[
j
]
+=
ebondhalf
;
}
}
else
{
if
(
bond
->
eflag_global
)
{
if
(
i
<
nlocal
)
thr
->
eng_bond
+=
ebondhalf
;
if
(
j
<
nlocal
)
thr
->
eng_bond
+=
ebondhalf
;
}
if
(
bond
->
eflag_atom
)
{
if
(
i
<
nlocal
)
thr
->
eatom_bond
[
i
]
+=
ebondhalf
;
if
(
j
<
nlocal
)
thr
->
eatom_bond
[
j
]
+=
ebondhalf
;
}
}
}
if
(
bond
->
vflag_either
)
{
double
v
[
6
];
v
[
0
]
=
delx
*
delx
*
fbond
;
v
[
1
]
=
dely
*
dely
*
fbond
;
v
[
2
]
=
delz
*
delz
*
fbond
;
v
[
3
]
=
delx
*
dely
*
fbond
;
v
[
4
]
=
delx
*
delz
*
fbond
;
v
[
5
]
=
dely
*
delz
*
fbond
;
if
(
bond
->
vflag_global
)
{
if
(
newton_bond
)
v_tally
(
thr
->
virial_bond
,
v
);
else
{
if
(
i
<
nlocal
)
v_tally
(
thr
->
virial_bond
,
0.5
,
v
);
if
(
j
<
nlocal
)
v_tally
(
thr
->
virial_bond
,
0.5
,
v
);
}
}
if
(
bond
->
vflag_atom
)
{
v
[
0
]
*=
0.5
;
v
[
1
]
*=
0.5
;
v
[
2
]
*=
0.5
;
v
[
3
]
*=
0.5
;
v
[
4
]
*=
0.5
;
v
[
5
]
*=
0.5
;
if
(
newton_bond
)
{
v_tally
(
thr
->
vatom_bond
[
i
],
v
);
v_tally
(
thr
->
vatom_bond
[
j
],
v
);
}
else
{
if
(
j
<
nlocal
)
v_tally
(
thr
->
vatom_bond
[
i
],
v
);
if
(
j
<
nlocal
)
v_tally
(
thr
->
vatom_bond
[
j
],
v
);
}
}
}
}
/* ----------------------------------------------------------------------
tally energy and virial into global and per-atom accumulators
virial = r1F1 + r2F2 + r3F3 = (r1-r2) F1 + (r3-r2) F3 = del1*f1 + del2*f3
------------------------------------------------------------------------- */
void
ThrOMP
::
ev_tally_thr
(
Angle
*
const
angle
,
const
int
i
,
const
int
j
,
const
int
k
,
const
int
nlocal
,
const
int
newton_bond
,
const
double
eangle
,
const
double
*
const
f1
,
const
double
*
const
f3
,
const
double
delx1
,
const
double
dely1
,
const
double
delz1
,
const
double
delx2
,
const
double
dely2
,
const
double
delz2
,
ThrData
*
const
thr
)
{
if
(
angle
->
eflag_either
)
{
const
double
eanglethird
=
THIRD
*
eangle
;
if
(
newton_bond
)
{
if
(
angle
->
eflag_global
)
thr
->
eng_angle
+=
eangle
;
if
(
angle
->
eflag_atom
)
{
thr
->
eatom_angle
[
i
]
+=
eanglethird
;
thr
->
eatom_angle
[
j
]
+=
eanglethird
;
thr
->
eatom_angle
[
k
]
+=
eanglethird
;
}
}
else
{
if
(
angle
->
eflag_global
)
{
if
(
i
<
nlocal
)
thr
->
eng_angle
+=
eanglethird
;
if
(
j
<
nlocal
)
thr
->
eng_angle
+=
eanglethird
;
if
(
k
<
nlocal
)
thr
->
eng_angle
+=
eanglethird
;
}
if
(
angle
->
eflag_atom
)
{
if
(
i
<
nlocal
)
thr
->
eatom_angle
[
i
]
+=
eanglethird
;
if
(
j
<
nlocal
)
thr
->
eatom_angle
[
j
]
+=
eanglethird
;
if
(
k
<
nlocal
)
thr
->
eatom_angle
[
k
]
+=
eanglethird
;
}
}
}
if
(
angle
->
vflag_either
)
{
double
v
[
6
];
v
[
0
]
=
delx1
*
f1
[
0
]
+
delx2
*
f3
[
0
];
v
[
1
]
=
dely1
*
f1
[
1
]
+
dely2
*
f3
[
1
];
v
[
2
]
=
delz1
*
f1
[
2
]
+
delz2
*
f3
[
2
];
v
[
3
]
=
delx1
*
f1
[
1
]
+
delx2
*
f3
[
1
];
v
[
4
]
=
delx1
*
f1
[
2
]
+
delx2
*
f3
[
2
];
v
[
5
]
=
dely1
*
f1
[
2
]
+
dely2
*
f3
[
2
];
if
(
angle
->
vflag_global
)
{
if
(
newton_bond
)
{
v_tally
(
thr
->
virial_angle
,
v
);
}
else
{
int
cnt
=
0
;
if
(
i
<
nlocal
)
++
cnt
;
if
(
j
<
nlocal
)
++
cnt
;
if
(
k
<
nlocal
)
++
cnt
;
v_tally
(
thr
->
virial_angle
,
cnt
*
THIRD
,
v
);
}
}
if
(
angle
->
vflag_atom
)
{
v
[
0
]
*=
THIRD
;
v
[
1
]
*=
THIRD
;
v
[
2
]
*=
THIRD
;
v
[
3
]
*=
THIRD
;
v
[
4
]
*=
THIRD
;
v
[
5
]
*=
THIRD
;
if
(
newton_bond
)
{
v_tally
(
thr
->
vatom_angle
[
i
],
v
);
v_tally
(
thr
->
vatom_angle
[
j
],
v
);
v_tally
(
thr
->
vatom_angle
[
k
],
v
);
}
else
{
if
(
j
<
nlocal
)
v_tally
(
thr
->
vatom_angle
[
i
],
v
);
if
(
j
<
nlocal
)
v_tally
(
thr
->
vatom_angle
[
j
],
v
);
if
(
k
<
nlocal
)
v_tally
(
thr
->
vatom_angle
[
k
],
v
);
}
}
}
}
/* ----------------------------------------------------------------------
tally energy and virial from 1-3 repulsion of SDK angle into accumulators
------------------------------------------------------------------------- */
void
ThrOMP
::
ev_tally13_thr
(
Angle
*
const
angle
,
const
int
i1
,
const
int
i3
,
const
int
nlocal
,
const
int
newton_bond
,
const
double
epair
,
const
double
fpair
,
const
double
delx
,
const
double
dely
,
const
double
delz
,
ThrData
*
const
thr
)
{
if
(
angle
->
eflag_either
)
{
const
double
epairhalf
=
0.5
*
epair
;
if
(
angle
->
eflag_global
)
{
if
(
newton_bond
||
i1
<
nlocal
)
thr
->
eng_angle
+=
epairhalf
;
if
(
newton_bond
||
i3
<
nlocal
)
thr
->
eng_angle
+=
epairhalf
;
}
if
(
angle
->
eflag_atom
)
{
if
(
newton_bond
||
i1
<
nlocal
)
thr
->
eatom_angle
[
i1
]
+=
epairhalf
;
if
(
newton_bond
||
i3
<
nlocal
)
thr
->
eatom_angle
[
i3
]
+=
epairhalf
;
}
}
if
(
angle
->
vflag_either
)
{
double
v
[
6
];
v
[
0
]
=
delx
*
delx
*
fpair
;
v
[
1
]
=
dely
*
dely
*
fpair
;
v
[
2
]
=
delz
*
delz
*
fpair
;
v
[
3
]
=
delx
*
dely
*
fpair
;
v
[
4
]
=
delx
*
delz
*
fpair
;
v
[
5
]
=
dely
*
delz
*
fpair
;
if
(
angle
->
vflag_global
)
{
double
*
const
va
=
thr
->
virial_angle
;
if
(
newton_bond
||
i1
<
nlocal
)
v_tally
(
va
,
0.5
,
v
);
if
(
newton_bond
||
i3
<
nlocal
)
v_tally
(
va
,
0.5
,
v
);
}
if
(
angle
->
vflag_atom
)
{
if
(
newton_bond
||
i1
<
nlocal
)
{
double
*
const
va
=
thr
->
vatom_angle
[
i1
];
v_tally
(
va
,
0.5
,
v
);
}
if
(
newton_bond
||
i3
<
nlocal
)
{
double
*
const
va
=
thr
->
vatom_angle
[
i3
];
v_tally
(
va
,
0.5
,
v
);
}
}
}
}
/* ----------------------------------------------------------------------
tally energy and virial into global and per-atom accumulators
virial = r1F1 + r2F2 + r3F3 + r4F4 = (r1-r2) F1 + (r3-r2) F3 + (r4-r2) F4
= (r1-r2) F1 + (r3-r2) F3 + (r4-r3 + r3-r2) F4
= vb1*f1 + vb2*f3 + (vb3+vb2)*f4
------------------------------------------------------------------------- */
void
ThrOMP
::
ev_tally_thr
(
Dihedral
*
const
dihed
,
const
int
i1
,
const
int
i2
,
const
int
i3
,
const
int
i4
,
const
int
nlocal
,
const
int
newton_bond
,
const
double
edihedral
,
const
double
*
const
f1
,
const
double
*
const
f3
,
const
double
*
const
f4
,
const
double
vb1x
,
const
double
vb1y
,
const
double
vb1z
,
const
double
vb2x
,
const
double
vb2y
,
const
double
vb2z
,
const
double
vb3x
,
const
double
vb3y
,
const
double
vb3z
,
ThrData
*
const
thr
)
{
if
(
dihed
->
eflag_either
)
{
if
(
dihed
->
eflag_global
)
{
if
(
newton_bond
)
{
thr
->
eng_dihed
+=
edihedral
;
}
else
{
const
double
edihedralquarter
=
0.25
*
edihedral
;
int
cnt
=
0
;
if
(
i1
<
nlocal
)
++
cnt
;
if
(
i2
<
nlocal
)
++
cnt
;
if
(
i3
<
nlocal
)
++
cnt
;
if
(
i4
<
nlocal
)
++
cnt
;
thr
->
eng_dihed
+=
static_cast
<
double
>
(
cnt
)
*
edihedralquarter
;
}
}
if
(
dihed
->
eflag_atom
)
{
const
double
edihedralquarter
=
0.25
*
edihedral
;
if
(
newton_bond
)
{
thr
->
eatom_dihed
[
i1
]
+=
edihedralquarter
;
thr
->
eatom_dihed
[
i2
]
+=
edihedralquarter
;
thr
->
eatom_dihed
[
i3
]
+=
edihedralquarter
;
thr
->
eatom_dihed
[
i4
]
+=
edihedralquarter
;
}
else
{
if
(
i1
<
nlocal
)
thr
->
eatom_dihed
[
i1
]
+=
edihedralquarter
;
if
(
i2
<
nlocal
)
thr
->
eatom_dihed
[
i2
]
+=
edihedralquarter
;
if
(
i3
<
nlocal
)
thr
->
eatom_dihed
[
i3
]
+=
edihedralquarter
;
if
(
i4
<
nlocal
)
thr
->
eatom_dihed
[
i4
]
+=
edihedralquarter
;
}
}
}
if
(
dihed
->
vflag_either
)
{
double
v
[
6
];
v
[
0
]
=
vb1x
*
f1
[
0
]
+
vb2x
*
f3
[
0
]
+
(
vb3x
+
vb2x
)
*
f4
[
0
];
v
[
1
]
=
vb1y
*
f1
[
1
]
+
vb2y
*
f3
[
1
]
+
(
vb3y
+
vb2y
)
*
f4
[
1
];
v
[
2
]
=
vb1z
*
f1
[
2
]
+
vb2z
*
f3
[
2
]
+
(
vb3z
+
vb2z
)
*
f4
[
2
];
v
[
3
]
=
vb1x
*
f1
[
1
]
+
vb2x
*
f3
[
1
]
+
(
vb3x
+
vb2x
)
*
f4
[
1
];
v
[
4
]
=
vb1x
*
f1
[
2
]
+
vb2x
*
f3
[
2
]
+
(
vb3x
+
vb2x
)
*
f4
[
2
];
v
[
5
]
=
vb1y
*
f1
[
2
]
+
vb2y
*
f3
[
2
]
+
(
vb3y
+
vb2y
)
*
f4
[
2
];
if
(
dihed
->
vflag_global
)
{
if
(
newton_bond
)
{
v_tally
(
thr
->
virial_dihed
,
v
);
}
else
{
int
cnt
=
0
;
if
(
i1
<
nlocal
)
++
cnt
;
if
(
i2
<
nlocal
)
++
cnt
;
if
(
i3
<
nlocal
)
++
cnt
;
if
(
i4
<
nlocal
)
++
cnt
;
v_tally
(
thr
->
virial_dihed
,
0.25
*
static_cast
<
double
>
(
cnt
),
v
);
}
}
v
[
0
]
*=
0.25
;
v
[
1
]
*=
0.25
;
v
[
2
]
*=
0.25
;
v
[
3
]
*=
0.25
;
v
[
4
]
*=
0.25
;
v
[
5
]
*=
0.25
;
if
(
dihed
->
vflag_atom
)
{
if
(
newton_bond
)
{
v_tally
(
thr
->
vatom_dihed
[
i1
],
v
);
v_tally
(
thr
->
vatom_dihed
[
i2
],
v
);
v_tally
(
thr
->
vatom_dihed
[
i3
],
v
);
v_tally
(
thr
->
vatom_dihed
[
i4
],
v
);
}
else
{
if
(
i1
<
nlocal
)
v_tally
(
thr
->
vatom_dihed
[
i1
],
v
);
if
(
i2
<
nlocal
)
v_tally
(
thr
->
vatom_dihed
[
i2
],
v
);
if
(
i3
<
nlocal
)
v_tally
(
thr
->
vatom_dihed
[
i3
],
v
);
if
(
i4
<
nlocal
)
v_tally
(
thr
->
vatom_dihed
[
i4
],
v
);
}
}
}
}
/* ----------------------------------------------------------------------
tally energy and virial into global and per-atom accumulators
virial = r1F1 + r2F2 + r3F3 + r4F4 = (r1-r2) F1 + (r3-r2) F3 + (r4-r2) F4
= (r1-r2) F1 + (r3-r2) F3 + (r4-r3 + r3-r2) F4
= vb1*f1 + vb2*f3 + (vb3+vb2)*f4
------------------------------------------------------------------------- */
void
ThrOMP
::
ev_tally_thr
(
Improper
*
const
imprp
,
const
int
i1
,
const
int
i2
,
const
int
i3
,
const
int
i4
,
const
int
nlocal
,
const
int
newton_bond
,
const
double
eimproper
,
const
double
*
const
f1
,
const
double
*
const
f3
,
const
double
*
const
f4
,
const
double
vb1x
,
const
double
vb1y
,
const
double
vb1z
,
const
double
vb2x
,
const
double
vb2y
,
const
double
vb2z
,
const
double
vb3x
,
const
double
vb3y
,
const
double
vb3z
,
ThrData
*
const
thr
)
{
if
(
imprp
->
eflag_either
)
{
if
(
imprp
->
eflag_global
)
{
if
(
newton_bond
)
{
thr
->
eng_imprp
+=
eimproper
;
}
else
{
const
double
eimproperquarter
=
0.25
*
eimproper
;
int
cnt
=
0
;
if
(
i1
<
nlocal
)
++
cnt
;
if
(
i2
<
nlocal
)
++
cnt
;
if
(
i3
<
nlocal
)
++
cnt
;
if
(
i4
<
nlocal
)
++
cnt
;
thr
->
eng_imprp
+=
static_cast
<
double
>
(
cnt
)
*
eimproperquarter
;
}
}
if
(
imprp
->
eflag_atom
)
{
const
double
eimproperquarter
=
0.25
*
eimproper
;
if
(
newton_bond
)
{
thr
->
eatom_imprp
[
i1
]
+=
eimproperquarter
;
thr
->
eatom_imprp
[
i2
]
+=
eimproperquarter
;
thr
->
eatom_imprp
[
i3
]
+=
eimproperquarter
;
thr
->
eatom_imprp
[
i4
]
+=
eimproperquarter
;
}
else
{
if
(
i1
<
nlocal
)
thr
->
eatom_imprp
[
i1
]
+=
eimproperquarter
;
if
(
i2
<
nlocal
)
thr
->
eatom_imprp
[
i2
]
+=
eimproperquarter
;
if
(
i3
<
nlocal
)
thr
->
eatom_imprp
[
i3
]
+=
eimproperquarter
;
if
(
i4
<
nlocal
)
thr
->
eatom_imprp
[
i4
]
+=
eimproperquarter
;
}
}
}
if
(
imprp
->
vflag_either
)
{
double
v
[
6
];
v
[
0
]
=
vb1x
*
f1
[
0
]
+
vb2x
*
f3
[
0
]
+
(
vb3x
+
vb2x
)
*
f4
[
0
];
v
[
1
]
=
vb1y
*
f1
[
1
]
+
vb2y
*
f3
[
1
]
+
(
vb3y
+
vb2y
)
*
f4
[
1
];
v
[
2
]
=
vb1z
*
f1
[
2
]
+
vb2z
*
f3
[
2
]
+
(
vb3z
+
vb2z
)
*
f4
[
2
];
v
[
3
]
=
vb1x
*
f1
[
1
]
+
vb2x
*
f3
[
1
]
+
(
vb3x
+
vb2x
)
*
f4
[
1
];
v
[
4
]
=
vb1x
*
f1
[
2
]
+
vb2x
*
f3
[
2
]
+
(
vb3x
+
vb2x
)
*
f4
[
2
];
v
[
5
]
=
vb1y
*
f1
[
2
]
+
vb2y
*
f3
[
2
]
+
(
vb3y
+
vb2y
)
*
f4
[
2
];
if
(
imprp
->
vflag_global
)
{
if
(
newton_bond
)
{
v_tally
(
thr
->
virial_imprp
,
v
);
}
else
{
int
cnt
=
0
;
if
(
i1
<
nlocal
)
++
cnt
;
if
(
i2
<
nlocal
)
++
cnt
;
if
(
i3
<
nlocal
)
++
cnt
;
if
(
i4
<
nlocal
)
++
cnt
;
v_tally
(
thr
->
virial_imprp
,
0.25
*
static_cast
<
double
>
(
cnt
),
v
);
}
}
v
[
0
]
*=
0.25
;
v
[
1
]
*=
0.25
;
v
[
2
]
*=
0.25
;
v
[
3
]
*=
0.25
;
v
[
4
]
*=
0.25
;
v
[
5
]
*=
0.25
;
if
(
imprp
->
vflag_atom
)
{
if
(
newton_bond
)
{
v_tally
(
thr
->
vatom_imprp
[
i1
],
v
);
v_tally
(
thr
->
vatom_imprp
[
i2
],
v
);
v_tally
(
thr
->
vatom_imprp
[
i3
],
v
);
v_tally
(
thr
->
vatom_imprp
[
i4
],
v
);
}
else
{
if
(
i1
<
nlocal
)
v_tally
(
thr
->
vatom_imprp
[
i1
],
v
);
if
(
i2
<
nlocal
)
v_tally
(
thr
->
vatom_imprp
[
i2
],
v
);
if
(
i3
<
nlocal
)
v_tally
(
thr
->
vatom_imprp
[
i3
],
v
);
if
(
i4
<
nlocal
)
v_tally
(
thr
->
vatom_imprp
[
i4
],
v
);
}
}
}
}
/* ----------------------------------------------------------------------
tally virial into per-atom accumulators
called by AIREBO potential, newton_pair is always on
fpair is magnitude of force on atom I
------------------------------------------------------------------------- */
void
ThrOMP
::
v_tally2_thr
(
const
int
i
,
const
int
j
,
const
double
fpair
,
const
double
*
const
drij
,
ThrData
*
const
thr
)
{
double
v
[
6
];
v
[
0
]
=
0.5
*
drij
[
0
]
*
drij
[
0
]
*
fpair
;
v
[
1
]
=
0.5
*
drij
[
1
]
*
drij
[
1
]
*
fpair
;
v
[
2
]
=
0.5
*
drij
[
2
]
*
drij
[
2
]
*
fpair
;
v
[
3
]
=
0.5
*
drij
[
0
]
*
drij
[
1
]
*
fpair
;
v
[
4
]
=
0.5
*
drij
[
0
]
*
drij
[
2
]
*
fpair
;
v
[
5
]
=
0.5
*
drij
[
1
]
*
drij
[
2
]
*
fpair
;
v_tally
(
thr
->
vatom_pair
[
i
],
v
);
v_tally
(
thr
->
vatom_pair
[
j
],
v
);
}
/* ----------------------------------------------------------------------
tally virial into per-atom accumulators
called by AIREBO and Tersoff potential, newton_pair is always on
------------------------------------------------------------------------- */
void
ThrOMP
::
v_tally3_thr
(
const
int
i
,
const
int
j
,
const
int
k
,
const
double
*
const
fi
,
const
double
*
const
fj
,
const
double
*
const
drik
,
const
double
*
const
drjk
,
ThrData
*
const
thr
)
{
double
v
[
6
];
v
[
0
]
=
THIRD
*
(
drik
[
0
]
*
fi
[
0
]
+
drjk
[
0
]
*
fj
[
0
]);
v
[
1
]
=
THIRD
*
(
drik
[
1
]
*
fi
[
1
]
+
drjk
[
1
]
*
fj
[
1
]);
v
[
2
]
=
THIRD
*
(
drik
[
2
]
*
fi
[
2
]
+
drjk
[
2
]
*
fj
[
2
]);
v
[
3
]
=
THIRD
*
(
drik
[
0
]
*
fi
[
1
]
+
drjk
[
0
]
*
fj
[
1
]);
v
[
4
]
=
THIRD
*
(
drik
[
0
]
*
fi
[
2
]
+
drjk
[
0
]
*
fj
[
2
]);
v
[
5
]
=
THIRD
*
(
drik
[
1
]
*
fi
[
2
]
+
drjk
[
1
]
*
fj
[
2
]);
v_tally
(
thr
->
vatom_pair
[
i
],
v
);
v_tally
(
thr
->
vatom_pair
[
j
],
v
);
v_tally
(
thr
->
vatom_pair
[
k
],
v
);
}
/* ----------------------------------------------------------------------
tally virial into per-atom accumulators
called by AIREBO potential, newton_pair is always on
------------------------------------------------------------------------- */
void
ThrOMP
::
v_tally4_thr
(
const
int
i
,
const
int
j
,
const
int
k
,
const
int
m
,
const
double
*
const
fi
,
const
double
*
const
fj
,
const
double
*
const
fk
,
const
double
*
const
drim
,
const
double
*
const
drjm
,
const
double
*
const
drkm
,
ThrData
*
const
thr
)
{
double
v
[
6
];
v
[
0
]
=
0.25
*
(
drim
[
0
]
*
fi
[
0
]
+
drjm
[
0
]
*
fj
[
0
]
+
drkm
[
0
]
*
fk
[
0
]);
v
[
1
]
=
0.25
*
(
drim
[
1
]
*
fi
[
1
]
+
drjm
[
1
]
*
fj
[
1
]
+
drkm
[
1
]
*
fk
[
1
]);
v
[
2
]
=
0.25
*
(
drim
[
2
]
*
fi
[
2
]
+
drjm
[
2
]
*
fj
[
2
]
+
drkm
[
2
]
*
fk
[
2
]);
v
[
3
]
=
0.25
*
(
drim
[
0
]
*
fi
[
1
]
+
drjm
[
0
]
*
fj
[
1
]
+
drkm
[
0
]
*
fk
[
1
]);
v
[
4
]
=
0.25
*
(
drim
[
0
]
*
fi
[
2
]
+
drjm
[
0
]
*
fj
[
2
]
+
drkm
[
0
]
*
fk
[
2
]);
v
[
5
]
=
0.25
*
(
drim
[
1
]
*
fi
[
2
]
+
drjm
[
1
]
*
fj
[
2
]
+
drkm
[
1
]
*
fk
[
2
]);
v_tally
(
thr
->
vatom_pair
[
i
],
v
);
v_tally
(
thr
->
vatom_pair
[
j
],
v
);
v_tally
(
thr
->
vatom_pair
[
k
],
v
);
v_tally
(
thr
->
vatom_pair
[
m
],
v
);
}
/* ---------------------------------------------------------------------- */
double
ThrOMP
::
memory_usage_thr
()
{
double
bytes
=
0.0
;
return
bytes
;
}
Event Timeline
Log In to Comment