Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F64179826
pair_hybrid.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
Sat, May 25, 03:07
Size
32 KB
Mime Type
text/x-c
Expires
Mon, May 27, 03:07 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
17864319
Attached To
rLAMMPS lammps
pair_hybrid.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.
------------------------------------------------------------------------- */
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include "pair_hybrid.h"
#include "atom.h"
#include "force.h"
#include "pair.h"
#include "neighbor.h"
#include "neigh_request.h"
#include "update.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
#include "respa.h"
using
namespace
LAMMPS_NS
;
/* ---------------------------------------------------------------------- */
PairHybrid
::
PairHybrid
(
LAMMPS
*
lmp
)
:
Pair
(
lmp
),
styles
(
NULL
),
keywords
(
NULL
),
multiple
(
NULL
),
nmap
(
NULL
),
map
(
NULL
),
special_lj
(
NULL
),
special_coul
(
NULL
)
{
nstyles
=
0
;
outerflag
=
0
;
respaflag
=
0
;
}
/* ---------------------------------------------------------------------- */
PairHybrid
::~
PairHybrid
()
{
if
(
nstyles
)
{
for
(
int
m
=
0
;
m
<
nstyles
;
m
++
)
{
delete
styles
[
m
];
delete
[]
keywords
[
m
];
if
(
special_lj
[
m
])
delete
[]
special_lj
[
m
];
if
(
special_coul
[
m
])
delete
[]
special_coul
[
m
];
}
}
delete
[]
styles
;
delete
[]
keywords
;
delete
[]
multiple
;
delete
[]
special_lj
;
delete
[]
special_coul
;
delete
[]
svector
;
if
(
allocated
)
{
memory
->
destroy
(
setflag
);
memory
->
destroy
(
cutsq
);
memory
->
destroy
(
cutghost
);
memory
->
destroy
(
nmap
);
memory
->
destroy
(
map
);
}
}
/* ----------------------------------------------------------------------
call each sub-style's compute() or compute_outer() function
accumulate sub-style global/peratom energy/virial in hybrid
for global vflag = 1:
each sub-style computes own virial[6]
sum sub-style virial[6] to hybrid's virial[6]
for global vflag = 2:
call sub-style with adjusted vflag to prevent it calling
virial_fdotr_compute()
hybrid calls virial_fdotr_compute() on final accumulated f
------------------------------------------------------------------------- */
void
PairHybrid
::
compute
(
int
eflag
,
int
vflag
)
{
int
i
,
j
,
m
,
n
;
// if no_virial_fdotr_compute is set and global component of
// incoming vflag = 2, then
// reset vflag as if global component were 1
// necessary since one or more sub-styles cannot compute virial as F dot r
if
(
no_virial_fdotr_compute
&&
vflag
%
4
==
2
)
vflag
=
1
+
vflag
/
4
*
4
;
if
(
eflag
||
vflag
)
ev_setup
(
eflag
,
vflag
);
else
evflag
=
vflag_fdotr
=
eflag_global
=
vflag_global
=
eflag_atom
=
vflag_atom
=
0
;
// check if global component of incoming vflag = 2
// if so, reset vflag passed to substyle as if it were 0
// necessary so substyle will not invoke virial_fdotr_compute()
int
vflag_substyle
;
if
(
vflag
%
4
==
2
)
vflag_substyle
=
vflag
/
4
*
4
;
else
vflag_substyle
=
vflag
;
double
*
saved_special
=
save_special
();
// check if we are running with r-RESPA using the hybrid keyword
Respa
*
respa
=
NULL
;
respaflag
=
0
;
if
(
strstr
(
update
->
integrate_style
,
"respa"
))
{
respa
=
(
Respa
*
)
update
->
integrate
;
if
(
respa
->
nhybrid_styles
>
0
)
respaflag
=
1
;
}
for
(
m
=
0
;
m
<
nstyles
;
m
++
)
{
set_special
(
m
);
if
(
!
respaflag
||
(
respaflag
&&
respa
->
hybrid_compute
[
m
]))
{
// invoke compute() unless compute flag is turned off or
// outerflag is set and sub-style has a compute_outer() method
if
(
styles
[
m
]
->
compute_flag
==
0
)
continue
;
if
(
outerflag
&&
styles
[
m
]
->
respa_enable
)
styles
[
m
]
->
compute_outer
(
eflag
,
vflag_substyle
);
else
styles
[
m
]
->
compute
(
eflag
,
vflag_substyle
);
}
restore_special
(
saved_special
);
// jump to next sub-style if r-RESPA does not want global accumulated data
if
(
respaflag
&&
!
respa
->
tally_global
)
continue
;
if
(
eflag_global
)
{
eng_vdwl
+=
styles
[
m
]
->
eng_vdwl
;
eng_coul
+=
styles
[
m
]
->
eng_coul
;
}
if
(
vflag_global
)
{
for
(
n
=
0
;
n
<
6
;
n
++
)
virial
[
n
]
+=
styles
[
m
]
->
virial
[
n
];
}
if
(
eflag_atom
)
{
n
=
atom
->
nlocal
;
if
(
force
->
newton_pair
)
n
+=
atom
->
nghost
;
double
*
eatom_substyle
=
styles
[
m
]
->
eatom
;
for
(
i
=
0
;
i
<
n
;
i
++
)
eatom
[
i
]
+=
eatom_substyle
[
i
];
}
if
(
vflag_atom
)
{
n
=
atom
->
nlocal
;
if
(
force
->
newton_pair
)
n
+=
atom
->
nghost
;
double
**
vatom_substyle
=
styles
[
m
]
->
vatom
;
for
(
i
=
0
;
i
<
n
;
i
++
)
for
(
j
=
0
;
j
<
6
;
j
++
)
vatom
[
i
][
j
]
+=
vatom_substyle
[
i
][
j
];
}
}
delete
[]
saved_special
;
if
(
vflag_fdotr
)
virial_fdotr_compute
();
}
/* ---------------------------------------------------------------------- */
void
PairHybrid
::
compute_inner
()
{
for
(
int
m
=
0
;
m
<
nstyles
;
m
++
)
if
(
styles
[
m
]
->
respa_enable
)
styles
[
m
]
->
compute_inner
();
}
/* ---------------------------------------------------------------------- */
void
PairHybrid
::
compute_middle
()
{
for
(
int
m
=
0
;
m
<
nstyles
;
m
++
)
if
(
styles
[
m
]
->
respa_enable
)
styles
[
m
]
->
compute_middle
();
}
/* ---------------------------------------------------------------------- */
void
PairHybrid
::
compute_outer
(
int
eflag
,
int
vflag
)
{
outerflag
=
1
;
compute
(
eflag
,
vflag
);
outerflag
=
0
;
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void
PairHybrid
::
allocate
()
{
allocated
=
1
;
int
n
=
atom
->
ntypes
;
memory
->
create
(
setflag
,
n
+
1
,
n
+
1
,
"pair:setflag"
);
for
(
int
i
=
1
;
i
<=
n
;
i
++
)
for
(
int
j
=
i
;
j
<=
n
;
j
++
)
setflag
[
i
][
j
]
=
0
;
memory
->
create
(
cutsq
,
n
+
1
,
n
+
1
,
"pair:cutsq"
);
memory
->
create
(
cutghost
,
n
+
1
,
n
+
1
,
"pair:cutghost"
);
memory
->
create
(
nmap
,
n
+
1
,
n
+
1
,
"pair:nmap"
);
memory
->
create
(
map
,
n
+
1
,
n
+
1
,
nstyles
,
"pair:map"
);
for
(
int
i
=
1
;
i
<=
n
;
i
++
)
for
(
int
j
=
i
;
j
<=
n
;
j
++
)
nmap
[
i
][
j
]
=
0
;
}
/* ----------------------------------------------------------------------
create one pair style for each arg in list
------------------------------------------------------------------------- */
void
PairHybrid
::
settings
(
int
narg
,
char
**
arg
)
{
if
(
narg
<
1
)
error
->
all
(
FLERR
,
"Illegal pair_style command"
);
// delete old lists, since cannot just change settings
if
(
nstyles
)
{
for
(
int
m
=
0
;
m
<
nstyles
;
m
++
)
delete
styles
[
m
];
delete
[]
styles
;
for
(
int
m
=
0
;
m
<
nstyles
;
m
++
)
delete
[]
keywords
[
m
];
delete
[]
keywords
;
}
if
(
allocated
)
{
memory
->
destroy
(
setflag
);
memory
->
destroy
(
cutsq
);
memory
->
destroy
(
cutghost
);
memory
->
destroy
(
nmap
);
memory
->
destroy
(
map
);
}
allocated
=
0
;
// allocate list of sub-styles as big as possibly needed if no extra args
styles
=
new
Pair
*
[
narg
];
keywords
=
new
char
*
[
narg
];
multiple
=
new
int
[
narg
];
special_lj
=
new
double
*
[
narg
];
special_coul
=
new
double
*
[
narg
];
// allocate each sub-style
// allocate uses suffix, but don't store suffix version in keywords,
// else syntax in coeff() will not match
// call settings() with set of args that are not pair style names
// use force->pair_map to determine which args these are
int
iarg
,
jarg
,
dummy
;
iarg
=
0
;
nstyles
=
0
;
while
(
iarg
<
narg
)
{
if
(
strcmp
(
arg
[
iarg
],
"hybrid"
)
==
0
)
error
->
all
(
FLERR
,
"Pair style hybrid cannot have hybrid as an argument"
);
if
(
strcmp
(
arg
[
iarg
],
"none"
)
==
0
)
error
->
all
(
FLERR
,
"Pair style hybrid cannot have none as an argument"
);
styles
[
nstyles
]
=
force
->
new_pair
(
arg
[
iarg
],
1
,
dummy
);
force
->
store_style
(
keywords
[
nstyles
],
arg
[
iarg
],
0
);
special_lj
[
nstyles
]
=
special_coul
[
nstyles
]
=
NULL
;
jarg
=
iarg
+
1
;
while
(
jarg
<
narg
&&
!
force
->
pair_map
->
count
(
arg
[
jarg
]))
jarg
++
;
styles
[
nstyles
]
->
settings
(
jarg
-
iarg
-
1
,
&
arg
[
iarg
+
1
]);
iarg
=
jarg
;
nstyles
++
;
}
// multiple[i] = 1 to M if sub-style used multiple times, else 0
for
(
int
i
=
0
;
i
<
nstyles
;
i
++
)
{
int
count
=
0
;
for
(
int
j
=
0
;
j
<
nstyles
;
j
++
)
{
if
(
strcmp
(
keywords
[
j
],
keywords
[
i
])
==
0
)
count
++
;
if
(
j
==
i
)
multiple
[
i
]
=
count
;
}
if
(
count
==
1
)
multiple
[
i
]
=
0
;
}
// set pair flags from sub-style flags
flags
();
}
/* ----------------------------------------------------------------------
set top-level pair flags from sub-style flags
------------------------------------------------------------------------- */
void
PairHybrid
::
flags
()
{
int
m
;
// set comm_forward, comm_reverse, comm_reverse_off to max of any sub-style
for
(
m
=
0
;
m
<
nstyles
;
m
++
)
{
if
(
styles
[
m
])
comm_forward
=
MAX
(
comm_forward
,
styles
[
m
]
->
comm_forward
);
if
(
styles
[
m
])
comm_reverse
=
MAX
(
comm_reverse
,
styles
[
m
]
->
comm_reverse
);
if
(
styles
[
m
])
comm_reverse_off
=
MAX
(
comm_reverse_off
,
styles
[
m
]
->
comm_reverse_off
);
}
// single_enable = 1 if any sub-style is set
// respa_enable = 1 if any sub-style is set
// manybody_flag = 1 if any sub-style is set
// no_virial_fdotr_compute = 1 if any sub-style is set
// ghostneigh = 1 if any sub-style is set
// ewaldflag, pppmflag, msmflag, dipoleflag, dispersionflag, tip4pflag = 1
// if any sub-style is set
// compute_flag = 1 if any sub-style is set
single_enable
=
0
;
compute_flag
=
0
;
for
(
m
=
0
;
m
<
nstyles
;
m
++
)
{
if
(
styles
[
m
]
->
single_enable
)
single_enable
=
1
;
if
(
styles
[
m
]
->
respa_enable
)
respa_enable
=
1
;
if
(
styles
[
m
]
->
manybody_flag
)
manybody_flag
=
1
;
if
(
styles
[
m
]
->
no_virial_fdotr_compute
)
no_virial_fdotr_compute
=
1
;
if
(
styles
[
m
]
->
ghostneigh
)
ghostneigh
=
1
;
if
(
styles
[
m
]
->
ewaldflag
)
ewaldflag
=
1
;
if
(
styles
[
m
]
->
pppmflag
)
pppmflag
=
1
;
if
(
styles
[
m
]
->
msmflag
)
msmflag
=
1
;
if
(
styles
[
m
]
->
dipoleflag
)
dipoleflag
=
1
;
if
(
styles
[
m
]
->
dispersionflag
)
dispersionflag
=
1
;
if
(
styles
[
m
]
->
tip4pflag
)
tip4pflag
=
1
;
if
(
styles
[
m
]
->
compute_flag
)
compute_flag
=
1
;
}
// single_extra = min of all sub-style single_extra
// allocate svector
single_extra
=
styles
[
0
]
->
single_extra
;
for
(
m
=
1
;
m
<
nstyles
;
m
++
)
single_extra
=
MIN
(
single_extra
,
styles
[
m
]
->
single_extra
);
if
(
single_extra
)
{
delete
[]
svector
;
svector
=
new
double
[
single_extra
];
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void
PairHybrid
::
coeff
(
int
narg
,
char
**
arg
)
{
if
(
narg
<
3
)
error
->
all
(
FLERR
,
"Incorrect args for pair 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
);
// 3rd arg = pair sub-style name
// 4th arg = pair sub-style index if name used multiple times
// allow for "none" as valid sub-style name
int
multflag
;
int
m
;
for
(
m
=
0
;
m
<
nstyles
;
m
++
)
{
multflag
=
0
;
if
(
strcmp
(
arg
[
2
],
keywords
[
m
])
==
0
)
{
if
(
multiple
[
m
])
{
multflag
=
1
;
if
(
narg
<
4
)
error
->
all
(
FLERR
,
"Incorrect args for pair coefficients"
);
if
(
!
isdigit
(
arg
[
3
][
0
]))
error
->
all
(
FLERR
,
"Incorrect args for pair coefficients"
);
int
index
=
force
->
inumeric
(
FLERR
,
arg
[
3
]);
if
(
index
==
multiple
[
m
])
break
;
else
continue
;
}
else
break
;
}
}
int
none
=
0
;
if
(
m
==
nstyles
)
{
if
(
strcmp
(
arg
[
2
],
"none"
)
==
0
)
none
=
1
;
else
error
->
all
(
FLERR
,
"Pair coeff for hybrid has invalid style"
);
}
// move 1st/2nd args to 2nd/3rd args
// if multflag: move 1st/2nd args to 3rd/4th args
// just copy ptrs, since arg[] points into original input line
arg
[
2
+
multflag
]
=
arg
[
1
];
arg
[
1
+
multflag
]
=
arg
[
0
];
// invoke sub-style coeff() starting with 1st remaining arg
if
(
!
none
)
styles
[
m
]
->
coeff
(
narg
-
1
-
multflag
,
&
arg
[
1
+
multflag
]);
// if sub-style only allows one pair coeff call (with * * and type mapping)
// then unset setflag/map assigned to that style before setting it below
// in case pair coeff for this sub-style is being called for 2nd time
if
(
!
none
&&
styles
[
m
]
->
one_coeff
)
for
(
int
i
=
1
;
i
<=
atom
->
ntypes
;
i
++
)
for
(
int
j
=
i
;
j
<=
atom
->
ntypes
;
j
++
)
if
(
nmap
[
i
][
j
]
&&
map
[
i
][
j
][
0
]
==
m
)
{
setflag
[
i
][
j
]
=
0
;
nmap
[
i
][
j
]
=
0
;
}
// set setflag and which type pairs map to which sub-style
// if sub-style is none: set hybrid setflag, wipe out map
// else: set hybrid setflag & map only if substyle setflag is set
// previous mappings are wiped out
int
count
=
0
;
for
(
int
i
=
ilo
;
i
<=
ihi
;
i
++
)
{
for
(
int
j
=
MAX
(
jlo
,
i
);
j
<=
jhi
;
j
++
)
{
if
(
none
)
{
setflag
[
i
][
j
]
=
1
;
nmap
[
i
][
j
]
=
0
;
count
++
;
}
else
if
(
styles
[
m
]
->
setflag
[
i
][
j
])
{
setflag
[
i
][
j
]
=
1
;
nmap
[
i
][
j
]
=
1
;
map
[
i
][
j
][
0
]
=
m
;
count
++
;
}
}
}
if
(
count
==
0
)
error
->
all
(
FLERR
,
"Incorrect args for pair coefficients"
);
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void
PairHybrid
::
init_style
()
{
int
i
,
m
,
itype
,
jtype
,
used
,
istyle
,
skip
;
// error if a sub-style is not used
int
ntypes
=
atom
->
ntypes
;
for
(
istyle
=
0
;
istyle
<
nstyles
;
istyle
++
)
{
used
=
0
;
for
(
itype
=
1
;
itype
<=
ntypes
;
itype
++
)
for
(
jtype
=
itype
;
jtype
<=
ntypes
;
jtype
++
)
for
(
m
=
0
;
m
<
nmap
[
itype
][
jtype
];
m
++
)
if
(
map
[
itype
][
jtype
][
m
]
==
istyle
)
used
=
1
;
if
(
used
==
0
)
error
->
all
(
FLERR
,
"Pair hybrid sub-style is not used"
);
}
// check if special_lj/special_coul overrides are compatible
for
(
istyle
=
0
;
istyle
<
nstyles
;
istyle
++
)
{
if
(
special_lj
[
istyle
])
{
for
(
i
=
1
;
i
<
4
;
++
i
)
{
if
(((
force
->
special_lj
[
i
]
==
0.0
)
||
(
force
->
special_lj
[
i
]
==
1.0
))
&&
(
force
->
special_lj
[
i
]
!=
special_lj
[
istyle
][
i
]))
error
->
all
(
FLERR
,
"Pair_modify special setting for pair hybrid "
"incompatible with global special_bonds setting"
);
}
}
if
(
special_coul
[
istyle
])
{
for
(
i
=
1
;
i
<
4
;
++
i
)
{
if
(((
force
->
special_coul
[
i
]
==
0.0
)
||
(
force
->
special_coul
[
i
]
==
1.0
))
&&
(
force
->
special_coul
[
i
]
!=
special_coul
[
istyle
][
i
]))
error
->
all
(
FLERR
,
"Pair_modify special setting for pair hybrid "
"incompatible with global special_bonds setting"
);
}
}
}
// each sub-style makes its neighbor list request(s)
for
(
istyle
=
0
;
istyle
<
nstyles
;
istyle
++
)
styles
[
istyle
]
->
init_style
();
// create skip lists inside each pair neigh request
// any kind of list can have its skip flag set at this stage
for
(
i
=
0
;
i
<
neighbor
->
nrequest
;
i
++
)
{
if
(
!
neighbor
->
requests
[
i
]
->
pair
)
continue
;
// istyle = associated sub-style for that request
for
(
istyle
=
0
;
istyle
<
nstyles
;
istyle
++
)
if
(
styles
[
istyle
]
==
neighbor
->
requests
[
i
]
->
requestor
)
break
;
// allocate iskip and ijskip
// initialize so as to skip all pair types
// set ijskip = 0 if type pair matches any entry in sub-style map
// set ijskip = 0 if mixing will assign type pair to this sub-style
// will occur if type pair is currently unassigned
// and both I,I and J,J are assigned to single sub-style
// and sub-style for both I,I and J,J match istyle
// set iskip = 1 only if all ijskip for itype are 1
int
*
iskip
=
new
int
[
ntypes
+
1
];
int
**
ijskip
;
memory
->
create
(
ijskip
,
ntypes
+
1
,
ntypes
+
1
,
"pair_hybrid:ijskip"
);
for
(
itype
=
1
;
itype
<=
ntypes
;
itype
++
)
for
(
jtype
=
1
;
jtype
<=
ntypes
;
jtype
++
)
ijskip
[
itype
][
jtype
]
=
1
;
for
(
itype
=
1
;
itype
<=
ntypes
;
itype
++
)
for
(
jtype
=
itype
;
jtype
<=
ntypes
;
jtype
++
)
{
for
(
m
=
0
;
m
<
nmap
[
itype
][
jtype
];
m
++
)
if
(
map
[
itype
][
jtype
][
m
]
==
istyle
)
ijskip
[
itype
][
jtype
]
=
ijskip
[
jtype
][
itype
]
=
0
;
if
(
nmap
[
itype
][
jtype
]
==
0
&&
nmap
[
itype
][
itype
]
==
1
&&
map
[
itype
][
itype
][
0
]
==
istyle
&&
nmap
[
jtype
][
jtype
]
==
1
&&
map
[
jtype
][
jtype
][
0
]
==
istyle
)
ijskip
[
itype
][
jtype
]
=
ijskip
[
jtype
][
itype
]
=
0
;
}
for
(
itype
=
1
;
itype
<=
ntypes
;
itype
++
)
{
iskip
[
itype
]
=
1
;
for
(
jtype
=
1
;
jtype
<=
ntypes
;
jtype
++
)
if
(
ijskip
[
itype
][
jtype
]
==
0
)
iskip
[
itype
]
=
0
;
}
// if any skipping occurs
// set request->skip and copy iskip and ijskip into request
// else delete iskip and ijskip
// no skipping if pair style assigned to all type pairs
skip
=
0
;
for
(
itype
=
1
;
itype
<=
ntypes
;
itype
++
)
for
(
jtype
=
1
;
jtype
<=
ntypes
;
jtype
++
)
if
(
ijskip
[
itype
][
jtype
]
==
1
)
skip
=
1
;
if
(
skip
)
{
neighbor
->
requests
[
i
]
->
skip
=
1
;
neighbor
->
requests
[
i
]
->
iskip
=
iskip
;
neighbor
->
requests
[
i
]
->
ijskip
=
ijskip
;
}
else
{
delete
[]
iskip
;
memory
->
destroy
(
ijskip
);
}
}
// combine sub-style neigh list requests and create new ones if needed
modify_requests
();
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double
PairHybrid
::
init_one
(
int
i
,
int
j
)
{
// if I,J is not set explicitly:
// perform mixing only if I,I sub-style = J,J sub-style
// also require I,I and J,J are both assigned to single sub-style
if
(
setflag
[
i
][
j
]
==
0
)
{
if
(
nmap
[
i
][
i
]
!=
1
||
nmap
[
j
][
j
]
!=
1
||
map
[
i
][
i
][
0
]
!=
map
[
j
][
j
][
0
])
error
->
one
(
FLERR
,
"All pair coeffs are not set"
);
nmap
[
i
][
j
]
=
1
;
map
[
i
][
j
][
0
]
=
map
[
i
][
i
][
0
];
}
// call init/mixing for all sub-styles of I,J
// set cutsq in sub-style just as Pair::init() does via call to init_one()
// set cutghost for I,J and J,I just as sub-style does
// sum tail corrections for I,J
// return max cutoff of all sub-styles assigned to I,J
// if no sub-styles assigned to I,J (pair_coeff none), cutmax = 0.0 returned
double
cutmax
=
0.0
;
cutghost
[
i
][
j
]
=
cutghost
[
j
][
i
]
=
0.0
;
if
(
tail_flag
)
etail_ij
=
ptail_ij
=
0.0
;
nmap
[
j
][
i
]
=
nmap
[
i
][
j
];
for
(
int
k
=
0
;
k
<
nmap
[
i
][
j
];
k
++
)
{
map
[
j
][
i
][
k
]
=
map
[
i
][
j
][
k
];
double
cut
=
styles
[
map
[
i
][
j
][
k
]]
->
init_one
(
i
,
j
);
styles
[
map
[
i
][
j
][
k
]]
->
cutsq
[
i
][
j
]
=
styles
[
map
[
i
][
j
][
k
]]
->
cutsq
[
j
][
i
]
=
cut
*
cut
;
if
(
styles
[
map
[
i
][
j
][
k
]]
->
ghostneigh
)
cutghost
[
i
][
j
]
=
cutghost
[
j
][
i
]
=
MAX
(
cutghost
[
i
][
j
],
styles
[
map
[
i
][
j
][
k
]]
->
cutghost
[
i
][
j
]);
if
(
tail_flag
)
{
etail_ij
+=
styles
[
map
[
i
][
j
][
k
]]
->
etail_ij
;
ptail_ij
+=
styles
[
map
[
i
][
j
][
k
]]
->
ptail_ij
;
}
cutmax
=
MAX
(
cutmax
,
cut
);
}
return
cutmax
;
}
/* ----------------------------------------------------------------------
invoke setup for each sub-style
------------------------------------------------------------------------- */
void
PairHybrid
::
setup
()
{
for
(
int
m
=
0
;
m
<
nstyles
;
m
++
)
styles
[
m
]
->
setup
();
}
/* ----------------------------------------------------------------------
examine sub-style neigh list requests
create new parent requests if needed, to derive sub-style requests from
------------------------------------------------------------------------- */
void
PairHybrid
::
modify_requests
()
{
int
i
,
j
;
NeighRequest
*
irq
,
*
jrq
;
// loop over pair requests only, including those added during looping
int
nrequest_original
=
neighbor
->
nrequest
;
for
(
i
=
0
;
i
<
neighbor
->
nrequest
;
i
++
)
{
if
(
!
neighbor
->
requests
[
i
]
->
pair
)
continue
;
// nothing more to do if this request:
// is not a skip list
// is a copy or half_from_full or granhistory list
// copy list setup is from pair style = hybrid/overlay
// which invokes this method at end of its modify_requests()
// if granhistory, turn off skip, since each gran sub-style
// its own history list, parent gran list does not have history
// if half_from_full, turn off skip, since it will derive
// from its full parent and its skip status
irq
=
neighbor
->
requests
[
i
];
if
(
irq
->
skip
==
0
)
continue
;
if
(
irq
->
copy
)
continue
;
if
(
irq
->
granhistory
||
irq
->
half_from_full
)
{
irq
->
skip
=
0
;
continue
;
}
// look for another list that matches via same_kind() and is not a skip list
// if one exists, point at that one via otherlist
// else make new parent request via copy_request() and point at that one
// new parent list is not a skip list
// parent does not need its ID set, since pair hybrid does not use it
for
(
j
=
0
;
j
<
neighbor
->
nrequest
;
j
++
)
{
if
(
!
neighbor
->
requests
[
j
]
->
pair
)
continue
;
jrq
=
neighbor
->
requests
[
j
];
if
(
irq
->
same_kind
(
jrq
)
&&
jrq
->
skip
==
0
)
break
;
}
if
(
j
<
neighbor
->
nrequest
)
irq
->
otherlist
=
j
;
else
{
int
newrequest
=
neighbor
->
request
(
this
,
instance_me
);
neighbor
->
requests
[
newrequest
]
->
copy_request
(
irq
);
irq
->
otherlist
=
newrequest
;
}
// for rRESPA inner/middle lists,
// which just created or set otherlist to parent:
// unset skip flag and otherlist
// this prevents neighbor from treating them as skip lists
if
(
irq
->
respainner
||
irq
->
respamiddle
)
{
irq
->
skip
=
0
;
irq
->
otherlist
=
-
1
;
}
}
// adjustments to newly added granular parent requests (gran = 1)
// set parent newton = 2 if has children with granonesided = 0 and 1
// else newton = 0 = setting of children
// if 2, also set child off2on for both granonesided kinds of children
// set parent gran onesided = 0 if has children with granonesided = 0 and 1
// else onesided = setting of children
for
(
i
=
nrequest_original
;
i
<
neighbor
->
nrequest
;
i
++
)
{
if
(
!
neighbor
->
requests
[
i
]
->
pair
)
continue
;
if
(
!
neighbor
->
requests
[
i
]
->
gran
)
continue
;
irq
=
neighbor
->
requests
[
i
];
int
onesided
=
-
1
;
for
(
j
=
0
;
j
<
nrequest_original
;
j
++
)
{
if
(
!
neighbor
->
requests
[
j
]
->
pair
)
continue
;
if
(
!
neighbor
->
requests
[
j
]
->
gran
)
continue
;
if
(
neighbor
->
requests
[
j
]
->
otherlist
!=
i
)
continue
;
jrq
=
neighbor
->
requests
[
j
];
if
(
onesided
<
0
)
onesided
=
jrq
->
granonesided
;
else
if
(
onesided
!=
jrq
->
granonesided
)
onesided
=
2
;
if
(
onesided
==
2
)
break
;
}
if
(
onesided
==
2
)
{
irq
->
newton
=
2
;
irq
->
granonesided
=
0
;
for
(
j
=
0
;
j
<
nrequest_original
;
j
++
)
{
if
(
!
neighbor
->
requests
[
j
]
->
pair
)
continue
;
if
(
!
neighbor
->
requests
[
j
]
->
gran
)
continue
;
if
(
neighbor
->
requests
[
j
]
->
otherlist
!=
i
)
continue
;
jrq
=
neighbor
->
requests
[
j
];
jrq
->
off2on
=
1
;
}
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void
PairHybrid
::
write_restart
(
FILE
*
fp
)
{
fwrite
(
&
nstyles
,
sizeof
(
int
),
1
,
fp
);
// each sub-style writes its settings, but no coeff info
int
n
;
for
(
int
m
=
0
;
m
<
nstyles
;
m
++
)
{
n
=
strlen
(
keywords
[
m
])
+
1
;
fwrite
(
&
n
,
sizeof
(
int
),
1
,
fp
);
fwrite
(
keywords
[
m
],
sizeof
(
char
),
n
,
fp
);
styles
[
m
]
->
write_restart_settings
(
fp
);
// write out per style special settings, if present
n
=
(
special_lj
[
m
]
==
NULL
)
?
0
:
1
;
fwrite
(
&
n
,
sizeof
(
int
),
1
,
fp
);
if
(
n
)
fwrite
(
special_lj
[
m
],
sizeof
(
double
),
4
,
fp
);
n
=
(
special_coul
[
m
]
==
NULL
)
?
0
:
1
;
fwrite
(
&
n
,
sizeof
(
int
),
1
,
fp
);
if
(
n
)
fwrite
(
special_coul
[
m
],
sizeof
(
double
),
4
,
fp
);
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void
PairHybrid
::
read_restart
(
FILE
*
fp
)
{
int
me
=
comm
->
me
;
if
(
me
==
0
)
fread
(
&
nstyles
,
sizeof
(
int
),
1
,
fp
);
MPI_Bcast
(
&
nstyles
,
1
,
MPI_INT
,
0
,
world
);
// allocate list of sub-styles
styles
=
new
Pair
*
[
nstyles
];
keywords
=
new
char
*
[
nstyles
];
multiple
=
new
int
[
nstyles
];
special_lj
=
new
double
*
[
nstyles
];
special_coul
=
new
double
*
[
nstyles
];
// each sub-style is created via new_pair()
// each reads its settings, but no coeff info
int
n
,
dummy
;
for
(
int
m
=
0
;
m
<
nstyles
;
m
++
)
{
if
(
me
==
0
)
fread
(
&
n
,
sizeof
(
int
),
1
,
fp
);
MPI_Bcast
(
&
n
,
1
,
MPI_INT
,
0
,
world
);
keywords
[
m
]
=
new
char
[
n
];
if
(
me
==
0
)
fread
(
keywords
[
m
],
sizeof
(
char
),
n
,
fp
);
MPI_Bcast
(
keywords
[
m
],
n
,
MPI_CHAR
,
0
,
world
);
styles
[
m
]
=
force
->
new_pair
(
keywords
[
m
],
0
,
dummy
);
styles
[
m
]
->
read_restart_settings
(
fp
);
// read back per style special settings, if present
special_lj
[
m
]
=
special_coul
[
m
]
=
NULL
;
if
(
me
==
0
)
fread
(
&
n
,
sizeof
(
int
),
1
,
fp
);
MPI_Bcast
(
&
n
,
1
,
MPI_INT
,
0
,
world
);
if
(
n
>
0
)
{
special_lj
[
m
]
=
new
double
[
4
];
if
(
me
==
0
)
fread
(
special_lj
[
m
],
sizeof
(
double
),
4
,
fp
);
MPI_Bcast
(
special_lj
[
m
],
4
,
MPI_DOUBLE
,
0
,
world
);
}
if
(
me
==
0
)
fread
(
&
n
,
sizeof
(
int
),
1
,
fp
);
MPI_Bcast
(
&
n
,
1
,
MPI_INT
,
0
,
world
);
if
(
n
>
0
)
{
special_coul
[
m
]
=
new
double
[
4
];
if
(
me
==
0
)
fread
(
special_coul
[
m
],
sizeof
(
double
),
4
,
fp
);
MPI_Bcast
(
special_coul
[
m
],
4
,
MPI_DOUBLE
,
0
,
world
);
}
}
// multiple[i] = 1 to M if sub-style used multiple times, else 0
for
(
int
i
=
0
;
i
<
nstyles
;
i
++
)
{
int
count
=
0
;
for
(
int
j
=
0
;
j
<
nstyles
;
j
++
)
{
if
(
strcmp
(
keywords
[
j
],
keywords
[
i
])
==
0
)
count
++
;
if
(
j
==
i
)
multiple
[
i
]
=
count
;
}
if
(
count
==
1
)
multiple
[
i
]
=
0
;
}
// set pair flags from sub-style flags
flags
();
}
/* ----------------------------------------------------------------------
call sub-style to compute single interaction
error if sub-style does not support single() call
since overlay could have multiple sub-styles, sum results explicitly
------------------------------------------------------------------------- */
double
PairHybrid
::
single
(
int
i
,
int
j
,
int
itype
,
int
jtype
,
double
rsq
,
double
factor_coul
,
double
factor_lj
,
double
&
fforce
)
{
if
(
nmap
[
itype
][
jtype
]
==
0
)
error
->
one
(
FLERR
,
"Invoked pair single on pair style none"
);
double
fone
;
fforce
=
0.0
;
double
esum
=
0.0
;
for
(
int
m
=
0
;
m
<
nmap
[
itype
][
jtype
];
m
++
)
{
if
(
rsq
<
styles
[
map
[
itype
][
jtype
][
m
]]
->
cutsq
[
itype
][
jtype
])
{
if
(
styles
[
map
[
itype
][
jtype
][
m
]]
->
single_enable
==
0
)
error
->
one
(
FLERR
,
"Pair hybrid sub-style does not support single call"
);
if
((
special_lj
[
map
[
itype
][
jtype
][
m
]]
!=
NULL
)
||
(
special_coul
[
map
[
itype
][
jtype
][
m
]]
!=
NULL
))
error
->
one
(
FLERR
,
"Pair hybrid single calls do not support"
" per sub-style special bond values"
);
esum
+=
styles
[
map
[
itype
][
jtype
][
m
]]
->
single
(
i
,
j
,
itype
,
jtype
,
rsq
,
factor_coul
,
factor_lj
,
fone
);
fforce
+=
fone
;
// copy substyle extra values into hybrid's svector
if
(
single_extra
&&
styles
[
map
[
itype
][
jtype
][
m
]]
->
single_extra
)
for
(
m
=
0
;
m
<
single_extra
;
m
++
)
svector
[
m
]
=
styles
[
map
[
itype
][
jtype
][
m
]]
->
svector
[
m
];
}
}
return
esum
;
}
/* ----------------------------------------------------------------------
modify parameters of the pair style and its sub-styles
------------------------------------------------------------------------- */
void
PairHybrid
::
modify_params
(
int
narg
,
char
**
arg
)
{
if
(
narg
==
0
)
error
->
all
(
FLERR
,
"Illegal pair_modify command"
);
// if 1st keyword is pair, apply other keywords to one sub-style
if
(
strcmp
(
arg
[
0
],
"pair"
)
==
0
)
{
if
(
narg
<
2
)
error
->
all
(
FLERR
,
"Illegal pair_modify command"
);
int
m
;
for
(
m
=
0
;
m
<
nstyles
;
m
++
)
if
(
strcmp
(
arg
[
1
],
keywords
[
m
])
==
0
)
break
;
if
(
m
==
nstyles
)
error
->
all
(
FLERR
,
"Unknown pair_modify hybrid sub-style"
);
int
iarg
=
2
;
if
(
multiple
[
m
])
{
if
(
narg
<
3
)
error
->
all
(
FLERR
,
"Illegal pair_modify command"
);
int
multiflag
=
force
->
inumeric
(
FLERR
,
arg
[
2
]);
for
(
m
=
0
;
m
<
nstyles
;
m
++
)
if
(
strcmp
(
arg
[
1
],
keywords
[
m
])
==
0
&&
multiflag
==
multiple
[
m
])
break
;
if
(
m
==
nstyles
)
error
->
all
(
FLERR
,
"Unknown pair_modify hybrid sub-style"
);
iarg
=
3
;
}
// if 2nd keyword (after pair) is special:
// invoke modify_special() for the sub-style
if
(
iarg
<
narg
&&
strcmp
(
arg
[
iarg
],
"special"
)
==
0
)
{
if
(
narg
<
iarg
+
5
)
error
->
all
(
FLERR
,
"Illegal pair_modify special command"
);
modify_special
(
m
,
narg
-
iarg
,
&
arg
[
iarg
+
1
]);
iarg
+=
5
;
}
// apply the remaining keywords to the base pair style itself and the
// sub-style except for "pair" and "special".
// the former is important for some keywords like "tail" or "compute"
if
(
narg
-
iarg
>
0
)
{
Pair
::
modify_params
(
narg
-
iarg
,
&
arg
[
iarg
]);
styles
[
m
]
->
modify_params
(
narg
-
iarg
,
&
arg
[
iarg
]);
}
// apply all keywords to pair hybrid itself and every sub-style
}
else
{
Pair
::
modify_params
(
narg
,
arg
);
for
(
int
m
=
0
;
m
<
nstyles
;
m
++
)
styles
[
m
]
->
modify_params
(
narg
,
arg
);
}
}
/* ----------------------------------------------------------------------
store a local per pair style override for special_lj and special_coul
------------------------------------------------------------------------- */
void
PairHybrid
::
modify_special
(
int
m
,
int
narg
,
char
**
arg
)
{
double
special
[
4
];
int
i
;
special
[
0
]
=
1.0
;
special
[
1
]
=
force
->
numeric
(
FLERR
,
arg
[
1
]);
special
[
2
]
=
force
->
numeric
(
FLERR
,
arg
[
2
]);
special
[
3
]
=
force
->
numeric
(
FLERR
,
arg
[
3
]);
if
(
strcmp
(
arg
[
0
],
"lj/coul"
)
==
0
)
{
if
(
!
special_lj
[
m
])
special_lj
[
m
]
=
new
double
[
4
];
if
(
!
special_coul
[
m
])
special_coul
[
m
]
=
new
double
[
4
];
for
(
i
=
0
;
i
<
4
;
++
i
)
special_lj
[
m
][
i
]
=
special_coul
[
m
][
i
]
=
special
[
i
];
}
else
if
(
strcmp
(
arg
[
0
],
"lj"
)
==
0
)
{
if
(
!
special_lj
[
m
])
special_lj
[
m
]
=
new
double
[
4
];
for
(
i
=
0
;
i
<
4
;
++
i
)
special_lj
[
m
][
i
]
=
special
[
i
];
}
else
if
(
strcmp
(
arg
[
0
],
"coul"
)
==
0
)
{
if
(
!
special_coul
[
m
])
special_coul
[
m
]
=
new
double
[
4
];
for
(
i
=
0
;
i
<
4
;
++
i
)
special_coul
[
m
][
i
]
=
special
[
i
];
}
else
error
->
all
(
FLERR
,
"Illegal pair_modify special command"
);
}
/* ----------------------------------------------------------------------
override global special bonds settings with per substyle values
------------------------------------------------------------------------- */
void
PairHybrid
::
set_special
(
int
m
)
{
int
i
;
if
(
special_lj
[
m
])
for
(
i
=
0
;
i
<
4
;
++
i
)
force
->
special_lj
[
i
]
=
special_lj
[
m
][
i
];
if
(
special_coul
[
m
])
for
(
i
=
0
;
i
<
4
;
++
i
)
force
->
special_coul
[
i
]
=
special_coul
[
m
][
i
];
}
/* ----------------------------------------------------------------------
store global special settings
------------------------------------------------------------------------- */
double
*
PairHybrid
::
save_special
()
{
double
*
saved
=
new
double
[
8
];
for
(
int
i
=
0
;
i
<
4
;
++
i
)
{
saved
[
i
]
=
force
->
special_lj
[
i
];
saved
[
i
+
4
]
=
force
->
special_coul
[
i
];
}
return
saved
;
}
/* ----------------------------------------------------------------------
restore global special settings from saved data
------------------------------------------------------------------------- */
void
PairHybrid
::
restore_special
(
double
*
saved
)
{
for
(
int
i
=
0
;
i
<
4
;
++
i
)
{
force
->
special_lj
[
i
]
=
saved
[
i
];
force
->
special_coul
[
i
]
=
saved
[
i
+
4
];
}
}
/* ----------------------------------------------------------------------
extract a ptr to a particular quantity stored by pair
pass request thru to sub-styles
return first non-NULL result except for cut_coul request
for cut_coul, insure all non-NULL results are equal since required by Kspace
------------------------------------------------------------------------- */
void
*
PairHybrid
::
extract
(
const
char
*
str
,
int
&
dim
)
{
void
*
cutptr
=
NULL
;
void
*
ptr
;
double
cutvalue
=
0.0
;
for
(
int
m
=
0
;
m
<
nstyles
;
m
++
)
{
ptr
=
styles
[
m
]
->
extract
(
str
,
dim
);
if
(
ptr
&&
strcmp
(
str
,
"cut_coul"
)
==
0
)
{
double
*
p_newvalue
=
(
double
*
)
ptr
;
double
newvalue
=
*
p_newvalue
;
if
(
cutptr
&&
newvalue
!=
cutvalue
)
error
->
all
(
FLERR
,
"Coulomb cutoffs of pair hybrid sub-styles do not match"
);
cutptr
=
ptr
;
cutvalue
=
newvalue
;
}
else
if
(
ptr
)
return
ptr
;
}
if
(
strcmp
(
str
,
"cut_coul"
)
==
0
)
return
cutptr
;
return
NULL
;
}
/* ---------------------------------------------------------------------- */
void
PairHybrid
::
reset_dt
()
{
for
(
int
m
=
0
;
m
<
nstyles
;
m
++
)
styles
[
m
]
->
reset_dt
();
}
/* ----------------------------------------------------------------------
check if itype,jtype maps to sub-style
------------------------------------------------------------------------- */
int
PairHybrid
::
check_ijtype
(
int
itype
,
int
jtype
,
char
*
substyle
)
{
for
(
int
m
=
0
;
m
<
nmap
[
itype
][
jtype
];
m
++
)
if
(
strcmp
(
keywords
[
map
[
itype
][
jtype
][
m
]],
substyle
)
==
0
)
return
1
;
return
0
;
}
/* ----------------------------------------------------------------------
memory usage of each sub-style
------------------------------------------------------------------------- */
double
PairHybrid
::
memory_usage
()
{
double
bytes
=
maxeatom
*
sizeof
(
double
);
bytes
+=
maxvatom
*
6
*
sizeof
(
double
);
for
(
int
m
=
0
;
m
<
nstyles
;
m
++
)
bytes
+=
styles
[
m
]
->
memory_usage
();
return
bytes
;
}
Event Timeline
Log In to Comment