Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F80867540
pair_gw_zbl.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
Tue, Sep 3, 00:06
Size
9 KB
Mime Type
text/x-c
Expires
Thu, Sep 5, 00:06 (2 d)
Engine
blob
Format
Raw Data
Handle
20395178
Attached To
rLAMMPS lammps
pair_gw_zbl.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: German Samolyuk (ORNL)
Based on PairTersoffZBL by Aidan Thompson (SNL) and David Farrell (NWU)
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pair_gw_zbl.h"
#include "atom.h"
#include "update.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "force.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
#include "math_const.h"
using
namespace
LAMMPS_NS
;
using
namespace
MathConst
;
#define MAXLINE 1024
#define DELTA 4
/* ---------------------------------------------------------------------- */
PairGWZBL
::
PairGWZBL
(
LAMMPS
*
lmp
)
:
PairGW
(
lmp
)
{
// hard-wired constants in metal or real units
// a0 = Bohr radius
// epsilon0 = permittivity of vacuum = q / energy-distance units
// e = unit charge
// 1 Kcal/mole = 0.043365121 eV
if
(
strcmp
(
update
->
unit_style
,
"metal"
)
==
0
)
{
global_a_0
=
0.529
;
global_epsilon_0
=
0.00552635
;
global_e
=
1.0
;
}
else
if
(
strcmp
(
update
->
unit_style
,
"real"
)
==
0
)
{
global_a_0
=
0.529
;
global_epsilon_0
=
0.00552635
*
0.043365121
;
global_e
=
1.0
;
}
else
error
->
all
(
FLERR
,
"Pair gw/zbl requires metal or real units"
);
}
/* ---------------------------------------------------------------------- */
void
PairGWZBL
::
read_file
(
char
*
file
)
{
int
params_per_line
=
21
;
char
**
words
=
new
char
*
[
params_per_line
+
1
];
memory
->
sfree
(
params
);
params
=
NULL
;
nparams
=
maxparam
=
0
;
// open file on proc 0
FILE
*
fp
;
if
(
comm
->
me
==
0
)
{
fp
=
force
->
open_potential
(
file
);
if
(
fp
==
NULL
)
{
char
str
[
128
];
sprintf
(
str
,
"Cannot open GW potential file %s"
,
file
);
error
->
one
(
FLERR
,
str
);
}
}
// read each line out of file, skipping blank lines or leading '#'
// store line of params if all 3 element tags are in element list
int
n
,
nwords
,
ielement
,
jelement
,
kelement
;
char
line
[
MAXLINE
],
*
ptr
;
int
eof
=
0
;
while
(
1
)
{
if
(
comm
->
me
==
0
)
{
ptr
=
fgets
(
line
,
MAXLINE
,
fp
);
if
(
ptr
==
NULL
)
{
eof
=
1
;
fclose
(
fp
);
}
else
n
=
strlen
(
line
)
+
1
;
}
MPI_Bcast
(
&
eof
,
1
,
MPI_INT
,
0
,
world
);
if
(
eof
)
break
;
MPI_Bcast
(
&
n
,
1
,
MPI_INT
,
0
,
world
);
MPI_Bcast
(
line
,
n
,
MPI_CHAR
,
0
,
world
);
// strip comment, skip line if blank
if
((
ptr
=
strchr
(
line
,
'#'
)))
*
ptr
=
'\0'
;
nwords
=
atom
->
count_words
(
line
);
if
(
nwords
==
0
)
continue
;
// concatenate additional lines until have params_per_line words
while
(
nwords
<
params_per_line
)
{
n
=
strlen
(
line
);
if
(
comm
->
me
==
0
)
{
ptr
=
fgets
(
&
line
[
n
],
MAXLINE
-
n
,
fp
);
if
(
ptr
==
NULL
)
{
eof
=
1
;
fclose
(
fp
);
}
else
n
=
strlen
(
line
)
+
1
;
}
MPI_Bcast
(
&
eof
,
1
,
MPI_INT
,
0
,
world
);
if
(
eof
)
break
;
MPI_Bcast
(
&
n
,
1
,
MPI_INT
,
0
,
world
);
MPI_Bcast
(
line
,
n
,
MPI_CHAR
,
0
,
world
);
if
((
ptr
=
strchr
(
line
,
'#'
)))
*
ptr
=
'\0'
;
nwords
=
atom
->
count_words
(
line
);
}
if
(
nwords
!=
params_per_line
)
error
->
all
(
FLERR
,
"Incorrect format in GW potential file"
);
// words = ptrs to all words in line
nwords
=
0
;
words
[
nwords
++
]
=
strtok
(
line
,
"
\t\n\r\f
"
);
while
((
words
[
nwords
++
]
=
strtok
(
NULL
,
"
\t\n\r\f
"
)))
continue
;
// ielement,jelement,kelement = 1st args
// if all 3 args are in element list, then parse this line
// else skip to next line
for
(
ielement
=
0
;
ielement
<
nelements
;
ielement
++
)
if
(
strcmp
(
words
[
0
],
elements
[
ielement
])
==
0
)
break
;
if
(
ielement
==
nelements
)
continue
;
for
(
jelement
=
0
;
jelement
<
nelements
;
jelement
++
)
if
(
strcmp
(
words
[
1
],
elements
[
jelement
])
==
0
)
break
;
if
(
jelement
==
nelements
)
continue
;
for
(
kelement
=
0
;
kelement
<
nelements
;
kelement
++
)
if
(
strcmp
(
words
[
2
],
elements
[
kelement
])
==
0
)
break
;
if
(
kelement
==
nelements
)
continue
;
// load up parameter settings and error check their values
if
(
nparams
==
maxparam
)
{
maxparam
+=
DELTA
;
params
=
(
Param
*
)
memory
->
srealloc
(
params
,
maxparam
*
sizeof
(
Param
),
"pair:params"
);
}
params
[
nparams
].
ielement
=
ielement
;
params
[
nparams
].
jelement
=
jelement
;
params
[
nparams
].
kelement
=
kelement
;
params
[
nparams
].
powerm
=
atof
(
words
[
3
]);
params
[
nparams
].
gamma
=
atof
(
words
[
4
]);
params
[
nparams
].
lam3
=
atof
(
words
[
5
]);
params
[
nparams
].
c
=
atof
(
words
[
6
]);
params
[
nparams
].
d
=
atof
(
words
[
7
]);
params
[
nparams
].
h
=
atof
(
words
[
8
]);
params
[
nparams
].
powern
=
atof
(
words
[
9
]);
params
[
nparams
].
beta
=
atof
(
words
[
10
]);
params
[
nparams
].
lam2
=
atof
(
words
[
11
]);
params
[
nparams
].
bigb
=
atof
(
words
[
12
]);
params
[
nparams
].
bigr
=
atof
(
words
[
13
]);
params
[
nparams
].
bigd
=
atof
(
words
[
14
]);
params
[
nparams
].
lam1
=
atof
(
words
[
15
]);
params
[
nparams
].
biga
=
atof
(
words
[
16
]);
params
[
nparams
].
Z_i
=
atof
(
words
[
17
]);
params
[
nparams
].
Z_j
=
atof
(
words
[
18
]);
params
[
nparams
].
ZBLcut
=
atof
(
words
[
19
]);
params
[
nparams
].
ZBLexpscale
=
atof
(
words
[
20
]);
// currently only allow m exponent of 1 or 3
params
[
nparams
].
powermint
=
int
(
params
[
nparams
].
powerm
);
if
(
params
[
nparams
].
lam3
<
0.0
||
params
[
nparams
].
c
<
0.0
||
params
[
nparams
].
d
<
0.0
||
params
[
nparams
].
powern
<
0.0
||
params
[
nparams
].
beta
<
0.0
||
params
[
nparams
].
lam2
<
0.0
||
params
[
nparams
].
bigb
<
0.0
||
params
[
nparams
].
bigr
<
0.0
||
params
[
nparams
].
bigd
<
0.0
||
params
[
nparams
].
bigd
>
params
[
nparams
].
bigr
||
params
[
nparams
].
lam3
<
0.0
||
params
[
nparams
].
biga
<
0.0
||
params
[
nparams
].
powerm
-
params
[
nparams
].
powermint
!=
0.0
||
(
params
[
nparams
].
powermint
!=
3
&&
params
[
nparams
].
powermint
!=
1
)
||
params
[
nparams
].
gamma
<
0.0
||
params
[
nparams
].
Z_i
<
1.0
||
params
[
nparams
].
Z_j
<
1.0
||
params
[
nparams
].
ZBLcut
<
0.0
||
params
[
nparams
].
ZBLexpscale
<
0.0
)
error
->
all
(
FLERR
,
"Illegal GW parameter"
);
nparams
++
;
}
delete
[]
words
;
}
/* ---------------------------------------------------------------------- */
void
PairGWZBL
::
repulsive
(
Param
*
param
,
double
rsq
,
double
&
fforce
,
int
eflag
,
double
&
eng
)
{
double
r
,
tmp_fc
,
tmp_fc_d
,
tmp_exp
;
// GW repulsive portion
r
=
sqrt
(
rsq
);
tmp_fc
=
gw_fc
(
r
,
param
);
tmp_fc_d
=
gw_fc_d
(
r
,
param
);
tmp_exp
=
exp
(
-
param
->
lam1
*
r
);
double
fforce_gw
=
param
->
biga
*
tmp_exp
*
(
tmp_fc_d
-
tmp_fc
*
param
->
lam1
);
double
eng_gw
=
tmp_fc
*
param
->
biga
*
tmp_exp
;
// ZBL repulsive portion
double
esq
=
pow
(
global_e
,
2.0
);
double
a_ij
=
(
0.8854
*
global_a_0
)
/
(
pow
(
param
->
Z_i
,
0.23
)
+
pow
(
param
->
Z_j
,
0.23
));
double
premult
=
(
param
->
Z_i
*
param
->
Z_j
*
esq
)
/
(
4.0
*
MY_PI
*
global_epsilon_0
);
double
r_ov_a
=
r
/
a_ij
;
double
phi
=
0.1818
*
exp
(
-
3.2
*
r_ov_a
)
+
0.5099
*
exp
(
-
0.9423
*
r_ov_a
)
+
0.2802
*
exp
(
-
0.4029
*
r_ov_a
)
+
0.02817
*
exp
(
-
0.2016
*
r_ov_a
);
double
dphi
=
(
1.0
/
a_ij
)
*
(
-
3.2
*
0.1818
*
exp
(
-
3.2
*
r_ov_a
)
-
0.9423
*
0.5099
*
exp
(
-
0.9423
*
r_ov_a
)
-
0.4029
*
0.2802
*
exp
(
-
0.4029
*
r_ov_a
)
-
0.2016
*
0.02817
*
exp
(
-
0.2016
*
r_ov_a
));
double
fforce_ZBL
=
premult
*-
phi
/
rsq
+
premult
*
dphi
/
r
;
double
eng_ZBL
=
premult
*
(
1.0
/
r
)
*
phi
;
// combine two parts with smoothing by Fermi-like function
fforce
=
-
(
-
F_fermi_d
(
r
,
param
)
*
eng_ZBL
+
(
1.0
-
F_fermi
(
r
,
param
))
*
fforce_ZBL
+
F_fermi_d
(
r
,
param
)
*
eng_gw
+
F_fermi
(
r
,
param
)
*
fforce_gw
)
/
r
;
if
(
eflag
)
eng
=
(
1.0
-
F_fermi
(
r
,
param
))
*
eng_ZBL
+
F_fermi
(
r
,
param
)
*
eng_gw
;
}
/* ---------------------------------------------------------------------- */
double
PairGWZBL
::
gw_fa
(
double
r
,
Param
*
param
)
{
if
(
r
>
param
->
bigr
+
param
->
bigd
)
return
0.0
;
return
-
param
->
bigb
*
exp
(
-
param
->
lam2
*
r
)
*
gw_fc
(
r
,
param
)
*
F_fermi
(
r
,
param
);
}
/* ---------------------------------------------------------------------- */
double
PairGWZBL
::
gw_fa_d
(
double
r
,
Param
*
param
)
{
if
(
r
>
param
->
bigr
+
param
->
bigd
)
return
0.0
;
return
param
->
bigb
*
exp
(
-
param
->
lam2
*
r
)
*
(
param
->
lam2
*
gw_fc
(
r
,
param
)
*
F_fermi
(
r
,
param
)
-
gw_fc_d
(
r
,
param
)
*
F_fermi
(
r
,
param
)
-
gw_fc
(
r
,
param
)
*
F_fermi_d
(
r
,
param
));
}
/* ----------------------------------------------------------------------
Fermi-like smoothing function
------------------------------------------------------------------------- */
double
PairGWZBL
::
F_fermi
(
double
r
,
Param
*
param
)
{
return
1.0
/
(
1.0
+
exp
(
-
param
->
ZBLexpscale
*
(
r
-
param
->
ZBLcut
)));
}
/* ----------------------------------------------------------------------
Fermi-like smoothing function derivative with respect to r
------------------------------------------------------------------------- */
double
PairGWZBL
::
F_fermi_d
(
double
r
,
Param
*
param
)
{
return
param
->
ZBLexpscale
*
exp
(
-
param
->
ZBLexpscale
*
(
r
-
param
->
ZBLcut
))
/
pow
(
1.0
+
exp
(
-
param
->
ZBLexpscale
*
(
r
-
param
->
ZBLcut
)),
2.0
);
}
Event Timeline
Log In to Comment