Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F65049917
pair_multi_lucy.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, May 31, 09:27
Size
24 KB
Mime Type
text/x-c
Expires
Sun, Jun 2, 09:27 (2 d)
Engine
blob
Format
Raw Data
Handle
17995342
Attached To
rLAMMPS lammps
pair_multi_lucy.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:
James Larentzos and Joshua Moore (U.S. Army Research Laboratory)
Please cite the related publications:
J.D. Moore, B.C. Barnes, S. Izvekov, M. Lisal, M.S. Sellers, D.E. Taylor & J.K. Brennan
"A coarse-grain force field for RDX: Density dependent and energy conserving"
The Journal of Chemical Physics, 2016, 144, 104501.
------------------------------------------------------------------------------------------- */
#include <mpi.h>
#include <math.h>
#include "math_const.h"
#include <stdlib.h>
#include <string.h>
#include "pair_multi_lucy.h"
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
#include "citeme.h"
using
namespace
LAMMPS_NS
;
enum
{
NONE
,
RLINEAR
,
RSQ
};
#define MAXLINE 1024
static
const
char
cite_pair_multi_lucy
[]
=
"pair_style multi/lucy command:
\n\n
"
"@Article{Moore16,
\n
"
" author = {J.D. Moore, B.C. Barnes, S. Izvekov, M. Lisal, M.S. Sellers, D.E. Taylor and J. K. Brennan},
\n
"
" title = {A coarse-grain force field for RDX: Density dependent and energy conserving},
\n
"
" journal = {J. Chem. Phys.},
\n
"
" year = 2016,
\n
"
" volume = 144
\n
"
" pages = {104501}
\n
"
"}
\n\n
"
;
/* ---------------------------------------------------------------------- */
PairMultiLucy
::
PairMultiLucy
(
LAMMPS
*
lmp
)
:
Pair
(
lmp
),
ntables
(
0
),
tables
(
NULL
),
tabindex
(
NULL
)
{
if
(
lmp
->
citeme
)
lmp
->
citeme
->
add
(
cite_pair_multi_lucy
);
if
(
atom
->
rho_flag
!=
1
)
error
->
all
(
FLERR
,
"Pair multi/lucy command requires atom_style with density (e.g. dpd, meso)"
);
ntables
=
0
;
tables
=
NULL
;
comm_forward
=
1
;
comm_reverse
=
1
;
}
/* ---------------------------------------------------------------------- */
PairMultiLucy
::~
PairMultiLucy
()
{
for
(
int
m
=
0
;
m
<
ntables
;
m
++
)
free_table
(
&
tables
[
m
]);
memory
->
sfree
(
tables
);
if
(
allocated
)
{
memory
->
destroy
(
setflag
);
memory
->
destroy
(
cutsq
);
memory
->
destroy
(
tabindex
);
}
}
/* ---------------------------------------------------------------------- */
void
PairMultiLucy
::
compute
(
int
eflag
,
int
vflag
)
{
int
i
,
j
,
ii
,
jj
,
inum
,
jnum
,
itype
,
jtype
,
itable
;
double
xtmp
,
ytmp
,
ztmp
,
delx
,
dely
,
delz
,
evdwl
,
fpair
,
rsq
;
int
*
ilist
,
*
jlist
,
*
numneigh
,
**
firstneigh
;
Table
*
tb
;
int
tlm1
=
tablength
-
1
;
evdwl
=
0.0
;
if
(
eflag
||
vflag
)
ev_setup
(
eflag
,
vflag
);
else
evflag
=
vflag_fdotr
=
0
;
double
**
x
=
atom
->
x
;
double
**
f
=
atom
->
f
;
int
*
type
=
atom
->
type
;
int
nlocal
=
atom
->
nlocal
;
int
newton_pair
=
force
->
newton_pair
;
double
pi
=
MathConst
::
MY_PI
;
double
A_i
;
double
A_j
;
double
fraction_i
,
fraction_j
;
int
jtable
;
double
*
rho
=
atom
->
rho
;
inum
=
list
->
inum
;
ilist
=
list
->
ilist
;
numneigh
=
list
->
numneigh
;
firstneigh
=
list
->
firstneigh
;
computeLocalDensity
();
// loop over neighbors of my atoms
for
(
ii
=
0
;
ii
<
inum
;
ii
++
)
{
i
=
ilist
[
ii
];
xtmp
=
x
[
i
][
0
];
ytmp
=
x
[
i
][
1
];
ztmp
=
x
[
i
][
2
];
itype
=
type
[
i
];
jlist
=
firstneigh
[
i
];
jnum
=
numneigh
[
i
];
for
(
jj
=
0
;
jj
<
jnum
;
jj
++
)
{
j
=
jlist
[
jj
];
j
&=
NEIGHMASK
;
delx
=
xtmp
-
x
[
j
][
0
];
dely
=
ytmp
-
x
[
j
][
1
];
delz
=
ztmp
-
x
[
j
][
2
];
rsq
=
delx
*
delx
+
dely
*
dely
+
delz
*
delz
;
jtype
=
type
[
j
];
if
(
rsq
<
cutsq
[
itype
][
jtype
])
{
tb
=
&
tables
[
tabindex
[
itype
][
jtype
]];
if
(
rho
[
i
]
*
rho
[
i
]
<
tb
->
innersq
||
rho
[
j
]
*
rho
[
j
]
<
tb
->
innersq
)
error
->
one
(
FLERR
,
"Density < table inner cutoff"
);
if
(
tabstyle
==
LOOKUP
)
{
itable
=
static_cast
<
int
>
(((
rho
[
i
]
*
rho
[
i
])
-
tb
->
innersq
)
*
tb
->
invdelta
);
jtable
=
static_cast
<
int
>
(((
rho
[
j
]
*
rho
[
j
])
-
tb
->
innersq
)
*
tb
->
invdelta
);
if
(
itable
>=
tlm1
||
jtable
>=
tlm1
)
error
->
one
(
FLERR
,
"Density > table outer cutoff"
);
A_i
=
tb
->
f
[
itable
];
A_j
=
tb
->
f
[
jtable
];
fpair
=
0.5
*
(
A_i
+
A_j
)
*
(
1.0
+
3.0
*
sqrt
(
rsq
)
/
sqrt
(
cutsq
[
itype
][
jtype
]))
*
(
1.0
-
sqrt
(
rsq
)
/
sqrt
(
cutsq
[
itype
][
jtype
]))
*
(
1.0
-
sqrt
(
rsq
)
/
sqrt
(
cutsq
[
itype
][
jtype
]))
*
(
1.0
-
sqrt
(
rsq
)
/
sqrt
(
cutsq
[
itype
][
jtype
]));
fpair
=
fpair
/
sqrt
(
rsq
);
}
else
if
(
tabstyle
==
LINEAR
)
{
itable
=
static_cast
<
int
>
((
rho
[
i
]
*
rho
[
i
]
-
tb
->
innersq
)
*
tb
->
invdelta
);
jtable
=
static_cast
<
int
>
(((
rho
[
j
]
*
rho
[
j
])
-
tb
->
innersq
)
*
tb
->
invdelta
);
if
(
itable
>=
tlm1
||
jtable
>=
tlm1
)
error
->
one
(
FLERR
,
"Density > table outer cutoff"
);
fraction_i
=
(((
rho
[
i
]
*
rho
[
i
])
-
tb
->
rsq
[
itable
])
*
tb
->
invdelta
);
fraction_j
=
(((
rho
[
j
]
*
rho
[
j
])
-
tb
->
rsq
[
jtable
])
*
tb
->
invdelta
);
A_i
=
tb
->
f
[
itable
]
+
fraction_i
*
tb
->
df
[
itable
];
A_j
=
tb
->
f
[
jtable
]
+
fraction_j
*
tb
->
df
[
jtable
];
fpair
=
0.5
*
(
A_i
+
A_j
)
*
(
1.0
+
3.0
*
sqrt
(
rsq
)
/
sqrt
(
cutsq
[
itype
][
jtype
]))
*
(
1.0
-
sqrt
(
rsq
)
/
sqrt
(
cutsq
[
itype
][
jtype
]))
*
(
1.0
-
sqrt
(
rsq
)
/
sqrt
(
cutsq
[
itype
][
jtype
]))
*
(
1.0
-
sqrt
(
rsq
)
/
sqrt
(
cutsq
[
itype
][
jtype
]));
fpair
=
fpair
/
sqrt
(
rsq
);
}
else
error
->
one
(
FLERR
,
"Only LOOKUP and LINEAR table styles have been implemented for pair multi/lucy"
);
f
[
i
][
0
]
+=
delx
*
fpair
;
f
[
i
][
1
]
+=
dely
*
fpair
;
f
[
i
][
2
]
+=
delz
*
fpair
;
if
(
newton_pair
||
j
<
nlocal
)
{
f
[
j
][
0
]
-=
delx
*
fpair
;
f
[
j
][
1
]
-=
dely
*
fpair
;
f
[
j
][
2
]
-=
delz
*
fpair
;
}
if
(
evflag
)
ev_tally
(
i
,
j
,
nlocal
,
newton_pair
,
0.0
,
0.0
,
fpair
,
delx
,
dely
,
delz
);
}
}
tb
=
&
tables
[
tabindex
[
itype
][
itype
]];
if
(
rho
[
i
]
*
rho
[
i
]
<
tb
->
innersq
||
rho
[
j
]
*
rho
[
j
]
<
tb
->
innersq
)
error
->
one
(
FLERR
,
"Density < table inner cutoff"
);
itable
=
static_cast
<
int
>
(((
rho
[
i
]
*
rho
[
i
])
-
tb
->
innersq
)
*
tb
->
invdelta
);
if
(
tabstyle
==
LOOKUP
)
evdwl
=
tb
->
e
[
itable
];
else
if
(
tabstyle
==
LINEAR
){
if
(
itable
>=
tlm1
)
error
->
one
(
FLERR
,
"Density > table outer cutoff"
);
if
(
itable
==
0
)
fraction_i
=
0.0
;
else
fraction_i
=
(((
rho
[
i
]
*
rho
[
i
])
-
tb
->
rsq
[
itable
])
*
tb
->
invdelta
);
evdwl
=
tb
->
e
[
itable
]
+
fraction_i
*
tb
->
de
[
itable
];
}
else
error
->
one
(
FLERR
,
"Only LOOKUP and LINEAR table styles have been implemented for pair multi/lucy"
);
evdwl
*=
(
pi
*
cutsq
[
itype
][
itype
]
*
cutsq
[
itype
][
itype
])
/
84.0
;
if
(
evflag
)
ev_tally
(
0
,
0
,
nlocal
,
newton_pair
,
evdwl
,
0.0
,
0.0
,
0.0
,
0.0
,
0.0
);
}
if
(
vflag_fdotr
)
virial_fdotr_compute
();
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void
PairMultiLucy
::
allocate
()
{
allocated
=
1
;
const
int
nt
=
atom
->
ntypes
+
1
;
memory
->
create
(
setflag
,
nt
,
nt
,
"pair:setflag"
);
memory
->
create
(
cutsq
,
nt
,
nt
,
"pair:cutsq"
);
memory
->
create
(
tabindex
,
nt
,
nt
,
"pair:tabindex"
);
memset
(
&
setflag
[
0
][
0
],
0
,
nt
*
nt
*
sizeof
(
int
));
memset
(
&
cutsq
[
0
][
0
],
0
,
nt
*
nt
*
sizeof
(
double
));
memset
(
&
tabindex
[
0
][
0
],
0
,
nt
*
nt
*
sizeof
(
int
));
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void
PairMultiLucy
::
settings
(
int
narg
,
char
**
arg
)
{
if
(
narg
!=
2
)
error
->
all
(
FLERR
,
"Illegal pair_style command"
);
// new settings
if
(
strcmp
(
arg
[
0
],
"lookup"
)
==
0
)
tabstyle
=
LOOKUP
;
else
if
(
strcmp
(
arg
[
0
],
"linear"
)
==
0
)
tabstyle
=
LINEAR
;
else
error
->
all
(
FLERR
,
"Unknown table style in pair_style command"
);
tablength
=
force
->
inumeric
(
FLERR
,
arg
[
1
]);
if
(
tablength
<
2
)
error
->
all
(
FLERR
,
"Illegal number of pair table entries"
);
// delete old tables, since cannot just change settings
for
(
int
m
=
0
;
m
<
ntables
;
m
++
)
free_table
(
&
tables
[
m
]);
memory
->
sfree
(
tables
);
if
(
allocated
)
{
memory
->
destroy
(
setflag
);
memory
->
destroy
(
cutsq
);
memory
->
destroy
(
tabindex
);
}
allocated
=
0
;
ntables
=
0
;
tables
=
NULL
;
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void
PairMultiLucy
::
coeff
(
int
narg
,
char
**
arg
)
{
if
(
narg
!=
4
&&
narg
!=
5
)
error
->
all
(
FLERR
,
"Illegal pair_coeff command"
);
if
(
!
allocated
)
allocate
();
int
ilo
,
ihi
,
jlo
,
jhi
;
force
->
bounds
(
FLERR
,
arg
[
0
],
atom
->
ntypes
,
ilo
,
ihi
);
force
->
bounds
(
FLERR
,
arg
[
1
],
atom
->
ntypes
,
jlo
,
jhi
);
int
me
;
MPI_Comm_rank
(
world
,
&
me
);
tables
=
(
Table
*
)
memory
->
srealloc
(
tables
,(
ntables
+
1
)
*
sizeof
(
Table
),
"pair:tables"
);
Table
*
tb
=
&
tables
[
ntables
];
null_table
(
tb
);
if
(
me
==
0
)
read_table
(
tb
,
arg
[
2
],
arg
[
3
]);
bcast_table
(
tb
);
// set table cutoff
if
(
narg
==
5
)
tb
->
cut
=
force
->
numeric
(
FLERR
,
arg
[
4
]);
else
if
(
tb
->
rflag
)
tb
->
cut
=
tb
->
rhi
;
else
tb
->
cut
=
tb
->
rfile
[
tb
->
ninput
-
1
];
// error check on table parameters
// insure cutoff is within table
if
(
tb
->
ninput
<=
1
)
error
->
one
(
FLERR
,
"Invalid pair table length"
);
double
rlo
;
if
(
tb
->
rflag
==
0
)
{
rlo
=
tb
->
rfile
[
0
];
}
else
{
rlo
=
tb
->
rlo
;
}
rho_0
=
rlo
;
tb
->
match
=
0
;
if
(
tabstyle
==
LINEAR
&&
tb
->
ninput
==
tablength
&&
tb
->
rflag
==
RSQ
)
tb
->
match
=
1
;
// spline read-in values and compute r,e,f vectors within table
if
(
tb
->
match
==
0
)
spline_table
(
tb
);
compute_table
(
tb
);
// store ptr to table in tabindex
int
count
=
0
;
for
(
int
i
=
ilo
;
i
<=
ihi
;
i
++
)
{
for
(
int
j
=
MAX
(
jlo
,
i
);
j
<=
jhi
;
j
++
)
{
tabindex
[
i
][
j
]
=
ntables
;
setflag
[
i
][
j
]
=
1
;
count
++
;
}
}
if
(
count
==
0
)
error
->
all
(
FLERR
,
"Illegal pair_coeff command"
);
ntables
++
;
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double
PairMultiLucy
::
init_one
(
int
i
,
int
j
)
{
if
(
setflag
[
i
][
j
]
==
0
)
error
->
all
(
FLERR
,
"All pair coeffs are not set"
);
tabindex
[
j
][
i
]
=
tabindex
[
i
][
j
];
return
tables
[
tabindex
[
i
][
j
]].
cut
;
}
/* ----------------------------------------------------------------------
read a table section from a tabulated potential file
only called by proc 0
this function sets these values in Table:
ninput,rfile,efile,ffile,rflag,rlo,rhi,fpflag,fplo,fphi
------------------------------------------------------------------------- */
void
PairMultiLucy
::
read_table
(
Table
*
tb
,
char
*
file
,
char
*
keyword
)
{
char
line
[
MAXLINE
];
// open file
FILE
*
fp
=
fopen
(
file
,
"r"
);
if
(
fp
==
NULL
)
{
char
str
[
128
];
sprintf
(
str
,
"Cannot open file %s"
,
file
);
error
->
one
(
FLERR
,
str
);
}
// loop until section found with matching keyword
while
(
1
)
{
if
(
fgets
(
line
,
MAXLINE
,
fp
)
==
NULL
)
error
->
one
(
FLERR
,
"Did not find keyword in table file"
);
if
(
strspn
(
line
,
"
\t\n\r
"
)
==
strlen
(
line
))
continue
;
// blank line
if
(
line
[
0
]
==
'#'
)
continue
;
// comment
char
*
word
=
strtok
(
line
,
"
\t\n\r
"
);
if
(
strcmp
(
word
,
keyword
)
==
0
)
break
;
// matching keyword
fgets
(
line
,
MAXLINE
,
fp
);
// no match, skip section
param_extract
(
tb
,
line
);
fgets
(
line
,
MAXLINE
,
fp
);
for
(
int
i
=
0
;
i
<
tb
->
ninput
;
i
++
)
fgets
(
line
,
MAXLINE
,
fp
);
}
// read args on 2nd line of section
// allocate table arrays for file values
fgets
(
line
,
MAXLINE
,
fp
);
param_extract
(
tb
,
line
);
memory
->
create
(
tb
->
rfile
,
tb
->
ninput
,
"pair:rfile"
);
memory
->
create
(
tb
->
efile
,
tb
->
ninput
,
"pair:efile"
);
memory
->
create
(
tb
->
ffile
,
tb
->
ninput
,
"pair:ffile"
);
// read r,e,f table values from file
// if rflag set, compute r
// if rflag not set, use r from file
int
itmp
;
double
rtmp
;
fgets
(
line
,
MAXLINE
,
fp
);
for
(
int
i
=
0
;
i
<
tb
->
ninput
;
i
++
)
{
fgets
(
line
,
MAXLINE
,
fp
);
sscanf
(
line
,
"%d %lg %lg %lg"
,
&
itmp
,
&
rtmp
,
&
tb
->
efile
[
i
],
&
tb
->
ffile
[
i
]);
if
(
tb
->
rflag
==
RLINEAR
)
rtmp
=
tb
->
rlo
+
(
tb
->
rhi
-
tb
->
rlo
)
*
i
/
(
tb
->
ninput
-
1
);
else
if
(
tb
->
rflag
==
RSQ
)
{
rtmp
=
tb
->
rlo
*
tb
->
rlo
+
(
tb
->
rhi
*
tb
->
rhi
-
tb
->
rlo
*
tb
->
rlo
)
*
i
/
(
tb
->
ninput
-
1
);
rtmp
=
sqrt
(
rtmp
);
}
tb
->
rfile
[
i
]
=
rtmp
;
}
// close file
fclose
(
fp
);
}
/* ----------------------------------------------------------------------
broadcast read-in table info from proc 0 to other procs
this function communicates these values in Table:
ninput,rfile,efile,ffile,rflag,rlo,rhi,fpflag,fplo,fphi
------------------------------------------------------------------------- */
void
PairMultiLucy
::
bcast_table
(
Table
*
tb
)
{
MPI_Bcast
(
&
tb
->
ninput
,
1
,
MPI_INT
,
0
,
world
);
int
me
;
MPI_Comm_rank
(
world
,
&
me
);
if
(
me
>
0
)
{
memory
->
create
(
tb
->
rfile
,
tb
->
ninput
,
"pair:rfile"
);
memory
->
create
(
tb
->
efile
,
tb
->
ninput
,
"pair:efile"
);
memory
->
create
(
tb
->
ffile
,
tb
->
ninput
,
"pair:ffile"
);
}
MPI_Bcast
(
tb
->
rfile
,
tb
->
ninput
,
MPI_DOUBLE
,
0
,
world
);
MPI_Bcast
(
tb
->
efile
,
tb
->
ninput
,
MPI_DOUBLE
,
0
,
world
);
MPI_Bcast
(
tb
->
ffile
,
tb
->
ninput
,
MPI_DOUBLE
,
0
,
world
);
MPI_Bcast
(
&
tb
->
rflag
,
1
,
MPI_INT
,
0
,
world
);
if
(
tb
->
rflag
)
{
MPI_Bcast
(
&
tb
->
rlo
,
1
,
MPI_DOUBLE
,
0
,
world
);
MPI_Bcast
(
&
tb
->
rhi
,
1
,
MPI_DOUBLE
,
0
,
world
);
}
MPI_Bcast
(
&
tb
->
fpflag
,
1
,
MPI_INT
,
0
,
world
);
if
(
tb
->
fpflag
)
{
MPI_Bcast
(
&
tb
->
fplo
,
1
,
MPI_DOUBLE
,
0
,
world
);
MPI_Bcast
(
&
tb
->
fphi
,
1
,
MPI_DOUBLE
,
0
,
world
);
}
}
/* ----------------------------------------------------------------------
build spline representation of e,f over entire range of read-in table
this function sets these values in Table: e2file,f2file
------------------------------------------------------------------------- */
void
PairMultiLucy
::
spline_table
(
Table
*
tb
)
{
memory
->
create
(
tb
->
e2file
,
tb
->
ninput
,
"pair:e2file"
);
memory
->
create
(
tb
->
f2file
,
tb
->
ninput
,
"pair:f2file"
);
double
ep0
=
-
tb
->
ffile
[
0
];
double
epn
=
-
tb
->
ffile
[
tb
->
ninput
-
1
];
spline
(
tb
->
rfile
,
tb
->
efile
,
tb
->
ninput
,
ep0
,
epn
,
tb
->
e2file
);
if
(
tb
->
fpflag
==
0
)
{
tb
->
fplo
=
(
tb
->
ffile
[
1
]
-
tb
->
ffile
[
0
])
/
(
tb
->
rfile
[
1
]
-
tb
->
rfile
[
0
]);
tb
->
fphi
=
(
tb
->
ffile
[
tb
->
ninput
-
1
]
-
tb
->
ffile
[
tb
->
ninput
-
2
])
/
(
tb
->
rfile
[
tb
->
ninput
-
1
]
-
tb
->
rfile
[
tb
->
ninput
-
2
]);
}
double
fp0
=
tb
->
fplo
;
double
fpn
=
tb
->
fphi
;
spline
(
tb
->
rfile
,
tb
->
ffile
,
tb
->
ninput
,
fp0
,
fpn
,
tb
->
f2file
);
}
/* ----------------------------------------------------------------------
extract attributes from parameter line in table section
format of line: N value R/RSQ lo hi FP fplo fphi
N is required, other params are optional
------------------------------------------------------------------------- */
void
PairMultiLucy
::
param_extract
(
Table
*
tb
,
char
*
line
)
{
tb
->
ninput
=
0
;
tb
->
rflag
=
NONE
;
tb
->
fpflag
=
0
;
char
*
word
=
strtok
(
line
,
"
\t\n\r\f
"
);
while
(
word
)
{
if
(
strcmp
(
word
,
"N"
)
==
0
)
{
word
=
strtok
(
NULL
,
"
\t\n\r\f
"
);
tb
->
ninput
=
atoi
(
word
);
}
else
if
(
strcmp
(
word
,
"R"
)
==
0
||
strcmp
(
word
,
"RSQ"
)
==
0
)
{
if
(
strcmp
(
word
,
"R"
)
==
0
)
tb
->
rflag
=
RLINEAR
;
else
if
(
strcmp
(
word
,
"RSQ"
)
==
0
)
tb
->
rflag
=
RSQ
;
word
=
strtok
(
NULL
,
"
\t\n\r\f
"
);
tb
->
rlo
=
atof
(
word
);
word
=
strtok
(
NULL
,
"
\t\n\r\f
"
);
tb
->
rhi
=
atof
(
word
);
}
else
if
(
strcmp
(
word
,
"FP"
)
==
0
)
{
tb
->
fpflag
=
1
;
word
=
strtok
(
NULL
,
"
\t\n\r\f
"
);
tb
->
fplo
=
atof
(
word
);
word
=
strtok
(
NULL
,
"
\t\n\r\f
"
);
tb
->
fphi
=
atof
(
word
);
}
else
{
printf
(
"WORD: %s
\n
"
,
word
);
error
->
one
(
FLERR
,
"Invalid keyword in pair table parameters"
);
}
word
=
strtok
(
NULL
,
"
\t\n\r\f
"
);
}
if
(
tb
->
ninput
==
0
)
error
->
one
(
FLERR
,
"Pair table parameters did not set N"
);
}
/* ----------------------------------------------------------------------
compute r,e,f vectors from splined values
------------------------------------------------------------------------- */
void
PairMultiLucy
::
compute_table
(
Table
*
tb
)
{
int
tlm1
=
tablength
-
1
;
// inner = inner table bound
// cut = outer table bound
// delta = table spacing in rsq for N-1 bins
double
inner
;
if
(
tb
->
rflag
)
inner
=
tb
->
rlo
;
else
inner
=
tb
->
rfile
[
0
];
tb
->
innersq
=
inner
*
inner
;
tb
->
delta
=
(
tb
->
rhi
*
tb
->
rhi
-
tb
->
innersq
)
/
tlm1
;
tb
->
invdelta
=
1.0
/
tb
->
delta
;
// direct lookup tables
// N-1 evenly spaced bins in rsq from inner to cut
// e,f = value at midpt of bin
// e,f are N-1 in length since store 1 value at bin midpt
// f is converted to f/r when stored in f[i]
// e,f are never a match to read-in values, always computed via spline interp
if
(
tabstyle
==
LOOKUP
)
{
memory
->
create
(
tb
->
e
,
tlm1
,
"pair:e"
);
memory
->
create
(
tb
->
f
,
tlm1
,
"pair:f"
);
double
r
,
rsq
;
for
(
int
i
=
0
;
i
<
tlm1
;
i
++
)
{
rsq
=
tb
->
innersq
+
(
i
+
0.5
)
*
tb
->
delta
;
r
=
sqrt
(
rsq
);
tb
->
e
[
i
]
=
splint
(
tb
->
rfile
,
tb
->
efile
,
tb
->
e2file
,
tb
->
ninput
,
r
);
tb
->
f
[
i
]
=
splint
(
tb
->
rfile
,
tb
->
ffile
,
tb
->
f2file
,
tb
->
ninput
,
r
);
}
}
// linear tables
// N-1 evenly spaced bins in rsq from inner to cut
// rsq,e,f = value at lower edge of bin
// de,df values = delta from lower edge to upper edge of bin
// rsq,e,f are N in length so de,df arrays can compute difference
// f is converted to f/r when stored in f[i]
// e,f can match read-in values, else compute via spline interp
if
(
tabstyle
==
LINEAR
)
{
memory
->
create
(
tb
->
rsq
,
tablength
,
"pair:rsq"
);
memory
->
create
(
tb
->
e
,
tablength
,
"pair:e"
);
memory
->
create
(
tb
->
f
,
tablength
,
"pair:f"
);
memory
->
create
(
tb
->
de
,
tlm1
,
"pair:de"
);
memory
->
create
(
tb
->
df
,
tlm1
,
"pair:df"
);
double
r
,
rsq
;
for
(
int
i
=
0
;
i
<
tablength
;
i
++
)
{
rsq
=
tb
->
innersq
+
i
*
tb
->
delta
;
r
=
sqrt
(
rsq
);
tb
->
rsq
[
i
]
=
rsq
;
if
(
tb
->
match
)
{
tb
->
e
[
i
]
=
tb
->
efile
[
i
];
tb
->
f
[
i
]
=
tb
->
ffile
[
i
];
}
else
{
tb
->
e
[
i
]
=
splint
(
tb
->
rfile
,
tb
->
efile
,
tb
->
e2file
,
tb
->
ninput
,
r
);
tb
->
f
[
i
]
=
splint
(
tb
->
rfile
,
tb
->
ffile
,
tb
->
f2file
,
tb
->
ninput
,
r
);
}
}
for
(
int
i
=
0
;
i
<
tlm1
;
i
++
)
{
tb
->
de
[
i
]
=
tb
->
e
[
i
+
1
]
-
tb
->
e
[
i
];
tb
->
df
[
i
]
=
tb
->
f
[
i
+
1
]
-
tb
->
f
[
i
];
}
}
}
/* ----------------------------------------------------------------------
set all ptrs in a table to NULL, so can be freed safely
------------------------------------------------------------------------- */
void
PairMultiLucy
::
null_table
(
Table
*
tb
)
{
tb
->
rfile
=
tb
->
efile
=
tb
->
ffile
=
NULL
;
tb
->
e2file
=
tb
->
f2file
=
NULL
;
tb
->
rsq
=
tb
->
drsq
=
tb
->
e
=
tb
->
de
=
NULL
;
tb
->
f
=
tb
->
df
=
tb
->
e2
=
tb
->
f2
=
NULL
;
}
/* ----------------------------------------------------------------------
free all arrays in a table
------------------------------------------------------------------------- */
void
PairMultiLucy
::
free_table
(
Table
*
tb
)
{
memory
->
destroy
(
tb
->
rfile
);
memory
->
destroy
(
tb
->
efile
);
memory
->
destroy
(
tb
->
ffile
);
memory
->
destroy
(
tb
->
e2file
);
memory
->
destroy
(
tb
->
f2file
);
memory
->
destroy
(
tb
->
rsq
);
memory
->
destroy
(
tb
->
drsq
);
memory
->
destroy
(
tb
->
e
);
memory
->
destroy
(
tb
->
de
);
memory
->
destroy
(
tb
->
f
);
memory
->
destroy
(
tb
->
df
);
memory
->
destroy
(
tb
->
e2
);
memory
->
destroy
(
tb
->
f2
);
}
/* ----------------------------------------------------------------------
spline and splint routines modified from Numerical Recipes
------------------------------------------------------------------------- */
void
PairMultiLucy
::
spline
(
double
*
x
,
double
*
y
,
int
n
,
double
yp1
,
double
ypn
,
double
*
y2
)
{
int
i
,
k
;
double
p
,
qn
,
sig
,
un
;
double
*
u
=
new
double
[
n
];
if
(
yp1
>
0.99e30
)
y2
[
0
]
=
u
[
0
]
=
0.0
;
else
{
y2
[
0
]
=
-
0.5
;
u
[
0
]
=
(
3.0
/
(
x
[
1
]
-
x
[
0
]))
*
((
y
[
1
]
-
y
[
0
])
/
(
x
[
1
]
-
x
[
0
])
-
yp1
);
}
for
(
i
=
1
;
i
<
n
-
1
;
i
++
)
{
sig
=
(
x
[
i
]
-
x
[
i
-
1
])
/
(
x
[
i
+
1
]
-
x
[
i
-
1
]);
p
=
sig
*
y2
[
i
-
1
]
+
2.0
;
y2
[
i
]
=
(
sig
-
1.0
)
/
p
;
u
[
i
]
=
(
y
[
i
+
1
]
-
y
[
i
])
/
(
x
[
i
+
1
]
-
x
[
i
])
-
(
y
[
i
]
-
y
[
i
-
1
])
/
(
x
[
i
]
-
x
[
i
-
1
]);
u
[
i
]
=
(
6.0
*
u
[
i
]
/
(
x
[
i
+
1
]
-
x
[
i
-
1
])
-
sig
*
u
[
i
-
1
])
/
p
;
}
if
(
ypn
>
0.99e30
)
qn
=
un
=
0.0
;
else
{
qn
=
0.5
;
un
=
(
3.0
/
(
x
[
n
-
1
]
-
x
[
n
-
2
]))
*
(
ypn
-
(
y
[
n
-
1
]
-
y
[
n
-
2
])
/
(
x
[
n
-
1
]
-
x
[
n
-
2
]));
}
y2
[
n
-
1
]
=
(
un
-
qn
*
u
[
n
-
2
])
/
(
qn
*
y2
[
n
-
2
]
+
1.0
);
for
(
k
=
n
-
2
;
k
>=
0
;
k
--
)
y2
[
k
]
=
y2
[
k
]
*
y2
[
k
+
1
]
+
u
[
k
];
delete
[]
u
;
}
/* ---------------------------------------------------------------------- */
double
PairMultiLucy
::
splint
(
double
*
xa
,
double
*
ya
,
double
*
y2a
,
int
n
,
double
x
)
{
int
klo
,
khi
,
k
;
double
h
,
b
,
a
,
y
;
klo
=
0
;
khi
=
n
-
1
;
while
(
khi
-
klo
>
1
)
{
k
=
(
khi
+
klo
)
>>
1
;
if
(
xa
[
k
]
>
x
)
khi
=
k
;
else
klo
=
k
;
}
h
=
xa
[
khi
]
-
xa
[
klo
];
a
=
(
xa
[
khi
]
-
x
)
/
h
;
b
=
(
x
-
xa
[
klo
])
/
h
;
y
=
a
*
ya
[
klo
]
+
b
*
ya
[
khi
]
+
((
a
*
a
*
a
-
a
)
*
y2a
[
klo
]
+
(
b
*
b
*
b
-
b
)
*
y2a
[
khi
])
*
(
h
*
h
)
/
6.0
;
return
y
;
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void
PairMultiLucy
::
write_restart
(
FILE
*
fp
)
{
write_restart_settings
(
fp
);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void
PairMultiLucy
::
read_restart
(
FILE
*
fp
)
{
read_restart_settings
(
fp
);
allocate
();
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void
PairMultiLucy
::
write_restart_settings
(
FILE
*
fp
)
{
fwrite
(
&
tabstyle
,
sizeof
(
int
),
1
,
fp
);
fwrite
(
&
tablength
,
sizeof
(
int
),
1
,
fp
);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void
PairMultiLucy
::
read_restart_settings
(
FILE
*
fp
)
{
if
(
comm
->
me
==
0
)
{
fread
(
&
tabstyle
,
sizeof
(
int
),
1
,
fp
);
fread
(
&
tablength
,
sizeof
(
int
),
1
,
fp
);
}
MPI_Bcast
(
&
tabstyle
,
1
,
MPI_INT
,
0
,
world
);
MPI_Bcast
(
&
tablength
,
1
,
MPI_INT
,
0
,
world
);
}
/* ---------------------------------------------------------------------- */
void
PairMultiLucy
::
computeLocalDensity
()
{
int
i
,
j
,
m
,
ii
,
jj
,
inum
,
jnum
,
itype
,
jtype
;
double
xtmp
,
ytmp
,
ztmp
,
delx
,
dely
,
delz
;
double
rsq
;
int
*
ilist
,
*
jlist
,
*
numneigh
,
**
firstneigh
;
double
**
x
=
atom
->
x
;
int
*
type
=
atom
->
type
;
int
nlocal
=
atom
->
nlocal
;
int
newton_pair
=
force
->
newton_pair
;
double
factor
;
inum
=
list
->
inum
;
ilist
=
list
->
ilist
;
numneigh
=
list
->
numneigh
;
firstneigh
=
list
->
firstneigh
;
double
pi
=
MathConst
::
MY_PI
;
double
*
rho
=
atom
->
rho
;
// zero out density
if
(
newton_pair
)
{
m
=
nlocal
+
atom
->
nghost
;
for
(
i
=
0
;
i
<
m
;
i
++
)
rho
[
i
]
=
0.0
;
}
else
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
rho
[
i
]
=
0.0
;
// rho = density at each atom
// loop over neighbors of my atoms
for
(
ii
=
0
;
ii
<
inum
;
ii
++
)
{
i
=
ilist
[
ii
];
xtmp
=
x
[
i
][
0
];
ytmp
=
x
[
i
][
1
];
ztmp
=
x
[
i
][
2
];
itype
=
type
[
i
];
jlist
=
firstneigh
[
i
];
jnum
=
numneigh
[
i
];
for
(
jj
=
0
;
jj
<
jnum
;
jj
++
)
{
j
=
jlist
[
jj
];
j
&=
NEIGHMASK
;
jtype
=
type
[
j
];
delx
=
xtmp
-
x
[
j
][
0
];
dely
=
ytmp
-
x
[
j
][
1
];
delz
=
ztmp
-
x
[
j
][
2
];
rsq
=
delx
*
delx
+
dely
*
dely
+
delz
*
delz
;
if
(
rsq
<
cutsq
[
itype
][
jtype
])
{
double
rcut
=
sqrt
(
cutsq
[
itype
][
jtype
]);
factor
=
(
84.0
/
(
5.0
*
pi
*
rcut
*
rcut
*
rcut
))
*
(
1.0
+
3.0
*
sqrt
(
rsq
)
/
(
2.0
*
rcut
))
*
(
1.0
-
sqrt
(
rsq
)
/
rcut
)
*
(
1.0
-
sqrt
(
rsq
)
/
rcut
)
*
(
1.0
-
sqrt
(
rsq
)
/
rcut
)
*
(
1.0
-
sqrt
(
rsq
)
/
rcut
);
rho
[
i
]
+=
factor
;
if
(
newton_pair
||
j
<
nlocal
)
{
rho
[
j
]
+=
factor
;
}
}
}
}
if
(
newton_pair
)
comm
->
reverse_comm_pair
(
this
);
comm
->
forward_comm_pair
(
this
);
}
/* ---------------------------------------------------------------------- */
int
PairMultiLucy
::
pack_forward_comm
(
int
n
,
int
*
list
,
double
*
buf
,
int
pbc_flag
,
int
*
pbc
)
{
int
i
,
j
,
m
;
double
*
rho
=
atom
->
rho
;
m
=
0
;
for
(
i
=
0
;
i
<
n
;
i
++
)
{
j
=
list
[
i
];
buf
[
m
++
]
=
rho
[
j
];
}
return
m
;
}
/* ---------------------------------------------------------------------- */
void
PairMultiLucy
::
unpack_forward_comm
(
int
n
,
int
first
,
double
*
buf
)
{
int
i
,
m
,
last
;
double
*
rho
=
atom
->
rho
;
m
=
0
;
last
=
first
+
n
;
for
(
i
=
first
;
i
<
last
;
i
++
)
rho
[
i
]
=
buf
[
m
++
];
}
/* ---------------------------------------------------------------------- */
int
PairMultiLucy
::
pack_reverse_comm
(
int
n
,
int
first
,
double
*
buf
)
{
int
i
,
m
,
last
;
double
*
rho
=
atom
->
rho
;
m
=
0
;
last
=
first
+
n
;
for
(
i
=
first
;
i
<
last
;
i
++
)
buf
[
m
++
]
=
rho
[
i
];
return
m
;
}
/* ---------------------------------------------------------------------- */
void
PairMultiLucy
::
unpack_reverse_comm
(
int
n
,
int
*
list
,
double
*
buf
)
{
int
i
,
j
,
m
;
double
*
rho
=
atom
->
rho
;
m
=
0
;
for
(
i
=
0
;
i
<
n
;
i
++
)
{
j
=
list
[
i
];
rho
[
j
]
+=
buf
[
m
++
];
}
}
Event Timeline
Log In to Comment