Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F69895550
pair_lj_sdk_coul_msm.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, Jul 4, 02:12
Size
10 KB
Mime Type
text/x-c++
Expires
Sat, Jul 6, 02:12 (2 d)
Engine
blob
Format
Raw Data
Handle
18746636
Attached To
rLAMMPS lammps
pair_lj_sdk_coul_msm.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: Axel Kohlmeyer (Temple U)
This style is a simplified re-implementation of the CG/CMM pair style
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pair_lj_sdk_coul_msm.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "kspace.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "update.h"
#include "integrate.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "lj_sdk_common.h"
using
namespace
LAMMPS_NS
;
using
namespace
MathConst
;
using
namespace
LJSDKParms
;
/* ---------------------------------------------------------------------- */
PairLJSDKCoulMSM
::
PairLJSDKCoulMSM
(
LAMMPS
*
lmp
)
:
PairLJSDKCoulLong
(
lmp
)
{
ewaldflag
=
pppmflag
=
0
;
msmflag
=
1
;
respa_enable
=
0
;
ftable
=
NULL
;
}
/* ---------------------------------------------------------------------- */
void
PairLJSDKCoulMSM
::
compute
(
int
eflag
,
int
vflag
)
{
if
(
force
->
kspace
->
scalar_pressure_flag
)
error
->
all
(
FLERR
,
"Must use 'kspace_modify pressure/scalar no' with Pair style"
);
if
(
eflag
||
vflag
)
{
ev_setup
(
eflag
,
vflag
);
}
else
evflag
=
vflag_fdotr
=
0
;
if
(
evflag
)
{
if
(
eflag
)
{
if
(
force
->
newton_pair
)
eval_msm
<
1
,
1
,
1
>
();
else
eval_msm
<
1
,
1
,
0
>
();
}
else
{
if
(
force
->
newton_pair
)
eval_msm
<
1
,
0
,
1
>
();
else
eval_msm
<
1
,
0
,
0
>
();
}
}
else
{
if
(
force
->
newton_pair
)
eval_msm
<
0
,
0
,
1
>
();
else
eval_msm
<
0
,
0
,
0
>
();
}
if
(
vflag_fdotr
)
virial_fdotr_compute
();
}
/* ---------------------------------------------------------------------- */
template
<
int
EVFLAG
,
int
EFLAG
,
int
NEWTON_PAIR
>
void
PairLJSDKCoulMSM
::
eval_msm
()
{
int
i
,
ii
,
j
,
jj
,
jtype
,
itable
;
double
qtmp
,
xtmp
,
ytmp
,
ztmp
,
delx
,
dely
,
delz
,
evdwl
,
ecoul
,
fpair
;
double
fraction
,
table
;
double
r
,
rsq
,
r2inv
,
forcecoul
,
forcelj
,
factor_coul
,
factor_lj
;
double
egamma
,
fgamma
,
prefactor
;
const
double
*
const
*
const
x
=
atom
->
x
;
double
*
const
*
const
f
=
atom
->
f
;
const
double
*
const
q
=
atom
->
q
;
const
int
*
const
type
=
atom
->
type
;
const
int
nlocal
=
atom
->
nlocal
;
const
double
*
const
special_coul
=
force
->
special_coul
;
const
double
*
const
special_lj
=
force
->
special_lj
;
const
double
qqrd2e
=
force
->
qqrd2e
;
double
fxtmp
,
fytmp
,
fztmp
;
const
int
inum
=
list
->
inum
;
const
int
*
const
ilist
=
list
->
ilist
;
const
int
*
const
numneigh
=
list
->
numneigh
;
const
int
*
const
*
const
firstneigh
=
list
->
firstneigh
;
// loop over neighbors of my atoms
for
(
ii
=
0
;
ii
<
inum
;
ii
++
)
{
i
=
ilist
[
ii
];
qtmp
=
q
[
i
];
xtmp
=
x
[
i
][
0
];
ytmp
=
x
[
i
][
1
];
ztmp
=
x
[
i
][
2
];
fxtmp
=
fytmp
=
fztmp
=
0.0
;
const
int
itype
=
type
[
i
];
const
int
*
const
jlist
=
firstneigh
[
i
];
const
int
jnum
=
numneigh
[
i
];
for
(
jj
=
0
;
jj
<
jnum
;
jj
++
)
{
forcecoul
=
forcelj
=
evdwl
=
ecoul
=
0.0
;
j
=
jlist
[
jj
];
factor_lj
=
special_lj
[
sbmask
(
j
)];
factor_coul
=
special_coul
[
sbmask
(
j
)];
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
])
{
r2inv
=
1.0
/
rsq
;
const
int
ljt
=
lj_type
[
itype
][
jtype
];
if
(
rsq
<
cut_coulsq
)
{
if
(
!
ncoultablebits
||
rsq
<=
tabinnersq
)
{
r
=
sqrt
(
rsq
);
prefactor
=
qqrd2e
*
qtmp
*
q
[
j
]
/
r
;
fgamma
=
1.0
+
(
rsq
/
cut_coulsq
)
*
force
->
kspace
->
dgamma
(
r
/
cut_coul
);
forcecoul
=
prefactor
*
fgamma
;
if
(
EFLAG
)
{
egamma
=
1.0
-
(
r
/
cut_coul
)
*
force
->
kspace
->
gamma
(
r
/
cut_coul
);
ecoul
=
prefactor
*
egamma
;
}
if
(
factor_coul
<
1.0
)
{
forcecoul
-=
(
1.0
-
factor_coul
)
*
prefactor
;
if
(
EFLAG
)
ecoul
-=
(
1.0
-
factor_coul
)
*
prefactor
;
}
}
else
{
union_int_float_t
rsq_lookup
;
rsq_lookup
.
f
=
rsq
;
itable
=
rsq_lookup
.
i
&
ncoulmask
;
itable
>>=
ncoulshiftbits
;
fraction
=
(
rsq_lookup
.
f
-
rtable
[
itable
])
*
drtable
[
itable
];
table
=
ftable
[
itable
]
+
fraction
*
dftable
[
itable
];
forcecoul
=
qtmp
*
q
[
j
]
*
table
;
if
(
EFLAG
)
ecoul
=
qtmp
*
q
[
j
]
*
(
etable
[
itable
]
+
fraction
*
detable
[
itable
]);
if
(
factor_coul
<
1.0
)
{
table
=
ctable
[
itable
]
+
fraction
*
dctable
[
itable
];
prefactor
=
qtmp
*
q
[
j
]
*
table
;
forcecoul
-=
(
1.0
-
factor_coul
)
*
prefactor
;
if
(
EFLAG
)
ecoul
-=
(
1.0
-
factor_coul
)
*
prefactor
;
}
}
}
else
forcecoul
=
0.0
;
if
(
rsq
<
cut_ljsq
[
itype
][
jtype
])
{
if
(
ljt
==
LJ12_4
)
{
const
double
r4inv
=
r2inv
*
r2inv
;
forcelj
=
r4inv
*
(
lj1
[
itype
][
jtype
]
*
r4inv
*
r4inv
-
lj2
[
itype
][
jtype
]);
if
(
EFLAG
)
evdwl
=
r4inv
*
(
lj3
[
itype
][
jtype
]
*
r4inv
*
r4inv
-
lj4
[
itype
][
jtype
])
-
offset
[
itype
][
jtype
];
}
else
if
(
ljt
==
LJ9_6
)
{
const
double
r3inv
=
r2inv
*
sqrt
(
r2inv
);
const
double
r6inv
=
r3inv
*
r3inv
;
forcelj
=
r6inv
*
(
lj1
[
itype
][
jtype
]
*
r3inv
-
lj2
[
itype
][
jtype
]);
if
(
EFLAG
)
evdwl
=
r6inv
*
(
lj3
[
itype
][
jtype
]
*
r3inv
-
lj4
[
itype
][
jtype
])
-
offset
[
itype
][
jtype
];
}
else
if
(
ljt
==
LJ12_6
)
{
const
double
r6inv
=
r2inv
*
r2inv
*
r2inv
;
forcelj
=
r6inv
*
(
lj1
[
itype
][
jtype
]
*
r6inv
-
lj2
[
itype
][
jtype
]);
if
(
EFLAG
)
evdwl
=
r6inv
*
(
lj3
[
itype
][
jtype
]
*
r6inv
-
lj4
[
itype
][
jtype
])
-
offset
[
itype
][
jtype
];
}
forcelj
*=
factor_lj
;
if
(
EFLAG
)
evdwl
*=
factor_lj
;
}
fpair
=
(
forcecoul
+
forcelj
)
*
r2inv
;
fxtmp
+=
delx
*
fpair
;
fytmp
+=
dely
*
fpair
;
fztmp
+=
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
,
evdwl
,
ecoul
,
fpair
,
delx
,
dely
,
delz
);
}
}
f
[
i
][
0
]
+=
fxtmp
;
f
[
i
][
1
]
+=
fytmp
;
f
[
i
][
2
]
+=
fztmp
;
}
}
/* ---------------------------------------------------------------------- */
double
PairLJSDKCoulMSM
::
single
(
int
i
,
int
j
,
int
itype
,
int
jtype
,
double
rsq
,
double
factor_coul
,
double
factor_lj
,
double
&
fforce
)
{
double
r2inv
,
r
,
fgamma
,
egamma
,
prefactor
;
double
fraction
,
table
,
forcecoul
,
forcelj
,
phicoul
,
philj
;
int
itable
;
forcecoul
=
forcelj
=
phicoul
=
philj
=
0.0
;
r2inv
=
1.0
/
rsq
;
if
(
rsq
<
cut_coulsq
)
{
if
(
!
ncoultablebits
||
rsq
<=
tabinnersq
)
{
r
=
sqrt
(
rsq
);
prefactor
=
force
->
qqrd2e
*
atom
->
q
[
i
]
*
atom
->
q
[
j
]
/
r
;
egamma
=
1.0
-
(
r
/
cut_coul
)
*
force
->
kspace
->
gamma
(
r
/
cut_coul
);
fgamma
=
1.0
+
(
rsq
/
cut_coulsq
)
*
force
->
kspace
->
dgamma
(
r
/
cut_coul
);
forcecoul
=
prefactor
*
fgamma
;
phicoul
=
prefactor
*
egamma
;
if
(
factor_coul
<
1.0
)
{
forcecoul
-=
(
1.0
-
factor_coul
)
*
prefactor
;
phicoul
-=
(
1.0
-
factor_coul
)
*
prefactor
;
}
}
else
{
union_int_float_t
rsq_lookup_single
;
rsq_lookup_single
.
f
=
rsq
;
itable
=
rsq_lookup_single
.
i
&
ncoulmask
;
itable
>>=
ncoulshiftbits
;
fraction
=
(
rsq_lookup_single
.
f
-
rtable
[
itable
])
*
drtable
[
itable
];
table
=
ftable
[
itable
]
+
fraction
*
dftable
[
itable
];
forcecoul
=
atom
->
q
[
i
]
*
atom
->
q
[
j
]
*
table
;
table
=
etable
[
itable
]
+
fraction
*
detable
[
itable
];
phicoul
=
atom
->
q
[
i
]
*
atom
->
q
[
j
]
*
table
;
if
(
factor_coul
<
1.0
)
{
table
=
ctable
[
itable
]
+
fraction
*
dctable
[
itable
];
prefactor
=
atom
->
q
[
i
]
*
atom
->
q
[
j
]
*
table
;
forcecoul
-=
(
1.0
-
factor_coul
)
*
prefactor
;
phicoul
-=
(
1.0
-
factor_coul
)
*
prefactor
;
}
}
}
if
(
rsq
<
cut_ljsq
[
itype
][
jtype
])
{
const
int
ljt
=
lj_type
[
itype
][
jtype
];
const
double
ljpow1
=
lj_pow1
[
ljt
];
const
double
ljpow2
=
lj_pow2
[
ljt
];
const
double
ljpref
=
lj_prefact
[
ljt
];
const
double
ratio
=
sigma
[
itype
][
jtype
]
/
sqrt
(
rsq
);
const
double
eps
=
epsilon
[
itype
][
jtype
];
forcelj
=
factor_lj
*
ljpref
*
eps
*
(
ljpow1
*
pow
(
ratio
,
ljpow1
)
-
ljpow2
*
pow
(
ratio
,
ljpow2
))
/
rsq
;
philj
=
factor_lj
*
(
ljpref
*
eps
*
(
pow
(
ratio
,
ljpow1
)
-
pow
(
ratio
,
ljpow2
))
-
offset
[
itype
][
jtype
]);
}
fforce
=
(
forcecoul
+
forcelj
)
*
r2inv
;
return
phicoul
+
philj
;
}
/* ---------------------------------------------------------------------- */
void
*
PairLJSDKCoulMSM
::
extract
(
const
char
*
str
,
int
&
dim
)
{
dim
=
2
;
if
(
strcmp
(
str
,
"epsilon"
)
==
0
)
return
(
void
*
)
epsilon
;
if
(
strcmp
(
str
,
"sigma"
)
==
0
)
return
(
void
*
)
sigma
;
if
(
strcmp
(
str
,
"lj_type"
)
==
0
)
return
(
void
*
)
lj_type
;
if
(
strcmp
(
str
,
"lj1"
)
==
0
)
return
(
void
*
)
lj1
;
if
(
strcmp
(
str
,
"lj2"
)
==
0
)
return
(
void
*
)
lj2
;
if
(
strcmp
(
str
,
"lj3"
)
==
0
)
return
(
void
*
)
lj3
;
if
(
strcmp
(
str
,
"lj4"
)
==
0
)
return
(
void
*
)
lj4
;
if
(
strcmp
(
str
,
"rminsq"
)
==
0
)
return
(
void
*
)
rminsq
;
if
(
strcmp
(
str
,
"emin"
)
==
0
)
return
(
void
*
)
emin
;
dim
=
0
;
if
(
strcmp
(
str
,
"cut_coul"
)
==
0
)
return
(
void
*
)
&
cut_coul
;
if
(
strcmp
(
str
,
"cut_msm"
)
==
0
)
return
(
void
*
)
&
cut_coul
;
return
NULL
;
}
Event Timeline
Log In to Comment