added LD potential and wrote html-style doc
This commit is contained in:
tanmoy.7989
BIN
doc/src/Eqs/pair_local_density_energy.jpg
Normal file
BIN
doc/src/Eqs/pair_local_density_energy.jpg
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 3.0 KiB |
11
doc/src/Eqs/pair_local_density_energy.tex
Normal file
11
doc/src/Eqs/pair_local_density_energy.tex
Normal file
@ -0,0 +1,11 @@
|
|||||||
|
\documentclass[12pt]{article}
|
||||||
|
|
||||||
|
\begin{document}
|
||||||
|
|
||||||
|
$$
|
||||||
|
U_{LD} = \sum_i F(\rho_i)
|
||||||
|
$$
|
||||||
|
|
||||||
|
|
||||||
|
\end{document}
|
||||||
|
~
|
||||||
BIN
doc/src/Eqs/pair_local_density_energy_implement.jpg
Normal file
BIN
doc/src/Eqs/pair_local_density_energy_implement.jpg
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 7.8 KiB |
9
doc/src/Eqs/pair_local_density_energy_implement.tex
Normal file
9
doc/src/Eqs/pair_local_density_energy_implement.tex
Normal file
@ -0,0 +1,9 @@
|
|||||||
|
\documentclass[12pt]{article}
|
||||||
|
|
||||||
|
\begin{document}
|
||||||
|
|
||||||
|
$$
|
||||||
|
U_{LD} = \sum_k U_{LD}^{(k)} = \sum_i \left[ \sum_k a_\alpha^{(k)} F^{(k)} \left(\rho_i^{(k)}\right) \right]
|
||||||
|
$$
|
||||||
|
|
||||||
|
\end{document}
|
||||||
BIN
doc/src/Eqs/pair_local_density_energy_multi.jpg
Normal file
BIN
doc/src/Eqs/pair_local_density_energy_multi.jpg
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 3.4 KiB |
9
doc/src/Eqs/pair_local_density_energy_multi.tex
Normal file
9
doc/src/Eqs/pair_local_density_energy_multi.tex
Normal file
@ -0,0 +1,9 @@
|
|||||||
|
\documentclass[12pt]{article}
|
||||||
|
|
||||||
|
\begin{document}
|
||||||
|
|
||||||
|
$$
|
||||||
|
U_{LD} = \sum_i a_\alpha F(\rho_i)
|
||||||
|
$$
|
||||||
|
|
||||||
|
\end{document}
|
||||||
BIN
doc/src/Eqs/pair_local_density_indicator_func.jpg
Normal file
BIN
doc/src/Eqs/pair_local_density_indicator_func.jpg
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 8.8 KiB |
16
doc/src/Eqs/pair_local_density_indicator_func.tex
Normal file
16
doc/src/Eqs/pair_local_density_indicator_func.tex
Normal file
@ -0,0 +1,16 @@
|
|||||||
|
\documentclass[12pt]{article}
|
||||||
|
\usepackage[utf8]{inputenc}
|
||||||
|
\usepackage{amsmath}
|
||||||
|
\usepackage{amsfonts}
|
||||||
|
|
||||||
|
\begin{document}
|
||||||
|
\[
|
||||||
|
\varphi(r) =
|
||||||
|
\begin{cases}
|
||||||
|
1 & r \le R_1 \\
|
||||||
|
c_0 + c_2r^2 + c_4r^4 + c_6r^6 & r \in (R_1, R_2) \\
|
||||||
|
0 & r \ge R_2
|
||||||
|
\end{cases}
|
||||||
|
\]
|
||||||
|
|
||||||
|
\end{document}
|
||||||
BIN
doc/src/Eqs/pair_local_density_ld.jpg
Normal file
BIN
doc/src/Eqs/pair_local_density_ld.jpg
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 3.0 KiB |
10
doc/src/Eqs/pair_local_density_ld.tex
Normal file
10
doc/src/Eqs/pair_local_density_ld.tex
Normal file
@ -0,0 +1,10 @@
|
|||||||
|
\documentclass[12pt]{article}
|
||||||
|
|
||||||
|
\begin{document}
|
||||||
|
|
||||||
|
|
||||||
|
$$
|
||||||
|
\rho_i = \sum_{j \neq i} \varphi(r_{ij})
|
||||||
|
$$
|
||||||
|
|
||||||
|
\end{document}
|
||||||
BIN
doc/src/Eqs/pair_local_density_ld_implement.jpg
Normal file
BIN
doc/src/Eqs/pair_local_density_ld_implement.jpg
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 4.2 KiB |
10
doc/src/Eqs/pair_local_density_ld_implement.tex
Normal file
10
doc/src/Eqs/pair_local_density_ld_implement.tex
Normal file
@ -0,0 +1,10 @@
|
|||||||
|
\documentstyle[12pt]{article}
|
||||||
|
|
||||||
|
\begin{document}
|
||||||
|
|
||||||
|
|
||||||
|
$$
|
||||||
|
\rho_i^{(k)} = \sum_j b_\beta^{(k)} \varphi^{(k)} (r_{ij})
|
||||||
|
$$
|
||||||
|
|
||||||
|
\end{document}
|
||||||
BIN
doc/src/Eqs/pair_local_density_ld_multi.jpg
Normal file
BIN
doc/src/Eqs/pair_local_density_ld_multi.jpg
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 3.4 KiB |
10
doc/src/Eqs/pair_local_density_ld_multi.tex
Normal file
10
doc/src/Eqs/pair_local_density_ld_multi.tex
Normal file
@ -0,0 +1,10 @@
|
|||||||
|
\documentclass[12pt]{article}
|
||||||
|
|
||||||
|
\begin{document}
|
||||||
|
|
||||||
|
|
||||||
|
$$
|
||||||
|
\rho_i = \sum_{j \neq i} b_\beta \varphi(r_{ij})
|
||||||
|
$$
|
||||||
|
|
||||||
|
\end{document}
|
||||||
207
doc/src/pair_local_density.txt
Normal file
207
doc/src/pair_local_density.txt
Normal file
@ -0,0 +1,207 @@
|
|||||||
|
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
|
||||||
|
|
||||||
|
:link(lws,http://lammps.sandia.gov)
|
||||||
|
:link(ld,Manual.html)
|
||||||
|
:link(lc,Commands_all.html)
|
||||||
|
|
||||||
|
:line
|
||||||
|
|
||||||
|
pair_style local/density command :h3
|
||||||
|
|
||||||
|
[Syntax:]
|
||||||
|
|
||||||
|
pair_style style arg :pre
|
||||||
|
|
||||||
|
style = {local/density}
|
||||||
|
arg = name of file containing tabulated values of local density and the potential :ul
|
||||||
|
|
||||||
|
[Examples:]
|
||||||
|
|
||||||
|
pair_style local/density benzene_water.localdensity.table :pre
|
||||||
|
|
||||||
|
pair_style hybrid/overlay table spline 500 local/density
|
||||||
|
pair_coeff * * local/density benzene_water.localdensity.table :pre
|
||||||
|
|
||||||
|
[Description:]
|
||||||
|
|
||||||
|
The local density (LD) potential is a new potential style that is, in some
|
||||||
|
sense,a generalization of embedded atom models (EAM). The name "local density
|
||||||
|
potential" arises from the fact that it assigns an energy to an atom depending
|
||||||
|
on the number of neighboring atoms of given type around it within a predefined
|
||||||
|
spherical volume (i.e., within a cutoff). The bottom-up coarse-graining (CG)
|
||||||
|
literature sugggests that such potentials can be widely useful in capturing
|
||||||
|
effective multibody forces in a computationally efficient manner so as to
|
||||||
|
improve the quality of CG models of implicit solvation"(Sanyal1)"_#Sanyal1 and
|
||||||
|
phase-segregation in liquid mixtures"(Sanyal2)"_#Sanyal2, and provide guidelines
|
||||||
|
to determine the extent of manybody correlations present in a CG
|
||||||
|
model."(Rosenberger)"_#Rosenberger The LD potential in LAMMPS is primarily
|
||||||
|
intended to be used as a corrective potential over traditional pair potentials
|
||||||
|
in bottom-up CG models, i.e., as a hybrid pair style with
|
||||||
|
other explicit pair interaction terms (e.g., table spline, Lennard Jones, etc.).
|
||||||
|
Because the LD potential is not a pair potential per se, it is implemented
|
||||||
|
simply as a single auxiliary file with all specifications that will be read
|
||||||
|
upon initialization.
|
||||||
|
|
||||||
|
NOTE: Thus when used as the only interaction in the system, there is no
|
||||||
|
corresponding pair_coeff command and when used with other pair styles using the
|
||||||
|
hybrid/overlay option, the corresponding pair_coeff command must be supplied
|
||||||
|
* * as placeholders for the atomtypes.
|
||||||
|
|
||||||
|
:line
|
||||||
|
|
||||||
|
[System with a single CG atom type:]
|
||||||
|
|
||||||
|
A system of a single atom type (e.g., LJ argon) with a single local density (LD)
|
||||||
|
potential would have an energy given by:
|
||||||
|
|
||||||
|
:c,image(Eqs/pair_local_density_energy.jpg)
|
||||||
|
|
||||||
|
where rho_i is the LD at atom i and F(rho) is similar in spirit to the
|
||||||
|
embedding function used in EAM potentials. The LD at atom i is given by the sum
|
||||||
|
|
||||||
|
:c,image(Eqs/pair_local_density_ld.jpg)
|
||||||
|
|
||||||
|
where phi is an indicator function that is one at r=0 and zero beyond a cutoff
|
||||||
|
distance R2. The choice of the functional form of phi is somewhat arbitrary,
|
||||||
|
but the following piecewise cubic function has proven sufficiently general:
|
||||||
|
"(Sanyal1)"_#Sanyal1, "(Sanyal2)"_#Sanyal2 "(Rosenberger)"_#Rosenberger
|
||||||
|
|
||||||
|
:c,image(Eqs/pair_local_density_indicator_func.jpg)
|
||||||
|
|
||||||
|
The constants {c} are chosen so that the indicator function smoothly
|
||||||
|
interpolates between 1 and 0 between the distances R1 and R2, which are
|
||||||
|
called the inner and outer cutoffs, respectively. Thus phi satisfies
|
||||||
|
phi(R1) = 1, phi(R2) = dphi/dr @ (r=R1) = dphi/dr @ (r=R2) = 0. The embedding
|
||||||
|
function F(rho) may or may not have a closed-form expression. To maintain
|
||||||
|
generality, it is practically represented with a spline-interpolated table
|
||||||
|
over a predetermined range of rho. Outside of that range it simply adopts zero
|
||||||
|
values at the endpoints.
|
||||||
|
|
||||||
|
It can be shown that the total force between two atoms due to the LD potential
|
||||||
|
takes the form of a pair force, which motivates its designation as a LAMMPS
|
||||||
|
pair style. Please see "(Sanyal1)"_#Sanyal1 for details of the derivation.
|
||||||
|
|
||||||
|
:line
|
||||||
|
|
||||||
|
[Systems with arbitrary numbers of atom types:]
|
||||||
|
|
||||||
|
The potential is easily generalized to systems involving multiple atom types:
|
||||||
|
|
||||||
|
:c,image(Eqs/pair_local_density_energy_multi.jpg)
|
||||||
|
|
||||||
|
with the LD expressed as
|
||||||
|
|
||||||
|
:c,image(Eqs/pair_local_density_ld_multi.jpg)
|
||||||
|
|
||||||
|
where alpha gives the type of atom i, beta the type of atom j, and the
|
||||||
|
coefficients a and b filter for atom types as specified by the user. a is
|
||||||
|
called the central atom filter as it determines to which atoms the
|
||||||
|
potential applies; a_alpha = 1 if the LD potential applies to atom type alpha
|
||||||
|
else zero. On the other hand, b is called the neighbor atom filter because it
|
||||||
|
specifies which atom types to use in the calculation of the LD; b_beta = 1 if
|
||||||
|
atom type beta contributes to the LD and zero otherwise.
|
||||||
|
|
||||||
|
NOTE: Note that the potentials need not be symmetric with respect to atom types,
|
||||||
|
which is the reason for two distinct sets of coefficients a and b. An atom type
|
||||||
|
may contribute to the LD but not the potential, or to the potential but not the
|
||||||
|
LD. Such decisions are made by the user and should (ideally) be motivated on
|
||||||
|
physical grounds for the problem at hand.
|
||||||
|
|
||||||
|
:line
|
||||||
|
|
||||||
|
[General form for implementation in LAMMPS:]
|
||||||
|
|
||||||
|
Of course, a system with many atom types may have many different possible LD
|
||||||
|
potentials, each with their own atom type filters, cutoffs, and embedding
|
||||||
|
functions. The most general form of this potential as implemented in the
|
||||||
|
pair_style local/density is:
|
||||||
|
|
||||||
|
:c,image(Eqs/pair_local_density_energy_implement.jpg)
|
||||||
|
|
||||||
|
where, k is an index that spans the (arbitrary) number of applied LD potentials
|
||||||
|
N_LD. Each LD is calculated as before with:
|
||||||
|
|
||||||
|
:c,image(Eqs/pair_local_density_ld_implement.jpg)
|
||||||
|
|
||||||
|
The superscript on the indicator function phi simply indicates that it is
|
||||||
|
associated with specific values of the cutoff distances R1(k) and R2(k). In
|
||||||
|
summary, there may be N_LD distinct LD potentials. With each potential type (k),
|
||||||
|
one must specify:
|
||||||
|
|
||||||
|
the inner and outer cutoffs as R1 and R2
|
||||||
|
the central type filter a(k), where k = 1,2,...N_LD
|
||||||
|
the neighbor type filter b(k), where k = 1,2,...N_LD
|
||||||
|
the LD potential function F(k)(rho), typically as a table that is later spline-interpolated :ul
|
||||||
|
|
||||||
|
:line
|
||||||
|
|
||||||
|
[Tabulated input file format:]
|
||||||
|
|
||||||
|
Line 1: comment or blank (ignored)
|
||||||
|
Line 2: comment or blank (ignored)
|
||||||
|
Line 3: N_LD N_rho (# of LD potentials and # of tabulated values, single space separated)
|
||||||
|
Line 4: blank (ignored)
|
||||||
|
Line 5: R1(k) R2(k) (lower and upper cutoffs, single space separated)
|
||||||
|
Line 6: central-types (central atom types, single space separated)
|
||||||
|
Line 7: neighbor-types (neighbor atom types single space separated)
|
||||||
|
Line 8: rho_min rho_max drho (min, max and diff. in tabulated rho values, single space separated)
|
||||||
|
Line 9: F(k)(rho_min + 0.drho)
|
||||||
|
Line 10: F(k)(rho_min + 1.drho)
|
||||||
|
Line 11: F(k)(rho_min + 2.drho)
|
||||||
|
............
|
||||||
|
Line 9+N_rho: F(k)(rho_min + N_rho . drho)
|
||||||
|
Line 10+N_rho: blank (ignored) :ul
|
||||||
|
|
||||||
|
Block 2 :ul
|
||||||
|
|
||||||
|
Block 3 :ul
|
||||||
|
|
||||||
|
Block N_LD :ul
|
||||||
|
|
||||||
|
Lines 5 to 9+N_rho constitute the first block. Thus the input file is separated
|
||||||
|
(by blank lines) into N_LD blocks each representing a separate LD potential and
|
||||||
|
each specifying its own upper and lower cutoffs, central and neighbor atoms,
|
||||||
|
and potential. In general, blank lines anywhere are ignored.
|
||||||
|
|
||||||
|
:line
|
||||||
|
|
||||||
|
[Mixing, shift, table, tail correction, restart, info]:
|
||||||
|
This pair style does not support automatic mixing. For atom type pairs alpha,
|
||||||
|
beta and alpha != beta, even if LD potentials of type (alpha, alpha) and
|
||||||
|
(beta, beta) are provided, you will need to explicitly provide LD potential
|
||||||
|
types (alpha, beta) and (beta, alpha) if need be (Here, the notation (alpha,
|
||||||
|
beta) means that alpha is the central atom to which the LD potential is applied
|
||||||
|
and beta is the neighbor atom which contributes to the LD potential on alpha).
|
||||||
|
|
||||||
|
This pair style does not support the "pair_modify"_pair_modify.html
|
||||||
|
shift, table, and tail options.
|
||||||
|
|
||||||
|
The local/density pair style does not write its information to "binary restart
|
||||||
|
files"_restart.html, since it is stored in tabulated potential files.
|
||||||
|
Thus, you need to re-specify the pair_style and pair_coeff commands in
|
||||||
|
an input script that reads a restart file.
|
||||||
|
|
||||||
|
:line
|
||||||
|
|
||||||
|
[Restrictions:]
|
||||||
|
|
||||||
|
The local/density pair style is a part of the USER-MISC package. It is only
|
||||||
|
enabled if LAMMPS was built with that package. See the "Build
|
||||||
|
package"_Build_package.html doc page for more info.
|
||||||
|
|
||||||
|
[Related commands:]
|
||||||
|
|
||||||
|
"pair_coeff"_pair_coeff.html
|
||||||
|
|
||||||
|
[Default:] none
|
||||||
|
|
||||||
|
:line
|
||||||
|
|
||||||
|
|
||||||
|
:link(Sanyal1)
|
||||||
|
[(Sanyal1)] Sanyal and Shell, Journal of Chemical Physics, 2016, 145 (3), 034109.
|
||||||
|
:link(Sanyal2)
|
||||||
|
[(Sanyal2)] Sanyal and Shell, Journal of Physical Chemistry B, 122 (21), 5678-5693.
|
||||||
|
|
||||||
|
:link(Rosenberger)
|
||||||
|
[(Rosenberger)] Rosenberger, Sanyal, Shell and van der Vegt, Journal of Chemical Physics, 2019, 151 (4), 044111.
|
||||||
1406
examples/USER/misc/local_density/benzene_water/benzene_water.data
Normal file
1406
examples/USER/misc/local_density/benzene_water/benzene_water.data
Normal file
File diff suppressed because it is too large
Load Diff
@ -0,0 +1,62 @@
|
|||||||
|
# LAMMPS input file for 26.5% benzene mole fraction solution
|
||||||
|
# with 380 benzene and 1000 water molecules,
|
||||||
|
# using all possible local density potentials
|
||||||
|
# between benzene and water
|
||||||
|
#
|
||||||
|
# Author: Tanmoy Sanyal, Shell Group, UC Santa Barbara
|
||||||
|
#
|
||||||
|
# Refer: Sanyal and Shell, JPC-B, 2018, 122 (21), 5678-5693
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
# Initialize simulation box
|
||||||
|
dimension 3
|
||||||
|
boundary p p p
|
||||||
|
units real
|
||||||
|
atom_style molecular
|
||||||
|
|
||||||
|
# Set potential styles
|
||||||
|
pair_style hybrid/overlay table spline 500 local/density
|
||||||
|
|
||||||
|
# Read molecule data and set initial velocities
|
||||||
|
read_data benzene_water.data
|
||||||
|
velocity all create 3.0000e+02 16611 rot yes dist gaussian
|
||||||
|
|
||||||
|
# Assign potentials
|
||||||
|
pair_coeff 1 1 table benzene_water.pair.table PairBB
|
||||||
|
pair_coeff 1 2 table benzene_water.pair.table PairWW
|
||||||
|
pair_coeff 2 2 table benzene_water.pair.table PairBW
|
||||||
|
pair_coeff * * local/density benzene_water.localdensity.table
|
||||||
|
|
||||||
|
# Recentering during minimization and equilibration
|
||||||
|
fix recentering all recenter 0.0 0.0 0.0 units box
|
||||||
|
|
||||||
|
# Thermostat & time integration
|
||||||
|
timestep 2.0
|
||||||
|
thermo 100
|
||||||
|
thermo_style custom temp ke pe etotal ebond eangle edihed evdwl
|
||||||
|
|
||||||
|
# Minimization
|
||||||
|
minimize 1.e-4 0.0 10000 10000
|
||||||
|
|
||||||
|
# Set up integration parameters
|
||||||
|
fix timeintegration all nve
|
||||||
|
fix thermostat all langevin 3.0000e+02 3.0000e+02 1.0000e+02 81890
|
||||||
|
|
||||||
|
# Equilibration (for realistic results, run for 5000000 steps)
|
||||||
|
reset_timestep 0
|
||||||
|
run 5000
|
||||||
|
|
||||||
|
# Turn off recentering during production phase
|
||||||
|
unfix recentering
|
||||||
|
|
||||||
|
# Setup trajectory output
|
||||||
|
dump myDump all custom 100 benzene_water.lammpstrj.gz id type x y z element
|
||||||
|
dump_modify myDump element B W
|
||||||
|
dump_modify myDump sort id
|
||||||
|
|
||||||
|
# Production (for realistic results, run for 10000000 steps)
|
||||||
|
reset_timestep 0
|
||||||
|
run 1000
|
||||||
|
|
||||||
|
|
||||||
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
267
examples/USER/misc/local_density/benzene_water/log.04Sep19.g++.1
Normal file
267
examples/USER/misc/local_density/benzene_water/log.04Sep19.g++.1
Normal file
@ -0,0 +1,267 @@
|
|||||||
|
LAMMPS (7 Aug 2019)
|
||||||
|
# LAMMPS input file for 26.5% benzene mole fraction solution
|
||||||
|
# with 380 benzene and 1000 water molecules,
|
||||||
|
# using all possible local density potentials
|
||||||
|
# between benzene and water
|
||||||
|
#
|
||||||
|
# Author: Tanmoy Sanyal, Shell Group, UC Santa Barbara
|
||||||
|
#
|
||||||
|
# Refer: Sanyal and Shell, JPC-B, 2018, 122 (21), 5678-5693
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
# Initialize simulation box
|
||||||
|
dimension 3
|
||||||
|
boundary p p p
|
||||||
|
units real
|
||||||
|
atom_style molecular
|
||||||
|
|
||||||
|
# Set potential styles
|
||||||
|
pair_style hybrid/overlay table spline 500 local/density
|
||||||
|
|
||||||
|
# Read molecule data and set initial velocities
|
||||||
|
read_data benzene_water.data
|
||||||
|
orthogonal box = (-12.865 -12.865 -64.829) to (12.865 12.865 64.829)
|
||||||
|
1 by 1 by 8 MPI processor grid
|
||||||
|
reading atoms ...
|
||||||
|
1380 atoms
|
||||||
|
0 = max # of 1-2 neighbors
|
||||||
|
0 = max # of 1-3 neighbors
|
||||||
|
0 = max # of 1-4 neighbors
|
||||||
|
1 = max # of special neighbors
|
||||||
|
special bonds CPU = 0.000566959 secs
|
||||||
|
read_data CPU = 0.00661397 secs
|
||||||
|
velocity all create 3.0000e+02 16611 rot yes dist gaussian
|
||||||
|
|
||||||
|
# Assign potentials
|
||||||
|
pair_coeff 1 1 table benzene_water.pair.table PairBB
|
||||||
|
WARNING: 33 of 500 force values in table are inconsistent with -dE/dr.
|
||||||
|
Should only be flagged at inflection points (../pair_table.cpp:483)
|
||||||
|
WARNING: 150 of 500 distance values in table with relative error
|
||||||
|
over 1e-06 to re-computed values (../pair_table.cpp:492)
|
||||||
|
pair_coeff 1 2 table benzene_water.pair.table PairWW
|
||||||
|
WARNING: 61 of 500 force values in table are inconsistent with -dE/dr.
|
||||||
|
Should only be flagged at inflection points (../pair_table.cpp:483)
|
||||||
|
WARNING: 90 of 500 distance values in table with relative error
|
||||||
|
over 1e-06 to re-computed values (../pair_table.cpp:492)
|
||||||
|
pair_coeff 2 2 table benzene_water.pair.table PairBW
|
||||||
|
WARNING: 108 of 500 force values in table are inconsistent with -dE/dr.
|
||||||
|
Should only be flagged at inflection points (../pair_table.cpp:483)
|
||||||
|
WARNING: 135 of 500 distance values in table with relative error
|
||||||
|
over 1e-06 to re-computed values (../pair_table.cpp:492)
|
||||||
|
pair_coeff * * local/density benzene_water.localdensity.table
|
||||||
|
|
||||||
|
# Recentering during minimization and equilibration
|
||||||
|
fix recentering all recenter 0.0 0.0 0.0 units box
|
||||||
|
|
||||||
|
# Thermostat & time integration
|
||||||
|
timestep 2.0
|
||||||
|
thermo 100
|
||||||
|
thermo_style custom temp ke pe etotal ebond eangle edihed evdwl
|
||||||
|
|
||||||
|
# Minimization
|
||||||
|
minimize 1.e-4 0.0 10000 10000
|
||||||
|
WARNING: Using 'neigh_modify every 1 delay 0 check yes' setting during minimization (../min.cpp:168)
|
||||||
|
Neighbor list info ...
|
||||||
|
update every 1 steps, delay 0 steps, check yes
|
||||||
|
max neighbors/atom: 2000, page size: 100000
|
||||||
|
master list distance cutoff = 15.25
|
||||||
|
ghost atom cutoff = 15.25
|
||||||
|
binsize = 7.625, bins = 4 4 18
|
||||||
|
2 neighbor lists, perpetual/occasional/extra = 2 0 0
|
||||||
|
(1) pair table, perpetual
|
||||||
|
attributes: half, newton on
|
||||||
|
pair build: half/bin/newton
|
||||||
|
stencil: half/bin/3d/newton
|
||||||
|
bin: standard
|
||||||
|
(2) pair local/density, perpetual, copy from (1)
|
||||||
|
attributes: half, newton on
|
||||||
|
pair build: copy
|
||||||
|
stencil: none
|
||||||
|
bin: none
|
||||||
|
Per MPI rank memory allocation (min/avg/max) = 8.061 | 8.32 | 8.674 Mbytes
|
||||||
|
Temp KinEng PotEng TotEng E_bond E_angle E_dihed E_vdwl
|
||||||
|
300 1233.1611 4162.3053 5395.4665 0 0 0 4162.3053
|
||||||
|
300 1233.1611 2275.526 3508.6871 0 0 0 2275.526
|
||||||
|
Loop time of 0.352822 on 8 procs for 40 steps with 1380 atoms
|
||||||
|
|
||||||
|
71.3% CPU use with 8 MPI tasks x no OpenMP threads
|
||||||
|
|
||||||
|
Minimization stats:
|
||||||
|
Stopping criterion = linesearch alpha is zero
|
||||||
|
Energy initial, next-to-last, final =
|
||||||
|
4162.30533361 2208.86525108 2275.52597861
|
||||||
|
Force two-norm initial, final = 259.364 69.3915
|
||||||
|
Force max component initial, final = 22.2077 8.31436
|
||||||
|
Final line search alpha, max atom move = 2.90022e-12 2.41135e-11
|
||||||
|
Iterations, force evaluations = 40 110
|
||||||
|
|
||||||
|
MPI task timing breakdown:
|
||||||
|
Section | min time | avg time | max time |%varavg| %total
|
||||||
|
---------------------------------------------------------------
|
||||||
|
Pair | 0.053192 | 0.23903 | 0.32779 | 17.2 | 67.75
|
||||||
|
Bond | 9.0599e-06 | 1.6302e-05 | 2.5272e-05 | 0.0 | 0.00
|
||||||
|
Neigh | 0.00044513 | 0.0023614 | 0.0063851 | 5.1 | 0.67
|
||||||
|
Comm | 0.015469 | 0.090432 | 0.20295 | 20.0 | 25.63
|
||||||
|
Output | 0 | 0 | 0 | 0.0 | 0.00
|
||||||
|
Modify | 0 | 0 | 0 | 0.0 | 0.00
|
||||||
|
Other | | 0.02098 | | | 5.95
|
||||||
|
|
||||||
|
Nlocal: 172.5 ave 348 max 72 min
|
||||||
|
Histogram: 5 0 0 0 0 0 0 0 1 2
|
||||||
|
Nghost: 2193.62 ave 4352 max 932 min
|
||||||
|
Histogram: 3 0 0 2 0 0 2 0 0 1
|
||||||
|
Neighs: 9700.5 ave 20535 max 3685 min
|
||||||
|
Histogram: 5 0 0 0 0 0 0 1 0 2
|
||||||
|
|
||||||
|
Total # of neighbors = 77604
|
||||||
|
Ave neighs/atom = 56.2348
|
||||||
|
Ave special neighs/atom = 0
|
||||||
|
Neighbor list builds = 2
|
||||||
|
Dangerous builds = 0
|
||||||
|
|
||||||
|
# Set up integration parameters
|
||||||
|
fix timeintegration all nve
|
||||||
|
fix thermostat all langevin 3.0000e+02 3.0000e+02 1.0000e+02 81890
|
||||||
|
|
||||||
|
# Equilibration (for realistic results, run for 5000000 steps)
|
||||||
|
reset_timestep 0
|
||||||
|
run 5000
|
||||||
|
WARNING: Fix recenter should come after all other integration fixes (../fix_recenter.cpp:131)
|
||||||
|
Per MPI rank memory allocation (min/avg/max) = 6.936 | 7.195 | 7.552 Mbytes
|
||||||
|
Temp KinEng PotEng TotEng E_bond E_angle E_dihed E_vdwl
|
||||||
|
300 1233.1611 2866.9109 4100.0721 0 0 0 2866.9109
|
||||||
|
273.33541 1123.5553 3983.2007 5106.756 0 0 0 3983.2007
|
||||||
|
293.68078 1207.1857 3319.6601 4526.8458 0 0 0 3319.6601
|
||||||
|
314.21462 1291.5908 3389.2178 4680.8086 0 0 0 3389.2178
|
||||||
|
323.77563 1330.8917 3332.9828 4663.8745 0 0 0 3332.9828
|
||||||
|
302.5902 1243.8082 3461.7692 4705.5774 0 0 0 3461.7692
|
||||||
|
295.39324 1214.2249 3411.5727 4625.7976 0 0 0 3411.5727
|
||||||
|
320.52341 1317.5234 3453.1931 4770.7164 0 0 0 3453.1931
|
||||||
|
312.00777 1282.5195 3403.3443 4685.8638 0 0 0 3403.3443
|
||||||
|
307.96774 1265.9128 3429.7809 4695.6937 0 0 0 3429.7809
|
||||||
|
294.75922 1211.6187 3388.8404 4600.4591 0 0 0 3388.8404
|
||||||
|
311.24567 1279.3869 3514.9603 4794.3472 0 0 0 3514.9603
|
||||||
|
306.6152 1260.3531 3447.2011 4707.5542 0 0 0 3447.2011
|
||||||
|
305.23306 1254.6718 3375.5092 4630.181 0 0 0 3375.5092
|
||||||
|
321.62889 1322.0675 3460.2581 4782.3256 0 0 0 3460.2581
|
||||||
|
316.37725 1300.4804 3437.0312 4737.5116 0 0 0 3437.0312
|
||||||
|
322.90522 1327.3139 3389.1262 4716.44 0 0 0 3389.1262
|
||||||
|
307.57893 1264.3146 3359.8491 4624.1637 0 0 0 3359.8491
|
||||||
|
302.22607 1242.3115 3406.1711 4648.4826 0 0 0 3406.1711
|
||||||
|
302.73997 1244.4239 3220.2582 4464.6821 0 0 0 3220.2582
|
||||||
|
303.66194 1248.2137 3318.4629 4566.6765 0 0 0 3318.4629
|
||||||
|
308.73862 1269.0815 3369.5894 4638.671 0 0 0 3369.5894
|
||||||
|
315.60294 1297.2976 3411.2405 4708.5381 0 0 0 3411.2405
|
||||||
|
310.0113 1274.3129 3360.1054 4634.4183 0 0 0 3360.1054
|
||||||
|
302.36229 1242.8714 3326.9845 4569.8559 0 0 0 3326.9845
|
||||||
|
317.78659 1306.2735 3355.4976 4661.7711 0 0 0 3355.4976
|
||||||
|
302.50479 1243.4571 3317.6846 4561.1417 0 0 0 3317.6846
|
||||||
|
304.29249 1250.8056 3423.5068 4674.3124 0 0 0 3423.5068
|
||||||
|
305.99948 1257.8222 3432.9395 4690.7617 0 0 0 3432.9395
|
||||||
|
309.93363 1273.9937 3393.657 4667.6506 0 0 0 3393.657
|
||||||
|
316.14884 1299.5415 3463.0636 4762.6051 0 0 0 3463.0636
|
||||||
|
300.38817 1234.7567 3309.2495 4544.0062 0 0 0 3309.2495
|
||||||
|
311.05735 1278.6128 3304.4418 4583.0546 0 0 0 3304.4418
|
||||||
|
311.11872 1278.865 3291.1891 4570.0542 0 0 0 3291.1891
|
||||||
|
315.74338 1297.8749 3341.3063 4639.1812 0 0 0 3341.3063
|
||||||
|
297.5658 1223.1552 3316.3862 4539.5414 0 0 0 3316.3862
|
||||||
|
311.79033 1281.6257 3357.4556 4639.0813 0 0 0 3357.4556
|
||||||
|
310.93666 1278.1167 3414.7694 4692.8861 0 0 0 3414.7694
|
||||||
|
307.37298 1263.468 3337.3889 4600.8569 0 0 0 3337.3889
|
||||||
|
298.84185 1228.4005 3329.6173 4558.0178 0 0 0 3329.6173
|
||||||
|
310.54684 1276.5143 3351.0852 4627.5995 0 0 0 3351.0852
|
||||||
|
300.0871 1233.5191 3302.2315 4535.7506 0 0 0 3302.2315
|
||||||
|
304.69078 1252.4427 3324.2508 4576.6935 0 0 0 3324.2508
|
||||||
|
313.50714 1288.6827 3330.4088 4619.0915 0 0 0 3330.4088
|
||||||
|
329.80018 1355.6559 3301.86 4657.5159 0 0 0 3301.86
|
||||||
|
304.57609 1251.9713 3365.2938 4617.2652 0 0 0 3365.2938
|
||||||
|
308.73584 1269.0701 3344.4155 4613.4856 0 0 0 3344.4155
|
||||||
|
306.90951 1261.5629 3304.4698 4566.0327 0 0 0 3304.4698
|
||||||
|
308.85761 1269.5707 3392.1511 4661.7218 0 0 0 3392.1511
|
||||||
|
302.78788 1244.6208 3317.0849 4561.7057 0 0 0 3317.0849
|
||||||
|
321.68092 1322.2813 3321.5755 4643.8568 0 0 0 3321.5755
|
||||||
|
Loop time of 16.3061 on 8 procs for 5000 steps with 1380 atoms
|
||||||
|
|
||||||
|
Performance: 52.986 ns/day, 0.453 hours/ns, 306.634 timesteps/s
|
||||||
|
69.6% CPU use with 8 MPI tasks x no OpenMP threads
|
||||||
|
|
||||||
|
MPI task timing breakdown:
|
||||||
|
Section | min time | avg time | max time |%varavg| %total
|
||||||
|
---------------------------------------------------------------
|
||||||
|
Pair | 2.1872 | 10.542 | 14.607 | 116.7 | 64.65
|
||||||
|
Bond | 0.00044084 | 0.00069669 | 0.00095081 | 0.0 | 0.00
|
||||||
|
Neigh | 0.026948 | 0.15225 | 0.44344 | 42.0 | 0.93
|
||||||
|
Comm | 0.63452 | 4.2953 | 9.49 | 133.9 | 26.34
|
||||||
|
Output | 0.0016391 | 0.012378 | 0.050919 | 13.9 | 0.08
|
||||||
|
Modify | 0.45894 | 1.2107 | 4.4629 | 116.4 | 7.42
|
||||||
|
Other | | 0.09292 | | | 0.57
|
||||||
|
|
||||||
|
Nlocal: 172.5 ave 380 max 70 min
|
||||||
|
Histogram: 5 0 0 0 0 0 0 1 1 1
|
||||||
|
Nghost: 2213 ave 4440 max 903 min
|
||||||
|
Histogram: 3 0 0 2 0 0 2 0 0 1
|
||||||
|
Neighs: 10042.5 ave 24051 max 3500 min
|
||||||
|
Histogram: 5 0 0 0 0 0 0 1 1 1
|
||||||
|
|
||||||
|
Total # of neighbors = 80340
|
||||||
|
Ave neighs/atom = 58.2174
|
||||||
|
Ave special neighs/atom = 0
|
||||||
|
Neighbor list builds = 123
|
||||||
|
Dangerous builds = 1
|
||||||
|
|
||||||
|
# Turn off recentering during production phase
|
||||||
|
unfix recentering
|
||||||
|
|
||||||
|
# Setup trajectory output
|
||||||
|
dump myDump all custom 100 benzene_water.lammpstrj.gz id type x y z element
|
||||||
|
dump_modify myDump element B W
|
||||||
|
dump_modify myDump sort id
|
||||||
|
|
||||||
|
# Production (for realistic results, run for 10000000 steps)
|
||||||
|
reset_timestep 0
|
||||||
|
run 1000
|
||||||
|
Per MPI rank memory allocation (min/avg/max) = 8.232 | 8.492 | 8.851 Mbytes
|
||||||
|
Temp KinEng PotEng TotEng E_bond E_angle E_dihed E_vdwl
|
||||||
|
321.68092 1322.2813 3784.0834 5106.3647 0 0 0 3784.0834
|
||||||
|
310.59763 1276.7231 3318.3283 4595.0513 0 0 0 3318.3283
|
||||||
|
303.39445 1247.1141 3324.1191 4571.2332 0 0 0 3324.1191
|
||||||
|
311.37275 1279.9092 3305.0901 4584.9993 0 0 0 3305.0901
|
||||||
|
311.29071 1279.572 3248.216 4527.788 0 0 0 3248.216
|
||||||
|
314.53456 1292.906 3283.4563 4576.3623 0 0 0 3283.4563
|
||||||
|
316.52595 1301.0916 3258.9171 4560.0087 0 0 0 3258.9171
|
||||||
|
318.92447 1310.9509 3235.6256 4546.5765 0 0 0 3235.6256
|
||||||
|
311.79212 1281.6331 3308.099 4589.7321 0 0 0 3308.099
|
||||||
|
305.52477 1255.8709 3267.6907 4523.5616 0 0 0 3267.6907
|
||||||
|
301.07457 1237.5782 3206.3997 4443.9779 0 0 0 3206.3997
|
||||||
|
Loop time of 4.44139 on 8 procs for 1000 steps with 1380 atoms
|
||||||
|
|
||||||
|
Performance: 38.907 ns/day, 0.617 hours/ns, 225.155 timesteps/s
|
||||||
|
60.8% CPU use with 8 MPI tasks x no OpenMP threads
|
||||||
|
|
||||||
|
MPI task timing breakdown:
|
||||||
|
Section | min time | avg time | max time |%varavg| %total
|
||||||
|
---------------------------------------------------------------
|
||||||
|
Pair | 0.656 | 2.5078 | 3.5775 | 57.7 | 56.46
|
||||||
|
Bond | 0.00013375 | 0.0001854 | 0.0002377 | 0.0 | 0.00
|
||||||
|
Neigh | 0.0048757 | 0.029188 | 0.090432 | 18.9 | 0.66
|
||||||
|
Comm | 0.51836 | 1.4427 | 2.6285 | 56.9 | 32.48
|
||||||
|
Output | 0.083084 | 0.089199 | 0.10333 | 2.3 | 2.01
|
||||||
|
Modify | 0.0087376 | 0.019705 | 0.038437 | 8.4 | 0.44
|
||||||
|
Other | | 0.3526 | | | 7.94
|
||||||
|
|
||||||
|
Nlocal: 172.5 ave 388 max 69 min
|
||||||
|
Histogram: 5 0 0 0 0 0 0 2 0 1
|
||||||
|
Nghost: 2207.88 ave 4429 max 896 min
|
||||||
|
Histogram: 3 0 0 2 0 0 2 0 0 1
|
||||||
|
Neighs: 10094.1 ave 24847 max 3403 min
|
||||||
|
Histogram: 5 0 0 0 0 0 1 1 0 1
|
||||||
|
|
||||||
|
Total # of neighbors = 80753
|
||||||
|
Ave neighs/atom = 58.5167
|
||||||
|
Ave special neighs/atom = 0
|
||||||
|
Neighbor list builds = 23
|
||||||
|
Dangerous builds = 0
|
||||||
|
|
||||||
|
|
||||||
|
Total wall time: 0:00:21
|
||||||
@ -0,0 +1,226 @@
|
|||||||
|
LAMMPS (7 Aug 2019)
|
||||||
|
# LAMMPS input file for 50.0% methanol mole fraction solution
|
||||||
|
# with 2500 methanol molecules in implicit water.
|
||||||
|
#
|
||||||
|
#
|
||||||
|
# Author: David Rosenberger, van der Vegt Group, TU Darmstadt
|
||||||
|
#
|
||||||
|
# Refer: Rosenberger, Sanyal, Shell, van der Vegt, J. Chem. Theory Comput. 15, 2881-2895 (2019)
|
||||||
|
|
||||||
|
|
||||||
|
# Initialize simulation box
|
||||||
|
dimension 3
|
||||||
|
boundary p p p
|
||||||
|
units real
|
||||||
|
atom_style molecular
|
||||||
|
|
||||||
|
# Set potential styles
|
||||||
|
pair_style hybrid/overlay table spline 500 local/density
|
||||||
|
|
||||||
|
# Read molecule data and set initial velocities
|
||||||
|
read_data methanol_implicit_water.data
|
||||||
|
orthogonal box = (-31.123 -31.123 -31.123) to (31.123 31.123 31.123)
|
||||||
|
2 by 2 by 2 MPI processor grid
|
||||||
|
reading atoms ...
|
||||||
|
2500 atoms
|
||||||
|
0 = max # of 1-2 neighbors
|
||||||
|
0 = max # of 1-3 neighbors
|
||||||
|
0 = max # of 1-4 neighbors
|
||||||
|
1 = max # of special neighbors
|
||||||
|
special bonds CPU = 0.00063014 secs
|
||||||
|
read_data CPU = 0.00599909 secs
|
||||||
|
velocity all create 3.0000e+02 12142 rot yes dist gaussian
|
||||||
|
|
||||||
|
# Assign potentials
|
||||||
|
pair_coeff 1 1 table methanol_implicit_water.pair.table PairMM
|
||||||
|
WARNING: 93 of 500 force values in table are inconsistent with -dE/dr.
|
||||||
|
Should only be flagged at inflection points (../pair_table.cpp:483)
|
||||||
|
WARNING: 254 of 500 distance values in table with relative error
|
||||||
|
over 1e-06 to re-computed values (../pair_table.cpp:492)
|
||||||
|
pair_coeff * * local/density methanol_implicit_water.localdensity.table
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
#Recentering during minimization and equilibration
|
||||||
|
fix recentering all recenter 0.0 0.0 0.0 units box
|
||||||
|
|
||||||
|
#Thermostat & time integration
|
||||||
|
timestep 1.0
|
||||||
|
thermo 100
|
||||||
|
thermo_style custom etotal ke pe temp evdwl
|
||||||
|
|
||||||
|
#minimization
|
||||||
|
minimize 1.e-4 0.0 1000 1000
|
||||||
|
WARNING: Using 'neigh_modify every 1 delay 0 check yes' setting during minimization (../min.cpp:168)
|
||||||
|
Neighbor list info ...
|
||||||
|
update every 1 steps, delay 0 steps, check yes
|
||||||
|
max neighbors/atom: 2000, page size: 100000
|
||||||
|
master list distance cutoff = 17
|
||||||
|
ghost atom cutoff = 17
|
||||||
|
binsize = 8.5, bins = 8 8 8
|
||||||
|
2 neighbor lists, perpetual/occasional/extra = 2 0 0
|
||||||
|
(1) pair table, perpetual
|
||||||
|
attributes: half, newton on
|
||||||
|
pair build: half/bin/newton
|
||||||
|
stencil: half/bin/3d/newton
|
||||||
|
bin: standard
|
||||||
|
(2) pair local/density, perpetual, copy from (1)
|
||||||
|
attributes: half, newton on
|
||||||
|
pair build: copy
|
||||||
|
stencil: none
|
||||||
|
bin: none
|
||||||
|
Per MPI rank memory allocation (min/avg/max) = 7.411 | 7.411 | 7.412 Mbytes
|
||||||
|
TotEng KinEng PotEng Temp E_vdwl
|
||||||
|
1470.3564 2234.7133 -764.35689 300 -764.35689
|
||||||
|
46.496766 2234.7133 -2188.2165 300 -2188.2165
|
||||||
|
7.9030246 2234.7133 -2226.8103 300 -2226.8103
|
||||||
|
Loop time of 0.463996 on 8 procs for 121 steps with 2500 atoms
|
||||||
|
|
||||||
|
91.4% CPU use with 8 MPI tasks x no OpenMP threads
|
||||||
|
|
||||||
|
Minimization stats:
|
||||||
|
Stopping criterion = linesearch alpha is zero
|
||||||
|
Energy initial, next-to-last, final =
|
||||||
|
-764.356892369 -2227.85589084 -2226.81026984
|
||||||
|
Force two-norm initial, final = 134.911 3.83896
|
||||||
|
Force max component initial, final = 14.1117 1.07422
|
||||||
|
Final line search alpha, max atom move = 5.06747e-10 5.44356e-10
|
||||||
|
Iterations, force evaluations = 121 154
|
||||||
|
|
||||||
|
MPI task timing breakdown:
|
||||||
|
Section | min time | avg time | max time |%varavg| %total
|
||||||
|
---------------------------------------------------------------
|
||||||
|
Pair | 0.41442 | 0.41976 | 0.42434 | 0.5 | 90.47
|
||||||
|
Bond | 1.1683e-05 | 2.0713e-05 | 3.5048e-05 | 0.0 | 0.00
|
||||||
|
Neigh | 0.0084722 | 0.0090862 | 0.010038 | 0.5 | 1.96
|
||||||
|
Comm | 0.022712 | 0.028157 | 0.034072 | 1.9 | 6.07
|
||||||
|
Output | 3.1948e-05 | 3.6925e-05 | 6.6996e-05 | 0.0 | 0.01
|
||||||
|
Modify | 0 | 0 | 0 | 0.0 | 0.00
|
||||||
|
Other | | 0.006937 | | | 1.50
|
||||||
|
|
||||||
|
Nlocal: 312.5 ave 333 max 299 min
|
||||||
|
Histogram: 2 2 0 0 1 0 2 0 0 1
|
||||||
|
Nghost: 2546 ave 2580 max 2517 min
|
||||||
|
Histogram: 1 1 0 3 0 1 0 0 0 2
|
||||||
|
Neighs: 33215.4 ave 37251 max 29183 min
|
||||||
|
Histogram: 1 0 0 1 2 2 0 1 0 1
|
||||||
|
|
||||||
|
Total # of neighbors = 265723
|
||||||
|
Ave neighs/atom = 106.289
|
||||||
|
Ave special neighs/atom = 0
|
||||||
|
Neighbor list builds = 6
|
||||||
|
Dangerous builds = 0
|
||||||
|
|
||||||
|
#set up integration parameters
|
||||||
|
fix timeintegration all nve
|
||||||
|
fix thermostat all langevin 3.0000e+02 3.0000e+02 1.0000e+02 59915
|
||||||
|
|
||||||
|
#Equilibration (for realistic results, run for 2000000 steps)
|
||||||
|
reset_timestep 0
|
||||||
|
thermo 200
|
||||||
|
thermo_style custom etotal ke pe temp evdwl
|
||||||
|
|
||||||
|
#run equilibration
|
||||||
|
run 2000
|
||||||
|
WARNING: Fix recenter should come after all other integration fixes (../fix_recenter.cpp:131)
|
||||||
|
Per MPI rank memory allocation (min/avg/max) = 6.286 | 6.286 | 6.287 Mbytes
|
||||||
|
TotEng KinEng PotEng Temp E_vdwl
|
||||||
|
177.26822 2234.7133 -2057.4451 300 -2057.4451
|
||||||
|
736.24287 2151.2608 -1415.0179 288.79688 -1415.0179
|
||||||
|
963.07617 2090.6433 -1127.5671 280.65926 -1127.5671
|
||||||
|
1148.9049 2173.1327 -1024.2279 291.73309 -1024.2279
|
||||||
|
1303.6409 2279.8586 -976.21767 306.06055 -976.21767
|
||||||
|
1355.42 2281.0383 -925.61826 306.21892 -925.61826
|
||||||
|
1394.5206 2276.2093 -881.68863 305.57064 -881.68863
|
||||||
|
1346.9764 2215.2973 -868.32091 297.3935 -868.32091
|
||||||
|
1381.3654 2248.8061 -867.44063 301.89189 -867.44063
|
||||||
|
1315.8059 2189.3193 -873.51332 293.90606 -873.51332
|
||||||
|
1314.4456 2209.7431 -895.29752 296.64787 -895.29752
|
||||||
|
Loop time of 6.38989 on 8 procs for 2000 steps with 2500 atoms
|
||||||
|
|
||||||
|
Performance: 27.043 ns/day, 0.887 hours/ns, 312.994 timesteps/s
|
||||||
|
80.5% CPU use with 8 MPI tasks x no OpenMP threads
|
||||||
|
|
||||||
|
MPI task timing breakdown:
|
||||||
|
Section | min time | avg time | max time |%varavg| %total
|
||||||
|
---------------------------------------------------------------
|
||||||
|
Pair | 5.2693 | 5.3572 | 5.457 | 2.1 | 83.84
|
||||||
|
Bond | 0.00028825 | 0.00033835 | 0.00039148 | 0.0 | 0.01
|
||||||
|
Neigh | 0.0296 | 0.032337 | 0.035071 | 0.9 | 0.51
|
||||||
|
Comm | 0.64679 | 0.73397 | 0.80847 | 5.2 | 11.49
|
||||||
|
Output | 0.00033498 | 0.00051582 | 0.0015228 | 0.0 | 0.01
|
||||||
|
Modify | 0.16395 | 0.18919 | 0.21056 | 3.9 | 2.96
|
||||||
|
Other | | 0.07636 | | | 1.19
|
||||||
|
|
||||||
|
Nlocal: 312.5 ave 337 max 295 min
|
||||||
|
Histogram: 2 2 0 1 0 0 0 1 1 1
|
||||||
|
Nghost: 2551.62 ave 2582 max 2525 min
|
||||||
|
Histogram: 2 1 0 0 1 1 1 0 1 1
|
||||||
|
Neighs: 33241.8 ave 37659 max 29705 min
|
||||||
|
Histogram: 2 0 0 2 2 0 0 0 1 1
|
||||||
|
|
||||||
|
Total # of neighbors = 265934
|
||||||
|
Ave neighs/atom = 106.374
|
||||||
|
Ave special neighs/atom = 0
|
||||||
|
Neighbor list builds = 21
|
||||||
|
Dangerous builds = 0
|
||||||
|
|
||||||
|
#turn off recentering during production run
|
||||||
|
unfix recentering
|
||||||
|
|
||||||
|
|
||||||
|
#setup trajectory output
|
||||||
|
dump myDump all custom 100 methanol_implicit_water.lammpstrj.gz id type x y z element
|
||||||
|
dump_modify myDump element M
|
||||||
|
dump_modify myDump sort id
|
||||||
|
|
||||||
|
#run production (for realistic results, run for 10000000 steps)
|
||||||
|
reset_timestep 0
|
||||||
|
thermo 1000
|
||||||
|
thermo_style custom etotal ke pe temp evdwl
|
||||||
|
run 10000
|
||||||
|
Per MPI rank memory allocation (min/avg/max) = 7.588 | 7.589 | 7.589 Mbytes
|
||||||
|
TotEng KinEng PotEng Temp E_vdwl
|
||||||
|
1442.5428 2209.7431 -767.20027 296.64787 -767.20027
|
||||||
|
1391.8624 2262.6889 -870.82656 303.7556 -870.82656
|
||||||
|
1375.914 2244.6176 -868.7036 301.3296 -868.7036
|
||||||
|
1345.9064 2227.2324 -881.32599 298.99573 -881.32599
|
||||||
|
1379.2334 2278.1156 -898.88222 305.82657 -898.88222
|
||||||
|
1389.7928 2255.8062 -866.01341 302.83163 -866.01341
|
||||||
|
1380.4549 2258.2108 -877.75582 303.15443 -877.75582
|
||||||
|
1380.8489 2256.9432 -876.09428 302.98426 -876.09428
|
||||||
|
1326.5151 2225.7408 -899.22577 298.79549 -899.22577
|
||||||
|
1376.6025 2253.0128 -876.41028 302.45662 -876.41028
|
||||||
|
1331.0008 2218.1033 -887.10258 297.77019 -887.10258
|
||||||
|
Loop time of 25.4591 on 8 procs for 10000 steps with 2500 atoms
|
||||||
|
|
||||||
|
Performance: 33.937 ns/day, 0.707 hours/ns, 392.787 timesteps/s
|
||||||
|
89.3% CPU use with 8 MPI tasks x no OpenMP threads
|
||||||
|
|
||||||
|
MPI task timing breakdown:
|
||||||
|
Section | min time | avg time | max time |%varavg| %total
|
||||||
|
---------------------------------------------------------------
|
||||||
|
Pair | 21.635 | 21.916 | 22.237 | 3.9 | 86.08
|
||||||
|
Bond | 0.0011308 | 0.0013149 | 0.0016932 | 0.5 | 0.01
|
||||||
|
Neigh | 0.14593 | 0.15675 | 0.16667 | 1.9 | 0.62
|
||||||
|
Comm | 1.3789 | 1.7502 | 1.9558 | 13.7 | 6.87
|
||||||
|
Output | 0.34664 | 0.82927 | 1.2013 | 32.8 | 3.26
|
||||||
|
Modify | 0.24904 | 0.25842 | 0.26907 | 1.2 | 1.02
|
||||||
|
Other | | 0.5475 | | | 2.15
|
||||||
|
|
||||||
|
Nlocal: 312.5 ave 327 max 298 min
|
||||||
|
Histogram: 2 0 0 1 1 0 1 1 1 1
|
||||||
|
Nghost: 2575 ave 2601 max 2559 min
|
||||||
|
Histogram: 2 0 3 1 0 0 0 0 1 1
|
||||||
|
Neighs: 33223.2 ave 35920 max 30303 min
|
||||||
|
Histogram: 1 1 1 1 0 1 0 0 0 3
|
||||||
|
|
||||||
|
Total # of neighbors = 265786
|
||||||
|
Ave neighs/atom = 106.314
|
||||||
|
Ave special neighs/atom = 0
|
||||||
|
Neighbor list builds = 103
|
||||||
|
Dangerous builds = 0
|
||||||
|
|
||||||
|
|
||||||
|
Total wall time: 0:00:32
|
||||||
File diff suppressed because it is too large
Load Diff
@ -0,0 +1,68 @@
|
|||||||
|
# LAMMPS input file for 50.0% methanol mole fraction solution
|
||||||
|
# with 2500 methanol molecules in implicit water.
|
||||||
|
#
|
||||||
|
#
|
||||||
|
# Author: David Rosenberger, van der Vegt Group, TU Darmstadt
|
||||||
|
#
|
||||||
|
# Refer: Rosenberger, Sanyal, Shell, van der Vegt, J. Chem. Theory Comput. 15, 2881-2895 (2019)
|
||||||
|
|
||||||
|
|
||||||
|
# Initialize simulation box
|
||||||
|
dimension 3
|
||||||
|
boundary p p p
|
||||||
|
units real
|
||||||
|
atom_style molecular
|
||||||
|
|
||||||
|
# Set potential styles
|
||||||
|
pair_style hybrid/overlay table spline 500 local/density
|
||||||
|
|
||||||
|
# Read molecule data and set initial velocities
|
||||||
|
read_data methanol_implicit_water.data
|
||||||
|
velocity all create 3.0000e+02 12142 rot yes dist gaussian
|
||||||
|
|
||||||
|
# Assign potentials
|
||||||
|
pair_coeff 1 1 table methanol_implicit_water.pair.table PairMM
|
||||||
|
pair_coeff * * local/density methanol_implicit_water.localdensity.table
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
#Recentering during minimization and equilibration
|
||||||
|
fix recentering all recenter 0.0 0.0 0.0 units box
|
||||||
|
|
||||||
|
#Thermostat & time integration
|
||||||
|
timestep 1.0
|
||||||
|
thermo 100
|
||||||
|
thermo_style custom etotal ke pe temp evdwl
|
||||||
|
|
||||||
|
#minimization
|
||||||
|
minimize 1.e-4 0.0 1000 1000
|
||||||
|
|
||||||
|
#set up integration parameters
|
||||||
|
fix timeintegration all nve
|
||||||
|
fix thermostat all langevin 3.0000e+02 3.0000e+02 1.0000e+02 59915
|
||||||
|
|
||||||
|
#Equilibration (for realistic results, run for 2000000 steps)
|
||||||
|
reset_timestep 0
|
||||||
|
thermo 200
|
||||||
|
thermo_style custom etotal ke pe temp evdwl
|
||||||
|
|
||||||
|
#run equilibration
|
||||||
|
run 2000
|
||||||
|
|
||||||
|
#turn off recentering during production run
|
||||||
|
unfix recentering
|
||||||
|
|
||||||
|
|
||||||
|
#setup trajectory output
|
||||||
|
dump myDump all custom 100 methanol_implicit_water.lammpstrj.gz id type x y z element
|
||||||
|
dump_modify myDump element M
|
||||||
|
dump_modify myDump sort id
|
||||||
|
|
||||||
|
#run production (for realistic results, run for 10000000 steps)
|
||||||
|
reset_timestep 0
|
||||||
|
thermo 1000
|
||||||
|
thermo_style custom etotal ke pe temp evdwl
|
||||||
|
run 10000
|
||||||
|
|
||||||
|
|
||||||
@ -0,0 +1,509 @@
|
|||||||
|
#LOCAL DENSITY POTENTIALS
|
||||||
|
|
||||||
|
1 500
|
||||||
|
|
||||||
|
5.3000000e+00 6.3000000e+00
|
||||||
|
1
|
||||||
|
1
|
||||||
|
0.0000000e+00 2.6000000e+01 5.2104208e-02
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4810000e-01
|
||||||
|
1.4807157e-01
|
||||||
|
1.4782582e-01
|
||||||
|
1.4711763e-01
|
||||||
|
1.4570179e-01
|
||||||
|
1.4333312e-01
|
||||||
|
1.3976643e-01
|
||||||
|
1.3478059e-01
|
||||||
|
1.2856173e-01
|
||||||
|
1.2163552e-01
|
||||||
|
1.1453802e-01
|
||||||
|
1.0780525e-01
|
||||||
|
1.0197328e-01
|
||||||
|
9.7575837e-02
|
||||||
|
9.4875548e-02
|
||||||
|
9.3613063e-02
|
||||||
|
9.3469690e-02
|
||||||
|
9.4126738e-02
|
||||||
|
9.5265515e-02
|
||||||
|
9.6567329e-02
|
||||||
|
9.7735007e-02
|
||||||
|
9.8575495e-02
|
||||||
|
9.8927186e-02
|
||||||
|
9.8628481e-02
|
||||||
|
9.7517779e-02
|
||||||
|
9.5433481e-02
|
||||||
|
9.2235018e-02
|
||||||
|
8.8072568e-02
|
||||||
|
8.3308496e-02
|
||||||
|
7.8309990e-02
|
||||||
|
7.3444241e-02
|
||||||
|
6.9078438e-02
|
||||||
|
6.5577180e-02
|
||||||
|
6.3110699e-02
|
||||||
|
6.1523109e-02
|
||||||
|
6.0627357e-02
|
||||||
|
6.0236386e-02
|
||||||
|
6.0163144e-02
|
||||||
|
6.0220573e-02
|
||||||
|
6.0233006e-02
|
||||||
|
6.0072080e-02
|
||||||
|
5.9621717e-02
|
||||||
|
5.8765838e-02
|
||||||
|
5.7388366e-02
|
||||||
|
5.5373224e-02
|
||||||
|
5.2623498e-02
|
||||||
|
4.9261717e-02
|
||||||
|
4.5550390e-02
|
||||||
|
4.1754290e-02
|
||||||
|
3.8138193e-02
|
||||||
|
3.4966871e-02
|
||||||
|
3.2501662e-02
|
||||||
|
3.0825931e-02
|
||||||
|
2.9762256e-02
|
||||||
|
2.9112455e-02
|
||||||
|
2.8678347e-02
|
||||||
|
2.8261751e-02
|
||||||
|
2.7664487e-02
|
||||||
|
2.6737788e-02
|
||||||
|
2.5509284e-02
|
||||||
|
2.4045951e-02
|
||||||
|
2.2414767e-02
|
||||||
|
2.0682707e-02
|
||||||
|
1.8916748e-02
|
||||||
|
1.7179645e-02
|
||||||
|
1.5493687e-02
|
||||||
|
1.3858641e-02
|
||||||
|
1.2274032e-02
|
||||||
|
1.0739385e-02
|
||||||
|
9.2542252e-03
|
||||||
|
7.8179601e-03
|
||||||
|
6.4255437e-03
|
||||||
|
5.0662231e-03
|
||||||
|
3.7288715e-03
|
||||||
|
2.4023618e-03
|
||||||
|
1.0755673e-03
|
||||||
|
-2.6263394e-04
|
||||||
|
-1.6141074e-03
|
||||||
|
-2.9522803e-03
|
||||||
|
-4.2451362e-03
|
||||||
|
-5.4606586e-03
|
||||||
|
-6.5668312e-03
|
||||||
|
-7.5316377e-03
|
||||||
|
-8.3294239e-03
|
||||||
|
-8.9860017e-03
|
||||||
|
-9.5521117e-03
|
||||||
|
-1.0078658e-02
|
||||||
|
-1.0616544e-02
|
||||||
|
-1.1216675e-02
|
||||||
|
-1.1929199e-02
|
||||||
|
-1.2782684e-02
|
||||||
|
-1.3781467e-02
|
||||||
|
-1.4928602e-02
|
||||||
|
-1.6227139e-02
|
||||||
|
-1.7680132e-02
|
||||||
|
-1.9290577e-02
|
||||||
|
-2.1031059e-02
|
||||||
|
-2.2793537e-02
|
||||||
|
-2.4456753e-02
|
||||||
|
-2.5899451e-02
|
||||||
|
-2.7000374e-02
|
||||||
|
-2.7638267e-02
|
||||||
|
-2.7719868e-02
|
||||||
|
-2.7344330e-02
|
||||||
|
-2.6691680e-02
|
||||||
|
-2.5942229e-02
|
||||||
|
-2.5276286e-02
|
||||||
|
-2.4874159e-02
|
||||||
|
-2.4909370e-02
|
||||||
|
-2.5403835e-02
|
||||||
|
-2.6230347e-02
|
||||||
|
-2.7255409e-02
|
||||||
|
-2.8345523e-02
|
||||||
|
-2.9367192e-02
|
||||||
|
-3.0187085e-02
|
||||||
|
-3.0712590e-02
|
||||||
|
-3.0944560e-02
|
||||||
|
-3.0896910e-02
|
||||||
|
-3.0583557e-02
|
||||||
|
-3.0018416e-02
|
||||||
|
-2.9215405e-02
|
||||||
|
-2.8195478e-02
|
||||||
|
-2.7020910e-02
|
||||||
|
-2.5768997e-02
|
||||||
|
-2.4517057e-02
|
||||||
|
-2.3342408e-02
|
||||||
|
-2.2322368e-02
|
||||||
|
-2.1532406e-02
|
||||||
|
-2.1015034e-02
|
||||||
|
-2.0784355e-02
|
||||||
|
-2.0853543e-02
|
||||||
|
-2.1235771e-02
|
||||||
|
-2.1944214e-02
|
||||||
|
-2.2991215e-02
|
||||||
|
-2.4278363e-02
|
||||||
|
-2.5486458e-02
|
||||||
|
-2.6270119e-02
|
||||||
|
-2.6283964e-02
|
||||||
|
-2.5182614e-02
|
||||||
|
-2.2620686e-02
|
||||||
|
-1.8367122e-02
|
||||||
|
-1.2765600e-02
|
||||||
|
-6.3400224e-03
|
||||||
|
3.8564733e-04
|
||||||
|
6.8874449e-03
|
||||||
|
1.2641406e-02
|
||||||
|
1.7151899e-02
|
||||||
|
2.0334733e-02
|
||||||
|
2.2416173e-02
|
||||||
|
2.3630118e-02
|
||||||
|
2.4210466e-02
|
||||||
|
2.4391115e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
2.4405000e-02
|
||||||
|
|
||||||
File diff suppressed because it is too large
Load Diff
@ -84,6 +84,7 @@ pair_style lebedeva/z, Zbigniew Koziol (National Center for Nuclear Research), s
|
|||||||
pair_style lennard/mdf, Paolo Raiteri, p.raiteri at curtin.edu.au, 2 Dec 15
|
pair_style lennard/mdf, Paolo Raiteri, p.raiteri at curtin.edu.au, 2 Dec 15
|
||||||
pair_style list, Axel Kohlmeyer (Temple U), akohlmey at gmail.com, 1 Jun 13
|
pair_style list, Axel Kohlmeyer (Temple U), akohlmey at gmail.com, 1 Jun 13
|
||||||
pair_style lj/mdf, Paolo Raiteri, p.raiteri at curtin.edu.au, 2 Dec 15
|
pair_style lj/mdf, Paolo Raiteri, p.raiteri at curtin.edu.au, 2 Dec 15
|
||||||
|
pair_style local/density, Tanmoy Sanyal (tanmoy dot 7989 at gmail.com) and M. Scott Shell (UCSB), and David Rosenberger (TU Darmstadt), 9 Sept 19
|
||||||
pair_style kolmogorov/crespi/full, Wengen Ouyang (Tel Aviv University), w.g.ouyang at gmail dot com, 30 Mar 18
|
pair_style kolmogorov/crespi/full, Wengen Ouyang (Tel Aviv University), w.g.ouyang at gmail dot com, 30 Mar 18
|
||||||
pair_style kolmogorov/crespi/z, Jaap Kroes (Radboud U), jaapkroes at gmail dot com, 28 Feb 17
|
pair_style kolmogorov/crespi/z, Jaap Kroes (Radboud U), jaapkroes at gmail dot com, 28 Feb 17
|
||||||
pair_style meam/spline, Alexander Stukowski (LLNL), alex at stukowski.com, 1 Feb 12
|
pair_style meam/spline, Alexander Stukowski (LLNL), alex at stukowski.com, 1 Feb 12
|
||||||
|
|||||||
871
src/USER-MISC/pair_local_density.cpp
Normal file
871
src/USER-MISC/pair_local_density.cpp
Normal file
@ -0,0 +1,871 @@
|
|||||||
|
/* -*- c++ -*- ----------------------------------------------------------
|
||||||
|
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:
|
||||||
|
Tanmoy Sanyal, M.Scott Shell, UC Santa Barbara
|
||||||
|
David Rosenberger, TU Darmstadt
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
#include "pair_local_density.h"
|
||||||
|
#include <mpi.h>
|
||||||
|
#include <cmath>
|
||||||
|
#include <cstdlib>
|
||||||
|
#include <cstring>
|
||||||
|
#include <string>
|
||||||
|
#include "atom.h"
|
||||||
|
#include "force.h"
|
||||||
|
#include "comm.h"
|
||||||
|
#include "neighbor.h"
|
||||||
|
#include "neigh_list.h"
|
||||||
|
#include "neigh_request.h"
|
||||||
|
#include "memory.h"
|
||||||
|
#include "error.h"
|
||||||
|
#include "domain.h"
|
||||||
|
|
||||||
|
using namespace LAMMPS_NS;
|
||||||
|
|
||||||
|
#define MAXLINE 1024
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
PairLocalDensity::PairLocalDensity(LAMMPS *lmp) : Pair(lmp)
|
||||||
|
{
|
||||||
|
restartinfo = 0;
|
||||||
|
one_coeff = 1;
|
||||||
|
single_enable = 1;
|
||||||
|
|
||||||
|
// stuff read from tabulated file
|
||||||
|
nLD = 0;
|
||||||
|
nrho = 0;
|
||||||
|
rho_min = NULL;
|
||||||
|
rho_max = NULL;
|
||||||
|
a = NULL;
|
||||||
|
b = NULL;
|
||||||
|
c0 = NULL;
|
||||||
|
c2 = NULL;
|
||||||
|
c4 = NULL;
|
||||||
|
c6 = NULL;
|
||||||
|
uppercut = NULL;
|
||||||
|
lowercut = NULL;
|
||||||
|
uppercutsq = NULL;
|
||||||
|
lowercutsq = NULL;
|
||||||
|
frho = NULL;
|
||||||
|
rho = NULL;
|
||||||
|
|
||||||
|
// splined arrays
|
||||||
|
frho_spline = NULL;
|
||||||
|
|
||||||
|
// per-atom arrays
|
||||||
|
nmax = 0;
|
||||||
|
fp = NULL;
|
||||||
|
localrho = NULL;
|
||||||
|
|
||||||
|
// set comm size needed by this pair
|
||||||
|
comm_forward = 1;
|
||||||
|
comm_reverse = 1;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
check if allocated, since class can be destructed when incomplete
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
PairLocalDensity::~PairLocalDensity()
|
||||||
|
{
|
||||||
|
|
||||||
|
memory->destroy(localrho);
|
||||||
|
memory->destroy(fp);
|
||||||
|
|
||||||
|
if (allocated) {
|
||||||
|
memory->destroy(setflag);
|
||||||
|
memory->destroy(cutsq);
|
||||||
|
}
|
||||||
|
|
||||||
|
memory->destroy(frho_spline);
|
||||||
|
|
||||||
|
memory->destroy(rho_min);
|
||||||
|
memory->destroy(rho_max);
|
||||||
|
memory->destroy(delta_rho);
|
||||||
|
memory->destroy(c0);
|
||||||
|
memory->destroy(c2);
|
||||||
|
memory->destroy(c4);
|
||||||
|
memory->destroy(c6);
|
||||||
|
memory->destroy(uppercut);
|
||||||
|
memory->destroy(lowercut);
|
||||||
|
memory->destroy(uppercutsq);
|
||||||
|
memory->destroy(lowercutsq);
|
||||||
|
memory->destroy(frho);
|
||||||
|
memory->destroy(rho);
|
||||||
|
|
||||||
|
memory->destroy(a);
|
||||||
|
memory->destroy(b);
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void PairLocalDensity::compute(int eflag, int vflag)
|
||||||
|
{
|
||||||
|
|
||||||
|
int i,j,ii,jj,m,k,inum,jnum,itype,jtype;
|
||||||
|
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
|
||||||
|
double rsqinv, phi, uLD, dphi, evdwl,fpair;
|
||||||
|
double p, *coeff;
|
||||||
|
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||||
|
|
||||||
|
phi = uLD = evdwl = fpair = rsqinv = 0.0;
|
||||||
|
|
||||||
|
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||||
|
else evflag = vflag_fdotr = eflag_global = eflag_atom = 0;
|
||||||
|
|
||||||
|
/* localrho = LD at each atom
|
||||||
|
fp = derivative of embedding energy at each atom for each LD potential
|
||||||
|
uLD = embedding energy of each atom due to each LD potential*/
|
||||||
|
|
||||||
|
// grow LD and fp arrays if necessary
|
||||||
|
// need to be atom->nmax in length
|
||||||
|
|
||||||
|
if (atom->nmax > nmax) {
|
||||||
|
memory->destroy(localrho);
|
||||||
|
memory->destroy(fp);
|
||||||
|
nmax = atom->nmax;
|
||||||
|
memory->create(localrho, nLD, nmax, "pairLD:localrho");
|
||||||
|
memory->create(fp, nLD, nmax, "pairLD:fp");
|
||||||
|
}
|
||||||
|
|
||||||
|
double **x = atom->x;
|
||||||
|
double **f = atom->f;
|
||||||
|
int *type = atom->type;
|
||||||
|
int nlocal = atom->nlocal;
|
||||||
|
int newton_pair = force->newton_pair;
|
||||||
|
|
||||||
|
inum = list->inum;
|
||||||
|
ilist = list->ilist;
|
||||||
|
numneigh = list->numneigh;
|
||||||
|
firstneigh = list->firstneigh;
|
||||||
|
|
||||||
|
// zero out LD and fp
|
||||||
|
|
||||||
|
if (newton_pair) {
|
||||||
|
m = nlocal + atom->nghost;
|
||||||
|
for (k = 0; k < nLD; k++) {
|
||||||
|
for (i = 0; i < m; i++) {
|
||||||
|
localrho[k][i] = 0.0;
|
||||||
|
fp[k][i] = 0.0;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
else {
|
||||||
|
for (k = 0; k < nLD; k++){
|
||||||
|
for (i = 0; i < nlocal; i++) {
|
||||||
|
localrho[k][i] = 0.0;
|
||||||
|
fp[k][i] = 0.0;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// loop over neighs of central atoms and types of LDs
|
||||||
|
|
||||||
|
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];
|
||||||
|
|
||||||
|
// calculate distance-squared between i,j atom-types
|
||||||
|
|
||||||
|
delx = xtmp - x[j][0];
|
||||||
|
dely = ytmp - x[j][1];
|
||||||
|
delz = ztmp - x[j][2];
|
||||||
|
rsq = delx*delx + dely*dely + delz*delz;
|
||||||
|
|
||||||
|
// calculating LDs based on central and neigh filters
|
||||||
|
|
||||||
|
for (k = 0; k < nLD; k++) {
|
||||||
|
if (rsq < lowercutsq[k]) {
|
||||||
|
phi = 1.0;
|
||||||
|
}
|
||||||
|
else if (rsq > uppercutsq[k]) {
|
||||||
|
phi = 0.0;
|
||||||
|
}
|
||||||
|
else {
|
||||||
|
phi = c0[k] + rsq * (c2[k] + rsq * (c4[k] + c6[k]*rsq));
|
||||||
|
}
|
||||||
|
localrho[k][i] += (phi * b[k][jtype]);
|
||||||
|
|
||||||
|
/*checking for both i,j is necessary
|
||||||
|
since a half neighbor list is processed.*/
|
||||||
|
|
||||||
|
if (newton_pair || j<nlocal) {
|
||||||
|
localrho[k][j] += (phi * b[k][itype]);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// communicate and sum LDs over all procs
|
||||||
|
if (newton_pair) comm->reverse_comm_pair(this);
|
||||||
|
|
||||||
|
//
|
||||||
|
|
||||||
|
for (ii = 0; ii < inum; ii++) {
|
||||||
|
i = ilist[ii];
|
||||||
|
itype = type[i];
|
||||||
|
uLD = 0.0;
|
||||||
|
|
||||||
|
for (k = 0; k < nLD; k++) {
|
||||||
|
|
||||||
|
/*skip over this loop if the LD potential
|
||||||
|
is not intendend for central atomtype <itype>*/
|
||||||
|
if (!(a[k][itype])) continue;
|
||||||
|
|
||||||
|
// linear extrapolation at rho_min and rho_max
|
||||||
|
|
||||||
|
if (localrho[k][i] <= rho_min[k]) {
|
||||||
|
coeff = frho_spline[k][0];
|
||||||
|
fp[k][i] = coeff[2];
|
||||||
|
uLD += a[k][itype] * ( coeff[6] + fp[k][i]*(localrho[k][i] - rho_min[k]) );
|
||||||
|
}
|
||||||
|
else if (localrho[k][i] >= rho_max[k]) {
|
||||||
|
coeff = frho_spline[k][nrho-2];
|
||||||
|
fp[k][i] = coeff[0] + coeff[1] + coeff[2];
|
||||||
|
uLD += a[k][itype] * ( (coeff[3] + coeff[4] + coeff[5] + coeff[6]) + fp[k][i]*(localrho[k][i] - rho_max[k]) );
|
||||||
|
}
|
||||||
|
else {
|
||||||
|
p = (localrho[k][i] - rho_min[k]) / delta_rho[k];
|
||||||
|
m = static_cast<int> (p);
|
||||||
|
m = MAX(0, MIN(m, nrho-2));
|
||||||
|
p -= m;
|
||||||
|
p = MIN(p, 1.0);
|
||||||
|
coeff = frho_spline[k][m];
|
||||||
|
fp[k][i] = (coeff[0]*p + coeff[1])*p + coeff[2];
|
||||||
|
uLD += a[k][itype] * (((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
if (eflag) {
|
||||||
|
if (eflag_global) eng_vdwl += uLD;
|
||||||
|
if (eflag_atom) eatom[i] += uLD;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// communicate LD and fp to all procs
|
||||||
|
|
||||||
|
comm->forward_comm_pair(this);
|
||||||
|
|
||||||
|
// compute forces on 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];
|
||||||
|
|
||||||
|
// calculate square of distance between i,j atoms
|
||||||
|
|
||||||
|
delx = xtmp - x[j][0];
|
||||||
|
dely = ytmp - x[j][1];
|
||||||
|
delz = ztmp - x[j][2];
|
||||||
|
rsq = delx*delx + dely*dely + delz*delz;
|
||||||
|
|
||||||
|
// calculate force between two atoms
|
||||||
|
fpair = 0.0;
|
||||||
|
if (rsq < cutforcesq) { // global cutoff check
|
||||||
|
rsqinv = 1.0/rsq;
|
||||||
|
for (k = 0; k < nLD; k++) {
|
||||||
|
if (rsq >= lowercutsq[k] && rsq < uppercutsq[k]) {
|
||||||
|
dphi = rsq * (2.0*c2[k] + rsq * (4.0*c4[k] + 6.0*c6[k]*rsq));
|
||||||
|
fpair += -(a[k][itype]*b[k][jtype]*fp[k][i] + a[k][jtype]*b[k][itype]*fp[k][j]) * dphi;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
fpair *= rsqinv;
|
||||||
|
|
||||||
|
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;
|
||||||
|
}
|
||||||
|
|
||||||
|
/*eng_vdwl has already been completely built,
|
||||||
|
so no need to add anything here*/
|
||||||
|
|
||||||
|
if (eflag) evdwl = 0.0;
|
||||||
|
|
||||||
|
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 PairLocalDensity::allocate()
|
||||||
|
{
|
||||||
|
allocated = 1;
|
||||||
|
int n = atom->ntypes;
|
||||||
|
|
||||||
|
memory->create(cutsq,n+1,n+1,"pair:cutsq");
|
||||||
|
|
||||||
|
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;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
global settings
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void PairLocalDensity::settings(int narg, char **arg)
|
||||||
|
{
|
||||||
|
if (narg > 0) error->all(FLERR,"Illegal pair_style command");
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
set coeffs for all type pairs
|
||||||
|
read tabulated LD input file
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void PairLocalDensity::coeff(int narg, char **arg)
|
||||||
|
{
|
||||||
|
int i, j;
|
||||||
|
if (!allocated) allocate();
|
||||||
|
|
||||||
|
if (narg != 3) error->all(FLERR,"Incorrect args for pair coefficients");
|
||||||
|
|
||||||
|
// insure I,J args are * *
|
||||||
|
|
||||||
|
if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
|
||||||
|
error->all(FLERR,"Incorrect args for pair coefficients");
|
||||||
|
|
||||||
|
// parse LD file
|
||||||
|
|
||||||
|
parse_file(arg[2]);
|
||||||
|
|
||||||
|
|
||||||
|
// clear setflag since coeff() called once with I,J = * *
|
||||||
|
|
||||||
|
for (i = 1; i <= atom->ntypes; i++)
|
||||||
|
for (j = i; j <= atom->ntypes; j++)
|
||||||
|
setflag[i][j] = 0;
|
||||||
|
|
||||||
|
// set setflag for all i,j type pairs
|
||||||
|
|
||||||
|
int count = 0;
|
||||||
|
for (i = 1; i <= atom->ntypes; i++) {
|
||||||
|
for (j = i; j <= atom->ntypes; j++) {
|
||||||
|
setflag[i][j] = 1;
|
||||||
|
count++;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
init specific to this pair style
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void PairLocalDensity::init_style()
|
||||||
|
{
|
||||||
|
// spline rho and frho arrays
|
||||||
|
// request half neighbor list
|
||||||
|
|
||||||
|
array2spline();
|
||||||
|
|
||||||
|
// half neighbor request
|
||||||
|
neighbor->request(this);
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
init for one type pair i,j and corresponding j,i
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
double PairLocalDensity::init_one(int i, int j)
|
||||||
|
{
|
||||||
|
// single global cutoff = max of all uppercuts read in from LD file
|
||||||
|
|
||||||
|
cutmax = 0.0;
|
||||||
|
for (int k = 0; k < nLD; k++)
|
||||||
|
cutmax = MAX(cutmax,uppercut[k]);
|
||||||
|
|
||||||
|
cutforcesq = cutmax*cutmax;
|
||||||
|
|
||||||
|
return cutmax;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
/*--------------------------------------------------------------------------
|
||||||
|
pair_write functionality for this pair style that gives just a snap-shot
|
||||||
|
of the LD potential without doing an actual MD run
|
||||||
|
---------------------------------------------------------------------------*/
|
||||||
|
|
||||||
|
double PairLocalDensity::single(int i, int j, int itype, int jtype, double rsq,
|
||||||
|
double factor_coul, double factor_lj,
|
||||||
|
double &fforce)
|
||||||
|
{
|
||||||
|
int m, k, index;
|
||||||
|
double rsqinv, p, uLD;
|
||||||
|
double *coeff, **LD;
|
||||||
|
double dFdrho, phi, dphi;
|
||||||
|
|
||||||
|
uLD = dFdrho = dphi = 0.0;
|
||||||
|
|
||||||
|
memory->create(LD, nLD, 3, "pairLD:LD");
|
||||||
|
for (k = 0; k < nLD; k++) {
|
||||||
|
LD[k][1] = 0.0; // itype:- 1
|
||||||
|
LD[k][2] = 0.0; // jtype:- 2
|
||||||
|
}
|
||||||
|
|
||||||
|
rsqinv = 1.0/rsq;
|
||||||
|
for (k = 0; k < nLD; k++) {
|
||||||
|
if (rsq < lowercutsq[k]) {
|
||||||
|
phi = 1.0;
|
||||||
|
}
|
||||||
|
else if (rsq > uppercutsq[k]) {
|
||||||
|
phi = 0.0;
|
||||||
|
}
|
||||||
|
else {
|
||||||
|
phi = c0[k] + rsq * (c2[k] + rsq * (c4[k] + c6[k]*rsq));
|
||||||
|
}
|
||||||
|
LD[k][1] += (phi * b[k][jtype]);
|
||||||
|
LD[k][2] += (phi * b[k][itype]);
|
||||||
|
}
|
||||||
|
|
||||||
|
for (k = 0; k < nLD; k++) {
|
||||||
|
if (a[k][itype]) index = 1;
|
||||||
|
if (a[k][jtype]) index = 2;
|
||||||
|
|
||||||
|
if (LD[k][index] <= rho_min[k]) {
|
||||||
|
coeff = frho_spline[k][0];
|
||||||
|
dFdrho = coeff[2];
|
||||||
|
uLD += a[k][itype] * ( coeff[6] + dFdrho*(LD[k][index] - rho_min[k]) );
|
||||||
|
}
|
||||||
|
else if (LD[k][index] >= rho_max[k]) {
|
||||||
|
coeff = frho_spline[k][nrho-1];
|
||||||
|
dFdrho = coeff[0] + coeff[1] + coeff[2];
|
||||||
|
uLD += a[k][itype] * ( (coeff[3] + coeff[4] + coeff[5] + coeff[6]) + dFdrho*(LD[k][index] - rho_max[k]) );
|
||||||
|
}
|
||||||
|
else {
|
||||||
|
p = (LD[k][index] - rho_min[k]) / delta_rho[k];
|
||||||
|
m = static_cast<int> (p);
|
||||||
|
m = MAX(0, MIN(m, nrho-2));
|
||||||
|
p -= m;
|
||||||
|
p = MIN(p, 1.0);
|
||||||
|
coeff = frho_spline[k][m];
|
||||||
|
dFdrho = (coeff[0]*p + coeff[1])*p + coeff[2];
|
||||||
|
uLD += a[k][itype] * (((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]);
|
||||||
|
}
|
||||||
|
|
||||||
|
if (rsq < lowercutsq[k]) {
|
||||||
|
dphi = 0.0;
|
||||||
|
}
|
||||||
|
else if (rsq > uppercutsq[k]) {
|
||||||
|
dphi = 0.0;
|
||||||
|
}
|
||||||
|
else {
|
||||||
|
dphi = rsq * (2.0*c2[k] + rsq * (4.0*c4[k] + 6.0*c6[k]*rsq));
|
||||||
|
}
|
||||||
|
fforce += -(a[k][itype]*b[k][jtype]*dFdrho + a[k][jtype]*b[k][itype]*dFdrho) * dphi *rsqinv;
|
||||||
|
}
|
||||||
|
memory->destroy(LD);
|
||||||
|
|
||||||
|
return uLD;
|
||||||
|
}
|
||||||
|
|
||||||
|
/*--------------------------------------------------------------------
|
||||||
|
spline the array frho read in from the file to create
|
||||||
|
frho_spline
|
||||||
|
---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void PairLocalDensity::array2spline() {
|
||||||
|
memory->destroy(frho_spline);
|
||||||
|
memory->create(frho_spline, nLD, nrho, 7, "pairLD:frho_spline");
|
||||||
|
|
||||||
|
for (int k = 0; k < nLD; k++)
|
||||||
|
interpolate_cbspl(nrho, delta_rho[k], frho[k], frho_spline[k]);
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
(one-dimensional) cubic spline interpolation sub-routine,
|
||||||
|
which determines the coeffs for a clamped cubic spline
|
||||||
|
given tabulated data
|
||||||
|
------------------------------------------------------------------------*/
|
||||||
|
|
||||||
|
void PairLocalDensity::interpolate_cbspl(int n, double delta,
|
||||||
|
double *f, double **spline)
|
||||||
|
{
|
||||||
|
/* inputs:
|
||||||
|
n number of interpolating points
|
||||||
|
|
||||||
|
f array containing function values to
|
||||||
|
be interpolated; f[i] is the function
|
||||||
|
value corresponding to x[i]
|
||||||
|
('x' refers to the independent var)
|
||||||
|
|
||||||
|
delta difference in tabulated values of x
|
||||||
|
|
||||||
|
outputs: (packaged as columns of the coeff matrix)
|
||||||
|
coeff_b coeffs of linear terms
|
||||||
|
coeff_c coeffs of quadratic terms
|
||||||
|
coeff_d coeffs of cubic terms
|
||||||
|
spline matrix that collects b,c,d
|
||||||
|
|
||||||
|
|
||||||
|
other parameters:
|
||||||
|
fpa derivative of function at x=a
|
||||||
|
fpb derivative of function at x=b
|
||||||
|
*/
|
||||||
|
|
||||||
|
double *dl, *dd, *du;
|
||||||
|
double *coeff_b, *coeff_c, *coeff_d;
|
||||||
|
double fpa, fpb;
|
||||||
|
|
||||||
|
int i;
|
||||||
|
|
||||||
|
coeff_b = new double [n];
|
||||||
|
coeff_c = new double [n];
|
||||||
|
coeff_d = new double [n];
|
||||||
|
dl = new double [n];
|
||||||
|
dd = new double [n];
|
||||||
|
du = new double [n];
|
||||||
|
|
||||||
|
// initialize values
|
||||||
|
for ( i = 0; i<n; i++) {
|
||||||
|
coeff_b[i] = coeff_c[i] = coeff_d[i] = 0.;
|
||||||
|
dl[i] = dd[i] = du[i] = 0.;
|
||||||
|
}
|
||||||
|
|
||||||
|
// set slopes at beginning and end
|
||||||
|
fpa = 0.;
|
||||||
|
fpb = 0.;
|
||||||
|
|
||||||
|
for ( i = 0; i < n-1; i++ ) {
|
||||||
|
dl[i] = du[i] = delta;
|
||||||
|
}
|
||||||
|
|
||||||
|
dd[0] = 2.0 * delta;
|
||||||
|
dd[n-1] = 2.0 * delta;
|
||||||
|
coeff_c[0] = ( 3.0 / delta ) * ( f[1] - f[0] ) - 3.0 * fpa;
|
||||||
|
coeff_c[n-1] = 3.0 * fpb - ( 3.0 / delta ) * ( f[n-1] - f[n-2] );
|
||||||
|
for ( i = 0; i < n-2; i++ ) {
|
||||||
|
dd[i+1] = 4.0 * delta;
|
||||||
|
coeff_c[i+1] = ( 3.0 / delta ) * ( f[i+2] - f[i+1] ) -
|
||||||
|
( 3.0 / delta ) * ( f[i+1] - f[i] );
|
||||||
|
}
|
||||||
|
|
||||||
|
// tridiagonal solver
|
||||||
|
for ( i = 0; i < n-1; i++ ) {
|
||||||
|
du[i] /= dd[i];
|
||||||
|
dd[i+1] -= dl[i]*du[i];
|
||||||
|
}
|
||||||
|
|
||||||
|
coeff_c[0] /= dd[0];
|
||||||
|
for ( i = 1; i < n; i++ )
|
||||||
|
coeff_c[i] = ( coeff_c[i] - dl[i-1] * coeff_c[i-1] ) / dd[i];
|
||||||
|
|
||||||
|
for ( i = n-2; i >= 0; i-- )
|
||||||
|
coeff_c[i] -= coeff_c[i+1] * du[i];
|
||||||
|
|
||||||
|
for ( i = 0; i < n-1; i++ ) {
|
||||||
|
coeff_d[i] = ( coeff_c[i+1] - coeff_c[i] ) / ( 3.0 * delta );
|
||||||
|
coeff_b[i] = ( f[i+1] - f[i] ) / delta - delta * ( coeff_c[i+1] + 2.0*coeff_c[i] ) / 3.0;
|
||||||
|
}
|
||||||
|
|
||||||
|
// normalize
|
||||||
|
for ( i = 0; i < n-1; i++ ) {
|
||||||
|
coeff_b[i] = coeff_b[i] * delta ;
|
||||||
|
coeff_c[i] = coeff_c[i] * delta*delta ;
|
||||||
|
coeff_d[i] = coeff_d[i] * delta*delta*delta;
|
||||||
|
}
|
||||||
|
|
||||||
|
//copy to coefficient matrix
|
||||||
|
for ( i = 0; i < n; i++) {
|
||||||
|
spline[i][3] = coeff_d[i];
|
||||||
|
spline[i][4] = coeff_c[i];
|
||||||
|
spline[i][5] = coeff_b[i];
|
||||||
|
spline[i][6] = f[i];
|
||||||
|
spline[i][2] = spline[i][5]/delta;
|
||||||
|
spline[i][1] = 2.0*spline[i][4]/delta;
|
||||||
|
spline[i][0] = 3.0*spline[i][3]/delta;
|
||||||
|
}
|
||||||
|
|
||||||
|
delete [] coeff_b;
|
||||||
|
delete [] coeff_c;
|
||||||
|
delete [] coeff_d;
|
||||||
|
delete [] du;
|
||||||
|
delete [] dd;
|
||||||
|
delete [] dl;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
read potential values from tabulated LD input file
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void PairLocalDensity::parse_file(char *filename) {
|
||||||
|
|
||||||
|
int k, n;
|
||||||
|
int me = comm->me;
|
||||||
|
FILE *fptr;
|
||||||
|
char line[MAXLINE];
|
||||||
|
double ratio, lc2, uc2, denom;
|
||||||
|
|
||||||
|
|
||||||
|
if (me == 0) {
|
||||||
|
fptr = fopen(filename, "r");
|
||||||
|
if (fptr == NULL) {
|
||||||
|
char str[128];
|
||||||
|
sprintf(str,"Cannot open Local Density potential file %s",filename);
|
||||||
|
error->one(FLERR,str);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
double *ftmp; // tmp var to extract the complete 2D frho array from file
|
||||||
|
|
||||||
|
// broadcast number of LD potentials and number of (rho,frho) pairs
|
||||||
|
if (me == 0) {
|
||||||
|
|
||||||
|
// first 2 comment lines ignored
|
||||||
|
fgets(line,MAXLINE,fptr);
|
||||||
|
fgets(line,MAXLINE,fptr);
|
||||||
|
|
||||||
|
// extract number of potentials and number of (frho, rho) points
|
||||||
|
fgets(line,MAXLINE,fptr);
|
||||||
|
sscanf(line, "%d %d", &nLD, &nrho);
|
||||||
|
fgets(line,MAXLINE,fptr);
|
||||||
|
}
|
||||||
|
|
||||||
|
MPI_Bcast(&nLD,1,MPI_INT,0,world);
|
||||||
|
MPI_Bcast(&nrho,1,MPI_INT,0,world);
|
||||||
|
|
||||||
|
// setting up all arrays to be read from files and broadcasted
|
||||||
|
memory->create(uppercut, nLD, "pairLD:uppercut");
|
||||||
|
memory->create(lowercut, nLD, "pairLD:lowercut");
|
||||||
|
memory->create(uppercutsq, nLD, "pairLD:uppercutsq");
|
||||||
|
memory->create(lowercutsq, nLD, "pairLD:lowercutsq");
|
||||||
|
memory->create(c0, nLD, "pairLD:c0");
|
||||||
|
memory->create(c2, nLD, "pairLD:c2");
|
||||||
|
memory->create(c4, nLD, "pairLD:c4");
|
||||||
|
memory->create(c6, nLD, "pairLD:c6");
|
||||||
|
memory->create(rho_min, nLD, "pairLD:rho_min");
|
||||||
|
memory->create(rho_max, nLD, "pairLD:rho_max");
|
||||||
|
memory->create(delta_rho, nLD,"pairLD:delta_rho");
|
||||||
|
memory->create(ftmp, nrho*nLD, "pairLD:ftmp");
|
||||||
|
|
||||||
|
// setting up central and neighbor atom filters
|
||||||
|
memory->create(a, nLD, atom->ntypes+1 , "pairLD:a");
|
||||||
|
memory->create(b, nLD, atom->ntypes+1, "pairLD:b");
|
||||||
|
if (me == 0) {
|
||||||
|
for (n = 1; n <= atom->ntypes; n++){
|
||||||
|
for (k = 0; k < nLD; k++) {
|
||||||
|
a[k][n] = 0;
|
||||||
|
b[k][n] = 0;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// read file block by block
|
||||||
|
|
||||||
|
if (me == 0) {
|
||||||
|
for (k = 0; k < nLD; k++) {
|
||||||
|
|
||||||
|
// parse upper and lower cut values
|
||||||
|
if (fgets(line,MAXLINE,fptr)==NULL) break;
|
||||||
|
sscanf(line, "%lf %lf", &lowercut[k], &uppercut[k]);
|
||||||
|
|
||||||
|
// parse and broadcast central atom filter
|
||||||
|
fgets(line, MAXLINE, fptr);
|
||||||
|
char *tmp = strtok(line, " /t/n/r/f");
|
||||||
|
while (tmp != NULL) {
|
||||||
|
a[k][atoi(tmp)] = 1;
|
||||||
|
tmp = strtok(NULL, " /t/n/r/f");
|
||||||
|
}
|
||||||
|
|
||||||
|
// parse neighbor atom filter
|
||||||
|
fgets(line, MAXLINE, fptr);
|
||||||
|
tmp = strtok(line, " /t/n/r/f");
|
||||||
|
while (tmp != NULL) {
|
||||||
|
b[k][atoi(tmp)] = 1;
|
||||||
|
tmp = strtok(NULL, " /t/n/r/f");
|
||||||
|
}
|
||||||
|
|
||||||
|
// parse min, max and delta rho values
|
||||||
|
fgets(line, MAXLINE, fptr);
|
||||||
|
sscanf(line, "%lf %lf %lf", &rho_min[k], &rho_max[k], &delta_rho[k]);
|
||||||
|
// recompute delta_rho from scratch for precision
|
||||||
|
delta_rho[k] = (rho_max[k] - rho_min[k]) / (nrho - 1);
|
||||||
|
|
||||||
|
// parse tabulated frho values from each line into temporary array
|
||||||
|
for (n = 0; n < nrho; n++) {
|
||||||
|
fgets(line,MAXLINE,fptr);
|
||||||
|
sscanf(line, "%lf", &ftmp[k*nrho + n]);
|
||||||
|
}
|
||||||
|
|
||||||
|
// ignore blank line at the end of every block
|
||||||
|
fgets(line,MAXLINE,fptr);
|
||||||
|
|
||||||
|
// set coefficients for local density indicator function
|
||||||
|
uc2 = uppercut[k] * uppercut[k];
|
||||||
|
uppercutsq[k] = uc2;
|
||||||
|
lc2 = lowercut[k] * lowercut[k];
|
||||||
|
lowercutsq[k] = lc2;
|
||||||
|
ratio = lc2/uc2;
|
||||||
|
denom = 1.0 - ratio;
|
||||||
|
denom = denom*denom*denom;
|
||||||
|
c0[k] = (1 - 3.0 * ratio) / denom;
|
||||||
|
c2[k] = (6.0 * ratio) / (uc2 * denom);
|
||||||
|
c4[k] = -(3.0 + 3.0*ratio) / (uc2*uc2 * denom);
|
||||||
|
c6[k] = 2.0 / (uc2*uc2*uc2 * denom);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// Broadcast all parsed arrays
|
||||||
|
MPI_Bcast(&lowercut[0], nLD, MPI_DOUBLE, 0, world);
|
||||||
|
MPI_Bcast(&uppercut[0], nLD, MPI_DOUBLE, 0, world);
|
||||||
|
MPI_Bcast(&lowercutsq[0], nLD, MPI_DOUBLE, 0, world);
|
||||||
|
MPI_Bcast(&uppercutsq[0], nLD, MPI_DOUBLE, 0, world);
|
||||||
|
MPI_Bcast(&c0[0], nLD, MPI_DOUBLE, 0, world);
|
||||||
|
MPI_Bcast(&c2[0], nLD, MPI_DOUBLE, 0, world);
|
||||||
|
MPI_Bcast(&c4[0], nLD, MPI_DOUBLE, 0, world);
|
||||||
|
MPI_Bcast(&c6[0], nLD, MPI_DOUBLE, 0, world);
|
||||||
|
for (k = 0; k < nLD; k++) {
|
||||||
|
MPI_Bcast(&a[k][1], atom->ntypes, MPI_INT, 0, world);
|
||||||
|
MPI_Bcast(&b[k][1], atom->ntypes, MPI_INT, 0, world);
|
||||||
|
}
|
||||||
|
MPI_Bcast(&rho_min[0], nLD, MPI_DOUBLE, 0, world);
|
||||||
|
MPI_Bcast(&rho_max[0], nLD, MPI_DOUBLE, 0, world);
|
||||||
|
MPI_Bcast(&delta_rho[0], nLD, MPI_DOUBLE, 0, world);
|
||||||
|
MPI_Bcast(&ftmp[0], nLD*nrho, MPI_DOUBLE, 0, world);
|
||||||
|
|
||||||
|
if (me == 0) fclose(fptr);
|
||||||
|
|
||||||
|
// set up rho and frho arrays
|
||||||
|
memory->create(rho, nLD, nrho, "pairLD:rho");
|
||||||
|
memory->create(frho, nLD, nrho, "pairLD:frho");
|
||||||
|
|
||||||
|
for (k = 0; k < nLD; k++) {
|
||||||
|
for (n = 0; n < nrho; n++) {
|
||||||
|
rho[k][n] = rho_min[k] + n*delta_rho[k];
|
||||||
|
frho[k][n] = ftmp[k*nrho + n];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// delete temporary array
|
||||||
|
memory->destroy(ftmp);
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
communication routines
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
int PairLocalDensity::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) {
|
||||||
|
int i,j,k;
|
||||||
|
int m;
|
||||||
|
|
||||||
|
m = 0;
|
||||||
|
for (i = 0; i < n; i++) {
|
||||||
|
j = list[i];
|
||||||
|
for (k = 0; k < nLD; k++) {
|
||||||
|
buf[m++] = fp[k][j];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
return nLD;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void PairLocalDensity::unpack_comm(int n, int first, double *buf) {
|
||||||
|
|
||||||
|
int i,k,m,last;
|
||||||
|
|
||||||
|
m = 0;
|
||||||
|
last = first + n;
|
||||||
|
for (i = first; i < last; i++) {
|
||||||
|
for (k = 0; k < nLD; k++) {
|
||||||
|
fp[k][i] = buf[m++];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
int PairLocalDensity::pack_reverse_comm(int n, int first, double *buf) {
|
||||||
|
|
||||||
|
int i,k,m,last;
|
||||||
|
|
||||||
|
m = 0;
|
||||||
|
last = first + n;
|
||||||
|
for (i = first; i < last; i++) {
|
||||||
|
for (k = 0; k < nLD; k++) {
|
||||||
|
buf[m++] = localrho[k][i];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return nLD;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ---------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
void PairLocalDensity::unpack_reverse_comm(int n, int *list, double *buf) {
|
||||||
|
|
||||||
|
int i,j,k;
|
||||||
|
int m;
|
||||||
|
|
||||||
|
m = 0;
|
||||||
|
for (i = 0; i < n; i++) {
|
||||||
|
j = list[i];
|
||||||
|
for (k = 0; k < nLD; k++) {
|
||||||
|
localrho[k][j] += buf[m++];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
memory usage of local atom-based arrays
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
double PairLocalDensity::memory_usage()
|
||||||
|
{
|
||||||
|
double bytes = maxeatom * sizeof(double);
|
||||||
|
bytes += maxvatom*6 * sizeof(double);
|
||||||
|
bytes += 2 * (nmax*nLD) * sizeof(double);
|
||||||
|
return bytes;
|
||||||
|
}
|
||||||
|
|
||||||
88
src/USER-MISC/pair_local_density.h
Normal file
88
src/USER-MISC/pair_local_density.h
Normal file
@ -0,0 +1,88 @@
|
|||||||
|
/* -*- c++ -*- ----------------------------------------------------------
|
||||||
|
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.
|
||||||
|
-------------------------------------------------------------------------
|
||||||
|
pair_LocalDensity written by:
|
||||||
|
Tanmoy Sanyal and M. Scott Shell from UC Santa Barbara
|
||||||
|
David Rosenberger: TU Darmstadt
|
||||||
|
-------------------------------------------------------------------------*/
|
||||||
|
|
||||||
|
|
||||||
|
#ifdef PAIR_CLASS
|
||||||
|
|
||||||
|
PairStyle(local/density,PairLocalDensity)
|
||||||
|
|
||||||
|
#else
|
||||||
|
|
||||||
|
#ifndef LMP_PAIR_LOCAL_DENSITY_H
|
||||||
|
#define LMP_PAIR_LOCAL_DENSITY_H
|
||||||
|
|
||||||
|
#include "pair.h"
|
||||||
|
|
||||||
|
|
||||||
|
namespace LAMMPS_NS {
|
||||||
|
|
||||||
|
class PairLocalDensity : public Pair {
|
||||||
|
public:
|
||||||
|
PairLocalDensity(class LAMMPS *);
|
||||||
|
virtual ~PairLocalDensity();
|
||||||
|
virtual void compute(int, int);
|
||||||
|
void settings(int, char **);
|
||||||
|
virtual void coeff(int, char **);
|
||||||
|
void init_style();
|
||||||
|
double init_one(int, int);
|
||||||
|
double single(int, int, int, int, double, double, double, double &);
|
||||||
|
|
||||||
|
virtual int pack_comm(int, int *, double *, int, int *);
|
||||||
|
virtual void unpack_comm(int, int, double *);
|
||||||
|
int pack_reverse_comm(int, int, double *);
|
||||||
|
void unpack_reverse_comm(int, int *, double *);
|
||||||
|
double memory_usage();
|
||||||
|
|
||||||
|
|
||||||
|
protected:
|
||||||
|
//------------------------------------------------------------------------
|
||||||
|
//This information is read from the tabulated input file
|
||||||
|
|
||||||
|
int nLD, nrho; // number of LD types
|
||||||
|
int **a, **b; // central and neigh atom filters
|
||||||
|
double *uppercut, *lowercut; // upper and lower cutoffs
|
||||||
|
double *uppercutsq, *lowercutsq; // square of above cutoffs
|
||||||
|
double *c0, *c2, *c4, *c6; // coeffs for indicator function
|
||||||
|
double *rho_min, *rho_max, *delta_rho; // min, max & grid-size for LDs
|
||||||
|
double **rho, **frho; // LD and LD function tables
|
||||||
|
|
||||||
|
//------------------------------------------------------------------------
|
||||||
|
|
||||||
|
double ***frho_spline; // splined LD potentials
|
||||||
|
double cutmax; // max cutoff for all elements
|
||||||
|
double cutforcesq; // square of global upper cutoff
|
||||||
|
|
||||||
|
int nmax; // max size of per-atom arrays
|
||||||
|
double **localrho; // per-atom LD
|
||||||
|
double **fp; // per-atom LD potential function derivative
|
||||||
|
|
||||||
|
void allocate();
|
||||||
|
|
||||||
|
// read tabulated input file
|
||||||
|
void parse_file(char *);
|
||||||
|
|
||||||
|
// convert array to spline
|
||||||
|
void array2spline();
|
||||||
|
|
||||||
|
// cubic spline interpolation
|
||||||
|
void interpolate_cbspl(int, double, double *, double **);
|
||||||
|
};
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
#endif
|
||||||
|
#endif
|
||||||
Reference in New Issue
Block a user