Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F65868491
pair_tdpd.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
Thu, Jun 6, 18:28
Size
14 KB
Mime Type
text/x-c
Expires
Sat, Jun 8, 18:28 (2 d)
Engine
blob
Format
Raw Data
Handle
18139098
Attached To
rLAMMPS lammps
pair_tdpd.cpp
View Options
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Zhen Li (Brown University)
Email: zhen_li@brown.edu
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "pair_tdpd.h"
#include "atom.h"
#include "atom_vec.h"
#include "comm.h"
#include "update.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "random_mars.h"
#include "citeme.h"
#include "memory.h"
#include "error.h"
#include <time.h>
using
namespace
LAMMPS_NS
;
#define MIN(A,B) ((A) < (B) ? (A) : (B))
#define MAX(A,B) ((A) > (B) ? (A) : (B))
#define EPSILON 1.0e-10
static
const
char
cite_pair_tdpd
[]
=
"pair tdpd command:
\n\n
"
"@Article{ZLi2015_JCP,
\n
"
" author = {Li, Z. and Yazdani, A. and Tartakovsky, A. and Karniadakis, G.E.},
\n
"
" title = {Transport dissipative particle dynamics model for mesoscopic advection-diffusion-reaction problems},
\n
"
" journal = {The Journal of Chemical Physics},
\n
"
" year = {2015},
\n
"
" volume = {143},
\n
"
" pages = {014101}
\n
"
"}
\n\n
"
;
/* ---------------------------------------------------------------------- */
PairTDPD
::
PairTDPD
(
LAMMPS
*
lmp
)
:
Pair
(
lmp
)
{
if
(
lmp
->
citeme
)
lmp
->
citeme
->
add
(
cite_pair_tdpd
);
cc_species
=
atom
->
cc_species
;
writedata
=
1
;
random
=
NULL
;
}
/* ---------------------------------------------------------------------- */
PairTDPD
::~
PairTDPD
()
{
if
(
allocated
)
{
memory
->
destroy
(
setflag
);
memory
->
destroy
(
cutsq
);
memory
->
destroy
(
cut
);
memory
->
destroy
(
cutcc
);
memory
->
destroy
(
a0
);
memory
->
destroy
(
gamma
);
memory
->
destroy
(
sigma
);
memory
->
destroy
(
power
);
memory
->
destroy
(
kappa
);
memory
->
destroy
(
epsilon
);
memory
->
destroy
(
powercc
);
}
if
(
random
)
delete
random
;
}
/* ---------------------------------------------------------------------- */
void
PairTDPD
::
compute
(
int
eflag
,
int
vflag
)
{
double
evdwl
=
0.0
;
if
(
eflag
||
vflag
)
ev_setup
(
eflag
,
vflag
);
else
evflag
=
vflag_fdotr
=
0
;
double
**
x
=
atom
->
x
;
double
**
v
=
atom
->
v
;
double
**
f
=
atom
->
f
;
double
**
cc
=
atom
->
cc
;
double
**
cc_flux
=
atom
->
cc_flux
;
int
*
type
=
atom
->
type
;
int
nlocal
=
atom
->
nlocal
;
double
*
special_lj
=
force
->
special_lj
;
int
newton_pair
=
force
->
newton_pair
;
double
dtinvsqrt
=
1.0
/
sqrt
(
update
->
dt
);
int
inum
=
list
->
inum
;
int
*
ilist
=
list
->
ilist
;
int
*
numneigh
=
list
->
numneigh
;
int
**
firstneigh
=
list
->
firstneigh
;
// loop over neighbors of my atoms
for
(
int
ii
=
0
;
ii
<
inum
;
ii
++
)
{
int
i
=
ilist
[
ii
];
double
xtmp
=
x
[
i
][
0
];
double
ytmp
=
x
[
i
][
1
];
double
ztmp
=
x
[
i
][
2
];
double
vxtmp
=
v
[
i
][
0
];
double
vytmp
=
v
[
i
][
1
];
double
vztmp
=
v
[
i
][
2
];
int
itype
=
type
[
i
];
int
*
jlist
=
firstneigh
[
i
];
int
jnum
=
numneigh
[
i
];
for
(
int
jj
=
0
;
jj
<
jnum
;
jj
++
)
{
int
j
=
jlist
[
jj
];
double
factor_dpd
=
special_lj
[
sbmask
(
j
)];
j
&=
NEIGHMASK
;
double
delx
=
xtmp
-
x
[
j
][
0
];
double
dely
=
ytmp
-
x
[
j
][
1
];
double
delz
=
ztmp
-
x
[
j
][
2
];
double
rsq
=
delx
*
delx
+
dely
*
dely
+
delz
*
delz
;
int
jtype
=
type
[
j
];
if
(
rsq
<
cutsq
[
itype
][
jtype
])
{
double
r
=
sqrt
(
rsq
);
if
(
r
<
EPSILON
)
continue
;
// r can be 0.0 in DPD systems
double
rinv
=
1.0
/
r
;
double
delvx
=
vxtmp
-
v
[
j
][
0
];
double
delvy
=
vytmp
-
v
[
j
][
1
];
double
delvz
=
vztmp
-
v
[
j
][
2
];
double
dot
=
delx
*
delvx
+
dely
*
delvy
+
delz
*
delvz
;
double
wc
=
1.0
-
r
/
cut
[
itype
][
jtype
];
wc
=
MAX
(
0
,
MIN
(
1.0
,
wc
));
double
wr
=
pow
(
wc
,
0.5
*
power
[
itype
][
jtype
]);
double
randnum
=
random
->
gaussian
();
// conservative force = a0 * wc
// drag force = -gamma * wr^2 * (delx dot delv) / r
// random force = sigma * wr^(power/2) * rnd * dtinvsqrt;
double
fpair
=
a0
[
itype
][
jtype
]
*
wc
;
fpair
-=
gamma
[
itype
][
jtype
]
*
wr
*
wr
*
dot
*
rinv
;
fpair
+=
sigma
[
itype
][
jtype
]
*
wr
*
randnum
*
dtinvsqrt
;
fpair
*=
factor_dpd
*
rinv
;
f
[
i
][
0
]
+=
delx
*
fpair
;
f
[
i
][
1
]
+=
dely
*
fpair
;
f
[
i
][
2
]
+=
delz
*
fpair
;
// chemical concentration transport
if
(
r
<
cutcc
[
itype
][
jtype
])
{
for
(
int
k
=
0
;
k
<
cc_species
;
k
++
)
{
double
wcr
=
1.0
-
r
/
cutcc
[
itype
][
jtype
];
wcr
=
MAX
(
0
,
wcr
);
wcr
=
pow
(
wcr
,
0.5
*
powercc
[
itype
][
jtype
][
k
]);
double
randnum
=
random
->
gaussian
();
randnum
=
MAX
(
-
5.0
,
MIN
(
randnum
,
5.0
));
double
dQc
=
-
kappa
[
itype
][
jtype
][
k
]
*
wcr
*
wcr
*
(
cc
[
i
][
k
]
-
cc
[
j
][
k
]);
double
dQr
=
epsilon
[
itype
][
jtype
][
k
]
*
wcr
*
(
cc
[
i
][
k
]
+
cc
[
j
][
k
])
*
randnum
*
dtinvsqrt
;
cc_flux
[
i
][
k
]
+=
(
dQc
+
dQr
);
if
(
newton_pair
||
j
<
nlocal
)
cc_flux
[
j
][
k
]
-=
(
dQc
+
dQr
);
}
}
//-----------------------------------------------------------
if
(
newton_pair
||
j
<
nlocal
)
{
f
[
j
][
0
]
-=
delx
*
fpair
;
f
[
j
][
1
]
-=
dely
*
fpair
;
f
[
j
][
2
]
-=
delz
*
fpair
;
}
if
(
eflag
)
{
evdwl
=
0.5
*
a0
[
itype
][
jtype
]
*
cut
[
itype
][
jtype
]
*
wc
*
wc
;
evdwl
*=
factor_dpd
;
}
if
(
evflag
)
ev_tally
(
i
,
j
,
nlocal
,
newton_pair
,
evdwl
,
0.0
,
fpair
,
delx
,
dely
,
delz
);
}
}
}
if
(
vflag_fdotr
)
virial_fdotr_compute
();
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void
PairTDPD
::
allocate
()
{
int
i
,
j
;
allocated
=
1
;
int
n
=
atom
->
ntypes
;
memory
->
create
(
setflag
,
n
+
1
,
n
+
1
,
"pair:setflag"
);
for
(
i
=
1
;
i
<=
n
;
i
++
)
for
(
j
=
i
;
j
<=
n
;
j
++
)
setflag
[
i
][
j
]
=
0
;
memory
->
create
(
cutsq
,
n
+
1
,
n
+
1
,
"pair:cutsq"
);
memory
->
create
(
cut
,
n
+
1
,
n
+
1
,
"pair:cut"
);
memory
->
create
(
cutcc
,
n
+
1
,
n
+
1
,
"pair:cutcc"
);
memory
->
create
(
a0
,
n
+
1
,
n
+
1
,
"pair:a0"
);
memory
->
create
(
gamma
,
n
+
1
,
n
+
1
,
"pair:gamma"
);
memory
->
create
(
sigma
,
n
+
1
,
n
+
1
,
"pair:sigma"
);
memory
->
create
(
power
,
n
+
1
,
n
+
1
,
"pair:power"
);
memory
->
create
(
kappa
,
n
+
1
,
n
+
1
,
cc_species
,
"pair:kappa"
);
memory
->
create
(
epsilon
,
n
+
1
,
n
+
1
,
cc_species
,
"pair:epsilon"
);
memory
->
create
(
powercc
,
n
+
1
,
n
+
1
,
cc_species
,
"pair:powercc"
);
for
(
i
=
0
;
i
<=
atom
->
ntypes
;
i
++
)
for
(
j
=
0
;
j
<=
atom
->
ntypes
;
j
++
)
sigma
[
i
][
j
]
=
gamma
[
i
][
j
]
=
0.0
;
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void
PairTDPD
::
settings
(
int
narg
,
char
**
arg
)
{
if
(
narg
!=
3
)
error
->
all
(
FLERR
,
"Illegal pair_style command"
);
temperature
=
force
->
numeric
(
FLERR
,
arg
[
0
]);
cut_global
=
force
->
numeric
(
FLERR
,
arg
[
1
]);
seed
=
force
->
inumeric
(
FLERR
,
arg
[
2
]);
// initialize Marsaglia RNG with processor-unique seed
if
(
seed
<=
0
)
{
struct
timespec
time
;
clock_gettime
(
CLOCK_REALTIME
,
&
time
);
seed
=
time
.
tv_nsec
;
// if seed is non-positive, get the current time as the seed
}
delete
random
;
random
=
new
RanMars
(
lmp
,(
seed
+
comm
->
me
)
%
900000000
);
// reset cutoffs that have been explicitly set
if
(
allocated
)
{
int
i
,
j
;
for
(
i
=
1
;
i
<=
atom
->
ntypes
;
i
++
)
for
(
j
=
i
+
1
;
j
<=
atom
->
ntypes
;
j
++
)
if
(
setflag
[
i
][
j
])
cut
[
i
][
j
]
=
cut_global
;
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void
PairTDPD
::
coeff
(
int
narg
,
char
**
arg
)
{
if
(
narg
!=
7
+
3
*
cc_species
)
error
->
all
(
FLERR
,
"Incorrect args for pair tdpd coefficients"
);
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
);
double
a0_one
=
force
->
numeric
(
FLERR
,
arg
[
2
]);
double
gamma_one
=
force
->
numeric
(
FLERR
,
arg
[
3
]);
double
power_one
=
force
->
numeric
(
FLERR
,
arg
[
4
]);
double
cut_one
=
force
->
numeric
(
FLERR
,
arg
[
5
]);
double
cutcc_one
=
force
->
numeric
(
FLERR
,
arg
[
6
]);
double
kappa_one
[
cc_species
],
epsilon_one
[
cc_species
],
powercc_one
[
cc_species
];
for
(
int
k
=
0
;
k
<
cc_species
;
k
++
)
{
kappa_one
[
k
]
=
force
->
numeric
(
FLERR
,
arg
[
7
+
3
*
k
]);
epsilon_one
[
k
]
=
force
->
numeric
(
FLERR
,
arg
[
8
+
3
*
k
]);
powercc_one
[
k
]
=
force
->
numeric
(
FLERR
,
arg
[
9
+
3
*
k
]);
}
int
count
=
0
;
for
(
int
i
=
ilo
;
i
<=
ihi
;
i
++
)
for
(
int
j
=
MAX
(
jlo
,
i
);
j
<=
jhi
;
j
++
)
{
a0
[
i
][
j
]
=
a0_one
;
gamma
[
i
][
j
]
=
gamma_one
;
power
[
i
][
j
]
=
power_one
;
cut
[
i
][
j
]
=
cut_one
;
cutcc
[
i
][
j
]
=
cutcc_one
;
for
(
int
k
=
0
;
k
<
cc_species
;
k
++
)
{
kappa
[
i
][
j
][
k
]
=
kappa_one
[
k
];
epsilon
[
i
][
j
][
k
]
=
epsilon_one
[
k
];
powercc
[
i
][
j
][
k
]
=
powercc_one
[
k
];
}
setflag
[
i
][
j
]
=
1
;
count
++
;
}
if
(
count
==
0
)
error
->
all
(
FLERR
,
"Incorrect args for pair coefficients"
);
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void
PairTDPD
::
init_style
()
{
if
(
comm
->
ghost_velocity
==
0
)
error
->
all
(
FLERR
,
"Pair tdpd requires ghost atoms store velocity"
);
// if newton off, forces between atoms ij will be double computed
// using different random numbers
if
(
force
->
newton_pair
==
0
&&
comm
->
me
==
0
)
error
->
warning
(
FLERR
,
"Pair tdpd needs newton pair on for momentum conservation"
);
neighbor
->
request
(
this
,
instance_me
);
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double
PairTDPD
::
init_one
(
int
i
,
int
j
)
{
if
(
setflag
[
i
][
j
]
==
0
)
error
->
all
(
FLERR
,
"All pair coeffs are not set"
);
sigma
[
i
][
j
]
=
sqrt
(
2.0
*
force
->
boltz
*
temperature
*
gamma
[
i
][
j
]);
cut
[
j
][
i
]
=
cut
[
i
][
j
];
cutcc
[
j
][
i
]
=
cutcc
[
i
][
j
];
a0
[
j
][
i
]
=
a0
[
i
][
j
];
gamma
[
j
][
i
]
=
gamma
[
i
][
j
];
sigma
[
j
][
i
]
=
sigma
[
i
][
j
];
power
[
j
][
i
]
=
power
[
i
][
j
];
for
(
int
k
=
0
;
k
<
cc_species
;
k
++
)
{
kappa
[
j
][
i
][
k
]
=
kappa
[
i
][
j
][
k
];
epsilon
[
j
][
i
][
k
]
=
epsilon
[
i
][
j
][
k
];
powercc
[
j
][
i
][
k
]
=
powercc
[
i
][
j
][
k
];
}
return
cut
[
i
][
j
];
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void
PairTDPD
::
write_restart
(
FILE
*
fp
)
{
write_restart_settings
(
fp
);
for
(
int
i
=
1
;
i
<=
atom
->
ntypes
;
i
++
)
for
(
int
j
=
i
;
j
<=
atom
->
ntypes
;
j
++
)
{
fwrite
(
&
setflag
[
i
][
j
],
sizeof
(
int
),
1
,
fp
);
if
(
setflag
[
i
][
j
])
{
fwrite
(
&
a0
[
i
][
j
],
sizeof
(
double
),
1
,
fp
);
fwrite
(
&
gamma
[
i
][
j
],
sizeof
(
double
),
1
,
fp
);
fwrite
(
&
power
[
i
][
j
],
sizeof
(
double
),
1
,
fp
);
fwrite
(
&
cut
[
i
][
j
],
sizeof
(
double
),
1
,
fp
);
fwrite
(
&
cutcc
[
i
][
j
],
sizeof
(
double
),
1
,
fp
);
for
(
int
k
=
0
;
k
<
cc_species
;
k
++
)
{
fwrite
(
&
kappa
[
i
][
j
][
k
],
sizeof
(
double
),
1
,
fp
);
fwrite
(
&
epsilon
[
i
][
j
][
k
],
sizeof
(
double
),
1
,
fp
);
fwrite
(
&
powercc
[
i
][
j
][
k
],
sizeof
(
double
),
1
,
fp
);
}
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void
PairTDPD
::
read_restart
(
FILE
*
fp
)
{
read_restart_settings
(
fp
);
allocate
();
int
me
=
comm
->
me
;
for
(
int
i
=
1
;
i
<=
atom
->
ntypes
;
i
++
)
for
(
int
j
=
i
;
j
<=
atom
->
ntypes
;
j
++
)
{
if
(
me
==
0
)
fread
(
&
setflag
[
i
][
j
],
sizeof
(
int
),
1
,
fp
);
MPI_Bcast
(
&
setflag
[
i
][
j
],
1
,
MPI_INT
,
0
,
world
);
if
(
setflag
[
i
][
j
])
{
if
(
me
==
0
)
{
fread
(
&
a0
[
i
][
j
],
sizeof
(
double
),
1
,
fp
);
fread
(
&
gamma
[
i
][
j
],
sizeof
(
double
),
1
,
fp
);
fread
(
&
power
[
i
][
j
],
sizeof
(
double
),
1
,
fp
);
fread
(
&
cut
[
i
][
j
],
sizeof
(
double
),
1
,
fp
);
fread
(
&
cutcc
[
i
][
j
],
sizeof
(
double
),
1
,
fp
);
for
(
int
k
=
0
;
k
<
cc_species
;
k
++
)
{
fread
(
&
kappa
[
i
][
j
][
k
],
sizeof
(
double
),
1
,
fp
);
fread
(
&
epsilon
[
i
][
j
][
k
],
sizeof
(
double
),
1
,
fp
);
fread
(
&
powercc
[
i
][
j
][
k
],
sizeof
(
double
),
1
,
fp
);
}
}
MPI_Bcast
(
&
a0
[
i
][
j
],
1
,
MPI_DOUBLE
,
0
,
world
);
MPI_Bcast
(
&
gamma
[
i
][
j
],
1
,
MPI_DOUBLE
,
0
,
world
);
MPI_Bcast
(
&
power
[
i
][
j
],
1
,
MPI_DOUBLE
,
0
,
world
);
MPI_Bcast
(
&
cut
[
i
][
j
],
1
,
MPI_DOUBLE
,
0
,
world
);
MPI_Bcast
(
&
cutcc
[
i
][
j
],
1
,
MPI_DOUBLE
,
0
,
world
);
for
(
int
k
=
0
;
k
<
cc_species
;
k
++
)
{
MPI_Bcast
(
&
kappa
[
i
][
j
][
k
],
1
,
MPI_DOUBLE
,
0
,
world
);
MPI_Bcast
(
&
epsilon
[
i
][
j
][
k
],
1
,
MPI_DOUBLE
,
0
,
world
);
MPI_Bcast
(
&
powercc
[
i
][
j
][
k
],
1
,
MPI_DOUBLE
,
0
,
world
);
}
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void
PairTDPD
::
write_restart_settings
(
FILE
*
fp
)
{
fwrite
(
&
temperature
,
sizeof
(
double
),
1
,
fp
);
fwrite
(
&
cut_global
,
sizeof
(
double
),
1
,
fp
);
fwrite
(
&
seed
,
sizeof
(
int
),
1
,
fp
);
fwrite
(
&
mix_flag
,
sizeof
(
int
),
1
,
fp
);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void
PairTDPD
::
read_restart_settings
(
FILE
*
fp
)
{
if
(
comm
->
me
==
0
)
{
fread
(
&
temperature
,
sizeof
(
double
),
1
,
fp
);
fread
(
&
cut_global
,
sizeof
(
double
),
1
,
fp
);
fread
(
&
seed
,
sizeof
(
int
),
1
,
fp
);
fread
(
&
mix_flag
,
sizeof
(
int
),
1
,
fp
);
}
MPI_Bcast
(
&
temperature
,
1
,
MPI_DOUBLE
,
0
,
world
);
MPI_Bcast
(
&
cut_global
,
1
,
MPI_DOUBLE
,
0
,
world
);
MPI_Bcast
(
&
seed
,
1
,
MPI_INT
,
0
,
world
);
MPI_Bcast
(
&
mix_flag
,
1
,
MPI_INT
,
0
,
world
);
// initialize Marsaglia RNG with processor-unique seed
// same seed that pair_style command initially specified
if
(
random
)
delete
random
;
random
=
new
RanMars
(
lmp
,
seed
+
comm
->
me
);
}
/* ---------------------------------------------------------------------- */
double
PairTDPD
::
single
(
int
i
,
int
j
,
int
itype
,
int
jtype
,
double
rsq
,
double
factor_coul
,
double
factor_dpd
,
double
&
fforce
)
{
double
r
,
rinv
,
wc
,
phi
;
r
=
sqrt
(
rsq
);
if
(
r
<
EPSILON
)
{
fforce
=
0.0
;
return
0.0
;
}
rinv
=
1.0
/
r
;
wc
=
1.0
-
r
/
cut
[
itype
][
jtype
];
fforce
=
a0
[
itype
][
jtype
]
*
wc
*
factor_dpd
*
rinv
;
phi
=
0.5
*
a0
[
itype
][
jtype
]
*
cut
[
itype
][
jtype
]
*
wc
*
wc
;
return
factor_dpd
*
phi
;
}
Event Timeline
Log In to Comment