Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F78853630
pair_mgpt.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
Fri, Aug 23, 08:38
Size
80 KB
Mime Type
text/x-c
Expires
Sun, Aug 25, 08:38 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20108075
Attached To
rLAMMPS lammps
pair_mgpt.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 authors: Tomas Oppelstrup, LLNL (oppelstrup2@llnl.gov)
and John Moriarty, LLNL (moriarty2@llnl.gov)
Fast MGPT algorithm developed by Tomas Oppelstrup (2015) based on the
matrix MGPT v4.4 FORTRAN routine of John Moriarty (2006) as converted
to C++ for LAMMPS application by Jamie Marian and Alexander Stukowski
(2011). See LLNL copyright notice at bottom of this file.
------------------------------------------------------------------------- */
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cassert>
#include "pair_mgpt.h"
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "memory.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "memory.h"
#include "error.h"
using
namespace
LAMMPS_NS
;
//#define TIMING_ON
#ifdef TIMING_ON
#include <sys/time.h>
#include <time.h>
//#include "rdtsc.h"
#ifdef __bgq__
#include <hwi/include/bqc/A2_inlines.h>
#endif
static
double
gettime
(
int
x
=
0
)
{
if
(
1
)
{
/*
struct timeval tv;
gettimeofday(&tv,NULL);
return tv.tv_sec + 1e-6 * tv.tv_usec;
*/
/*
const double x = 1.0 / CLOCKS_PER_SEC;
return clock() * x;
*/
//const double invfreq = 1.0 / 2394.108e6;
/*
const double invfreq = 1.0 / 700e6;
unsigned long long int x = rdtsc();
return x*invfreq;
*/
const
double
invfreq
=
1.0
/
1.6e9
;
unsigned
long
long
int
x
=
GetTimeBase
();
return
x
*
invfreq
;
}
else
return
0.0
;
}
#else
static
double
gettime
(
int
x
=
0
)
{
return
0.0
;
}
#endif
/* ---------------------------------------------------------------------- */
PairMGPT
::
PairMGPT
(
LAMMPS
*
lmp
)
:
Pair
(
lmp
)
{
single_enable
=
0
;
one_coeff
=
1
;
ghostneigh
=
1
;
}
PairMGPT
::~
PairMGPT
()
{
if
(
allocated
)
{
memory
->
destroy
(
setflag
);
memory
->
destroy
(
cutsq
);
memory
->
destroy
(
cutghost
);
}
}
/* ---------------------------------------------------------------------- */
static
double
t_make_b2
=
0.0
,
n_make_b2
=
0.0
;
template
<
typename
intype
,
typename
outtype
,
int
ni
,
int
nj
>
void
fmatconv
(
intype
*
array
)
{
outtype
*
cast
=
(
outtype
*
)
array
;
for
(
int
i
=
0
;
i
<
ni
;
i
++
)
for
(
int
j
=
0
;
j
<
nj
;
j
++
)
cast
[
i
*
nj
+
j
]
=
array
[
i
*
nj
+
j
];
}
void
PairMGPT
::
make_bond
(
const
double
xx
[][
3
],
int
i
,
int
j
,
bond_data
*
bptr
)
{
double
rrij
[
3
],
rij
;
int
p
;
double
t0
,
t1
;
/* Check that alignment requirements for SIMD code are fullfilled */
assert
(
(((
unsigned
long
long
int
)
(
bptr
->
H
.
m
))
&
31
)
==
0
);
assert
(
(((
unsigned
long
long
int
)
(
bptr
->
Hx
.
m
))
&
31
)
==
0
);
assert
(
(((
unsigned
long
long
int
)
(
bptr
->
Hy
.
m
))
&
31
)
==
0
);
assert
(
(((
unsigned
long
long
int
)
(
bptr
->
Hz
.
m
))
&
31
)
==
0
);
rij
=
0.0
;
for
(
p
=
0
;
p
<
3
;
p
++
)
{
rrij
[
p
]
=
xx
[
i
][
p
]
-
xx
[
j
][
p
];
rij
=
rij
+
rrij
[
p
]
*
rrij
[
p
];
}
/* Zero all matrix elements */
for
(
i
=
0
;
i
<
8
;
i
++
)
for
(
j
=
0
;
j
<
8
;
j
++
)
{
bptr
->
H
.
m
[
i
][
j
]
=
0.0
;
bptr
->
Hx
.
m
[
i
][
j
]
=
0.0
;
bptr
->
Hy
.
m
[
i
][
j
]
=
0.0
;
bptr
->
Hz
.
m
[
i
][
j
]
=
0.0
;
bptr
->
Hz
.
m
[
j
][
i
]
=
0.0
;
}
if
(
rij
<=
rcrit
*
rcrit
)
{
t0
=
gettime
();
if
(
lang
==
3
)
{
hamltn_5_raw
(
rrij
[
0
],
rrij
[
1
],
rrij
[
2
],
bptr
->
H
.
m
,
bptr
->
Hx
.
m
,
bptr
->
Hy
.
m
,
bptr
->
Hz
.
m
,
&
bptr
->
fl_deriv_sum
);
}
else
{
hamltn_7_raw
(
rrij
[
0
],
rrij
[
1
],
rrij
[
2
],
bptr
->
H
.
m
,
bptr
->
Hx
.
m
,
bptr
->
Hy
.
m
,
bptr
->
Hz
.
m
,
&
bptr
->
fl_deriv_sum
);
}
t1
=
gettime
();
t_make_b2
+=
t1
-
t0
;
n_make_b2
++
;
}
else
{
bptr
->
fl_deriv_sum
=
0.0
;
}
if
(
linalg
.
single
)
{
fmatconv
<
double
,
float
,
7
,
8
>
(
&
(
bptr
->
H
.
m
[
1
][
0
]));
fmatconv
<
double
,
float
,
7
,
8
>
(
&
(
bptr
->
Hx
.
m
[
1
][
0
]));
fmatconv
<
double
,
float
,
7
,
8
>
(
&
(
bptr
->
Hy
.
m
[
1
][
0
]));
fmatconv
<
double
,
float
,
7
,
8
>
(
&
(
bptr
->
Hz
.
m
[
1
][
0
]));
}
}
static
double
t_trace
=
0.0
,
n_trace
=
0.0
;
/*
static inline double mtrace(int n,double A[8][8],double B[8][8]) {
double t0,t1;
double s;
t0 = gettime();
if(n == 5) s = mtrace_5(A,B);
else if(n == 7) s = mtrace_7(A,B);
else {
s = 0.0;
for(int i = 1; i<=n; i++)
for(int j = 1; j<=n; j++)
s = s + A[i][j]*B[i][j];
}
t1 = gettime();
t_trace += t1-t0;
n_trace++;
return s;
}
*/
void
PairMGPT
::
make_triplet
(
bond_data
*
ij_bond
,
bond_data
*
ik_bond
,
triplet_data
*
triptr
)
{
if
(
1
)
{
const
trmul_fun
tr_mul
=
linalg
.
tr_mul
;
tr_mul
(
&
(
ij_bond
->
H
.
m
[
1
][
0
]),
&
(
ik_bond
->
H
.
m
[
1
][
0
])
,
&
(
triptr
->
H1H2
.
m
[
1
][
0
])
);
tr_mul
(
&
(
ij_bond
->
Hx
.
m
[
1
][
0
]),
&
(
ik_bond
->
H
.
m
[
1
][
0
])
,
&
(
triptr
->
H1xH2
.
m
[
1
][
0
]));
tr_mul
(
&
(
ij_bond
->
Hy
.
m
[
1
][
0
]),
&
(
ik_bond
->
H
.
m
[
1
][
0
])
,
&
(
triptr
->
H1yH2
.
m
[
1
][
0
]));
tr_mul
(
&
(
ij_bond
->
Hz
.
m
[
1
][
0
]),
&
(
ik_bond
->
H
.
m
[
1
][
0
])
,
&
(
triptr
->
H1zH2
.
m
[
1
][
0
]));
tr_mul
(
&
(
ij_bond
->
H
.
m
[
1
][
0
])
,
&
(
ik_bond
->
Hx
.
m
[
1
][
0
]),
&
(
triptr
->
H1H2x
.
m
[
1
][
0
]));
tr_mul
(
&
(
ij_bond
->
H
.
m
[
1
][
0
])
,
&
(
ik_bond
->
Hy
.
m
[
1
][
0
]),
&
(
triptr
->
H1H2y
.
m
[
1
][
0
]));
tr_mul
(
&
(
ij_bond
->
H
.
m
[
1
][
0
])
,
&
(
ik_bond
->
Hz
.
m
[
1
][
0
]),
&
(
triptr
->
H1H2z
.
m
[
1
][
0
]));
}
else
{
transprod
(
ij_bond
->
H
,
ik_bond
->
H
,
triptr
->
H1H2
);
transprod
(
ij_bond
->
Hx
,
ik_bond
->
H
,
triptr
->
H1xH2
);
transprod
(
ij_bond
->
Hy
,
ik_bond
->
H
,
triptr
->
H1yH2
);
transprod
(
ij_bond
->
Hz
,
ik_bond
->
H
,
triptr
->
H1zH2
);
transprod
(
ij_bond
->
H
,
ik_bond
->
Hx
,
triptr
->
H1H2x
);
transprod
(
ij_bond
->
H
,
ik_bond
->
Hy
,
triptr
->
H1H2y
);
transprod
(
ij_bond
->
H
,
ik_bond
->
Hz
,
triptr
->
H1H2z
);
}
}
static
double
t_make_t
=
0.0
,
t_make_b
=
0.0
,
n_make
=
0.0
;
PairMGPT
::
triplet_data
*
PairMGPT
::
get_triplet
(
const
double
xx
[][
3
],
int
i
,
int
j
,
int
k
,
Hash
<
bond_data
,
Doublet
>
*
bhash
,
triplet_data
*
twork
,
double
*
dvir_ij_p
,
double
*
dvir_ik_p
)
{
const
int
recompute
=
0
;
static
bond_data
bij_work
,
bik_work
;
double
t0
,
t1
;
bond_data
*
bij
=
0
,
*
bik
=
0
;
triplet_data
*
tptr
=
0
;
t0
=
gettime
();
if
(
recompute
==
0
)
{
bij
=
bhash
->
Lookup
(
Doublet
(
i
,
j
));
bik
=
bhash
->
Lookup
(
Doublet
(
i
,
k
));
}
if
(
bij
==
0
)
{
if
(
recompute
==
0
)
bij
=
bhash
->
Insert
(
Doublet
(
i
,
j
));
else
bij
=
&
bij_work
;
if
(
i
<
j
)
make_bond
(
xx
,
i
,
j
,
bij
);
else
make_bond
(
xx
,
j
,
i
,
bij
);
}
if
(
bik
==
0
)
{
if
(
recompute
==
0
)
bik
=
bhash
->
Insert
(
Doublet
(
i
,
k
));
else
bik
=
&
bik_work
;
if
(
i
<
k
)
make_bond
(
xx
,
i
,
k
,
bik
);
else
make_bond
(
xx
,
k
,
i
,
bik
);
}
t1
=
gettime
();
t_make_b
+=
t1
-
t0
;
t0
=
gettime
();
if
(
bij
!=
0
&&
bik
!=
0
)
{
tptr
=
twork
;
make_triplet
(
bij
,
bik
,
tptr
);
*
dvir_ij_p
=
bij
->
fl_deriv_sum
;
*
dvir_ik_p
=
bik
->
fl_deriv_sum
;
}
else
{
*
dvir_ij_p
=
0.0
;
*
dvir_ik_p
=
0.0
;
}
t1
=
gettime
();
t_make_t
+=
t1
-
t0
;
n_make
++
;
return
tptr
;
}
double
PairMGPT
::
numderiv3t
(
double
xx
[][
3
],
int
i
,
int
j
,
int
k
,
int
p
)
{
static
bond_data
Bij
,
Bjk
,
Bki
;
const
double
delta
=
1e-5
;
const
double
xsave
=
xx
[
i
][
p
];
double
e1
,
e2
;
const
double
vc
=
splinepot
.
vc
;
xx
[
i
][
p
]
=
xsave
+
delta
;
make_bond
(
xx
,
i
,
j
,
&
Bij
);
make_bond
(
xx
,
j
,
k
,
&
Bjk
);
make_bond
(
xx
,
k
,
i
,
&
Bki
);
e1
=
trace
(
prodmat
(
Bij
.
H
,
Bjk
.
H
),
Bki
.
H
)
*
(
vc
/
anorm3
);
xx
[
i
][
p
]
=
xsave
-
delta
;
make_bond
(
xx
,
i
,
j
,
&
Bij
);
if
(
0
)
{
/* This bond doesn't change when i is perturbed */
make_bond
(
xx
,
j
,
k
,
&
Bjk
);
}
make_bond
(
xx
,
k
,
i
,
&
Bki
);
e2
=
trace
(
prodmat
(
Bij
.
H
,
Bjk
.
H
),
Bki
.
H
)
*
(
vc
/
anorm3
);
xx
[
i
][
p
]
=
xsave
;
return
(
e1
-
e2
)
/
(
2.0
*
delta
);
}
double
PairMGPT
::
numderiv3v
(
double
xx
[][
3
],
int
i
,
int
j
,
int
k
,
int
p
,
int
ipert
)
{
static
bond_data
Bij
,
Bik
;
const
double
delta
=
1e-5
;
const
double
xsave
=
xx
[
ipert
][
p
];
double
e1
,
e2
;
const
double
vd
=
splinepot
.
vd
;
xx
[
ipert
][
p
]
=
xsave
+
delta
;
make_bond
(
xx
,
i
,
j
,
&
Bij
);
make_bond
(
xx
,
i
,
k
,
&
Bik
);
e1
=
trace
(
prodmat
(
Bij
.
H
,
Bij
.
H
),
prodmat
(
Bik
.
H
,
Bik
.
H
))
*
(
vd
/
anorm4
);
xx
[
ipert
][
p
]
=
xsave
-
delta
;
make_bond
(
xx
,
i
,
j
,
&
Bij
);
make_bond
(
xx
,
i
,
k
,
&
Bik
);
e2
=
trace
(
prodmat
(
Bij
.
H
,
Bij
.
H
),
prodmat
(
Bik
.
H
,
Bik
.
H
))
*
(
vd
/
anorm4
);
xx
[
ipert
][
p
]
=
xsave
;
return
(
e1
-
e2
)
/
(
2.0
*
delta
);
}
double
PairMGPT
::
numderiv4
(
double
xx
[][
3
],
int
i
,
int
j
,
int
k
,
int
m
,
int
p
)
{
static
bond_data
Bij
,
Bjk
,
Bkm
,
Bmi
;
const
double
delta
=
1e-5
;
const
double
xsave
=
xx
[
i
][
p
];
double
e1
,
e2
;
const
double
ve
=
splinepot
.
ve
;
xx
[
i
][
p
]
=
xsave
+
delta
;
make_bond
(
xx
,
i
,
j
,
&
Bij
);
make_bond
(
xx
,
j
,
k
,
&
Bjk
);
make_bond
(
xx
,
k
,
m
,
&
Bkm
);
make_bond
(
xx
,
m
,
i
,
&
Bmi
);
e1
=
trace
(
prodmat
(
Bij
.
H
,
Bjk
.
H
),
prodmat
(
Bkm
.
H
,
Bmi
.
H
))
*
(
ve
/
anorm4
);
xx
[
i
][
p
]
=
xsave
-
delta
;
make_bond
(
xx
,
i
,
j
,
&
Bij
);
if
(
0
)
{
/* Only the i coordinates changed... */
make_bond
(
xx
,
j
,
k
,
&
Bjk
);
make_bond
(
xx
,
k
,
m
,
&
Bkm
);
}
make_bond
(
xx
,
m
,
i
,
&
Bmi
);
e2
=
trace
(
prodmat
(
Bij
.
H
,
Bjk
.
H
),
prodmat
(
Bkm
.
H
,
Bmi
.
H
))
*
(
ve
/
anorm4
);
xx
[
i
][
p
]
=
xsave
;
return
(
e1
-
e2
)
/
(
2.0
*
delta
);
}
static
double
dtol
=
1e-6
;
void
PairMGPT
::
force_debug_3t
(
double
xx
[][
3
],
int
i0
,
int
j0
,
int
k0
,
int
i
,
int
j
,
int
k
,
double
dfix
,
double
dfiy
,
double
dfiz
,
double
dfjx
,
double
dfjy
,
double
dfjz
,
double
dfkx
,
double
dfky
,
double
dfkz
)
{
double
dfi
[
3
],
dfj
[
3
],
dfk
[
3
];
dfi
[
0
]
=
dfix
;
dfi
[
1
]
=
dfiy
;
dfi
[
2
]
=
dfiz
;
dfj
[
0
]
=
dfjx
;
dfj
[
1
]
=
dfjy
;
dfj
[
2
]
=
dfjz
;
dfk
[
0
]
=
dfkx
;
dfk
[
1
]
=
dfky
;
dfk
[
2
]
=
dfkz
;
for
(
int
p
=
0
;
p
<
3
;
p
++
)
{
/* Compute numerical derivatives by displacing atoms i,j,k */
double
ndfi
,
ndfj
,
ndfk
;
ndfi
=
-
numderiv3t
(
xx
,
i
,
j
,
k
,
p
);
ndfj
=
-
numderiv3t
(
xx
,
j
,
k
,
i
,
p
);
ndfk
=
-
numderiv3t
(
xx
,
k
,
i
,
j
,
p
);
if
((
fabs
(
dfi
[
p
]
-
ndfi
)
>
dtol
&&
fabs
(
dfi
[
p
]
-
ndfi
)
>
dtol
*
fabs
(
ndfi
))
||
(
fabs
(
dfj
[
p
]
-
ndfj
)
>
dtol
&&
fabs
(
dfj
[
p
]
-
ndfj
)
>
dtol
*
fabs
(
ndfj
))
||
(
fabs
(
dfk
[
p
]
-
ndfk
)
>
dtol
&&
fabs
(
dfk
[
p
]
-
ndfk
)
>
dtol
*
fabs
(
ndfk
)))
{
printf
(
"Force error in T12 & T23 & T31 :: i,j,k = %d,%d,%d
\n
"
,
i0
,
j0
,
k0
);
printf
(
" dE/d%c[i] = %20.10e %20.10e
\n
"
,
'x'
+
p
,
ndfi
,
dfi
[
p
]);
printf
(
" dE/d%c[j] = %20.10e %20.10e
\n
"
,
'x'
+
p
,
ndfj
,
dfj
[
p
]);
printf
(
" dE/d%c[k] = %20.10e %20.10e
\n
"
,
'x'
+
p
,
ndfk
,
dfk
[
p
]);
printf
(
"
\n
"
);
}
}
}
void
PairMGPT
::
force_debug_3v
(
double
xx
[][
3
],
int
i0
,
int
j0
,
int
k0
,
int
i
,
int
j
,
int
k
,
double
dfix
,
double
dfiy
,
double
dfiz
,
double
dfjx
,
double
dfjy
,
double
dfjz
,
double
dfkx
,
double
dfky
,
double
dfkz
)
{
double
dfi
[
3
],
dfj
[
3
],
dfk
[
3
];
dfi
[
0
]
=
dfix
;
dfi
[
1
]
=
dfiy
;
dfi
[
2
]
=
dfiz
;
dfj
[
0
]
=
dfjx
;
dfj
[
1
]
=
dfjy
;
dfj
[
2
]
=
dfjz
;
dfk
[
0
]
=
dfkx
;
dfk
[
1
]
=
dfky
;
dfk
[
2
]
=
dfkz
;
for
(
int
p
=
0
;
p
<
3
;
p
++
)
{
/* Compute numerical derivatives by displacing atoms i,j,k */
double
ndfi
,
ndfj
,
ndfk
;
ndfi
=
-
numderiv3v
(
xx
,
i
,
j
,
k
,
p
,
i0
);
ndfj
=
-
numderiv3v
(
xx
,
i
,
j
,
k
,
p
,
j0
);
ndfk
=
-
numderiv3v
(
xx
,
i
,
j
,
k
,
p
,
k0
);
if
((
fabs
(
dfi
[
p
]
-
ndfi
)
>
dtol
&&
fabs
(
dfi
[
p
]
-
ndfi
)
>
dtol
*
fabs
(
ndfi
))
||
(
fabs
(
dfj
[
p
]
-
ndfj
)
>
dtol
&&
fabs
(
dfj
[
p
]
-
ndfj
)
>
dtol
*
fabs
(
ndfj
))
||
(
fabs
(
dfk
[
p
]
-
ndfk
)
>
dtol
&&
fabs
(
dfk
[
p
]
-
ndfk
)
>
dtol
*
fabs
(
ndfk
)))
{
printf
(
"Force error in T12 :: i,j,k = %d,%d,%d
\n
"
,
i0
,
j0
,
k0
);
printf
(
" dE/d%c[i] = %20.10e %20.10e
\n
"
,
'x'
+
p
,
ndfi
,
dfi
[
p
]);
printf
(
" dE/d%c[j] = %20.10e %20.10e
\n
"
,
'x'
+
p
,
ndfj
,
dfj
[
p
]);
printf
(
" dE/d%c[k] = %20.10e %20.10e
\n
"
,
'x'
+
p
,
ndfk
,
dfk
[
p
]);
printf
(
"
\n
"
);
}
}
}
void
PairMGPT
::
force_debug_4
(
double
xx
[][
3
],
int
i0
,
int
j0
,
int
k0
,
int
m0
,
int
i
,
int
j
,
int
k
,
int
m
,
double
dfix
,
double
dfiy
,
double
dfiz
,
double
dfjx
,
double
dfjy
,
double
dfjz
,
double
dfkx
,
double
dfky
,
double
dfkz
,
double
dfmx
,
double
dfmy
,
double
dfmz
)
{
double
dfi
[
3
],
dfj
[
3
],
dfk
[
3
],
dfm
[
3
];
dfi
[
0
]
=
dfix
;
dfi
[
1
]
=
dfiy
;
dfi
[
2
]
=
dfiz
;
dfj
[
0
]
=
dfjx
;
dfj
[
1
]
=
dfjy
;
dfj
[
2
]
=
dfjz
;
dfk
[
0
]
=
dfkx
;
dfk
[
1
]
=
dfky
;
dfk
[
2
]
=
dfkz
;
dfm
[
0
]
=
dfmx
;
dfm
[
1
]
=
dfmy
;
dfm
[
2
]
=
dfmz
;
const
int
ii0
[]
=
{
i0
,
j0
,
k0
,
m0
},
ii
[]
=
{
i
,
j
,
k
,
m
,
i
,
j
,
k
};
for
(
int
p
=
0
;
p
<
3
;
p
++
)
{
/* Compute numerical derivatives by displacing atoms i,j,k,m */
double
ndfi
,
ndfj
,
ndfk
,
ndfm
;
if
(
1
)
{
double
ndf
[]
=
{
0.0
,
0.0
,
0.0
,
0.0
};
for
(
int
s
=
0
;
s
<
4
;
s
++
)
for
(
int
t
=
0
;
t
<
4
;
t
++
)
if
(
ii
[
s
]
==
ii0
[
t
])
ndf
[
t
]
=
-
numderiv4
(
xx
,
ii
[
s
],
ii
[
s
+
1
],
ii
[
s
+
2
],
ii
[
s
+
3
],
p
);
ndfi
=
ndf
[
0
];
ndfj
=
ndf
[
1
];
ndfk
=
ndf
[
2
];
ndfm
=
ndf
[
3
];
}
else
{
ndfi
=
-
numderiv4
(
xx
,
i
,
j
,
k
,
m
,
p
);
ndfj
=
-
numderiv4
(
xx
,
j
,
k
,
m
,
i
,
p
);
ndfk
=
-
numderiv4
(
xx
,
k
,
m
,
i
,
j
,
p
);
ndfm
=
-
numderiv4
(
xx
,
m
,
i
,
j
,
k
,
p
);
}
if
((
fabs
(
dfi
[
p
]
-
ndfi
)
>
dtol
&&
fabs
(
dfi
[
p
]
-
ndfi
)
>
dtol
*
fabs
(
ndfi
))
||
(
fabs
(
dfj
[
p
]
-
ndfj
)
>
dtol
&&
fabs
(
dfj
[
p
]
-
ndfj
)
>
dtol
*
fabs
(
ndfj
))
||
(
fabs
(
dfk
[
p
]
-
ndfk
)
>
dtol
&&
fabs
(
dfk
[
p
]
-
ndfk
)
>
dtol
*
fabs
(
ndfk
))
||
(
fabs
(
dfm
[
p
]
-
ndfm
)
>
dtol
&&
fabs
(
dfm
[
p
]
-
ndfm
)
>
dtol
*
fabs
(
ndfm
)))
{
printf
(
"Force error in T31 & T64 :: i,j,k,m = %d,%d,%d,%d
\n
"
,
i0
,
j0
,
k0
,
m0
);
printf
(
" dE/d%c[i] = %20.10e %20.10e
\n
"
,
'x'
+
p
,
ndfi
,
dfi
[
p
]);
printf
(
" dE/d%c[j] = %20.10e %20.10e
\n
"
,
'x'
+
p
,
ndfj
,
dfj
[
p
]);
printf
(
" dE/d%c[k] = %20.10e %20.10e
\n
"
,
'x'
+
p
,
ndfk
,
dfk
[
p
]);
printf
(
" dE/d%c[m] = %20.10e %20.10e
\n
"
,
'x'
+
p
,
ndfm
,
dfm
[
p
]);
printf
(
"
\n
"
);
}
}
}
/*
#define trd_update_4(T12,T45,coord) \
do { \
trd1 = transtrace(T12->H1##coord##H2,T45->H1H2 ); \
trd2 = transtrace(T12->H1H2##coord,T45->H1H2 ); \
trd3 = transtrace(T12->H1H2 ,T45->H1##coord##H2); \
trd4 = transtrace(T12->H1H2 ,T45->H1H2##coord ); \
} while(0)
*/
#define trd_update_4(T12,T45) \
do { \
tr_trace3(&(T45->H1H2.m[1][0]), \
&(T12->H1xH2.m[1][0]),&utr1x.d, \
&(T12->H1yH2.m[1][0]),&utr1y.d, \
&(T12->H1zH2.m[1][0]),&utr1z.d); \
tr_trace3(&(T45->H1H2.m[1][0]), \
&(T12->H1H2x.m[1][0]),&utr2x.d, \
&(T12->H1H2y.m[1][0]),&utr2y.d, \
&(T12->H1H2z.m[1][0]),&utr2z.d); \
tr_trace3(&(T12->H1H2.m[1][0]), \
&(T45->H1xH2.m[1][0]),&utr3x.d, \
&(T45->H1yH2.m[1][0]),&utr3y.d, \
&(T45->H1zH2.m[1][0]),&utr3z.d); \
tr_trace3(&(T12->H1H2.m[1][0]), \
&(T45->H1H2x.m[1][0]),&utr4x.d, \
&(T45->H1H2y.m[1][0]),&utr4y.d, \
&(T45->H1H2z.m[1][0]),&utr4z.d); \
if(linalg.single) { \
trd1x = utr1x.f; trd2x = utr2x.f; trd3x = utr3x.f; trd4x = utr4x.f; \
trd1y = utr1y.f; trd2y = utr2y.f; trd3y = utr3y.f; trd4y = utr4y.f; \
trd1z = utr1z.f; trd2z = utr2z.f; trd3z = utr3z.f; trd4z = utr4z.f; \
} else { \
trd1x = utr1x.d; trd2x = utr2x.d; trd3x = utr3x.d; trd4x = utr4x.d; \
trd1y = utr1y.d; trd2y = utr2y.d; trd3y = utr3y.d; trd4y = utr4y.d; \
trd1z = utr1z.d; trd2z = utr2z.d; trd3z = utr3z.d; trd4z = utr4z.d; \
} \
} while(0)
#define dfix_update_4a(coord) \
do { \
dfi##coord = ( (-sij)*trd1##coord + (-sim)*trd3##coord ) * (ve / anorm4); \
dfj##coord = ( ( sij)*trd1##coord + (-sjk)*trd2##coord ) * (ve / anorm4); \
dfk##coord = ( ( sjk)*trd2##coord + (-skm)*trd4##coord ) * (ve / anorm4); \
dfm##coord = ( ( sim)*trd3##coord + ( skm)*trd4##coord ) * (ve / anorm4); \
} while(0)
#define dfix_update_4b(coord) \
do { \
dfi##coord = ( ( ski)*trd1##coord + (-sim)*trd3##coord ) * (ve / anorm4); \
dfj##coord = ( (-sjk)*trd2##coord + (-sjm)*trd4##coord ) * (ve / anorm4); \
dfk##coord = ( (-ski)*trd1##coord + ( sjk)*trd2##coord ) * (ve / anorm4); \
dfm##coord = ( ( sim)*trd3##coord + ( sjm)*trd4##coord ) * (ve / anorm4); \
} while(0);
#define dfix_update_4c(coord) \
do { \
dfi##coord = ( (-sij)*trd1##coord + ( ski)*trd2##coord ) * (ve / anorm4); \
dfj##coord = ( ( sij)*trd1##coord + (-sjm)*trd3##coord ) * (ve / anorm4); \
dfk##coord = ( (-ski)*trd2##coord + (-skm)*trd4##coord ) * (ve / anorm4); \
dfm##coord = ( ( sjm)*trd3##coord + ( skm)*trd4##coord ) * (ve / anorm4); \
} while(0);
#define accumulate_forces_2(w) \
do { \
fix = fix + dfix*(w); \
fiy = fiy + dfiy*(w); \
fiz = fiz + dfiz*(w); \
\
fjx = fjx + dfjx*(w); \
fjy = fjy + dfjy*(w); \
fjz = fjz + dfjz*(w); \
} while(0)
#define accumulate_forces_3(w) \
do { \
accumulate_forces_2(w); \
fkx = fkx + dfkx*(w); \
fky = fky + dfky*(w); \
fkz = fkz + dfkz*(w); \
} while(0)
#define accumulate_forces_4(w) \
do { \
accumulate_forces_3(w); \
fmx = fmx + dfmx*(w); \
fmy = fmy + dfmy*(w); \
fmz = fmz + dfmz*(w); \
} while(0)
#define restrict __restrict__
#ifdef __bg__
#define const
#endif
static
int
ntr_calls
=
0
;
static
trtrace3_fun
tr_internal
;
static
void
tr_count
(
const
double
*
restrict
A
,
const
double
*
restrict
B1
,
double
*
restrict
t1
,
const
double
*
restrict
B2
,
double
*
restrict
t2
,
const
double
*
restrict
B3
,
double
*
restrict
t3
)
{
tr_internal
(
A
,
B1
,
t1
,
B2
,
t2
,
B3
,
t3
);
ntr_calls
++
;
}
#ifdef __bg__
#undef const
#endif
#undef restrict
int
PairMGPT
::
Matrix
::
sz
;
void
PairMGPT
::
compute_x
(
const
int
*
nnei
,
const
int
*
const
*
nlist
,
double
*
e_s
,
double
*
e_p
,
double
*
e_t
,
double
*
e_q
,
int
evflag
,
int
newton_pair
)
{
Hash
<
bond_data
,
Doublet
>
bond_hash
(
100000
);
int
i
,
j
,
k
,
m
,
ix
,
jx
,
kx
,
p
;
double
e_single
,
e_pair
,
e_triplet
,
e_triplet_c
,
e_quad
;
double
volvir2
;
double
nbc
=
0.0
,
tbl
=
0.0
,
tbm
=
0.0
;
const
int
lmax_local
=
lmax
;
//if(evflag) printf("##### ev flag is set... wasting cycles...\n");
*
e_s
=
-
99.0
;
*
e_p
=
-
99.0
;
*
e_t
=
-
99.0
;
*
e_q
=
-
99.0
;
double
t0
,
t1
;
t0
=
gettime
(
1
);
e_single
=
e_pair
=
e_triplet
=
e_triplet_c
=
e_quad
=
0.0
;
volvir2
=
0.0
;
t_make_t
=
t_make_b
=
t_make_b2
=
t_trace
=
0.0
;
n_make
=
n_make_b2
=
n_trace
=
0.0
;
double
tx0
,
tx1
,
tsort
=
0.0
,
tpair
=
0.0
,
tlookup
=
0.0
;
double
ttriplet
=
0.0
,
tquad
=
0.0
,
tmem
=
0.0
;
double
ntsort
=
0.0
,
ntpair
=
0.0
,
ntlookup
=
0.0
;
double
nttriplet
=
0.0
,
ntquad
=
0.0
,
ntmem
=
0.0
,
ntquaditer
=
0.0
;
double
mcount
=
0.0
,
mcount2
=
0.0
,
qcount
=
0.0
;
double
fix
,
fjx
,
fkx
,
fmx
,
dfix
,
dfjx
,
dfkx
,
dfmx
;
double
fiy
,
fjy
,
fky
,
fmy
,
dfiy
,
dfjy
,
dfky
,
dfmy
;
double
fiz
,
fjz
,
fkz
,
fmz
,
dfiz
,
dfjz
,
dfkz
,
dfmz
;
double
fsave
[
4
][
3
]
=
{
{
0.0
}
}
/* {{0.0}} is to get rid of uninitialized use warning */
;
//const int numerical_pair_forces = (nbody_flag/16)%2;
const
int
pair_forces
=
(
nbody_flag
/
2
)
%
2
,
three_body_forces
=
(
nbody_flag
/
4
)
%
2
,
four_body_forces
=
(
nbody_flag
/
8
)
%
2
;
const
int
pair_energies
=
(
nbody_flag
/
2
)
%
2
,
three_body_energies
=
(
nbody_flag
/
4
)
%
2
,
four_body_energies
=
(
nbody_flag
/
8
)
%
2
;
const
int
single_energies
=
nbody_flag
%
2
;
const
int
triplet_debug
=
0
,
quad_debug
=
0
;
/* Energy and force scale factor for unit conversion. */
const
double
e_scale
=
0.5
;
#ifdef NEIGHMASK
#define NIDX(x) (x)
#else
#define NIDX(x) ((x) & NEIGHMASK)
#endif
int
nneitot
,
*
first
,
*
nlist_short
;
double
w2
,
w3
,
w4
;
triplet_data
T12work
,
T23work
,
T31work
,
T45work
,
T56work
,
T64work
;
triplet_data
*
T12
,
*
T23
,
*
T31
,
*
T45
,
*
T56
,
*
T64
;
int
c_ij
,
c_jk
,
c_ki
,
c_im
,
c_jm
,
c_km
;
int
mi
,
mj
,
mk
;
double
tr0
,
tr1
,
tr2
,
tr3
;
double
v33
,
v43
;
double
rcut2_pair
=
rmax
*
rmax
,
rcut2_bond
=
rcrit
*
rcrit
,
rij2
;
int
ntot
,
nloc
;
double
dvir_ij
,
dvir_jk
,
dvir_ki
,
dvir_im
,
dvir_jm
,
dvir_km
;
double
vir3t
=
0.0
,
vir3v
=
0.0
,
vir4
=
0.0
;
double
(
*
xx
)[
3
],(
*
ff
)[
3
],(
*
ss
)[
3
];
#ifdef TIMING_ON
tr_internal
=
linalg
.
tr_trace
;
ntr_calls
=
0
;
const
trtrace3_fun
tr_trace3
=
tr_count
;
#else
const
trtrace3_fun
tr_trace3
=
linalg
.
tr_trace
;
#endif
union
{
double
d
;
float
f
;
}
utr1x
,
utr2x
,
utr3x
,
utr4x
,
utr1y
,
utr2y
,
utr3y
,
utr4y
,
utr1z
,
utr2z
,
utr3z
,
utr4z
;
double
trd1x
,
trd2x
,
trd3x
,
trd4x
;
double
trd1y
,
trd2y
,
trd3y
,
trd4y
;
double
trd1z
,
trd2z
,
trd3z
,
trd4z
;
tx0
=
gettime
();
double
rhoinv
;
{
double
vtot
=
1.0
;
double
ntot
=
atom
->
natoms
;
for
(
i
=
0
;
i
<
3
;
i
++
)
vtot
=
vtot
*
(
domain
->
boxhi
[
i
]
-
domain
->
boxlo
[
i
]);
rhoinv
=
vtot
/
ntot
;
}
/* Make sure triplet data work area is aligned and zeroed out. */
{
assert
(
T12work
.
align_check
()
==
0
);
assert
(
T23work
.
align_check
()
==
0
);
assert
(
T31work
.
align_check
()
==
0
);
assert
(
T45work
.
align_check
()
==
0
);
assert
(
T56work
.
align_check
()
==
0
);
assert
(
T64work
.
align_check
()
==
0
);
T12work
.
zero
();
T23work
.
zero
();
T31work
.
zero
();
T45work
.
zero
();
T56work
.
zero
();
T64work
.
zero
();
}
ntot
=
atom
->
nlocal
+
atom
->
nghost
;
nloc
=
atom
->
nlocal
;
//printf("[%3d] Allocating local array, size is %d atoms...\n",comm->me,j);
xx
=
(
double
(
*
)[
3
])
memory
->
smalloc
(
sizeof
(
double
[
3
])
*
ntot
,
"mgpt: local position vector."
);
ff
=
(
double
(
*
)[
3
])
memory
->
smalloc
(
sizeof
(
double
[
3
])
*
ntot
,
"mgpt: local force vector."
);
//printf("[%3d] Initializing arrays...\n",comm->me);
const
int
triclinic
=
domain
->
triclinic
;
double
alpha
[
3
]
=
{
0.0
,
0.0
,
0.0
};
if
(
triclinic
)
{
double
E
[
3
][
3
],
EX
[
3
][
3
];
int
cyc
[]
=
{
0
,
1
,
2
,
0
,
1
};
ss
=
(
double
(
*
)[
3
])
memory
->
smalloc
(
sizeof
(
double
[
3
])
*
ntot
,
"mgpt: local reduced coordinate vector."
);
for
(
i
=
0
;
i
<
3
;
i
++
)
{
for
(
j
=
0
;
j
<
3
;
j
++
)
E
[
i
][
j
]
=
0.0
;
E
[
i
][
i
]
=
domain
->
subhi_lamda
[
i
]
-
domain
->
sublo_lamda
[
i
];
domain
->
lamda2x
(
E
[
i
],
EX
[
i
]);
}
for
(
i
=
0
;
i
<
3
;
i
++
)
{
int
i1
=
cyc
[
i
+
1
],
i2
=
cyc
[
i
+
2
];
double
dot
=
0.0
,
ns2
=
0.0
;
for
(
j
=
0
;
j
<
3
;
j
++
)
{
int
j1
=
cyc
[
j
+
1
],
j2
=
cyc
[
j
+
2
];
double
cj
=
EX
[
i1
][
j1
]
*
EX
[
i2
][
j2
]
-
EX
[
i1
][
j2
]
*
EX
[
i2
][
j1
];
ns2
=
ns2
+
cj
*
cj
;
dot
=
dot
+
EX
[
i
][
j
]
*
cj
;
}
alpha
[
i
]
=
E
[
i
][
i
]
/
(
dot
/
sqrt
(
ns2
));
if
(
comm
->
me
==
0
)
{
static
int
count
=
0
;
if
(
count
<
3
)
printf
(
"@@@ alpha(%d) = %15.5e
\n
"
,
i
+
1
,
alpha
[
i
]);
count
++
;
}
if
(
alpha
[
i
]
<
0.0
)
alpha
[
i
]
=
-
alpha
[
i
];
}
}
else
ss
=
xx
;
nneitot
=
0
;
for
(
ix
=
0
;
ix
<
ntot
;
ix
++
)
{
for
(
p
=
0
;
p
<
3
;
p
++
)
{
xx
[
ix
][
p
]
=
atom
->
x
[
ix
][
p
];
ff
[
ix
][
p
]
=
0.0
;
}
if
(
triclinic
)
domain
->
x2lamda
(
xx
[
ix
],
ss
[
ix
]);
nneitot
=
nneitot
+
nnei
[
ix
];
}
first
=
(
int
*
)
memory
->
smalloc
(
sizeof
(
int
)
*
(
ntot
+
1
),
"mgpt: first"
);
nlist_short
=
(
int
*
)
memory
->
smalloc
(
sizeof
(
int
)
*
nneitot
,
"mgpt: nlist_short"
);
tx1
=
gettime
();
tmem
+=
tx1
-
tx0
;
ntmem
++
;
//printf("[%3d] Starting calculation...\n",comm->me);
fix
=
fjx
=
fkx
=
fmx
=
0.0
;
fiy
=
fjy
=
fky
=
fmy
=
0.0
;
fiz
=
fjz
=
fkz
=
fmz
=
0.0
;
int
c_p
=
0
,
c_t
=
0
,
c_q
=
0
;
if
(
0
)
if
(
domain
->
triclinic
)
{
if
(
comm
->
me
==
0
)
printf
(
"Can not handle triclinic box yet
\n
"
);
error
->
all
(
__FILE__
,
__LINE__
,
"Can not handle triclinic cell with mgpt yet."
);
}
/*
for(i = 0; i<nloc; i++) {
printf("Atom %3d:: %10.3f %10.3f %10.3f\n",
i,xx[i][0],xx[i][1],xx[i][2]);
}
*/
first
[
0
]
=
0
;
for
(
i
=
0
;
i
<
ntot
;
i
++
)
{
fix
=
fiy
=
fiz
=
0.0
;
first
[
i
+
1
]
=
first
[
i
];
const
int
c1
=
c1_outside
(
ss
[
i
],
triclinic
,
alpha
);
tx0
=
gettime
();
for
(
jx
=
0
;
jx
<
nnei
[
i
];
jx
++
)
{
fjx
=
fjy
=
fjz
=
0.0
;
j
=
NIDX
(
nlist
[
i
][
jx
]
);
rij2
=
0.0
;
for
(
p
=
0
;
p
<
3
;
p
++
)
{
double
t
=
xx
[
i
][
p
]
-
xx
[
j
][
p
];
rij2
=
rij2
+
t
*
t
;
}
if
(
c1
==
0
&&
rij2
<
rcut2_pair
)
{
if
(
j
<
i
)
{
w2
=
get_weight
(
triclinic
,
ss
[
i
],
ss
[
j
]);
if
(
w2
>
0.0
)
{
/*
Compute pair energy/force
*/
double
de_pair
,
df
,
rij
=
sqrt
(
rij2
);
splinepot
.
eval_pot
(
rij
,
&
de_pair
,
&
df
);
de_pair
=
de_pair
*
e_scale
*
w2
;
df
=
df
/
rij
*
w2
;
if
(
pair_energies
==
0
)
de_pair
=
0.0
;
e_pair
=
e_pair
+
de_pair
;
c_p
++
;
if
(
pair_forces
==
0
)
df
=
0.0
;
if
(
volpres_flag
&&
pair_energies
)
{
double
dvir
;
splinepot
.
eval_vir
(
rij
,
&
dvir
);
volvir2
=
volvir2
-
dvir
*
w2
;
/* Per-atom virial contribution of volumetric energy term */
if
(
vflag_atom
)
for
(
int
pp
=
0
;
pp
<
3
;
pp
++
)
{
//virial[i] = virial[i] + rhoinv*e_scale*volvir2;
vatom
[
i
][
pp
]
-=
0.5
*
rhoinv
*
e_scale
*
dvir
*
w2
;
vatom
[
j
][
pp
]
-=
0.5
*
rhoinv
*
e_scale
*
dvir
*
w2
;
}
}
double
drijx
=
xx
[
j
][
0
]
-
xx
[
i
][
0
];
double
drijy
=
xx
[
j
][
1
]
-
xx
[
i
][
1
];
double
drijz
=
xx
[
j
][
2
]
-
xx
[
i
][
2
];
fix
=
fix
+
df
*
drijx
;
fjx
=
fjx
-
df
*
drijx
;
fiy
=
fiy
+
df
*
drijy
;
fjy
=
fjy
-
df
*
drijy
;
fiz
=
fiz
+
df
*
drijz
;
fjz
=
fjz
-
df
*
drijz
;
if
(
evflag
)
{
//ev_tally(i,j,nloc,newton_pair,de_pair,0.0,df,-drijx,-drijy,-drijz);
/* To fix stress-per-atom scaling, and sign */
ev_tally
(
i
,
j
,
nloc
,
newton_pair
,
de_pair
,
0.0
,
-
df
*
e_scale
,
-
drijx
,
-
drijy
,
-
drijz
);
}
ff
[
j
][
0
]
+=
fjx
*
e_scale
;
ff
[
j
][
1
]
+=
fjy
*
e_scale
;
ff
[
j
][
2
]
+=
fjz
*
e_scale
;
}
}
}
if
(
rij2
<
rcut2_bond
&&
c2_outside
(
ss
[
i
],
ss
[
j
],
triclinic
,
alpha
)
==
0
)
{
/*
Add j to short neighbor list for i.
Insert j to keep list sorted.
*/
p
=
first
[
i
+
1
]
-
1
;
while
(
p
>=
first
[
i
]
&&
nlist_short
[
p
]
>
j
)
{
nlist_short
[
p
+
1
]
=
nlist_short
[
p
];
p
=
p
-
1
;
}
nlist_short
[
p
+
1
]
=
j
;
first
[
i
+
1
]
=
first
[
i
+
1
]
+
1
;
if
(
first
[
i
+
1
]
>
nneitot
)
{
printf
(
"nneitot = %d, short list full. i=%d
\n
"
,
nneitot
,
i
);
error
->
one
(
__FILE__
,
__LINE__
,
"Shit! Short list full
\n
"
);
}
}
}
ff
[
i
][
0
]
+=
fix
*
e_scale
;
ff
[
i
][
1
]
+=
fiy
*
e_scale
;
ff
[
i
][
2
]
+=
fiz
*
e_scale
;
tx1
=
gettime
();
tpair
+=
tx1
-
tx0
;
ntpair
+=
nnei
[
i
];
}
for
(
i
=
0
;
i
<
ntot
;
i
++
)
{
fix
=
fiy
=
fiz
=
0.0
;
/*
Use short lists for triplets and quadruplets.
For open (2-bonded) triplets, can only use k<j, but not k<i.
For closed (3-bonded) triplets, we can assume k<j<i.
Quadruplets:
Always use k<j<i, and require m<i.
If 5-bonded with im bond, ignore the quadruplet.
If 6-bonded, require m<k.
For 4-bonded quadruplets, we can still use k<j, but also
assume max(m,j)<i
*/
if
(
three_body_energies
||
three_body_forces
||
four_body_energies
||
four_body_forces
)
for
(
jx
=
first
[
i
];
jx
<
first
[
i
+
1
];
jx
++
)
{
fjx
=
fjy
=
fjz
=
0.0
;
j
=
nlist_short
[
jx
];
for
(
kx
=
first
[
i
];
kx
<
jx
;
kx
++
)
{
fkx
=
fky
=
fkz
=
0.0
;
k
=
nlist_short
[
kx
];
/*
Search lists of j and k, and see if
1) j is in k-list (closed triplet)
2) j and k have a common neighbor (closed quadruplet)
*/
c_ij
=
c_ki
=
1
;
const
int
sij
=
(
i
<
j
)
?
1
:
-
1
;
const
int
sjk
=
(
j
<
k
)
?
1
:
-
1
;
const
int
ski
=
(
k
<
i
)
?
1
:
-
1
;
T12
=
T23
=
T31
=
0
;
mj
=
first
[
j
];
/*
Since i is in the j-list, and i > k and the list
is sorted, the loop below terminates:-)
*/
while
(
mj
<
first
[
j
+
1
]
&&
nlist_short
[
mj
]
<
k
)
mj
=
mj
+
1
;
if
(
mj
<
first
[
j
+
1
]
&&
nlist_short
[
mj
]
==
k
)
{
/* Closed triplet */
c_jk
=
1
;
if
(
j
>
i
)
continue
;
/* Require k<j<i for closed triplets */
}
else
{
/* Open triplet */
c_jk
=
0
;
}
tx0
=
gettime
();
w3
=
get_weight
(
triclinic
,
ss
[
i
],
ss
[
j
],
ss
[
k
]);
int
triplet_defer
;
if
(
w3
>
0.0
)
{
triplet_defer
=
0
;
dvir_ij
=
dvir_jk
=
dvir_ki
=
0.0
;
if
(
c_ij
&&
c_jk
)
T12
=
get_triplet
(
xx
,
j
,
i
,
k
,
&
bond_hash
,
&
T12work
,
&
dvir_ij
,
&
dvir_jk
);
if
(
c_ki
&&
c_jk
)
T23
=
get_triplet
(
xx
,
k
,
i
,
j
,
&
bond_hash
,
&
T23work
,
&
dvir_ki
,
&
dvir_jk
);
if
(
c_ij
&&
c_ki
)
T31
=
get_triplet
(
xx
,
i
,
j
,
k
,
&
bond_hash
,
&
T31work
,
&
dvir_ij
,
&
dvir_ki
);
if
(
evflag
)
{
fsave
[
0
][
0
]
=
fix
;
fsave
[
0
][
1
]
=
fiy
;
fsave
[
0
][
2
]
=
fiz
;
fsave
[
1
][
0
]
=
fjx
;
fsave
[
1
][
1
]
=
fjy
;
fsave
[
1
][
2
]
=
fjz
;
fsave
[
2
][
0
]
=
fkx
;
fsave
[
2
][
1
]
=
fky
;
fsave
[
2
][
2
]
=
fkz
;
fix
=
fiy
=
fiz
=
0.0
;
fjx
=
fjy
=
fjz
=
0.0
;
fkx
=
fky
=
fkz
=
0.0
;
}
tr0
=
tr1
=
tr2
=
tr3
=
0.0
;
double
xvir3t
,
xvir3v
;
xvir3t
=
xvir3v
=
0.0
;
if
(
T12
&&
T23
)
{
bond_data
*
bki
=
bond_hash
.
Lookup
(
Doublet
(
k
,
i
));
if
(
three_body_energies
&&
evflag
)
{
tr0
=
transtrace
(
T12
->
H1H2
,
bki
->
H
);
double
dvir
=
((
dvir_ij
+
dvir_jk
+
bki
->
fl_deriv_sum
)
*
splinepot
.
vc
+
splinepot
.
dvc
)
*
tr0
*
w3
/
anorm3
;
vir3t
=
vir3t
+
dvir
;
xvir3t
=
xvir3t
+
dvir
;
}
mcount2
++
;
{
const
double
vc
=
splinepot
.
vc
;
tr_trace3
(
&
(
bki
->
H
.
m
[
1
][
0
]),
&
(
T12
->
H1xH2
.
m
[
1
][
0
]),
&
utr1x
.
d
,
&
(
T12
->
H1yH2
.
m
[
1
][
0
]),
&
utr1y
.
d
,
&
(
T12
->
H1zH2
.
m
[
1
][
0
]),
&
utr1z
.
d
);
tr_trace3
(
&
(
bki
->
H
.
m
[
1
][
0
]),
&
(
T12
->
H1H2x
.
m
[
1
][
0
]),
&
utr2x
.
d
,
&
(
T12
->
H1H2y
.
m
[
1
][
0
]),
&
utr2y
.
d
,
&
(
T12
->
H1H2z
.
m
[
1
][
0
]),
&
utr2z
.
d
);
tr_trace3
(
&
(
T12
->
H1H2
.
m
[
1
][
0
]),
&
(
bki
->
Hx
.
m
[
1
][
0
]),
&
utr3x
.
d
,
&
(
bki
->
Hy
.
m
[
1
][
0
]),
&
utr3y
.
d
,
&
(
bki
->
Hz
.
m
[
1
][
0
]),
&
utr3z
.
d
);
if
(
linalg
.
single
)
{
trd1x
=
utr1x
.
f
;
trd2x
=
utr2x
.
f
;
trd3x
=
utr3x
.
f
;
trd1y
=
utr1y
.
f
;
trd2y
=
utr2y
.
f
;
trd3y
=
utr3y
.
f
;
trd1z
=
utr1z
.
f
;
trd2z
=
utr2z
.
f
;
trd3z
=
utr3z
.
f
;
}
else
{
trd1x
=
utr1x
.
d
;
trd2x
=
utr2x
.
d
;
trd3x
=
utr3x
.
d
;
trd1y
=
utr1y
.
d
;
trd2y
=
utr2y
.
d
;
trd3y
=
utr3y
.
d
;
trd1z
=
utr1z
.
d
;
trd2z
=
utr2z
.
d
;
trd3z
=
utr3z
.
d
;
}
dfix
=
(
(
-
sij
)
*
trd1x
+
(
ski
)
*
trd3x
)
*
(
vc
/
anorm3
);
dfjx
=
(
(
sij
)
*
trd1x
+
(
-
sjk
)
*
trd2x
)
*
(
vc
/
anorm3
);
dfkx
=
(
(
sjk
)
*
trd2x
+
(
-
ski
)
*
trd3x
)
*
(
vc
/
anorm3
);
dfiy
=
(
(
-
sij
)
*
trd1y
+
(
ski
)
*
trd3y
)
*
(
vc
/
anorm3
);
dfjy
=
(
(
sij
)
*
trd1y
+
(
-
sjk
)
*
trd2y
)
*
(
vc
/
anorm3
);
dfky
=
(
(
sjk
)
*
trd2y
+
(
-
ski
)
*
trd3y
)
*
(
vc
/
anorm3
);
dfiz
=
(
(
-
sij
)
*
trd1z
+
(
ski
)
*
trd3z
)
*
(
vc
/
anorm3
);
dfjz
=
(
(
sij
)
*
trd1z
+
(
-
sjk
)
*
trd2z
)
*
(
vc
/
anorm3
);
dfkz
=
(
(
sjk
)
*
trd2z
+
(
-
ski
)
*
trd3z
)
*
(
vc
/
anorm3
);
}
if
(
triplet_debug
)
force_debug_3t
(
xx
,
i
,
j
,
k
,
i
,
j
,
k
,
dfix
,
dfiy
,
dfiz
,
dfjx
,
dfjy
,
dfjz
,
dfkx
,
dfky
,
dfkz
);
if
(
three_body_forces
)
accumulate_forces_3
(
w3
);
}
if
(
T12
!=
0
)
{
//printf("T12 i,j,k = %d,%d,%d\n",i,j,k);
mcount
++
;
if
(
three_body_energies
&&
evflag
)
{
tr1
=
transtrace
(
T12
->
H1H2
,
T12
->
H1H2
);
double
dvir
=
(
2.0
*
(
dvir_ij
+
dvir_jk
)
*
splinepot
.
vd
+
splinepot
.
dvd
)
*
tr1
*
w3
/
anorm4
;
vir3v
=
vir3v
+
dvir
;
xvir3v
=
xvir3v
+
dvir
;
}
{
const
double
vd
=
splinepot
.
vd
;
tr_trace3
(
&
(
T12
->
H1H2
.
m
[
1
][
0
]),
&
(
T12
->
H1xH2
.
m
[
1
][
0
]),
&
utr1x
.
d
,
&
(
T12
->
H1yH2
.
m
[
1
][
0
]),
&
utr1y
.
d
,
&
(
T12
->
H1zH2
.
m
[
1
][
0
]),
&
utr1z
.
d
);
tr_trace3
(
&
(
T12
->
H1H2
.
m
[
1
][
0
]),
&
(
T12
->
H1H2x
.
m
[
1
][
0
]),
&
utr2x
.
d
,
&
(
T12
->
H1H2y
.
m
[
1
][
0
]),
&
utr2y
.
d
,
&
(
T12
->
H1H2z
.
m
[
1
][
0
]),
&
utr2z
.
d
);
if
(
linalg
.
single
)
{
trd1x
=
utr1x
.
f
;
trd2x
=
utr2x
.
f
;
trd1y
=
utr1y
.
f
;
trd2y
=
utr2y
.
f
;
trd1z
=
utr1z
.
f
;
trd2z
=
utr2z
.
f
;
}
else
{
trd1x
=
utr1x
.
d
;
trd2x
=
utr2x
.
d
;
trd1y
=
utr1y
.
d
;
trd2y
=
utr2y
.
d
;
trd1z
=
utr1z
.
d
;
trd2z
=
utr2z
.
d
;
}
dfix
=
2.0
*
(
-
sij
)
*
trd1x
*
(
vd
/
anorm4
);
dfkx
=
2.0
*
(
sjk
)
*
trd2x
*
(
vd
/
anorm4
);
dfjx
=
-
(
dfix
+
dfkx
);
dfiy
=
2.0
*
(
-
sij
)
*
trd1y
*
(
vd
/
anorm4
);
dfky
=
2.0
*
(
sjk
)
*
trd2y
*
(
vd
/
anorm4
);
dfjy
=
-
(
dfiy
+
dfky
);
dfiz
=
2.0
*
(
-
sij
)
*
trd1z
*
(
vd
/
anorm4
);
dfkz
=
2.0
*
(
sjk
)
*
trd2z
*
(
vd
/
anorm4
);
dfjz
=
-
(
dfiz
+
dfkz
);
}
if
(
triplet_debug
)
/* Compare forces to numerical derivatives */
force_debug_3v
(
xx
,
i
,
j
,
k
,
j
,
i
,
k
,
dfix
,
dfiy
,
dfiz
,
dfjx
,
dfjy
,
dfjz
,
dfkx
,
dfky
,
dfkz
);
if
(
three_body_forces
)
accumulate_forces_3
(
w3
);
}
if
(
T23
!=
0
)
{
//printf("T23 i,j,k = %d,%d,%d\n",i,j,k);
mcount
++
;
if
(
three_body_energies
&&
evflag
)
{
tr2
=
transtrace
(
T23
->
H1H2
,
T23
->
H1H2
);
double
dvir
=
(
2.0
*
(
dvir_jk
+
dvir_ki
)
*
splinepot
.
vd
+
splinepot
.
dvd
)
*
tr2
*
w3
/
anorm4
;
vir3v
=
vir3v
+
dvir
;
xvir3v
=
xvir3v
+
dvir
;
}
{
const
double
vd
=
splinepot
.
vd
;
tr_trace3
(
&
(
T23
->
H1H2
.
m
[
1
][
0
]),
&
(
T23
->
H1xH2
.
m
[
1
][
0
]),
&
utr1x
.
d
,
&
(
T23
->
H1yH2
.
m
[
1
][
0
]),
&
utr1y
.
d
,
&
(
T23
->
H1zH2
.
m
[
1
][
0
]),
&
utr1z
.
d
);
tr_trace3
(
&
(
T23
->
H1H2
.
m
[
1
][
0
]),
&
(
T23
->
H1H2x
.
m
[
1
][
0
]),
&
utr2x
.
d
,
&
(
T23
->
H1H2y
.
m
[
1
][
0
]),
&
utr2y
.
d
,
&
(
T23
->
H1H2z
.
m
[
1
][
0
]),
&
utr2z
.
d
);
if
(
linalg
.
single
)
{
trd1x
=
utr1x
.
f
;
trd2x
=
utr2x
.
f
;
trd1y
=
utr1y
.
f
;
trd2y
=
utr2y
.
f
;
trd1z
=
utr1z
.
f
;
trd2z
=
utr2z
.
f
;
}
else
{
trd1x
=
utr1x
.
d
;
trd2x
=
utr2x
.
d
;
trd1y
=
utr1y
.
d
;
trd2y
=
utr2y
.
d
;
trd1z
=
utr1z
.
d
;
trd2z
=
utr2z
.
d
;
}
dfix
=
2.0
*
(
ski
)
*
trd1x
*
(
vd
/
anorm4
);
dfjx
=
2.0
*
(
-
sjk
)
*
trd2x
*
(
vd
/
anorm4
);
dfkx
=
-
(
dfix
+
dfjx
);
dfiy
=
2.0
*
(
ski
)
*
trd1y
*
(
vd
/
anorm4
);
dfjy
=
2.0
*
(
-
sjk
)
*
trd2y
*
(
vd
/
anorm4
);
dfky
=
-
(
dfiy
+
dfjy
);
dfiz
=
2.0
*
(
ski
)
*
trd1z
*
(
vd
/
anorm4
);
dfjz
=
2.0
*
(
-
sjk
)
*
trd2z
*
(
vd
/
anorm4
);
dfkz
=
-
(
dfiz
+
dfjz
);
}
if
(
triplet_debug
)
/* Compare forces to numerical derivatives */
force_debug_3v
(
xx
,
i
,
j
,
k
,
k
,
i
,
j
,
dfix
,
dfiy
,
dfiz
,
dfjx
,
dfjy
,
dfjz
,
dfkx
,
dfky
,
dfkz
);
if
(
three_body_forces
)
accumulate_forces_3
(
w3
);
}
if
(
T31
!=
0
)
{
//printf("T31 i,j,k = %d,%d,%d\n",i,j,k);
mcount
++
;
if
(
three_body_energies
&&
evflag
)
{
tr3
=
transtrace
(
T31
->
H1H2
,
T31
->
H1H2
);
double
dvir
=
(
2.0
*
(
dvir_ki
+
dvir_ij
)
*
splinepot
.
vd
+
splinepot
.
dvd
)
*
tr3
*
w3
/
anorm4
;
vir3v
=
vir3v
+
dvir
;
xvir3v
=
xvir3v
+
dvir
;
}
{
const
double
vd
=
splinepot
.
vd
;
tr_trace3
(
&
(
T31
->
H1H2
.
m
[
1
][
0
]),
&
(
T31
->
H1xH2
.
m
[
1
][
0
]),
&
utr1x
.
d
,
&
(
T31
->
H1yH2
.
m
[
1
][
0
]),
&
utr1y
.
d
,
&
(
T31
->
H1zH2
.
m
[
1
][
0
]),
&
utr1z
.
d
);
tr_trace3
(
&
(
T31
->
H1H2
.
m
[
1
][
0
]),
&
(
T31
->
H1H2x
.
m
[
1
][
0
]),
&
utr2x
.
d
,
&
(
T31
->
H1H2y
.
m
[
1
][
0
]),
&
utr2y
.
d
,
&
(
T31
->
H1H2z
.
m
[
1
][
0
]),
&
utr2z
.
d
);
if
(
linalg
.
single
)
{
trd1x
=
utr1x
.
f
;
trd2x
=
utr2x
.
f
;
trd1y
=
utr1y
.
f
;
trd2y
=
utr2y
.
f
;
trd1z
=
utr1z
.
f
;
trd2z
=
utr2z
.
f
;
}
else
{
trd1x
=
utr1x
.
d
;
trd2x
=
utr2x
.
d
;
trd1y
=
utr1y
.
d
;
trd2y
=
utr2y
.
d
;
trd1z
=
utr1z
.
d
;
trd2z
=
utr2z
.
d
;
}
dfjx
=
2.0
*
(
sij
)
*
trd1x
*
(
vd
/
anorm4
);
dfkx
=
2.0
*
(
-
ski
)
*
trd2x
*
(
vd
/
anorm4
);
dfix
=
-
(
dfjx
+
dfkx
);
dfjy
=
2.0
*
(
sij
)
*
trd1y
*
(
vd
/
anorm4
);
dfky
=
2.0
*
(
-
ski
)
*
trd2y
*
(
vd
/
anorm4
);
dfiy
=
-
(
dfjy
+
dfky
);
dfjz
=
2.0
*
(
sij
)
*
trd1z
*
(
vd
/
anorm4
);
dfkz
=
2.0
*
(
-
ski
)
*
trd2z
*
(
vd
/
anorm4
);
dfiz
=
-
(
dfjz
+
dfkz
);
}
if
(
triplet_debug
)
/* Compare forces to numerical derivatives */
force_debug_3v
(
xx
,
i
,
j
,
k
,
i
,
j
,
k
,
dfix
,
dfiy
,
dfiz
,
dfjx
,
dfjy
,
dfjz
,
dfkx
,
dfky
,
dfkz
);
if
(
three_body_forces
)
accumulate_forces_3
(
w3
);
}
v33
=
tr0
/
anorm3
;
v43
=
(
tr1
+
tr2
+
tr3
)
/
anorm4
;
double
de_triplet
=
(
splinepot
.
vc
*
v33
+
splinepot
.
vd
*
v43
)
*
e_scale
*
w3
;
e_triplet
=
e_triplet
+
de_triplet
;
e_triplet_c
=
e_triplet_c
+
splinepot
.
vc
*
v33
*
e_scale
*
w3
;
c_t
++
;
//printf("xxxx %6d %6d %6d :: %20.10e\n",1,2,3,de_triplet);
if
(
evflag
)
{
double
drji
[
3
],
drki
[
3
];
double
fj
[
3
]
=
{
fjx
,
fjy
,
fjz
},
fk
[
3
]
=
{
fkx
,
fky
,
fkz
};
for
(
int
p
=
0
;
p
<
3
;
p
++
)
{
drji
[
p
]
=
xx
[
j
][
p
]
-
xx
[
i
][
p
];
drki
[
p
]
=
xx
[
k
][
p
]
-
xx
[
i
][
p
];
/* To fix stress-per-atom scaling. */
fj
[
p
]
*=
e_scale
;
fk
[
p
]
*=
e_scale
;
}
ev_tally3
(
i
,
j
,
k
,
de_triplet
,
0.0
,
fj
,
fk
,
drji
,
drki
);
if
(
volpres_flag
&&
vflag_atom
)
{
//virial[i] = virial[i] - (vir3v + vir3t) * rhoinv*e_scale;
double
dvir
=
-
(
xvir3v
+
xvir3t
)
*
rhoinv
*
e_scale
*
(
1.0
/
3.0
);
for
(
int
pp
=
0
;
pp
<
3
;
pp
++
)
{
vatom
[
i
][
pp
]
+=
dvir
;
vatom
[
j
][
pp
]
+=
dvir
;
vatom
[
k
][
pp
]
+=
dvir
;
}
}
fix
=
fix
+
fsave
[
0
][
0
];
fiy
=
fiy
+
fsave
[
0
][
1
];
fiz
=
fiz
+
fsave
[
0
][
2
];
fjx
=
fjx
+
fsave
[
1
][
0
];
fjy
=
fjy
+
fsave
[
1
][
1
];
fjz
=
fjz
+
fsave
[
1
][
2
];
fkx
=
fkx
+
fsave
[
2
][
0
];
fky
=
fky
+
fsave
[
2
][
1
];
fkz
=
fkz
+
fsave
[
2
][
2
];
}
tx1
=
gettime
();
ttriplet
+=
tx1
-
tx0
;
nttriplet
++
;
}
else
{
triplet_defer
=
1
;
}
if
(
four_body_energies
||
four_body_forces
)
if
(
j
<
i
)
{
/* Search for quadruplet */
tx0
=
gettime
();
mj
=
first
[
j
];
mk
=
first
[
k
];
/*
i is in both the j-list and the k-list, and i > k,
and lists are sorted, so the loop terminates.
*/
while
(
nlist_short
[
mj
]
<
i
&&
nlist_short
[
mk
]
<
i
)
{
if
(
mj
>=
first
[
j
+
1
]
||
mk
>=
first
[
k
+
1
])
{
printf
(
"Illegal quad...
\n
"
" j=%d first[j]=%d first[j+1]=%d mj=%d
\n
"
" k=%d first[k]=%d first[k+1]=%d mk=%d
\n
"
,
j
,
first
[
j
],
first
[
j
+
1
],
mj
,
k
,
first
[
k
],
first
[
k
+
1
],
mk
);
error
->
one
(
__FILE__
,
__LINE__
,
"Shit, brkoen quad loop"
);
}
if
(
nlist_short
[
mj
]
==
nlist_short
[
mk
])
{
/* Closed quadruplet */
m
=
nlist_short
[
mj
];
c_jm
=
c_km
=
1
;
const
int
sim
=
(
i
<
m
)
?
1
:
-
1
;
const
int
sjm
=
(
j
<
m
)
?
1
:
-
1
;
const
int
skm
=
(
k
<
m
)
?
1
:
-
1
;
w4
=
get_weight
(
triclinic
,
ss
[
i
],
ss
[
j
],
ss
[
k
],
ss
[
m
]);
if
(
w4
>
0.0
)
{
/* Alrady know ij,jk,ki,jm,km bonds. Look for im bond. */
mi
=
first
[
i
];
while
(
mi
<
first
[
i
+
1
]
&&
nlist_short
[
mi
]
<
m
)
mi
=
mi
+
1
;
if
(
mi
<
first
[
i
+
1
]
&&
nlist_short
[
mi
]
==
m
)
c_im
=
1
;
else
c_im
=
0
;
if
(
c_im
==
0
||
c_jk
==
0
||
(
c_jk
&&
c_im
&&
m
<
k
))
{
if
(
triplet_defer
)
{
dvir_ij
=
dvir_jk
=
dvir_ki
=
0.0
;
if
(
c_ij
&&
c_jk
)
T12
=
get_triplet
(
xx
,
j
,
i
,
k
,
&
bond_hash
,
&
T12work
,
&
dvir_ij
,
&
dvir_jk
);
if
(
c_ki
&&
c_jk
)
T23
=
get_triplet
(
xx
,
k
,
i
,
j
,
&
bond_hash
,
&
T23work
,
&
dvir_ki
,
&
dvir_jk
);
if
(
c_ij
&&
c_ki
)
T31
=
get_triplet
(
xx
,
i
,
j
,
k
,
&
bond_hash
,
&
T31work
,
&
dvir_ij
,
&
dvir_ki
);
triplet_defer
=
0
;
}
fmx
=
fmy
=
fmz
=
0.0
;
double
xvir4
=
0.0
;
if
(
evflag
)
{
fsave
[
0
][
0
]
=
fix
;
fsave
[
0
][
1
]
=
fiy
;
fsave
[
0
][
2
]
=
fiz
;
fsave
[
1
][
0
]
=
fjx
;
fsave
[
1
][
1
]
=
fjy
;
fsave
[
1
][
2
]
=
fjz
;
fsave
[
2
][
0
]
=
fkx
;
fsave
[
2
][
1
]
=
fky
;
fsave
[
2
][
2
]
=
fkz
;
fsave
[
3
][
0
]
=
fmx
;
fsave
[
3
][
1
]
=
fmy
;
fsave
[
3
][
2
]
=
fmz
;
fix
=
fiy
=
fiz
=
0.0
;
fjx
=
fjy
=
fjz
=
0.0
;
fkx
=
fky
=
fkz
=
0.0
;
fmx
=
fmy
=
fmz
=
0.0
;
}
tr1
=
tr2
=
tr3
=
0.0
;
dvir_im
=
dvir_jm
=
dvir_km
=
0.0
;
T45
=
T56
=
T64
=
0
;
if
(
T12
!=
0
&&
c_km
&&
c_im
)
T45
=
get_triplet
(
xx
,
m
,
i
,
k
,
&
bond_hash
,
&
T45work
,
&
dvir_im
,
&
dvir_km
);
if
(
T23
!=
0
&&
c_im
&&
c_jm
)
T56
=
get_triplet
(
xx
,
m
,
i
,
j
,
&
bond_hash
,
&
T56work
,
&
dvir_im
,
&
dvir_jm
);
if
(
T31
!=
0
&&
c_jm
&&
c_km
)
T64
=
get_triplet
(
xx
,
m
,
j
,
k
,
&
bond_hash
,
&
T64work
,
&
dvir_jm
,
&
dvir_km
);
if
(
T12
!=
0
&&
T45
!=
0
)
{
if
(
four_body_energies
&&
evflag
)
{
tr1
=
transtrace
(
T12
->
H1H2
,
T45
->
H1H2
);
double
dvir
=
(
(
dvir_ij
+
dvir_jk
+
dvir_im
+
dvir_km
)
*
splinepot
.
ve
+
splinepot
.
dve
)
*
tr1
*
w4
/
anorm4
;
vir4
=
vir4
+
dvir
;
xvir4
=
xvir4
+
dvir
;
}
qcount
++
;
{
const
double
ve
=
splinepot
.
ve
;
trd_update_4
(
T12
,
T45
);
dfix_update_4a
(
x
);
dfix_update_4a
(
y
);
dfix_update_4a
(
z
);
}
if
(
quad_debug
)
/* Compare forces to numerical derivatives */
force_debug_4
(
xx
,
i
,
j
,
k
,
m
,
i
,
j
,
k
,
m
,
dfix
,
dfiy
,
dfiz
,
dfjx
,
dfjy
,
dfjz
,
dfkx
,
dfky
,
dfkz
,
dfmx
,
dfmy
,
dfmz
);
if
(
four_body_forces
)
accumulate_forces_4
(
w4
);
}
if
(
T23
!=
0
&&
T56
!=
0
)
{
if
(
four_body_energies
&&
evflag
)
{
tr2
=
transtrace
(
T23
->
H1H2
,
T56
->
H1H2
);
double
dvir
=
(
(
dvir_ki
+
dvir_jk
+
dvir_im
+
dvir_jm
)
*
splinepot
.
ve
+
splinepot
.
dve
)
*
tr2
*
w4
/
anorm4
;
vir4
=
vir4
+
dvir
;
xvir4
=
xvir4
+
dvir
;
}
qcount
++
;
{
const
double
ve
=
splinepot
.
ve
;
trd_update_4
(
T23
,
T56
);
dfix_update_4b
(
x
);
dfix_update_4b
(
y
);
dfix_update_4b
(
z
);
}
if
(
quad_debug
)
/* Compare forces to numerical derivatives */
force_debug_4
(
xx
,
i
,
j
,
k
,
m
,
i
,
m
,
j
,
k
,
dfix
,
dfiy
,
dfiz
,
dfjx
,
dfjy
,
dfjz
,
dfkx
,
dfky
,
dfkz
,
dfmx
,
dfmy
,
dfmz
);
if
(
four_body_forces
)
accumulate_forces_4
(
w4
);
}
if
(
T31
!=
0
&&
T64
!=
0
)
{
if
(
four_body_energies
&&
evflag
)
{
tr3
=
transtrace
(
T31
->
H1H2
,
T64
->
H1H2
);
double
dvir
=
(
(
dvir_ki
+
dvir_ij
+
dvir_jm
+
dvir_km
)
*
splinepot
.
ve
+
splinepot
.
dve
)
*
tr3
*
w4
/
anorm4
;
vir4
=
vir4
+
dvir
;
xvir4
=
xvir4
+
dvir
;
}
qcount
++
;
{
const
double
ve
=
splinepot
.
ve
;
/* X */
trd_update_4
(
T31
,
T64
);
dfix_update_4c
(
x
);
dfix_update_4c
(
y
);
dfix_update_4c
(
z
);
}
if
(
quad_debug
)
/* Compare forces to numerical derivatives */
force_debug_4
(
xx
,
i
,
j
,
k
,
m
,
i
,
j
,
m
,
k
,
dfix
,
dfiy
,
dfiz
,
dfjx
,
dfjy
,
dfjz
,
dfkx
,
dfky
,
dfkz
,
dfmx
,
dfmy
,
dfmz
);
if
(
four_body_forces
)
accumulate_forces_4
(
w4
);
}
double
de_quad
=
splinepot
.
ve
*
(
tr1
+
tr2
+
tr3
)
/
anorm4
*
e_scale
*
w4
;
e_quad
=
e_quad
+
de_quad
;
if
((
T12
&&
T45
)
||
(
T23
&&
T56
)
||
(
T31
&&
T64
))
{
c_q
++
;
}
if
(
evflag
)
{
double
drim
[
3
],
drjm
[
3
],
drkm
[
3
];
double
fi
[
3
]
=
{
fix
,
fiy
,
fiz
};
double
fj
[
3
]
=
{
fjx
,
fjy
,
fjz
};
double
fk
[
3
]
=
{
fkx
,
fky
,
fkz
};
for
(
int
p
=
0
;
p
<
3
;
p
++
)
{
drim
[
p
]
=
xx
[
i
][
p
]
-
xx
[
m
][
p
];
drjm
[
p
]
=
xx
[
j
][
p
]
-
xx
[
m
][
p
];
drkm
[
p
]
=
xx
[
k
][
p
]
-
xx
[
m
][
p
];
fi
[
p
]
*=
e_scale
;
fj
[
p
]
*=
e_scale
;
fk
[
p
]
*=
e_scale
;
}
ev_tally4
(
i
,
j
,
k
,
m
,
de_quad
,
fi
,
fj
,
fk
,
drim
,
drjm
,
drkm
);
if
(
volpres_flag
&&
vflag_atom
)
{
//virial[i] = virial[i] - vir4 * rhoinv*e_scale;
double
dvir
=
-
xvir4
*
rhoinv
*
e_scale
*
(
1.0
/
4.0
);
for
(
int
pp
=
0
;
pp
<
3
;
pp
++
)
{
vatom
[
i
][
pp
]
+=
dvir
;
vatom
[
j
][
pp
]
+=
dvir
;
vatom
[
k
][
pp
]
+=
dvir
;
vatom
[
m
][
pp
]
+=
dvir
;
}
}
fix
=
fix
+
fsave
[
0
][
0
];
fiy
=
fiy
+
fsave
[
0
][
1
];
fiz
=
fiz
+
fsave
[
0
][
2
];
fjx
=
fjx
+
fsave
[
1
][
0
];
fjy
=
fjy
+
fsave
[
1
][
1
];
fjz
=
fjz
+
fsave
[
1
][
2
];
fkx
=
fkx
+
fsave
[
2
][
0
];
fky
=
fky
+
fsave
[
2
][
1
];
fkz
=
fkz
+
fsave
[
2
][
2
];
fmx
=
fmx
+
fsave
[
3
][
0
];
fmy
=
fmy
+
fsave
[
3
][
1
];
fmz
=
fmz
+
fsave
[
3
][
2
];
}
ff
[
m
][
0
]
+=
fmx
*
e_scale
;
ff
[
m
][
1
]
+=
fmy
*
e_scale
;
ff
[
m
][
2
]
+=
fmz
*
e_scale
;
}
}
mj
=
mj
+
1
;
mk
=
mk
+
1
;
}
else
if
(
nlist_short
[
mj
]
<
nlist_short
[
mk
])
{
mj
=
mj
+
1
;
}
else
{
mk
=
mk
+
1
;
}
}
tx1
=
gettime
();
tquad
+=
tx1
-
tx0
;
ntquad
++
;
ntquaditer
++
;
}
ff
[
k
][
0
]
+=
fkx
*
e_scale
;
ff
[
k
][
1
]
+=
fky
*
e_scale
;
ff
[
k
][
2
]
+=
fkz
*
e_scale
;
}
#undef transtrace
ff
[
j
][
0
]
+=
fjx
*
e_scale
;
ff
[
j
][
1
]
+=
fjy
*
e_scale
;
ff
[
j
][
2
]
+=
fjz
*
e_scale
;
}
ff
[
i
][
0
]
+=
fix
*
e_scale
;
ff
[
i
][
1
]
+=
fiy
*
e_scale
;
ff
[
i
][
2
]
+=
fiz
*
e_scale
;
if
(
single_energies
==
1
&&
i
<
nloc
)
{
const
double
evol0
=
splinepot
.
evol0
;
if
(
eflag_global
)
{
e_single
=
e_single
+
evol0
*
e_scale
;
eng_vdwl
=
eng_vdwl
+
evol0
*
e_scale
;
}
if
(
eflag_atom
)
eatom
[
i
]
=
eatom
[
i
]
+
evol0
*
e_scale
;
if
(
volpres_flag
&&
vflag_atom
)
{
for
(
int
pp
=
0
;
pp
<
3
;
pp
++
)
vatom
[
i
][
pp
]
=
vatom
[
i
][
pp
]
-
rhoinv
*
splinepot
.
devol0
*
e_scale
;
}
}
}
tx0
=
gettime
();
for
(
i
=
0
;
i
<
ntot
;
i
++
)
for
(
p
=
0
;
p
<
3
;
p
++
)
atom
->
f
[
i
][
p
]
=
atom
->
f
[
i
][
p
]
+
ff
[
i
][
p
];
memory
->
sfree
(
nlist_short
);
memory
->
sfree
(
first
);
if
(
ss
!=
xx
)
memory
->
sfree
(
ss
);
memory
->
sfree
(
ff
);
memory
->
sfree
(
xx
);
tx1
=
gettime
();
tmem
+=
tx1
-
tx0
;
ntmem
++
;
t1
=
gettime
(
1
);
//printf("compute_x: c_p = %d c_t = %d c_q = %d\n",c_p,c_t,c_q);
#ifdef TIMING_ON
if
(
comm
->
me
==
0
)
{
double
tsum
=
(
tmem
+
tsort
+
tpair
+
tlookup
+
ttriplet
+
tquad
);
double
nsum
=
(
ntmem
+
ntsort
+
ntpair
+
ntlookup
+
nttriplet
+
ntquad
);
//double adj = ((t1-t0)-tsum)/nsum;
/* Use adj = 6ns for RDTSC, and 58ns for gettimeofday,
on monkfish.llnl.gov, 2.4GHz Intel
Use adj = 35.945ns for RDTSC on uBGL (assumed rate set to 700MHz)
*/
double
adj
=
35.945e-9
;
double
memadj
=
tmem
-
adj
*
ntmem
,
sortadj
=
tsort
-
adj
*
ntsort
,
pairadj
=
tpair
-
adj
*
ntpair
,
lookupadj
=
tlookup
-
adj
*
ntlookup
,
tripletadj
=
ttriplet
-
adj
*
nttriplet
,
quadadj
=
tquad
-
adj
*
ntquad
,
make_b_adj
=
t_make_b
-
adj
*
n_make
,
make_t_adj
=
t_make_t
-
adj
*
n_make
,
make_b2_adj
=
t_make_b2
-
adj
*
n_make_b2
,
trace_adj
=
t_trace
-
adj
*
n_trace
;
printf
(
"mgpt engy = %10.3fms
\n
"
,(
t1
-
t0
)
*
1e3
);
printf
(
" mem = %10.3fms n = %8.0f adj = %10.3fms one = %10.3fns
\n
"
,
tmem
*
1e3
,
ntmem
,
memadj
*
1e3
,
memadj
/
ntmem
*
1e9
);
printf
(
" sort = %10.3fms n = %8.0f adj = %10.3fms one = %10.3fns
\n
"
,
tsort
*
1e3
,
ntsort
,
sortadj
*
1e3
,
sortadj
/
ntsort
*
1e9
);
printf
(
" pair = %10.3fms n = %8.0f adj = %10.3fms one = %10.3fns
\n
"
,
tpair
*
1e3
,
ntpair
,
pairadj
*
1e3
,
pairadj
/
ntpair
*
1e9
);
printf
(
" lookup = %10.3fms n = %8.0f adj = %10.3fms one = %10.3fns
\n
"
,
tlookup
*
1e3
,
ntlookup
,
lookupadj
*
1e3
,
lookupadj
/
ntlookup
*
1e9
);
printf
(
" triplet = %10.3fms n = %8.0f adj = %10.3fms one = %10.3fns
\n
"
,
ttriplet
*
1e3
,
nttriplet
,
tripletadj
*
1e3
,
tripletadj
/
nttriplet
*
1e9
);
printf
(
" quad = %10.3fms n = %8.0f adj = %10.3fms one = %10.3fns
\n
"
,
tquad
*
1e3
,
ntquaditer
,
quadadj
*
1e3
,
quadadj
/
ntquaditer
*
1e9
);
printf
(
" sum = %10.3fms adj = %10.3fms
\n
"
,
tsum
*
1e3
,(
tsum
-
adj
*
nsum
)
*
1e3
);
printf
(
"
\n
make_b = %10.3fms n = %8.0f adj = %10.3fms one = %10.3fns
\n
"
,
t_make_b
*
1e3
,
n_make
,
make_b_adj
*
1e3
,
make_b_adj
/
n_make
*
1e9
);
printf
(
" make_b2 = %10.3fms n = %8.0f adj = %10.3fms one = %10.3fns
\n
"
,
t_make_b2
*
1e3
,
n_make_b2
,
make_b2_adj
*
1e3
,
make_b2_adj
/
n_make_b2
*
1e9
);
printf
(
" make_t = %10.3fms n = %8.0f adj = %10.3fms one = %10.3fns
\n\n
"
,
t_make_t
*
1e3
,
n_make
,
make_t_adj
*
1e3
,
make_t_adj
/
n_make
*
1e9
);
printf
(
" trace = %10.3fms n = %8.0f adj = %10.3fms one = %10.3fns
\n\n
"
,
t_trace
*
1e3
,
n_trace
,
trace_adj
*
1e3
,
trace_adj
/
n_trace
*
1e9
);
printf
(
"mcount (transpose + trace for triplet) = %.0f , %.0f qcount = %.0f lmax = %d
\n
"
,
mcount
,
mcount2
,
qcount
,
lmax
);
printf
(
"nbc=%.0f tbl=%.3fms tbm=%.3fms one tbl=%.3fns one tbm=%.3fns
\n
"
,
nbc
,(
tbl
-
adj
*
nbc
)
*
1e3
,(
tbm
-
adj
*
nbc
)
*
1e3
,(
tbl
/
nbc
-
adj
)
*
1e9
,
(
tbm
/
nbc
-
adj
)
*
1e9
);
printf
(
"
\n\n
Forces:
\n
"
);
printf
(
"fix = %.3f fiy=%.3f fiz=%.3f
\n
"
,
fix
,
fiy
,
fiz
);
printf
(
"fjx = %.3f fjy=%.3f fjz=%.3f
\n
"
,
fjx
,
fjy
,
fjz
);
printf
(
"fkx = %.3f fky=%.3f fkz=%.3f
\n
"
,
fkx
,
fky
,
fkz
);
printf
(
"
\n
"
);
printf
(
"Bonds : nsearch=%d maxlen=%d avg.len=%.3f
\n
"
,
bond_hash
.
NSearch
(),
bond_hash
.
MaxLength
(),
bond_hash
.
NStep
()
/
(
double
)
bond_hash
.
NSearch
());
printf
(
"compute_x: c_p = %d c_t = %d c_q = %d
\n
"
,
c_p
,
c_t
,
c_q
);
printf
(
"@@ Total number of trace3 calls is %d, total number of make_triplet is %.1f
\n
"
,
ntr_calls
,
n_make
);
{
Hash
<
bond_data
,
Doublet
>::
Iterator
iter
=
bond_hash
.
begin
();
int
nitem
=
0
,
nhit
=
0
;
while
(
iter
!=
bond_hash
.
end
())
{
nitem
++
;
nhit
+=
iter
.
link
()
->
hits
;
iter
.
next
();
}
printf
(
"bond_hash hits: nitems=%d nhits=%d hits/item = %.3f
\n
"
,
nitem
,
nhit
,
nhit
/
(
double
)
nitem
);
}
}
#endif
if
(
volpres_flag
)
{
/*
Include contributions to the pressure due to derivatines
of the energy with respect to the potential input volume.
*/
/* The following lines have moved to beginning of functions,
since they are used in calculating per-atom virial contributions */
/*
double vtot = 1.0;
double ntot = atom->natoms;
for(i = 0; i<3; i++)
vtot = vtot * (domain->boxhi[i] - domain->boxlo[i]);
double rhoinv = vtot / ntot;
*/
if
(
single_energies
)
// Virial correction for self energy
for
(
i
=
0
;
i
<
3
;
i
++
)
{
//virial[i] = virial[i] + nloc*pot_input_vol*pvol0*e_scale;
virial
[
i
]
=
virial
[
i
]
-
nloc
*
rhoinv
*
splinepot
.
devol0
*
e_scale
;
}
if
(
pair_energies
)
// Virial correction for pair energy
for
(
i
=
0
;
i
<
3
;
i
++
)
virial
[
i
]
=
virial
[
i
]
+
rhoinv
*
e_scale
*
volvir2
;
if
(
three_body_energies
)
// Virial correction for three body enegries
for
(
i
=
0
;
i
<
3
;
i
++
)
{
//virial[i] = virial[i] - pot_input_vol*(e_triplet_c*pc + (e_triplet-e_triplet_c)*pd);
virial
[
i
]
=
virial
[
i
]
-
(
vir3v
+
vir3t
)
*
rhoinv
*
e_scale
;
}
if
(
four_body_energies
)
// Virial correction for four body enegries
for
(
i
=
0
;
i
<
3
;
i
++
)
{
//virial[i] = virial[i] - pot_input_vol*e_quad*pe;
virial
[
i
]
=
virial
[
i
]
-
vir4
*
rhoinv
*
e_scale
;
}
}
*
e_s
=
e_single
;
*
e_p
=
e_pair
;
*
e_t
=
e_triplet
;
*
e_q
=
e_quad
;
}
void
PairMGPT
::
compute
(
int
eflag
,
int
vflag
)
{
if
(
eflag
||
vflag
)
ev_setup
(
eflag
,
vflag
);
else
evflag
=
vflag_fdotr
=
eflag_global
=
vflag_global
=
eflag_atom
=
vflag_atom
=
0
;
int
newton_pair
=
force
->
newton_pair
;
double
e_s
,
e_p
,
e_t
,
e_q
;
//printf("newton_pair = %d, newton = %d, tag_enable = %d\n",force->newton_pair,force->newton,atom->tag_enable);
if
(
newton_pair
==
0
)
{
printf
(
"This is a problem. MGPT requires newton_pair flag to be on. Exiting...
\n
"
);
exit
(
1
);
}
if
(
atom
->
tag_enable
==
0
)
{
printf
(
"This is a problem. MGPT requires tag_enable flag to be on. Exiting...
\n
"
);
exit
(
1
);
}
compute_x
(
listfull
->
numneigh
,
listfull
->
firstneigh
,
&
e_s
,
&
e_p
,
&
e_t
,
&
e_q
,
evflag
,
newton_pair
);
if
(
0
)
{
// Stupid force calculation / verification
int
ii
,
nmax
=-
1
;
for
(
ii
=
0
;
ii
<
listfull
->
inum
+
listfull
->
gnum
;
ii
++
)
{
int
i
=
listfull
->
ilist
[
ii
];
if
(
i
>
nmax
)
nmax
=
i
;
}
nmax
++
;
double
*
ffwork
=
new
double
[
3
*
nmax
];
double
*
ffloc
=
new
double
[
3
*
listfull
->
inum
];
double
*
ffloc2
=
new
double
[
3
*
listfull
->
inum
];
double
**
ffptr
=
new
double
*
[
nmax
];
for
(
ii
=
0
;
ii
<
listfull
->
inum
+
listfull
->
gnum
;
ii
++
)
ffptr
[
ii
]
=
&
ffwork
[
3
*
ii
];
printf
(
"Computing boundary forces
\n
"
);
for
(
ii
=
0
;
ii
<
listfull
->
inum
;
ii
++
)
{
ffloc2
[
3
*
ii
]
=
0.0
;
ffloc2
[
3
*
ii
+
1
]
=
0.0
;
ffloc2
[
3
*
ii
+
2
]
=
0.0
;
int
i
=
listfull
->
ilist
[
ii
];
for
(
int
jj
=
0
;
jj
<
listfull
->
inum
+
listfull
->
gnum
;
jj
++
)
{
int
j
=
listfull
->
ilist
[
jj
];
if
(
atom
->
tag
[
i
]
==
atom
->
tag
[
j
])
for
(
int
p
=
0
;
p
<
3
;
p
++
)
ffloc2
[
3
*
ii
+
p
]
+=
atom
->
f
[
j
][
p
];
}
}
printf
(
"Starting main displacement force calculation
\n
"
);
for
(
ii
=
0
;
ii
<
listfull
->
inum
;
ii
++
)
{
int
i
=
listfull
->
ilist
[
ii
];
double
**
atom_f_save
=
atom
->
f
;
atom
->
f
=
ffptr
;
for
(
int
p
=
0
;
p
<
3
;
p
++
)
{
double
xsave
=
atom
->
x
[
i
][
p
];
const
double
delta
=
1e-3
;
atom
->
x
[
i
][
p
]
=
xsave
+
delta
;
for
(
int
jj
=
0
;
jj
<
3
*
nmax
;
jj
++
)
ffwork
[
jj
]
=
0.0
;
compute_x
(
listfull
->
numneigh
,
listfull
->
firstneigh
,
&
e_s
,
&
e_p
,
&
e_t
,
&
e_q
,
evflag
,
newton_pair
);
double
e1
=
e_s
+
e_p
+
e_t
+
e_q
;
atom
->
x
[
i
][
p
]
=
xsave
-
delta
;
for
(
int
jj
=
0
;
jj
<
3
*
nmax
;
jj
++
)
ffwork
[
jj
]
=
0.0
;
compute_x
(
listfull
->
numneigh
,
listfull
->
firstneigh
,
&
e_s
,
&
e_p
,
&
e_t
,
&
e_q
,
evflag
,
newton_pair
);
double
e2
=
e_s
+
e_p
+
e_t
+
e_q
;
ffloc
[
3
*
ii
+
p
]
=
-
(
e1
-
e2
)
/
(
2
*
delta
);
atom
->
x
[
i
][
p
]
=
xsave
;
}
atom
->
f
=
atom_f_save
;
printf
(
"Force on i=%4d:
\n
"
,
i
);
printf
(
" Position %20.10e %20.10e %20.10e
\n
"
,
atom
->
x
[
i
][
0
],
atom
->
x
[
i
][
1
],
atom
->
x
[
i
][
2
]);
printf
(
" Exact %20.10e %20.10e %20.10e
\n
"
,
atom
->
f
[
i
][
0
],
atom
->
f
[
i
][
1
],
atom
->
f
[
i
][
2
]);
printf
(
" Numerical %20.10e %20.10e %20.10e
\n
"
,
ffloc
[
3
*
ii
+
0
],
ffloc
[
3
*
ii
+
1
],
ffloc
[
3
*
ii
+
2
]);
printf
(
" Boundary %20.10e %20.10e %20.10e
\n
"
,
ffloc2
[
3
*
ii
+
0
],
ffloc2
[
3
*
ii
+
1
],
ffloc2
[
3
*
ii
+
2
]);
}
delete
[]
ffloc2
;
delete
[]
ffloc
;
delete
[]
ffptr
;
delete
[]
ffwork
;
}
if
(
0
)
{
printf
(
"
\n
Forces MGPT:
\n
"
);
const
int
iimax
=
(
listfull
->
inum
<
10
)
?
listfull
->
inum
:
10
;
for
(
int
ii
=
0
;
ii
<
iimax
;
ii
++
)
{
int
i
=
listfull
->
ilist
[
ii
];
printf
(
"%4d = %20.10e %20.10e %20.10e
\n
"
,
i
,
atom
->
f
[
i
][
0
],
atom
->
f
[
i
][
1
],
atom
->
f
[
i
][
2
]);
}
printf
(
"
\n\n
"
);
}
if
(
vflag_fdotr
)
{
//printf("##### Using virial_compute!!!\n");
virial_fdotr_compute
();
}
}
void
PairMGPT
::
allocate
()
{
allocated
=
1
;
int
n
=
atom
->
ntypes
;
memory
->
create
(
setflag
,
n
+
1
,
n
+
1
,
"pair:setflag"
);
for
(
int
i
=
0
;
i
<=
n
;
i
++
)
for
(
int
j
=
0
;
j
<=
n
;
j
++
)
setflag
[
i
][
j
]
=
0
;
memory
->
create
(
cutsq
,
n
+
1
,
n
+
1
,
"pair:cutsq"
);
memory
->
create
(
cutghost
,
n
+
1
,
n
+
1
,
"pair:cutsq"
);
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void
PairMGPT
::
settings
(
int
narg
,
char
**
arg
)
{
if
(
narg
!=
0
)
error
->
all
(
__FILE__
,
__LINE__
,
"Illegal pair_style command"
);
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void
PairMGPT
::
coeff
(
int
narg
,
char
**
arg
)
{
int
single_precision
=
0
;
if
(
narg
<
5
)
error
->
all
(
__FILE__
,
__LINE__
,
"Not enough arguments for mgpt (MGPT) pair coefficients."
);
if
(
!
allocated
)
allocate
();
// Make sure I,J args are * *
if
(
strcmp
(
arg
[
0
],
"*"
)
!=
0
||
strcmp
(
arg
[
1
],
"*"
)
!=
0
)
error
->
all
(
__FILE__
,
__LINE__
,
"Incorrect args for pair coefficients"
);
double
vol
;
if
(
sscanf
(
arg
[
4
],
"%lg"
,
&
vol
)
!=
1
||
vol
<=
0.0
)
error
->
all
(
__FILE__
,
__LINE__
,
"Invalid volume in mgpt (MGPT) pair coefficients."
);
volpres_flag
=
1
;
single_precision
=
0
;
/* Parse arguments */
{
int
volpres_tag
=
0
,
precision_tag
=
0
,
nbody_tag
=
0
;
int
iarg
=
5
;
while
(
iarg
<
narg
)
{
if
(
strcmp
(
arg
[
iarg
],
"volpress"
)
==
0
)
{
/* Volumetric pressure flag */
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Incorrect args for pair coefficients"
);
if
(
strcmp
(
arg
[
iarg
+
1
],
"yes"
)
==
0
)
volpres_flag
=
1
;
else
if
(
strcmp
(
arg
[
iarg
+
1
],
"no"
)
==
0
)
volpres_flag
=
0
;
else
{
char
line
[
1024
];
sprintf
(
line
,
"(In %s:%d) Invalid value for volumetric pressure argument.
\n
"
"It should be
\"
volpress yes
\"
or
\"
volpress no
\"
.
\n
"
"The value is
\"
%s
\"
.
\n
"
,
__FILE__
,
__LINE__
,
arg
[
iarg
+
1
]);
error
->
all
(
__FILE__
,
__LINE__
,
line
);
}
volpres_tag
=
1
;
iarg
+=
2
;
if
(
comm
->
me
==
0
)
printf
(
"* volpress: volpres_flag = %d [%s %s]
\n
"
,
volpres_flag
,
arg
[
iarg
-
2
],
arg
[
iarg
-
1
]);
}
else
if
(
strcmp
(
arg
[
iarg
],
"nbody"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Incorrect args for pair coefficients"
);
if
(
strspn
(
arg
[
iarg
+
1
],
"1234"
)
==
strlen
(
arg
[
iarg
+
1
]))
{
nbody_flag
=
0
;
for
(
int
i
=
0
;
i
<
4
;
i
++
)
if
(
strchr
(
arg
[
iarg
+
1
],
'1'
+
i
)
!=
NULL
)
{
nbody_flag
=
nbody_flag
+
(
1
<<
i
);
if
(
comm
->
me
==
0
)
printf
(
"Explicitly adding %d-tuple forces.
\n
"
,
i
+
1
);
}
}
else
{
char
line
[
1024
];
sprintf
(
line
,
"(In %s:%d) Invalid value for nbody flag.
\n
"
"It should be e.g.
\"
nbody=1234
\"
(for single, pair, triple, and quad forces/energiers)
\n
"
"For e.g. only pair and triple forces/energies, use
\"
nbody=23
\"
.
\n
"
"The default is
\"
nbody=1234
\"
.
\n
"
"The current value is
\"
%s
\"
.
\n
"
,
__FILE__
,
__LINE__
,
arg
[
iarg
+
1
]);
error
->
all
(
__FILE__
,
__LINE__
,
line
);
}
nbody_tag
=
1
;
iarg
+=
2
;
}
else
if
(
strcmp
(
arg
[
iarg
],
"precision"
)
==
0
)
{
if
(
iarg
+
2
>
narg
)
error
->
all
(
FLERR
,
"Incorrect args for pair coefficients"
);
if
(
strcmp
(
arg
[
iarg
+
1
],
"single"
)
==
0
)
single_precision
=
1
;
else
if
(
strcmp
(
arg
[
iarg
+
1
],
"double"
)
==
0
)
single_precision
=
0
;
else
{
char
line
[
1024
];
sprintf
(
line
,
"(In %s:%d) Invalid value for precision argument.
\n
"
"It should be
\"
precision single
\"
or
\"
precision double
\"
.
\n
"
"The value is
\"
%s
\"
.
\n
"
,
__FILE__
,
__LINE__
,
arg
[
iarg
+
1
]);
error
->
all
(
__FILE__
,
__LINE__
,
line
);
}
precision_tag
=
1
;
iarg
+=
2
;
if
(
comm
->
me
==
0
)
printf
(
"* precision: single_flag = %d [%s %s]
\n
"
,
single_precision
,
arg
[
iarg
-
2
],
arg
[
iarg
-
1
]);
}
else
{
char
line
[
1024
];
sprintf
(
line
,
"(In %s:%d) Invalid argument. Allowed arguments are:
\n
"
" volpress {yes|no} , default = yes
\n
"
" precision {single|double} , default = double
\n
"
" nbody {[1234,]*} , default = whichever terms potential require
\n
"
"The invalid argument is
\"
%s
\"
.
\n
"
,
__FILE__
,
__LINE__
,
arg
[
iarg
]);
error
->
all
(
__FILE__
,
__LINE__
,
line
);
}
}
if
(
comm
->
me
==
0
)
printf
(
"Volumetric pressure is %s.
\n
"
,
volpres_flag
?
"on"
:
"off"
);
if
(
comm
->
me
==
0
)
{
FILE
*
parmin_fp
=
force
->
open_potential
(
arg
[
2
]);
FILE
*
potin_fp
=
force
->
open_potential
(
arg
[
3
]);
if
(
parmin_fp
==
NULL
||
potin_fp
==
NULL
)
{
char
str
[
128
];
sprintf
(
str
,
"Cannot open MGPT potential files %s %s"
,
arg
[
2
],
arg
[
3
]);
error
->
one
(
FLERR
,
str
);
}
fclose
(
parmin_fp
);
fclose
(
potin_fp
);
splinepot
.
readpot
(
arg
[
2
],
arg
[
3
],
vol
);
printf
(
"evol0 = %.10e
\n
"
,
splinepot
.
evol0
);
/* Set up default and requested nbody forces to include */
{
int
nbody_default
=
(
1
<<
0
)
+
(
1
<<
1
)
+
(
1
<<
2
)
+
(
1
<<
3
);
if
(
splinepot
.
vd
==
0.0
&&
splinepot
.
dvd
==
0.0
)
nbody_default
-=
(
1
<<
2
);
// No 3-body contributions
if
(
splinepot
.
ve
==
0.0
&&
splinepot
.
dve
==
0.0
)
nbody_default
-=
(
1
<<
3
);
// No 4-body contributions
if
(
nbody_tag
==
0
)
nbody_flag
=
nbody_default
;
if
(
nbody_flag
!=
nbody_default
)
{
printf
(
"Warning: nbody=%d (suggested=%d) set to disregard multibody-forces in potential.
\n
"
,
nbody_flag
,
nbody_default
);
}
}
}
}
MPI_Bcast
(
&
nbody_flag
,
sizeof
(
nbody_flag
),
MPI_BYTE
,
0
,
world
);
/*
Broadcast structure to all processes. In receiving
processes, pointes will be screwed up. We allocate
memory, and then broadcast contents of arrays.
*/
MPI_Bcast
(
&
splinepot
,
sizeof
(
splinepot
),
MPI_BYTE
,
0
,
world
);
if
(
comm
->
me
!=
0
)
{
splinepot
.
vpair_spline
=
new
double
[
splinepot
.
nr
-
1
][
4
];
splinepot
.
dvpair_spline
=
new
double
[
splinepot
.
nr
-
1
][
4
];
}
MPI_Bcast
(
splinepot
.
vpair_spline
,
4
*
(
splinepot
.
nr
-
1
),
MPI_DOUBLE
,
0
,
world
);
MPI_Bcast
(
splinepot
.
dvpair_spline
,
4
*
(
splinepot
.
nr
-
1
),
MPI_DOUBLE
,
0
,
world
);
anorm3
=
splinepot
.
anorm3
;
anorm4
=
splinepot
.
anorm4
;
lmax
=
splinepot
.
lmax
;
lang
=
splinepot
.
lang
;
//ipot = splinepot.ipot;
for
(
int
i
=
0
;
i
<
(
int
)
(
sizeof
(
ddl
)
/
sizeof
(
double
));
i
++
)
ddl
[
i
]
=
splinepot
.
ddl
[
i
];
for
(
int
i
=
0
;
i
<
lmax
;
i
++
)
{
for
(
int
j
=
0
;
j
<
lmax
;
j
++
)
del0
.
m
[
i
+
1
][
j
+
1
]
=
0.0
;
del0
.
m
[
i
+
1
][
i
+
1
]
=
1.0
;
}
/* Set matrix param, cutoff, LAMMPS param */
Matrix
::
sz
=
lmax
;
rcrit
=
splinepot
.
rcrit
;
rmax
=
splinepot
.
rmax
;
cutoff
=
rmax
;
if
(
rcrit
>
rmax
)
cutoff
=
rcrit
;
// Set LAMMPS pair interaction flags.
for
(
int
i
=
1
;
i
<=
atom
->
ntypes
;
i
++
)
{
for
(
int
j
=
1
;
j
<=
atom
->
ntypes
;
j
++
)
{
setflag
[
i
][
j
]
=
1
;
cutsq
[
i
][
j
]
=
cutoff
;
cutghost
[
i
][
j
]
=
cutoff
;
}
}
// Set atomic mass.
for
(
int
i
=
1
;
i
<=
atom
->
ntypes
;
i
++
)
atom
->
set_mass
(
i
,
splinepot
.
mass
);
// Initialize linear algebra routines.
linalg
=
mgpt_linalg
(
lmax
,
single_precision
);
if
(
comm
->
me
==
0
)
printf
(
"%s"
,
linalg
.
msg
);
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void
PairMGPT
::
init_style
()
{
if
(
force
->
newton_pair
==
0
)
error
->
all
(
__FILE__
,
__LINE__
,
"Pair style mgpt requires newton pair on."
);
// Need full neighbor list.
int
irequest_full
=
neighbor
->
request
(
this
);
neighbor
->
requests
[
irequest_full
]
->
id
=
1
;
neighbor
->
requests
[
irequest_full
]
->
half
=
0
;
neighbor
->
requests
[
irequest_full
]
->
full
=
1
;
neighbor
->
requests
[
irequest_full
]
->
ghost
=
1
;
// Also need half neighbor list.
int
irequest_half
=
neighbor
->
request
(
this
);
neighbor
->
requests
[
irequest_half
]
->
id
=
2
;
neighbor
->
requests
[
irequest_half
]
->
half
=
0
;
neighbor
->
requests
[
irequest_half
]
->
half_from_full
=
1
;
neighbor
->
requests
[
irequest_half
]
->
otherlist
=
irequest_full
;
}
/* ----------------------------------------------------------------------
neighbor callback to inform pair style of neighbor list to use
half or full
------------------------------------------------------------------------- */
void
PairMGPT
::
init_list
(
int
id
,
NeighList
*
ptr
)
{
if
(
id
==
1
)
listfull
=
ptr
;
else
if
(
id
==
2
)
listhalf
=
ptr
;
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double
PairMGPT
::
init_one
(
int
i
,
int
j
)
{
return
cutoff
;
}
/************************************************************************
**** REIMPLEMENTATION OF FL AND HAMLTN WITH ANALYTICAL DERIVATIVES ****
************************************************************************/
/*
Reimplementation of bond length potential, including
derivatives with respect to x,y, and z.
*/
void
PairMGPT
::
fl_deriv_new
(
double
r
,
double
ri
,
double
xhat
,
double
yhat
,
double
zhat
,
double
&
fl_0
,
double
&
fl_x
,
double
&
fl_y
,
double
&
fl_z
,
double
&
fl_rp
,
double
&
fl_p1
,
double
&
fl_r0
,
double
&
fl_al
)
{
const
double
rp
=
splinepot
.
rp
,
p1
=
splinepot
.
p1
,
r0
=
splinepot
.
r00
,
al
=
splinepot
.
al
;
const
int
mode
=
splinepot
.
mode
;
const
double
pn
=
splinepot
.
pn
;
double
t
,
tx
,
ty
,
tz
,
t_rp_ti
,
t_p1_ti
;
double
s
;
/*
// Original code
double term;
double pn=1.0;
if (mode <= 4)
term = pow(rp/r, p1);
else
term = exp(-p1*(pow(r/rp, pn) - 1.0)/pn);
*/
double
rpi
=
1.0
/
rp
;
if
(
mode
<=
4
)
{
t
=
pow
(
rp
*
ri
,
p1
);
s
=
-
p1
*
t
*
ri
;
t_rp_ti
=
p1
*
rpi
;
t_p1_ti
=
log
(
rp
*
ri
);
}
else
{
if
(
pn
==
1.0
)
{
double
p1_rpi
=
-
p1
*
rpi
;
t
=
exp
(
p1
+
r
*
p1_rpi
);
s
=
p1_rpi
*
t
;
t_rp_ti
=
-
r
*
p1_rpi
*
rpi
;
t_p1_ti
=
1.0
-
r
*
rpi
;
}
else
{
double
pni
=
1.0
/
pn
;
double
rprpn
=
pow
(
r
*
rpi
,
pn
);
t
=
exp
(
-
p1
*
pni
*
(
rprpn
-
1.0
));
s
=
-
p1
*
rprpn
*
ri
*
t
;
t_rp_ti
=
p1
*
rprpn
*
rpi
;
t_p1_ti
=
pni
-
pni
*
rprpn
;
// -pni*(rprpn - 1.0);
}
}
tx
=
s
*
xhat
;
ty
=
s
*
yhat
;
tz
=
s
*
zhat
;
fl_rp
=
t_rp_ti
;
fl_p1
=
t_p1_ti
;
if
(
r
<=
r0
)
{
fl_0
=
t
;
fl_x
=
tx
;
fl_y
=
ty
;
fl_z
=
tz
;
fl_r0
=
0.0
;
fl_al
=
0.0
;
}
else
{
double
q
,
qx
,
qy
,
qz
,
exp_q
,
q_r0
,
q_al
;
double
r0i
,
u
;
r0i
=
1.0
/
r0
;
u
=
r
*
r0i
-
1.0
;
q
=
al
*
u
*
u
;
s
=
2
*
al
*
u
*
r0i
;
qx
=
s
*
xhat
;
qy
=
s
*
yhat
;
qz
=
s
*
zhat
;
q_r0
=
-
2.0
*
al
*
u
*
r
*
r0i
*
r0i
;
q_al
=
u
*
u
;
exp_q
=
exp
(
-
q
);
if
(
mode
<=
2
)
{
fl_0
=
exp_q
*
t
;
fl_x
=
exp_q
*
(
tx
-
t
*
qx
);
fl_y
=
exp_q
*
(
ty
-
t
*
qy
);
fl_z
=
exp_q
*
(
tz
-
t
*
qz
);
fl_r0
=
-
q_r0
;
fl_al
=
-
q_al
;
}
else
{
fl_0
=
exp_q
*
(
1.0
+
q
)
*
t
;
fl_x
=
exp_q
*
(
tx
+
q
*
(
tx
-
t
*
qx
));
fl_y
=
exp_q
*
(
ty
+
q
*
(
ty
-
t
*
qy
));
fl_z
=
exp_q
*
(
tz
+
q
*
(
tz
-
t
*
qz
));
fl_r0
=
-
q_r0
*
q
/
(
1.0
+
q
);
fl_al
=
-
q_al
*
q
/
(
1.0
+
q
);
}
}
}
/*
Macros to build elements of the bond matrix, and also
its derivatives with repsect to x,y, and z.
*/
#define MAKE_ELEMENT_5(i,j) \
do { \
const double dl0 = del0.m[i][j]; \
const double dl4 = gsl_##i * gsl_##j; \
const double dl4x = gsl_##i##x * gsl_##j + gsl_##i * gsl_##j##x; \
const double dl4y = gsl_##i##y * gsl_##j + gsl_##i * gsl_##j##y; \
const double dl4z = gsl_##i##z * gsl_##j + gsl_##i * gsl_##j##z; \
\
const double tmp = w4*dl4 + w2*dl2 + w0*dl0; \
const double tmpx = w4*dl4x + w2*dl2x; \
const double tmpy = w4*dl4y + w2*dl2y; \
const double tmpz = w4*dl4z + w2*dl2z; \
const double tmpsum = tmpx*x + tmpy*y + tmpz*z; \
M [i][j] = fl *tmp; \
Mx[i][j] = fl_x*tmp + fl_ri*(tmpx - x*tmpsum); \
My[i][j] = fl_y*tmp + fl_ri*(tmpy - y*tmpsum); \
Mz[i][j] = fl_z*tmp + fl_ri*(tmpz - z*tmpsum); \
if (i != j) { \
M[j][i] = M[i][j]; Mx[j][i] = Mx[i][j]; \
My[j][i] = My[i][j]; Mz[j][i] = Mz[i][j]; \
} \
} while(0)
#define MAKE_ELEMENT_7(i,j) \
do { \
const double dl0 = del0.m[i][j]; \
const double dl6 = gsl_##i * gsl_##j; \
const double dl6x = gsl_##i##x * gsl_##j + gsl_##i * gsl_##j##x; \
const double dl6y = gsl_##i##y * gsl_##j + gsl_##i * gsl_##j##y; \
const double dl6z = gsl_##i##z * gsl_##j + gsl_##i * gsl_##j##z; \
\
const double tmp = w6*dl6 + w4*dl4 + w2*dl2 + w0*dl0; \
const double tmpx = w6*dl6x + w4*dl4x + w2*dl2x; \
const double tmpy = w6*dl6y + w4*dl4y + w2*dl2y; \
const double tmpz = w6*dl6z + w4*dl4z + w2*dl2z; \
const double tmpsum = tmpx*x + tmpy*y + tmpz*z; \
M [i][j] = fl *tmp; \
Mx[i][j] = fl_x*tmp + fl_ri*(tmpx - x*tmpsum); \
My[i][j] = fl_y*tmp + fl_ri*(tmpy - y*tmpsum); \
Mz[i][j] = fl_z*tmp + fl_ri*(tmpz - z*tmpsum); \
if (i != j) { \
M[j][i] = M[i][j]; Mx[j][i] = Mx[i][j]; \
My[j][i] = My[i][j]; Mz[j][i] = Mz[i][j]; \
} \
} while(0)
/* End of bond matrix macros */
/*
Construction of bond matrix, and its derivatives
with respect to the coordinates
*/
void
PairMGPT
::
hamltn_5_raw
(
const
double
xin
,
const
double
yin
,
const
double
zin
,
double
M
[
8
][
8
],
double
Mx
[
8
][
8
],
double
My
[
8
][
8
],
double
Mz
[
8
][
8
],
double
*
fl_deriv_sum_p
)
{
const
double
r
=
sqrt
(
xin
*
xin
+
yin
*
yin
+
zin
*
zin
),
ri
=
1.0
/
r
;
const
double
x
=
xin
*
ri
,
y
=
yin
*
ri
,
z
=
zin
*
ri
;
// d-d
// call delndd(x,y,z)
const
double
x2
=
x
*
x
,
y2
=
y
*
y
,
z2
=
z
*
z
;
const
double
xy
=
x
*
y
,
xz
=
x
*
z
,
yz
=
y
*
z
;
const
double
sr3
=
sqrt
(
3.0
),
sr3i
=
1.0
/
sr3
;
const
double
frac_1_3
=
1.0
/
3.0
,
frac_2_3
=
2.0
/
3.0
,
frac_4_3
=
4.0
/
3.0
;
const
double
ddl_1
=
ddl
[
1
],
ddl_2
=
ddl
[
2
],
ddl_3
=
ddl
[
3
];
const
double
w4
=
ddl_1
-
frac_4_3
*
ddl_2
+
frac_1_3
*
ddl_3
;
const
double
w2
=
ddl_2
-
ddl_3
;
const
double
w0
=
ddl_2
;
//del4
double
gsl_1
,
gsl_2
,
gsl_3
,
gsl_4
,
gsl_5
;
double
gsl_1x
,
gsl_2x
,
gsl_3x
,
gsl_4x
,
gsl_5x
;
double
gsl_1y
,
gsl_2y
,
gsl_3y
,
gsl_4y
,
gsl_5y
;
double
gsl_1z
,
gsl_2z
,
gsl_3z
,
gsl_4z
,
gsl_5z
;
double
dl2
,
dl2x
,
dl2y
,
dl2z
;
double
fl
,
fl_x
,
fl_y
,
fl_z
,
fl_ri
;
double
fl_rp
,
fl_p1
,
fl_r0
,
fl_al
;
gsl_1
=
0.5
*
(
3.0
*
z2
-
1.0
);
gsl_1x
=
0.0
;
gsl_1y
=
0.0
;
gsl_1z
=
3.0
*
z
;
gsl_2
=
sr3
*
xz
;
gsl_2x
=
sr3
*
z
;
gsl_2y
=
0.0
;
gsl_2z
=
sr3
*
x
;
gsl_3
=
sr3
*
yz
;
gsl_3x
=
0.0
;
gsl_3y
=
sr3
*
z
;
gsl_3z
=
sr3
*
y
;
gsl_4
=
sr3
*
(
x2
-
y2
)
*
0.5
;
gsl_4x
=
sr3
*
x
;
gsl_4y
=
-
sr3
*
y
;
gsl_4z
=
0.0
;
gsl_5
=
sr3
*
xy
;
gsl_5x
=
sr3
*
y
;
gsl_5y
=
sr3
*
x
;
gsl_5z
=
0.0
;
// Compute bond length potential
fl_deriv_new
(
r
,
ri
,
x
,
y
,
z
,
fl
,
fl_x
,
fl_y
,
fl_z
,
fl_rp
,
fl_p1
,
fl_r0
,
fl_al
);
fl_ri
=
fl
*
ri
;
*
fl_deriv_sum_p
=
fl_rp
*
splinepot
.
drp
+
fl_p1
*
splinepot
.
dp1
+
fl_r0
*
splinepot
.
dr00
+
fl_al
*
splinepot
.
dal
;
// del2
//del2.m[1][1] = z2 - 2.0/3.0;
dl2
=
z2
-
frac_2_3
;
dl2x
=
0.0
;
dl2y
=
0.0
;
dl2z
=
2
*
z
;
MAKE_ELEMENT_5
(
1
,
1
);
//del2.m[1][2] = xz/sr3;
dl2
=
xz
*
sr3i
;
dl2x
=
z
*
sr3i
;
dl2y
=
0.0
;
dl2z
=
x
*
sr3i
;
MAKE_ELEMENT_5
(
1
,
2
);
//del2.m[1][3] = yz/sr3;
dl2
=
yz
*
sr3i
;
dl2x
=
0.0
;
dl2y
=
z
*
sr3i
;
dl2z
=
y
*
sr3i
;
MAKE_ELEMENT_5
(
1
,
3
);
//del2.m[1][4] = -(x2 - y2)*sr3i;
dl2
=
-
(
x2
-
y2
)
*
sr3i
;
dl2x
=
-
2.0
*
sr3i
*
x
;
dl2y
=
2.0
*
sr3i
*
y
;
dl2z
=
0.0
;
MAKE_ELEMENT_5
(
1
,
4
);
//del2.m[1][5] = -2.0*xy*sr3i;
dl2
=
-
2.0
*
xy
*
sr3i
;
dl2x
=
-
2.0
*
y
*
sr3i
;
dl2y
=
-
2.0
*
x
*
sr3i
;
dl2z
=
0.0
;
MAKE_ELEMENT_5
(
1
,
5
);
//del2.m[2][2] = -y2;
dl2
=
-
y2
;
dl2x
=
0.0
;
dl2y
=
-
2.0
*
y
;
dl2z
=
0.0
;
MAKE_ELEMENT_5
(
2
,
2
);
//del2.m[2][3] = xy;
dl2
=
xy
;
dl2x
=
y
;
dl2y
=
x
;
dl2z
=
0.0
;
MAKE_ELEMENT_5
(
2
,
3
);
//del2.m[2][4] = xz;
dl2
=
xz
;
dl2x
=
z
;
dl2y
=
0.0
;
dl2z
=
x
;
MAKE_ELEMENT_5
(
2
,
4
);
//del2.m[2][5] = yz;
dl2
=
yz
;
dl2x
=
0.0
;
dl2y
=
z
;
dl2z
=
y
;
MAKE_ELEMENT_5
(
2
,
5
);
//del2.m[3][3] = -x2;
dl2
=
-
x2
;
dl2x
=
-
2.0
*
x
;
dl2y
=
0.0
;
dl2z
=
0.0
;
MAKE_ELEMENT_5
(
3
,
3
);
//del2.m[3][4] = -yz;
dl2
=
-
yz
;
dl2x
=
0.0
;
dl2y
=
-
z
;
dl2z
=
-
y
;
MAKE_ELEMENT_5
(
3
,
4
);
//del2.m[3][5] = xz;
dl2
=
xz
;
dl2x
=
z
;
dl2y
=
0.0
;
dl2z
=
x
;
MAKE_ELEMENT_5
(
3
,
5
);
//del2.m[4][4] = -z2;
dl2
=
-
z2
;
dl2x
=
0.0
;
dl2y
=
0.0
;
dl2z
=
-
2.0
*
z
;
MAKE_ELEMENT_5
(
4
,
4
);
//del2.m[4][5] = 0.0;
dl2
=
0.0
;
dl2x
=
0.0
;
dl2y
=
0.0
;
dl2z
=
0.0
;
MAKE_ELEMENT_5
(
4
,
5
);
//del2.m[5][5] = -z2;
dl2
=
-
z2
;
dl2x
=
0.0
;
dl2y
=
0.0
;
dl2z
=
-
2.0
*
z
;
MAKE_ELEMENT_5
(
5
,
5
);
}
void
PairMGPT
::
hamltn_7_raw
(
const
double
xin
,
const
double
yin
,
const
double
zin
,
double
M
[
8
][
8
],
double
Mx
[
8
][
8
],
double
My
[
8
][
8
],
double
Mz
[
8
][
8
],
double
*
fl_deriv_sum_p
)
{
const
double
r
=
sqrt
(
xin
*
xin
+
yin
*
yin
+
zin
*
zin
),
ri
=
1.0
/
r
;
const
double
x
=
xin
*
ri
,
y
=
yin
*
ri
,
z
=
zin
*
ri
;
// d-d
// call delndd(x,y,z)
const
double
x2
=
x
*
x
,
y2
=
y
*
y
,
z2
=
z
*
z
;
const
double
xy
=
x
*
y
,
xz
=
x
*
z
,
yz
=
y
*
z
;
const
double
x4
=
x2
*
x2
,
y4
=
y2
*
y2
;
//const double sr3 = sqrt(3.0);//,sr3i = 1.0/sr3;
//const double frac_1_3 = 1.0/3.0,frac_2_3 = 2.0/3.0,frac_4_3 = 4.0/3.0;
const
double
sr01
=
sqrt
(
0.1
);
const
double
sr015
=
sqrt
(
0.15
);
const
double
sr024
=
sqrt
(
0.24
);
const
double
sr0375
=
sqrt
(
0.375
);
const
double
sr06
=
sqrt
(
0.6
);
const
double
sr0625
=
sqrt
(
0.625
);
const
double
sr09
=
sqrt
(
0.9
);
const
double
sr15
=
sqrt
(
1.5
);
const
double
sr24
=
sqrt
(
2.4
);
const
double
sr36
=
sqrt
(
3.6
);
const
double
sr375
=
sqrt
(
3.75
);
const
double
sr96
=
sqrt
(
9.6
);
const
double
sr150
=
sqrt
(
15.0
);
const
double
ddl_1
=
ddl
[
1
],
ddl_2
=
ddl
[
2
],
ddl_3
=
ddl
[
3
],
ddl_4
=
ddl
[
4
];
const
double
w6
=
ddl_1
-
1.5
*
ddl_2
+
0.6
*
ddl_3
-
0.1
*
ddl_4
;
const
double
w4
=
0.625
*
ddl_2
-
ddl_3
+
0.375
*
ddl_4
;
const
double
w2
=
0.625
*
(
ddl_2
-
ddl_4
);
const
double
w0
=
0.625
*
ddl_2
+
0.375
*
ddl_4
;
//del6
double
gsl_1
,
gsl_2
,
gsl_3
,
gsl_4
,
gsl_5
,
gsl_6
,
gsl_7
;
double
gsl_1x
,
gsl_2x
,
gsl_3x
,
gsl_4x
,
gsl_5x
,
gsl_6x
,
gsl_7x
;
double
gsl_1y
,
gsl_2y
,
gsl_3y
,
gsl_4y
,
gsl_5y
,
gsl_6y
,
gsl_7y
;
double
gsl_1z
,
gsl_2z
,
gsl_3z
,
gsl_4z
,
gsl_5z
,
gsl_6z
,
gsl_7z
;
double
dl2
,
dl2x
,
dl2y
,
dl2z
;
double
dl4
,
dl4x
,
dl4y
,
dl4z
;
double
t1
;
double
fl
,
fl_x
,
fl_y
,
fl_z
,
fl_ri
;
double
fl_rp
,
fl_p1
,
fl_r0
,
fl_al
;
//gslf[1] = 0.5*(5.0*n2 - 3.0)*n;
gsl_1
=
0.5
*
(
5.0
*
z2
-
3.0
)
*
z
;
gsl_1x
=
0.0
;
gsl_1y
=
0.0
;
gsl_1z
=
7.5
*
z2
-
1.5
;
//gslf[2] = sr0375*(5.0*n2 - 1.0)*l;
gsl_2
=
sr0375
*
(
5.0
*
z2
-
1.0
)
*
x
;
gsl_2x
=
sr0375
*
(
5.0
*
z2
-
1.0
);
gsl_2y
=
0.0
;
gsl_2z
=
sr0375
*
10.0
*
xz
;
//gslf[3] = sr0375*(5.0*n2 - 1.0)*m;
gsl_3
=
sr0375
*
(
5.0
*
z2
-
1.0
)
*
y
;
gsl_3x
=
0.0
;
gsl_3y
=
sr0375
*
(
5.0
*
z2
-
1.0
);
gsl_3z
=
sr0375
*
10.0
*
yz
;
//gslf[4] = sr375*(l2 - m2)*n;
gsl_4
=
sr375
*
(
x2
-
y2
)
*
z
;
gsl_4x
=
2.0
*
sr375
*
xz
;
gsl_4y
=
-
2.0
*
sr375
*
yz
;
gsl_4z
=
sr375
*
(
x2
-
y2
);
//gslf[5] = sr150*lm*n;
gsl_5
=
sr150
*
xy
*
z
;
gsl_5x
=
sr150
*
yz
;
gsl_5y
=
sr150
*
xz
;
gsl_5z
=
sr150
*
xy
;
//gslf[6] = sr0625*(l2 - 3.0*m2)*l;
gsl_6
=
sr0625
*
(
x2
-
3.0
*
y2
)
*
x
;
gsl_6x
=
3.0
*
sr0625
*
(
x2
-
y2
);
gsl_6y
=
-
6.0
*
sr0625
*
xy
;
gsl_6z
=
0.0
;
//gslf[7] = sr0625*(3.0*l2 - m2)*m;
gsl_7
=
sr0625
*
(
3.0
*
x2
-
y2
)
*
y
;
gsl_7x
=
6.0
*
sr0625
*
xy
;
gsl_7y
=
3.0
*
sr0625
*
(
x2
-
y2
);
gsl_7z
=
0.0
;
// Compute bond length potential
fl_deriv_new
(
r
,
ri
,
x
,
y
,
z
,
fl
,
fl_x
,
fl_y
,
fl_z
,
fl_rp
,
fl_p1
,
fl_r0
,
fl_al
);
fl_ri
=
fl
*
ri
;
*
fl_deriv_sum_p
=
fl_rp
*
splinepot
.
drp
+
fl_p1
*
splinepot
.
dp1
+
fl_r0
*
splinepot
.
dr00
+
fl_al
*
splinepot
.
dal
;
// del2f
//del2f.m[1][1] = 0.4*(3.0*n2 - 1.0);
dl2
=
0.4
*
(
3.0
*
z2
-
1.0
);
dl2x
=
0.0
;
dl2y
=
0.0
;
dl2z
=
2.4
*
z
;
//del4f.m[1][1] = 0.60*(5.0*n2 - 4.0)*n2;
dl4
=
0.60
*
(
5.0
*
z2
-
4.0
)
*
z2
;
dl4x
=
0.0
;
dl4y
=
0.0
;
dl4z
=
0.60
*
(
20.0
*
z2
-
8.0
)
*
z
;
MAKE_ELEMENT_7
(
1
,
1
);
//del2f.m[1][2] = sr024*ln;
dl2
=
sr024
*
xz
;
dl2x
=
sr024
*
z
;
dl2y
=
0.0
;
dl2z
=
sr024
*
x
;
//del4f.m[1][2] = sr024*(5.0*n2 - 2.0)*ln;
dl4
=
sr024
*
(
5.0
*
z2
-
2.0
)
*
xz
;
dl4x
=
sr024
*
(
5.0
*
z2
-
2.0
)
*
z
;
dl4y
=
0.0
;
dl4z
=
sr024
*
(
15.0
*
z2
-
2.0
)
*
x
;
MAKE_ELEMENT_7
(
1
,
2
);
//del2f.m[1][3] = sr024*mn;
dl2
=
sr024
*
yz
;
dl2x
=
0.0
;
dl2y
=
sr024
*
z
;
dl2z
=
sr024
*
y
;
//del4f.m[1][3] = sr024*(5.0*n2 - 2.0)*mn;
dl4
=
sr024
*
(
5.0
*
z2
-
2.0
)
*
yz
;
dl4x
=
0.0
;
dl4y
=
sr024
*
(
5.0
*
z2
-
2.0
)
*
z
;
dl4z
=
sr024
*
(
15.0
*
z2
-
2.0
)
*
y
;
MAKE_ELEMENT_7
(
1
,
3
);
//del2f.m[1][4] = -sr06*(l2 - m2);
dl2
=
-
sr06
*
(
x2
-
y2
);
dl2x
=
-
2.0
*
sr06
*
x
;
dl2y
=
2.0
*
sr06
*
y
;
dl2z
=
0.0
;
//del4f.m[1][4] = -sr06*(l2 - m2)*n2;
dl4
=
-
sr06
*
(
x2
-
y2
)
*
z2
;
dl4x
=
-
2.0
*
sr06
*
x
*
z2
;
dl4y
=
2.0
*
sr06
*
y
*
z2
;
dl4z
=
-
2.0
*
sr06
*
(
x2
-
y2
)
*
z
;
MAKE_ELEMENT_7
(
1
,
4
);
//del2f.m[1][5] = -sr24*lm;
dl2
=
-
sr24
*
xy
;
dl2x
=
-
sr24
*
y
;
dl2y
=
-
sr24
*
x
;
dl2z
=
0.0
;
//del4f.m[1][5] = -sr24*lm*n2;
dl4
=
-
sr24
*
xy
*
z2
;
dl4x
=
-
sr24
*
y
*
z2
;
dl4y
=
-
sr24
*
x
*
z2
;
dl4z
=
-
2.0
*
sr24
*
xy
*
z
;
MAKE_ELEMENT_7
(
1
,
5
);
//del2f.m[1][6] = 0.0;
dl2
=
0.0
;
dl2x
=
0.0
;
dl2y
=
0.0
;
dl2z
=
0.0
;
//del4f.m[1][6] = sr36*(3.0*m2 - l2)*ln;
dl4
=
sr36
*
(
3.0
*
y2
-
x2
)
*
xz
;
dl4x
=
3.0
*
sr36
*
(
y2
-
x2
)
*
z
;
dl4y
=
6.0
*
sr36
*
y
*
xz
;
dl4z
=
sr36
*
(
3.0
*
y2
-
x2
)
*
x
;
MAKE_ELEMENT_7
(
1
,
6
);
//del2f.m[1][7] = 0.0;
dl2
=
0.0
;
dl2x
=
0.0
;
dl2y
=
0.0
;
dl2z
=
0.0
;
//del4f.m[1][7] = -sr36*(3.0*l2 - m2)*mn;
dl4
=
-
sr36
*
(
3.0
*
x2
-
y2
)
*
yz
;
dl4x
=
-
6.0
*
sr36
*
x
*
yz
;
dl4y
=
-
3.0
*
sr36
*
(
x2
-
y2
)
*
z
;
dl4z
=
-
sr36
*
(
3.0
*
x2
-
y2
)
*
y
;
MAKE_ELEMENT_7
(
1
,
7
);
//del2f.m[2][2] = 0.3*(1.0 - 4.0*m2 + n2);
dl2
=
0.3
-
1.2
*
y2
+
0.3
*
z2
;
dl2x
=
0.0
;
dl2y
=
-
2.4
*
y
;
dl2z
=
0.6
*
z
;
//del4f.m[2][2] = -0.4*l*l - 2.5*(m2 - 0.6*l2)*n2;
dl4
=
-
0.4
*
x2
-
2.5
*
(
y2
-
0.6
*
x2
)
*
z2
;
dl4x
=
-
0.8
*
x
+
3.0
*
x
*
z2
;
dl4y
=
-
5.0
*
y
*
z2
;
dl4z
=
-
5.0
*
(
y2
-
0.6
*
x2
)
*
z
;
MAKE_ELEMENT_7
(
2
,
2
);
//del2f.m[2][3] = 1.2*lm;
dl2
=
1.2
*
xy
;
dl2x
=
1.2
*
y
;
dl2y
=
1.2
*
x
;
dl2z
=
0.0
;
//del4f.m[2][3] = 0.4*(10.0*n2 - 1.0)*lm;
dl4
=
(
4.0
*
z2
-
0.4
)
*
xy
;
dl4x
=
(
4.0
*
z2
-
0.4
)
*
y
;
dl4y
=
(
4.0
*
z2
-
0.4
)
*
x
;
dl4z
=
8.0
*
z
*
xy
;
MAKE_ELEMENT_7
(
2
,
3
);
//del2f.m[2][4] = sr09*ln;
dl2
=
sr09
*
xz
;
dl2x
=
sr09
*
z
;
dl2y
=
0.0
;
dl2z
=
sr09
*
x
;
//del4f.m[2][4] = sr01*(6.0*n2 - 8.0*m2 - 1.0)*ln;
dl4
=
sr01
*
(
6.0
*
z2
-
8.0
*
y2
-
1.0
)
*
xz
;
dl4x
=
sr01
*
(
6.0
*
z2
-
8.0
*
y2
-
1.0
)
*
z
;
dl4y
=
-
16.0
*
sr01
*
y
*
xz
;
dl4z
=
sr01
*
(
18.0
*
z2
-
8.0
*
y2
-
1.0
)
*
x
;
MAKE_ELEMENT_7
(
2
,
4
);
//del2f.m[2][5] = sr09*mn;
dl2
=
sr09
*
yz
;
dl2x
=
0.0
;
dl2y
=
sr09
*
z
;
dl2z
=
sr09
*
y
;
//del4f.m[2][5] = sr01*(2.0*n2 - 8.0*m2 + 3.0)*mn;
dl4
=
sr01
*
(
2.0
*
z2
-
8.0
*
y2
+
3.0
)
*
yz
;
dl4x
=
0.0
;
dl4y
=
sr01
*
(
2.0
*
z2
-
24.0
*
y2
+
3.0
)
*
z
;
dl4z
=
sr01
*
(
6.0
*
z2
-
8.0
*
y2
+
3.0
)
*
y
;
MAKE_ELEMENT_7
(
2
,
5
);
//del2f.m[2][6] = -sr015*(l2 - m2);
dl2
=
-
sr015
*
(
x2
-
y2
);
dl2x
=
-
2.0
*
sr015
*
x
;
dl2y
=
2.0
*
sr015
*
y
;
dl2z
=
0.0
;
//del4f.m[2][6] = sr375*(l2 - m2 - 1.4*l4 + 1.2*l2*m2 + m4);
dl4
=
sr375
*
(
x2
-
y2
-
1.4
*
x4
+
1.2
*
x2
*
y2
+
y4
);
dl4x
=
sr375
*
(
2.0
-
5.6
*
x2
+
2.4
*
y2
)
*
x
;
dl4y
=
sr375
*
(
-
2.0
+
2.4
*
x2
+
4.0
*
y2
)
*
y
;
dl4z
=
0.0
;
MAKE_ELEMENT_7
(
2
,
6
);
//del2f.m[2][7] = -sr06*lm;
dl2
=
-
sr06
*
xy
;
dl2x
=
-
sr06
*
y
;
dl2y
=
-
sr06
*
x
;
dl2z
=
0.0
;
//del4f.m[2][7] = sr96*(n2 - l2 + 0.25)*lm;
dl4
=
sr96
*
(
z2
-
x2
+
0.25
)
*
xy
;
dl4x
=
sr96
*
(
z2
-
3.0
*
x2
+
0.25
)
*
y
;
dl4y
=
sr96
*
(
z2
-
x2
+
0.25
)
*
x
;
dl4z
=
2.0
*
sr96
*
z
*
xy
;
MAKE_ELEMENT_7
(
2
,
7
);
//del2f.m[3][3] = 0.30*(1.0 - 4.0*l2 + n2);
dl2
=
0.3
-
1.2
*
x2
+
0.3
*
z2
;
dl2x
=
-
2.4
*
x
;
dl2y
=
0.0
;
dl2z
=
0.6
*
z
;
//del4f.m[3][3] = -0.4*m2 - 2.5*(l2 - 0.6*m2)*n2;
dl4
=
-
0.4
*
y2
-
2.5
*
(
x2
-
0.6
*
y2
)
*
z2
;
dl4x
=
-
5.0
*
x
*
z2
;
dl4y
=
y
*
(
3.0
*
z2
-
0.8
);
dl4z
=
-
5.0
*
(
x2
-
0.6
*
y2
)
*
z
;
MAKE_ELEMENT_7
(
3
,
3
);
//del2f.m[3][4] = -sr09*mn;
dl2
=
-
sr09
*
yz
;
dl2x
=
0.0
;
dl2y
=
-
sr09
*
z
;
dl2z
=
-
sr09
*
y
;
//del4f.m[3][4] = -sr01*(6.0*n2 - 8.0*l2 - 1.0)*mn;
dl4
=
-
sr01
*
(
6.0
*
z2
-
8.0
*
x2
-
1.0
)
*
yz
;
dl4x
=
16.0
*
sr01
*
x
*
yz
;
dl4y
=
-
sr01
*
(
6.0
*
z2
-
8.0
*
x2
-
1.0
)
*
z
;
dl4z
=
-
sr01
*
(
18.0
*
z2
-
8.0
*
x2
-
1.0
)
*
y
;
MAKE_ELEMENT_7
(
3
,
4
);
//del2f.m[3][5] = sr09*ln;
dl2
=
sr09
*
xz
;
dl2x
=
sr09
*
z
;
dl2y
=
0.0
;
dl2z
=
sr09
*
x
;
//del4f.m[3][5] = sr01*(2.0*n2 - 8.0*l2 + 3.0)*ln;
dl4
=
sr01
*
(
2.0
*
z2
-
8.0
*
x2
+
3.0
)
*
xz
;
dl4x
=
sr01
*
(
2.0
*
z2
-
24.0
*
x2
+
3.0
)
*
z
;
dl4y
=
0.0
;
dl4z
=
sr01
*
(
6.0
*
z2
-
8.0
*
x2
+
3.0
)
*
x
;
MAKE_ELEMENT_7
(
3
,
5
);
//del2f.m[3][6] = sr06*lm;
dl2
=
sr06
*
xy
;
dl2x
=
sr06
*
y
;
dl2y
=
sr06
*
x
;
dl2z
=
0.0
;
//del4f.m[3][6] = sr96*(m2 - n2 - 0.25)*lm;
dl4
=
sr96
*
(
y2
-
z2
-
0.25
)
*
xy
;
dl4x
=
sr96
*
(
y2
-
z2
-
0.25
)
*
y
;
dl4y
=
sr96
*
(
3.0
*
y2
-
z2
-
0.25
)
*
x
;
dl4z
=
-
2.0
*
sr96
*
z
*
xy
;
MAKE_ELEMENT_7
(
3
,
6
);
//del2f.m[3][7] = -sr015*(l2 - m2);
dl2
=
-
sr015
*
(
x2
-
y2
);
dl2x
=
-
2.0
*
sr015
*
x
;
dl2y
=
2.0
*
sr015
*
y
;
dl2z
=
0.0
;
//del4f.m[3][7] = sr375*(l2 - m2 + 1.4*m4 - 1.2*l2*m2 - l4);
dl4
=
sr375
*
(
x2
-
y2
+
1.4
*
y4
-
1.2
*
x2
*
y2
-
x4
);
dl4x
=
sr375
*
(
2.0
-
2.4
*
y2
-
4.0
*
x2
)
*
x
;
dl4y
=
sr375
*
(
-
2.0
+
5.6
*
y2
-
2.4
*
x2
)
*
y
;
dl4z
=
0.0
;
MAKE_ELEMENT_7
(
3
,
7
);
//del2f.m[4][4] = 0.0;
dl2
=
0.0
;
dl2x
=
0.0
;
dl2y
=
0.0
;
dl2z
=
0.0
;
//del4f.m[4][4] = (2.0 - 3.0*n2)*n2 - 4.0*l2*m2;
dl4
=
(
2.0
-
3.0
*
z2
)
*
z2
-
4.0
*
x2
*
y2
;
dl4x
=
-
8.0
*
x
*
y2
;
dl4y
=
-
8.0
*
x2
*
y
;
dl4z
=
(
4.0
-
12.0
*
z2
)
*
z
;
MAKE_ELEMENT_7
(
4
,
4
);
//del2f.m[4][5] = 0.0;
dl2
=
0.0
;
dl2x
=
0.0
;
dl2y
=
0.0
;
dl2z
=
0.0
;
//del4f.m[4][5] = 2.0*(l2 - m2)*lm;
dl4
=
2.0
*
(
x2
-
y2
)
*
xy
;
dl4x
=
2.0
*
(
3.0
*
x2
-
y2
)
*
y
;
dl4y
=
2.0
*
(
x2
-
3.0
*
y2
)
*
x
;
dl4z
=
0.0
;
MAKE_ELEMENT_7
(
4
,
5
);
//del2f.m[4][6] = sr15*ln;
dl2
=
sr15
*
xz
;
dl2x
=
sr15
*
z
;
dl2y
=
0.0
;
dl2z
=
sr15
*
x
;
//del4f.m[4][6] = -sr15*(2.0*n2 - 1.0)*ln;
dl4
=
-
sr15
*
(
2.0
*
z2
-
1.0
)
*
xz
;
dl4x
=
-
sr15
*
(
2.0
*
z2
-
1.0
)
*
z
;
dl4y
=
0.0
;
dl4z
=
-
sr15
*
(
6.0
*
z2
-
1.0
)
*
x
;
MAKE_ELEMENT_7
(
4
,
6
);
//del2f.m[4][7] = sr15*mn;
dl2
=
sr15
*
yz
;
dl2x
=
0.0
;
dl2y
=
sr15
*
z
;
dl2z
=
sr15
*
y
;
//del4f.m[4][7] = -sr15*(2.0*n2 - 1.0)*mn;
dl4
=
-
sr15
*
(
2.0
*
z2
-
1.0
)
*
yz
;
dl4x
=
0.0
;
dl4y
=
-
sr15
*
(
2.0
*
z2
-
1.0
)
*
z
;
dl4z
=
-
sr15
*
(
6.0
*
z2
-
1.0
)
*
y
;
MAKE_ELEMENT_7
(
4
,
7
);
//del2f.m[5][5] = 0.0;
dl2
=
0.0
;
dl2x
=
0.0
;
dl2y
=
0.0
;
dl2z
=
0.0
;
//del4f.m[5][5] = -pow((2.0*n2 - 1.0),2) + 4.0*l2*m2;
t1
=
2.0
*
z2
-
1.0
;
dl4
=
-
t1
*
t1
+
4.0
*
x2
*
y2
;
dl4x
=
8.0
*
x
*
y2
;
dl4y
=
8.0
*
x2
*
y
;
dl4z
=
-
8.0
*
t1
*
z
;
MAKE_ELEMENT_7
(
5
,
5
);
//del2f.m[5][6] = -sr15*mn;
dl2
=
-
sr15
*
yz
;
dl2x
=
0.0
;
dl2y
=
-
sr15
*
z
;
dl2z
=
-
sr15
*
y
;
//del4f.m[5][6] = sr15*(2.0*n2 - 1.0)*mn;
dl4
=
sr15
*
(
2.0
*
z2
-
1.0
)
*
yz
;
dl4x
=
0.0
;
dl4y
=
sr15
*
(
2.0
*
z2
-
1.0
)
*
z
;
dl4z
=
sr15
*
(
6.0
*
z2
-
1.0
)
*
y
;
MAKE_ELEMENT_7
(
5
,
6
);
//del2f.m[5][7] = sr15*ln;
dl2
=
sr15
*
xz
;
dl2x
=
sr15
*
z
;
dl2y
=
0.0
;
dl2z
=
sr15
*
x
;
//del4f.m[5][7] = -sr15*(2.0*n2 - 1.0)*ln;
dl4
=
-
sr15
*
(
2.0
*
z2
-
1.0
)
*
xz
;
dl4x
=
-
sr15
*
(
2.0
*
z2
-
1.0
)
*
z
;
dl4y
=
0.0
;
dl4z
=
-
sr15
*
(
6.0
*
z2
-
1.0
)
*
x
;
MAKE_ELEMENT_7
(
5
,
7
);
//del2f.m[6][6] = -(3.0*n2 - 1.0)/2.0;
dl2
=
0.5
-
1.5
*
z2
;
dl2x
=
0.0
;
dl2y
=
0.0
;
dl2z
=
-
3.0
*
z
;
//del4f.m[6][6] = 1.5*(n2 - 1.0)*n2;
dl4
=
(
1.5
*
z2
-
1.5
)
*
z2
;
dl4x
=
0.0
;
dl4y
=
0.0
;
dl4z
=
(
6.0
*
z2
-
3.0
)
*
z
;
MAKE_ELEMENT_7
(
6
,
6
);
//del2f.m[6][7] = 0.0;
dl2
=
0.0
;
dl2x
=
0.0
;
dl2y
=
0.0
;
dl2z
=
0.0
;
//del4f.m[6][7] = 0.0;
dl4
=
0.0
;
dl4x
=
0.0
;
dl4y
=
0.0
;
dl4z
=
0.0
;
MAKE_ELEMENT_7
(
6
,
7
);
//del2f.m[7][7] = -(3.0*n2 - 1.0)/2.0;
dl2
=
0.5
-
1.5
*
z2
;
dl2x
=
0.0
;
dl2y
=
0.0
;
dl2z
=
-
3.0
*
z
;
//del4f.m[7][7] = 1.5*(n2 - 1.0)*n2;
dl4
=
(
1.5
*
z2
-
1.5
)
*
z2
;
dl4x
=
0.0
;
dl4y
=
0.0
;
dl4z
=
(
6.0
*
z2
-
3.0
)
*
z
;
MAKE_ELEMENT_7
(
7
,
7
);
}
/************************************************************************/
/* ----------------------------------------------------------------------
* Fast Model Generalized Pseudopotential Theory (MGPT) interatomic
* potential routine.
*
* Copyright (2015) Lawrence Livermore National Security, LLC.
* Produced at the Lawrence Livermore National Laboratory.
* Written by Tomas Oppelstrup (oppelstrup2@llnl.gov) and John Moriarty
* (moriarty2@llnl.gov)
* LLNL-CODE-674031 All rights reserved.
*
* 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) version 2, dated June 1991.
*
* 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 terms and conditions of the
* GNU General Public License for more details.
*
* LLNL Preamble Notice
* A. This notice is required to be provided under our contract with the
* U.S. Department of Energy (DOE). This work was performed under the auspices
* of the DOE by Lawrence Livermore National Laboratory under Contract No.
* DE-AC52-07NA27344.
*
* B. Neither the United States Government nor Lawrence Livermore National
* Security, LLC nor any of their employees, makes any warranty, express or
* implied, or assumes any liability or responsibility for the accuracy,
* completeness, or usefulness of any information, apparatus, product, or
* process disclosed, or represents that its use would not infringe
* privately-owned rights.
*
* C. Also, reference herein to any specific commercial products, process,
* or services by trade name, trademark, manufacturer or otherwise does not
* necessarily constitute or imply its endorsement, recommendation, or
* favoring by the United States Government or Lawrence Livermore National
* Security, LLC. The views and opinions of authors expressed herein do not
* necessarily state or reflect those of the United States Government or
* Lawrence Livermore National Security, LLC, and shall not be used for
* advertising or product endorsement purposes.
------------------------------------------------------------------------- */
Event Timeline
Log In to Comment