Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F80613515
reaxc_valence_angles_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
Sun, Sep 1, 16:24
Size
23 KB
Mime Type
text/x-c
Expires
Tue, Sep 3, 16:24 (2 d)
Engine
blob
Format
Raw Data
Handle
20404240
Attached To
rLAMMPS lammps
reaxc_valence_angles_omp.cpp
View Options
/*----------------------------------------------------------------------
PuReMD - Purdue ReaxFF Molecular Dynamics Program
Website: https://www.cs.purdue.edu/puremd
Copyright (2010) Purdue University
Contributing authors:
H. M. Aktulga, J. Fogarty, S. Pandit, A. Grama
Corresponding author:
Hasan Metin Aktulga, Michigan State University, hma@cse.msu.edu
Please cite the related publication:
H. M. Aktulga, J. C. Fogarty, S. A. Pandit, A. Y. Grama,
"Parallel Reactive Molecular Dynamics: Numerical Methods and
Algorithmic Techniques", Parallel Computing, 38 (4-5), 245-259
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License as
published by the Free Software Foundation; either version 2 of
the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU General Public License for more details:
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "pair_reaxc_omp.h"
#include "thr_data.h"
#include "reaxc_types.h"
#include "reaxc_valence_angles.h"
#include "reaxc_valence_angles_omp.h"
#include "reaxc_bond_orders_omp.h"
#include "reaxc_list.h"
#include "reaxc_vector.h"
#if defined(_OPENMP)
#include <omp.h>
#endif
using
namespace
LAMMPS_NS
;
/* ---------------------------------------------------------------------- */
void
Calculate_dCos_ThetaOMP
(
rvec
dvec_ji
,
double
d_ji
,
rvec
dvec_jk
,
double
d_jk
,
rvec
*
dcos_theta_di
,
rvec
*
dcos_theta_dj
,
rvec
*
dcos_theta_dk
)
{
double
sqr_d_ji
=
SQR
(
d_ji
);
double
sqr_d_jk
=
SQR
(
d_jk
);
double
inv_dists
=
1.0
/
(
d_ji
*
d_jk
);
double
inv_dists3
=
inv_dists
*
inv_dists
*
inv_dists
;
double
dot_dvecs
=
dvec_ji
[
0
]
*
dvec_jk
[
0
]
+
dvec_ji
[
1
]
*
dvec_jk
[
1
]
+
dvec_ji
[
2
]
*
dvec_jk
[
2
];
double
Cdot_inv3
=
dot_dvecs
*
inv_dists3
;
double
csqr_jk
=
Cdot_inv3
*
sqr_d_jk
;
double
csqr_ji
=
Cdot_inv3
*
sqr_d_ji
;
// Try to help compiler out by unrolling
// x-component
double
dinv_jk
=
dvec_jk
[
0
]
*
inv_dists
;
double
dinv_ji
=
dvec_ji
[
0
]
*
inv_dists
;
double
cdev_ji
=
csqr_jk
*
dvec_ji
[
0
];
double
cdev_jk
=
csqr_ji
*
dvec_jk
[
0
];
(
*
dcos_theta_di
)[
0
]
=
dinv_jk
-
cdev_ji
;
(
*
dcos_theta_dj
)[
0
]
=
-
(
dinv_jk
+
dinv_ji
)
+
cdev_ji
+
cdev_jk
;
(
*
dcos_theta_dk
)[
0
]
=
dinv_ji
-
cdev_jk
;
// y-component
dinv_jk
=
dvec_jk
[
1
]
*
inv_dists
;
dinv_ji
=
dvec_ji
[
1
]
*
inv_dists
;
cdev_ji
=
csqr_jk
*
dvec_ji
[
1
];
cdev_jk
=
csqr_ji
*
dvec_jk
[
1
];
(
*
dcos_theta_di
)[
1
]
=
dinv_jk
-
cdev_ji
;
(
*
dcos_theta_dj
)[
1
]
=
-
(
dinv_jk
+
dinv_ji
)
+
cdev_ji
+
cdev_jk
;
(
*
dcos_theta_dk
)[
1
]
=
dinv_ji
-
cdev_jk
;
// z-component
dinv_jk
=
dvec_jk
[
2
]
*
inv_dists
;
dinv_ji
=
dvec_ji
[
2
]
*
inv_dists
;
cdev_ji
=
csqr_jk
*
dvec_ji
[
2
];
cdev_jk
=
csqr_ji
*
dvec_jk
[
2
];
(
*
dcos_theta_di
)[
2
]
=
dinv_jk
-
cdev_ji
;
(
*
dcos_theta_dj
)[
2
]
=
-
(
dinv_jk
+
dinv_ji
)
+
cdev_ji
+
cdev_jk
;
(
*
dcos_theta_dk
)[
2
]
=
dinv_ji
-
cdev_jk
;
}
/* ---------------------------------------------------------------------- */
/* this is a 3-body interaction in which the main role is
played by j which sits in the middle of the other two. */
void
Valence_AnglesOMP
(
reax_system
*
system
,
control_params
*
control
,
simulation_data
*
data
,
storage
*
workspace
,
reax_list
**
lists
,
output_controls
*
out_control
)
{
#ifdef OMP_TIMING
double
endTimeBase
,
startTimeBase
;
startTimeBase
=
MPI_Wtime
();
#endif
reax_list
*
bonds
=
(
*
lists
)
+
BONDS
;
reax_list
*
thb_intrs
=
(
*
lists
)
+
THREE_BODIES
;
// Precompute and store valence_angle offsets for OpenMP code.
int
*
_my_offset
=
workspace
->
valence_angle_atom_myoffset
;
/* global parameters used in these calculations */
double
p_val6
=
system
->
reax_param
.
gp
.
l
[
14
];
double
p_val8
=
system
->
reax_param
.
gp
.
l
[
33
];
double
p_val9
=
system
->
reax_param
.
gp
.
l
[
16
];
double
p_val10
=
system
->
reax_param
.
gp
.
l
[
17
];
double
total_Eang
=
0
;
double
total_Epen
=
0
;
double
total_Ecoa
=
0
;
int
nthreads
=
control
->
nthreads
;
int
num_thb_intrs
=
0
;
int
TWICE
=
2
;
#if defined(_OPENMP)
#pragma omp parallel default(shared) reduction(+:total_Eang, total_Epen, total_Ecoa, num_thb_intrs)
#endif
{
int
i
,
j
,
pi
,
k
,
pk
,
t
;
int
type_i
,
type_j
,
type_k
;
int
start_j
,
end_j
,
start_pk
,
end_pk
;
int
cnt
,
my_offset
;
double
temp
,
temp_bo_jt
,
pBOjt7
;
double
p_val1
,
p_val2
,
p_val3
,
p_val4
,
p_val5
,
p_val7
;
double
p_pen1
,
p_pen2
,
p_pen3
,
p_pen4
;
double
p_coa1
,
p_coa2
,
p_coa3
,
p_coa4
;
double
trm8
,
expval6
,
expval7
,
expval2theta
,
expval12theta
,
exp3ij
,
exp3jk
;
double
exp_pen2ij
,
exp_pen2jk
,
exp_pen3
,
exp_pen4
,
trm_pen34
,
exp_coa2
;
double
dSBO1
,
dSBO2
,
SBO
,
SBO2
,
CSBO2
,
SBOp
,
prod_SBO
,
vlpadj
;
double
CEval1
,
CEval2
,
CEval3
,
CEval4
,
CEval5
,
CEval6
,
CEval7
,
CEval8
;
double
CEpen1
,
CEpen2
,
CEpen3
;
double
e_ang
,
e_coa
,
e_pen
;
double
CEcoa1
,
CEcoa2
,
CEcoa3
,
CEcoa4
,
CEcoa5
;
double
Cf7ij
,
Cf7jk
,
Cf8j
,
Cf9j
;
double
f7_ij
,
f7_jk
,
f8_Dj
,
f9_Dj
;
double
Ctheta_0
,
theta_0
,
theta_00
,
theta
,
cos_theta
,
sin_theta
;
double
BOA_ij
,
BOA_jk
;
rvec
force
,
ext_press
;
// rtensor temp_rtensor, total_rtensor;
// Tallying variables
double
eng_tmp
,
fi_tmp
[
3
],
fj_tmp
[
3
],
fk_tmp
[
3
];
double
delij
[
3
],
delkj
[
3
];
three_body_header
*
thbh
;
three_body_parameters
*
thbp
;
three_body_interaction_data
*
p_ijk
,
*
p_kji
;
bond_data
*
pbond_ij
,
*
pbond_jk
,
*
pbond_jt
;
bond_order_data
*
bo_ij
,
*
bo_jk
,
*
bo_jt
;
#if defined(_OPENMP)
int
tid
=
omp_get_thread_num
();
#else
int
tid
=
0
;
#endif
long
reductionOffset
=
(
system
->
N
*
tid
);
class
PairReaxCOMP
*
pair_reax_ptr
;
pair_reax_ptr
=
static_cast
<
class
PairReaxCOMP
*>
(
system
->
pair_ptr
);
class
ThrData
*
thr
=
pair_reax_ptr
->
getFixOMP
()
->
get_thr
(
tid
);
pair_reax_ptr
->
ev_setup_thr_proxy
(
system
->
pair_ptr
->
eflag_either
,
system
->
pair_ptr
->
vflag_either
,
system
->
N
,
system
->
pair_ptr
->
eatom
,
system
->
pair_ptr
->
vatom
,
thr
);
// Run through a minimal for(j<N) loop once to precompute offsets with safe number of threads
const
int
per_thread
=
thb_intrs
->
num_intrs
/
nthreads
;
#if defined(_OPENMP)
#pragma omp for schedule(dynamic,50)
#endif
for
(
j
=
0
;
j
<
system
->
N
;
++
j
)
{
type_j
=
system
->
my_atoms
[
j
].
type
;
_my_offset
[
j
]
=
0
;
if
(
type_j
<
0
)
continue
;
start_j
=
Start_Index
(
j
,
bonds
);
end_j
=
End_Index
(
j
,
bonds
);
// Always point to start of workspace to count angles
my_offset
=
tid
*
per_thread
;
for
(
pi
=
start_j
;
pi
<
end_j
;
++
pi
)
{
Set_Start_Index
(
pi
,
my_offset
,
thb_intrs
);
pbond_ij
=
&
(
bonds
->
select
.
bond_list
[
pi
]);
bo_ij
=
&
(
pbond_ij
->
bo_data
);
BOA_ij
=
bo_ij
->
BO
-
control
->
thb_cut
;
if
(
BOA_ij
>
0.0
)
{
i
=
pbond_ij
->
nbr
;
/* first copy 3-body intrs from previously computed ones where i>k.
in the second for-loop below,
we compute only new 3-body intrs where i < k */
for
(
pk
=
start_j
;
pk
<
pi
;
++
pk
)
{
start_pk
=
Start_Index
(
pk
,
thb_intrs
);
end_pk
=
End_Index
(
pk
,
thb_intrs
);
for
(
t
=
start_pk
;
t
<
end_pk
;
++
t
)
if
(
thb_intrs
->
select
.
three_body_list
[
t
].
thb
==
i
)
{
p_ijk
=
&
(
thb_intrs
->
select
.
three_body_list
[
my_offset
]
);
p_ijk
->
thb
=
bonds
->
select
.
bond_list
[
pk
].
nbr
;
++
my_offset
;
break
;
}
}
// for(pk)
/* and this is the second for loop mentioned above */
for
(
pk
=
pi
+
1
;
pk
<
end_j
;
++
pk
)
{
pbond_jk
=
&
(
bonds
->
select
.
bond_list
[
pk
]);
k
=
pbond_jk
->
nbr
;
if
(
j
>=
system
->
n
&&
i
>=
system
->
n
&&
k
>=
system
->
n
)
continue
;
p_ijk
=
&
(
thb_intrs
->
select
.
three_body_list
[
my_offset
]
);
p_ijk
->
thb
=
k
;
++
my_offset
;
// add this to the list of 3-body interactions
}
// for(pk)
}
// if()
Set_End_Index
(
pi
,
my_offset
,
thb_intrs
);
}
// for(pi)
// Confirm that thb_intrs->num_intrs / nthreads is enough to hold all angles from a single atom
if
(
my_offset
>=
(
tid
+
1
)
*
per_thread
)
{
int
me
;
MPI_Comm_rank
(
MPI_COMM_WORLD
,
&
me
);
fprintf
(
stderr
,
"step%d-ran out of space on angle_list on proc %i for atom %i:"
,
data
->
step
,
me
,
j
);
fprintf
(
stderr
,
" nthreads= %d, tid=%d, my_offset=%d, per_thread=%d
\n
"
,
nthreads
,
tid
,
my_offset
,
per_thread
);
fprintf
(
stderr
,
" num_intrs= %i N= %i
\n
"
,
thb_intrs
->
num_intrs
,
system
->
N
);
MPI_Abort
(
MPI_COMM_WORLD
,
INSUFFICIENT_MEMORY
);
}
// Number of angles owned by this atom
_my_offset
[
j
]
=
my_offset
-
tid
*
per_thread
;
}
// for(j)
// Wait for all threads to finish counting angles
#if defined(_OPENMP) && !defined(__NVCC__)
#pragma omp barrier
#endif
// Master thread uses angle counts to compute offsets
// This can be threaded
#if defined(_OPENMP) && !defined(__NVCC__)
#pragma omp master
#endif
{
int
current_count
=
0
;
int
m
=
_my_offset
[
0
];
_my_offset
[
0
]
=
current_count
;
for
(
j
=
1
;
j
<
system
->
N
;
j
++
)
{
current_count
+=
m
;
m
=
_my_offset
[
j
];
_my_offset
[
j
]
=
current_count
;
}
_my_offset
[
system
->
N
]
=
current_count
+
m
;
// Used to test if last particle has any angles
}
// All threads wait till master thread finished computing offsets
#if defined(_OPENMP) && !defined(__NVCC__)
#pragma omp barrier
#endif
// Original loop, but now using precomputed offsets
// Safe to use all threads available, regardless of threads tasked above
// We also now skip over atoms that have no angles assigned
#if defined(_OPENMP)
#pragma omp for schedule(dynamic,50)
//(dynamic,chunksize)//(guided)
#endif
for
(
j
=
0
;
j
<
system
->
N
;
++
j
)
{
// Ray: the first one with system->N
type_j
=
system
->
my_atoms
[
j
].
type
;
if
(
type_j
<
0
)
continue
;
// Skip if no angles for this atom
if
(
_my_offset
[
j
]
==
_my_offset
[
j
+
1
])
continue
;
start_j
=
Start_Index
(
j
,
bonds
);
end_j
=
End_Index
(
j
,
bonds
);
type_j
=
system
->
my_atoms
[
j
].
type
;
my_offset
=
_my_offset
[
j
];
p_val3
=
system
->
reax_param
.
sbp
[
type_j
].
p_val3
;
p_val5
=
system
->
reax_param
.
sbp
[
type_j
].
p_val5
;
SBOp
=
0
,
prod_SBO
=
1
;
for
(
t
=
start_j
;
t
<
end_j
;
++
t
)
{
bo_jt
=
&
(
bonds
->
select
.
bond_list
[
t
].
bo_data
);
SBOp
+=
(
bo_jt
->
BO_pi
+
bo_jt
->
BO_pi2
);
temp
=
SQR
(
bo_jt
->
BO
);
temp
*=
temp
;
temp
*=
temp
;
prod_SBO
*=
exp
(
-
temp
);
}
// modifications to match Adri's code - 09/01/09
if
(
workspace
->
vlpex
[
j
]
>=
0
){
vlpadj
=
0
;
dSBO2
=
prod_SBO
-
1
;
}
else
{
vlpadj
=
workspace
->
nlp
[
j
];
dSBO2
=
(
prod_SBO
-
1
)
*
(
1
-
p_val8
*
workspace
->
dDelta_lp
[
j
]);
}
SBO
=
SBOp
+
(
1
-
prod_SBO
)
*
(
-
workspace
->
Delta_boc
[
j
]
-
p_val8
*
vlpadj
);
dSBO1
=
-
8
*
prod_SBO
*
(
workspace
->
Delta_boc
[
j
]
+
p_val8
*
vlpadj
);
if
(
SBO
<=
0
)
SBO2
=
0
,
CSBO2
=
0
;
else
if
(
SBO
>
0
&&
SBO
<=
1
)
{
SBO2
=
pow
(
SBO
,
p_val9
);
CSBO2
=
p_val9
*
pow
(
SBO
,
p_val9
-
1
);
}
else
if
(
SBO
>
1
&&
SBO
<
2
)
{
SBO2
=
2
-
pow
(
2
-
SBO
,
p_val9
);
CSBO2
=
p_val9
*
pow
(
2
-
SBO
,
p_val9
-
1
);
}
else
SBO2
=
2
,
CSBO2
=
0
;
expval6
=
exp
(
p_val6
*
workspace
->
Delta_boc
[
j
]
);
for
(
pi
=
start_j
;
pi
<
end_j
;
++
pi
)
{
Set_Start_Index
(
pi
,
my_offset
,
thb_intrs
);
pbond_ij
=
&
(
bonds
->
select
.
bond_list
[
pi
]);
bo_ij
=
&
(
pbond_ij
->
bo_data
);
BOA_ij
=
bo_ij
->
BO
-
control
->
thb_cut
;
if
(
BOA_ij
>
0.0
)
{
i
=
pbond_ij
->
nbr
;
type_i
=
system
->
my_atoms
[
i
].
type
;
/* first copy 3-body intrs from previously computed ones where i>k.
in the second for-loop below,
we compute only new 3-body intrs where i < k */
for
(
pk
=
start_j
;
pk
<
pi
;
++
pk
)
{
start_pk
=
Start_Index
(
pk
,
thb_intrs
);
end_pk
=
End_Index
(
pk
,
thb_intrs
);
for
(
t
=
start_pk
;
t
<
end_pk
;
++
t
)
if
(
thb_intrs
->
select
.
three_body_list
[
t
].
thb
==
i
)
{
p_ijk
=
&
(
thb_intrs
->
select
.
three_body_list
[
my_offset
]
);
p_kji
=
&
(
thb_intrs
->
select
.
three_body_list
[
t
]);
p_ijk
->
thb
=
bonds
->
select
.
bond_list
[
pk
].
nbr
;
p_ijk
->
pthb
=
pk
;
p_ijk
->
theta
=
p_kji
->
theta
;
rvec_Copy
(
p_ijk
->
dcos_di
,
p_kji
->
dcos_dk
);
rvec_Copy
(
p_ijk
->
dcos_dj
,
p_kji
->
dcos_dj
);
rvec_Copy
(
p_ijk
->
dcos_dk
,
p_kji
->
dcos_di
);
++
my_offset
;
++
num_thb_intrs
;
break
;
}
}
// for(pk)
/* and this is the second for loop mentioned above */
for
(
pk
=
pi
+
1
;
pk
<
end_j
;
++
pk
)
{
pbond_jk
=
&
(
bonds
->
select
.
bond_list
[
pk
]);
bo_jk
=
&
(
pbond_jk
->
bo_data
);
BOA_jk
=
bo_jk
->
BO
-
control
->
thb_cut
;
k
=
pbond_jk
->
nbr
;
type_k
=
system
->
my_atoms
[
k
].
type
;
p_ijk
=
&
(
thb_intrs
->
select
.
three_body_list
[
my_offset
]
);
// Fix by Sudhir
// if (BOA_jk <= 0) continue;
if
(
j
>=
system
->
n
&&
i
>=
system
->
n
&&
k
>=
system
->
n
)
continue
;
Calculate_Theta
(
pbond_ij
->
dvec
,
pbond_ij
->
d
,
pbond_jk
->
dvec
,
pbond_jk
->
d
,
&
theta
,
&
cos_theta
);
Calculate_dCos_ThetaOMP
(
pbond_ij
->
dvec
,
pbond_ij
->
d
,
pbond_jk
->
dvec
,
pbond_jk
->
d
,
&
(
p_ijk
->
dcos_di
),
&
(
p_ijk
->
dcos_dj
),
&
(
p_ijk
->
dcos_dk
)
);
p_ijk
->
thb
=
k
;
p_ijk
->
pthb
=
pk
;
p_ijk
->
theta
=
theta
;
sin_theta
=
sin
(
theta
);
if
(
sin_theta
<
1.0e-5
)
sin_theta
=
1.0e-5
;
++
my_offset
;
// add this to the list of 3-body interactions
++
num_thb_intrs
;
if
((
j
<
system
->
n
)
&&
(
BOA_jk
>
0.0
)
&&
(
bo_ij
->
BO
>
control
->
thb_cut
)
&&
(
bo_jk
->
BO
>
control
->
thb_cut
)
&&
(
bo_ij
->
BO
*
bo_jk
->
BO
>
control
->
thb_cutsq
))
{
thbh
=
&
(
system
->
reax_param
.
thbp
[
type_i
][
type_j
][
type_k
]
);
for
(
cnt
=
0
;
cnt
<
thbh
->
cnt
;
++
cnt
)
{
if
(
fabs
(
thbh
->
prm
[
cnt
].
p_val1
)
>
0.001
)
{
thbp
=
&
(
thbh
->
prm
[
cnt
]
);
/* ANGLE ENERGY */
p_val1
=
thbp
->
p_val1
;
p_val2
=
thbp
->
p_val2
;
p_val4
=
thbp
->
p_val4
;
p_val7
=
thbp
->
p_val7
;
theta_00
=
thbp
->
theta_00
;
exp3ij
=
exp
(
-
p_val3
*
pow
(
BOA_ij
,
p_val4
)
);
f7_ij
=
1.0
-
exp3ij
;
Cf7ij
=
p_val3
*
p_val4
*
pow
(
BOA_ij
,
p_val4
-
1.0
)
*
exp3ij
;
exp3jk
=
exp
(
-
p_val3
*
pow
(
BOA_jk
,
p_val4
)
);
f7_jk
=
1.0
-
exp3jk
;
Cf7jk
=
p_val3
*
p_val4
*
pow
(
BOA_jk
,
p_val4
-
1.0
)
*
exp3jk
;
expval7
=
exp
(
-
p_val7
*
workspace
->
Delta_boc
[
j
]
);
trm8
=
1.0
+
expval6
+
expval7
;
f8_Dj
=
p_val5
-
(
(
p_val5
-
1.0
)
*
(
2.0
+
expval6
)
/
trm8
);
Cf8j
=
(
(
1.0
-
p_val5
)
/
SQR
(
trm8
)
)
*
(
p_val6
*
expval6
*
trm8
-
(
2.0
+
expval6
)
*
(
p_val6
*
expval6
-
p_val7
*
expval7
)
);
theta_0
=
180.0
-
theta_00
*
(
1.0
-
exp
(
-
p_val10
*
(
2.0
-
SBO2
)));
theta_0
=
DEG2RAD
(
theta_0
);
expval2theta
=
exp
(
-
p_val2
*
SQR
(
theta_0
-
theta
)
);
if
(
p_val1
>=
0
)
expval12theta
=
p_val1
*
(
1.0
-
expval2theta
);
else
// To avoid linear Me-H-Me angles (6/6/06)
expval12theta
=
p_val1
*
-
expval2theta
;
CEval1
=
Cf7ij
*
f7_jk
*
f8_Dj
*
expval12theta
;
CEval2
=
Cf7jk
*
f7_ij
*
f8_Dj
*
expval12theta
;
CEval3
=
Cf8j
*
f7_ij
*
f7_jk
*
expval12theta
;
CEval4
=
-
2.0
*
p_val1
*
p_val2
*
f7_ij
*
f7_jk
*
f8_Dj
*
expval2theta
*
(
theta_0
-
theta
);
Ctheta_0
=
p_val10
*
DEG2RAD
(
theta_00
)
*
exp
(
-
p_val10
*
(
2.0
-
SBO2
)
);
CEval5
=
-
CEval4
*
Ctheta_0
*
CSBO2
;
CEval6
=
CEval5
*
dSBO1
;
CEval7
=
CEval5
*
dSBO2
;
CEval8
=
-
CEval4
/
sin_theta
;
total_Eang
+=
e_ang
=
f7_ij
*
f7_jk
*
f8_Dj
*
expval12theta
;
/* END ANGLE ENERGY*/
/* PENALTY ENERGY */
p_pen1
=
thbp
->
p_pen1
;
p_pen2
=
system
->
reax_param
.
gp
.
l
[
19
];
p_pen3
=
system
->
reax_param
.
gp
.
l
[
20
];
p_pen4
=
system
->
reax_param
.
gp
.
l
[
21
];
exp_pen2ij
=
exp
(
-
p_pen2
*
SQR
(
BOA_ij
-
2.0
)
);
exp_pen2jk
=
exp
(
-
p_pen2
*
SQR
(
BOA_jk
-
2.0
)
);
exp_pen3
=
exp
(
-
p_pen3
*
workspace
->
Delta
[
j
]
);
exp_pen4
=
exp
(
p_pen4
*
workspace
->
Delta
[
j
]
);
trm_pen34
=
1.0
+
exp_pen3
+
exp_pen4
;
f9_Dj
=
(
2.0
+
exp_pen3
)
/
trm_pen34
;
Cf9j
=
(
-
p_pen3
*
exp_pen3
*
trm_pen34
-
(
2.0
+
exp_pen3
)
*
(
-
p_pen3
*
exp_pen3
+
p_pen4
*
exp_pen4
)
)
/
SQR
(
trm_pen34
);
total_Epen
+=
e_pen
=
p_pen1
*
f9_Dj
*
exp_pen2ij
*
exp_pen2jk
;
CEpen1
=
e_pen
*
Cf9j
/
f9_Dj
;
temp
=
-
2.0
*
p_pen2
*
e_pen
;
CEpen2
=
temp
*
(
BOA_ij
-
2.0
);
CEpen3
=
temp
*
(
BOA_jk
-
2.0
);
/* END PENALTY ENERGY */
/* COALITION ENERGY */
p_coa1
=
thbp
->
p_coa1
;
p_coa2
=
system
->
reax_param
.
gp
.
l
[
2
];
p_coa3
=
system
->
reax_param
.
gp
.
l
[
38
];
p_coa4
=
system
->
reax_param
.
gp
.
l
[
30
];
exp_coa2
=
exp
(
p_coa2
*
workspace
->
Delta_val
[
j
]
);
total_Ecoa
+=
e_coa
=
p_coa1
/
(
1.
+
exp_coa2
)
*
exp
(
-
p_coa3
*
SQR
(
workspace
->
total_bond_order
[
i
]
-
BOA_ij
)
)
*
exp
(
-
p_coa3
*
SQR
(
workspace
->
total_bond_order
[
k
]
-
BOA_jk
)
)
*
exp
(
-
p_coa4
*
SQR
(
BOA_ij
-
1.5
)
)
*
exp
(
-
p_coa4
*
SQR
(
BOA_jk
-
1.5
)
);
CEcoa1
=
-
2
*
p_coa4
*
(
BOA_ij
-
1.5
)
*
e_coa
;
CEcoa2
=
-
2
*
p_coa4
*
(
BOA_jk
-
1.5
)
*
e_coa
;
CEcoa3
=
-
p_coa2
*
exp_coa2
*
e_coa
/
(
1
+
exp_coa2
);
CEcoa4
=
-
2
*
p_coa3
*
(
workspace
->
total_bond_order
[
i
]
-
BOA_ij
)
*
e_coa
;
CEcoa5
=
-
2
*
p_coa3
*
(
workspace
->
total_bond_order
[
k
]
-
BOA_jk
)
*
e_coa
;
/* END COALITION ENERGY */
/* FORCES */
bo_ij
->
Cdbo
+=
(
CEval1
+
CEpen2
+
(
CEcoa1
-
CEcoa4
));
bo_jk
->
Cdbo
+=
(
CEval2
+
CEpen3
+
(
CEcoa2
-
CEcoa5
));
workspace
->
CdDelta
[
j
]
+=
((
CEval3
+
CEval7
)
+
CEpen1
+
CEcoa3
);
workspace
->
CdDeltaReduction
[
reductionOffset
+
i
]
+=
CEcoa4
;
workspace
->
CdDeltaReduction
[
reductionOffset
+
k
]
+=
CEcoa5
;
for
(
t
=
start_j
;
t
<
end_j
;
++
t
)
{
pbond_jt
=
&
(
bonds
->
select
.
bond_list
[
t
]
);
bo_jt
=
&
(
pbond_jt
->
bo_data
);
temp_bo_jt
=
bo_jt
->
BO
;
temp
=
CUBE
(
temp_bo_jt
);
pBOjt7
=
temp
*
temp
*
temp_bo_jt
;
bo_jt
->
Cdbo
+=
(
CEval6
*
pBOjt7
);
bo_jt
->
Cdbopi
+=
CEval5
;
bo_jt
->
Cdbopi2
+=
CEval5
;
}
if
(
control
->
virial
==
0
)
{
rvec_ScaledAdd
(
workspace
->
f
[
j
],
CEval8
,
p_ijk
->
dcos_dj
);
rvec_ScaledAdd
(
workspace
->
forceReduction
[
reductionOffset
+
i
],
CEval8
,
p_ijk
->
dcos_di
);
rvec_ScaledAdd
(
workspace
->
forceReduction
[
reductionOffset
+
k
],
CEval8
,
p_ijk
->
dcos_dk
);
}
else
{
/* terms not related to bond order derivatives are
added directly into forces and pressure vector/tensor */
rvec_Scale
(
force
,
CEval8
,
p_ijk
->
dcos_di
);
rvec_Add
(
workspace
->
forceReduction
[
reductionOffset
+
i
],
force
);
rvec_iMultiply
(
ext_press
,
pbond_ij
->
rel_box
,
force
);
rvec_Add
(
workspace
->
my_ext_pressReduction
[
tid
],
ext_press
);
rvec_ScaledAdd
(
workspace
->
f
[
j
],
CEval8
,
p_ijk
->
dcos_dj
);
rvec_Scale
(
force
,
CEval8
,
p_ijk
->
dcos_dk
);
rvec_Add
(
workspace
->
forceReduction
[
reductionOffset
+
k
],
force
);
rvec_iMultiply
(
ext_press
,
pbond_jk
->
rel_box
,
force
);
rvec_Add
(
workspace
->
my_ext_pressReduction
[
tid
],
ext_press
);
}
/* tally into per-atom virials */
if
(
system
->
pair_ptr
->
vflag_atom
||
system
->
pair_ptr
->
evflag
)
{
/* Acquire vectors */
rvec_ScaledSum
(
delij
,
1.
,
system
->
my_atoms
[
i
].
x
,
-
1.
,
system
->
my_atoms
[
j
].
x
);
rvec_ScaledSum
(
delkj
,
1.
,
system
->
my_atoms
[
k
].
x
,
-
1.
,
system
->
my_atoms
[
j
].
x
);
rvec_Scale
(
fi_tmp
,
-
CEval8
,
p_ijk
->
dcos_di
);
rvec_Scale
(
fj_tmp
,
-
CEval8
,
p_ijk
->
dcos_dj
);
rvec_Scale
(
fk_tmp
,
-
CEval8
,
p_ijk
->
dcos_dk
);
eng_tmp
=
e_ang
+
e_pen
+
e_coa
;
if
(
system
->
pair_ptr
->
evflag
)
pair_reax_ptr
->
ev_tally_thr_proxy
(
system
->
pair_ptr
,
j
,
j
,
system
->
N
,
1
,
eng_tmp
,
0.0
,
0.0
,
0.0
,
0.0
,
0.0
,
thr
);
if
(
system
->
pair_ptr
->
vflag_atom
)
// NEED TO MAKE AN OMP VERSION OF THIS CALL!
system
->
pair_ptr
->
v_tally3
(
i
,
j
,
k
,
fi_tmp
,
fk_tmp
,
delij
,
delkj
);
}
}
// if(p_val1>0.001)
}
// for(cnt)
}
// if(j<n && BOA_jk>0)
}
// for(pk)
}
// if(BOA_ij>0)
Set_End_Index
(
pi
,
my_offset
,
thb_intrs
);
}
// for(pi)
}
// for(j)
pair_reax_ptr
->
reduce_thr_proxy
(
system
->
pair_ptr
,
system
->
pair_ptr
->
eflag_either
,
system
->
pair_ptr
->
vflag_either
,
thr
);
}
// end omp parallel
data
->
my_en
.
e_ang
=
total_Eang
;
data
->
my_en
.
e_pen
=
total_Epen
;
data
->
my_en
.
e_coa
=
total_Ecoa
;
if
(
num_thb_intrs
>=
thb_intrs
->
num_intrs
*
DANGER_ZONE
)
{
workspace
->
realloc
.
num_3body
=
num_thb_intrs
*
TWICE
;
if
(
num_thb_intrs
>
thb_intrs
->
num_intrs
)
{
fprintf
(
stderr
,
"step%d-ran out of space on angle_list: top=%d, max=%d"
,
data
->
step
,
num_thb_intrs
,
thb_intrs
->
num_intrs
);
MPI_Abort
(
MPI_COMM_WORLD
,
INSUFFICIENT_MEMORY
);
}
}
#ifdef OMP_TIMING
endTimeBase
=
MPI_Wtime
();
ompTimingData
[
COMPUTEVALENCEANGLESBOINDEX
]
+=
(
endTimeBase
-
startTimeBase
);
#endif
}
Event Timeline
Log In to Comment